Michel DE LARA

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

• stabilization at given depth ;
• ascent at fixed speed ;
• maximum sojourn time at given depth ;
• minimum upward time.

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

• mass ;
• volume given, with obvious notations, by
 (2.1)

• gravity center located at depth ( is the surface).

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

The impulsion is

 (2.2)

## 2.1.2 Newton's law of motion

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

 (2.3)

## 2.1.6 Drag force:

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

 (2.4)

where the coefficients and 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:

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

## 2.2.2 Lungs volume: (stationary model)

The simplest lungs volume evolution model consists in assuming that is fixed, taken equal to zero (meaning thus that this volume is included in the diver's body volume).

## 2.2.3 Lungs volume: (periodic model)

### Amplitude: volume 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 ). We denote by the volume of air breathed per unit of time:

 (2.5)

### Periodicity: expiration and inspiration

On the other hand, the volume in lungs is periodical due to expiration and inspiration. We denote by the period ( ).

### A more elaborate lungs volume evolution model

For both above reasons, we assume that

 (2.6)

A typical shape for is given on picture 2.1. Normalized volume variation in lungsf:lungs

# 2.3 Pressure: hydrostatic distribution

The pressure is supposed to follow the hydrostatic distribution, given here by

 (2.7)

is the atmospheric pressure and 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:

The diving suit contains a number of moles of gas. By the perfect-gas law (A.10), the volume of the gas captured in the diving suit satisfies

 (2.8)

Note that the number of moles may be estimated by
 (2.9)

where is the volume of the gas captured in the diving suit at atmospheric pressure and at temperature ( ).

## 2.4.2 Jacket volume:

The diving suit contains a number of moles of gas. By the perfect-gas law (A.10), the volume of the gas captured in the diving suit satisfies

 (2.10)

The volume of the air-jacket is related to the number of moles of air in the jacket by

 (2.11)

## 2.4.3 Diving suite and jacket moles of gas:

The diving suit and jacket volume evolutions may be treated jointly by introducing the number of moles of gas in the diving suite (fixed) and in the jacket (variable):

 (2.12)

# 2.5 Mass evolution

With obvious notations, the mass of the equipped diver is given by

 (2.13)

The diver's body mass (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:

The lungs air mass is related to the air volume in lungs:

 (2.14)

The simplest lungs air mass evolution model consists in assuming that 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:

The cylinder's air mass decreases with time, both for lung and jacket filling. The rate of cylinder's air mass variations is related to the variations of volume in lungs:

 (2.15)

where is the rate of air mass transfer from cylinder to jacket.

### Decrease for lungs filling: a simple model

We assume here that

 (2.16)

With obvious notations, we also have

 (2.17)

where
 (2.18)

# 3.1 A model of compartments for dissolved nitrogen

Decompression theory - neo-Haldane models Let denote the percentage of nitrogen in air. The partial pressure , to which the diver is submitted at ambient pressure , derives from Dalton's law:

 (3.1)

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

• a variable , the quantity of dissolved in blood (unit is as pressure),
• a parameter , the period,
• a parameter , the over saturation critical coefficient.
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 (dimensionless) 2.72 2.54 2.38 2.2 2.04 1.82 1.68 1.61 1.58 1.56 1.55 1.54

By Henry's law, equilibrium between dissolved in blood and gazeous in lungs is reached when in compartment . The dynamics of dissolution in blood is given by the differential system:

 (3.2)

 (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:

 (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 be bounded. In practice, this requires that the diver ascending speed be bounded by ( ):

 (3.5)

or
 (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.

 (3.7)

 (3.8)

 (4.1)

 (4.2)

 (4.3)

 (4.4)

 (4.5)

 (4.6)

 (4.7)

 (4.8)

 (4.9)

 (4.10)

# 5. Control problems

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

# 5.1 A global model

With

 (5.1)

we have a controlled state model with constraints as follows:
 (5.2)

 (5.3)

 (5.4)

 (5.5)

 (5.6)

 (5.7)

measured

# A.1 Universal constants

 atmospheric pressure = water mass per unit of volume = gravitation acceleration = air molecular weight = Boltzmann constant = perfect-gas constant = Avogadro constant = percentage of nitrogen in air =

# A.2 Other constants

 volume of air breathed per unit of time = volume of gas in the diving suit = drag coefficient of the diver 1 area of the diver projection orthogonal to speed 1

# A.3 Basic physical laws for diver's mechanics

## A.3.1 Newton's law

A punctual mass located at depth , moving vertically in a fluid, has an impulsion

 (A.1)

and satisfies
 (A.2)

This remains valid for a rigid body, where 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 ( ), the pressure such that if is any surface of the fluid with a unit normal vector , the force of stress exerced by the fluid per unit of area on is given by:

 (A.3)

By Stoke's law (divergence theorem), we get that
 (A.4)

### Application: hydrostatic distribution

For water at rest ( ), Newton's law (A.2) reduces to

 (A.5)

The resulting pressure is the hydrostatic distribution, given here by:
 (A.6)

### Application: The two laws of buyoancy of Archimedes

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

 (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 .

### At low speed

 (A.8)

where depends on the geometry and has the dimension of a length ( for a sphere of radius ), and is the viscosity coefficient ( at ).

### At higher speeds

 (A.9)

where the drag coefficient is a dimensionless number depending on the solid shape (), is the area of the solid projection orthogonal to the speed vector ().

## A.3.4 Perfect-gas law

Let a gas occupy the volume , be submitted to the pressure , be at temperature , consist of moles (and molecules). Then, perfect-gas law states that

 (A.10)

If we introduce the gas molecular weight , this can be rephrased in term of gas density (mass per unit of volume)
 (A.11)

## A.3.5 Bernouilli's law

In a non turbulent fluid with volumic mass , at pressure , with speed , the quantity is uniform in the fluid.

## A.3.6 Cylinder pressure evolution

Let 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:

 (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 given by Bernouilli's law:

 (A.13)

The rate of air mass transfer from cylinder to jacket through a surface is
 (A.14)

We find
 (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 of a gas in air derives from Dalton's law:

 (A.16)

## A.4.2 Gas tension, Henry's law

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

 (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 , and by the law:

 (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:

 (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