J.C. Kantor ([email protected])

# Energy Dispersing Air Bag Landing Systems¶

In [26]:
#Initializations
from IPython.core.display import HTML

Out[26]:
In [27]:
%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt

import math
pi = math.pi

from pint import UnitRegistry
ur = UnitRegistry()


## Air Bag Landing Systems¶

Air bag landing systems have been adapted to a wide range of uses including space vehicle recovery, planetary exploration, humanitarian supply missions, and field delivery of military supplies. The purpose of the airbags is to dissipate the landing energy and thereby safely and securely deliver cargo to the ground. A typical mission profile is shown in the following diagram (Figure 1 from Wen, et. al, 2010).

The air bags are stored in a folded position underneath a payload platform. Upon exiting the aircraft, the parachute deploys and the air bags arefilled by ram air entering through a one-way vent. After initial contact with the ground, the lower vent closes and bag pressure builds due to the payload momentum. A membrane covering the upper vent bursts at a pre-determined pressure after which the airbag continues to collapse until the payload finally and safely touches down. This image shows deployment of Russian towed artillary piece.

A typical airbag design for delivering military cargo is shown in these figures (courtesy of Duramold, Inc.). In this configuration, the airbags are constructed in pairs, and two pairs (a total of four airbags) are installed beneath each 2 feet by nine feet section of a segmented payload platform. A typical platorm (Type V platform manufactured by Airlift Technologies, Inc.) is shown in the accompanying figure.

Typical operating parameters:

• Parachute descent velocity: 28 ft/sec
• Maximum g-force: 8g
• No bouncing

## Preliminary Analysis¶

This is a preliminary analysis to obtain an approximate model for the design analysis of an air bag system. For this purpose, the assumptions include ideal operation of all components, an airbag system of constant cross-sectional area, non-elastic materials, and the ideal gas law.

### Air Bag Dimensions¶

The dimensions of a typical airbag are taken from drawings provided by Duramold, Inc.

In [28]:
# specifications
height = 36.0 * ur.inches
diameter = 36.0 * ur.inches

# calculated
area = 0.25*pi*diameter**2
volume = height*area

print "Airbag area = {0:8.3f}".format(area.to(ur.m**2))
print "Airbag volume = {0:8.3f}".format(volume.to(ur.L))

Airbag area =    0.657 meter ** 2
Airbag volume =  600.480 liter


We'll assume the initial velocity is the standard descent velocity of U.S. Army parachutes (28 ft/sec), and the acceptable landing velocity is 5 ft/sec.

In [29]:
# average loading per airbag
v_initial = 28 * ur.feet/ur.sec
v_final = 5 * ur.feet/ur.sec

print "Initial Descent Velocity = {0:8.4f}".format(v_initial.to(ur.m/ur.sec))
print "Final Descent Velocity = {0:8.4f}".format(v_final.to(ur.m/ur.sec))

Payload Mass = 680.3886 kilogram
Initial Descent Velocity =   8.5344 meter / second
Final Descent Velocity =   1.5240 meter / second


### Energy Analysis¶

We can get a preliminary estimate of the forces and accelerations by assuming uniform negative acceleration of the payload. This spreads the landing event over a maximum period, and gives a lower bound on the peak g-force to be experienced by the payload. We'll neglect the force of the parachute.

In [30]:
KE_initial = 0.5*m_payload*v_initial**2

print "Initial KE = {0:8.3f}".format(KE_initial.to(ur.kJ))
print "  Final KE = {0:8.3f}".format(KE_final.to(ur.kJ))

Initial KE =   24.778 kilojoule
Final KE =    0.790 kilojoule

In [31]:
# Average force
force = (KE_initial-KE_final)/height
print "Average Force = {0:8.3f}".format(force.to(ur.kN))
print "Average Force = {0:8.3f}".format(force.to(ur.lbf))

# Average pressure required
pressure = force/area
print "\nAverage Overpressure = {0:8.3f}".format(pressure.to(ur.kPa))
print "Average Overpressure = {0:8.3f}".format(pressure.to(ur.psi))

# Average decleration
print "Average Acceleration = {0:8.3f}".format(acceleration.to(ur.m/ur.s**2))

# Duration
t = (v_initial-v_final)/acceleration
print "Time = {0:8.3f}".format(t.to(ur.s))

Average Force =   26.234 kilonewton
Average Force = 5897.610 force_pound

