🔬
Solver Purpose & Physical Scope
Size converging-diverging nozzles using isentropic compressible relations. Compute throat/exit areas, expansion ratio, thrust, and specific impulse.
📂 Discipline: Fluid
⚡ Precision: IEEE-754 64-bit Real(`real(8)`)
📥 Total Downloads: 397 times
📄 Source File:
rocket_nozzle.f90
📁 calcul/FluidMechanics /
rocket_nozzle.f90
program rocket_nozzle
implicit none
integer :: i, iostat_val
double precision :: Pc, Tc, gamma, R_gas, M_mol, Pe, Pa
double precision :: At, Ae, Dt, De, Me, Te, Ve, rho_e, mdot
double precision :: CF, Isp, F_thrust, c_star, expansion_ratio
double precision :: Pt_ratio, Tt_ratio, rho_ratio, a_ratio
double precision :: M_i, AR_i, P_i, T_i, V_i, F_i, CF_i
double precision :: R_specific, cp, a_throat, rho_c, V_throat
double precision, parameter :: PI = 3.141592653589793d0
double precision, parameter :: g0 = 9.80665d0
read(*,*,iostat=iostat_val) Pc
if (iostat_val /= 0) then
write(*,*) 'ERROR: Invalid chamber pressure input.'
stop
end if
read(*,*,iostat=iostat_val) Tc
read(*,*,iostat=iostat_val) gamma
read(*,*,iostat=iostat_val) M_mol
read(*,*,iostat=iostat_val) Pe
read(*,*,iostat=iostat_val) Pa
read(*,*,iostat=iostat_val) Dt
if (iostat_val /= 0) then
write(*,*) 'ERROR: Failed to read all nozzle inputs.'
stop
end if
if (Pc <= 0.0d0 .or. Tc <= 0.0d0 .or. gamma <= 1.0d0) then
write(*,*) 'ERROR: Pc, Tc must be positive and gamma > 1.'
stop
end if
if (M_mol <= 0.0d0) then
write(*,*) 'ERROR: Molar mass must be positive.'
stop
end if
if (Dt <= 0.0d0) then
write(*,*) 'ERROR: Throat diameter must be positive.'
stop
end if
if (Pe <= 0.0d0) Pe = 101325.0d0
if (Pa < 0.0d0) Pa = 0.0d0
R_specific = 8314.46d0 / M_mol ! J/(kg K)
cp = gamma * R_specific / (gamma - 1.0d0)
At = PI / 4.0d0 * Dt**2
! Characteristic velocity c*
c_star = sqrt(gamma * R_specific * Tc) / gamma / &
sqrt((2.0d0/(gamma+1.0d0))**((gamma+1.0d0)/(gamma-1.0d0)))
! Mass flow rate (choked throat)
rho_c = Pc / (R_specific * Tc)
a_throat = sqrt(gamma * R_specific * Tc * &
(2.0d0/(gamma+1.0d0)))
mdot = Pc * At / c_star
! Exit Mach number from pressure ratio Pe/Pc
Me = exit_mach(Pc, Pe, gamma)
! Expansion ratio Ae/At from exit Mach
expansion_ratio = area_ratio(Me, gamma)
Ae = At * expansion_ratio
De = sqrt(4.0d0 * Ae / PI)
! Exit conditions
Te = Tc / (1.0d0 + (gamma-1.0d0)/2.0d0 * Me**2)
Ve = Me * sqrt(gamma * R_specific * Te)
rho_e = Pe / (R_specific * Te)
! Thrust
F_thrust = mdot * Ve + (Pe - Pa) * Ae
! Thrust coefficient CF = F / (Pc At)
CF = F_thrust / (Pc * At)
! Specific impulse
Isp = F_thrust / (mdot * g0)
! Isentropic ratios at throat (M=1)
Pt_ratio = (1.0d0 + (gamma-1.0d0)/2.0d0)**(-gamma/(gamma-1.0d0))
Tt_ratio = 2.0d0 / (gamma + 1.0d0)
V_throat = sqrt(gamma * R_specific * Tc * Tt_ratio)
write(*,'(A)') '============================================================'
write(*,'(A)') ' ROCKET NOZZLE DESIGN (DE LAVAL) ENGINE'
write(*,'(A)') '============================================================'
write(*,*)
write(*,'(A)') '--- INPUTS --------------------------------------------------'
write(*,'(A,ES12.4,A)') ' Chamber Pressure Pc = ', Pc, ' Pa'
write(*,'(A,ES12.4,A)') ' Chamber Temperature Tc = ', Tc, ' K'
write(*,'(A,ES12.4)') ' Specific Heat Ratio gamma = ', gamma
write(*,'(A,ES12.4,A)') ' Molar Mass = ', M_mol, ' kg/kmol'
write(*,'(A,ES12.4,A)') ' R_specific = ', R_specific, ' J/(kg.K)'
write(*,'(A,ES12.4,A)') ' Exit Pressure Pe = ', Pe, ' Pa'
write(*,'(A,ES12.4,A)') ' Ambient Pressure Pa = ', Pa, ' Pa'
write(*,'(A,ES12.4,A)') ' Throat Diameter Dt = ', Dt, ' m'
write(*,*)
write(*,'(A)') '--- NOZZLE GEOMETRY -----------------------------------------'
write(*,'(A,ES12.4,A)') ' Throat Area At = ', At, ' m2'
write(*,'(A,ES12.4,A)') ' Exit Area Ae = ', Ae, ' m2'
write(*,'(A,ES12.4,A)') ' Exit Diameter De = ', De, ' m'
write(*,'(A,ES12.4)') ' Expansion Ratio Ae/At = ', expansion_ratio
write(*,*)
write(*,'(A)') '--- EXIT CONDITIONS -----------------------------------------'
write(*,'(A,ES12.4)') ' Exit Mach Number Me = ', Me
write(*,'(A,ES12.4,A)') ' Exit Temperature Te = ', Te, ' K'
write(*,'(A,ES12.4,A)') ' Exit Velocity Ve = ', Ve, ' m/s'
write(*,'(A,ES12.4,A)') ' Exit Density = ', rho_e, ' kg/m3'
write(*,*)
write(*,'(A)') '--- PERFORMANCE ---------------------------------------------'
write(*,'(A,ES12.4,A)') ' Mass Flow Rate = ', mdot, ' kg/s'
write(*,'(A,ES12.4,A)') ' Characteristic Velocity = ', c_star, ' m/s'
write(*,'(A,ES12.4,A)') ' Thrust = ', F_thrust, ' N'
write(*,'(A,ES12.4,A)') ' Thrust = ', F_thrust/1000.0d0, ' kN'
write(*,'(A,ES12.4)') ' Thrust Coefficient CF = ', CF
write(*,'(A,ES12.4,A)') ' Specific Impulse Isp = ', Isp, ' s'
write(*,'(A,ES12.4,A)') ' Effective Exhaust Vel = ', Isp*g0, ' m/s'
write(*,*)
write(*,'(A)') '--- THROAT CONDITIONS (M=1) ---------------------------------'
write(*,'(A,ES12.4)') ' P_throat / Pc = ', Pt_ratio
write(*,'(A,ES12.4)') ' T_throat / Tc = ', Tt_ratio
write(*,'(A,ES12.4,A)') ' Throat Velocity = ', V_throat, ' m/s'
write(*,*)
! Nozzle profile: Mach, P, T, V vs A/At
write(*,'(A)') '--- NOZZLE PROFILE ------------------------------------------'
write(*,'(A)') ' A/At Mach P[Pa] T[K] V[m/s]'
write(*,'(A)') ' ----------------------------------------------------------------'
do i = 1, 80
AR_i = 1.0d0 + (expansion_ratio - 1.0d0) * (dble(i-1)/79.0d0)**1.5
M_i = mach_from_area_ratio(AR_i, gamma)
T_i = Tc / (1.0d0 + (gamma-1.0d0)/2.0d0 * M_i**2)
P_i = Pc * (T_i/Tc)**(gamma/(gamma-1.0d0))
V_i = M_i * sqrt(gamma * R_specific * T_i)
write(*,'(F10.4,2X,F10.5,2X,ES12.4,2X,F10.2,2X,F12.2)') &
AR_i, M_i, P_i, T_i, V_i
end do
write(*,*)
! Thrust vs expansion ratio sweep
write(*,'(A)') '--- THRUST VS EXPANSION RATIO SWEEP -------------------------'
write(*,'(A)') ' Ae/At Me CF F[kN] Isp[s]'
write(*,'(A)') ' ----------------------------------------------------------------'
do i = 1, 50
AR_i = 1.0d0 + 99.0d0 * dble(i-1) / 49.0d0
M_i = mach_from_area_ratio(AR_i, gamma)
T_i = Tc / (1.0d0 + (gamma-1.0d0)/2.0d0 * M_i**2)
P_i = Pc * (T_i/Tc)**(gamma/(gamma-1.0d0))
V_i = M_i * sqrt(gamma * R_specific * T_i)
F_i = mdot * V_i + (P_i - Pa) * At * AR_i
CF_i = F_i / (Pc * At)
write(*,'(F10.2,2X,F10.5,2X,F10.5,2X,F12.4,2X,F10.2)') &
AR_i, M_i, CF_i, F_i/1000.0d0, F_i/(mdot*g0)
end do
write(*,*)
write(*,'(A)') '--- CORRELATIONS USED ---------------------------------------'
write(*,'(A)') ' Isentropic flow: T/Tc = 1/(1+(gamma-1)/2 M^2).'
write(*,'(A)') ' Area-Mach: A/At = (1/M)((2+(gamma-1)M^2)/(gamma+1))^((gamma+1)/(2(gamma-1))).'
write(*,'(A)') ' Thrust: F = mdot Ve + (Pe-Pa) Ae.'
write(*,'(A)') ' c* = sqrt(gamma R Tc)/gamma / sqrt((2/(gamma+1))^((gamma+1)/(gamma-1))).'
write(*,'(A)') ' CF = F/(Pc At); Isp = F/(mdot g0).'
contains
double precision function exit_mach(Pcin, Pein, gam)
implicit none
double precision, intent(in) :: Pcin, Pein, gam
double precision :: ratio, lo, hi, mid, P_mid
integer :: it
ratio = Pein / Pcin
if (ratio >= 1.0d0) then
exit_mach = 0.0d0
return
end if
lo = 1.0d0
hi = 50.0d0
do it = 1, 200
mid = 0.5d0*(lo+hi)
P_mid = (1.0d0 + (gam-1.0d0)/2.0d0*mid**2)**(- gam/(gam-1.0d0))
if (P_mid > ratio) then
lo = mid
else
hi = mid
end if
end do
exit_mach = 0.5d0*(lo+hi)
end function exit_mach
double precision function area_ratio(Min, gam)
implicit none
double precision, intent(in) :: Min, gam
double precision :: term
if (Min < 1.0d-10) then
area_ratio = 1.0d30
return
end if
term = (2.0d0 + (gam-1.0d0)*Min**2) / (gam+1.0d0)
area_ratio = (1.0d0/Min) * term**((gam+1.0d0)/(2.0d0*(gam-1.0d0)))
end function area_ratio
double precision function mach_from_area_ratio(AR_in, gam)
implicit none
double precision, intent(in) :: AR_in, gam
double precision :: lo, hi, mid, ar_mid
integer :: it
! Supersonic branch
lo = 1.0d0
hi = 50.0d0
do it = 1, 200
mid = 0.5d0*(lo+hi)
ar_mid = area_ratio(mid, gam)
if (ar_mid < AR_in) then
lo = mid
else
hi = mid
end if
end do
mach_from_area_ratio = 0.5d0*(lo+hi)
end function mach_from_area_ratio
end program rocket_nozzle
💻 How to Compile & Run Locally
1. Compilation (GNU Fortran / Intel oneAPI):
gfortran -O3 rocket_nozzle.f90 -o rocket_nozzle
2. Execution with input.txt redirection:
rocket_nozzle < input.txt
📄 Sample input.txt File Structure
Sample Data:
0.0 0.0 0.0 0.0 0.0 0.0 0.0
Parameter Description:
Chamber pressure Pc\nChamber temperature Tc\nSpecific heat ratio γ\nMolar mass [kg/kmol]\nExit pressure Pe\nAmbient pressure Pa\nThroat diameter Dt