! ============================================================================ ! ThermoFluidCalc — Fanno Flow Solver ! Reference: Shapiro, Dynamics and Thermodynamics of Compressible Fluid Flow ! ============================================================================ program fanno_flow implicit none ! Input variables double precision :: M1 ! Inlet Mach number double precision :: f ! Darcy friction factor double precision :: D ! Pipe diameter [mm] double precision :: L ! Pipe length [m] double precision :: gamma ! Specific heat ratio (default 1.4) double precision :: T0 ! Inlet Stagnation Temp [K] (0 = skip) double precision :: P0 ! Inlet Stagnation Pressure [kPa] (0 = skip) ! Constants double precision, parameter :: R_air = 287.0d0 ! Gas constant [J/kg-K] ! Intermediate and output variables double precision :: D_m ! Diameter in meters double precision :: fLstar_over_D_inlet, fL_over_D double precision :: Lstar ! Choking length [m] double precision :: Lstar_max_shock double precision :: M2 ! Exit Mach number double precision :: T2_T1, P2_P1, rho2_rho1, P02_P01 double precision :: T1, P1, rho1, V1, T2, P2, rho2, V2 double precision :: P02, T02 ! exit stagnation properties double precision :: T_Tstar1, P_Pstar1, rho_rhostar1, V_Vstar1, P0_P0star1 double precision :: T_Tstar2, P_Pstar2, rho_rhostar2, V_Vstar2, P0_P0star2 double precision :: delta_s_over_R double precision :: M_shock1, M_shock2, x_shock integer :: status ! 1=Subsonic, 2=Supersonic, 3=Choked Subsonic, 4=Choked Shock in Duct, 5=Choked Shock Pushed Upstream ! Read input from stdin read(*,*) M1 read(*,*) f read(*,*) D read(*,*) L read(*,*) gamma read(*,*) T0 read(*,*) P0 ! Set defaults and bounds if (gamma <= 1.0d0) gamma = 1.4d0 if (f <= 0.0d0) f = 0.02d0 if (D <= 0.0d0) D = 50.0d0 if (L < 0.0d0) L = 0.0d0 if (M1 <= 0.0d0) M1 = 0.5d0 D_m = D / 1000.0d0 fLstar_over_D_inlet = fanno_param(M1, gamma) Lstar = fLstar_over_D_inlet * D_m / f fL_over_D = f * L / D_m ! Initialize shock variables M_shock1 = 0.0d0 M_shock2 = 0.0d0 x_shock = 0.0d0 ! ── DETERMINE FLOW STATUS AND SOLVE EXIT MACH ──────────── if (M1 < 1.0d0) then ! SUBSONIC INLET FLOW if (L <= Lstar) then status = 1 ! Subsonic unchoked M2 = solve_mach_from_fanno(fLstar_over_D_inlet - fL_over_D, gamma, M1, 1.0d0) else status = 3 ! Choked Subsonic (Inlet Mach must decrease) ! Solve new inlet Mach such that Lstar = L M1 = solve_mach_from_fanno(fL_over_D, gamma, 0.001d0, 1.0d0) fLstar_over_D_inlet = fanno_param(M1, gamma) Lstar = L M2 = 1.0d0 end if else ! SUPERSONIC INLET FLOW if (L <= Lstar) then status = 2 ! Supersonic unchoked M2 = solve_mach_from_fanno(fLstar_over_D_inlet - fL_over_D, gamma, 1.0d0, M1) else ! Check if shock can fit inside the duct ! Shock at inlet means flow becomes subsonic at inlet. ! Post-shock Mach of M1: M_shock2 = sqrt((2.0d0 + (gamma - 1.0d0) * M1**2) / (2.0d0 * gamma * M1**2 - (gamma - 1.0d0))) Lstar_max_shock = Lstar + fanno_param(M_shock2, gamma) * D_m / f if (L <= Lstar_max_shock) then status = 4 ! Shock inside the duct M2 = 1.0d0 ! Solve for shock Mach number M_shock1 in range [1.0, M1] call solve_shock_position(M1, f, D_m, L, gamma, M_shock1, M_shock2, x_shock) else status = 5 ! Shock pushed back upstream ! Flow is completely subsonic in the duct, choked at exit. ! Inlet Mach is reduced to subsonic value M1 = solve_mach_from_fanno(fL_over_D, gamma, 0.001d0, 1.0d0) fLstar_over_D_inlet = fanno_param(M1, gamma) Lstar = L M2 = 1.0d0 end if end if end if ! ── COMPUTE FANNO RATIOS AT INLET AND EXIT ─────────────── T_Tstar1 = (gamma + 1.0d0) / (2.0d0 + (gamma - 1.0d0) * M1**2) P_Pstar1 = (1.0d0 / M1) * sqrt(T_Tstar1) rho_rhostar1 = (1.0d0 / M1) * sqrt((2.0d0 + (gamma - 1.0d0) * M1**2) / (gamma + 1.0d0)) V_Vstar1 = 1.0d0 / rho_rhostar1 P0_P0star1 = (1.0d0 / M1) * ((2.0d0 + (gamma - 1.0d0) * M1**2) / (gamma + 1.0d0))**((gamma + 1.0d0) / (2.0d0 * (gamma - 1.0d0))) T_Tstar2 = (gamma + 1.0d0) / (2.0d0 + (gamma - 1.0d0) * M2**2) P_Pstar2 = (1.0d0 / M2) * sqrt(T_Tstar2) rho_rhostar2 = (1.0d0 / M2) * sqrt((2.0d0 + (gamma - 1.0d0) * M2**2) / (gamma + 1.0d0)) V_Vstar2 = 1.0d0 / rho_rhostar2 P0_P0star2 = (1.0d0 / M2) * ((2.0d0 + (gamma - 1.0d0) * M2**2) / (gamma + 1.0d0))**((gamma + 1.0d0) / (2.0d0 * (gamma - 1.0d0))) ! Property changes across duct T2_T1 = T_Tstar2 / T_Tstar1 P2_P1 = P_Pstar2 / P_Pstar1 rho2_rho1 = rho_rhostar2 / rho_rhostar1 P02_P01 = P0_P0star2 / P0_P0star1 delta_s_over_R = log(P0_P0star1 / P0_P0star2) ! ── COMPUTE OPTIONAL ACTUAL PROPERTIES ─────────────────── T1 = 0.0d0; P1 = 0.0d0; rho1 = 0.0d0; V1 = 0.0d0 T2 = 0.0d0; P2 = 0.0d0; rho2 = 0.0d0; V2 = 0.0d0 P02 = 0.0d0; T02 = 0.0d0 if (T0 > 0.0d0) then T1 = T0 / (1.0d0 + (gamma - 1.0d0) / 2.0d0 * M1**2) V1 = M1 * sqrt(gamma * R_air * T1) T2 = T1 * T2_T1 V2 = M2 * sqrt(gamma * R_air * T2) T02 = T0 ! Stagnation temperature is constant in adiabatic Fanno flow end if if (P0 > 0.0d0) then P1 = P0 / (1.0d0 + (gamma - 1.0d0) / 2.0d0 * M1**2)**(gamma / (gamma - 1.0d0)) P2 = P1 * P2_P1 P02 = P0 * P02_P01 if (T0 > 0.0d0) then rho1 = P1 / (R_air * T1 / 1000.0d0) ! P1 is in kPa, convert R to kJ/kg-K rho2 = P2 / (R_air * T2 / 1000.0d0) end if end if ! ── OUTPUT RESULTS IN KEY-VALUE FORMAT ─────────────────── write(*, '(A, I2)') "Status Code = ", status select case (status) case (1) write(*, '(A)') "Status = Subsonic Flow" case (2) write(*, '(A)') "Status = Supersonic Flow" case (3) write(*, '(A)') "Status = Choked Subsonic Flow (Inlet Mach Reduced)" case (4) write(*, '(A)') "Status = Choked Supersonic Flow (Normal Shock in Duct)" case (5) write(*, '(A)') "Status = Choked Supersonic Flow (Shock Pushed Upstream)" end select write(*, '(A, F14.6)') "Inlet Mach (M1) = ", M1 write(*, '(A, F14.6)') "Exit Mach (M2) = ", M2 write(*, '(A, F14.6)') "Friction Factor (f) = ", f write(*, '(A, F14.6)') "Pipe Diameter (D) = ", D write(*, '(A, F14.6)') "Pipe Length (L) = ", L write(*, '(A, F14.6)') "Choking Length (Lstar) = ", Lstar write(*, '(A, F14.6)') "Specific Heat Ratio (gamma) = ", gamma write(*, '(A, F14.6)') "Inlet T/Tstar = ", T_Tstar1 write(*, '(A, F14.6)') "Inlet P/Pstar = ", P_Pstar1 write(*, '(A, F14.6)') "Inlet rho/rhostar = ", rho_rhostar1 write(*, '(A, F14.6)') "Inlet V/Vstar = ", V_Vstar1 write(*, '(A, F14.6)') "Inlet P0/P0star = ", P0_P0star1 write(*, '(A, F14.6)') "Exit T/Tstar = ", T_Tstar2 write(*, '(A, F14.6)') "Exit P/Pstar = ", P_Pstar2 write(*, '(A, F14.6)') "Exit rho/rhostar = ", rho_rhostar2 write(*, '(A, F14.6)') "Exit V/Vstar = ", V_Vstar2 write(*, '(A, F14.6)') "Exit P0/P0star = ", P0_P0star2 write(*, '(A, F14.6)') "Temperature Ratio (T2/T1) = ", T2_T1 write(*, '(A, F14.6)') "Pressure Ratio (P2/P1) = ", P2_P1 write(*, '(A, F14.6)') "Density Ratio (rho2/rho1) = ", rho2_rho1 write(*, '(A, F14.6)') "Stagnation Pressure Ratio (P02/P01) = ", P02_P01 write(*, '(A, F14.6)') "Entropy Change (delta_s/R) = ", delta_s_over_R if (status == 4) then write(*, '(A, F14.6)') "Shock Pre-Mach (Ms1) = ", M_shock1 write(*, '(A, F14.6)') "Shock Post-Mach (Ms2) = ", M_shock2 write(*, '(A, F14.6)') "Shock Position (x_shock) = ", x_shock end if if (T0 > 0.0d0) then write(*, '(A, F14.4)') "Inlet Temperature (T1) = ", T1 write(*, '(A, F14.4)') "Exit Temperature (T2) = ", T2 write(*, '(A, F14.2)') "Inlet Velocity (V1) = ", V1 write(*, '(A, F14.2)') "Exit Velocity (V2) = ", V2 write(*, '(A, F14.4)') "Inlet Stagnation Temp (T01) = ", T0 write(*, '(A, F14.4)') "Exit Stagnation Temp (T02) = ", T02 end if if (P0 > 0.0d0) then write(*, '(A, F14.4)') "Inlet Pressure (P1) = ", P1 write(*, '(A, F14.4)') "Exit Pressure (P2) = ", P2 write(*, '(A, F14.4)') "Inlet Stagnation Pres (P01) = ", P0 write(*, '(A, F14.4)') "Exit Stagnation Pres (P02) = ", P02 if (T0 > 0.0d0) then write(*, '(A, F14.6)') "Inlet Density (rho1) = ", rho1 write(*, '(A, F14.6)') "Exit Density (rho2) = ", rho2 end if end if contains ! Function to compute fL*/D for a given Mach number double precision function fanno_param(M, g) double precision, intent(in) :: M, g double precision :: Term1, Term2 if (M <= 0.0d0) then fanno_param = 1.0d10 else Term1 = (1.0d0 - M**2) / (g * M**2) Term2 = (g + 1.0d0) / (2.0d0 * g) * log(((g + 1.0d0) * M**2) / (2.0d0 + (g - 1.0d0) * M**2)) fanno_param = Term1 + Term2 end if end function fanno_param ! Function derivative of Fanno parameter dfL*/D / dM double precision function fanno_df(M, g) double precision, intent(in) :: M, g fanno_df = -2.0d0 * (1.0d0 - M**2) / (g * M**3 * (1.0d0 + (g - 1.0d0)/2.0d0 * M**2)) end function fanno_df ! Newton-Raphson solver to find Mach number from a target Fanno parameter value double precision function solve_mach_from_fanno(target_val, g, lower_bound, upper_bound) double precision, intent(in) :: target_val, g, lower_bound, upper_bound double precision :: M, diff, f_val, df_val integer :: i if (target_val <= 0.0d0) then solve_mach_from_fanno = 1.0d0 return end if ! Initial guess M = 0.5d0 * (lower_bound + upper_bound) if (lower_bound < 1.0d0 .and. upper_bound <= 1.0d0) then ! Subsonic guess M = 0.5d0 else if (lower_bound >= 1.0d0) then ! Supersonic guess M = 2.0d0 end if do i = 1, 100 f_val = fanno_param(M, g) - target_val df_val = fanno_df(M, g) diff = f_val / df_val M = M - diff ! Keep M in bounds if (M < lower_bound) M = lower_bound + 1.0d-5 if (M > upper_bound) M = upper_bound - 1.0d-5 if (abs(diff) < 1.0d-12) exit end do solve_mach_from_fanno = M end function solve_mach_from_fanno ! Subroutine to solve for shock Mach number and position subroutine solve_shock_position(M1, f_fac, D_val, L_val, g, Ms1, Ms2, xs) double precision, intent(in) :: M1, f_fac, D_val, L_val, g double precision, intent(out) :: Ms1, Ms2, xs double precision :: low_M, high_M, mid_M, fL_D, test_val integer :: iter fL_D = f_fac * L_val / D_val ! We use bisection solver to find Ms1 in [1.0, M1] low_M = 1.0d0 high_M = M1 do iter = 1, 100 mid_M = 0.5d0 * (low_M + high_M) ! Post-shock Mach number of mid_M Ms2 = sqrt((2.0d0 + (g - 1.0d0) * mid_M**2) / (2.0d0 * g * mid_M**2 - (g - 1.0d0))) ! Evaluate constraint: L_shock_Lstar = L*(M1) - L*(Ms1) + L*(Ms2) test_val = fanno_param(M1, g) - fanno_param(mid_M, g) + fanno_param(Ms2, g) if (test_val < fL_D) then ! Need higher mid_M (less shock distance) low_M = mid_M else ! Need lower mid_M high_M = mid_M end if if (abs(high_M - low_M) < 1.0d-12) exit end do Ms1 = mid_M Ms2 = sqrt((2.0d0 + (g - 1.0d0) * Ms1**2) / (2.0d0 * g * Ms1**2 - (g - 1.0d0))) xs = (fanno_param(M1, g) - fanno_param(Ms1, g)) * D_val / f_fac end subroutine solve_shock_position end program fanno_flow