💻 Fortran Source Code Library
Run calculations locally on your own machine. View code structure, read technical explanations, and download compilation packages including sample input files.
Natural Convection Buoyancy Solver
Core Numerical Engine in Fortran 90 • 0 total downloads
natural_convection.f90
! =========================================================================
! Source File: 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
Solver Description
Analyzes natural convection heat transfer. Solves Grashof ($Gr$) and Rayleigh ($Ra$) numbers, selects appropriate Nusselt correlations (such as Churchill-Chu for vertical plates), and calculates convective coefficient ($h$) and buoyancy-induced heat loss ($Q$).
Key Numerical Methods & Architecture
- Input Redirection: Reads parameters sequentially from standard input (`stdin`) using Fortran sequential read (`read(*,*)`), ensuring modular integration.
- Modular Design: Formulated using pure mathematical routines, separation of equations from output formatting, and precise numerical solvers (e.g. bisection, Newton-Raphson).
- Standard Compliant: Written in clean, standards-compliant Fortran 90 to ensure cross-compiler compatibility.
🛠️ Local Compilation
To test this code on your machine, compile the source code file(s) using a standard Fortran compiler (e.g., `gfortran`).
Compilation Command:
gfortran -ffree-line-length-none natural_convection.f90 -o natural_convection
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
natural_convection < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Geometry Code (1=Vertical Plate, 2=Horiz Upper, 3=Horiz Lower, 4=Horiz Cylinder, 5=Sphere)
1
! Characteristic Dimension ($L$ or $D$) [m]
0.3
! Surface Temp ($T_s$) [°C]
60.0
! Ambient Temp ($T_\infty$) [°C]
20.0
! Fluid Density ($\rho$) [kg/m³]
1.164
! Fluid Viscosity ($\mu$) [Pa-s]
1.9e-5
! Fluid Conductivity ($k$) [W/m-K]
0.026
! Prandtl Number ($Pr$)
0.707
! Specific Heat ($C_p$) [J/kg-K]
1007.0
! Volumetric Expansion Coefficient ($\beta$) [1/K]
0.0031
1
! Characteristic Dimension ($L$ or $D$) [m]
0.3
! Surface Temp ($T_s$) [°C]
60.0
! Ambient Temp ($T_\infty$) [°C]
20.0
! Fluid Density ($\rho$) [kg/m³]
1.164
! Fluid Viscosity ($\mu$) [Pa-s]
1.9e-5
! Fluid Conductivity ($k$) [W/m-K]
0.026
! Prandtl Number ($Pr$)
0.707
! Specific Heat ($C_p$) [J/kg-K]
1007.0
! Volumetric Expansion Coefficient ($\beta$) [1/K]
0.0031