💻 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.
Forced Tube Flow Internal Convection
Core Numerical Engine in Fortran 90 • 0 total downloads
! =========================================================================
! Source File: internal_flow_convection.f90
! =========================================================================
program internal_flow_convection
implicit none
double precision :: D, L, T_s, T_in, m_dot, rho, mu, nu, k_f, Pr, Cp
double precision :: A_c, P_w, U_m, Re_D, f
double precision :: x_fd_h, x_fd_t
double precision :: Nu_DB, Nu_Gn, Nu_used, h_avg
double precision :: T_out, T_lm, Q_total
double precision :: dP, Pump_power
character(len=30) :: flow_regime, correlation_used
double precision :: pi
integer :: i, n_points
double precision :: x, dx, T_x, T_mx
pi = 3.14159265358979d0
! Read inputs
read(*,*) D ! Tube inner diameter [m]
read(*,*) L ! Tube length [m]
read(*,*) T_s ! Surface temperature [°C] (constant wall temp)
read(*,*) T_in ! Inlet fluid temperature [°C]
read(*,*) m_dot ! Mass flow rate [kg/s]
read(*,*) rho ! Fluid density [kg/m³]
read(*,*) mu ! Dynamic viscosity [Pa·s]
read(*,*) k_f ! Thermal conductivity [W/m·K]
read(*,*) Pr ! Prandtl number
read(*,*) Cp ! Specific heat [J/kg·K]
! Computed properties
nu = mu / rho
A_c = pi * D**2 / 4.0d0
P_w = pi * D
U_m = m_dot / (rho * A_c)
Re_D = rho * U_m * D / mu
! Hydrodynamic entry length
if (Re_D < 2300.0d0) then
flow_regime = 'Laminar'
x_fd_h = 0.05d0 * Re_D * D
x_fd_t = 0.05d0 * Re_D * Pr * D
else if (Re_D < 4000.0d0) then
flow_regime = 'Transitional'
x_fd_h = 10.0d0 * D
x_fd_t = 10.0d0 * D
else
flow_regime = 'Turbulent'
x_fd_h = 10.0d0 * D
x_fd_t = 10.0d0 * D
end if
! ============================================
! FRICTION FACTOR
! ============================================
if (Re_D < 2300.0d0) then
! Laminar: f = 64/Re
f = 64.0d0 / Re_D
else
! Turbulent: Petukhov correlation
f = (0.790d0 * log(Re_D) - 1.64d0)**(-2.0d0)
end if
! ============================================
! NUSSELT NUMBER
! ============================================
if (Re_D < 2300.0d0) then
! Laminar, fully developed, constant wall temperature
Nu_DB = 3.66d0
Nu_Gn = 3.66d0
Nu_used = 3.66d0
correlation_used = 'Nu = 3.66 (const T_s)'
else
! Dittus-Boelter: Nu = 0.023 Re^0.8 Pr^n (n=0.4 for heating)
if (T_s > T_in) then
Nu_DB = 0.023d0 * Re_D**0.8d0 * Pr**0.4d0
else
Nu_DB = 0.023d0 * Re_D**0.8d0 * Pr**0.3d0
end if
! Gnielinski: Nu = (f/8)(Re-1000)Pr / [1 + 12.7(f/8)^0.5 (Pr^(2/3)-1)]
Nu_Gn = ((f/8.0d0) * (Re_D - 1000.0d0) * Pr) / &
(1.0d0 + 12.7d0 * sqrt(f/8.0d0) * (Pr**(2.0d0/3.0d0) - 1.0d0))
! Use Gnielinski (more accurate for transition)
if (Re_D >= 3000.0d0) then
Nu_used = Nu_Gn
correlation_used = 'Gnielinski'
else
Nu_used = Nu_DB
correlation_used = 'Dittus-Boelter'
end if
end if
! Average heat transfer coefficient
h_avg = Nu_used * k_f / D
! ============================================
! OUTLET TEMPERATURE (constant wall temp)
! ============================================
! T_out = T_s - (T_s - T_in) * exp(-P_w * L * h_avg / (m_dot * Cp))
T_out = T_s - (T_s - T_in) * exp(-pi * D * L * h_avg / (m_dot * Cp))
! Log-mean temperature difference
if (abs(T_s - T_out) > 1.0d-6 .and. abs(T_s - T_in) > 1.0d-6) then
T_lm = ((T_s - T_in) - (T_s - T_out)) / log((T_s - T_in) / (T_s - T_out))
else
T_lm = abs(T_s - T_in)
end if
! Total heat transfer
Q_total = m_dot * Cp * (T_out - T_in)
! Pressure drop
dP = f * (L / D) * rho * U_m**2 / 2.0d0
! Pumping power
Pump_power = dP * m_dot / rho
! ============================================
! OUTPUT
! ============================================
write(*,'(A)') '============================================================'
write(*,'(A)') ' FORCED CONVECTION - INTERNAL FLOW IN CIRCULAR TUBE'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- INPUT PARAMETERS ---'
write(*,'(A,F12.4,A)') ' Tube Diameter (D) = ', D*1000, ' mm'
write(*,'(A,F12.4,A)') ' Tube Length (L) = ', L, ' m'
write(*,'(A,F12.2,A)') ' Wall Temperature (T_s) = ', T_s, ' °C'
write(*,'(A,F12.2,A)') ' Inlet Temperature = ', T_in, ' °C'
write(*,'(A,ES12.4,A)') ' Mass Flow Rate = ', m_dot, ' kg/s'
write(*,*)
write(*,'(A)') '--- FLUID PROPERTIES ---'
write(*,'(A,F12.2,A)') ' Density = ', rho, ' kg/m³'
write(*,'(A,ES12.4,A)') ' Dynamic Viscosity = ', mu, ' Pa·s'
write(*,'(A,F12.4,A)') ' Thermal Conductivity = ', k_f, ' W/m·K'
write(*,'(A,F12.4)') ' Prandtl Number = ', Pr
write(*,'(A,F12.2,A)') ' Specific Heat = ', Cp, ' J/kg·K'
write(*,*)
write(*,'(A)') '--- FLOW ANALYSIS ---'
write(*,'(A,F12.4,A)') ' Mean Velocity = ', U_m, ' m/s'
write(*,'(A,ES12.4)') ' Reynolds Number (Re_D) = ', Re_D
write(*,'(A,A)') ' Flow Regime = ', trim(flow_regime)
write(*,'(A,ES12.4)') ' Friction Factor (f) = ', f
write(*,'(A,F12.4,A)') ' Entry Length (hydro) = ', x_fd_h, ' m'
write(*,'(A,F12.4,A)') ' Entry Length (thermal) = ', x_fd_t, ' m'
write(*,*)
write(*,'(A)') '--- HEAT TRANSFER RESULTS ---'
write(*,'(A,A)') ' Correlation Used = ', trim(correlation_used)
if (Re_D >= 2300.0d0) then
write(*,'(A,F12.4)') ' Nu (Dittus-Boelter) = ', Nu_DB
write(*,'(A,F12.4)') ' Nu (Gnielinski) = ', Nu_Gn
end if
write(*,'(A,F12.4)') ' Nusselt Number Used = ', Nu_used
write(*,'(A,F12.4,A)') ' Average h = ', h_avg, ' W/m²·K'
write(*,'(A,F12.2,A)') ' Outlet Temperature = ', T_out, ' °C'
write(*,'(A,F12.4,A)') ' LMTD = ', T_lm, ' °C'
write(*,'(A,F12.4,A)') ' Total Heat Transfer (Q) = ', Q_total, ' W'
write(*,*)
write(*,'(A)') '--- PRESSURE DROP ---'
write(*,'(A,F12.4,A)') ' Pressure Drop (ΔP) = ', dP, ' Pa'
write(*,'(A,F12.4,A)') ' Pressure Drop = ', dP/1000.0d0, ' kPa'
write(*,'(A,F12.6,A)') ' Pumping Power = ', Pump_power, ' W'
write(*,*)
! Temperature profile along tube
write(*,'(A)') '--- TEMPERATURE PROFILE ALONG TUBE ---'
write(*,'(A)') ' x [m] T_bulk [°C] T_wall [°C] ΔT [°C]'
write(*,'(A)') ' ---'
n_points = 30
dx = L / dble(n_points)
do i = 0, n_points
x = dble(i) * dx
! T_m(x) = T_s - (T_s - T_in) * exp(-P_w * x * h_avg / (m_dot * Cp))
T_mx = T_s - (T_s - T_in) * exp(-pi * D * x * h_avg / (m_dot * Cp))
write(*,'(F8.3,4X,F10.2,6X,F10.2,6X,F10.2)') x, T_mx, T_s, abs(T_s - T_mx)
end do
write(*,*)
write(*,'(A)') '--- CORRELATIONS USED ---'
write(*,'(A)') ' Laminar: Nu = 3.66 (const T_s, fully developed)'
write(*,'(A)') ' Dittus-B: Nu = 0.023 Re^0.8 Pr^n (n=0.4 heat, 0.3 cool)'
write(*,'(A)') ' Gnielinski: Nu = (f/8)(Re-1000)Pr/[1+12.7√(f/8)(Pr^(2/3)-1)]'
write(*,'(A)') ' Valid for: 0.5 < Pr < 2000, 3000 < Re < 5×10⁶'
end program internal_flow_convection
Solver Description
Models forced internal convection inside circular tubes. Calculates Reynolds number ($Re_D$), determines if flow is laminar ($Re < 2300$) or turbulent, computes friction factor ($f$), applies Dittus-Boelter or Gnielinski correlations for Nusselt number ($Nu$), and calculates outlet temperature ($T_{out}$) using log-mean temperature difference ($LMTD$) equations.
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:
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
📥 Downloads & Local Files
Preview of the required input file (input.txt):
0.025
! Tube Length ($L$) [m]
6.0
! Wall Temperature ($T_s$) [°C]
100.0
! Inlet Fluid Temperature ($T_{in}$) [°C]
20.0
! Mass Flow Rate (m) [kg/s]
0.15
! Fluid Density ($\rho$) [kg/m³]
997.0
! Fluid Viscosity ($\mu$) [Pa-s]
0.00089
! Fluid Conductivity ($k$) [W/m-K]
0.607
! Prandtl Number ($Pr$)
6.2
! Specific Heat ($C_p$) [J/kg-K]
4180.0