Average Overpressure =   39.948 kilopascal
Average Overpressure =    5.794 psi
Average Acceleration =   38.557 meter / second ** 2
Time =    0.182 second


## Descent Phase¶

Starting at rest, to achieve a velocity $v$ in free-fall with gravitational acceleration $g$ with no other significant drags, the period of time required is

$$t_f = \frac{v}{g}$$

and the distance required is

$$y_f = \frac{v^2}{2 g}$$
In [32]:
gravity = 9.81 * ur.m/ur.sec**2

tf = v_initial/gravity
print "Free Fall Time = {0:8.4f}".format(tf.to(ur.sec))

yf = (v_initial**2)/2.0/gravity
print "Free Fall Distance = {0:8.4f}".format(yf.to(ur.m))
print "Free Fall Distance = {0:8.4f}".format(yf.to(ur.ft))

Free Fall Time =   0.8700 second
Free Fall Distance =   3.7123 meter
Free Fall Distance =  12.1796 foot


The parachute descent is slow enough that the air flow is well approximated as an incompressible flow. Then from Bernoulli's law, the stagnation pressure at the bottom vent of the air bag is

$$P_{total} = P_{ambient} + \frac{1}{2}\rho v^2$$

For these low pressures the ideal gas law provides an accurate equation of state. Then

$$\rho = MW_{air}\frac{P_{ambient}}{R T_{ambient}}$$

Assuming the air bag fills at the stagnation pressure and equilibrates with the ambient temperature, thn mass of air in the air bag at the moment it touches down is given by

$$m_{air} = MW_{air}\frac{P_{total} V}{R T_{ambient}}$$
In [33]:
P_ambient = 1.0 * ur.atm
T_ambient = ur.Quantity(15.0,ur.degC).to(ur.degK)

MW = 28.966 * ur.g/ur.mol
R = 8.314 * ur.J/ur.mol/ur.degK

rho = MW*P_ambient/(R*T_ambient)
rho.ito(ur.kg/ur.m**3)

print "Air density = {0:8.3f}".format(rho)

Air density =    1.225 kilogram / meter ** 3

In [34]:
# for low velocities, assume essentially incompressible flow
P_stagnation = 0.5*rho*v_initial**2
P_stagnation.ito(ur.Pa)

print "Stagnation Pressure = {0:8.5f}".format(P_stagnation)

Stagnation Pressure = 44.61619 pascal

In [35]:
P_initial = P_ambient + P_stagnation

print "Total Pressure = {0:8.5f}".format(P_initial)

Total Pressure =  1.00044 atmosphere

In [36]:
m_air = MW*P_initial*volume/(R*T_ambient)
print "Air mass = {0:8.3f}".format(m_air.to(ur.g))

Air mass =  735.980 gram


## Airbag Dynamics without Vent¶

Let $h$ denote the height of the cargo platform above the ground, and let $v$ denote velocity. At the moment the bottom of the air bag touches down, the dynamics become

\begin{align*} \frac{dh}{dt} & = v \\ \frac{dv}{dt} & = \frac{(P_{bag}-P_{ambient})A_{bag}}{m_{payload}} - g \end{align*}

where

\begin{align*} h(0) & = H_{bag} \\ v(0) & = -v_{parachute} \end{align*}

$v_{parachute}$ is descent velocity from the parachute deployment phase, and $H_{bag}$ is the height of the bag when fully extended. We assume the parachute drag is immediately released upon touchdown.

The bottom inlet vent immediately closes on touchdown. The air mass in the bag is constant until the upper vent membrane bursts. For this first analysis, we'll assume the vent membrane remains intact.

$$\frac{dE}{dt}=\dot{Q}+\dot{W}$$

where $\dot{Q}$ is the heat transferred to the system from the surroundings and $\dot{W}$ is the work done on the system. The entire compression phase is only a few hundred milliseconds, so to a good approximation the adiabatic assumption $\dot{Q} = 0$ holds.

The specific energy of an ideal gas is a function of temperature alone, so

\begin{aligned} \frac{dE}{dt} & =\frac{d(n_{air}\hat{U})}{dt}\\ & =n_{air}\frac{d\hat{U}}{dt}\\ & =n_{air}\frac{\partial\hat{U}}{\partial T}\frac{dT}{dt}\\ & =n_{air}C_{v}\frac{dT}{dt} \end{aligned}

On the work side of the equation, the rate of work done on the system ($\dot{W}$)

\begin{aligned} \dot{W} & = - P \frac{dV}{dt} \end{aligned}

