💻 Fortran Source Code Library
We currently offer 172 open-source, production-grade Fortran codes for offline testing. Run calculations locally on your own machine, view code structure, read technical explanations, and download compilation packages including sample input files.
1D Transient Conduction FDM Solver
Core Numerical Engine in Fortran 90 • 58 total downloads
transient_fdm_solver.f90
! =========================================================================
! Source File: transient_fdm_solver.f90
! =========================================================================
program transient_fdm_solver
implicit none
! Variable declarations
integer :: geom_type, n_nodes, method
integer :: bc_left_type, bc_right_type, init_type
real(8) :: length, dt_choice, t_sim
real(8) :: k_0, rho, cp_0, beta_k, beta_cp
real(8) :: bc_left_p1, bc_left_p2, bc_right_p1, bc_right_p2
real(8) :: init_p1, init_p2
! Grid and solver variables
real(8), allocatable, dimension(:) :: T, T_old, r_nodes, V_nodes
real(8), allocatable, dimension(:) :: a_coeff, b_coeff, c_coeff, d_coeff
real(8), allocatable, dimension(:) :: k_nodes, cp_nodes
real(8) :: dx, dt, alpha, current_time
integer :: i, step, total_steps, print_interval
character(len=20) :: geom_name, method_name
real(8) :: E_in_total, E_out_total, E_st_initial, E_st_current, E_st_change
real(8) :: balance_err, heat_flux_left, heat_flux_right
real(8) :: dt_max_stable, Fo_max, Fo_node, Bi_node
logical :: is_stable
! Read parameters from stdin
read *, geom_type ! 1=Plane, 2=Cylinder, 3=Sphere
read *, length ! Total thickness L (m) or outer radius r0 (m)
read *, n_nodes ! Number of nodes
read *, dt_choice ! <= 0 for auto stability, > 0 for manual time step (s)
read *, method ! 1=Explicit, 2=Implicit BE, 3=Crank-Nicolson
read *, k_0 ! Thermal conductivity k0 [W/m-K]
read *, rho ! Density rho [kg/m3]
read *, cp_0 ! Specific heat cp0 [J/kg-K]
read *, beta_k ! Temperature coeff for conductivity
read *, beta_cp ! Temperature coeff for specific heat
read *, bc_left_type ! Left BC: 1=Const T, 2=Const q", 3=Convection, 4=Insulated
read *, bc_left_p1 ! T_s, q", or h
read *, bc_left_p2 ! T_inf (if convection)
read *, bc_right_type ! Right BC: 1=Const T, 2=Const q", 3=Convection, 4=Insulated
read *, bc_right_p1 ! T_s, q", or h
read *, bc_right_p2 ! T_inf (if convection)
read *, init_type ! 1=Uniform, 2=Linear, 3=Custom (parabolic)
read *, init_p1 ! Ti or T_left or T_center
read *, init_p2 ! T_right or T_outer
read *, t_sim ! Total simulation time [s]
! Allocations
allocate(T(n_nodes))
allocate(T_old(n_nodes))
allocate(r_nodes(n_nodes))
allocate(V_nodes(n_nodes))
allocate(a_coeff(n_nodes))
allocate(b_coeff(n_nodes))
allocate(c_coeff(n_nodes))
allocate(d_coeff(n_nodes))
allocate(k_nodes(n_nodes))
allocate(cp_nodes(n_nodes))
! Grid setup
dx = length / dble(n_nodes - 1)
do i = 1, n_nodes
r_nodes(i) = dble(i - 1) * dx
end do
! Calculate volumes for nodes (for energy balance check)
select case (geom_type)
case (1)
geom_name = "1D Plane Wall"
do i = 1, n_nodes
if (i == 1 .or. i == n_nodes) then
V_nodes(i) = dx / 2.0d0
else
V_nodes(i) = dx
end if
end do
case (2)
geom_name = "1D Cylinder"
do i = 1, n_nodes
if (i == 1) then
V_nodes(i) = 3.14159265358979d0 * (dx / 2.0d0)**2
elseif (i == n_nodes) then
V_nodes(i) = 3.14159265358979d0 * (length**2 - (length - dx/2.0d0)**2)
else
V_nodes(i) = 3.14159265358979d0 * ((r_nodes(i) + dx/2.0d0)**2 - (r_nodes(i) - dx/2.0d0)**2)
end if
end do
case (3)
geom_name = "1D Sphere"
do i = 1, n_nodes
if (i == 1) then
V_nodes(i) = (4.0d0/3.0d0) * 3.14159265358979d0 * (dx / 2.0d0)**3
elseif (i == n_nodes) then
V_nodes(i) = (4.0d0/3.0d0) * 3.14159265358979d0 * (length**3 - (length - dx/2.0d0)**3)
else
V_nodes(i) = (4.0d0/3.0d0) * 3.14159265358979d0 * ((r_nodes(i) + dx/2.0d0)**3 - (r_nodes(i) - dx/2.0d0)**3)
end if
end do
end select
! Initial temperature distribution
do i = 1, n_nodes
select case (init_type)
case (1) ! Uniform
T(i) = init_p1
case (2) ! Linear
T(i) = init_p1 + (init_p2 - init_p1) * (r_nodes(i) / length)
case (3) ! Parabolic
T(i) = init_p1 - (init_p1 - init_p2) * (r_nodes(i) / length)**2
end select
end do
! Estimate stability limit for Explicit method
! We evaluate at initial temperature (which is safe)
call get_properties(T(1), k_0, cp_0, beta_k, beta_cp, k_nodes(1), cp_nodes(1))
alpha = k_nodes(1) / (rho * cp_nodes(1))
select case (geom_type)
case (1)
dt_max_stable = (dx**2) / (2.0d0 * alpha)
! Adjust for convection if present
if (bc_left_type == 3) then
dt_max_stable = min(dt_max_stable, (dx**2) / (2.0d0 * alpha * (1.0d0 + bc_left_p1 * dx / k_nodes(1))))
endif
if (bc_right_type == 3) then
dt_max_stable = min(dt_max_stable, (dx**2) / (2.0d0 * alpha * (1.0d0 + bc_right_p1 * dx / k_nodes(1))))
endif
case (2)
! Cylinder center node has 4*Fo limit
dt_max_stable = (dx**2) / (4.0d0 * alpha)
case (3)
! Sphere center node has 6*Fo limit
dt_max_stable = (dx**2) / (6.0d0 * alpha)
end select
! Time step selection
if (dt_choice <= 0.0d0) then
dt = 0.95d0 * dt_max_stable
else
dt = dt_choice
end if
! Determine method name
select case (method)
case (1)
method_name = "Explicit (Euler)"
case (2)
method_name = "Implicit (Backward)"
case (3)
method_name = "Crank-Nicolson"
end select
total_steps = ceiling(t_sim / dt)
if (total_steps < 1) total_steps = 1
! Adjust print interval to output ~10 profiles in text report
print_interval = max(1, total_steps / 10)
! Initial energy stored
E_st_initial = 0.0d0
do i = 1, n_nodes
call get_properties(T(i), k_0, cp_0, beta_k, beta_cp, k_nodes(i), cp_nodes(i))
E_st_initial = E_st_initial + rho * cp_nodes(i) * T(i) * V_nodes(i)
end do
E_in_total = 0.0d0
E_out_total = 0.0d0
current_time = 0.0d0
! ==========================================
! OUTPUT HEADER REPORT
! ==========================================
print *, '==================================================='
print *, ' 1D TRANSIENT CONDUCTION NUMERICAL SOLVER REPORT'
print *, '==================================================='
print *, ''
print '(A,A)', ' Geometry Configuration: ', trim(geom_name)
print '(A,F10.6,A)', ' Domain Length/Radius (L): ', length, ' m'
print '(A,I5)', ' Number of Nodal Points (N): ', n_nodes
print '(A,F10.6,A)', ' Nodal Grid Spacing (dx): ', dx, ' m'
print '(A,A)', ' Time-Marching Algorithm: ', trim(method_name)
print '(A,F12.6,A)', ' Time Step Size (dt): ', dt, ' s'
print '(A,I8)', ' Total Number of Steps: ', total_steps
print '(A,F12.4,A)', ' Total Simulation Time: ', t_sim, ' s'
print *, ''
print *, '==================================================='
print *, ' MATERIAL PROPERTIES & BOUNDARY CONDITIONS'
print *, '==================================================='
print '(A,F10.4,A)', ' Thermal Conductivity k0: ', k_0, ' W/(m-K)'
if (beta_k /= 0.0d0) then
print '(A,ES12.4,A)', ' Temperature Coeff (beta_k):', beta_k, ' 1/K'
endif
print '(A,F10.2,A)', ' Density (rho): ', rho, ' kg/m3'
print '(A,F10.2,A)', ' Specific Heat Cp0: ', cp_0, ' J/(kg-K)'
if (beta_cp /= 0.0d0) then
print '(A,ES12.4,A)', ' Temperature Coeff (beta_cp):', beta_cp, ' 1/K'
endif
print *, ''
! Boundary Condition displays
call print_bc_desc("Left/Inner (x=0)", bc_left_type, bc_left_p1, bc_left_p2)
call print_bc_desc("Right/Outer (x=L)", bc_right_type, bc_right_p1, bc_right_p2)
print *, ''
! Stability indicator
is_stable = .true.
if (method == 1 .and. dt > dt_max_stable) then
is_stable = .false.
print *, '⚠️ WARNING: Explicit stability limit exceeded!'
print '(A,F12.6,A)', ' Calculated stable dt limit is ', dt_max_stable, ' s'
print '(A,F12.6,A)', ' Your selected time step dt is ', dt, ' s'
else
print '(A,F12.6,A)', ' Maximum stable explicit dt limit: ', dt_max_stable, ' s'
print *, ' Solver stability status: OK'
endif
print *, ''
print *, '==================================================='
print *, ' TRANSIENT TEMPERATURE FIELD PROGRESSION'
print *, '==================================================='
print '(A,9A10)', ' Time [s]', 'Node 1', 'Node 25%', 'Node 50%', 'Node 75%', 'Node N', 'E_in[J]', 'E_st[J]', 'Err [%]'
! Print initial state
call compute_energy_change(T, T_old, V_nodes, rho, k_0, cp_0, beta_k, beta_cp, &
bc_left_type, bc_left_p1, bc_left_p2, &
bc_right_type, bc_right_p1, bc_right_p2, &
0.0d0, 0.0d0, E_st_initial, E_st_current, &
E_in_total, E_out_total, balance_err)
call print_state_line(0.0d0, T, n_nodes, E_in_total, E_out_total, E_st_current, E_st_initial, balance_err)
! Time Loop
do step = 1, total_steps
T_old = T
current_time = dble(step) * dt
! Update node properties
do i = 1, n_nodes
call get_properties(T_old(i), k_0, cp_0, beta_k, beta_cp, k_nodes(i), cp_nodes(i))
end do
select case (method)
case (1) ! Explicit Euler
call solve_explicit(T, T_old, n_nodes, geom_type, dx, dt, rho, k_nodes, cp_nodes, &
bc_left_type, bc_left_p1, bc_left_p2, &
bc_right_type, bc_right_p1, bc_right_p2)
case (2, 3) ! Implicit BE or Crank-Nicolson
call solve_implicit(T, T_old, n_nodes, geom_type, dx, dt, rho, k_nodes, cp_nodes, method, &
bc_left_type, bc_left_p1, bc_left_p2, &
bc_right_type, bc_right_p1, bc_right_p2)
end select
! Compute border heat fluxes and cumulative energy inputs (for 1 m2 surface)
call get_boundary_fluxes(T, T_old, n_nodes, geom_type, dx, k_nodes(1), k_nodes(n_nodes), &
bc_left_type, bc_left_p1, bc_left_p2, &
bc_right_type, bc_right_p1, bc_right_p2, &
heat_flux_left, heat_flux_right)
E_in_total = E_in_total + heat_flux_left * dt
E_out_total = E_out_total + heat_flux_right * dt
if (mod(step, print_interval) == 0 .or. step == total_steps) then
call compute_energy_change(T, T_old, V_nodes, rho, k_0, cp_0, beta_k, beta_cp, &
bc_left_type, bc_left_p1, bc_left_p2, &
bc_right_type, bc_right_p1, bc_right_p2, &
E_in_total, E_out_total, E_st_initial, E_st_current, &
E_in_total, E_out_total, balance_err)
call print_state_line(current_time, T, n_nodes, E_in_total, E_out_total, E_st_current, E_st_initial, balance_err)
end if
end do
print *, ''
print *, '==================================================='
print *, ' FINAL NODAL TEMPERATURE DISTRIBUTION'
print *, '==================================================='
print '(A12,A18,A20)', ' Node Index', 'Coordinate [m]', 'Temperature [C]'
do i = 1, n_nodes
if (i == 1 .or. i == n_nodes .or. mod(i, max(1, n_nodes/20)) == 0) then
print '(I10,F18.6,F20.4)', i, r_nodes(i), T(i)
end if
end do
print *, '==================================================='
contains
subroutine get_properties(Temp, k0, cp0, beta_k, beta_cp, k_out, cp_out)
real(8), intent(in) :: Temp, k0, cp0, beta_k, beta_cp
real(8), intent(out) :: k_out, cp_out
k_out = k0 * (1.0d0 + beta_k * (Temp - 20.0d0))
cp_out = cp0 * (1.0d0 + beta_cp * (Temp - 20.0d0))
end subroutine get_properties
subroutine print_bc_desc(side, bc_type, p1, p2)
character(len=*), intent(in) :: side
integer, intent(in) :: bc_type
real(8), intent(in) :: p1, p2
select case (bc_type)
case (1)
print '(A,A,F10.2,A)', ' BC at ', side, p1, ' deg-C (Constant Temperature)'
case (2)
print '(A,A,F10.2,A)', ' BC at ', side, p1, ' W/m2 (Constant Heat Flux)'
case (3)
print '(A,A,F10.2,A,F10.2,A)', ' BC at ', side, p1, ' W/m2-K (Convection), T_inf = ', p2, ' deg-C'
case (4)
print '(A,A,A)', ' BC at ', side, ' Insulated (Adiabatic)'
end select
end subroutine print_bc_desc
subroutine get_boundary_fluxes(T_now, T_prev, n, geom, dx_val, k_left, k_right, &
bc_l_type, bc_l_p1, bc_l_p2, &
bc_r_type, bc_r_p1, bc_r_p2, &
q_left, q_right)
real(8), dimension(n), intent(in) :: T_now, T_prev
integer, intent(in) :: n, geom, bc_l_type, bc_r_type
real(8), intent(in) :: dx_val, k_left, k_right, bc_l_p1, bc_l_p2, bc_r_p1, bc_r_p2
real(8), intent(out) :: q_left, q_right
! Left Boundary Inward Heat Flux
select case (bc_l_type)
case (1) ! Constant T
q_left = k_left * (T_now(1) - T_now(2)) / dx_val
case (2) ! Constant q"
q_left = bc_l_p1
case (3) ! Convection
q_left = bc_l_p1 * (bc_l_p2 - T_now(1))
case (4) ! Insulated
q_left = 0.0d0
end select
! Center symmetry check for radial domains
if (geom > 1) q_left = 0.0d0
! Right Boundary Outward Heat Flux (positive escaping)
select case (bc_r_type)
case (1) ! Constant T
q_right = k_right * (T_now(n) - T_now(n-1)) / dx_val
case (2) ! Constant q" (inward flux so escaping is -inward)
q_right = -bc_r_p1
case (3) ! Convection
q_right = bc_r_p1 * (T_now(n) - bc_r_p2)
case (4) ! Insulated
q_right = 0.0d0
end select
end subroutine get_boundary_fluxes
subroutine print_state_line(t_val, T_arr, n, E_in, E_out, E_st, E_st_init, err)
real(8), intent(in) :: t_val, E_in, E_out, E_st, E_st_init, err
integer, intent(in) :: n
real(8), dimension(n), intent(in) :: T_arr
print '(F10.3,5F10.3,3ES10.2)', t_val, T_arr(1), T_arr(max(1,n/4)), T_arr(max(1,n/2)), &
T_arr(max(1,3*n/4)), T_arr(n), (E_in - E_out), (E_st - E_st_init), err
end subroutine print_state_line
subroutine compute_energy_change(T_now, T_prev, V_arr, rho_val, k0, cp0, beta_k, beta_cp, &
bc_l_type, bc_l_p1, bc_l_p2, &
bc_r_type, bc_r_p1, bc_r_p2, &
E_in, E_out, E_st_init, E_st, &
E_in_out, E_out_out, err)
real(8), dimension(:), intent(in) :: T_now, T_prev, V_arr
real(8), intent(in) :: rho_val, k0, cp0, beta_k, beta_cp
integer, intent(in) :: bc_l_type, bc_r_type
real(8), intent(in) :: bc_l_p1, bc_l_p2, bc_r_p1, bc_r_p2, E_in, E_out, E_st_init
real(8), intent(out) :: E_st, E_in_out, E_out_out, err
integer :: idx
real(8) :: temp_k, temp_cp, dE_net
E_st = 0.0d0
do idx = 1, size(T_now)
call get_properties(T_now(idx), k0, cp0, beta_k, beta_cp, temp_k, temp_cp)
E_st = E_st + rho_val * temp_cp * T_now(idx) * V_arr(idx)
end do
dE_net = E_in - E_out
E_in_out = E_in
E_out_out = E_out
err = abs(dE_net - (E_st - E_st_init)) / (abs(E_st - E_st_init) + 1.0d0) * 100.0d0
end subroutine compute_energy_change
subroutine solve_explicit(T_new, T_now, n, geom, dx_val, dt_val, rho_val, k_arr, cp_arr, &
bc_l_type, bc_l_p1, bc_l_p2, &
bc_r_type, bc_r_p1, bc_r_p2)
integer, intent(in) :: n, geom, bc_l_type, bc_r_type
real(8), dimension(n), intent(in) :: T_now, k_arr, cp_arr
real(8), dimension(n), intent(out) :: T_new
real(8), intent(in) :: dx_val, dt_val, rho_val, bc_l_p1, bc_l_p2, bc_r_p1, bc_r_p2
integer :: idx
real(8) :: r_i, r_left, r_right, k_left, k_right, Fo_local
! Interior nodes update
do idx = 2, n-1
k_left = 0.5d0 * (k_arr(idx-1) + k_arr(idx))
k_right = 0.5d0 * (k_arr(idx) + k_arr(idx+1))
select case (geom)
case (1) ! Plane
T_new(idx) = T_now(idx) + (dt_val / (rho_val * cp_arr(idx) * dx_val**2)) * &
(k_left * (T_now(idx-1) - T_now(idx)) + k_right * (T_now(idx+1) - T_now(idx)))
case (2) ! Cylinder
r_i = dble(idx - 1) * dx_val
r_left = r_i - 0.5d0 * dx_val
r_right = r_i + 0.5d0 * dx_val
T_new(idx) = T_now(idx) + (dt_val / (rho_val * cp_arr(idx) * r_i * dx_val**2)) * &
(r_left * k_left * (T_now(idx-1) - T_now(idx)) + r_right * k_right * (T_now(idx+1) - T_now(idx)))
case (3) ! Sphere
r_i = dble(idx - 1) * dx_val
r_left = r_i - 0.5d0 * dx_val
r_right = r_i + 0.5d0 * dx_val
T_new(idx) = T_now(idx) + (dt_val / (rho_val * cp_arr(idx) * r_i**2 * dx_val**2)) * &
(r_left**2 * k_left * (T_now(idx-1) - T_now(idx)) + r_right**2 * k_right * (T_now(idx+1) - T_now(idx)))
end select
end do
! Left/Inner boundary node (idx=1)
if (geom > 1) then
! Center node (insulated symmetry boundary)
Fo_local = k_arr(1) * dt_val / (rho_val * cp_arr(1) * dx_val**2)
if (geom == 2) then
T_new(1) = T_now(1) + 4.0d0 * Fo_local * (T_now(2) - T_now(1))
else ! Sphere
T_new(1) = T_now(1) + 6.0d0 * Fo_local * (T_now(2) - T_now(1))
end if
else
select case (bc_l_type)
case (1) ! Constant T
T_new(1) = bc_l_p1
case (2) ! Constant q"
Fo_local = k_arr(1) * dt_val / (rho_val * cp_arr(1) * dx_val**2)
T_new(1) = T_now(1) + 2.0d0 * Fo_local * (T_now(2) - T_now(1) + bc_l_p1 * dx_val / k_arr(1))
case (3) ! Convection
Fo_local = k_arr(1) * dt_val / (rho_val * cp_arr(1) * dx_val**2)
T_new(1) = T_now(1) + 2.0d0 * Fo_local * (T_now(2) - T_now(1) + bc_l_p1 * dx_val * (bc_l_p2 - T_now(1)) / k_arr(1))
case (4) ! Insulated
Fo_local = k_arr(1) * dt_val / (rho_val * cp_arr(1) * dx_val**2)
T_new(1) = T_now(1) + 2.0d0 * Fo_local * (T_now(2) - T_now(1))
end select
end if
! Right/Outer boundary node (idx=n)
select case (bc_r_type)
case (1) ! Constant T
T_new(n) = bc_r_p1
case (2) ! Constant q"
Fo_local = k_arr(n) * dt_val / (rho_val * cp_arr(n) * dx_val**2)
select case (geom)
case (1)
T_new(n) = T_now(n) + 2.0d0 * Fo_local * (T_now(n-1) - T_now(n) + bc_r_p1 * dx_val / k_arr(n))
case (2)
T_new(n) = T_now(n) + 2.0d0 * Fo_local * ((1.0d0 - dx_val/(2.0d0*length)) * &
(T_now(n-1) - T_now(n)) + bc_r_p1 * dx_val / k_arr(n))
case (3)
T_new(n) = T_now(n) + 2.0d0 * Fo_local * ((1.0d0 - dx_val/length)**2 * &
(T_now(n-1) - T_now(n)) + bc_r_p1 * dx_val / k_arr(n))
end select
case (3) ! Convection
Fo_local = k_arr(n) * dt_val / (rho_val * cp_arr(n) * dx_val**2)
select case (geom)
case (1)
T_new(n) = T_now(n) + 2.0d0 * Fo_local * (T_now(n-1) - T_now(n) + bc_r_p1 * dx_val * (bc_r_p2 - T_now(n)) / k_arr(n))
case (2)
T_new(n) = T_now(n) + 2.0d0 * Fo_local * ((1.0d0 - dx_val/(2.0d0*length)) * &
(T_now(n-1) - T_now(n)) + bc_r_p1 * dx_val * (bc_r_p2 - T_now(n)) / k_arr(n))
case (3)
T_new(n) = T_now(n) + 2.0d0 * Fo_local * ((1.0d0 - dx_val/length)**2 * &
(T_now(n-1) - T_now(n)) + bc_r_p1 * dx_val * (bc_r_p2 - T_now(n)) / k_arr(n))
end select
case (4) ! Insulated
Fo_local = k_arr(n) * dt_val / (rho_val * cp_arr(n) * dx_val**2)
select case (geom)
case (1)
T_new(n) = T_now(n) + 2.0d0 * Fo_local * (T_now(n-1) - T_now(n))
case (2)
T_new(n) = T_now(n) + 2.0d0 * Fo_local * (1.0d0 - dx_val/(2.0d0*length)) * (T_now(n-1) - T_now(n))
case (3)
T_new(n) = T_now(n) + 2.0d0 * Fo_local * (1.0d0 - dx_val/length)**2 * (T_now(n-1) - T_now(n))
end select
end select
end subroutine solve_explicit
subroutine solve_implicit(T_new, T_now, n, geom, dx_val, dt_val, rho_val, k_arr, cp_arr, meth_idx, &
bc_l_type, bc_l_p1, bc_l_p2, &
bc_r_type, bc_r_p1, bc_r_p2)
integer, intent(in) :: n, geom, meth_idx, bc_l_type, bc_r_type
real(8), dimension(n), intent(in) :: T_now, k_arr, cp_arr
real(8), dimension(n), intent(out) :: T_new
real(8), intent(in) :: dx_val, dt_val, rho_val, bc_l_p1, bc_l_p2, bc_r_p1, bc_r_p2
real(8), dimension(n) :: a, b, c, d
integer :: idx
real(8) :: r_i, r_left, r_right, k_left, k_right, Fo_local
real(8) :: factor, Bi_local
factor = 1.0d0
if (meth_idx == 3) factor = 0.5d0 ! Crank-Nicolson uses 0.5 weight for both implicit/explicit
! Interior nodes coefficients
do idx = 2, n-1
k_left = 0.5d0 * (k_arr(idx-1) + k_arr(idx))
k_right = 0.5d0 * (k_arr(idx) + k_arr(idx+1))
select case (geom)
case (1) ! Plane
a(idx) = -factor * (dt_val * k_left / (rho_val * cp_arr(idx) * dx_val**2))
c(idx) = -factor * (dt_val * k_right / (rho_val * cp_arr(idx) * dx_val**2))
case (2) ! Cylinder
r_i = dble(idx - 1) * dx_val
r_left = r_i - 0.5d0 * dx_val
r_right = r_i + 0.5d0 * dx_val
a(idx) = -factor * (dt_val * r_left * k_left / (rho_val * cp_arr(idx) * r_i * dx_val**2))
c(idx) = -factor * (dt_val * r_right * k_right / (rho_val * cp_arr(idx) * r_i * dx_val**2))
case (3) ! Sphere
r_i = dble(idx - 1) * dx_val
r_left = r_i - 0.5d0 * dx_val
r_right = r_i + 0.5d0 * dx_val
a(idx) = -factor * (dt_val * r_left**2 * k_left / (rho_val * cp_arr(idx) * r_i**2 * dx_val**2))
c(idx) = -factor * (dt_val * r_right**2 * k_right / (rho_val * cp_arr(idx) * r_i**2 * dx_val**2))
end select
b(idx) = 1.0d0 - a(idx) - c(idx)
if (meth_idx == 2) then
d(idx) = T_now(idx)
else ! Crank-Nicolson RHS includes explicit part
d(idx) = T_now(idx) - a(idx) * T_now(idx-1) + (a(idx) + c(idx)) * T_now(idx) - c(idx) * T_now(idx+1)
end if
end do
! Left boundary (idx=1)
a(1) = 0.0d0
if (geom > 1) then
! Symmetry at center
Fo_local = k_arr(1) * dt_val / (rho_val * cp_arr(1) * dx_val**2)
if (geom == 2) then
b(1) = 1.0d0 + 4.0d0 * factor * Fo_local
c(1) = -4.0d0 * factor * Fo_local
else ! Sphere
b(1) = 1.0d0 + 6.0d0 * factor * Fo_local
c(1) = -6.0d0 * factor * Fo_local
end if
if (meth_idx == 2) then
d(1) = T_now(1)
else
d(1) = T_now(1) - c(1) * (T_now(2) - T_now(1))
end if
else
select case (bc_l_type)
case (1) ! Constant T
b(1) = 1.0d0
c(1) = 0.0d0
d(1) = bc_l_p1
case (2) ! Constant q"
Fo_local = k_arr(1) * dt_val / (rho_val * cp_arr(1) * dx_val**2)
b(1) = 1.0d0 + 2.0d0 * factor * Fo_local
c(1) = -2.0d0 * factor * Fo_local
if (meth_idx == 2) then
d(1) = T_now(1) + 2.0d0 * Fo_local * bc_l_p1 * dx_val / k_arr(1)
else
d(1) = T_now(1) - c(1) * (T_now(2) - T_now(1)) + 2.0d0 * Fo_local * bc_l_p1 * dx_val / k_arr(1)
end if
case (3) ! Convection
Fo_local = k_arr(1) * dt_val / (rho_val * cp_arr(1) * dx_val**2)
Bi_local = bc_l_p1 * dx_val / k_arr(1)
b(1) = 1.0d0 + 2.0d0 * factor * Fo_local * (1.0d0 + Bi_local)
c(1) = -2.0d0 * factor * Fo_local
if (meth_idx == 2) then
d(1) = T_now(1) + 2.0d0 * Fo_local * Bi_local * bc_l_p2
else
d(1) = T_now(1) - c(1) * (T_now(2) - T_now(1)) + &
Fo_local * Bi_local * (2.0d0 * bc_l_p2 - T_now(1))
end if
case (4) ! Insulated
Fo_local = k_arr(1) * dt_val / (rho_val * cp_arr(1) * dx_val**2)
b(1) = 1.0d0 + 2.0d0 * factor * Fo_local
c(1) = -2.0d0 * factor * Fo_local
if (meth_idx == 2) then
d(1) = T_now(1)
else
d(1) = T_now(1) - c(1) * (T_now(2) - T_now(1))
end if
end select
end if
! Right boundary (idx=n)
c(n) = 0.0d0
select case (bc_r_type)
case (1) ! Constant T
a(n) = 0.0d0
b(n) = 1.0d0
d(n) = bc_r_p1
case (2) ! Constant q"
Fo_local = k_arr(n) * dt_val / (rho_val * cp_arr(n) * dx_val**2)
select case (geom)
case (1)
a(n) = -2.0d0 * factor * Fo_local
b(n) = 1.0d0 + 2.0d0 * factor * Fo_local
if (meth_idx == 2) then
d(n) = T_now(n) + 2.0d0 * Fo_local * bc_r_p1 * dx_val / k_arr(n)
else
d(n) = T_now(n) - a(n) * (T_now(n-1) - T_now(n)) + 2.0d0 * Fo_local * bc_r_p1 * dx_val / k_arr(n)
end if
case (2)
a(n) = -2.0d0 * factor * Fo_local * (1.0d0 - dx_val/(2.0d0*length))
b(n) = 1.0d0 - a(n)
if (meth_idx == 2) then
d(n) = T_now(n) + 2.0d0 * Fo_local * bc_r_p1 * dx_val / k_arr(n)
else
d(n) = T_now(n) - a(n) * (T_now(n-1) - T_now(n)) + 2.0d0 * Fo_local * bc_r_p1 * dx_val / k_arr(n)
end if
case (3)
a(n) = -2.0d0 * factor * Fo_local * (1.0d0 - dx_val/length)**2
b(n) = 1.0d0 - a(n)
if (meth_idx == 2) then
d(n) = T_now(n) + 2.0d0 * Fo_local * bc_r_p1 * dx_val / k_arr(n)
else
d(n) = T_now(n) - a(n) * (T_now(n-1) - T_now(n)) + 2.0d0 * Fo_local * bc_r_p1 * dx_val / k_arr(n)
end if
end select
case (3) ! Convection
Fo_local = k_arr(n) * dt_val / (rho_val * cp_arr(n) * dx_val**2)
Bi_local = bc_r_p1 * dx_val / k_arr(n)
select case (geom)
case (1)
a(n) = -2.0d0 * factor * Fo_local
b(n) = 1.0d0 + 2.0d0 * factor * Fo_local * (1.0d0 + Bi_local)
if (meth_idx == 2) then
d(n) = T_now(n) + 2.0d0 * Fo_local * Bi_local * bc_r_p2
else
d(n) = T_now(n) - a(n) * (T_now(n-1) - T_now(n)) + &
Fo_local * Bi_local * (2.0d0 * bc_r_p2 - T_now(n))
end if
case (2)
a(n) = -2.0d0 * factor * Fo_local * (1.0d0 - dx_val/(2.0d0*length))
b(n) = 1.0d0 - a(n) + 2.0d0 * factor * Fo_local * Bi_local
if (meth_idx == 2) then
d(n) = T_now(n) + 2.0d0 * Fo_local * Bi_local * bc_r_p2
else
d(n) = T_now(n) - a(n) * (T_now(n-1) - T_now(n)) + &
Fo_local * Bi_local * (2.0d0 * bc_r_p2 - T_now(n))
end if
case (3)
a(n) = -2.0d0 * factor * Fo_local * (1.0d0 - dx_val/length)**2
b(n) = 1.0d0 - a(n) + 2.0d0 * factor * Fo_local * Bi_local
if (meth_idx == 2) then
d(n) = T_now(n) + 2.0d0 * Fo_local * Bi_local * bc_r_p2
else
d(n) = T_now(n) - a(n) * (T_now(n-1) - T_now(n)) + &
Fo_local * Bi_local * (2.0d0 * bc_r_p2 - T_now(n))
end if
end select
case (4) ! Insulated
Fo_local = k_arr(n) * dt_val / (rho_val * cp_arr(n) * dx_val**2)
select case (geom)
case (1)
a(n) = -2.0d0 * factor * Fo_local
b(n) = 1.0d0 + 2.0d0 * factor * Fo_local
case (2)
a(n) = -2.0d0 * factor * Fo_local * (1.0d0 - dx_val/(2.0d0*length))
b(n) = 1.0d0 - a(n)
case (3)
a(n) = -2.0d0 * factor * Fo_local * (1.0d0 - dx_val/length)**2
b(n) = 1.0d0 - a(n)
end select
if (meth_idx == 2) then
d(n) = T_now(n)
else
d(n) = T_now(n) - a(n) * (T_now(n-1) - T_now(n))
end if
end select
! Solve tridiagonal system
call thomas_solve(a, b, c, d, T_new, n)
end subroutine solve_implicit
subroutine thomas_solve(a, b, c, d, x, n)
integer, intent(in) :: n
real(8), dimension(n), intent(in) :: a, b, c, d
real(8), dimension(n), intent(out) :: x
real(8), dimension(n) :: c_prime, d_prime
integer :: idx
real(8) :: m
c_prime(1) = c(1) / b(1)
d_prime(1) = d(1) / b(1)
do idx = 2, n - 1
m = b(idx) - a(idx) * c_prime(idx-1)
c_prime(idx) = c(idx) / m
d_prime(idx) = (d(idx) - a(idx) * d_prime(idx-1)) / m
end do
m = b(n) - a(n) * c_prime(n-1)
d_prime(n) = (d(n) - a(n) * d_prime(n-1)) / m
x(n) = d_prime(n)
do idx = n - 1, 1, -1
x(idx) = d_prime(idx) - c_prime(idx) * x(idx+1)
end do
end subroutine thomas_solve
end program transient_fdm_solver
Solver Description
Simulate transient heat conduction in 1D plates, cylinders, or spheres. Compare Explicit, Implicit, and Crank-Nicolson numerical methods in real-time.
Key Numerical Methods & Architecture
- Input Redirection: Reads parameters sequentially from standard input (`stdin`) using Fortran sequential read (`read(*,*)`), ensuring modular integration.
- Modular Design: Formulated using pure mathematical routines, separation of equations from output formatting, and precise numerical solvers (e.g. bisection, Newton-Raphson).
- Standard Compliant: Written in clean, standards-compliant Fortran 90 to ensure cross-compiler compatibility.
🛠️ Local Compilation
To test this code on your machine, compile the source code file(s) using a standard Fortran compiler (e.g., `gfortran`).
Compilation Command:
gfortran -O3 transient_fdm_solver.f90 -o fdm_solver_calc
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
fdm_solver_calc < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Geometry (1=Slab, 2=Cylinder, 3=Sphere)
1
! Length/Radius L [m]
0.05
! Number of nodes N
21
! Time step dt [s]
0.0
! Solver method (1=Explicit, 2=Implicit Euler, 3=Crank-Nicolson)
3
! Thermal conductivity k0 [W/m-K]
43.0
! Density rho [kg/m3]
7850.0
! Specific heat Cp0 [J/kg-K]
475.0
! Temp-dependency beta_k [1/K]
0.0
! Temp-dependency beta_cp [1/K]
0.0
! Left/Inner BC type (1=Temp, 2=Flux, 3=Convection, 4=Insulated)
4
! Left BC param 1
0.0
! Left BC param 2
0.0
! Right/Outer BC type
3
! Right BC param 1
150.0
! Right BC param 2
20.0
! Initial temperature type (1=Uniform, 2=Linear)
1
! Initial temp param 1
180.0
! Initial temp param 2
20.0
! Simulation time t_sim [s]
600.0
1
! Length/Radius L [m]
0.05
! Number of nodes N
21
! Time step dt [s]
0.0
! Solver method (1=Explicit, 2=Implicit Euler, 3=Crank-Nicolson)
3
! Thermal conductivity k0 [W/m-K]
43.0
! Density rho [kg/m3]
7850.0
! Specific heat Cp0 [J/kg-K]
475.0
! Temp-dependency beta_k [1/K]
0.0
! Temp-dependency beta_cp [1/K]
0.0
! Left/Inner BC type (1=Temp, 2=Flux, 3=Convection, 4=Insulated)
4
! Left BC param 1
0.0
! Left BC param 2
0.0
! Right/Outer BC type
3
! Right BC param 1
150.0
! Right BC param 2
20.0
! Initial temperature type (1=Uniform, 2=Linear)
1
! Initial temp param 1
180.0
! Initial temp param 2
20.0
! Simulation time t_sim [s]
600.0