💻 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

internal_flow_convection.f90
! =========================================================================
! 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:

gfortran -ffree-line-length-none internal_flow_convection.f90 -o internal_flow_conv

Execution Command:

Execute the program by feeding the sample input file into the program using stdin redirection:

internal_flow_conv < input.txt

📥 Downloads & Local Files

Preview of the required input file (input.txt):

! Tube Inner Diameter ($D$) [m]
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