Giving an energy balance

$$n_{air}C_{v}\frac{dT}{dt} = -P\frac{dV}{dt}$$

Since it is pressure that appears in other parts of the model, we'd like to have the left side of this equation in terms of pressure. Substituting the ideal gas and taking a total derivative of temperature gives

$$\frac{C_{v}}{R}\left(V \frac{dP}{dt} + P\frac{dV}{dt}\right) = -P\frac{dV}{dt}$$

or

$$\frac{1}{P} \frac{dP}{dt} = -\left(1+\frac{R}{C_v}\right) \frac{1}{V} \frac{dV}{dt}$$

For an ideal gas, $C_p = R + C_v$. Rearranging we get

$$\frac{dP}{dt} = - \frac{P}{V} \left(\frac{C_p}{C_v}\right)\frac{dV}{dt}$$

Putting this all together for a constant area ($V = hA$) leaves a system of three differential equations with well-defined initial conditions for the initial compression phase of the air bag.

In [37]:
R = 8.314 * ur.J/(ur.mol*ur.degK)
MW = 28.966 * ur.grams/ur.mol

P_ambient = 1.0 * ur.atm

\begin{align*} \frac{dh}{dt} & = v \\ \frac{dv}{dt} & = \frac{(P - P_{ambient})A}{m_{payload}} - g \\ \frac{dP}{dt} & = -\frac{P}{h} \left(\frac{C_p}{C_v}\right)v \end{align*}

where

\begin{align*} h(0) & = H_{bag} \\ v(0) & = -v_{parachute} \\ P(0) & = P_{total} \end{align*}
In [38]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Convert all parameters to base units
gravity = 9.81 * ur.m/ur.sec**2
k = 1.4
area.ito(ur.m**2)
P_ambient.ito(ur.Pa)

# Verify parameter values, some calculated in earlier cells

print "gravity = {0:8.4f}".format(gravity)
print "k = {0:8.4f}".format(k)
print "area = {0:8.4f}".format(area)
print "P_ambient = {0:8.4f}".format(P_ambient)

# Set and verify initial conditions

h_ic = (height.to(ur.m)).magnitude
v_ic = (-v_initial.to(ur.m/ur.sec)).magnitude
P_ic = (P_initial.to(ur.Pa)).magnitude

x0 = [h_ic, v_ic, P_ic]

print "\nInitial Conditions for Initial Compression Phase"
print "h_ic = {0:8.4f}".format(h_ic)
print "v_ic = {0:8.4f}".format(v_ic)
print "P_ic = {0:8.1f}".format(P_ic)

def phase1(x,t):
h,v,P = x
dhdt = v
if h > h_ic:
dvdt = -gravity.magnitude
dPdt = 0
else:
dvdt = -gravity.magnitude + \
dPdt = -P*k*dhdt/h
return np.array([dhdt,dvdt,dPdt])

t = np.linspace(0,3.00,1000)
soln = odeint(phase1,x0,t);

h = soln[:,0] * ur.m
v = soln[:,1] * ur.m/ur.sec
P = soln[:,2] * ur.Pa

Pr = P/P_initial.to(ur.Pa)
Tr = Pr**(2.0/7.0)

T = Tr*T_ambient

gravity =   9.8100 meter / second ** 2
k =   1.4000
area =   0.6567 meter ** 2
P_ambient = 101325.0000 pascal

Initial Conditions for Initial Compression Phase
h_ic =   0.9144
v_ic =  -8.5344
P_ic = 101369.6

In [39]:
plt.figure(figsize=(10,4))
plt.plot(t,h,t,v)
plt.title('Height')
plt.ylabel('meters')
plt.legend(['Height','Velocity'])

plt.figure(figsize=(10,4))
plt.plot(t,v)
plt.title('Velocity')
plt.ylabel('meters/sec')

plt.figure(figsize=(10,4))
plt.plot(t,gForce)
plt.title('Braking G Force')
plt.ylabel("G's")

plt.figure(figsize=(10,4))
plt.plot(t,P.to(ur.kPa))
plt.title('Air Bag Pressure (absolute)')
plt.ylabel('kiloPascals')

plt.figure(figsize=(10,4))
plt.plot(t,T.to(ur.degC))
plt.title('Air Bag Temperature')
plt.ylabel('deg C')

Out[39]:
<matplotlib.text.Text at 0x10ac38f90>

## Venting Air Bag¶

