🔬
Solver Purpose & Physical Scope
Calculates combined natural and forced convection on vertical/horizontal surfaces. Uses Richardson number $Ri = Gr / Re^2$ to identify heat transfer regimes: pure forced ($Ri < 0.1$), mixed ($0.1 \le Ri \le 10$), or pure natural ($Ri > 10$). Computes combined Nusselt number $Nu^3 = Nu_{forced}^3 \pm Nu_{natural}^3$ based on alignment.
📂 Discipline: Convection
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 305 times
📄 Source File:
combined_convection.f90
📁 calcul/Convection /
combined_convection.f90
program combined_convection
implicit none
double precision :: D, V, Ti, Ts, Lc, rho, mu, kf, Pr, cp, beta
integer :: or_type, fl_type
double precision :: Tf, dT, nu, alpha, Re, Gr, Ra, Ri, Nuf, Nun, Nuc, h, g
double precision :: temp_val
integer :: iostat_val, i
double precision :: Vs, Res, Grs, Ris, Nufs, Nuns, Nucs, hs
g = 9.80665d0
! Read inputs sequentially
read(*,*,iostat=iostat_val) D
read(*,*,iostat=iostat_val) V
read(*,*,iostat=iostat_val) Ti
read(*,*,iostat=iostat_val) Ts
read(*,*,iostat=iostat_val) Lc
read(*,*,iostat=iostat_val) or_type
read(*,*,iostat=iostat_val) fl_type
read(*,*,iostat=iostat_val) rho
read(*,*,iostat=iostat_val) mu
read(*,*,iostat=iostat_val) kf
read(*,*,iostat=iostat_val) Pr
read(*,*,iostat=iostat_val) cp
read(*,*,iostat=iostat_val) temp_val ! read the 13th value which is beta from php, or custom beta
if (iostat_val /= 0) then
write(*,*) 'ERROR: Failed to read all input values.'
stop
end if
Tf = (Ti + Ts) / 2.0d0
dT = abs(Ts - Ti)
! Select defaults based on fluid type
if (fl_type == 1) then
! Air
rho = 1.177d0
mu = 1.85d-5
kf = 0.0263d0
Pr = 0.71d0
cp = 1007d0
beta = 1.0d0 / (Tf + 273.15d0)
else if (fl_type == 2) then
! Water
rho = 997d0
mu = 8.9d-4
kf = 0.613d0
Pr = 6.13d0
cp = 4180d0
beta = 2.1d-4
else
! Engine Oil
rho = 870d0
mu = 0.05d0
kf = 0.14d0
Pr = 500d0
cp = 2000d0
beta = 7.0d-4
end if
! Override if user provided custom beta
if (temp_val > 0.0d0) beta = temp_val
nu = mu / rho
alpha = kf / (rho * cp)
Re = rho * V * D / mu
Gr = g * beta * dT * Lc**3 / nu**2
Ra = Gr * Pr
if (Re > 1.0d-10) then
Ri = Gr / Re**2
else
Ri = 9999.0d0
end if
! Forced Convection Nusselt Nuf (Flat Plate model)
if (Re < 5.0d5) then
Nuf = 0.664d0 * Re**0.5d0 * Pr**(1.0d0/3.0d0)
else
Nuf = (0.037d0 * Re**0.8d0 - 871.0d0) * Pr**(1.0d0/3.0d0)
end if
if (Nuf < 0.0d0) Nuf = 0.0d0
! Natural Convection Nusselt Nun (Vertical Plate model)
Nun = (0.825d0 + (0.387d0 * Ra**(1.0d0/6.0d0)) / &
((1.0d0 + (0.492d0/Pr)**(9.0d0/16.0d0))**(8.0d0/27.0d0)))**2
! Combined Nusselt Nuc
if (or_type == 1) then
! Assisting
Nuc = (Nuf**3 + Nun**3)**(1.0d0/3.0d0)
else if (or_type == 2) then
! Opposing
Nuc = abs(Nuf**3 - Nun**3)**(1.0d0/3.0d0)
else
! Transverse / Cross Flow
Nuc = (Nuf**3 + Nun**3)**(1.0d0/3.0d0)
end if
h = Nuc * kf / Lc
write(*,'(A)') '============================================'
write(*,'(A)') ' COMBINED FREE & FORCED CONVECTION ENGINE'
write(*,'(A)') '============================================'
write(*,*)
write(*,'(A,ES14.4)') ' Reynolds Re = ', Re
write(*,'(A,ES14.4)') ' Grashof Number = ', Gr
write(*,'(A,ES14.4)') ' Richardson Number Ri = ', Ri
write(*,'(A,F12.4)') ' Nu Forced = ', Nuf
write(*,'(A,F12.4)') ' Nu Natural = ', Nun
write(*,'(A,F12.4)') ' Nu Combined = ', Nuc
write(*,'(A,F12.4,A)') ' Coeff h = ', h, ' W/m2K'
write(*,*)
write(*,'(A)') '--- VELOCITY SWEEP ---'
write(*,'(A)') ' V[m/s] Re Gr Ri Nu_f Nu_n Nu_c h[W/m2K]'
write(*,'(A)') ' ------------------------------------------------------------------------------------------------'
do i=1,25
Vs = 0.01d0 + (max(V,0.2d0)*3.0d0 - 0.01d0)*dble(i-1)/24.0d0
Res = rho * Vs * D / mu
Grs = g * beta * dT * Lc**3 / nu**2
if (Res > 1.0d-10) then
Ris = Grs / Res**2
else
Ris = 9999.0d0
end if
if (Res < 5.0d5) then
Nufs = 0.664d0 * Res**0.5d0 * Pr**(1.0d0/3.0d0)
else
Nufs = (0.037d0 * Res**0.8d0 - 871.0d0) * Pr**(1.0d0/3.0d0)
end if
if (Nufs < 0.0d0) Nufs = 0.0d0
Nuns = (0.825d0 + (0.387d0 * Ra**(1.0d0/6.0d0)) / &
((1.0d0 + (0.492d0/Pr)**(9.0d0/16.0d0))**(8.0d0/27.0d0)))**2
if (or_type == 1) then
Nucs = (Nufs**3 + Nuns**3)**(1.0d0/3.0d0)
else if (or_type == 2) then
Nucs = abs(Nufs**3 - Nuns**3)**(1.0d0/3.0d0)
else
Nucs = (Nufs**3 + Nuns**3)**(1.0d0/3.0d0)
end if
hs = Nucs * kf / Lc
write(*,'(2X,F8.3,2X,ES10.3,2X,ES10.3,2X,ES10.3,2X,F9.2,2X,F9.2,2X,F9.2,2X,F10.3)') Vs, Res, Grs, Ris, Nufs, Nuns, Nucs, hs
end do
write(*,*)
end program combined_convection
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 combined_convection.f90 -o combined_convection
2. Execution with input.txt redirection:
combined_convection < input.txt
📄 Sample input.txt File Structure
Sample Data:
0.025 0.5 25.0 100.0 0.5 1 1 0.0 0.0 0.0 0.0 0.0 0.0
Parameter Description:
Characteristic diameter D [m] Free stream velocity V [m/s] Free stream temperature Tinf [C] Surface temperature Ts [C] Characteristic length L [m] Orientation (1=Assisting, 2=Opposing, 3=Transverse) Fluid type (1=Air, 2=Water, 3=Oil) Custom density [kg/m3] (0=auto) Custom viscosity [Pa-s] (0=auto) Custom thermal conductivity [W/m-K] (0=auto) Custom Prandtl number (0=auto) Custom specific heat [J/kg-K] (0=auto) Custom thermal expansion coefficient beta [1/K] (0=auto)