💻 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.
Prandtl-Meyer Expansion Fan
Core Numerical Engine in Fortran 90 • 53 total downloads
prandtl_meyer.f90
! =========================================================================
! Source File: prandtl_meyer.f90
! =========================================================================
! ============================================================================
! ThermoFluidCalc — Prandtl-Meyer Expansion Fan Solver
! Reference: Anderson Modern Compressible Flow Ch. 4
! ============================================================================
program prandtl_meyer
implicit none
! Input variables
integer :: solve_for ! 1 = Given M1 & theta, solve M2; 2 = Given M1 & M2, solve theta; 3 = Given M2 & theta, solve M1
double precision :: M1 ! Upstream Mach number
double precision :: theta ! Deflection angle [deg]
double precision :: M2 ! Downstream Mach number
double precision :: gamma ! Specific heat ratio (default 1.4)
double precision :: T1 ! Upstream Temperature [K] (0 = skip)
double precision :: P1 ! Upstream Pressure [kPa] (0 = skip)
! Constants
double precision, parameter :: pi = 3.141592653589793d0
double precision, parameter :: rad = pi / 180.0d0
double precision, parameter :: R_air = 287.0d0 ! Gas constant for air [J/kg-K]
! Intermediate and output variables
double precision :: nu1_rad, nu2_rad, theta_rad
double precision :: mu1, mu2, nu_max_rad, nu_max_deg
double precision :: P2_P1, T2_T1, rho2_rho1
double precision :: T2, P2, a1, a2, V1, V2
! Read input from stdin
read(*,*) solve_for
read(*,*) M1
read(*,*) theta
read(*,*) M2
read(*,*) gamma
read(*,*) T1
read(*,*) P1
! Check inputs
if (gamma <= 1.0d0) gamma = 1.4d0
! Compute nu_max analytically for the given gamma
nu_max_rad = (sqrt((gamma + 1.0d0) / (gamma - 1.0d0)) - 1.0d0) * (pi / 2.0d0)
nu_max_deg = nu_max_rad / rad
! ── PERFORM SOLVING ACCORDING TO MODE ────────────────────
select case (solve_for)
case (1) ! Given M1 and theta, solve M2
if (M1 <= 1.0d0) then
write(*,*) "ERROR: Upstream Mach number M1 must be supersonic (> 1.0) for Prandtl-Meyer expansion."
stop
end if
nu1_rad = nu_mach(M1, gamma)
theta_rad = theta * rad
nu2_rad = nu1_rad + theta_rad
if (nu2_rad >= nu_max_rad) then
write(*,*) "ERROR: Deflection angle theta exceeds maximum possible turn angle for M1 = ", M1
stop
end if
M2 = solve_mach_from_nu(nu2_rad, gamma, M1)
case (2) ! Given M1 and M2, solve theta
if (M1 <= 1.0d0 .or. M2 <= 1.0d0) then
write(*,*) "ERROR: Mach numbers must be supersonic (> 1.0) for expansion waves."
stop
end if
if (M2 < M1) then
write(*,*) "ERROR: Downstream Mach number M2 must be greater than upstream Mach number M1."
stop
end if
nu1_rad = nu_mach(M1, gamma)
nu2_rad = nu_mach(M2, gamma)
theta_rad = nu2_rad - nu1_rad
theta = theta_rad / rad
case (3) ! Given M2 and theta, solve M1
if (M2 <= 1.0d0) then
write(*,*) "ERROR: Downstream Mach number M2 must be supersonic (> 1.0)."
stop
end if
nu2_rad = nu_mach(M2, gamma)
theta_rad = theta * rad
nu1_rad = nu2_rad - theta_rad
if (nu1_rad < 0.0d0) then
write(*,*) "ERROR: Deflection angle theta exceeds total downstream PM angle."
stop
end if
M1 = solve_mach_from_nu(nu1_rad, gamma, 1.5d0)
case default
write(*,*) "ERROR: Invalid solve target mode."
stop
end select
! ── COMPUTE DOWNSTREAM PROPERTIES ────────────────────────
nu1_rad = nu_mach(M1, gamma)
nu2_rad = nu_mach(M2, gamma)
mu1 = asin(1.0d0 / M1) / rad
mu2 = asin(1.0d0 / M2) / rad
! Isentropic ratios
T2_T1 = (1.0d0 + (gamma - 1.0d0)/2.0d0 * M1**2) / (1.0d0 + (gamma - 1.0d0)/2.0d0 * M2**2)
P2_P1 = (T2_T1)**(gamma / (gamma - 1.0d0))
rho2_rho1 = (T2_T1)**(1.0d0 / (gamma - 1.0d0))
! Optional static temperatures/pressures calculation
T2 = 0.0d0
P2 = 0.0d0
V1 = 0.0d0
V2 = 0.0d0
if (T1 > 0.0d0) then
T2 = T1 * T2_T1
a1 = sqrt(gamma * R_air * T1)
a2 = sqrt(gamma * R_air * T2)
V1 = M1 * a1
V2 = M2 * a2
end if
if (P1 > 0.0d0) then
P2 = P1 * P2_P1
end if
! ── OUTPUT RESULTS IN KEY-VALUE FORMAT ───────────────────
write(*, '(A, F14.6)') "Upstream Mach (M1) = ", M1
write(*, '(A, F14.6)') "Downstream Mach (M2) = ", M2
write(*, '(A, F14.6)') "Deflection Angle (theta) = ", theta
write(*, '(A, F14.6)') "PM Angle Upstream (nu1) = ", nu1_rad / rad
write(*, '(A, F14.6)') "PM Angle Downstream (nu2) = ", nu2_rad / rad
write(*, '(A, F14.6)') "Mach Angle Upstream (mu1) = ", mu1
write(*, '(A, F14.6)') "Mach Angle Downstream (mu2) = ", mu2
write(*, '(A, F14.6)') "Max Turning Angle (nu_max) = ", nu_max_deg
write(*, '(A, F14.6)') "Pressure Ratio (P2/P1) = ", P2_P1
write(*, '(A, F14.6)') "Temperature Ratio (T2/T1) = ", T2_T1
write(*, '(A, F14.6)') "Density Ratio (rho2/rho1) = ", rho2_rho1
write(*, '(A, F14.6)') "Specific Heat Ratio (gamma) = ", gamma
if (T1 > 0.0d0) then
write(*, '(A, F14.4)') "Upstream Temp (T1) = ", T1
write(*, '(A, F14.4)') "Downstream Temp (T2) = ", T2
write(*, '(A, F14.2)') "Upstream Velocity (V1) = ", V1
write(*, '(A, F14.2)') "Downstream Velocity (V2) = ", V2
end if
if (P1 > 0.0d0) then
write(*, '(A, F14.4)') "Upstream Pressure (P1) = ", P1
write(*, '(A, F14.4)') "Downstream Pressure (P2) = ", P2
end if
contains
! Prandtl-Meyer function calculation: nu(M) in radians
double precision function nu_mach(M, g)
double precision, intent(in) :: M, g
double precision :: factor
if (M <= 1.0d0) then
nu_mach = 0.0d0
else
factor = sqrt((g + 1.0d0) / (g - 1.0d0))
nu_mach = factor * atan(sqrt((g - 1.0d0) / (g + 1.0d0) * (M**2 - 1.0d0))) - &
atan(sqrt(M**2 - 1.0d0))
end if
end function nu_mach
! Derivative of the Prandtl-Meyer function: dnu/dM
double precision function dnu_dM(M, g)
double precision, intent(in) :: M, g
if (M <= 1.0d0) then
dnu_dM = 1.0d-10 ! Avoid division by zero
else
dnu_dM = sqrt(M**2 - 1.0d0) / (M * (1.0d0 + (g - 1.0d0) / 2.0d0 * M**2))
end if
end function dnu_dM
! Newton-Raphson solver to find M from a target nu in radians
double precision function solve_mach_from_nu(nu_target, g, start_M)
double precision, intent(in) :: nu_target, g, start_M
double precision :: M, diff, f_val, df_val
integer :: i
if (nu_target <= 1.0d-12) then
solve_mach_from_nu = 1.0d0
return
end if
! Initialize solver with starting Mach or standard guess
M = start_M
if (M <= 1.0d0) M = 2.0d0
do i = 1, 100
f_val = nu_mach(M, g) - nu_target
df_val = dnu_dM(M, g)
diff = f_val / df_val
M = M - diff
if (M < 1.0d0) M = 1.0001d0 ! Floor Mach to supersonic limit
if (abs(diff) < 1.0d-14) exit
end do
solve_mach_from_nu = M
end function solve_mach_from_nu
end program prandtl_meyer
Solver Description
Solves supersonic flow expansion around a convex corner. Using the Prandtl-Meyer function, it relates the deflection angle to the upstream and downstream Mach numbers. A high-precision Newton-Raphson iterative solver is used to find the exit Mach number. It also evaluates isentropic static pressure, temperature, and density ratios, as well as the orientations of the leading and trailing Mach waves.
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 prandtl_meyer.f90 -o prandtl_meyer_calc
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
prandtl_meyer_calc < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Solve mode (1=Solve M2, 2=Solve theta, 3=Solve M1)
1
! Upstream Mach M1
2.0
! Deflection angle theta [deg]
15.0
! Downstream Mach M2
0.0
! Specific heat ratio gamma
1.4
! Upstream Temp T1 [K]
288.15
! Upstream Pressure P1 [kPa]
101.325
1
! Upstream Mach M1
2.0
! Deflection angle theta [deg]
15.0
! Downstream Mach M2
0.0
! Specific heat ratio gamma
1.4
! Upstream Temp T1 [K]
288.15
! Upstream Pressure P1 [kPa]
101.325