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