🔬
Solver Purpose & Physical Scope
Models natural convection within enclosed cavities including vertical and horizontal rectangular gaps, concentric cylinders, and concentric spheres. Computes effective thermal conductivity ($k_{eff}$) and convection heat transfer rate using Catton and Globe-Dropkin correlations.
📂 Discipline: Convection
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 405 times
📄 Source File:
enclosure_convection.f90
📁 calcul/Convection /
enclosure_convection.f90
program enclosure_convection
implicit none
integer :: etype,ftype,i
double precision :: Lg,Hh,Ww,Di,Do,Thot,Tcold
double precision :: rho,mu,kf,Pr,cp,beta
double precision :: Tf,dT,nu_v,alp,Ra,Nu,h_eff,Q,keff,Lc,HL,g
double precision :: dTs,Ras,Nus,Qs
double precision, parameter :: pi=3.14159265358979d0
g=9.80665d0
read(*,*) etype; read(*,*) Lg; read(*,*) Hh; read(*,*) Ww
read(*,*) Di; read(*,*) Do
read(*,*) Thot; read(*,*) Tcold
read(*,*) ftype; read(*,*) rho; read(*,*) mu
read(*,*) kf; read(*,*) Pr; read(*,*) cp; read(*,*) beta
Tf=(Thot+Tcold)/2d0; dT=abs(Thot-Tcold)
if(ftype==1) then
rho=1.177d0;mu=1.85d-5;kf=0.0263d0;Pr=0.71d0;cp=1007d0;beta=1d0/(Tf+273.15d0)
elseif(ftype==2) then
rho=997d0;mu=8.9d-4;kf=0.613d0;Pr=6.13d0;cp=4180d0;beta=2.1d-4
elseif(ftype==3) then
rho=870d0;mu=0.05d0;kf=0.14d0;Pr=500d0;cp=2000d0;beta=7d-4
endif
nu_v=mu/rho; alp=kf/(rho*cp)
Ra=g*beta*dT*Lg**3/(nu_v*alp)
Nu=1d0; Q=0d0; keff=kf; Lc=Lg; HL=1d0
select case(etype)
case(1)
HL=Hh/max(Lg,1d-6)
if(HL<=2d0) then
Nu=0.18d0*(Pr*Ra/(0.2d0+Pr))**0.29d0
elseif(HL<=10d0) then
Nu=0.22d0*(Pr*Ra/(0.2d0+Pr))**0.28d0*(HL)**(-0.25d0)
else
Nu=0.42d0*Ra**0.25d0*Pr**0.012d0*(HL)**(-0.3d0)
endif
if(Nu<1d0) Nu=1d0
h_eff=Nu*kf/Lg; Q=h_eff*Hh*Ww*dT
case(2)
Nu=0.069d0*Ra**(1d0/3d0)*Pr**0.074d0
if(Nu<1d0) Nu=1d0
h_eff=Nu*kf/Lg; Q=h_eff*Hh*Ww*dT
case(3)
Nu=1d0; h_eff=kf/Lg; Q=h_eff*Hh*Ww*dT
case(4)
Lc=(Do-Di)/2d0
Ra=g*beta*dT*Lc**3/(nu_v*alp)
keff=0.386d0*kf*(Pr/(0.861d0+Pr))**0.25d0*Ra**0.25d0
if(keff<kf) keff=kf
Q=2d0*pi*keff*dT/log(Do/max(Di,1d-6))
Nu=keff/kf; h_eff=Nu*kf/Lc
case default
Lc=(Do-Di)/2d0
Ra=g*beta*dT*Lc**3/(nu_v*alp)
keff=0.74d0*kf*(Pr/(0.861d0+Pr))**0.25d0*Ra**0.25d0
if(keff<kf) keff=kf
Q=pi*Di*Do*keff*dT/Lc
Nu=keff/kf; h_eff=Nu*kf/Lc
end select
write(*,'(A)') '============================================'
write(*,'(A)') ' FREE CONVECTION IN ENCLOSURES'
write(*,'(A)') '============================================'
write(*,'(A)') ''
write(*,'(A)') '--- INPUTS ---'
write(*,'(A,I4)') ' Enclosure Type = ',etype
write(*,'(A,F10.4,A)') ' Gap L = ',Lg,' m'
write(*,'(A,F10.4,A)') ' Height H = ',Hh,' m'
write(*,'(A,F10.4,A)') ' Width W = ',Ww,' m'
write(*,'(A,F10.2,A)') ' T_hot = ',Thot,' C'
write(*,'(A,F10.2,A)') ' T_cold = ',Tcold,' C'
write(*,'(A,F10.2,A)') ' Delta T = ',dT,' C'
write(*,'(A)') ''
write(*,'(A)') '--- FLUID ---'
write(*,'(A,F10.4,A)') ' rho = ',rho,' kg/m3'
write(*,'(A,ES12.4,A)') ' mu = ',mu,' Pa.s'
write(*,'(A,F10.6,A)') ' k = ',kf,' W/mK'
write(*,'(A,F10.4)') ' Pr = ',Pr
write(*,'(A,ES12.4,A)') ' beta = ',beta,' 1/K'
write(*,'(A)') ''
write(*,'(A)') '--- RESULTS ---'
write(*,'(A,ES14.4)') ' Rayleigh Number Ra = ',Ra
write(*,'(A,F12.4)') ' Nusselt Number Nu = ',Nu
write(*,'(A,F12.4,A)') ' Effective h = ',h_eff,' W/m2K'
write(*,'(A,F12.4,A)') ' k_effective = ',keff,' W/mK'
write(*,'(A,F12.4,A)') ' Heat Transfer Q = ',Q,' W (or W/m)'
if(etype<=3) write(*,'(A,F8.2)') ' Aspect H/L = ',HL
write(*,'(A)') ''
write(*,'(A)') '--- DELTA-T SWEEP ---'
write(*,'(A)') ' dT[C] Ra Nu Q[W]'
write(*,'(A)') ' ---------------------------------------------------'
do i=1,25
dTs=1d0+(dT*3d0-1d0)*dble(i-1)/24d0
if(etype<=3) then; Ras=g*beta*dTs*Lg**3/(nu_v*alp)
else; Ras=g*beta*dTs*Lc**3/(nu_v*alp); endif
select case(etype)
case(1)
if(HL<=2d0) then; Nus=0.18d0*(Pr*Ras/(0.2d0+Pr))**0.29d0
elseif(HL<=10d0) then; Nus=0.22d0*(Pr*Ras/(0.2d0+Pr))**0.28d0*HL**(-0.25d0)
else; Nus=0.42d0*Ras**0.25d0*Pr**0.012d0*HL**(-0.3d0); endif
if(Nus<1d0) Nus=1d0; Qs=Nus*kf/Lg*Hh*Ww*dTs
case(2)
Nus=0.069d0*Ras**(1d0/3d0)*Pr**0.074d0
if(Nus<1d0) Nus=1d0; Qs=Nus*kf/Lg*Hh*Ww*dTs
case(3); Nus=1d0; Qs=kf/Lg*Hh*Ww*dTs
case(4)
keff=0.386d0*kf*(Pr/(0.861d0+Pr))**0.25d0*Ras**0.25d0
if(keff<kf) keff=kf; Nus=keff/kf
Qs=2d0*pi*keff*dTs/log(Do/max(Di,1d-6))
case default
keff=0.74d0*kf*(Pr/(0.861d0+Pr))**0.25d0*Ras**0.25d0
if(keff<kf) keff=kf; Nus=keff/kf
Qs=pi*Di*Do*keff*dTs/Lc
end select
write(*,'(2X,F8.2,2X,ES10.3,2X,F9.3,2X,F12.4)') dTs,Ras,Nus,Qs
enddo
write(*,'(A)') ''
write(*,'(A)') '--- CORRELATIONS ---'
write(*,'(A)') ' Vertical rect: Catton (1978) correlations'
write(*,'(A)') ' Horizontal hot bottom: Globe-Dropkin Nu=0.069*Ra^(1/3)*Pr^0.074'
write(*,'(A)') ' Concentric cyl: k_eff/k=0.386(Pr/(0.861+Pr))^0.25*Ra^0.25'
write(*,'(A)') ' Concentric sph: k_eff/k=0.74(Pr/(0.861+Pr))^0.25*Ra^0.25'
write(*,'(A)') ' Ref: Incropera Ch.9 Sec.9.8, Kothandaraman Ch.10'
end program enclosure_convection
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 enclosure_convection.f90 -o enclosure_convection
2. Execution with input.txt redirection:
enclosure_convection < input.txt
📄 Sample input.txt File Structure
Sample Data:
1 0.02 0.5 1.0 0.0 0.0 60.0 20.0 1 0.0
Parameter Description:
Enclosure type (1=Vert Rect, 2=Horiz hot bottom, 3=Horiz hot top, 4=Concentric Cyl, 5=Concentric Sph) Gap size L [m] Height H [m] Width W [m] Inner/First diameter Di [m] Outer/Second diameter Do [m] Hot surface temperature T_hot [C] Cold surface temperature T_cold [C] Fluid type (1=Air, 2=Water, 3=Oil) Custom thermal expansion coefficient beta [1/K] (0=auto)