Diving into Mathematics
or some Mathematics of Scuba Diving

Michel DE LARA


Table des matières

1. Introduction

A diver is submitted to kinematic and dynamic laws for its movements in water: gravity, water drag force, buyoancy of Archimedes, perfect-gas law are the essential physical elements which contribute to the diver's kinematics.

The diver's safety depends on the dynamics of dissolved nitrogen in blood: Dalton's law coupled with linear dynamics provide the so called Haldane model of compartments for dissolved nitrogen in blood. More recent and accurate models exist among which the famous Reduced Gradient Bubble Model. The diver safety depends on the realism of such models: indeed, they are coded in tables or in wrist computers which inform the diver and prescribe diving profiles (stops and ascents) during the dive.

We present here the basic elements to build models of diver which combine both type of physics. The diver state consists of position (depth), impulsion, mass, air-jacket mass and dissolved nitrogen in compartments. The diver has control variables to modify his or her state: body orientation (area of the solid projection orthogonal to the speed vector), palming, jacket filling and unfilling. The state follows a controlled ordinary differential equation. Controls belong to bounded sets and linear combinations of the state must satisfay given inequalities for safety and physical reasons.

Such controlled ordinary differential equation with contraints models allow us to formulate various control problems:

At this stage, we can present the buiding of the model, suggest control problems and give some hints as to their resolution.

2. The diver in movement:
kinematics and dynamics

We consider an equipped diver (cylinder, suit, jacket) identified with a rigid body with

2.1 Depth evolution:
impulsion, forces and Newton's law of motion

2.1.1 Impulsion

The impulsion is

\begin{displaymath} \textbf{q}= \textbf{m}\dot{\textbf{z}} . \end{displaymath} (2.2)

2.1.2 Newton's law of motion

The impulsion satisfies, by Newton's law (A.2), the differential equation

\begin{displaymath} \displaystyle \frac{d \textbf{q}}{dt} = \textbf{m}g - \rho_{w} \textbf{v}g+ \mbox{drag} (\dot{\textbf{z}}) + \textbf{f}\, . \end{displaymath} (2.3)

2.1.3 Gravity force: $\textbf{m}g$

2.1.4 Archimedes force: $- \rho _{w} \textbf {v}g$

2.1.5 Palming force: $\textbf{f}$

2.1.6 Drag force: $\mbox{drag}(\dot{\textbf{z}})$

The drag force is resulting from the equipped diver's shape, type of diving suit, etc.: we take

\begin{displaymath} \mbox{drag} (\dot{\textbf{z}}) = - \frac{1}{2} \rho_{w}\, C\!x\, S \, \mbox{sign}(\dot{\textbf{z}}) \dot{\textbf{z}}^2 \end{displaymath} (2.4)

where the coefficients $C\!x$ and $S$ are described in paragraph A.3.3.

2.2 Lungs and diver's body volume evolution:
physiological laws

2.2.1 Diver's body volume: $\textbf {v}_{body}$

We define the diver's body volume $\textbf {v}_{body}$ as the one with empty lungs and including equipment (cylinder, suit, empty jacket). We assume the diver's body to be incompressible and thus $\textbf {v}_{body}$ to have a fixed value.

2.2.2 Lungs volume: $\textbf {v}_{air-lung}$ (stationary model)

The simplest lungs volume evolution model consists in assuming that $\textbf {v}_{air-lung}$ is fixed, taken equal to zero (meaning thus that this volume is included in the diver's body volume).

2.2.3 Lungs volume: $\textbf {v}_{air-lung}$ (periodic model)

Amplitude: volume $c$ of air breathed per unit of time

A physiological property of the human body is the following: whatever the pressure, one breathes (and thus expires) the same volume of air per unit of time (approximately $20~l min^{-1}$). We denote by $c$ the volume of air breathed per unit of time:

\begin{displaymath} c\approx 0.33~ l s^{-1}= 3.3~ 10^{-4}\, m^3 s^{-1} \, . \end{displaymath} (2.5)

Periodicity: expiration and inspiration

On the other hand, the volume $\textbf {v}_{air-lung}$ in lungs is periodical due to expiration and inspiration. We denote by $\tau$ the period ( $\tau \approx 12 s$).

A more elaborate lungs volume evolution model

For both above reasons, we assume that

\begin{displaymath} \textbf{v}_{air-lung}(t) = c\upsilon(t) \quad\mbox{where}\qu... ...) &=&\upsilon(t) \, , \quad \forall t \, . \end{array} \right. \end{displaymath} (2.6)

A typical shape for $\upsilon (t)$ is given on picture 2.1. Normalized volume variation $\upsilon (t)$ in lungsf:lungs

2.3 Pressure: hydrostatic distribution

The pressure $\textbf{p}$ is supposed to follow the hydrostatic distribution, given here by

\begin{displaymath} \textbf{p}= p_{a} + \rho_{w} g\textbf{z} . \end{displaymath} (2.7)

$p_{a}$ is the atmospheric pressure and $\rho_{w}$ is the water mass per unit of volume.

