🔬
Solver Purpose & Physical Scope
Simulate transient heat conduction in 1D plates, cylinders, or spheres. Compare Explicit, Implicit, and Crank-Nicolson numerical methods in real-time.
📂 Discipline: Conduction
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 527 times
📄 Source File:
transient_fdm_solver.f90
📁 calcul/Conduction /
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
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 transient_fdm_solver.f90 -o fdm_solver_calc
2. Execution with input.txt redirection:
fdm_solver_calc < input.txt
📄 Sample input.txt File Structure
Sample Data:
1 0.05 21 0.0 3 43.0 7850.0 475.0 0.0 0.0 4 0.0 0.0 3 150.0 20.0 1 180.0 20.0 600.0
Parameter Description:
Geometry (1=Slab, 2=Cylinder, 3=Sphere) Length/Radius L [m] Number of nodes N Time step dt [s] Solver method (1=Explicit, 2=Implicit Euler, 3=Crank-Nicolson) Thermal conductivity k0 [W/m-K] Density rho [kg/m3] Specific heat Cp0 [J/kg-K] Temp-dependency beta_k [1/K] Temp-dependency beta_cp [1/K] Left/Inner BC type (1=Temp, 2=Flux, 3=Convection, 4=Insulated) Left BC param 1 Left BC param 2 Right/Outer BC type Right BC param 1 Right BC param 2 Initial temperature type (1=Uniform, 2=Linear) Initial temp param 1 Initial temp param 2 Simulation time t_sim [s]