💻 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.
Rocket Nozzle Design (De Laval)
Core Numerical Engine in Fortran 90 • 24 total downloads
rocket_nozzle.f90
! =========================================================================
! Source File: 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
Solver Description
Size converging-diverging nozzles using isentropic compressible relations. Compute throat/exit areas, expansion ratio, thrust, and specific impulse.
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 rocket_nozzle.f90 -o rocket_nozzle
Execution Command:
Execute the program by feeding the sample input file into the program using stdin redirection:
rocket_nozzle < input.txt
📥 Downloads & Local Files
Preview of the required input file (input.txt):
! Chamber pressure Pc\nChamber temperature Tc\nSpecific heat ratio γ\nMolar mass [kg/kmol]\nExit pressure Pe\nAmbient pressure Pa\nThroat diameter Dt
0.0
! Parameter 2
0.0
! Parameter 3
0.0
! Parameter 4
0.0
! Parameter 5
0.0
! Parameter 6
0.0
! Parameter 7
0.0
0.0
! Parameter 2
0.0
! Parameter 3
0.0
! Parameter 4
0.0
! Parameter 5
0.0
! Parameter 6
0.0
! Parameter 7
0.0