### Mass Balance¶

Venting the air bag provides a mechanism to release accumulated energy. The mass balance is given by

$$\frac{d\left(\rho V\right)}{dt} = - C_DA_o \rho_{o} q_{o}$$

where $\rho_o$ and $q_o$ refer to the density and velocity of the exit gas at the orifice conditions, and where $A_o$ is the area of the orifice. For unchoked flow, the density and pressure at the orifice are equal to ambient conditions, so $\rho_o = \rho_a$.

The air bag is compressed by the rapid deacceleration of the payload. From an adiabatic energy balance on an ideal gas, the pressure and density of the remaining air inside the bag increase follow the relationship

$$\frac{\rho}{\rho_a} = \left(\frac{P}{P_a}\right)^{\frac{1}{\gamma}}$$

where the subscript $a$ refers to ambient conditions. Applying the chain rule

$$\rho\frac{dV}{dt} + V\frac{d\rho}{dt} = - C_DA_o \rho_{a} q_{o}$$$$\frac{V\rho_a}{\gamma P_a}\left(\frac{P}{P_a}\right)^{\frac{1}{\gamma}-1}\frac{dP}{dt} = - C_DA_o \rho_a q_o - \rho_a\ \left(\frac{P}{P_a}\right)^{\frac{1}{\gamma}} \frac{dV}{dt}$$

Cancelling terms

$$\frac{V}{\gamma P_a}\left(\frac{P}{P_a}\right)^{\frac{1}{\gamma}-1}\frac{dP}{dt} = - C_DA_o q_o - \left(\frac{P}{P_a}\right)^{\frac{1}{\gamma}} \frac{dV}{dt}$$$$\frac{V}{\gamma P_a}\left(\frac{P}{P_a}\right)^{\frac{1}{\gamma}}\frac{P_a}{P}\frac{dP}{dt} = - C_DA_o q_o - \left(\frac{P}{P_a}\right)^{\frac{1}{\gamma}} \frac{dV}{dt}$$

Leaving

$$\frac{1}{P}\frac{dP}{dt} = - \frac{\gamma C_DA_o q_o}{V \left(\frac{P}{P_a}\right)^{\frac{1}{\gamma}}} - \frac{\gamma}{V}\frac{dV}{dt}$$

### Bernoulli's Equation¶

To a close approximation, energy balances on the streamlines from the interior of the airbag satisfy Bernoulli's law

$$\frac{\gamma}{\gamma-1} \frac{P_o}{\rho_o} + \frac{q_o^2}{2} + gz_o = \mbox{constant}$$

For unchoked flow, $P_o = P_a$ and $\rho_o = \rho_a$ at the orifice. Assuming a vent with a small cross-sectional area compared to the bag dimensions, and neglecting minor elevation effects,

$$\frac{q_o^2}{2} + \frac{\gamma}{\gamma-1} \frac{P_a}{\rho_a} = \frac{\gamma}{\gamma-1} \frac{P}{\rho}$$

so solving for $q_o$ and using the adiabatic relationship one more time

$$q_o =\sqrt{\frac{2\gamma}{\gamma-1}\frac{P_a}{\rho_a}\left[\left(\frac{P}{P_a}\right)^{1-\frac{1}{\gamma}} - 1\right]}$$

Using the ideal gas relationship

$$\frac{P_a}{\rho_a} = \frac{R T_a}{MW}$$

leaves an expression for orifice velocity

$$q_o =\sqrt{\frac{2\gamma}{\gamma-1}\frac{R T_a}{MW}\left[\left(\frac{P}{P_a}\right)^{1-\frac{1}{\gamma}} - 1\right]}$$
In [40]:
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# orifice parameters

Cd = 0.72
Ao = 20*8 * ur.cm**2

# air parameters
k = 1.4
gamma = k
MW = 0.028966 * ur.kg/ur.mol

# Convert all parameters to base units
gravity = 9.81 * ur.m/ur.sec**2
k = 1.4
area.ito(ur.m**2)
P_ambient.ito(ur.Pa)
T_ambient.ito(ur.degK)
Ao.ito(ur.m**2)

# Verify parameter values, some of which were calculated above

print "gravity = {0:8.4f}".format(gravity)
print "k = {0:8.4f}".format(k)
print "area = {0:8.4f}".format(area)
print "P_ambient = {0:8.4f}".format(P_ambient)
print "T_ambient = {0:8.4f}".format(T_ambient)

