! ============================================================================ ! 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