---
# System prepended metadata

title: 4.2 Characteristics in the internal gas dynamics
tags: [Mathematical modeling and optimization]

---

---
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\}
$$

<!--
<img src="https://encrypted-tbn0.gstatic.com/images?q=tbn:ANd9GcSd4xVNf8sr-wqar6UisRzDYl3GODo-yTkpGQ&s" width="65%"> 
--> 

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">

<br/>
<br/>

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.

<br/>
<br/>
<br/>

#### 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)