# Set and verify initial conditions

h_ic = (height.to(ur.m)).magnitude
v_ic = (-v_initial.to(ur.m/ur.sec)).magnitude
P_ic = (P_initial.to(ur.Pa)).magnitude

x0 = [h_ic, v_ic, P_ic]

print "\nInitial Conditions for Initial Compression Phase"
print "h_ic = {0:8.4f}".format(h_ic)
print "v_ic = {0:8.4f}".format(v_ic)
print "P_ic = {0:8.1f}".format(P_ic)

def phase1(x,t):
h = x[0]
v = x[1]
P = x[2]
pr = P/P_ambient.magnitude
cr = ((gamma+1)/2.0)**(gamma/(gamma-1.0))
#if pr > cr:
#q = np.sqrt(gamma*R.magnitude*T_ambient.magnitude/MW.magnitude)
#q = q*np.sqrt((2/(gamma+1))**((gamma+1)/(gamma-1)))
# else:
#q = np.sqrt(2*(gamma/(gamma-1))*(R.magnitude*T_ambient.magnitude/MW.magnitude))
#q = q*np.sqrt(pr**(1.0 - (1.0/gamma))-1.0)
if pr > cr:
pr = cr
q = np.sqrt(2*(gamma/(gamma-1))*(R.magnitude*T_ambient.magnitude/MW.magnitude))
q = q*np.sqrt(pr**(1.0 - (1.0/gamma))-1.0)
dhdt = v
if h > h_ic:
dvdt = -gravity.magnitude
dPdt = 0
# elif h <= 0.01:
#     dhdt = 0
#    dvdt = 0
#    dPdt = 0
else:
dvdt = -gravity.magnitude + \
dPdt = -P*k*dhdt/h - P*gamma*Cd*Ao.magnitude*q/(h*area.magnitude*pr**(1.0/gamma))
return np.array([dhdt,dvdt,dPdt])

t = np.linspace(0,0.3,1000)
soln = odeint(phase1,x0,t);

h = soln[:,0] * ur.m
v = soln[:,1] * ur.m/ur.sec
P = soln[:,2] * ur.Pa

Pr = P/P_initial.to(ur.Pa)
Tr = Pr**(2.0/7.0)

T = Tr*T_ambient

plt.figure(figsize=(10,4))
plt.plot(t,h,t,v)
plt.title('Height')
plt.ylabel('meters')
plt.legend(['Height','Velocity'])

plt.figure(figsize=(10,4))
plt.plot(t,v)
plt.title('Velocity')
plt.ylabel('meters/sec')

plt.figure(figsize=(10,4))
plt.plot(t,gForce)
plt.title('Braking G Force')
plt.ylabel("G's")

plt.figure(figsize=(10,4))
plt.plot(t,P.to(ur.kPa))
plt.title('Air Bag Pressure (absolute)')
plt.ylabel('kiloPascals')

plt.figure(figsize=(10,4))
plt.plot(t,T.to(ur.degC))
plt.title('Air Bag Temperature')
plt.ylabel('deg C')

gravity =   9.8100 meter / second ** 2
k =   1.4000
area =   0.6567 meter ** 2
P_ambient = 101325.0000 pascal
T_ambient = 288.1500 kelvin

Initial Conditions for Initial Compression Phase
h_ic =   0.9144
v_ic =  -8.5344
P_ic = 101369.6

Out[40]:
<matplotlib.text.Text at 0x108ffbc50>

## Crush Zone Dynamics¶

We'll assume that a 3 inch crush zone of cardboard honeycomb or other structure is the final cushion for the falling mass. We'll further assume the maximum tolerable de-acceleration rate is 8g. From this data we can determine the maximum permissable velocity of the payload at completion of the air bag phase.

In [41]:
max_accel = 8*g
h_crush = 3.0 * ur.inches

t_crush = np.sqrt(2.0*h_crush/max_accel)
print "Crush Zone Time = {0:8.4f}".format(t_crush.to_base_units())

v_crush = max_accel*t_crush
print "Crush Zone Velocity = {0:8.4f}".format(v_crush.to_base_units())

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-41-a695337b1fec> in <module>()
----> 1 max_accel = 8*g
2 h_crush = 3.0 * ur.inches
3
4 t_crush = np.sqrt(2.0*h_crush/max_accel)
5 print "Crush Zone Time = {0:8.4f}".format(t_crush.to_base_units())

NameError: name 'g' is not defined