🔬
Solver Purpose & Physical Scope
Analyze numerical stability limits for 1D/2D transient diffusion solvers using the FTCS scheme.
📂 Discipline: Cfd
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 248 times
📄 Source File:
heat_equation_stability.f90
📁 calcul/CFD /
heat_equation_stability.f90
program heat_equation_stability
implicit none
double precision :: alpha_in,dx,dy,dz,dt_in,k_cond,rho,cp_val,T_init,T_bc
integer :: scheme,ios,ndim,i,npts
double precision :: alpha,inv_dx2,Fo_max,rx,ry,rz,Fo,dt_use,dt_max
double precision :: t_cell,t_diff,pen_depth,margin
character(len=30) :: sch_name,stab_type,status
character(len=10) :: order_str
double precision :: tc,fc,rxc,pc,nc,Fo_ftcs
read(*,*,iostat=ios) alpha_in; if(ios/=0) stop 'ERROR: alpha'
read(*,*,iostat=ios) dx; if(ios/=0) stop 'ERROR: dx'
read(*,*,iostat=ios) dy; if(ios/=0) stop 'ERROR: dy'
read(*,*,iostat=ios) dz; if(ios/=0) stop 'ERROR: dz'
read(*,*,iostat=ios) dt_in; if(ios/=0) stop 'ERROR: dt'
read(*,*,iostat=ios) k_cond; if(ios/=0) stop 'ERROR: k'
read(*,*,iostat=ios) rho; if(ios/=0) stop 'ERROR: rho'
read(*,*,iostat=ios) cp_val; if(ios/=0) stop 'ERROR: cp'
read(*,*,iostat=ios) T_init; if(ios/=0) stop 'ERROR: T_init'
read(*,*,iostat=ios) T_bc; if(ios/=0) stop 'ERROR: T_bc'
read(*,*,iostat=ios) scheme; if(ios/=0) stop 'ERROR: scheme'
if(dx<=0d0) stop 'ERROR: dx must be positive'
if(alpha_in>0d0)then
alpha=alpha_in
else
if(k_cond<=0d0.or.rho<=0d0.or.cp_val<=0d0) stop 'ERROR: need alpha or k,rho,cp'
alpha=k_cond/(rho*cp_val)
end if
ndim=1
if(dy>0d0) ndim=2
if(dy>0d0.and.dz>0d0) ndim=3
if(dy<=0d0) dy=1d30
if(dz<=0d0) dz=1d30
inv_dx2=1d0/dx**2
if(ndim>=2) inv_dx2=inv_dx2+1d0/dy**2
if(ndim>=3) inv_dx2=inv_dx2+1d0/dz**2
select case(scheme)
case(1)
sch_name='FTCS Explicit'
stab_type='Conditionally Stable'
Fo_max=1d0/(2d0*dble(ndim))
order_str='O(dt,dx2)'
case(2)
sch_name='Implicit BE'
stab_type='Unconditionally Stable'
Fo_max=1d10
order_str='O(dt,dx2)'
case(3)
sch_name='Crank-Nicolson'
stab_type='Unconditionally Stable'
Fo_max=1d10
order_str='O(dt2,dx2)'
case(4)
sch_name='DuFort-Frankel'
stab_type='Unconditionally Stable'
Fo_max=1d10
order_str='O(dt2,dx2)'
case default
sch_name='FTCS Explicit'
stab_type='Conditionally Stable'
Fo_max=1d0/(2d0*dble(ndim))
order_str='O(dt,dx2)'
end select
dt_max=Fo_max/(alpha*inv_dx2)
if(dt_max>1d20) dt_max=1d10
if(dt_in>0d0)then
dt_use=dt_in
else
dt_use=min(dt_max*0.9d0,1d5)
end if
rx=alpha*dt_use/dx**2; ry=0d0; rz=0d0
if(ndim>=2) ry=alpha*dt_use/dy**2
if(ndim>=3) rz=alpha*dt_use/dz**2
Fo=rx+ry+rz
Fo_ftcs=1d0/(2d0*dble(ndim))
if(scheme>=2)then
status='UNCONDITIONAL'; margin=100d0
else if(Fo<=Fo_max)then
status='STABLE'; margin=(Fo_max-Fo)/max(Fo_max,1d-30)*100d0
else
status='UNSTABLE'; margin=(Fo_max-Fo)/max(Fo_max,1d-30)*100d0
end if
t_cell=dx**2/alpha; pen_depth=2d0*sqrt(alpha*dt_use); t_diff=(10d0*dx)**2/alpha
write(*,'(A)') '============================================================'
write(*,'(A)') ' HEAT EQUATION STABILITY CALCULATOR'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- INPUT CONDITIONS ----------------------------------------'
write(*,'(A,ES14.6,A)') ' alpha (diffusivity) = ', alpha, ' m2/s'
write(*,'(A,ES14.6,A)') ' k (conductivity) = ', k_cond, ' W/m.K'
write(*,'(A,F12.4,A)') ' rho = ', rho, ' kg/m3'
write(*,'(A,F12.4,A)') ' cp = ', cp_val, ' J/kg.K'
write(*,'(A,F12.2,A)') ' T_init = ', T_init, ' K'
write(*,'(A,F12.2,A)') ' T_bc = ', T_bc, ' K'
write(*,*)
write(*,'(A)') '--- GRID INFORMATION ----------------------------------------'
write(*,'(A,ES14.6,A)') ' dx = ', dx, ' m'
if(ndim>=2) write(*,'(A,ES14.6,A)') ' dy = ', dy, ' m'
if(ndim>=3) write(*,'(A,ES14.6,A)') ' dz = ', dz, ' m'
write(*,'(A,I2,A)') ' Dimensions = ', ndim, 'D'
write(*,*)
write(*,'(A)') '--- SCHEME INFORMATION --------------------------------------'
write(*,'(A,A)') ' Scheme = ', trim(sch_name)
write(*,'(A,A)') ' Stability type = ', trim(stab_type)
if(Fo_max<1d5) write(*,'(A,F12.6)') ' Fo_max (limit) = ', Fo_max
if(Fo_max>=1d5) write(*,'(A)') ' Fo_max (limit) = unlimited'
write(*,'(A,A)') ' Accuracy order = ', trim(order_str)
write(*,*)
write(*,'(A)') '--- TIME STEP -----------------------------------------------'
write(*,'(A,ES14.6,A)') ' dt (used) = ', dt_use, ' s'
write(*,'(A,ES14.6,A)') ' dt_max (FTCS stable) = ', min(dt_max,1d10), ' s'
write(*,*)
write(*,'(A)') '--- FOURIER NUMBERS -----------------------------------------'
write(*,'(A,F14.6)') ' Fo (total) = ', Fo
write(*,'(A,F14.6)') ' rx = alpha*dt/dx2 = ', rx
if(ndim>=2) write(*,'(A,F14.6)') ' ry = alpha*dt/dy2 = ', ry
if(ndim>=3) write(*,'(A,F14.6)') ' rz = alpha*dt/dz2 = ', rz
if(Fo_max<1d5) write(*,'(A,F14.6)') ' Fo_max (scheme) = ', Fo_max
write(*,*)
write(*,'(A)') '--- STABILITY STATUS ----------------------------------------'
write(*,'(A,A)') ' STATUS = ', trim(status)
write(*,'(A,F12.2,A)') ' Margin = ', margin, ' %'
write(*,*)
write(*,'(A)') '--- DIFFUSION CHARACTERISTICS -------------------------------'
write(*,'(A,ES14.6,A)') ' t_cell = dx2/alpha = ', t_cell, ' s'
write(*,'(A,ES14.6,A)') ' Penetration depth = ', pen_depth, ' m'
write(*,'(A,F14.4)') ' Pen. depth / dx = ', pen_depth/dx
write(*,'(A,ES14.6,A)') ' t_diff (10*dx domain) = ', t_diff, ' s'
if(dt_use>0d0) write(*,'(A,F14.1)') ' Steps to t_diff = ', t_diff/dt_use
write(*,*)
write(*,'(A)') '--- ALL SCHEMES COMPARISON ----------------------------------'
write(*,'(A)') ' Scheme Fo_max Stable? Order'
write(*,'(A)') ' --------------------------------------------------------'
if(Fo<=Fo_ftcs)then
write(*,'(A,F10.4,A)') ' FTCS ',Fo_ftcs,' YES O(dt,dx2)'
else
write(*,'(A,F10.4,A)') ' FTCS ',Fo_ftcs,' NO O(dt,dx2)'
end if
write(*,'(A)') ' Implicit BE unlimited YES O(dt,dx2)'
write(*,'(A)') ' Crank-Nicolson unlimited YES O(dt2,dx2)'
write(*,'(A)') ' DuFort-Frankel unlimited YES* O(dt2,dx2)'
write(*,*)
write(*,'(A)') '--- PROFILE vs dt ------------------------------------------'
write(*,'(A)') ' dt Fo rx Status pen_depth N_steps'
write(*,'(A)') ' --------------------------------------------------------------------------'
npts=40
do i=1,npts
if(dt_max<1d5)then
tc=dt_max/20d0+(dt_max*4d0-dt_max/20d0)*dble(i-1)/dble(npts-1)
else
tc=dt_use/20d0+(dt_use*40d0-dt_use/20d0)*dble(i-1)/dble(npts-1)
end if
if(tc<=0d0) tc=1d-15
fc=alpha*tc*inv_dx2; rxc=alpha*tc/dx**2; pc=2d0*sqrt(alpha*tc)
nc=t_diff/max(tc,1d-30)
if(fc<=Fo_ftcs)then
write(*,'(ES12.4,2X,F10.6,2X,F10.6,2X,A8,2X,ES12.4,2X,ES12.4)') tc,fc,rxc,'STABLE ',pc,nc
else
write(*,'(ES12.4,2X,F10.6,2X,F10.6,2X,A8,2X,ES12.4,2X,ES12.4)') tc,fc,rxc,'UNSTABLE',pc,nc
end if
end do
write(*,*)
write(*,'(A)') '--- EQUATIONS USED ------------------------------------------'
write(*,'(A)') ' Fo = alpha*dt/dx^2 (Fourier number, Eq. 2.80)'
write(*,'(A)') ' FTCS: Fo <= 1/(2*ndim)'
write(*,'(A)') ' 1D: Fo<=0.5, 2D: Fo<=0.25, 3D: Fo<=1/6'
write(*,'(A)') ' dt_max = Fo_max / (alpha * sum(1/dxi^2))'
write(*,'(A)') ' Penetration depth = 2*sqrt(alpha*t)'
write(*,'(A)') '============================================================'
end program heat_equation_stability
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 heat_equation_stability.f90 -o heat_equation_stability
2. Execution with input.txt redirection:
heat_equation_stability < input.txt
📄 Sample input.txt File Structure
Sample Data:
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 8
Parameter Description:
a_i\nx [m]\ny [m]\nz [m]\nt [s]\nk_i\nr_i\ncp [J/kgK]\nti_i\ntb_i\nsc_i