2.4 Diving suit and jacket volume evolution:
perfect-gas law

2.4.1 Diving suit volume: $\textbf {v}_{gas-suit}$

The diving suit contains a number $\textbf{n}_{gas-suit}$ of moles of gas. By the perfect-gas law (A.10), the volume of the gas captured in the diving suit satisfies

\begin{displaymath} \textbf{p}\textbf{v}_{gas-suit} = \textbf{n}_{gas-suit} R T\, . \end{displaymath} (2.8)

Note that the number of moles $\textbf{n}_{gas-suit}$ may be estimated by
\begin{displaymath} \textbf{n}_{gas-suit} = \frac{p_{a} \textbf{v}_{gas-suit}^0}{R T^0} \end{displaymath} (2.9)

where $\textbf{v}_{gas-suit}^0$ is the volume of the gas captured in the diving suit at atmospheric pressure $p_{a}$ and at temperature $T^0$ ( $\textbf{v}_{gas-suit}^0 \approx 8 l = 8.10^{-3} m^3$).

2.4.2 Jacket volume: $\textbf {v}_{air-jac}$

The diving suit contains a number $\textbf{n}_{gas-suit}$ of moles of gas. By the perfect-gas law (A.10), the volume of the gas captured in the diving suit satisfies

\begin{displaymath} \textbf{p}\textbf{v}_{gas-suit} = \textbf{n}_{gas-suit} R T\, . \end{displaymath} (2.10)

The volume $\textbf {v}_{air-jac}$ of the air-jacket is related to the number $\textbf{n}_{air-jac}$ of moles of air in the jacket by

\begin{displaymath} \textbf{p}\textbf{v}_{air-jac} = \textbf{n}_{air-jac} R T\, . \end{displaymath} (2.11)

2.4.3 Diving suite and jacket moles of gas: $\textbf {n}_{gas}$

The diving suit and jacket volume evolutions may be treated jointly by introducing the number $\textbf {n}_{gas}$ of moles of gas in the diving suite (fixed) and in the jacket (variable):

\begin{displaymath} \textbf{p}(\textbf{v}_{gas-suit} + \textbf{v}_{air-jac}) = \textbf{n}_{gas} R T\, . \end{displaymath} (2.12)

2.5 Mass evolution

With obvious notations, the mass $\textbf{m}$ of the equipped diver is given by

\begin{displaymath} \textbf{m}= \textbf{m}_{body} + \textbf{m}_{air-lung} + \textbf{m}_{air-cyl} . \end{displaymath} (2.13)

The diver's body mass $\textbf{m}_{body}$ (including everything, body, jacket, cylinder... but not the air in cylinder, neither in lungs, nor in jacket) is supposed to be stationary.

2.5.1 Lungs air mass: $\textbf{m}_{air-lung}$

The lungs air mass $\textbf{m}_{air-lung}$ is related to the air volume $\textbf {v}_{air-lung}$ in lungs:

\begin{displaymath} \textbf{m}_{air-lung} = \rho_{air} \textbf{v}_{air-lung} \, . \end{displaymath} (2.14)

The simplest lungs air mass evolution model consists in assuming that $\textbf{m}_{air-lung}$ is fixed, taken equal to zero (meaning thus that this mass is included in the diver's body mass).

2.5.2 Cylinder's air mass variation rate: $\dot{\textbf{m}}_{air-cyl}$

The cylinder's air mass $\textbf{m}_{air-cyl}$ decreases with time, both for lung and jacket filling. The rate $\dot{\textbf{m}}_{air-cyl}$ of cylinder's air mass variations is related to the variations $\dot{\textbf{v}}_{air-lung}$ of volume in lungs:

\begin{displaymath} \dot{\textbf{m}}_{air-cyl} = - \dot{\textbf{m}}_{air-lung} \... ...{\textbf{m}}_{air-lung} \geq 0} - \dot{\textbf{m}}_{cyl->jac} \end{displaymath} (2.15)

where $\dot{\textbf{m}}_{cyl->jac} \geq 0$ is the rate of air mass transfer from cylinder to jacket.

Decrease for lungs filling: a simple model

We assume here that

\begin{displaymath} \dot{\textbf{m}}_{air-cyl} = - c\rho_{air} = - \frac{M_{air}}{R T} (p_{a} + \rho_{w} g\textbf{z}) c\, . \end{displaymath} (2.16)

With obvious notations, we also have

\begin{displaymath} \dot{\textbf{m}}_{air-jac} = \dot{\textbf{m}}_{cyl->jac} - \dot{\textbf{m}}_{jac->wat} \end{displaymath} (2.17)

