🔬
Solver Purpose & Physical Scope
Natural convection heat transfer is driven by buoyancy forces. The strength of the natural convection is measured by the Rayleigh number ($Ra_L$), which dictates the Nusselt number ($Nu_L$) via empirical relations (like the Churchill and Chu correlation for vertical plates):
📂 Discipline: Convection
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 291 times
📄 Source File:
natural_convection.f90
📁 calcul/Convection /
natural_convection.f90
program natural_convection
implicit none
! Inputs
integer :: geom_type
double precision :: dim1, dim2
double precision :: T_s, T_inf
double precision :: rho, mu, k_f, Pr, Cp, beta
! Constants
double precision, parameter :: g = 9.80665d0
double precision, parameter :: PI = 3.141592653589793d0
! Physical properties
double precision :: nu, L_c, A_s, T_f, dT
double precision :: Gr, Ra, Nu_avg, h_avg, Q
! Local profile variables
integer :: i, n_points
double precision :: x, dx, local_Ra, local_Gr, local_Nu, local_h, local_delta
character(len=50) :: geom_name
integer :: iostat_val
! Read inputs
read(*,*,iostat=iostat_val) geom_type
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid geometry type input.'
stop
end if
read(*,*,iostat=iostat_val) dim1
read(*,*,iostat=iostat_val) dim2
read(*,*,iostat=iostat_val) T_s
read(*,*,iostat=iostat_val) T_inf
read(*,*,iostat=iostat_val) rho
read(*,*,iostat=iostat_val) mu
read(*,*,iostat=iostat_val) k_f
read(*,*,iostat=iostat_val) Pr
read(*,*,iostat=iostat_val) Cp
read(*,*,iostat=iostat_val) beta
if (iostat_val /= 0) then
write(*,*) 'ERROR: Failed to read all convection parameters.'
stop
end if
! Basic validations
if (dim1 <= 0.0d0) then
write(*,*) 'ERROR: Dimension 1 must be positive.'
stop
end if
if (rho <= 0.0d0 .or. mu <= 0.0d0 .or. k_f <= 0.0d0 .or. Pr <= 0.0d0 .or. Cp <= 0.0d0) then
write(*,*) 'ERROR: Fluid properties must be positive.'
stop
end if
! Compute film temperature and kinematic viscosity
T_f = (T_s + T_inf) / 2.0d0
nu = mu / rho
dT = abs(T_s - T_inf)
! Setup Geometry parameters
select case (geom_type)
case (1)
geom_name = "Vertical Plate"
L_c = dim1
A_s = dim1 * dim2
case (2)
geom_name = "Horizontal Plate (Upper Hot / Lower Cold)"
if (dim2 <= 0.0d0) dim2 = dim1 ! Default to square if W is zero
L_c = (dim1 * dim2) / (2.0d0 * (dim1 + dim2))
A_s = dim1 * dim2
case (3)
geom_name = "Horizontal Plate (Lower Hot / Upper Cold)"
if (dim2 <= 0.0d0) dim2 = dim1
L_c = (dim1 * dim2) / (2.0d0 * (dim1 + dim2))
A_s = dim1 * dim2
case (4)
geom_name = "Horizontal Cylinder"
if (dim2 <= 0.0d0) dim2 = 1.0d0 ! Default to 1m length
L_c = dim1
A_s = PI * dim1 * dim2
case (5)
geom_name = "Sphere"
L_c = dim1
A_s = PI * dim1**2
case default
write(*,*) 'ERROR: Invalid geometry type selected.'
stop
end select
! Grashof and Rayleigh numbers
Gr = (g * beta * dT * L_c**3) / (nu**2)
Ra = Gr * Pr
! Nusselt Number Calculations
select case (geom_type)
case (1)
! Churchill and Chu (1975) for Vertical Plate (All Ra)
Nu_avg = (0.825d0 + (0.387d0 * Ra**(1.0d0/6.0d0)) / &
((1.0d0 + (0.492d0/Pr)**(9.0d0/16.0d0))**(8.0d0/27.0d0)))**2
case (2)
! Horizontal Plate (Upper Hot / Lower Cold)
if (Ra < 1.0d4) then
Nu_avg = 0.54d0 * Ra**0.25d0 ! Extrapolated/lower-laminar
else if (Ra <= 1.0d7) then
Nu_avg = 0.54d0 * Ra**0.25d0
else if (Ra <= 1.0d11) then
Nu_avg = 0.15d0 * Ra**(1.0d0/3.0d0)
else
Nu_avg = 0.15d0 * Ra**(1.0d0/3.0d0) ! Extrapolated/upper-turbulent
end if
case (3)
! Horizontal Plate (Lower Hot / Upper Cold)
! Valid for 10^5 to 10^11
Nu_avg = 0.27d0 * Ra**0.25d0
case (4)
! Churchill and Chu (1975) for Horizontal Cylinder (Ra <= 10^12)
Nu_avg = (0.60d0 + (0.387d0 * Ra**(1.0d0/6.0d0)) / &
((1.0d0 + (0.559d0/Pr)**(9.0d0/16.0d0))**(8.0d0/27.0d0)))**2
case (5)
! Churchill (1983) for Sphere (Ra <= 10^11 and Pr >= 0.7)
Nu_avg = 2.0d0 + (0.589d0 * Ra**0.25d0) / &
((1.0d0 + (0.469d0/Pr)**(9.0d0/16.0d0))**(4.0d0/9.0d0))
end select
! Average heat transfer coefficient
h_avg = Nu_avg * k_f / L_c
! Heat transfer rate
Q = h_avg * A_s * dT
! Output results
write(*,'(A)') '============================================================'
write(*,'(A)') ' NATURAL CONVECTION HEAT TRANSFER ENGINE'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A,A)') ' Geometry Name = ', trim(geom_name)
write(*,'(A,I2)') ' Geometry Type Code = ', geom_type
write(*,'(A,F12.4,A)') ' Characteristic Length = ', L_c, ' m'
write(*,'(A,F12.4,A)') ' Surface Area (As) = ', A_s, ' m2'
write(*,'(A,F12.2,A)') ' Surface Temperature = ', T_s, ' deg-C'
write(*,'(A,F12.2,A)') ' Ambient Temperature = ', T_inf, ' deg-C'
write(*,'(A,F12.2,A)') ' Film Temperature = ', T_f, ' deg-C'
write(*,*)
write(*,'(A)') '--- FLUID PROPERTIES ----------------------------------------'
write(*,'(A,ES12.4,A)') ' Density (rho) = ', rho, ' kg/m3'
write(*,'(A,ES12.4,A)') ' Dynamic Viscosity (mu) = ', mu, ' Pa.s'
write(*,'(A,ES12.4,A)') ' Kinematic Viscosity (nu) = ', nu, ' m2/s'
write(*,'(A,F12.4,A)') ' Thermal Conductivity = ', k_f, ' W/m.K'
write(*,'(A,F12.4)') ' Prandtl Number (Pr) = ', Pr
write(*,'(A,F12.2,A)') ' Specific Heat (Cp) = ', Cp, ' J/kg.K'
write(*,'(A,ES12.4,A)') ' Expansion Coeff (beta) = ', beta, ' 1/K'
write(*,*)
write(*,'(A)') '--- CONVECTION RESULTS --------------------------------------'
write(*,'(A,ES12.4)') ' Grashof Number (Gr) = ', Gr
write(*,'(A,ES12.4)') ' Rayleigh Number (Ra) = ', Ra
write(*,'(A,F12.4)') ' Average Nusselt Number = ', Nu_avg
write(*,'(A,F12.4,A)') ' Average h = ', h_avg, ' W/m2.K'
write(*,'(A,F12.4,A)') ' Heat Transfer Rate (Q) = ', Q, ' W'
write(*,*)
! ============================================
! LOCAL PROFILE DATA (along coordinate x)
! ============================================
write(*,'(A)') '--- LOCAL PROFILE ALONG SURFACE ----------------------------'
write(*,'(A)') ' x [m] Ra_x Nu_x h_x [W/m2.K] d [mm]'
write(*,'(A)') ' ----------------------------------------------------------'
n_points = 40
dx = L_c / dble(n_points)
do i = 1, n_points
x = dble(i) * dx
local_Gr = (g * beta * dT * x**3) / (nu**2)
local_Ra = local_Gr * Pr
if (dT <= 1.0d-10 .or. local_Gr <= 1.0d-10) then
local_Nu = 0.0d0
local_delta = 0.0d0
else
if (local_Ra < 1.0d9) then
! Laminar boundary layer correlations
local_Nu = 0.508d0 * (Pr / (0.952d0 + Pr))**0.25d0 * local_Ra**0.25d0
local_delta = 3.93d0 * x * ((0.952d0 + Pr) / (Pr**2 * local_Gr))**0.25d0
else
! Turbulent boundary layer approximations
local_Nu = 0.0292d0 * local_Ra**0.4d0
local_delta = 0.15d0 * x * local_Ra**(-0.1d0)
end if
end if
if (x > 0.0d0) then
local_h = local_Nu * k_f / x
else
local_h = 0.0d0
end if
write(*,'(F8.4,2X,ES12.4,2X,F10.2,4X,F10.4,4X,F10.4)') &
x, local_Ra, local_Nu, local_h, local_delta*1000.0d0
end do
write(*,*)
write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
select case (geom_type)
case (1)
write(*,'(A)') ' Vertical Plate: Churchill and Chu (1975) Correlation'
case (2)
write(*,'(A)') ' Horizontal Plate (Upper Hot): Nu = 0.54 Ra^0.25 (lam), Nu = 0.15 Ra^(1/3) (turb)'
case (3)
write(*,'(A)') ' Horizontal Plate (Lower Hot): Nu = 0.27 Ra^0.25'
case (4)
write(*,'(A)') ' Horizontal Cylinder: Churchill and Chu (1975) Correlation'
case (5)
write(*,'(A)') ' Sphere: Churchill (1983) Correlation'
end select
write(*,'(A)') ' Boundary Layer: Laminar similarity solutions for natural convection.'
end program natural_convection
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 natural_convection.f90 -o natural_convection
2. Execution with input.txt redirection:
natural_convection < input.txt
📄 Sample input.txt File Structure
Sample Data:
1 1.0 0.5 80.0 20.0 1.161 1.9e-5 0.0263 0.71 1007.0 0.003
Parameter Description:
Geometry (1=Vertical Plate, 2=Horizontal Plate, 3=Vertical Cylinder, 4=Horizontal Cylinder) Characteristic dimension (Length/Height) [m] Secondary dimension (Width/Diameter) [m] Surface Temperature Ts [°C] Ambient Temperature Tinf [°C] Fluid Density [kg/m3] Dynamic Viscosity [Pa-s] Thermal Conductivity k [W/m-K] Prandtl Number Pr Specific Heat Cp [J/kg-K] Thermal Expansion Coefficient beta [1/K]