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