๐ป Fortran Source Code Library
Run calculations locally on your own machine. View code structure, read technical explanations, and download compilation packages including sample input files.
Convergent-Divergent Nozzle Flow Sizer
Core Numerical Engine in Fortran 90 โข 0 total downloads
nozzle_flow.f90
! =========================================================================
! Source File: nozzle_flow.f90
! =========================================================================
program nozzle_flow
implicit none
double precision :: gamma, T0, P0, A_throat, A_exit
double precision :: P_back, AR, R_gas
double precision :: pi
! Isentropic properties
double precision :: M_sub, M_sup, P_exit_sub, P_exit_sup
double precision :: T_exit, rho_exit, V_exit, snd_exit
double precision :: mdot, mdot_choked
double precision :: P_e, T_e, M_e
! Shock variables
double precision :: M_shock_up, M_shock_down
double precision :: P2_P1, T2_T1, P02_P01
double precision :: A_shock, AR_shock
double precision :: P_exit_shock
! Regime detection
integer :: regime ! 1-6
character(len=60) :: regime_name
double precision :: P_crit1, P_crit2, P_crit3
! Chart data
integer :: i, n_pts
double precision :: x, AR_x, M_x, P_x, T_x
double precision :: x_throat, x_shock
! Bisection variables
double precision :: M_lo, M_hi, M_mid, AR_mid
integer :: j
pi = 3.14159265358979d0
R_gas = 287.0d0 ! J/(kgยทK) for air
n_pts = 60
! Read inputs
read(*,*) gamma ! Specific heat ratio
read(*,*) T0 ! Stagnation temperature [K]
read(*,*) P0 ! Stagnation pressure [kPa]
read(*,*) A_throat ! Throat area [cm2]
read(*,*) A_exit ! Exit area [cm2]
read(*,*) P_back ! Back pressure [kPa]
AR = A_exit / A_throat
! ============================================
! FIND SUBSONIC AND SUPERSONIC EXIT MACH
! ============================================
! Subsonic solution (bisection: M from 0.01 to 1.0)
M_lo = 0.001d0; M_hi = 1.0d0
do j = 1, 200
M_mid = (M_lo + M_hi) / 2.0d0
AR_mid = area_mach_ratio(M_mid, gamma)
if (AR_mid > AR) then
M_lo = M_mid
else
M_hi = M_mid
end if
end do
M_sub = M_mid
! Supersonic solution (bisection: M from 1.0 to 20)
M_lo = 1.0d0; M_hi = 20.0d0
do j = 1, 200
M_mid = (M_lo + M_hi) / 2.0d0
AR_mid = area_mach_ratio(M_mid, gamma)
if (AR_mid < AR) then
M_lo = M_mid
else
M_hi = M_mid
end if
end do
M_sup = M_mid
! Exit pressures for subsonic and supersonic isentropic cases
P_exit_sub = P0 / isen_P_ratio(M_sub, gamma)
P_exit_sup = P0 / isen_P_ratio(M_sup, gamma)
! Normal shock at exit: compute post-shock pressure
M_shock_down = sqrt(((gamma-1.0d0)*M_sup**2 + 2.0d0) / &
(2.0d0*gamma*M_sup**2 - (gamma-1.0d0)))
P2_P1 = 1.0d0 + 2.0d0*gamma/(gamma+1.0d0) * (M_sup**2 - 1.0d0)
P_exit_shock = P_exit_sup * P2_P1
! Critical back pressures for regime boundaries
P_crit1 = P_exit_sub ! Below: choked flow begins
P_crit2 = P_exit_shock ! Below: shock moves into nozzle โ exits
P_crit3 = P_exit_sup ! Below: overexpanded โ design โ underexpanded
! ============================================
! DETERMINE OPERATING REGIME
! ============================================
if (P_back >= P0) then
regime = 1
regime_name = 'No flow (Pb >= P0)'
M_e = 0.0d0
P_e = P_back
T_e = T0
else if (P_back >= P_crit1) then
regime = 2
regime_name = 'Subsonic throughout (not choked)'
! Find actual exit Mach from P_back
M_e = mach_from_pressure(P_back, P0, gamma)
P_e = P_back
T_e = T0 / (1.0d0 + (gamma-1.0d0)/2.0d0 * M_e**2)
else if (P_back >= P_crit2) then
regime = 3
regime_name = 'Normal shock in diverging section'
! Need to find shock location
M_e = M_sub ! Will be adjusted
P_e = P_back
T_e = T0 / (1.0d0 + (gamma-1.0d0)/2.0d0 * M_e**2)
else if (P_back > P_crit3) then
regime = 4
regime_name = 'Overexpanded (oblique shocks at exit)'
M_e = M_sup
P_e = P_exit_sup
T_e = T0 / (1.0d0 + (gamma-1.0d0)/2.0d0 * M_e**2)
else if (abs(P_back - P_crit3) < 0.001d0) then
regime = 5
regime_name = 'Design condition (perfectly expanded)'
M_e = M_sup
P_e = P_exit_sup
T_e = T0 / (1.0d0 + (gamma-1.0d0)/2.0d0 * M_e**2)
else
regime = 6
regime_name = 'Underexpanded (expansion fans at exit)'
M_e = M_sup
P_e = P_exit_sup
T_e = T0 / (1.0d0 + (gamma-1.0d0)/2.0d0 * M_e**2)
end if
! Mass flow rate
if (P_back < P_crit1) then
! Choked
mdot_choked = P0*1000.0d0 * (A_throat*1.0d-4) * &
sqrt(gamma/(R_gas*T0)) * (2.0d0/(gamma+1.0d0))**((gamma+1.0d0)/(2.0d0*(gamma-1.0d0)))
mdot = mdot_choked
else
! Not choked โ subsonic throughout
mdot = P0*1000.0d0 * (A_throat*1.0d-4) * M_e * &
sqrt(gamma/(R_gas*T0)) / isen_P_ratio(M_e, gamma) * &
sqrt(1.0d0 + (gamma-1.0d0)/2.0d0*M_e**2)
mdot_choked = P0*1000.0d0 * (A_throat*1.0d-4) * &
sqrt(gamma/(R_gas*T0)) * (2.0d0/(gamma+1.0d0))**((gamma+1.0d0)/(2.0d0*(gamma-1.0d0)))
end if
! ============================================
! OUTPUT
! ============================================
write(*,'(A)') '============================================================'
write(*,'(A)') ' CONVERGING-DIVERGING NOZZLE FLOW ANALYSIS'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- INPUT PARAMETERS ---'
write(*,'(A,F12.4)') ' Specific heat ratio ฮณ = ', gamma
write(*,'(A,F12.2,A)') ' Stagnation temp Tโ = ', T0, ' K'
write(*,'(A,F12.4,A)') ' Stagnation pressure Pโ = ', P0, ' kPa'
write(*,'(A,F12.4,A)') ' Throat area A* = ', A_throat, ' cmยฒ'
write(*,'(A,F12.4,A)') ' Exit area A_e = ', A_exit, ' cmยฒ'
write(*,'(A,F12.4)') ' Area ratio A_e/A* = ', AR
write(*,'(A,F12.4,A)') ' Back pressure P_b = ', P_back, ' kPa'
write(*,*)
write(*,'(A)') '--- OPERATING REGIME ---'
write(*,'(A,I1,A,A)') ' Regime ', regime, ': ', trim(regime_name)
write(*,*)
write(*,'(A)') '--- CRITICAL BACK PRESSURES ---'
write(*,'(A,F12.4,A)') ' P_sub (choke onset) = ', P_crit1, ' kPa'
write(*,'(A,F12.4,A)') ' P_shock (shock at exit) = ', P_crit2, ' kPa'
write(*,'(A,F12.4,A)') ' P_design (isentropic) = ', P_crit3, ' kPa'
write(*,*)
write(*,'(A)') '--- ISENTROPIC EXIT SOLUTIONS ---'
write(*,'(A,F12.6)') ' M_exit (subsonic) = ', M_sub
write(*,'(A,F12.6)') ' M_exit (supersonic) = ', M_sup
write(*,'(A,F12.4,A)') ' P_exit (subsonic) = ', P_exit_sub, ' kPa'
write(*,'(A,F12.4,A)') ' P_exit (supersonic) = ', P_exit_sup, ' kPa'
write(*,*)
write(*,'(A)') '--- ACTUAL EXIT CONDITIONS ---'
write(*,'(A,F12.6)') ' Exit Mach number = ', M_e
write(*,'(A,F12.4,A)') ' Exit pressure = ', P_e, ' kPa'
write(*,'(A,F12.2,A)') ' Exit temperature = ', T_e, ' K'
if (M_e > 0.0d0) then
snd_exit = sqrt(gamma * R_gas * T_e)
V_exit = M_e * snd_exit
write(*,'(A,F12.2,A)') ' Exit velocity = ', V_exit, ' m/s'
write(*,'(A,F12.2,A)') ' Speed of sound (exit) = ', snd_exit, ' m/s'
end if
write(*,*)
write(*,'(A)') '--- MASS FLOW RATE ---'
write(*,'(A,F12.6,A)') ' แน (actual) = ', mdot, ' kg/s'
write(*,'(A,F12.6,A)') ' แน (choked max) = ', mdot_choked, ' kg/s'
if (P_back < P_crit1) then
write(*,'(A)') ' Status: CHOKED (แน = แน_max)'
else
write(*,'(A,F8.2,A)') ' Status: NOT CHOKED (', (mdot/mdot_choked)*100.0d0, '% of max)'
end if
write(*,*)
! Normal shock info if applicable
if (regime == 3) then
write(*,'(A)') '--- NORMAL SHOCK IN NOZZLE ---'
! Find shock location by bisection
call find_shock_location(AR, P_back, P0, gamma, A_throat, A_shock, AR_shock, &
M_shock_up, M_shock_down)
P2_P1 = 1.0d0 + 2.0d0*gamma/(gamma+1.0d0) * (M_shock_up**2 - 1.0d0)
T2_T1 = (1.0d0 + 2.0d0*gamma/(gamma+1.0d0)*(M_shock_up**2-1.0d0)) * &
((2.0d0+(gamma-1.0d0)*M_shock_up**2)/((gamma+1.0d0)*M_shock_up**2))
P02_P01 = (((gamma+1.0d0)*M_shock_up**2)/((gamma-1.0d0)*M_shock_up**2+2.0d0)) &
**(gamma/(gamma-1.0d0)) * &
(1.0d0/(1.0d0+2.0d0*gamma/(gamma+1.0d0)*(M_shock_up**2-1.0d0))) &
**(1.0d0/(gamma-1.0d0))
write(*,'(A,F12.4)') ' Shock at A/A* = ', AR_shock
write(*,'(A,F12.6)') ' M upstream of shock = ', M_shock_up
write(*,'(A,F12.6)') ' M downstream of shock = ', M_shock_down
write(*,'(A,F12.4)') ' P2/P1 across shock = ', P2_P1
write(*,'(A,F12.4)') ' T2/T1 across shock = ', T2_T1
write(*,'(A,F12.6)') ' P02/P01 (stag. loss) = ', P02_P01
write(*,'(A,F8.2,A)') ' Stag. pressure loss = ', (1.0d0-P02_P01)*100.0d0, ' %'
write(*,*)
end if
! ============================================
! NOZZLE PROFILE TABLE
! ============================================
write(*,'(A)') '--- NOZZLE FLOW PROFILE ---'
write(*,'(A)') ' x/L A/A* M P/P0 T/T0'
write(*,'(A)') ' ------ ------ ------ ------ ------'
x_throat = 0.5d0 ! Throat at mid-nozzle
do i = 1, n_pts
x = dble(i) / dble(n_pts)
! Area ratio: linear converging then diverging
if (x <= x_throat) then
AR_x = AR + (1.0d0 - AR) * (x / x_throat)
else
AR_x = 1.0d0 + (AR - 1.0d0) * ((x - x_throat) / (1.0d0 - x_throat))
end if
if (AR_x < 1.0d0) AR_x = 1.0d0
! Find Mach based on regime
if (regime <= 2) then
M_x = mach_from_area_sub(AR_x, gamma)
else if (x <= x_throat) then
M_x = mach_from_area_sub(AR_x, gamma)
else
if (regime >= 4) then
M_x = mach_from_area_sup(AR_x, gamma)
else
M_x = mach_from_area_sup(AR_x, gamma)
end if
end if
P_x = 1.0d0 / isen_P_ratio(M_x, gamma)
T_x = 1.0d0 / (1.0d0 + (gamma-1.0d0)/2.0d0 * M_x**2)
write(*,'(F7.4,3X,F8.4,3X,F8.4,3X,F8.4,3X,F8.4)') x, AR_x, M_x, P_x, T_x
end do
write(*,*)
! ============================================
! CHART DATA
! ============================================
write(*,'(A)') '--- CHART_DATA ---'
do i = 1, n_pts
x = dble(i) / dble(n_pts)
if (x <= x_throat) then
AR_x = AR + (1.0d0 - AR) * (x / x_throat)
else
AR_x = 1.0d0 + (AR - 1.0d0) * ((x - x_throat) / (1.0d0 - x_throat))
end if
if (AR_x < 1.0d0) AR_x = 1.0d0
if (regime <= 2) then
M_x = mach_from_area_sub(AR_x, gamma)
else if (x <= x_throat) then
M_x = mach_from_area_sub(AR_x, gamma)
else
if (regime >= 4) then
M_x = mach_from_area_sup(AR_x, gamma)
else
M_x = mach_from_area_sup(AR_x, gamma)
end if
end if
P_x = 1.0d0 / isen_P_ratio(M_x, gamma)
T_x = 1.0d0 / (1.0d0 + (gamma-1.0d0)/2.0d0 * M_x**2)
write(*,'(F10.6,A,F10.6,A,F10.6,A,F10.6,A,F10.6)') &
x, ',', AR_x, ',', M_x, ',', P_x, ',', T_x
end do
write(*,'(A)') '--- END_CHART_DATA ---'
write(*,*)
write(*,'(A)') '--- EQUATIONS USED ---'
write(*,'(A)') ' A/A* = (1/M)[(2/(g+1))(1+(g-1)/2 M^2)]^((g+1)/(2(g-1)))'
write(*,'(A)') ' T0/T = 1 + (g-1)/2 x M^2'
write(*,'(A)') ' P0/P = (T0/T)^(g/(g-1))'
write(*,'(A)') ' mdot_choked = P0*A*sqrt(g/(R*T0)) x (2/(g+1))^((g+1)/(2(g-1)))'
write(*,'(A)') ' Normal shock: M2^2 = [(g-1)M1^2+2]/[2g*M1^2-(g-1)]'
contains
! Area-Mach ratio function
function area_mach_ratio(M, g) result(ar)
double precision, intent(in) :: M, g
double precision :: ar
ar = (1.0d0/M) * ((2.0d0/(g+1.0d0)) * &
(1.0d0 + (g-1.0d0)/2.0d0 * M**2))**((g+1.0d0)/(2.0d0*(g-1.0d0)))
end function
! Isentropic pressure ratio P0/P
function isen_P_ratio(M, g) result(pr)
double precision, intent(in) :: M, g
double precision :: pr
pr = (1.0d0 + (g-1.0d0)/2.0d0 * M**2)**(g/(g-1.0d0))
end function
! Find Mach from pressure ratio (subsonic)
function mach_from_pressure(Pb, P0_val, g) result(M_out)
double precision, intent(in) :: Pb, P0_val, g
double precision :: M_out
M_out = sqrt(2.0d0/(g-1.0d0) * ((P0_val/Pb)**((g-1.0d0)/g) - 1.0d0))
end function
! Find subsonic Mach from area ratio (bisection)
function mach_from_area_sub(ar_target, g) result(M_out)
double precision, intent(in) :: ar_target, g
double precision :: M_out, ml, mh, mm, arm
integer :: jj
ml = 0.001d0; mh = 1.0d0
do jj = 1, 200
mm = (ml + mh) / 2.0d0
arm = area_mach_ratio(mm, g)
if (arm > ar_target) then
ml = mm
else
mh = mm
end if
end do
M_out = mm
end function
! Find supersonic Mach from area ratio (bisection)
function mach_from_area_sup(ar_target, g) result(M_out)
double precision, intent(in) :: ar_target, g
double precision :: M_out, ml, mh, mm, arm
integer :: jj
ml = 1.0d0; mh = 20.0d0
do jj = 1, 200
mm = (ml + mh) / 2.0d0
arm = area_mach_ratio(mm, g)
if (arm < ar_target) then
ml = mm
else
mh = mm
end if
end do
M_out = mm
end function
! Find shock location in diverging section
subroutine find_shock_location(AR_exit, Pb, P0_val, g, A_t_val, A_sh, AR_sh, M_up, M_dn)
double precision, intent(in) :: AR_exit, Pb, P0_val, g, A_t_val
double precision, intent(out) :: A_sh, AR_sh, M_up, M_dn
double precision :: ar_lo, ar_hi, ar_mid
double precision :: M1, M2, P02_01, P0_2, AR_downstream
double precision :: P_exit_test, M_exit_test
integer :: kk
ar_lo = 1.0d0
ar_hi = AR_exit
do kk = 1, 200
ar_mid = (ar_lo + ar_hi) / 2.0d0
! Supersonic Mach at shock location
M1 = mach_from_area_sup(ar_mid, g)
! Post-shock Mach
M2 = sqrt(((g-1.0d0)*M1**2 + 2.0d0) / (2.0d0*g*M1**2 - (g-1.0d0)))
! Stagnation pressure ratio across shock
P02_01 = (((g+1.0d0)*M1**2)/((g-1.0d0)*M1**2+2.0d0))**(g/(g-1.0d0)) * &
(1.0d0/(1.0d0+2.0d0*g/(g+1.0d0)*(M1**2-1.0d0)))**(1.0d0/(g-1.0d0))
P0_2 = P0_val * P02_01
! Post-shock A/A* for M2
AR_downstream = area_mach_ratio(M2, g)
! Exit ratio: A_exit/A*_new = (AR_exit/ar_mid) * AR_downstream
M_exit_test = mach_from_area_sub((AR_exit/ar_mid) * AR_downstream, g)
P_exit_test = P0_2 / isen_P_ratio(M_exit_test, g)
if (P_exit_test > Pb) then
ar_lo = ar_mid
else
ar_hi = ar_mid
end if
end do
AR_sh = ar_mid
A_sh = ar_mid * A_t_val
M_up = mach_from_area_sup(ar_mid, g)
M_dn = sqrt(((g-1.0d0)*M_up**2 + 2.0d0) / (2.0d0*g*M_up**2 - (g-1.0d0)))
end subroutine
end program nozzle_flow
Solver Description
Analyzes quasi-1D steady compressible flow through convergent-divergent nozzles. Finds critical pressures, evaluates if throat is choked ($M=1$), and computes exit Mach number ($M_e$) and nozzle state (underexpanded, overexpanded, normal shock).
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 -ffree-line-length-none nozzle_flow.f90 -o nozzle_flow_calc
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
nozzle_flow_calc < input.txt
๐ฅ Downloads & Local Files
Preview of the required input file (input.txt):
! Stagnation Temp ($T_0$) [K]
300.0
! Stagnation Pressure ($P_0$) [kPa]
100.0
! Throat Area ($A^*$) [cmยฒ]
10.0
! Exit Area ($A_e$) [cmยฒ]
20.0
! Specific Heat Ratio ($\gamma$)
1.4
! Back Pressure ($P_b$) [kPa]
80.0
300.0
! Stagnation Pressure ($P_0$) [kPa]
100.0
! Throat Area ($A^*$) [cmยฒ]
10.0
! Exit Area ($A_e$) [cmยฒ]
20.0
! Specific Heat Ratio ($\gamma$)
1.4
! Back Pressure ($P_b$) [kPa]
80.0