🔬
Solver Purpose & Physical Scope
Solves 1D compressible flow with pipe friction (Fanno Flow) for both subsonic and supersonic regimes. Computes Fanno line ratios (T/T*, P/P*, rho/rho*, V/V*, P0/P0*), choking pipe length, and entropy changes. Resolves exit Mach numbers for arbitrary duct lengths, determines choking limits, and performs numerical bisection to locate normal shock waves standing within supersonic ducts.
📂 Discipline: Fluid
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 410 times
📄 Source File:
fanno_flow.f90
📁 calcul/FluidMechanics /
fanno_flow.f90
! ============================================================================
! 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
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 fanno_flow.f90 -o fanno_flow_calc
2. Execution with input.txt redirection:
fanno_flow_calc < input.txt
📄 Sample input.txt File Structure
Sample Data:
2.0 0.02 50.0 1.0 1.4 0.0 0.0
Parameter Description:
Inlet Mach number M1 Friction factor f Pipe diameter D [mm] Pipe length L [m] Specific heat ratio gamma Inlet stagnation temp T0 [K] (optional) Inlet stagnation pressure P0 [kPa] (optional)