💻 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.
Rayleigh Flow Calculator
Core Numerical Engine in Fortran 90 • 31 total downloads
rayleigh_flow.f90
! =========================================================================
! Source File: rayleigh_flow.f90
! =========================================================================
! ============================================================================
! ThermoFluidCalc — Rayleigh Flow Solver
! Reference: Shapiro, Dynamics and Thermodynamics of Compressible Fluid Flow
! ============================================================================
program rayleigh_flow
implicit none
! Input variables
double precision :: M1 ! Inlet Mach number
double precision :: gamma ! Specific heat ratio (default 1.4)
double precision :: T01 ! Inlet Stagnation Temp [K] (0 = skip)
double precision :: P01 ! Inlet Stagnation Pressure [kPa] (0 = skip)
double precision :: q ! Heat added per unit mass [kJ/kg] (0 = skip)
! Constants
double precision, parameter :: R_air = 287.05d0 ! Gas constant [J/kg-K]
! Intermediate and output variables
double precision :: Cp, Cp_kJ
double precision :: T0star, q_max, q_Joules
double precision :: M2 ! Exit Mach number
double precision :: T2_T1, T02_T01, P2_P1, P02_P01, rho2_rho1
double precision :: T1, P1, rho1, V1, T2, P2, rho2, V2
double precision :: P02, T02 ! exit stagnation properties
double precision :: T_Tstar1, T0_T0star1, P_Pstar1, P0_P0star1, V_Vstar1, rho_rhostar1
double precision :: T_Tstar2, T0_T0star2, P_Pstar2, P0_P0star2, V_Vstar2, rho_rhostar2
double precision :: delta_s_over_R
double precision :: phi_target
double precision :: M1_orig
double precision :: P01_eff
integer :: status ! 1=Subsonic unchoked, 2=Supersonic unchoked, 3=Thermally choked
! Read input from stdin
read(*,*) M1
read(*,*) gamma
read(*,*) T01
read(*,*) P01
read(*,*) q
! Set defaults and bounds
if (gamma <= 1.0d0) gamma = 1.4d0
if (M1 <= 0.0d0) M1 = 0.5d0
if (q < 0.0d0) q = 0.0d0
! If heat is added, we must have reference stagnation temperature and pressure
! If the user did not provide them, we default them
if (T01 <= 0.0d0) T01 = 300.0d0
if (P01 <= 0.0d0) P01 = 101.325d0
M1_orig = M1
Cp = gamma * R_air / (gamma - 1.0d0) ! J/kg-K
Cp_kJ = Cp / 1000.0d0 ! kJ/kg-K
T0star = T01 / rayleigh_t0_ratio(M1, gamma)
q_max = Cp_kJ * (T0star - T01)
! Initialize shock effects for supersonic choking
P01_eff = P01
! ── DETERMINE FLOW STATUS AND SOLVE EXIT MACH ────────────
if (q <= q_max) then
! FLOW IS UNCHOKED
if (M1 < 1.0d0) then
status = 1 ! Subsonic unchoked
else
status = 2 ! Supersonic unchoked
end if
T02 = T01 + q / Cp_kJ
phi_target = T02 / T0star
if (M1 < 1.0d0) then
M2 = solve_mach_from_t0_ratio(phi_target, gamma, .true., M1)
else
M2 = solve_mach_from_t0_ratio(phi_target, gamma, .false., M1)
end if
else
! FLOW IS THERMALLY CHOKED
status = 3
M2 = 1.0d0
T02 = T01 + q / Cp_kJ ! Exit stagnation temperature
! The new inlet flow must become subsonic and adapt to a lower Mach number.
! This holds true for both subsonic and supersonic initial states.
phi_target = T01 / T02
M1 = solve_mach_from_t0_ratio(phi_target, gamma, .true., 0.5d0)
! Stagnation temperature at inlet remains constant
! Stagnation pressure at inlet remains constant if subsonic, or suffers shock loss if supersonic
if (M1_orig > 1.0d0) then
P01_eff = P01 * shock_p0_ratio(M1_orig, gamma)
else
P01_eff = P01
end if
! Re-calculate choking conditions for the new Rayleigh line
T0star = T01 / rayleigh_t0_ratio(M1, gamma)
q_max = Cp_kJ * (T0star - T01)
end if
! ── COMPUTE RAYLEIGH RATIOS AT INLET AND EXIT ────────────
T_Tstar1 = ((gamma + 1.0d0) * M1 / (1.0d0 + gamma * M1**2))**2
T0_T0star1 = rayleigh_t0_ratio(M1, gamma)
P_Pstar1 = (gamma + 1.0d0) / (1.0d0 + gamma * M1**2)
P0_P0star1 = P_Pstar1 * ((2.0d0 + (gamma - 1.0d0) * M1**2) / (gamma + 1.0d0))**(gamma / (gamma - 1.0d0))
V_Vstar1 = (gamma + 1.0d0) * M1**2 / (1.0d0 + gamma * M1**2)
rho_rhostar1 = 1.0d0 / V_Vstar1
T_Tstar2 = ((gamma + 1.0d0) * M2 / (1.0d0 + gamma * M2**2))**2
T0_T0star2 = rayleigh_t0_ratio(M2, gamma)
P_Pstar2 = (gamma + 1.0d0) / (1.0d0 + gamma * M2**2)
P0_P0star2 = P_Pstar2 * ((2.0d0 + (gamma - 1.0d0) * M2**2) / (gamma + 1.0d0))**(gamma / (gamma - 1.0d0))
V_Vstar2 = (gamma + 1.0d0) * M2**2 / (1.0d0 + gamma * M2**2)
rho_rhostar2 = 1.0d0 / V_Vstar2
! Property changes across duct
T2_T1 = T_Tstar2 / T_Tstar1
T02_T01 = T0_T0star2 / T0_T0star1
P2_P1 = P_Pstar2 / P_Pstar1
P02_P01 = (P01_eff / P01) * (P0_P0star2 / P0_P0star1)
rho2_rho1 = rho_rhostar2 / rho_rhostar1
delta_s_over_R = (gamma / (gamma - 1.0d0)) * log(T2_T1) - log(P2_P1)
! ── COMPUTE ACTUAL PHYSICAL PROPERTIES ───────────────────
T1 = T01 / (1.0d0 + (gamma - 1.0d0) / 2.0d0 * M1**2)
V1 = M1 * sqrt(gamma * R_air * T1)
P1 = P01_eff / (1.0d0 + (gamma - 1.0d0) / 2.0d0 * M1**2)**(gamma / (gamma - 1.0d0))
rho1 = P1 / (R_air * T1 / 1000.0d0)
T2 = T1 * T2_T1
V2 = M2 * sqrt(gamma * R_air * T2)
P2 = P1 * P2_P1
rho2 = P2 / (R_air * T2 / 1000.0d0)
P02 = P01_eff * (P0_P0star2 / P0_P0star1)
! ── 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)
if (M1_orig < 1.0d0) then
write(*, '(A)') "Status = Thermally Choked Subsonic Flow"
else
write(*, '(A)') "Status = Thermally Choked Supersonic Flow (Shocked)"
end if
end select
write(*, '(A, F14.6)') "Inlet Mach (M1) = ", M1
write(*, '(A, F14.6)') "Exit Mach (M2) = ", M2
write(*, '(A, F14.6)') "Original Inlet Mach = ", M1_orig
write(*, '(A, F14.6)') "Heat Added (q) = ", q
write(*, '(A, F14.6)') "Max Heat Added (q_max) = ", q_max
write(*, '(A, F14.6)') "Specific Heat Ratio (gamma) = ", gamma
write(*, '(A, F14.6)') "Inlet T/Tstar = ", T_Tstar1
write(*, '(A, F14.6)') "Inlet T0/T0star = ", T0_T0star1
write(*, '(A, F14.6)') "Inlet P/Pstar = ", P_Pstar1
write(*, '(A, F14.6)') "Inlet P0/P0star = ", P0_P0star1
write(*, '(A, F14.6)') "Inlet V/Vstar = ", V_Vstar1
write(*, '(A, F14.6)') "Inlet rho/rhostar = ", rho_rhostar1
write(*, '(A, F14.6)') "Exit T/Tstar = ", T_Tstar2
write(*, '(A, F14.6)') "Exit T0/T0star = ", T0_T0star2
write(*, '(A, F14.6)') "Exit P/Pstar = ", P_Pstar2
write(*, '(A, F14.6)') "Exit P0/P0star = ", P0_P0star2
write(*, '(A, F14.6)') "Exit V/Vstar = ", V_Vstar2
write(*, '(A, F14.6)') "Exit rho/rhostar = ", rho_rhostar2
write(*, '(A, F14.6)') "Temperature Ratio (T2/T1) = ", T2_T1
write(*, '(A, F14.6)') "Stagnation Temp Ratio (T02/T01) = ", T02_T01
write(*, '(A, F14.6)') "Pressure Ratio (P2/P1) = ", P2_P1
write(*, '(A, F14.6)') "Stagnation Pressure Ratio (P02/P01) = ", P02_P01
write(*, '(A, F14.6)') "Density Ratio (rho2/rho1) = ", rho2_rho1
write(*, '(A, F14.6)') "Entropy Change (delta_s/R) = ", delta_s_over_R
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 Pressure (P1) = ", P1
write(*, '(A, F14.4)') "Exit Pressure (P2) = ", P2
write(*, '(A, F14.4)') "Inlet Stagnation Temp (T01) = ", T01
write(*, '(A, F14.4)') "Exit Stagnation Temp (T02) = T02"
write(*, '(A, F14.4)') "Inlet Stagnation Pres (P01) = ", P01
write(*, '(A, F14.4)') "Exit Stagnation Pres (P02) = ", P02
write(*, '(A, F14.6)') "Inlet Density (rho1) = ", rho1
write(*, '(A, F14.6)') "Exit Density (rho2) = ", rho2
contains
! Stagnation temperature ratio function T0/T0*
double precision function rayleigh_t0_ratio(M, g)
double precision, intent(in) :: M, g
rayleigh_t0_ratio = (2.0d0 * (g + 1.0d0) * M**2 * (1.0d0 + (g - 1.0d0) / 2.0d0 * M**2)) / (1.0d0 + g * M**2)**2
end function rayleigh_t0_ratio
! Stagnation pressure ratio across a normal shock P02/P01
double precision function shock_p0_ratio(M, g)
double precision, intent(in) :: M, g
double precision :: t1, t2
t1 = (g + 1.0d0) / (2.0d0 * g * M**2 - (g - 1.0d0))
t2 = ((g + 1.0d0) * M**2) / (2.0d0 + (g - 1.0d0) * M**2)
shock_p0_ratio = t1**(1.0d0 / (g - 1.0d0)) * t2**(g / (g - 1.0d0))
end function shock_p0_ratio
! Solve Mach number from Stagnation Temperature Ratio using Newton-Raphson
double precision function solve_mach_from_t0_ratio(target_phi, g, is_sub, init_guess)
double precision, intent(in) :: target_phi, g, init_guess
logical, intent(in) :: is_sub
double precision :: M, M_next, f_val, df_val
double precision :: num, den, dnum, dden
integer :: iter
M = init_guess
if (is_sub .and. M >= 1.0d0) M = 0.5d0
if (.not. is_sub .and. M <= 1.0d0) M = 2.0d0
! Safeguard target_phi
if (target_phi >= 0.99999d0) then
solve_mach_from_t0_ratio = 1.0d0
return
end if
do iter = 1, 100
! We solve N(M) - target_phi * D(M) = 0
! N(M) = 2(g+1)M^2 + (g^2 - 1)M^4
! D(M) = (1 + g M^2)^2
num = 2.0d0 * (g + 1.0d0) * M**2 + (g**2 - 1.0d0) * M**4
den = (1.0d0 + g * M**2)**2
f_val = num - target_phi * den
! Derivatives:
! N'(M) = 4(g+1)M + 4(g^2 - 1)M^3
! D'(M) = 4 g M (1 + g M^2)
dnum = 4.0d0 * (g + 1.0d0) * M + 4.0d0 * (g**2 - 1.0d0) * M**3
dden = 4.0d0 * g * M * (1.0d0 + g * M**2)
df_val = dnum - target_phi * dden
if (abs(df_val) < 1.0d-12) exit
M_next = M - f_val / df_val
! Keep guess in bounds
if (is_sub) then
if (M_next <= 0.0d0) M_next = M / 2.0d0
if (M_next >= 1.0d0) M_next = (M + 1.0d0) / 2.0d0
else
if (M_next <= 1.0d0) M_next = (M + 1.0d0) / 2.0d0
end if
if (abs(M_next - M) < 1.0d-8) then
M = M_next
exit
end if
M = M_next
end do
solve_mach_from_t0_ratio = M
end function solve_mach_from_t0_ratio
end program rayleigh_flow
Solver Description
Solves 1D compressible flow with heat addition or extraction (Rayleigh Flow) for both subsonic and supersonic regimes. Computes Rayleigh line ratios (T/T*, T0/T0*, P/P*, P0/P0*, V/V*), maximum heat addition before choking, and entropy changes. Resolves exit Mach numbers for arbitrary heat inputs, determines thermal choking limits, and computes adjusted subsonic inlet states and shock losses for choked 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 rayleigh_flow.f90 -o rayleigh_flow_calc
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
rayleigh_flow_calc < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Inlet Mach number M1
0.5
! Specific heat ratio gamma
1.4
! Inlet stagnation temp T0 [K]
300.0
! Inlet stagnation pressure P0 [kPa]
101.325
! Heat added q [kJ/kg]
50.0
0.5
! Specific heat ratio gamma
1.4
! Inlet stagnation temp T0 [K]
300.0
! Inlet stagnation pressure P0 [kPa]
101.325
! Heat added q [kJ/kg]
50.0