🔬
Solver Purpose & Physical Scope
For forced convection inside tubes, flow is laminar for $Re_D 10,000$. The Nusselt number is determined by Dittus-Boelter correlation for turbulent flow:
📂 Discipline: Convection
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 311 times
📄 Source File:
internal_flow_convection.f90
📁 calcul/Convection /
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
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 internal_flow_convection.f90 -o internal_flow_conv
2. Execution with input.txt redirection:
internal_flow_conv < input.txt
📄 Sample input.txt File Structure
Sample Data:
0.025 5.0 100.0 20.0 0.1 998.0 1.002e-3 0.598 7.01 4182.0
Parameter Description:
Tube Diameter D [m] Tube Length L [m] Surface Temperature Ts [°C] Inlet Temperature Tin [°C] Mass Flow Rate mdot [kg/s] Fluid Density [kg/m3] Dynamic Viscosity [Pa-s] Thermal Conductivity k [W/m-K] Prandtl Number Pr Specific Heat Cp [J/kg-K]