💻 Fortran Source Code Library
We currently offer 172 open-source, production-grade Fortran codes for offline testing. Run calculations locally on your own machine, view code structure, read technical explanations, and download compilation packages including sample input files.
Fanno Flow Calculator
Core Numerical Engine in Fortran 90 • 36 total downloads
fanno_flow.f90
! =========================================================================
! Source File: 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
Solver Description
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.
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 -O3 fanno_flow.f90 -o fanno_flow_calc
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
fanno_flow_calc < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Inlet Mach number M1
2.0
! Friction factor f
0.02
! Pipe diameter D [mm]
50.0
! Pipe length L [m]
1.0
! Specific heat ratio gamma
1.4
! Inlet stagnation temp T0 [K] (optional)
0.0
! Inlet stagnation pressure P0 [kPa] (optional)
0.0
2.0
! Friction factor f
0.02
! Pipe diameter D [mm]
50.0
! Pipe length L [m]
1.0
! Specific heat ratio gamma
1.4
! Inlet stagnation temp T0 [K] (optional)
0.0
! Inlet stagnation pressure P0 [kPa] (optional)
0.0