where
\begin{displaymath} \left\{ \begin{array}{rcl} {\textbf{m}}_{air-jac} & \geq & ... ...dot{\textbf{m}}_{jac->wat} & \geq & 0 . \end{array} \right. \end{displaymath} (2.18)

2.6 Temperature variation with depth

3. The diver breathing: nitrogen saturation and desaturation

3.1 A model of compartments for dissolved nitrogen

Decompression theory - neo-Haldane models Let $\alpha_{N_2}$ denote the percentage of nitrogen in air. The partial $N_2$ pressure $\textbf{p}_{N_2}$, to which the diver is submitted at ambient pressure $\textbf{p}$, derives from Dalton's law:

\begin{displaymath} \textbf{p}_{N_2} = \alpha_{N_2} \textbf{p} . \end{displaymath} (3.1)

For nitrogen saturation and desaturation in the blood, the human body is schematically described by $n$ compartments, where each compartment $i$ is characterized by

For instance, official French tables are characterized by the following values:
period (minutes) 5 7 10 15 20 30 40 50 60 80 100 120
$S\!c$ (dimensionless) 2.72 2.54 2.38 2.20 2.04 1.82 1.68 1.61 1.58 1.56 1.55 1.54

By Henry's law, equilibrium between dissolved $N_2$ in blood and gazeous $N_2$ in lungs is reached when $\textbf{p}^{i}_{N_2}=\textbf{p}_{N_2}$ in compartment $i$. The dynamics of $N_2$ dissolution in blood is given by the differential system:

\begin{displaymath} \left\{ \begin{array}{rcl} \displaystyle \frac{d \textbf{p}^... ...textbf{p}^{n}_{N_2} - \textbf{p}_{N_2})\\ \end{array} \right. \end{displaymath} (3.2)


\begin{displaymath} \lambda_i \stackrel{\mbox{\tiny def}}{=}\frac{\log 2}{\theta_i} \end{displaymath} (3.3)

3.1.1 Safety values

Experiments have shown that nitrogen bubbles do not appear in blood as long as the over saturation critical coefficient is not reached, that is when:

\begin{displaymath} \textbf{p}^{i}_{N_2} \leq S\!c_i \textbf{p} . \end{displaymath} (3.4)

More elaborate models have been developed for the formation of nitrogen bubbles.

What is more, it is asked that, in the ascending phase, the decrease of the ambient pressure $\textbf{p}$ be bounded. In practice, this requires that the diver ascending speed be bounded by $s_{max}$ ( $s_{max}=12 m/min = 0.2 m s^{-1}$):

\begin{displaymath} \dot{\textbf{z}} < 0 \Rightarrow \vert \dot{\textbf{z}} \vert \leq s_{max} \end{displaymath} (3.5)

or
\begin{displaymath} \textbf{m}s_{max} + \textbf{q}\geq 0 . \end{displaymath} (3.6)

3.1.2 A dynamic model for nitrogen saturation, desaturation and bubbles formation

Summing up the above relationships, we get the following controlled linear state model with constraints as follows.


\begin{displaymath} \frac{d}{dt} \left( \begin{array}{c} \textbf{p}^{1}_{N_2} \\... ... \displaystyle -\frac{\log 2}{\theta_n} \end{array} \right) \end{displaymath} (3.7)


\begin{displaymath} \quad\mbox{Constraints: }\quad \left\{ \begin{array}{l} \fo... ... \textbf{m}s_{max} + \textbf{q}\geq 0 . \end{array} \right. \end{displaymath} (3.8)

3.2 The Reduced Gradient Bubble Model

Bubble Decompression Strategies

Bubble Decompression Strategies

4. Diver's equations

4.1 A simple model for local control


\begin{displaymath} \frac{d}{dt} \left( \begin{array}{c} \textbf{z}\\ [5mm] \do... ...^2 + \textbf{f}\\ [5mm] \dot{\textbf{n}} \end{array} \right) \end{displaymath} (4.1)


\begin{displaymath} \quad\mbox{Control: }\quad (\textbf{f}, \dot{\textbf{n}}) \end{displaymath} (4.2)


\begin{displaymath} \quad\mbox{Constraints: }\quad \left\{ \begin{array}{rcl} \t... ...x} & \geq & \mid \dot{\textbf{n}} \mid \, \end{array} \right. \end{displaymath} (4.3)

4.2 A first simple model for global control


\begin{displaymath} \frac{d}{dt} \left( \begin{array}{c} \textbf{z}\\ [5mm] \te... ...pha_{N_2} (p_{a} + \rho_{w} g\textbf{z}) ) \end{array} \right) \end{displaymath} (4.4)


\begin{displaymath} \quad\mbox{Control: }\quad (\textbf{f}, \dot{\textbf{n}}) \end{displaymath} (4.5)


\begin{displaymath} \quad\mbox{Constraints: }\quad \left\{ \begin{array}{rcl} \t... ... s_{max} + \dot{\textbf{z}} & \geq & 0 \, \end{array} \right. \end{displaymath} (4.6)

4.3 A general model


\begin{displaymath} \begin{array}{cl} & \displaystyle \quad\mbox{State equation:... ...(p_{a} + \rho_{w} g\textbf{z})) \end{array} \right) \end{array}\end{displaymath} (4.7)


\begin{displaymath} \quad\mbox{Control: }\quad (\textbf{f}, \dot{\textbf{m}}_{cyl->jac}, \dot{\textbf{m}}_{jac->wat} ) \end{displaymath} (4.8)


\begin{displaymath} \quad\mbox{State constraints: }\quad \left\{ \begin{array}{r... ...tbf{z}) - \textbf{p}^{n}_{N_2} & \geq & 0 \end{array} \right. \end{displaymath} (4.9)


\begin{displaymath} \quad\mbox{Control constraints: }\quad \left\{ \begin{array}... ...eq & \textbf{f}& \geq & -\textbf{f}_{max} \end{array} \right. \end{displaymath} (4.10)

5. Control problems

We try and formulate different control problems naturally arising in scuba diving.

5.1 A global model

With

\begin{displaymath} \left\{ \begin{array}{rcl} x & \stackrel{\mbox{\tiny def}}{=... ...}, \dot{\textbf{m}}_{jac->wat} ) \in \RR^2 \end{array} \right. \end{displaymath} (5.1)

we have a controlled state model with constraints as follows:
\begin{displaymath} \dot x = f(t,x,u) , \quad Rx -r \geq 0 , \quad u \geq 0 . \end{displaymath} (5.2)

5.2 Stabilization at given depth


\begin{displaymath} \forall t \in [t_1,t_2] , \quad \textbf{z}(t) = \textbf{z}_{max} \end{displaymath} (5.3)

5.3 Ascent at fixed speed


\begin{displaymath} \textbf{z}(T)=0 \quad\mbox{and}\quad \forall t \in [t_2,T] \, , \quad \dot{\textbf{z}}(t) = s_{max} \leq 0 \end{displaymath} (5.4)

5.4 Maximum sojourn time at given depth


\begin{displaymath} \left\{ \begin{array}{ll} \displaystyle \sup \int_{t_1}^{t_... ... t \in [t_1,t_2] \\ [4mm] S x(T)=s \, , & \end{array} \right. \end{displaymath} (5.5)


\begin{displaymath} \sup \int_0^{+\infty} \1_{ \{ \textbf{z}(t) \geq h \} } dt \end{displaymath} (5.6)

5.5 Minimum upward time


\begin{displaymath} \inf \int_0^{+\infty} \1_{ \{ \dot{\textbf{z}}(t) \leq 0 \} } dt \end{displaymath} (5.7)

5.6 With observation of one compartment

$\textbf{p}^{1}_{N_2}$ measured

A. Appendix

A.1 Universal constants

$p_{a}$ atmospheric pressure $p_{a}$ = $ 1.01 325 10^5 Pa$
$\rho_{w}$ water mass per unit of volume $\rho_{w}$ = $ 10^3 kg \; m^{-3}$
$g$ gravitation acceleration $g$ = $ 9.81 m s^{-2}$
$M_{air}$ air molecular weight $M_{air}$ = $ 2.91 10^{-2} kg \; mol^{-1}$
$k$ Boltzmann constant $k$ = $1.38 10^{-23} J K^{-1}$
$R$ perfect-gas constant $R$ = $8.314 J K^{-1} mol^{-1}$
${\cal N}$ Avogadro constant ${\cal N}$ = $6.022 10^{23} mol^{-1}$
$\alpha_{N_2}$ percentage of nitrogen in air $\alpha_{N_2}$ = $0.81$

A.2 Other constants

$c$ volume of air breathed per unit of time $c$ = $ 3.3 10^{-4} m^3 s^{-1}$
$\textbf{v}_{gas-suit}^0$ volume of gas in the diving suit $\textbf{v}_{gas-suit}^0$ = $ 8.10^{-3} m^3$
$C\!x$ drag coefficient of the diver $C\!x$ $\approx$ 1
$S$ area of the diver projection orthogonal to speed $S$ $\approx$ 1

A.3 Basic physical laws for diver's mechanics

A.3.1 Newton's law

A punctual mass $\textbf{m}$ located at depth $\textbf{z}$, moving vertically in a fluid, has an impulsion

\begin{displaymath} \textbf{q}= \textbf{m}\dot{\textbf{z}} \end{displaymath} (A.1)

and satisfies
\begin{displaymath} \frac{d \textbf{q}}{dt} = \mbox{gravitation force} - \mbox{... ...of stress} \pm \mbox{other forces (friction, palming...)} . \end{displaymath} (A.2)

This remains valid for a rigid body, where $\textbf{z}$ is the depth of the center of gravity.

A.3.2 Pressure and force of stress for an ideal fluid

For scuba diving problems, we make the (reasonable) assumption that water is an ideal fluid, homogeneous in the horizontal dimension. Thus, there exists a function of depth $\textbf{z}$ ( $\textbf{z}\geq 0$), the pressure $\textbf{p}=\textbf{p}(\textbf{z})$ such that if $S$ is any surface of the fluid with a unit normal vector $\stackrel{\to}{\textbf{n}}$, the force of stress exerced by the fluid per unit of area on $S$ is given by:

\begin{displaymath} \quad\mbox{force of stress across \emph{S} per unit area}\quad = \textbf{p}\stackrel{\to}{\textbf{n}} . \end{displaymath} (A.3)

By Stoke's law (divergence theorem), we get that
\begin{displaymath} \quad\mbox{force of stress per unit volume}\quad = - \stackrel{\rightarrow}{\nabla}\textbf{p} . \end{displaymath} (A.4)

Application: hydrostatic distribution

For water at rest ( $\dot{\textbf{z}} =0$), Newton's law (A.2) reduces to

\begin{displaymath} \stackrel{\rightarrow}{\nabla}\textbf{p}= \rho_{w} \stackrel{\to}{g} . \end{displaymath} (A.5)

The resulting pressure $\textbf{p}$ is the hydrostatic distribution, given here by:
\begin{displaymath} \textbf{p}= p_{a} + \rho_{w} g\textbf{z} . \end{displaymath} (A.6)

Application: The two laws of buyoancy of Archimedes

By equations (A.4) and (A.5), pressure exerts, on a volume $W$, a force

\begin{displaymath} \int_W (- \stackrel{\rightarrow}{\nabla}\textbf{p}) dw = - \rho_{w} \int_W \stackrel{\to}{g}dw . \end{displaymath} (A.7)

This can be formulated as Archimedes' laws.
  1. A body immersed in a fluid experiences a vertical buyoant force equal to the weight of the fluid it displaces, having origin at the center of buoyancy (center of mass computed as if with uniform density).
  2. A floating body displaces its own weight in the fluid in which it floats.


A.3.3 Drag force

The drag force is opposed to the speed $\dot{\textbf{z}}$.

At low speed


\begin{displaymath} \mbox{drag force} = - \kappa \eta \dot{\textbf{z}} \end{displaymath} (A.8)

where $\kappa$ depends on the geometry and has the dimension of a length ( $\kappa= 6 \pi r$ for a sphere of radius $r$), and $\eta$ is the viscosity coefficient ( $\eta_{water} \approx 10^{-3}$ at $20°C$).

At higher speeds


\begin{displaymath} \mbox{drag force} = - \frac{1}{2} \rho_{w} C\!x S \dot{\textbf{z}}^2 \end{displaymath} (A.9)

where the drag coefficient $C\!x$ is a dimensionless number depending on the solid shape ($C\!x\approx 1$), $S$ is the area of the solid projection orthogonal to the speed vector ($S \approx 1$).

A.3.4 Perfect-gas law

Let a gas occupy the volume $\textbf{v}$, be submitted to the pressure $\textbf{p}$, be at temperature $T$, consist of $\textbf{n}$ moles (and $N=\textbf{n}{\cal N}$ molecules). Then, perfect-gas law states that

\begin{displaymath} \left\{ \begin{array}{rcl} \textbf{p}\textbf{v}&=& \display... ...uad\mbox{is the Boltzmann constant.}\quad \end{array} \right. \end{displaymath} (A.10)

If we introduce the gas molecular weight $M_{gas}$, this can be rephrased in term of gas density $\rho_{gas}$ (mass per unit of volume)
\begin{displaymath} \rho_{gas} = \displaystyle \frac{\textbf{n}M_{gas}}{\textbf{v}} = \frac{M_{gas}}{R T} \textbf{p}\, . \end{displaymath} (A.11)

A.3.5 Bernouilli's law

In a non turbulent fluid with volumic mass $\rho$, at pressure $\textbf{p}$, with speed $u$, the quantity $\textbf{p}+ \frac{1}{2} \rho u^2$ is uniform in the fluid.

A.3.6 Cylinder pressure evolution

Let $\textbf{p}_{cylinder}$ be the pressure inside the cylinder. According to the perfect-gas law, the pressure inside the cylinder decreases in proportion to the decrease of air particles:

\begin{displaymath} \begin{array}{rcl} \dot{\textbf{p}}_{cylinder} &=& \displays... ...cylinder}M_{air}} \dot{\textbf{m}}_{cyl->jac} \, . \end{array}\end{displaymath} (A.12)

Rate of air mass transfer from cylinder to jacket

When we action the so called direct system, air particles flee from the cylinder with a speed $u_{cyl->jac}$ given by Bernouilli's law:

\begin{displaymath} \textbf{p}_{cylinder} = \textbf{p}+ \frac{1}{2} \rho_{air-cyl} u_{cyl->jac}^2 . \end{displaymath} (A.13)

The rate of air mass transfer from cylinder to jacket $\dot{\textbf{m}}_{cyl->jac}$ through a surface $\sigma$ is
\begin{displaymath} \dot{\textbf{m}}_{cyl->jac} = \sigma \rho_{air-cyl} u_{cyl->jac} . \end{displaymath} (A.14)

We find
\begin{displaymath} \dot{\textbf{m}}_{cylinder->jacket} = \sigma \sqrt{2 \frac{... ...xtbf{p}_{cylinder} (\textbf{p}_{cylinder} - \textbf{p})} . \end{displaymath} (A.15)

A.4 Basic (bio)physical laws for diver's nitrogen saturation and desaturation

A.4.1 Partial pressures and Dalton's law

The partial pressure $p_G$ of a gas $G$ in air derives from Dalton's law:

\begin{displaymath} \textbf{p}_G = (\% G \textrm{in air}) \textbf{p}_{air} . \end{displaymath} (A.16)

A.4.2 Gas tension, Henry's law

At equilibrium, the tension $T_G$ of a gas $G$ in a liquid (blood) is equal to the partial pressure $p_G$ of the gas $G$ at the surface of the liquid (Henry's law):

\begin{displaymath} T_G = \textbf{p}_G . \end{displaymath} (A.17)

A.4.3 Dynamics of saturation and desaturation

The dynamics of saturation and desaturation that leads to Henry's law's equilibrium is governed by a specific time, the period $\theta$, and by the law:

\begin{displaymath} \frac{d T_G}{dt} = - \frac{\log 2}{\theta} (T_G - \textbf{p}_G) . \end{displaymath} (A.18)

A.4.4 Bubbles formation

Bubbles do not appear as long as the over saturation critical coefficient is not reached, that is when:

\begin{displaymath} T_G \leq S\!c\textbf{p}= S\!c\frac{\textbf{p}_G}{\% G \textrm{ in air}} . \end{displaymath} (A.19)

A.5 Scilab code

A.5.1 diving4.sce 

//exec('diving4.sce')

getf('diving4.sci')
// charge le fichier de fonctions

//////////////////////////////////////////////////////////////////////
 //                   DONNEES
//////////////////////////////////////////////////////////////////////

Init_diving();
// data initialisation

// MASS
mass_body=70;
mass_suit=4;
mass_cylinder=12;
mass=mass_body+mass_suit+mass_cylinder;

// VOLUME 
vbody=7*10^(-2);

// MOLES
vgassuit=8*10^(-3);
ngassuit=patm*vgassuit/(R*Tatm);


//////////////////////////////////////////////////////////////////////
 //     CONDITIONS INITIALES D'UNE EDO    
//////////////////////////////////////////////////////////////////////

// état initial des compartiments
w=zeros(16,1);
w(2)=mass_body; w(5:16)=-patm;
//Quand on part a 0m de profondeur, on suppose que 
// tous les caissons sont a la pression atmospherique
w(2)=mass_body; w(5:16)=-patm*Sc';
//Quand on part du fond, on suppose que tous les caissons sont saturés

// initial volume in the air-jacket at depth 0
[air_mass0,volume0]=balance(0,mass,1)
// vérification
diving(0,[0,0,mass,0],[0,0,0])
diving(0,[0,0,mass,air_mass0],[0,0,0])

// volume in the air-jacket at depth "depth"
//depth=20;
[air_mass20,volume20]=balance(20,mass,0)
// vérification
diving(6,[20,0,mass,air_mass20],[0,0,0])
// diving(6... for breathing(6)=0, i.e. empty lungs

// initial acceleration
x0=[0,0,mass,air_mass0]';
// à la surface
x0=[20,0,mass,air_mass20]';
// à 20 mètres


//////////////////////////////////////////////////////////////////////
 //     RESOLUTION D'UNE EDO    
//////////////////////////////////////////////////////////////////////

// EDO
t=0:300;
//deff('y=f(t,x)','y=diving(t,x,0*air_mass0/120*feedback(x))');
deff('y=f(t,x)','y=diving_Haldane(t,x,[0,0,0])');
y=ode([x0; patm*alpha*ones(1,12)'],0,t,f);

deff('y=f(t,x)','y=diving(t,x,[0,0,0])');
y=ode(x0,0,t,f);



//////////////////////////////////////////////////////////////////////
 //     VISUALISATION DE TRAJECTOIRES
//////////////////////////////////////////////////////////////////////

//y(1,$+1)=-1000;
//[u,v]=mini(sign(y(1,2:$)));
// v est l'indice pour lequel la profondeur est negative !
// i.e. le plongeur sort de l'eau...
//tt=1:(v-1);
//xbasc();plot(t(tt),-y(1,tt));


// Afin de toujours tracer des vecteurs de même taille, 
// meme en cas de divergence de l'integration, 
// ne prendre que les valeurs de t pour lesquelles la solution y
// a effectivement ete calculee
s=size(y);
//xbasc();plot(t(1):t(s(2)),-y(1,:));
//attention: le pas de temps dans plot doit etre le meme que dans t

xdel(1);
xset("window",1);
//xtitle("profondeur du plongeur (mètres) en fonction du temps (secondes)");
xtitle("diver''s depth (meters) with time (seconds)");
plot2d(t(1):t(s(2)),-y(1,:));
// fenêtre 1 : profondeur du plongeur

xdel(2);
xset("window",2);
plot2d(t(1):t(s(2)),-60*y(2,:)./y(3,:));
// fenêtre 2 : vitesse du plongeur en m/min

xdel(3);
xset("window",3);
temps=t(1):t(s(2));
CU=diag([ones(1:4),10^(-5)*ones(5:16)]);
// changement d'unité pour les pressions
//plot2d(ones(1:16)'*temps,W*y-w*ones(y(1,:));
plot2d((ones(1:16)'*temps)',(CU*(W*y-w*ones(y(1,:))))');

//sign(W*y-w*ones(y(1,:)));
// le signe de la matrice W*y-w*ones(y(1,:)) nous informe sur 
// le premier moment ou les contraintes ne sont plus vérifiées (Haldane...)



//////////////////////////////////////////////////////////////////////
 //     COMMANDE
//////////////////////////////////////////////////////////////////////


[air_mole20]=balance_0(20,mass)
x0=[20,0,air_mole20]';
// à 20 mètres
diving_0(0,x0,[0,0,0])

deff('y=ff(t,x,u)','y=diving_0(t,x,[(-u)*Heavyside(-u),u*Heavyside(u),0])');
ff(0,x0,0)

[f,g,complin]=tangent('ff',x0,0);
contr(f,g);

k=ppol(f,g,[-1 -1 -1]);

function [xdot]=diving_0_cont(t,x)
xdot=ff(t,x,-k*(x-x0))
endfunction

t=0:300;
y=ode(x0,0,t,diving_0_cont);
s=size(y);
xdel(1);
xset("window",1);
//xtitle("profondeur du plongeur (mètres) en fonction du temps (secondes)");
xtitle("diver''s depth (meters) with time (seconds)");
xbasc();plot2d(t(1):t(s(2)),-y(1,:),rect=[t(1),-x0(1)-2,t(s(2)),2]);
// fenêtre 1 : profondeur du plongeur



[air_mass,volume]=balance(20,mass,0);
x0=[20,0,mass,air_mass]'

function [xdot]=diving_cont(t,x)
  u=-[k(1) k(2) 0 k(3)]*(x-x0);
xdot=diving(t,x,[(-u)*Heavyside(-u),u*Heavyside(u),0])
endfunction

t=0:300;
y=ode(x0,0,t,diving_cont);
s=size(y);
xdel(2);
xset("window",2);
//xtitle("profondeur du plongeur (mètres) en fonction du temps (secondes)");
xtitle("diver''s depth (meters) with time (seconds)");
xbasc();plot2d(t(1):t(s(2)),-y(1,:),rect=[t(1),-x0(1)-2,t(s(2)),2]);
// fenêtre 1 : profondeur du plongeur

A.5.2 diving4.sci 

function []=Init_diving()
// initialise les données

kk=1.38*10^{-23};
R=8.314;
T=273+18;
Tatm=273+25;
patm=1.01325*10^5;
gg=9.81;
rhow=10^3;
nu=6*%pi*10^{-3};
cc=3.3*10^{-4};
Mair=2.91*10^{-2};
S=0.5;
Cx=1;
theta=60*[5 7 10 15 20 30 40 50 60 80 100 120]';
alpha=0.81;
Sc=[2.72 2.54 2.38 2.20 2.04 1.82 1.68 1.61 1.58 1.56 1.55 1.54];
// coefficients de sécurité
smax=0.2;
W=zeros(16,16);
W(1,1)=1; W(2,3)=1; W(3,4)=1;W(4,2)=1;W(4,3)=smax;
W(5:16,5:16)=diag(-ones(5:16));
W(5:16,1)=rhow*gg*Sc';

[kk,R,T,Tatm,patm,gg,rhow,nu,cc,Mair,S,Cx,theta,alpha,Sc,smax,W]=resume...
(kk,R,T,Tatm,patm,gg,rhow,nu,cc,Mair,S,Cx,theta,alpha,Sc,smax,W);
// resume



//////////////////////////////////////////////////////////////////////
 //     FONCTIONS INTERMEDIAIRES
//////////////////////////////////////////////////////////////////////

function [vol]=breathing(t)
dd=pmodulo(t,12)
vol=abs(dd-6)*2*cc
// volume total inspiré sur 12 s = 12*cc
// soit bien cc par seconde

function [vol]=expiring(t)
dd=pmodulo(t,12)
vol=2*cc*Heavyside(6-dd)
// volume expiré par unité de temps

function [H] = Heavyside(x)
// Fonction de Heavyside
    H = max(sign(x+%eps),0);
// so that Heavyside(0)=1


//////////////////////////////////////////////////////////////////////
 //     SYSTÈME DYNAMIQUE
//////////////////////////////////////////////////////////////////////


// MODELE SANS DESATURATION SANS RESPIRATION

function [xdot]=diving_0(t,x,u)
// x(1) : depth z
// z>0 si le plongeur est dans l'eau
// x(2) : speed 
// x(3) : number of moles of air in the jacket + suit
// u(1) : rate of air mass transfer from cylinder to jacket (mdotcj)
// u(2) : rate of air mass transfer from jacket to water (mdotjw)
// u(3) : palming force (palm)
xdot(1)=x(2);
xdot(2)= Heavyside(x(1)).*( mass*gg - rhow*vbody*gg - ...
( (rhow*R*T) ./ (patm+rhow*gg.*x(1)) ) .* ( (ngassuit+x(3)) * gg )...
- 0.5*rhow*S*Cx.*sign(x(2)).*((x(2)).^2) + u(3) );
// hors de l'eau on n'a que dq/dt=m*g d'où le Heavyside au dessus
xdot(3)=Heavyside(x(1)).*( u(1)-u(2) ) ;
// hors de l'eau il n'y a plus d'échanges gazeux ( cyl->jac,cyl->wat...)
endfunction

function [air_mole]=balance_0(depth,mass)
// air moles in the air-jacket
// to ensure stability of the diver with mass "mass" at depth "depth"
air_mole= (mass-rhow*vbody ).* ...
(patm+rhow*gg.*depth)./(rhow*R*T)-ngassuit ; 
endfunction



// MODELE SANS DESATURATION AVEC RESPIRATION

function [xdot]=diving(t,x,u)
// x(1) : depth z
// z>0 si le plongeur est dans l'eau
// x(2) : impulsion q
// x(3) : mass m
// x(4) : mass of air in the jacket
// u(1) : rate of air mass transfer from cylinder to jacket (mdotcj)
// u(2) : rate of air mass transfer from jacket to water (mdotjw)
// u(3) : palming force (palm)
//u(3)=0.01*rand(x(1));
  xdot(1)=x(2)./x(3);
xdot(2)=  x(3)*gg + Heavyside(x(1)).*( ...
- rhow*(vbody+breathing(t))*gg - ...
(rhow*R*T)./(patm+rhow*gg.*x(1)).*(ngassuit+x(4)/Mair)*gg...
- 0.5*rhow*S*Cx.*sign(x(2)).*((x(2)./x(3))^2) + u(3) );
// + 0.1*rand(x(3));
// hors de l'eau on n'a que dq/dt=m*g d'où le Heavyside au dessus
xdot(3)=Heavyside(x(1)).*(...
-expiring(t)*Mair/(R*T)*(patm+rhow*gg.*x(1)) - u(1) ) ; 
xdot(4)=Heavyside(x(1)).*( u(1)-u(2) ) ;
// hors de l'eau il n'y a plus d'échanges gazeux ( cyl->jac,cyl->wat...)

function [air_mass,volume]=balance(depth,mass,breath)
// air mass and corresponding volume in the air-jacket
// to ensure stability of the diver with mass "mass" at depth "depth"
// when his lungs are full (breath=1) or empty (breath=0)
air_mass=Mair*((mass-rhow*(vbody+breathing(6*(1-breath)))).*...
(patm+rhow*gg.*depth)./(rhow*R*T)-ngassuit); 
volume=(air_mass/Mair*R*T)./(patm+rhow*gg.*depth);


// MODELE AVEC DESATURATION

function [xdot]=diving_Haldane(t,x,u)
// x(1) : depth z
// z>0 si le plongeur est dans l'eau
// x(2) : impulsion q
// x(3) : mass m
// x(4) : mass of air in the jacket
// x(5:16) : N2 dans compartiments
// u(1) : rate of air mass transfer from cylinder to jacket (mdotcj)
// u(2) : rate of air mass transfer from jacket to water (mdotjw)
// u(3) : palming force (palm)
xdot(1)=x(2)./x(3);
xdot(2)=  x(3)*gg + Heavyside(x(1)).*( ...
- rhow*(vbody+breathing(t))*gg - ...
(rhow*R*T)./(patm+rhow*gg.*x(1)).*(ngassuit+x(4)/Mair)*gg...
- 0.5*rhow*S*Cx.*sign(x(2)).*((x(2)./x(3))^2) + u(3) );
// hors de l'eau on n'a que dq/dt=m*g d'où le Heavyside au dessus
xdot(3)=Heavyside(x(1)).*(...
-expiring(t)*Mair/(R*T)*(patm+rhow*gg.*x(1)) - u(1) ) ; 
xdot(4)=Heavyside(x(1)).*( u(1)-u(2) ) ;
// hors de l'eau il n'y a plus d'échanges gazeux ( cyl->jac,cyl->wat...)
xdot(5:16)=-log(2)*diag(1 ./theta)*x(5:16)-...
           alpha*log(2)*(1 ./theta)*(patm+rhow*gg.*x(1));




//////////////////////////////////////////////////////////////////////
 //     COMMANDE
//////////////////////////////////////////////////////////////////////



function [u]=feedback(x)
u(1)=5*(x(1)-20)*Heavyside(x(1)-20);
// A z > 20m, le plongeur injecte une quantité d'air (heavyside(0)=1)
u(2)=(20-x(1))*Heavyside(20-x(1));
// when , eject air from the jacket. 
// Plus le plongeur remonte plus il vide le gilet
u(3)=0;
// pas de force de palmage