4.2 Characteristics in the internal gas dynamics

--- tags: Mathematical modeling and optimization title: 4.2 Characteristics in the internal gas dynamics --- ## 4.2. Characteristics in the internal gas dynamics Gas dyanmics refers to the discipline of compressible flows, which describes the behavior of gases when the density changes significantly in response to pressure variations. This is a common occurrence in e.g. propulsion systems, gas pipelines and pneumatic applications. Unlike incompressible flow, where the density is constant, compressible flow requires accounting for changes in density, pressure, and temperature. [__Speed of sound__](https://en.wikipedia.org/wiki/Speed_of_sound) ($a$) is a specific physical property in the compressible flow scop, which denotes the pressure change over density change in the isentropic condition, derived from the continuety and momentum conservation over a wave front. This yields to $$ c^2 = \left(\frac{\partial p}{\partial \rho}\right)_s $$ In the combination of isentropic process ($\ p v^\gamma = \text{const}$) and ideal gas law ($\ p = \rho \ R\ T$), speed of sound has a simple expression as $c = \sqrt{\gamma \ R\ T}$. The speed of sound concept allows the Mach number ($M$) to quantify its ratio to the gas (or related to gas) velocity, written as $$ M = \frac{u}{c}\quad . $$ [__Stagnation state__](https://en.wikipedia.org/wiki/Stagnation_temperature) (subscript $_0$) is another concept which well applied in the regime of gas dynamics. The stagnation state denotes the condition where no dynamics (velocity) is taken into account. Accordingly, the condition with dynamic expression is name as _static_ state. Thus, $$ \begin{align} \underbrace{h_0}_{\text{stagnation enthalpy}} = \underbrace{h}_{\text{static enthalpy}} +\underbrace{\frac{1}{2}\ \ u^2}_{\text{kinematic energy}} \end{align} $$ Involving the Mach number and speed of sound definition, the enthalpy can be represented in another form: $$ h_0 = h + \frac{\gamma -1}{2}M^2 \ c_p \ T = (1 + \frac{\gamma -1}{2} M^2 ) \ h $$ and further reduced due to calorific perfect gas assumption $(\ h = c_p \ T\ )$ to $$ \frac{T_0}{T} = (1 + \frac{\gamma -1}{2} M^2 ) $$ where $\gamma$ is the ratio of specific heat capacity ${c_p}/{c_v}$. This relation can be interpreted as: _for a given stagnation temperature $(T_0)$, the static temperature $(T)$ on any point is in a simple behavior to its Mach number incorporated with its themrmodynamic property $\gamma$ of the gas._ ### 4.2.1 Isentropic flow and de Laval nozzle The concept of isentropic process $(pv^\gamma = \text{const})$ futher brings up the $M$ representations between static and stagnation state to other phyiscal properties. Involving the ideal gas law, the pressure and density relation to stagnation states yield $$ \begin{align} \frac{p_0}{p} = \left(\frac{T_0}{T}\right)^{\gamma / (\gamma -1)}= \left(1 + \frac{\gamma -1}{2} M^2 \right)^{\gamma / (\gamma -1)}\\ \\ \frac{\rho_0}{\rho} = \left(\frac{T_0}{T}\right)^{1 / (\gamma -1)}= \left(1 + \frac{\gamma -1}{2} M^2 \right)^{1 / (\gamma -1)}\\ \end{align} $$ The starting points for internal compressible flow are the continuity and momentum equation for one-dimenional state. They present in the form of mass flow rate and Euler's equation, where $$ \begin{align} \dot{m} = \rho \ u\ A\ = & \ \text{const},\\ \\ \frac{1}{\rho} \ dp + \frac{1}{2} d u^2 \ &= \ 0 \end{align} $$ where, - $\rho$ is the gas density - $u$ is the flow velocity, and - $A$ is the cross section area The further handling takes the derivative on both sides, the governing equations become $$ \begin{align} d(\rho uA) =\ Au\ d\rho + \rho u \ dA+ \rho A \ du &= d(\text{const}) = \ 0 \\ \to \ \frac{d\rho}{\rho} + \frac{dA}{A} + \frac{du}{u} &= 0 \\ \\ d(\rho \ u^2 + p )= \ d(\text{const}) & = 0\\ \to du = -\frac{dp}{\rho\ u} \end{align} $$ In alignment to the speed of sound definition, the cominition of continuity-momentum equation results to $$ \begin{align} \frac{d\rho}{\rho}\ =& \ \frac{M^2}{1 - M^2} \frac{dA}{A}.\\ \\ \frac{du}{u}\ =& -\frac{1}{1-M^2} \frac{dA}{A}\\ \\ dp\ = & \ \rho u^2 \frac{1}{1-M^2} \frac{dA}{A} \end{align} $$ These relation provide the baseline rules for analysis in the 1-D dominant flow patterns, such as - For subsonic state $(M < 1)$, __decreasing of cross section area $(dA < 0)$__ leads to - increase of velocity $(du >0)$ and Mach number, - decrease of pressure $(dp < 0)$, - decrease of density $(d \rho < 0)$ and - decrease of temperature $(dT < 0)$, vice versa. - For supersonic state $(M > 1)$, __increase of cross section area $(dA > 0)$__ leads to - increase of velocity $(du > 0)$ and Mach number, - decrease of pressure $(dp < 0)$, - decrease of density $(d \rho < 0)$ and - decrease of temperature $(dT < 0)$, vice versa. It is found that the velocity reaches its maximum for a converging pipe at speed of sound, where $M = 1$. A further speed up can only be achieved in a diverging shape $(dA > 0)$. This is a common strategy to gain the kinetic energy out of the nozzle, well applied in rocket design. #### Sonic (critical) state when $M = 1$, flow _chokes_, denoting that the maximal mass flow is reached. The choke state is usually symbolized with the superscript $(^*)$. The physical relation to the stagnation state read $$ \begin{align} \frac{T^*}{T_0}\ =& \ \frac{2}{\gamma +1} \\ \frac{p^*}{p_0}\ =& \ \left(\frac{2}{\gamma +1} \right)^{\gamma / (\gamma -1)} \\ \frac{\rho^*}{\rho_0}\ =& \ \left(\frac{2}{\gamma +1} \right)^{1 / (\gamma -1)} \end{align} $$ and further, the throat area relation reveals the Mach number at the end of the nozzle, which $$ \frac{A}{A^*} \ = \frac{1}{M}\ \left( \frac{2}{\gamma +1}\left(1+\frac{\gamma -1}{2}M^2\right)\right)^{\frac{\gamma+1}{2(\gamma-1)}} $$ #### [De Laval (Convergent-Divergent) Nozzle](https://en.wikipedia.org/wiki/De_Laval_nozzle) <img style="float: right;" src="https://upload.wikimedia.org/wikipedia/commons/thumb/4/4c/Nozzle_de_Laval_diagram.svg/800px-Nozzle_de_Laval_diagram.svg.png" width = 40%> Based on this understanding, convergent-divergent nozzle (also called de Laval nozzle) is created to generate internal flows in supersonic state. The behavior of physical properties within the nozzle can be then described exactly based on the shape of nozzle, i.e. the nozzle cross area. During the whole process in the convergent-divergent nozzle, the pressure is descreasing. If the pressure level on the exit of nozzle does not match the environment pressure, [shock](https://en.wikipedia.org/wiki/Thermodynamic_relations_across_normal_shocks) occurs. Shock can be regarded as the non-kontineous interface between two pressure levels. The shock front can be in the divergent section, can also be outside of the nozzle. Physics and effects of the normal shock in or outside of the internal flow are one of the central topics inw the discipline of gas dynamics. Nevertheless, the pressure relation provides a usful function to determine the nozzle throat state and further evaluate the mass flow rate over the nozzle system. <img src="https://media.springernature.com/m685/springer-static/image/art%3A10.1038%2Fs41598-024-53680-2/MediaObjects/41598_2024_53680_Fig1_HTML.png" width="95%"> #### Convergent nozzle and back pressure impact The decreasing pressure also reveals significant technical insights. Say, if the pressure level outside of the nozzle ($p_b$, also called back pressure) is less than the critical pressure ($p^*$), the flow chokes ($M=1$) on the nozzle throat, e.g. state C, D, E, F on the upper diagram. Cutting the convergent-divergent nozzle on the nozzle troat, the situation becomes easier. The 6 scenarios are reduced to 3. For back pressure higher than scenario 3, i.e. $(p_b > p^*)$ scenario A, B, the Mach number on the nozzle throat, or better said nozzle exit (with subscript "$_e$"), $M_e$ can be calculated using the stagnation relation $$ \frac{p_0}{p_b} = \left(1 + \frac{\gamma -1}{2} M_e^2 \right)^{\gamma / (\gamma -1)}\qquad.\\ $$ For back pressure lower than scenario 3 $(p_b < p^*)$, $M_e = 1$. [Expansion waves](https://aerospaceweb.org/question/propulsion/q0224.shtml) occurs behind the nozzle exit. Summarizing the two scenarios, an unified formulation yields: $$ M_e = \text{min}\left\{ \ \left( \frac{2}{\gamma - 1} \right)^{\frac{1}{2}}\ \left[ \left( \frac{p_t}{p_b} \right)^{\frac{\gamma -1}{\gamma }} -1 \right]^{\frac{1}{2}}, \ 1.0 \right\} $$  The convergent nozzle can also be regarded as valve in the technical applications. By using this formulation, the mass flow releasing from a system can be described, as $$ \dot{m} = \rho_e \ A_e \ u_e= \rho_t \ A_e \ M_e \ c_e \qquad . $$ This is the counterpart description for compressible gas to the Hagen-Poiseuelli equation / Bernoulli equation of the incompressinle fluid. #### Scenario -- Cartridge Gas releasing mechanism Combining the energy representation within a gas system with the mass/energy releasing throughout the nozzle, we can describe a dynamic process by focuing the system as a controll volume. The mass and energy flow rate out of the nozzle/valve are the loss terms in the mass and energy conservation. <img src="https://live.staticflickr.com/65535/53904972165_23096376f7.jpg" width = 70%> $$ \begin{align} \frac{d}{d t} \ m (t) &= -\dot{m}_e(t) \qquad &\text{Mass conservation}\\ \\ \frac{d}{d t}\ \left( c_v\ m(t) \ T(t)\right) &= - \dot{m}_e c_p T_e(t) - \frac{1}{2} \dot{m}_e u_e^2 -\dot{Q}(T, t)\qquad &\text{Energy conservation}\\ \\ P(t) &= \rho(t) \ R \ T(t) \qquad &\text{Equation of state}\\ \\ \dot{m}_e &= \rho_e \ c_e\ M_e \ A \end{align} $$ with $\dot{m}_e$, $T_e$, and $u_e$ as the mass flow rate, temperature and the velocity outwards of the nozzle exit, as the $Q$ as the heat loss of the system. (mostly due to the heat transfer.) One can also involve the temperature dependent specific heat capacity ($c_p(T)$ and $c_v(T)$) towards the real gas correction. A daily example is the gas filling tank, where pressurized gas is filled in a container. Once the valve is opened, the high pressurized gas flees due to the pressure drop. We applied this mechanism to pump up the tire or inflate a ballon. Airbag also rely on such mechanism, the stored gas inflators are exactly designed in this way. Here is a scenario: _A stored gas inflator cartridge of 150cc volume is filled with Nitrogen of 50bar. The cartridge is equipped with a vent of 1.5mm diameter._ _Once the vent is opened,_ _1. how long does it take to release the stored gas completely? 2. what is the temperature level in the bottle at the end of inflation? 3. What if the gas is changed to He?_ A python implementation for this reads ``` #! PYTHON v3.7 SCRIPT from numpy import pi, sqrt, square from scipy.optimize import newton, fsolve import matplotlib.pyplot as plt #-------------------------------------------------------------------------# # Task: # A cartridge of 150cc volume is filled with Nitrogen of 50bar. The Cartridge is equipped with a vent of 1.5mm diameter. Once the vent is opened, # 1. how long does it take to release stored gas completely? # 2. what is the temperature level in the bottle at the end of inflation? # 3. What if the gas is changed to H2? #-------------------------------------------------------------------------# # Assumption: # - caloric perfect gas, ideal gas law # - constant cp, cv and gamma (no temperature dependency obtained) # - no discharge on vent, cross area perfectly used # - no heat loss considered #><><><><><><><><><><><><><>< Functions ><><><><><><><><><><><><><><# # ----------------------------------------------------------------------- # def calc_cross_area(diameter): area = 1.0/4.0*pi*diameter**2.0 return area # ----------------------------------------------------------------------- # def calc_Mach_number(stagnation_pressure, back_pressure, gamma): if gamma != 0.0 : Mach_formula = sqrt(max(2.0/(gamma -1.0)*((stagnation_pressure/ back_pressure)**((gamma -1.0)/gamma) - 1.0 ),0.0)) Mach = min(1.0, Mach_formula) else: Mach = 0.0 return Mach # ----------------------------------------------------------------------- # def calc_nozzle_properties(stagnation_pressure, stagnation_temerature, stagnation_density, stagnation_gamma, back_pressure): Mach = calc_Mach_number(stagnation_pressure, back_pressure, stagnation_gamma) temperature = stagnation_temerature/(((stagnation_gamma-1.0)*Mach**2.0)/2.0+1.0) density = stagnation_density/(((stagnation_gamma-1.0)*Mach**2.0)/2.0+1.0)**(1.0/(stagnation_gamma-1.0)) pressure = stagnation_pressure /(((stagnation_gamma-1.0)*Mach**2.0)/2.0+1.0)**(stagnation_gamma/(stagnation_gamma-1.0)) gamma = stagnation_gamma speedOfSound = sqrt(gamma* Ru/molar_weight* temperature) velocity = speedOfSound*Mach return Mach, temperature, density, pressure, velocity # ----------------------------------------------------------------------- # # ------------------- Equation solver ----------------# # ----------------------------------------------------------------------- # def solv_system_Temperature(Temp_system , internalEnergy_previous, energy_supply, mass_system): internal_energy_system = mass_system * cv * Temp_system error = internal_energy_system - (internalEnergy_previous - energy_supply) return (error) #><><><><><><><><><><><><><>< Gas Properties ><><><><><><><><><><><><><><# Ru = 8.314 #J/(K mole) universal gas constant #====================== Nitrogen =======================# molar_weight = 28*1e-3 #[kg] cp = 1040 #J/(kg K) ''' #====================== CO2 =======================# molar_weight = 44*1e-3 #[kg] cp = 846 #J/(kg K) ''' R = Ru/molar_weight #=======================================================================# cv = cp - Ru/molar_weight gamma = cp/cv #><><><><><><><><><><><><><>< Physical properties ><><><><><><><><><><><><><><# pressure_environment = 1.013*1e5 #[Pa] temperature_environment = 23+273 #[K] volume_system = 150*1e-6 #[m^3] pressure_system = 30*1e5 #[Pa] temperature_system = temperature_environment diameter_nozzle = 1.0*1e-3 #[m] discharge = 0.5 #[-] area_nozzle = discharge* calc_cross_area(diameter_nozzle) #><><><><><><><><><><><><><>< numerical properties ><><><><><><><><><><><><><><# ent_time = 10.0 #[s] timeStep = 1e-5 #[s] outputInterval =1e-2 #[s] status = "go" iterationNumber = int(ent_time/timeStep) outputIteration = int(outputInterval/timeStep) #%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Main %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%# #---------------------- initial condition ------------------------# pressure_initial = pressure_system temperature_initial = temperature_environment mass_initial = pressure_initial*volume_system*molar_weight/(Ru * temperature_initial) print(f"initial gas weight: {mass_initial*1e3:.3f} g") #---------------------- creating output file -------------------------# IO_time = [] IO_pressure = [] IO_temperature = [] IO_mass = [] IO_massFlowRate = [] IO_Mach = [] for i in range(0, iterationNumber): currentTime = i * timeStep density_initial = mass_initial/volume_system internal_energy_initial = mass_initial * cv * temperature_initial #--------------------- calculate exit mass and energy ------------------------# Mach_exit, temperature_exit, density_exit, pressure_exit, velocity_exit = \ calc_nozzle_properties(pressure_initial, temperature_initial, density_initial, gamma, pressure_environment) massFlowRate_exit = velocity_exit* density_exit* area_nozzle massAmount_exit = massFlowRate_exit* timeStep enthalpy_exit = 0.5 * massAmount_exit * square(velocity_exit) +\ massAmount_exit * cp * temperature_exit mass_initial = mass_initial - massAmount_exit temperature_system = newton(solv_system_Temperature, temperature_initial, maxiter = 300, tol=1.0e-10, \ args=(internal_energy_initial, enthalpy_exit, mass_initial)) #temperature_system = (internal_energy_initial - enthalpy_exit)/(cv*mass_initial) pressure_system = mass_initial/volume_system * R * temperature_system pressure_initial = pressure_system temperature_initial = temperature_system if pressure_system <= pressure_environment and status == "go" : status = "stop" print(f"complete released at {currentTime:.3f} s") print(f"temperature: {temperature_system-273:.3f} °C") break if i % outputIteration == 0 : IO_time.append(currentTime) IO_pressure.append((pressure_system-101300)/1e5) IO_temperature.append(temperature_initial-273) IO_mass.append(mass_initial*1e3) IO_massFlowRate.append(massFlowRate_exit) IO_Mach.append(Mach_exit) fig, (ax1, ax2) = plt.subplots(2) fig.suptitle('Gas release characteristics') ax1.plot(IO_time, IO_pressure) ax1.grid(True) ax1.set_ylabel("pressure [bar]") ax2.plot(IO_time, IO_temperature) ax2.grid(True) ax2.set_xlabel("time [ms]") ax2.set_ylabel("temperature [°C]") ax3.plot(IO_time, IO_Mach) ax3.grid(True) ax3.set_xlabel("time [ms]") ax3.set_ylabel("Mach number [-]") fig.show() ``` ### 4.2.2 Solid propellant deflagration and interior ballistics Interior ballistics is a discipline studying the dynamics of pressure condition in one or multiple closed systems, which further utilized to e.g. trigger projectiles in a gun/barrel or serve as propulsion system for rocket movement. <img src="https://d3i71xaburhd42.cloudfront.net/04e2d66d9d75b2d8575ee35bfb323710a2486888/14-Figure1.4-1.png" width=55%> <img src="https://www.researchgate.net/publication/348035213/figure/fig1/AS:11431281090011087@1665852823265/Basic-components-and-structure-of-the-solid-propellant-rocket-engine.jpg" width=35%> Imagine the high velocity of a projectile or the movement of a heavy rocket, all these applications require very high pressure level to push the objectives into demanded dynamics. The pressure level are commonly achieved by deflagration of propellent, which creates huge amount of gas and heat very rapidly (in millisecond scale). These mechanisms involve the chemical reactions, i.e. the fuel-oxidizer behavior, the propellant shape and the combination chamber design. #### Solid propellant deflagration The word _deflagration_ denotes rapid combustion, which requires the homogeneous mixing of fuel and oxidizers. If one thinks on the combustion process in e.g. in the internal combustion engine, the oxidizer is always taken from the environment, i.e. air. Differ to that, the characteristics of solid proepellant is the premixing of fuel and oxidizer in the solid form. Thus, the solid propellant devices can be deployed in deep see, as well as in outer space. The chemical reaction process during the deflagration process is very complex. It involves several intermediate steps to reach the equilibrium state. Fortunately, there are [_database_](http://www.jpyro.co.uk/wp-content/uploads/j01_11_htfsr.pdf) and [_calculatuon tools_](https://rocketworkbench.sourceforge.net/equil.phtml) for common composites. An example is shown in the figure below for balck powder composition. <img src="https://live.staticflickr.com/65535/53906981904_e76f0709fa_c.jpg"> Purposing in quantifying the system performance, the most important information are the _gas yield, species composition_ and the _heat of combustion_ at the final equilibrium state. These are the representative information serve as the source terms in the mass and energy conservation. _For the black powder example, the CO, CO2 and N2 are the main products obtaining ~70% gas yield, shown in the Figure. The combustion enthalpy is calculated as 3494 kJ/kg, as [exothermic process](https://en.wikipedia.org/wiki/Exothermic_reaction)._ #### Burn rate and Vieille's law The combustion of solid propellant always occurs on the surface of the solid body, stated in [__Piobert's law__](https://en.wikipedia.org/wiki/Piobert%27s_law). The propellant design is all based on this statement. Shape of the propellant affects to different burning dynamics, which supports diverse applications. <img src="https://www.nakka-rocketry.net/th_pix/grains1.gif"> [__Vieille's law__](https://www.nakka-rocketry.net/burnrate.html) is another puzzle to complete the mass and energy production rate in the transport process. This is a simplified way to determine the propellant burn rate, in the dimension of [mm/s]. As a rule of thumb, the burn rate of the solid propellant is exponentially proportional to the pressure state, which $$ \dot{\gamma} = \ a \ p^n $$ where - $\dot{\gamma}$ is the burn rate, dimensioned in [mm/s] - $p$ is the pressure, usually dimensioned in [MPa] - $a$, $n$ are the empirical material parameters, which are determined by e.g. [closed vessel tests](https://www.sciencedirect.com/science/article/abs/pii/0010218075900899) #### Solid rocket motor (SRM) application Solid rocket motor is one of the classical applications working with solid propellant. For which, PEP calculation tool also provides the rocket nozzle specific results such as nozzle throat behavior, $c^*$ and Isp, etc. <img src="https://live.staticflickr.com/65535/53907271670_68ea6c8f24_z.jpg" width="70%"> [__Thrust__](https://www.grc.nasa.gov/www/k-12/airplane/rockth.html) is the key performance characteristics for rocket relevant applications. The thrust Equation read: $$ F=\dot{m}\ u_e+(p_e−p_a)A_e $$ where - $F$ is the thrust. - $\dot{m}$ is the mass flow rate of the exhaust. - $u_e$, $p_e$ and $A_e$ are the velocity, pressure and area at the nozzle exit. - $P_a$ is the ambient pressure. Based on the thrust, some analysis indices are further defined. - $\pmb{A_e/A_t}$ ratio is the ratio of the exit area to the throat area of the nozzle, acquired by the isentropic behavior as listed in the last subsection. - $\pmb c^*$ is the characteristic exhaust velocity, defined as $$ c^∗=\dot{m}\ p_c\ A_t = \ \sqrt{\frac{R\ T_c}{\gamma}} \left( \frac{2}{\gamma + 1} \right)^{(\gamma + 1)/\left(2(\gamma -1)\right)} $$ where - $A_t$ is the throat Area (the narrowest part of the rocket nozzle, where the cross-sectional area is minimized). This is where the exhaust gases reach their highest velocity (assuming choked flow). - $T_c$ is the temperature in the combustion chamber, which is dominated by the chemical reaction, i.e. combustion temperature. - Specific Impulse [__*Isp*__](https://www.grc.nasa.gov/www/k-12/airplane/specimp.html) is the ratio of _total impulse_ to _total weight_, defined by $$ \textit{Isp}=\frac{\int F \ dt}{g_0 \ \int \dot{m}\ dt}= \frac{F}{\dot{m}\ g_0} $$ with $g_0$ as the gravitational acceleration. #### Numerical Approach Tracing back to describe a solid rocket propulsion system within a solid inner space as control volume. It results to the equation system: <img style="float: right;" src="https://live.staticflickr.com/65535/53907044578_0cd50fa4cc_w.jpg" width ="18%"> $$ \begin{align} &\text{Mass conservation (solid)}:\\ & \qquad \frac{d m_s}{d t} = -\rho_s \dot{\gamma} A(t) \\ \\ &\text{Mass conservation (gas)}:\\ & \qquad \frac{d m_g}{d t} = \frac{d}{dt} (\rho_g V_g) = \rho_s \dot{\gamma} A(t) \ \eta - \dot{m}_e \\ \\ &\text{Energy conservation}:\\ & \qquad \frac{d E_g}{d t} = \frac{d}{dt} (\rho_g V_g e_g) = \rho_s \dot{\gamma} A(t) \ \eta \ h_c - (\dot{m}_e h_e + \frac{1}{2} \dot{m}_e u_e^2) - \dot{Q}\\ \end{align} $$ with - Internal energy: $\displaystyle e_g = \sum_{i}^{n}\ \phi_i\ c_{v_i}\ T$ - Equation of state (ideal gas law): $\displaystyle p = \rho_g \ R \ T$ - Mass flow rate $\displaystyle \dot{m}_e = \rho_t \ c_t M_t \ A ,\qquad c_t = \sqrt{\gamma R T_t}$ - Gas yield ($\eta$), specific heat capacity in constant volume ($c_v$)/ pressure ($c_p$) - Thrust force $F=\dot{m_e}\ u_e+(p_e−p_a)A_e$ Example: _Configuration of an Ariane4 Solid Rocket Motor filled with Ammonium Nitrate (AN) in a hollow column propellant form. Nozzle of the rocket is applied in the expansion ratio Ae/At as 8.3:1 with the exit diameter 0.96 m. The configuration of Ariane 4 and AN properties reads:_ <img style="float: right;" src="https://live.staticflickr.com/65535/53907035113_a0e7a37ff9.jpg" width="70%"> - _Motor chamber inner-diameter: 1.07 [m]_ - _Motor chamber inner-height: 11.32 [m]_ - _Propellant hollow column inner-diameter: 0.718 [m]_ - _Propellant mass: 9640 [kg]_ - _Density AN: 1724 [kg/m^3]_ - _Burn rate AN:_ - _b0_AN = 0.05*25.4_ - _n0_AN = 0.88_ - _bor_AN = 30_ - _b1_AN = 0.32*25.4_ - _b1_AN = 0.33_ - _Gas yield AN: 100% [-]_ - _$c_p$ of AN combustion production: 1760 [J/(kg.K)]_ - _$c_v$ of AN combustion production: 1390 [J/(kg.K)]_ - _molar weight AN combustion production: 22.85 #[g/mole]_ - _combustion enthalpy: 1825 #[kJ/kg]_ _Assuming the combustion chamber has been heat up to 200°C for a proper ignition of AN, and no air portion in the expelling gas considered._ Assignments: - _The $A_e/A_t$ ratio is 8.3, what is the expected Mach number on the exit side?_ - _What is the pressure, temperature and thrust evolution in the combustion chamber during the entire flight time?_ - _How long does the combustion last?_ - _Is the temperature in the chamber identical to the combustion temperature? If not, why is it?_ A python implementation for this system can be: ``` #! PYTHON 3.7 SCRIPT #------------------------------------------------------# # Task 04 --Solid Rocket motor # Configuration of an Arian4 Solid Rocket Motor filled with Ammonium Nitrate in a hollow column propellant form. # Nozzle of the rocket is applied in the expansion ratio Ae/At as 8.3:1 with the exit diameter 0.96 m. # The configuration of Ariane 4 reads: # Motor chmaber inner-diameter: 1.07 #[m] # Motor chamber inner-height: 11.32 #[m] # Propellant hollow colume inner-diameter: 0.718 #[m] # Propellant mass: 9640 #[kg] # Density AN: 1724 [kg/m^3] # Burn rate AN: # a0_AN = 0.05*25.4 # n0_AN = 0.88 # bor_AN = 30 # a1_AN = 0.32*25.4 # n1_AN = 0.33 from numpy import pi, sqrt, square from scipy.optimize import newton, fsolve import time, sys Ru = 8.314 #universal gas constant #================== Air =================== molar_weight_air = 28.964*1e-3 #[kg/mole] cp_air =1000 #[J/kg.K] cv_air = 718 #[J/kg.K] gamma_air = 1.4 #%%%%%%%%%%%% Universal constants %%%%%%%%%%%%% #================== Propellant =================== density_propellant = 1.724*1e3 #[kg/m^3] gas_yield_weight = 1.0 gas_yield_volumetric = 1.0 diameter_outer_propellant = 1.07 #[m] diameter_inner_propellant = 0.718 #[m] height_propellant = 11.32 #[m] burnrate_AN = [1.27, 0.88, 30.0e6, 8.13, 0.33] cp_gas = 1760 #[J/kg.K] cv_gas = 1390 #[J/kg.K] gamma_gas = cp_gas/cv_gas molar_weight_gas = 22.85*1e-3 #[kg/mole] combustion_enthalpy = 1825.4*1e3 #[J/kg] temperature_initial = 200+273 #[K] pressure_initial = 1.597*1e5 #[Pa] ratio_nozzle = 8.3 #[-] crossArea_throat = 0.087 #[m^2] crossArea_exit = crossArea_throat*ratio_nozzle #================== Environment =================== volume_vessel_default = 10.17 #[m^3] temperature_environment = 23+273 #[K] pressure_environment = 1.013*1e5 #[Pa] area_vessel_surface = 38.033 #[m^2] heat_conduct_coefficient = 0#1.5*1e3 #[J/(K m^2)] #================== Numerical Setups ============ startTime = 0 #[s] timeStep = 1e-5 #[s] outputInterval = 0.5 #[s] endTime = 200.0 #[s] numberOfIteration = int((endTime-startTime)/timeStep) ouputIteration = int(outputInterval/timeStep) outputFile = "results.dat" writeOutput = open(outputFile, "w") # %%%%%%%%%%%%%%%%%%%%%%%%%% Functions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # ----------------------------------------------------------------------- # # ------------------- OS environmental ------------------------# # ----------------------------------------------------------------------- # def print_no_newline(string): sys.stdout.write(string) sys.stdout.flush() # ========================================================================# # ----------------------------------------------------------------------- # # --------------------- State functions ------------------------# # ----------------------------------------------------------------------- # def cal_pressure(mass_gas, mass_air, volume_vessel, temperature): # ideal gas following Dalton's law law pressure_air = mass_air/volume_vessel *Ru/molar_weight_air * temperature pressure_gas = mass_gas/volume_vessel *Ru/molar_weight_gas * temperature pressure = pressure_air + pressure_gas return pressure #---------------------------------------------- def cal_temperature(internal_energy, energy_supply, mass_gas, mass_air): # 1st law of thermodynamics temperature = (internal_energy + energy_supply)/(mass_gas* cv_gas + mass_air* cv_air) return temperature #---------------------------------------------- def cal_internalEnergy(mass_gas, mass_air, temperature): internalEnergy = (mass_gas* cv_gas + mass_air* cv_air)* temperature return internalEnergy # ========================================================================# # ----------------------------------------------------------------------- # # --------------------- Geometric functions ------------------------# # ----------------------------------------------------------------------- # def calc_area_type2(diameter, height): area_surface = pi*diameter*height return area_surface #---------------------------------------------- def calc_volume_type2(diameter_inner, diameter_outer, height): volume = pi/4.0*(diameter_outer**2.0 - diameter_inner**2.0) *height return volume #---------------------------------------------- def evol_geometric_type2(cylinder_geometry, burning_velocity, time_step): diameter_inner_before, diameter_outer,height = cylinder_geometry volume_before = calc_volume_type2(diameter_inner_before, diameter_outer, height) surfaceArea_before = calc_area_type2(diameter_inner_before, height) volume_burnt = min(surfaceArea_before * burning_velocity * time_step , volume_before) diameter_inner_after = max(diameter_inner_before + 2.0 * burning_velocity * time_step , 0.0) volume_after = calc_volume_type2(diameter_inner_after, diameter_outer, height) area_after = calc_area_type2(diameter_inner_after, height) cylinder_geometry = diameter_inner_after, diameter_outer, height return (volume_burnt, volume_after, cylinder_geometry) # ========================================================================# # ----------------------------------------------------------------------- # # --------------------- Propellant functions ------------------------# # ----------------------------------------------------------------------- # def calc_burning_rate(pressure, b_0, n_0, border_pressure, b_1, n_1): if pressure <= border_pressure: burning_rate = 1e-3*b_0 * (pressure/1e6)**n_0 #[m/s] else: burning_rate = 1e-3*b_1 * (pressure/1e6)**n_1 #[m/s] return burning_rate # ========================================================================# # ----------------------------------------------------------------------- # # ------------------- Transition properties ----------------# # ----------------------------------------------------------------------- # def calc_throat_Mach_number(stagnation_pressure, back_pressure, gamma): if gamma != 0.0 : Mach_formula = sqrt(max(2.0/(gamma -1.0)*((stagnation_pressure/ back_pressure)**((gamma -1.0)/gamma) - 1.0 ),0.0)) Mach = min(1.0, Mach_formula) else: Mach = 0.0 return Mach # ========================================================================# def calc_troat_properties(stagnation_pressure, stagnation_temperature, stagnation_density, stagnation_gamma, molar_weight, back_pressure): Mach = calc_throat_Mach_number(stagnation_pressure, back_pressure, stagnation_gamma) temperature = stagnation_temperature/(((stagnation_gamma-1.0)*Mach**2.0)/2.0+1.0) density = stagnation_density/(((stagnation_gamma-1.0)*Mach**2.0)/2.0+1.0)**(1.0/(stagnation_gamma-1.0)) pressure = stagnation_pressure /(((stagnation_gamma-1.0)*Mach**2.0)/2.0+1.0)**(stagnation_gamma/(stagnation_gamma-1.0)) speedOfSound = sqrt(stagnation_gamma* Ru/molar_weight* temperature) velocity = speedOfSound*Mach return Mach, temperature, density, pressure, velocity # ========================================================================# def solv_exit_Mach_number(Mach_exit, ratio_area, Mach_throat, gamma): rhs =Mach_exit/ Mach_throat * ((1.0+((gamma - 1.0)/2.0)*Mach_throat**2.0) / (1.0+((gamma - 1.0)/2.0)*Mach_exit**2.0) )**((gamma+1.0)/(2.0*(gamma-1.0))) lhs = 1.0/ratio_area error = lhs - rhs return error # ========================================================================# def calc_exit_velocity(Mach_throat, Mach_exit, temperature_throat, molar_weight, gamma): temperature_exit = temperature_throat * ((1.0+((gamma-1.0)/2.0)*Mach_throat**2.0)/(1.0+((gamma-1.0)/2.0)*Mach_exit**2.0)) speedOfSound_exit = sqrt(gamma* Ru/molar_weight* temperature_exit ) velocity_exit = Mach_exit*speedOfSound_exit return velocity_exit # ========================================================================# def calc_exit_presure(Mach_throat, Mach_exit, pressure_throat, gamma): pressure_exit = pressure_throat * ((1.0+((gamma-1.0)/2.0)*Mach_throat**2.0)/(1.0+((gamma-1.0)/2.0)*Mach_exit**2.0))**(gamma/(gamma-1.0)) return pressure_exit # %%%%%%%%%%%%%%%%%%%%%%%%%% Main %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% # ========================================================================# # ----------------------------------------------------------------------- # # ------------------- initial condition ----------------# # ----------------------------------------------------------------------- # volume_propellant = calc_volume_type2(diameter_inner_propellant, diameter_outer_propellant, height_propellant) volume_vessel_current = volume_vessel_default - volume_propellant temperature_current = temperature_initial pressure_current = pressure_initial density_air = pressure_current/(Ru/molar_weight_air * temperature_current) mass_air =density_air * volume_vessel_current mass_gas = 0.0 burnt_volume = 0.0 writeOutput.write("#time"+" "+"Pressure"+" "+"Temperature"+" "+"diameter"+" "+"height"+" "+"surface_area"+" "+\ "burnt_volume"+" "+"rest_volume"+"\n") writeOutput.write("#[s]"+" "+"[bar]"+" "+"[K]"+" "+"[mm]"+" "+"[mm]"+" "+"[mm^2]"+" "+\ "[mm^3]"+" "+"[mm^3]"+"\n") print("pressure",pressure_current/1e6, "MPa") for i in range(0, numberOfIteration): currentTime = i*timeStep + startTime internalEnergy = cal_internalEnergy(mass_gas, mass_air, temperature_current) gasDensity_current = mass_gas/volume_vessel_current mass_supply = 0.0 energy_supply = 0.0 #================================ transitional properties ==================================# #---------------------------- assuming all the expelled gas are propellant producted -----------# #-------------------------------------------------------------------------------------------------# [Mach_throat, temperature_throat, density_throat, pressure_throat, velocity_throat] = \ calc_troat_properties(pressure_current, temperature_current, gasDensity_current, gamma_gas, molar_weight_gas, pressure_environment) massFlowRate_throat = density_throat * velocity_throat* crossArea_throat massAmount_throat = min(massFlowRate_throat* timeStep, mass_gas) energyAmount_throat = 0.5 * massAmount_throat * square(velocity_throat) +\ massAmount_throat * cp_gas * temperature_throat mass_supply = mass_supply - massAmount_throat energy_supply = energy_supply - energyAmount_throat #--------------------------- exit properties --------------------------------------------------# if Mach_throat < 1.0: Mach_exit = newton(solv_exit_Mach_number, 0.7, args=(ratio_nozzle, Mach_throat, gamma_gas)) else: Mach_exit = newton(solv_exit_Mach_number, 3.0, args=(ratio_nozzle, Mach_throat, gamma_gas)) velocity_exit = calc_exit_velocity(Mach_throat, Mach_exit, temperature_throat, molar_weight_gas, gamma_gas) pressure_exit = calc_exit_presure(Mach_throat, Mach_exit, pressure_throat, gamma_gas) thrust = massFlowRate_throat * velocity_exit +max((pressure_exit - pressure_environment),0.0)*crossArea_exit #================================ propellant deflagration ==================================# burnrate = calc_burning_rate(pressure_current, burnrate_AN[0], burnrate_AN[1], burnrate_AN[2], burnrate_AN[3], burnrate_AN[4]) propellant_geometry = diameter_inner_propellant, diameter_outer_propellant, height_propellant volume_burnt_current, volume_rest_propellant, propellant_geometry_after = \ evol_geometric_type2(propellant_geometry, burnrate, timeStep) diameter_inner_propellant, diameter_outer_propellant, height_propellant = propellant_geometry_after #=================================================================================================# heatLoss = area_vessel_surface* (temperature_current - temperature_environment)*heat_conduct_coefficient *timeStep mass_supply = mass_supply+ volume_burnt_current* density_propellant* gas_yield_weight energy_supply = energy_supply + volume_burnt_current* density_propellant* gas_yield_weight *combustion_enthalpy energy_supply = energy_supply - heatLoss #================================= update state ===========================================# burnt_volume = burnt_volume +volume_burnt_current volume_vessel_current = volume_vessel_current + volume_burnt_current* gas_yield_volumetric mass_gas = mass_gas + mass_supply temperature_current = cal_temperature(internalEnergy, energy_supply, mass_gas, mass_air) pressure_current = cal_pressure(mass_gas, mass_air, volume_vessel_current, temperature_current) if i % ouputIteration == 0 : print_no_newline(">>"+str("%.2f" % (currentTime))+"s") writeOutput.write(str(currentTime)+" "+str(pressure_current/1e5)+" "+str(temperature_current)+" "+\ str(thrust/1e3)+" "+str(Mach_throat)+" "+str(Mach_exit)+" "+\ str(burnt_volume*1e6)+" "+str(volume_rest_propellant*1e6)+"\n") ``` [__BACK TO CONTENT__](https://hackmd.io/@SamuelChang/H1LvI_eYn)