1 Description

Le présent programme est un simulateur de comportement mécanique de corps solides indéformables dans un environnement fluide. Son objet principal est la résolution des équations de mouvements du ou des corps considérés, dans le domaine temporel. Ces équations de mouvement sont construites sur la base d’efforts extérieurs calculés par des modèles spécifiques inclus dans le logiciel, ou pouvant y être ajoutés. Chaque modèle est associé à un jeu de données spécifique.

Sur le principe, le simulateur permet de traiter n’importe quel problème de mouvement de corps solide. Les premiers modèles d’efforts implémentés et fournis concernent essentiellement le comportement de navires sur houle.

Concernant le mode de fonctionnement du logiciel, il s’agit d’un outil en ligne de commande, c’est-à-dire qu’on le lance depuis un terminal (ou une invite de commande MS DOS).

Le simulateur a été construit pour permettre une utilisation modulaire en ligne de commande. Pour ce faire, les différents domaines à configurer sont séparés : la description du problème physique est faite dans un ou plusieurs fichiers d’entrée et la description de la simulation (pas de temps, solveur…) est faite en ligne de commande. De cette manière, on peut facilement lancer des simulations du même problème avec des solveurs différents et sur des durées différentes, sans toucher au(x) fichier(s) de configuration du problème.

Le fichier de description du problème peut être séparé en autant de fichiers que l’on veut, ce qui permet par exemple de versionner la configuration de l’environnement physique séparément de celle du ou des solides.

En sortie, le simulateur génère un fichier de résultat au format CSV contenant tous les états du ou des solides. Il peut également générer un fichier YAML contenant des hauteurs de houle sur un maillage. La figure suivante résume les entrées et sorties du simulateur :

Parallèlement à la présente documentation utilisateur, il existe également une documentation d’implémentation décrivant l’architecture logicielle et détaillant les API pour faciliter la maintenance du code (correctifs et ajouts de fonctionnalités).

2 Ligne de commande

2.1 Liste des arguments

xdyn <yaml file> [-hd] [-y ARG] [-s ARG] [dt ARG] [--tstart ARG]
[--tend ARG] [-o ARG] [-w ARG]
Options:
  -h [ --help ]              Show this help message
  -y [ --yml ] arg           Name(s) of the YAML file(s)
  -s [ --solver ] arg (=rk4) Name of the solver: euler, rk4, rkck for
Euler,
                             Runge-Kutta 4 & Runge-Kutta-Cash-Karp
                             respectively.
  --dt arg                   Initial time step (or value of the fixed
time step
                             for fixed step solvers)
  --tstart arg (=0)          Date corresponding to the beginning of
the
                             simulation (in seconds)
  --tend arg                 Last time step
  -o [ --output ] arg        Name of the output file where all
computed data
                             will be exported.
                             Possible values/extensions are csv, tsv,
json,
                             hdf5, h5, ws
  -w [ --waves ] arg         Name of the output file where the wave
heights
                             will be stored ('output' section of the
YAML
                             file). In case output is made to a HDF5
file or
                             web sockets, this option appends the wave
height
                             to the main output
  -d [ --debug ]             Used by the application's support team to
help
                             error diagnosis. Allows us to pinpoint
the exact
                             location in code where the error occurred
(do not
                             catch exceptions), eg. for use in a
debugger.

2.2 Exemples

Les exemples suivants peuvent être lancés à partir du répertoire demos:

2.2.1 Simulation avec un solveur Runge-Kutta d’ordre 4 en commençant à t=0

./xdyn tutorial_01_falling_ball.yml --dt 0.1 --tend 1

ou si l’on utilise une invite de commande MS-DOS :

xdyn tutorial_01_falling_ball.yml --dt 0.1 --tend 1

En ajoutant l’option -o csv, les sorties se font sur la sortie standard (le terminal). Ceci permet de chaîner les traitements (pipe UNIX), par exemple :

./xdyn tutorial_01_falling_ball.yml --dt 0.1 --tend 1 -o csv| python plot.py test 0 3

La commande précédente lance la simulation et génère un tracé (format SVG) à l’aide du script python de post-traitement livré avec le simulateur. Des scripts de post-traitement sont livrés avec l’installation du simulateur, pour Matlab, Python.

2.2.2 Simulation avec un solveur Euler en commençant à t=1

./xdyn tutorial_01_falling_ball.yml -s euler --dt 0.1 --tstart 1 --tend 1.2

3 Documentations des données d’entrées du simulateur

Les données d’entrées du simulateur se basent sur un format de fichiers YAML qui fonctionne par clef-valeur.

Ce format présente l’avantage d’être facilement lisible et éditable. Des parsers YAML existent pour de nombreux langages de programmation.

3.1 Vue d’ensemble

3.1.1 Liste des sections

Le fichier YAML comprend quatre sections de haut niveau :

• la section rotations convention définit la convention d’angles utilisée,
• environmental constants donne les valeurs de la gravité et la densité de l’eau,
• les modèles environnementaux figurent dans environment models,
• bodies décrit les corps simulés.

3.1.2 Remarques sur les unités

On retrouve fréquemment dans le fichier YAML des lignes du type :

clef: {value: 65456, unit: km}

Les unités ne sont pas vérifiées par le système : le parser se contente de convertir toutes les entrées en unité du système international avant simulation, sans préjuger de l’homogénéité. Ce décodage est fait en convertissant l’unité en un facteur multiplicatif (en utilisant la liste des unités de l’utilitaire UNIX units) et en multipliant la valeur originale par ce facteur multiplicatif. Par exemple:

masse: {value: 10, unit: lb}

sera vu par xdyn comme:

masse = 4.5359237;

L’exemple précédent :

clef: {value: 65456, unit: km}

aurait tout aussi bien pu être écrit :

clef: {value: 65456, unit: kW}

et on aurait obtenu exactement la même valeur numérique, bien que la grandeur physique ne soit pas la même : on suppose donc que l’utilisateur renseigne des données de façon homogène. En interne, tous les calculs sont faits en unité du système international.

3.2 Sorties

La spécification des sorties se fait au moyen de la section output, à la racine du fichier YAML, dont voici un exemple :

output:
   - format: csv
     filename: test.csv
     data: [t, x(ball), 'Fx(gravity,ball)']

• format : csv pour un fichier texte dont les colonnes sont séparées par une virgule ou hdf5 pour le format des fichiers .mat de MatLab (HDF5)
• filename : nom du fichier de sortie
• data : liste des colonnes à écrire. Le temps est noté t, et les états sont x(body), y(body) z(body), u(body), v(body), w(body), p(body), q(body), r(body), qr(body), qi(body), qj(body), qk(body). body doit être remplacé par le nom du corps (ball dans l’exemple ci-dessus). Les sorties d’effort sont Fx(modèle,corps,repère), Fy(modèle,corps,repère), Fz(modèle,corps,repère), Mx(modèle,corps,repère), My(modèle,corps,repère), Mz(modèle,corps,repère) où modèle est le nom du modèle d’effort (renseigné dans la clef modèle de chaque modèle d’effort), corps est le nom du corps sur lequel agit l’effort et repère est le repère d’expression (qui ne peut être que NED ou le nom du corps). Les sorties de houle sont notées waves et leur contenu est décrit dans la section Modèle de houle/Sorties. La somme des efforts appliqués à un corps est accessible par Fx(sum of forces,corps,repère) (resp. Fy, Fz, Mx, My, Mz).

4 Interface MatLab

xdyn peut être appelé depuis le logiciel MatLab. Cela présente l’avantage de disposer dans la foulée d’un environnement graphique pour afficher les résultats de simulation.

Voici les fonctions de base pour travailler avec xdyn.

• xdyn_run exécute xdyn depuis MatLab.
• xdyn_loadResultsFromHdf5File permet de charger les résultats d’une simulation dans MatLab.
• xdyn_postProcess lance l’ensemble des post-traitements disponibles.
• xdyn_plotStates permet de tracer les états, à partir de résultats de simulations. Cela comprend les positions et les vitesses de chaque corps.
• xdyn_plotPositions permet de tracer les positions et les orientations de chaque corps.
• xdyn_plotVelocities permet de tracer les vitesses de translation et de rotation de chaque corps dans le repère de chaque repère.
• xdyn_animate3d permet de lancer une animation 3d d’une simulation avec les objets simulés et le champ de vagues lorsque celui-ci est exporté.
• xdyn_extractMatLabScript permet d’extraire les scripts MatLab de dépouillement intégrés dans le fichier de sortie hdf5.

Enfin, le fichier xdyn_demos.m permet de lancer les différents tutoriels.

Il est également possible d’utiliser l’environnement MatLab pour analyser des simulations déjà réalisées. Par exemple, le script xdyn_demos.m va lire le fichier tutorial_01_falling_ball.h5 qui aura préalablement été généré à partir de la commande :

./xdyn tutorial_01_falling_ball.yml --dt 0.1 --tend 1 -o tutorial_01_falling_ball.h5

5 Interface Docker et génération automatique de rapports

xdyn peut être utilisé sous la forme d’un conteneur Docker ce qui permet de l’exécuter sur n’importe quelle plateforme.

En outre, c’est actuellement le seul moyen d’utiliser le module de génération automatique de rapport décrit ci-après.

Pour utiliser l’image xdyn, il faut d’abord installer Docker puis récupérer l’image sirehna/xdyn en exécutant :

docker pull sirehna/xdyn

Pour utiliser l’image xdyn avec son environnement de post-traitement, on récupérera l’image sirehna/xweave en exécutant :

docker pull sirehna/xweave

5.1 Lancement d’xdyn via l’image docker

Une fois l’image chargée par la commande précédente, on peut lancer :

docker run -it --rm -v $(pwd):/work -w /work sirehna/xdyn

• le paramètre -v $(pwd):/work permet de faire correspondre le répertoire courant avec le répertoire /work du conteneur,
• le paramètre -w /work précise que le répertoire de travail du conteneur (celui depuis lequel sera exécuté xdyn) est /work,
• sirehna/xdyn correspond au nom de l’image Docker,

Si l’on souhaite lancer une simulation, on ajoute les arguments d’xdyn à la suite :

docker run -it --rm -v $(pwd):/work -w /work sirehna/xdyn tutorial_01_falling_ball.yml --dt 0.1 --tend 1 -o tsv

5.2 Génération automatique de rapports (X-Weave)

Cette fonctionnalité n’est accessible qu’en utilisant le conteneur Docker sirehna/xweave.

5.2.1 Principe de fonctionnement

On utilise Pweave pour inclure des balises Python dans du texte au format Markdown). Par exemple, la présente documentation est générée en utilisant ce système, y compris les tutoriels et leurs courbes de résultats.

Grâce à ce générateur (appelé “x-weave”), on peut obtenir des rapports dans le format que l’on souhaite (DOCX, PDF, HTML…). Ces rapports peuvent par exemple lancer des simulations xdyn et inclure des courbes à partir de ces résultats.

L’intérêt est que l’on peut gérer le code de ses rapports en configuration (puisque c’est du texte), regarder facilement les différences entre deux versions sans utiliser d’outil propriétaire et regénérer à la demande une nouvelle version du rapport lorsque l’on change les entrées.

X-Weave n’est ni plus ni moins que Pweave dans un containeur Docker avec xdyn pré-installé, ce qui permet de lancer xdyn depuis Pweave.

5.2.2 Ligne de commande

La ligne de commande à utiliser est :

docker run --rm -it -v $(pwd):/work -w /work -u $(id -u):$(id -g) --entrypoint /usr/bin/xdyn-weave sirehna/xweave -f markdown concatenated_doc.pmd -o concatenated_doc.md

• La première partie de la ligne de commande (docker run --rm -it -v $(pwd):/work -w /work -u $(id -u):$(id -g) --entrypoint /usr/bin/xdyn-weave sirehna/xweave) contient les arguments passés à Docker pour partager le dossier courant avec le containeur et créer les fichiers de sortie avec les permissions de l’utilisateur courant.
• La partie après sirehna/xweave contient les arguments passés à pweave. Une liste complète de ces arguments peut être obtenue en exécutant docker run --rm -it -v $(pwd):/work -w /work -u $(id -u):$(id -g) --entrypoint /usr/bin/xdyn-weave sirehna/xweave -h ou en consultant la documentation de Pweave.

On obtient alors un fichier Markdown et (le cas échéant) des images. On peut ensuite utiliser Pandoc pour convertir le Markdown au format de son choix en utilisant le Pandoc installé dans le containeur X-Weave :

docker run --rm -it -v $(shell pwd):/work -w /work -u $(shell id -u):$(shell id -g) --entrypoint /usr/bin/pandoc xweave concatenated_doc_pandoc.md -o doc.html

5.2.3 Commandes Python supplémentaires

Par rapport à Pweave “classique”, on a inclus quelques fonctions Python pour simplifier le travail avec xdyn.

5.2.3.1 Chargement du fichier YAML : load_yaml

L’insertion de la section suivante dans un document X-weave ne génèrera pas de sortie. Par contre, dans les sections python ultérieures, on pourra utiliser les données ainsi chargées.

```python
yaml_data = load_yaml('tutorial_06_1D_propulsion.yml')
```

5.2.3.2 Affichage de l’intégrabilité du YAML : print_yaml

Suite au chargement du YAML, on peut vouloir l’afficher dans le document. Pour cela, on utilise une section python de la forme :

```python
yaml_data = load_yaml('tutorial_06_1D_propulsion.yml')

# Pour afficher l'intégralité du YAML avec de la coloration syntaxique
print_yaml(yaml_data)

# Pour n'afficher qu'une sous-section du YAML
print_yaml(yaml_data, 'bodies/0/controlled forces')
```

5.2.3.3 Exécution d’une simulation

La commande suivante ne produira pas de contenu visible dans le document. En revanche, il sera possible de récupérer et afficher les données ainsi générées.

```python
execCommand('xdyn tutorial_01_falling_ball.yml --dt 0.01 --tend 1 -o out.csv')
```

5.2.3.4 Tracé des résultats

```python
# Récupération des données générées
data = csv('out.csv')

# Définition de la donnée à tracer
plot = prepare_plot_data(data, x = 't', y = 'z(ball)', name='Résultat')

# Définition du type de graph
g = cartesian_graph([plot], x='t (s)', y='Élévation (m)')

# Tracé de la planche
create_layout(g, title='Élévation au cours du temps')
```

5.3 Utilisation de xdyn en serveur websocket

5.3.1 Description

Il est possible de lancer xdyn en tant que serveur, afin de l’intégrer à un autre environnement de simulation. L’utilisation d’un serveur plutôt qu’une MEX-fonction (MATLAB executable - fonction) ou un FMU (Functional Mock-up Unit) permet d’exécuter xdyn sur une autre machine et de l’interroger par plusieurs clients.

Il s’agit d’une utilisation en “model exchange” (au sens de la spécification “Functional Mockup Interface”), par opposition à “xdyn for Co-simulation”. La différence se situe dans l’utilisation du solveur : dans le cas de la co-simulation, on utilise le solveur interne de xdyn (le serveur renvoie les états intégrés au pas suivant $X(t+dt)$). Dans le cas “model exchange”, le serveur renvoie la dérivée des états $\frac{dX}{dt}$.

5.3.2 Justification technique

L’utilisation des websockets permet des temps de réponse plus courts (puisque c’est un protocole assez léger, comparé au HTTP par exemple). Dans l’implémentation actuelle, les messages envoyés sont en JSON, pour offrir un bon compromis entre la verbosité (moins que du XML mais plus qu’un format binaire) et une utilisation plus aisée qu’un format type Protobuf ou Thrift, quitte à sacrifier un peu de performance (taille des messages, temps d’encodage/décodage).

5.3.3 Lancement du serveur “Model Exchange”

Le serveur est lancé grâce à l’exécutable xdyn-for-me (xdyn for Model Exchange) avec un ou plusieurs fichiers YAML en paramètre : une fois lancé, ce serveur ne simulera donc qu’une seule configuration de l’environnement et du navire. Concrètement, on lance le serveur comme suit :

./xdyn-for-me --port 9002 tutorial_01_falling_ball.yml

où --port sert à définir le port sur lequel écoute le serveur websocket.

La liste complète des options avec leur description est obtenue en lançant l’exécutable avec le flag -h.

Ensuite, on peut se connecter à l’adresse du serveur pour l’interroger.

5.3.4 Lancement du serveur “Co-Simulation”

Le serveur est lancé grâce à l’exécutable xdyn-for-cs (xdyn for Co-Simulation) avec un ou plusieurs fichiers YAML en paramètre, un pas de temps et un solveur : une fois lancé, ce serveur ne simulera donc qu’une seule configuration de l’environnement et du navire, en utilisant un solveur et un pas de temps. Concrètement, on lance le serveur comme suit :

./xdyn-for-cs --port 9002 tutorial_01_falling_ball.yml --dt 0.1

où --port sert à définir le port sur lequel écoute le serveur websocket.

La liste complète des options avec leur description est obtenue en lançant l’exécutable avec le flag -h.

Ensuite, on peut se connecter à l’adresse du serveur pour l’interroger.

5.3.5 Utilisation avec Chrome

Le navigateur Chrome dispose d’une extension websocket Simple Websocket Client qui permet de faire quelques tests de bon fonctionnement.

5.3.6 Utilisation avec Matlab

On initie une connexion websocket via MATLAB en utilisant par exemple MatlabWebSocket. Il faut également pouvoir encoder et décoder du JSON en MATLAB, par exemple en utilisant les fonctions MATLAB jsondecode et jsonencode.

5.3.7 Description des entrées/sorties pour une utilisation en “Model Exchange” (x -> dx/dt)

Dans ce mode, xdyn calcule uniquement la dérivée des 13 états navire mais n’effectue pas l’intégration numérique, ce qui permet d’utiliser un solveur externe, par exemple Matlab ou Simulink.

Chaque élément de type « État » est composé des éléments suivants:

Exemple d’entrée:

{
  "states": [
    {
      "t": 0,
      "x": 0,
      "y": 8,
      "z": 12,
      "u": 1,
      "v": 0,
      "w": 0,
      "p": 0,
      "q": 1,
      "r": 0,
      "qr": 1,
      "qi": 0,
      "qj": 0,
      "qk": 0
    },
    {
      "t": 0.1,
      "x": 0.1,
      "y": -0.156,
      "z": 10,
      "u": 0,
      "v": 0,
      "w": 0,
      "p": 0,
      "q": 0,
      "r": 0,
      "qr": 0,
      "qi": 0,
      "qj": 0,
      "qk": 0,
    }
  ],
  "commands": {
    "beta": 0.1
  }
}

La sortie du “Model Exchange” correspond à la dérivée des états à savoir dx/dt, dy/dt, dz/dt, du/dt, dv/dt, dw/dt, dp/dt, dq/dt, dr/dt, dqr/dt, dqi/dt, dqj/dt, dqk/dt.

Exemple de sortie:

{
  "dx_dt": 0.1,
  "dy_dt": -0.156,
  "dz_dt": 10,
  "du_dt": 0,
  "dv_dt": 0,
  "dw_dt": 0,
  "dp_dt": 0,
  "dq_dt": 0,
  "dr_dt": 0,
  "dqr_dt": 0,
  "dqi_dt": 0,
  "dqj_dt": 0,
  "dqk_dt": 0
}

La représentation textuelle de ces nombres flottants est faite de façon unique (et non exacte) : cela signifie que si les flottants sont différents (binairement), leur représentation textuelle sera différente. Cela signifie également que si xdyn lit la représentation textuelle et la traduit en binaire, on retrouvera bien la valeur binaire initiale. En d’autres termes, la fonction flottant -> texte est injective. Cela n’implique pas qu’elle soit bijective, puisque si l’on part d’une représentation textuelle, que l’on convertit en binaire pour reconvertir ensuite en texte on ne retrouvera pas nécessairement le texte initial.

5.3.8 Description des entrées/sorties pour une utilisation en “Co-Simulation” (x(t) -> [x(t), …,x(t+Dt)])

Chaque élément de type « État » est composé des éléments suivants:

Exemple d’entrée:

{
  "Dt": 10,
  "states": [
    {
      "t": 0,
      "x": 0,
      "y": 8,
      "z": 12,
      "u": 1,
      "v": 0,
      "w": 0,
      "p": 0,
      "q": 1,
      "r": 0,
      "qr": 1,
      "qi": 0,
      "qj": 0,
      "qk": 0
    },
    {
      "t": 0,
      "x": 0.1,
      "y": -0.156,
      "z": 10,
      "u": 0,
      "v": 0,
      "w": 0,
      "p": 0,
      "q": 0,
      "r": 0,
      "qr": 0,
      "qi": 0,
      "qj": 0,
      "qk": 0,
    }
  ],
  "commands": {
    "beta": 0.1
  }
}

La sortie est une liste d’éléments de type « État augmenté » contenant l’historique des états de t0 à t0 + Dt par pas de dt.

La signification exacte des angles d’Euler dépend de la convention d’angle choisie dans le fichier YAML d’entrée de xdyn (voir la section correspondante dans la documentation). Cette sortie est fournie pour faciliter le travail du client du serveur, mais n’est pas utilisée en interne par xdyn.

Exemple de sortie:

[
  {
    "t": 10,
    "x": 0,
    "y": 8,
    "z": 12,
    "u": 1,
    "v": 0,
    "w": 0,
    "p": 0,
    "q": 1,
    "r": 0,
    "qr": 1,
    "qi": 0,
    "qj": 0,
    "qk": 0,
    "phi": 0,
    "theta": 0,
    "psi": 0
  },
  {
    "t": 10.1,
    "x": 0.0999968,
    "y": 8,
    "z": 12.0491,
    "u": 0.897068,
    "v": 0,
    "w": 1.07593,
    "p": 0,
    "q": 1,
    "r": 0,
    "qr": 0.99875,
    "qi": 0,
    "qj": 0.0499792,
    "qk": 0,
    "phi": 0,
    "theta": 0.1,
    "psi": 0
  }
]

Comme pour le “model exchange”, la représentation textuelle de ces nombres flottants est faite de façon unique (et non exacte) : cela signifie que si les flottants sont différents (binairement), leur représentation textuelle sera différente. Cela signifie également que si xdyn lit la représentation textuelle et la traduit en binaire, on retrouvera bien la valeur binaire initiale. En d’autres termes, la fonction flottant -> texte est injective. Cela n’implique pas qu’elle soit bijective, puisque si l’on part d’une représentation textuelle, que l’on convertit en binaire pour reconvertir ensuite en texte on ne retrouvera pas nécessairement le texte initial.

6 Repères et conventions

6.1 Repères

Afin de connaître et décrire l’attitude d’un ou plusieurs corps dans l’espace, il est nécessaire de les placer par rapport à un repère de référence.

6.1.1 Repère de référence (NED)

Le repère NED (North-East-Down) est utilisé comme repère de référence, avec un point de référence $O$ et une base pointant les directions Nord-Est-Bas. Il sert à exprimer les déplacements des corps de la simulation.

6.1.2 Repère navire (mobile ou “body” ou “repère de résolution”)

Le repère navire correspond au repère attaché au navire lors de la simulation. Le point de référence de ce repère correspond généralement au centre de gravité du navire. Les axes du repère navire sont les suivants:

• $X$ vers l’avant
• $Y$ sur tribord
• $Z$ vers le bas

On résout les équations du mouvement (principe fondamental de la dynamique) à l’origine de ce repère, c’est-à-dire que tous les torseurs d’effort sont déplacés au point (0,0,0) du repère “body”. De même, les matrices d’inertie et de masse ajoutée sont déplacées au point de résolution.

6.1.3 Repère NED local

Lors des exports de champs de vagues, ni le repère NED ni le repère navire ne sont parfaitement adaptés : en effet, si le maillage sur lequel on calcule la houle est lié à NED, le navire finira par sortir de cette zone lorsqu’il se déplacera. Si l’on calcule les hauteurs de vague dans le repère navire, l’aire de la grille vue dans le repère NED va varier en fonction de l’attitude du navire et, à la limite, pour un navire vertical ($\theta=\pi/2$ par exemple), la projection de la grille est un segment.

On définit donc un NED “local”, c’est-à-dire un repère centré au même point que le repère navire mais simplement translaté par rapport à NED :

Ce repère est nommé “NED(body)”. Ainsi, si le navire s’appelle “nav1”, le repère NED local sera “NED(nav1)”.

6.1.4 Attitude navire

L’attitude d’un corps permet de connaître son orientation par rapport à un repère. La position est donnée par le triplet $\left(X,Y,Z\right)$ et l’orientation par un triplet d’angles $\left(\phi,\theta,\psi\right)$. L’interprétation de ce triplet en termes de rotations autour des axes $x$, $y$, $z$ dépend de la convention de rotation choisie. L’orientation peut également être exprimée de manière différente notamment avec des quaternions (c’est d’ailleurs ainsi qu’elle est exprimée dans le code d’xdyn).

6.2 Conventions d’orientations

Cette section présente les notations utilisées pour définir l’orientation d’un élément dans l’espace à partir d’un triplet d’angles $(\phi,\theta,\psi)$.

6.2.1 Définition d’une orientation

Pour définir la composition de rotations donnant l’orientation d’un élément dans l’espace à partir d’un triplet d’angles $\left(\phi,\theta,\psi\right)$, plusieurs éléments doivent être définis:

• une convention d’angles ou d’axes. Elle permet de définir si ce sont les angles ou les axes qui évoluent pour la notation utilisée.
• une composition des rotations interne ou externe. Elle définit si la composition des rotations s’effectue par rapport à un système d’axes fixes ou alors par rapport au système d’axes nouvellement modifiés.
• un ordre dans lequel est appliqué les rotations. Il permet de définir complètement la composition de rotations.

6.2.2 Énumération des conventions possibles

Si on choisit une convention d’angles, alors chaque angle du triplet définit respectivement une rotation autour d’un axe $X$, $Y$ ou $Z$. Les axes ne peuvent être répétés. Il est possible de définir 6 conventions d’angles, qui correspondent à la permutation des trois axes: $XYZ$ ,$XZY$ ,$YXZ$ ,$YZX$ ,$ZXY$ ,$ZYX$. Par exemple la rotation $R_{YZX}$ appliquée au triplet $\left(\phi,\theta,\psi\right)$ s’interprétera comme une rotation de $R_{Y}\left(\theta\right)$, suivie de la rotation $R_{Z}\left(\psi\right)$, et terminée par la rotation $R_{X}\left(\phi\right)$.

Si on choisit une convention d’axes, alors on modifie l’ordre des axes sur lesquels appliquer successivement les rotations. Des répétitions des axes sont alors possibles, si elles ne se suivent pas. Par exemple, $XYX$ sera valide, mais pas $XXY$. Par exemple, une convention ZXY définit une composition de rotations. Il est possible de définir 12 conventions d’axes: $XYX$, $XYZ$, $XZX$, $XZY$, $YXY$, $YXZ$, $YZX$, $YZY$, $ZXY$, $ZXZ$, $ZYX$, $ZYZ$. Par exemple la rotation $R_{YZX}$ appliquée au triplet $\left(\phi,\theta,\psi\right)$ s’interprétera comme une rotation de $R_{Y}\left(\phi\right)$, suivie de la rotation $R_{Z}\left(\theta\right)$, et terminée par la rotation $R_{X}\left(\psi\right)$.

Avec ces conventions d’angles et d’axes, il existe déjà 18 combinaisons. Ce nombre est doublé du fait que la composition de rotations peut être interne (intrinsic) ou externe (extrinsic). Si les rotations sont composées par rapport au repère fixe, on parle de composition externe. Si les rotations sont composées par rapport aux repères nouvellement créés, on parle de composition interne. C’est cette dernière qui est utilisée dans la majorité des cas. Au total, ce sont ainsi 36 conventions qu’il est possible de définir.

6.2.3 Synthèse des différentes conventions

Les deux tableaux suivants présentent les 36 conventions possibles :

où les matrices de rotation autour des trois axes $X$, $Y$ et $Z$ s’écrivent

$R_X \left( \alpha \right) = \left[\begin{array}{ccc} 1 & 0 & 0\\ 0 & +\cos \left( \alpha \right) & -\sin \left( \alpha \right)\\ 0 & +\sin \left( \alpha \right) & +\cos \left( \alpha \right) \end{array}\right]$

$R_Y \left( \alpha \right) = \left[\begin{array}{ccc} +\cos\left( \alpha \right) & 0 & +\sin \left( \alpha \right) \\ 0 & 1 & 0 \\ -\sin\left( \alpha \right) & 0 & +\cos \left( \alpha \right) \\ \end{array}\right]$

$R_Z \left( \alpha \right) = \left[\begin{array}{ccc} +\cos \left( \alpha \right) & -\sin \left( \alpha \right) & 0 \\ +\sin \left( \alpha \right) & +\cos \left( \alpha \right) & 0 \\ 0 & 0 & 1 \\ \end{array}\right]$

6.3 Conventions couramment utilisées

Parmi l’ensemble des conventions possibles, certaines sont plus utilisées que d’autres.

6.3.1 Convention des angles aéronautiques

La convention des angles aéronautiques (convention de type 2 de la norme AFNOR) exprimée par le triplet (Roulis, Tangage, Lacet) régulièrement utilisée est référencée id=12 dans le tableau ci-dessus.

Elle se comprend de la manière suivante, on effectue une rotation de l’angle de lacet autour de l’axe $Z$, suivie d’une rotation de l’angle d’assiette autour du nouvel axe $Y'$ suivie d’une rotation de l’angle de roulis autour du nouvel axe $X''$.

Si on exprime ce triplet de la manière suivante (Lacet, Roulis, Tangage), on obtient id=33 dans le tableau ci-dessus.

La composition de rotations de ces deux conventions est représentée sur la figure suivante

6.3.2 Convention ParaView

La convention d’orientation utilisée dans le logiciel ParaView est identifiée par id=11 dans le tableau ci-dessus.

Cette composition de rotation se comprend comme une rotation $\psi$ autour de l’axe $Z$, suivie d’une rotation $\theta$ autour du nouvel axe $X'$ et finalement d’une rotation $\phi$ autour du nouvel axe $Y''$.

6.4 Quaternions

L’utilisation des angles d’Euler pose deux problèmes principaux :

• Le blocage de Cardan, d’une part, qui se manifeste par la perte d’un degré de liberté lorsque les axes de deux des trois rotations définissant la transformation du repère NED au repère navire sont confondus,
• La périodicité des angles, d’autre part, qui introduit des discontinuités dans les états.

Une manière usuelle de contourner ces problèmes est d’utiliser des quaternions qui, au prix de l’ajout d’un état supplémentaire permettent de définir les rotations sans ambiguïté et de façon unique, quelle que soit la convention d’angle adoptée.

Lors des calculs de tenue à la mer (avec une formulation fréquentielle), on utilise souvent un repère linéarisé. Ce repère, qui peut être impliqué lors du lien avec les bases de données hydrodynamiques issues du fréquentiel, est calculé de la façon suivante. En faisant l’hypothèse que les mouvements sont faibles, on effectue un développement des rotations limité au premier ordre et ainsi elles peuvent être exprimées indépendamment par rapport aux axes principaux liés à la position moyenne du navire, dans n’importe quel ordre. Ce repère n’est pas utilisé dans la version actuelle d’xdyn.

6.5 États navires

Le simulateur est multi-corps en ce sens que plusieurs corps peuvent être simulés en même temps. Ainsi, on peut modéliser plusieurs corps mécaniquement indépendants, avec leurs interactions hydrodynamiques, pourvu que l’on implémente le modèle d’interaction (qui peut venir d’un fichier HDB multicorps). Actuellement, aucun effort d’interaction ni de liaison cinématique ne sont implémentés.

Chaque corps possède des états, permettant de reconstituer exactement son mouvement. Ces états sont les suivants :

• La position du corps par rapport à l’origine du NED projetée dans le repère corps est notée $p^n = [x,y,z]^T$ et est exprimée en mètres.
• L’attitude du corps est notée $\Theta = [\phi,\theta,\psi]^T$ et est définie au paragraphe précédent. En pratique, on utilise plutôt des quaternions $q = [q_r, q_i, q_j, q_k]$ en interne dans le simulateur pour l’intégration des équations du mouvement, mais les raisonnements présentés dans cette documentation se transposent aisément.
• La vitesse de translation du corps par rapport au repère fixe NED, projetée dans le repère du corps (ou “body”) est notée $v^b = [u,v,w]^T$ et s’exprime en m/s.
• Le vecteur vitesse de rotation du repère “body” par rapport au repère NED, projeté dans le repère “body”, est noté $\omega_{nb}^b = [p,q,r]^T$ et s’exprime en rad/s.
• Les efforts appliqués au navire et projetés dans le repère “body” sont notés : $f^b = [X,Y,Z]^T$. Ils s’expriment en N.
• Les moments appliqués au navire et projetés dans le repère “body” sont notés : $m^b = [K,M,N]^T$. Ils s’expriment en N.m.

6.6 Définition des corps simulés

La section bodies du fichier YAML contient une liste de corps, chacun commençant par un tiret :

bodies:
    - name: ....
      ...
    - name: ....
      ...

Chaque corps comprend :

• un nom (section name)
• éventuellement un maillage (section mesh)
• la position du repère body par rapport au maillage (section position of body frame relative to mesh)
• ses conditions initiales (sections initial position of body frame relative to NED et initial velocity of body frame relative to NED)
• des données nécessaires pour les calculs cinétiques (section dynamics)
• la liste des efforts auxquels il est soumis (sections external forces et controlled forces).
• de façon facultative, des états forcés.

6.6.1 Exemple complet

bodies: # All bodies have NED as parent frame
  - name: TestShip
    mesh: test_ship.stl
    position of body frame relative to mesh:
        frame: mesh
        x: {value: 9.355, unit: m}
        y: {value: 0, unit: m}
        z: {value: -3.21, unit: m}
        phi: {value: 0, unit: rad}
        theta: {value: 0, unit: rad}
        psi: {value: 0, unit: rad}
    initial position of body frame relative to NED:
        frame: NED
        x: {value: 0, unit: m}
        y: {value: 0, unit: m}
        z: {value: -5, unit: m}
        phi: {value: 0, unit: deg}
        theta: {value: -.0058, unit: rad}
        psi: {value: 0, unit: deg}
    initial velocity of body frame relative to NED:
        frame: TestShip
        u: {value: 0, unit: m/s}
        v: {value: 0, unit: m/s}
        w: {value: 0, unit: m/s}
        p: {value: 0, unit: rad/s}
        q: {value: 0, unit: rad/s}
        r: {value: 0, unit: rad/s}
    dynamics:
        hydrodynamic forces calculation point in body frame:
            x: {value: 0.696, unit: m}
            y: {value: 0, unit: m}
            z: {value: 1.418, unit: m}
        centre of inertia:
            frame: TestShip
            x: {value: 0.258, unit: m}
            y: {value: 0, unit: m}
            z: {value: 0.432, unit: m}
        mass: {value: 253.31, unit: tonne} # Caution: 'ton' is the british ton which is 907.185 kg
        rigid body inertia matrix at the center of gravity and projected in the body frame:
            row 1: [253310,0,0,0,0,0]
            row 2: [0,253310,0,0,0,0]
            row 3: [0,0,253310,0,0,0]
            row 4: [0,0,0,1.522e6,0,0]
            row 5: [0,0,0,0,8.279e6,0]
            row 6: [0,0,0,0,0,7.676e6]
        added mass matrix at the center of gravity and projected in the body frame:
            row 1: [3.519e4,0,0,0,0,0]
            row 2: [0,3.023e5,0,0,0,0]
            row 3: [0,0,1.980e5,0,0,0]
            row 4: [0,0,0,3.189e5,0,0]
            row 5: [0,0,0,0,8.866e6,0]
            row 6: [0,0,0,0,0,6.676e6]
    external forces:
      - model: gravity
      - model: non-linear hydrostatic (fast)
    blocked dof:
       from CSV:
         - state: u
           t: T
           value: PS
           interpolation: spline
           filename: test.csv
       from YAML:
         - state: p
           t: [4.2]
           value: [5]
           interpolation: piecewise constant

6.6.2 Nommage du solide

Le nom du solide a son importance puisqu’en définissant un solide, on définit implicitement le repère qui lui est attaché (le repère “body”, cf. documentation des repères).

On peut ensuite y faire référence, notamment pour les post-traitements.

6.6.3 Utilisation d’un maillage

Pour les efforts intégrés sur la carène (par exemple, les efforts hydrostatiques non-linéaires et les efforts de Froude-Krylov), il est nécessaire de définir un maillage.

6.6.3.1 Définition du fichier de maillage

La section mesh est optionnelle. Si l’on choisit de l’utiliser, elle doit contenir un nom de fichier STL contenant le maillage surfacique du navire. Les formats STL ASCII et STL binaire sont supportés. Ce chemin doit être donné relativement à l’endroit où l’on lance le simulateur.

Par exemple, si l’exécutable du simulateur est dans le répertoire A/B/C, que le maillage m.stl est dans A et qu’on lance l’exécutable depuis B, on écrira :

mesh: ../m.stl

6.6.3.2 Passage du repère maillage au repère body

L’origine du repère “body” (qui est le repère dans lequel est réalisé le bilan des efforts) est spécifiée par rapport au repère du maillage. En pratique, $x,y,z$ peuvent définir la position du centre de gravité dans le repère maillage et $\phi,\theta,\psi$ définissent la rotation permettant de passer du repère maillage au repère body (suivant la convention choisie dans la section rotations).

position of body frame relative to mesh:
    frame: mesh
    x: {value: 0, unit: m}
    y: {value: 0, unit: m}
    z: {value: -10, unit: m}
    phi: {value: 1, unit: rad}
    theta: {value: 3, unit: rad}
    psi: {value: 2, unit: rad}

6.6.4 Champs dynamics

La section dynamics permet de décrire l’inertie du solide. Elle est composée de cinq sous-sections :

• hydrodynamic forces calculation point in body frame est le point de calcul des efforts hydrodynamiques
• centre of inertia (si le repère “body” n’est pas au centre de masse)
• mass contenant la masse du corps considéré
• rigid body inertia matrix at the center of gravity and projected in the body frame définissant la matrice d’inertie
• added mass matrix at the center of gravity and projected in the body frame pour les masses ajoutées.

6.6.4.1 Position du centre d’inertie

La section centre of inertia a l’allure suivante :

centre of inertia:
    frame: TestShip
    x: {value: 0, unit: m}
    y: {value: 0, unit: m}
    z: {value: 0, unit: m}

Le champs frame doit contenir soit NED, soit un nom de solide, soit mesh(<nom de solide>) suivant le repère choisi pour exprimer les coordonnées du centre de gravité.

6.6.4.2 Définition de la masse du solide

La section mass contient la masse du solide, utilisée notamment pour le modèle de gravité :

mass: {value: 1000, unit: tonne}

Attention : si l’on spécifie ton, le système utilisera la tonne britannique qui vaut 907.185 kg.

6.6.4.3 Matrice d’inertie

La matrice d’inertie n’est pas normalisée et on n’effectue pas de changement de repère.

rigid body inertia matrix at the center of gravity and projected in the body frame:
    row 1: [253310,0,0,0,0,0]
    row 2: [0,253310,0,0,0,0]
    row 3: [0,0,253310,0,0,0]
    row 4: [0,0,0,1.522e6,0,0]
    row 5: [0,0,0,0,8.279e6,0]
    row 6: [0,0,0,0,0,7.676e6]

Les calculs sont réalisés en supposant que les termes $a_{ij}$ de la matrice d’inertie sont en unité SI, c’est-à-dire:

• en $\texttt{kg}$ pour les termes $a_{ij}$ avec $1\leq i,j \leq 3$,
• en $\texttt{kg}\times\texttt{m}^2$ pour les termes $a_{ij}$ avec $4\leq i,j \leq 6$,
• en $\texttt{kg}\times\texttt{m}$ pour les deux blocs 3x3 extra-diagonaux représentant les termes croisés de la matrice d’inertie.

6.6.4.4 Inerties ajoutées

La matrice de masse ajoutée est définie de la même façon :

added mass matrix at the center of gravity and projected in the body frame:
    row 1: [3.519e4,0,0,0,0,0]
    row 2: [0,3.023e5,0,0,0,0]
    row 3: [0,0,1.980e5,0,0,0]
    row 4: [0,0,0,3.189e5,0,0]
    row 5: [0,0,0,0,8.866e6,0]
    row 6: [0,0,0,0,0,6.676e6]

Elle figure dans la section dynamics et non dans la section external forces (bien qu’il s’agisse d’un modèle d’effort, proportionnel à l’accélération) parce que ce modèle d’effort fait l’objet d’un traitement particulier : il figure dans le membre de gauche de l’équation fondamentale de la dynamique

$$M\ddot{X} = \sum F_i$$(6.1) pour des raisons de stabilité numérique (l’effort dépend des accélérations qui doivent justement être calculées par la résolution de l’équation fondamentale de la dynamique). La matrice de masses ajoutées n’est cependant pas équivalente à une masse supplémentaire, car les termes de Coriolis et centripète qui correspondraient ne sont pas pris en compte.

Il est également possible d’extrapoler les masses ajoutées à pulsation infinie à partir d’un fichier HDB. Pour cela, on écrit (pour lire depuis le fichier test_ship.hdb) :

added mass matrix at the center of gravity and projected in the body frame:
    from hdb file: test_ship.hdb

Dans ce cas, il ne faut pas spécifier les clefs frame et row (le programme lance une exception si on le fait). Comme le fichier STL, le chemin du fichier HDB est relatif à l’endroit d’où on lance l’exécutable. La section correspondante dans le fichier HDB est Added_mass_Radiation_Damping. La valeur utilisée est la matrice de masse ajoutée à la période minimale définie dans le fichier HDB (aucune extrapolation n’est faite).

6.6.5 Forçage de degrés de liberté

Il est possible de forcer les valeurs des degrés de liberté suivant : U, V, W, P, Q, R pour chaque corps.

Pour forcer les degrés de liberté, on ajoute la section (facultative) suivante à la section body :

blocked dof:
   from CSV:
     - state: u
       t: T
       value: PS
       interpolation: spline
       filename: test.csv
   from YAML:
     - state: p
       t: [4.2]
       value: [5]
       interpolation: piecewise constant

Soit les états sont donnés directement dans le fichier YAML, soit ils sont lus depuis un fichier CSV.

• state: nom de l’état à forcer. u, v, w, p, q ou r.

Si les valeurs des états sont dans le YAML :

• value: Valeur de l’état forcé pour chaque instant (dans l’exemple ci-dessus, u=5 pour t $\geq$ 4.2).
• t: instants auxquels est défini l’état
• interpolation: type d’interpolation à réaliser. piecewise constant, linear ou spline.

Si les valeurs sont lues depuis un fichier CSV :

• filename: nom du fichier contenant les valeurs forcées. Le chemin du fichier s’entend relativement à l’endroit d’où est lancé le simulateur.
• value: nom de la colonne contenant les valeurs à lire
• t: nom de la colonne contenant le temps
• interpolation: type d’interpolation à réaliser. piecewise constant, linear ou spline.

Pour les valeurs de t hors de l’intervalle [tmin,tmax], l’état est supposé libre (non forcé).

Il est possible de récupérer dans les sorties l’écart entre l’effort réel et l’effort permettant de conserver les forçages, en d’autres termes il est possible de récupérer

$$(M+M_a)\dot{X_{\textrm{blocked}}} - \sum F_i$$(6.2)

Pour ce faire, on utilise dans la section ‘output’ les clefs suivantes (si le corps s’appelle ’TestShip):

• Fx(blocked states,TestShip,TestShip)
• Fy(blocked states,TestShip,TestShip)
• Fz(blocked states,TestShip,TestShip)
• Mx(blocked states,TestShip,TestShip)
• My(blocked states,TestShip,TestShip)
• Mz(blocked states,TestShip,TestShip)

Il est à noter que ces efforts sont exprimés dans le repère BODY.

6.7 Fichiers HDB issus d’un calcul fréquentiel

Le format HDB (Hydrodynamic DataBase) est le format standard du logiciel Diodore. Le logiciel AQUA+ (développé et utilisé en interne par l’École Centrale de Nantes et SIREHNA, dont des parties ont été reprises dans logiciel NEMOH) dispose également d’une sortie permettant de créer un fichier HDB. Ces fichiers peuvent être utilisés par xdyn pour calculer :

• les masses ajoutées (cf. paragraphe précédent)
• les amortissements de radiation
• les efforts de diffraction, calculés à partir de fonctions de transfert (RAO)

6.7.1 Conventions des fichiers HDB

Les fichiers HDB ne spécifient ni leur repère de calcul, ni les points d’expression. Plutôt que de laisser l’utilisateur la spécifier (avec les risques d’erreur que cela comporte) ou de la détecter automatiquement, on suppose pour xdyn que les fichiers HDB ont été générés avec les conventions du logiciel AQUA+ (convention “z vers le haut”). Le repère dans lequel sont exprimées toutes les matrices de tous les fichiers HDB lus par xdyn est donc : x longitudinal, y vers bâbord et z vers le haut.

Les RAO dépendent des conventions (repère, origine des phases, provenance vs. propagation, signe des angles). Or la convention utilisée ne figure pas dans le format HDB. Comme les RAO que l’on utilise actuellement proviennent exclusivement du logiciel AQUA+, dans la version actuelle xdyn suppose que la convention utilisée est celle d’AQUA+.

La différence entre la convention de xdyn et celle d’AQUA+ est que l’axe $z$ est ascendant pour AQUA+ et descendant pour xdyn. L’angle de la houle représente bien dans les deux cas une direction de propagation (et non de provenance). Le potentiel des vitesses de la houle est, dans xdyn comme dans AQUA+ :

$\phi(x,y,z,t) = -\frac{g\eta_a}{\omega}\frac{\cosh(k\cdot(h-z))} {\cosh(k\cdot h)}\cos(k\cdot(x\cdot \cos(\gamma)+ y\cdot \sin(\gamma))-\omega\cdot t+\phi)$

Toutes les données issues des fichiers HDB sont données en convention z vers le haut : par conséquent, il faut effectuer un changement de repère (rotation de $\pi$ autour de l’axe $X$) pour les mettre dans le repère d’xdyn (z vers le bas). La matrice de changement de base est :

$$R_X(\pi)=\left[\begin{array}{ccc} 1 & 0 &0\\0&-1&0\\0&0&-1\end{array}\right]$$(6.3)

• Pour les vecteurs

$$\left[\begin{array}{c}x_d\\y_d\\z_d\\\end{array}\right]_{\mbox{xdyn}}=\left[\begin{array}{ccc} 1 & 0 &0\\0&-1&0\\0&0&-1\end{array}\right]\left[\begin{array}{c}x_a\\y_a\\z_a\\\end{array}\right]_{\mbox{AQUA+}} =\left[\begin{array}{c}x_a\\-y_a\\-z_a\\\end{array}\right]_{\mbox{xdyn}}$$(6.4)

• Pour les matrices (masses ajoutées, amortissements…) Si $M=((m_{ij}))$ désigne une matrice exprimée dans le repère AQUA+ et $M_d$ la même matrice exprimée dans le repère xdyn, on a :

$$M_d = \left[\begin{array}{cc}R_X(\pi)&S(AB)R_X(\pi)\\0&R_X(\pi)\end{array}\right]^\top M \left[\begin{array}{cc}R_X(\pi)&S(AB)R_X(\pi)\\0&R_X(\pi)\end{array}\right]$$(6.5)

Par conséquent, toutes les matrices lues depuis le fichier HDB (masses ajoutées et amortissement de radiation) subissent un changement de repère décrit au paragraphe suivant.

6.7.2 Transport des matrices d’inertie et d’amortissement lues depuis le fichier HDB

La formule de Huygens permet de translater les inerties du centre de gravité vers un point quelconque. On s’intéresse ici au déplacement et au changement de repère d’une matrice d’inertie généralisée (inerties propres et inerties ajoutées). Par rapport à la formule de Huygens, la différence vient de ce que l’inertie n’est pas identique sur tous les axes (mathématiquement, la matrice peut donc être n’importe quelle forme bilinéaire symétrique, positive et définie).

Soit

$$S(\lambda)=\left[\begin{array}{ccc} 0&-\lambda_3&\lambda_2\\\lambda_3&0&-\lambda_1\\-\lambda_2&\lambda_1&0\end{array}\right]_{\mathcal{R_0}}$$(6.6)

la matrice du produit vectoriel par $\lambda^\top=[\lambda_1,\lambda_2,\lambda_3]$ : $\forall a\in\mathcal{R_0},\lambda\times a = S(\lambda) a$

Pour effectuer un transport d’une matrice 6x6 en coordonnées généralisées (masse, masse ajoutée ou amortissement) d’un point $A$ vers un point $B$, on utilise :

$$M_B=\left[\begin{array}{cc}I&S(AB)\\0&I\end{array}\right]^\top M_A \left[\begin{array}{cc}I&S(AB)\\0&I\end{array}\right]$$(6.7)

Cette formule est une généralisation de la formule de Huygens.

En combinant avec un changement de base (de la base $a$ vers la base $b$) par la matrice de rotation ${}^a R_b$ de $b$ vers $a$ on obtient l’expression plus générale :

$${}^bM_B=\left[\begin{array}{cc}{}^a R_b&S(AB){}^a R_b\\0&{}^a R_b\end{array}\right]^\top {}^aM_A \left[\begin{array}{cc}{}^a R_b&S(AB){}^a R_b\\0&{}^a R_b\end{array}\right]$$(6.8)

Cette formule permet d’effectuer à la fois le transport d’une matrice d’inertie généralisée 6x6 d’un point $A$ à un point $B$ et le changement de son repère d’expression de $a$ vers $b$.

La matrice d’inertie est le plus souvent exprimée au centre de gravité. Cependant, rien ne le garantit dans les fichiers HDB. Dans le cas général, il faut donc effectuer un transport, mais on suppose dans xdyn qu’il n’y a pas de transport à faire et donc que la seule opération pour la conversion est le passage d’un repère z vers le haut (HDB) à un repère z vers le bas (xdyn). La formule précédente se simplifie alors en :

$${}^bM_B=\left[ \begin{array}{cccccc} 1&0&0&0&0&0\\ 0&-1&0&0&0&0\\ 0&0&-1&0&0&0\\ 0&0&0&1&0&0\\ 0&0&0&0&-1&0\\ 0&0&0&0&0&-1\\ \end{array} \right]^\top {}^aM_A \left[ \begin{array}{cccccc} 1&0&0&0&0&0\\ 0&-1&0&0&0&0\\ 0&0&-1&0&0&0\\ 0&0&0&1&0&0\\ 0&0&0&0&-1&0\\ 0&0&0&0&0&-1\\  \end{array} \right]$$(6.9)

En prenant la notation

$$M_{\mbox{HDB}} = \left[ \begin{array}{cccccc} m_{11} & m_{12} & m_{13} & m_{14} & m_{15} & m_{16}\\ m_{21} & m_{22} & m_{23} & m_{24} & m_{25} & m_{26} \\ m_{31} & m_{32} & m_{33} & m_{34} & m_{35} & m_{36} \\ m_{41} & m_{42} & m_{43} & m_{44} & m_{45} & m_{46}\\ m_{51} & m_{52} & m_{53} & m_{54} & m_{55} & m_{56} \\ m_{61} & m_{62} & m_{63} & m_{64} & m_{65} & m_{66} \\ \end{array} \right]$$(6.10)

on obtient :

$$M_{\mbox{xdyn}} = \left[ \begin{array}{cccccc} m_{11} & -m_{12} & -m_{13} & m_{14} & -m_{15} & -m_{16}\\ -m_{21} & m_{22} & m_{23} & -m_{24} & m_{25} & m_{26} \\ -m_{31} & m_{32} & m_{33} & -m_{34} & m_{35} & m_{36} \\ m_{41} & -m_{42} & -m_{43} & m_{44} & -m_{45} & -m_{46}\\ -m_{51} & m_{52} & m_{53} & -m_{54} & m_{55} & m_{56} \\ -m_{61} & m_{62} & m_{63} & -m_{64} & m_{65} & m_{66} \\ \end{array} \right]$$(6.11)

Pour les torseurs, on obtient le changement suivant :

$$\tau_{\mbox{HDB}} = \left[ \begin{array}{c} F_X\\ F_Y\\ F_Z\\ M_X\\ M_Y\\ M_Z\\ \end{array} \right]$$(6.12)

$$\tau_{\mbox{xdyn}} = \left[ \begin{array}{c} F_X\\ -F_Y\\ -F_Z\\ M_X\\ -M_Y\\ -M_Z\\ \end{array} \right]$$(6.13)

Cf. SimBody Theory Manual, Release 3.1, March, 2013, page 137, §12.3.1, Rigid body shift of rigid body spatial inertia

6.7.3 Structure des fichiers HDB

Les fichiers HDB sont des sorties générées par un code Fortran 77. Aucune spécification formelle de ce format n’a été trouvée. Pour pouvoir les importer dans xdyn, leur grammaire EBNF a été déduite de l’analyse de plusieurs fichiers (reverse-engineering). Cette grammaire, qui permet de parser tous les fichiers HDB qui nous ont été fournis jusqu’ici, peut être représentée comme suit (sans garantie d’exhaustivité) :

ast                             = string-key
                                | value-key
                                | vector-section
                                | matrix-section
                                | list-of-matrix-sections
                                | list-of-matrix-sections-with-id;
str                             = character , +extended-char , *" "  , +extended-char;
header                          = "[" , str , "]";
string-key                      = header , str , -eol;
value-key                       = header , double , -eol;
values                          = double_ , +double , -eol;
vector-section                  = header , eol , +(double , eol) , -eol;
matrix-section                  = header , eol , +(*(eol, values)) , -eol;
list-of-matrix-sections         = header , eol , +matrix-section , -eol;
section-with-id                 = header , double , eol , *(eol, values) , -eol;
list-of-matrix-sections-with-id = header , eol , +(header , double , eol , +(*(eol, values)));
character                       = "A" | "B" | "C" | "D" | "E" | "F" | "G"
                                | "H" | "I" | "J" | "K" | "L" | "M" | "N"
                                | "O" | "P" | "Q" | "R" | "S" | "T" | "U"
                                | "V" | "W" | "X" | "Y" | "Z" | "a" | "b"
                                | "c" | "d" | "e" | "f" | "g" | "h" | "i"
                                | "j" | "k" | "l" | "m" | "n" | "O" | "p"
                                | "q" | "r" | "s" | "t" | "u" | "v" | "w"
                                | "x" | "y" | "z" | "_";
extended-char                   = character | "-" | "+"

Les fichiers HDB sont donc constitués d’une liste d’éléments, ces éléments étant pouvant être de type string-key, value-key, vector-section, matrix-section, list-of-matrix-sections ou list-of-matrix-sections-with-id. xdyn ne tient pas compte de l’ordre dans lequel sont ces éléments (mais aucune garantie ne peut être fournie pour d’éventuels autres outils utilisant les HDB).

6.7.3.1 Exemple de “string-key”

 [SOFT] AQUA+

6.7.3.2 Exemple de “value-key”

[FORWARD_SPEED]   0.00

6.7.3.3 Exemple de “vector-section”

 [List_calculated_periods]
   4.00
  64.00
 125.00

6.7.3.4 Exemple de “matrix-section”

[Mass_Inertia_matrix]
 8.635000E+06  0.000000E+00  0.000000E+00  0.000000E+00  0.000000E+000.000000E+00
 0.000000E+00  8.635000E+06  0.000000E+00  0.000000E+00  0.000000E+000.000000E+00
 0.000000E+00  0.000000E+00  8.635000E+06  0.000000E+00  0.000000E+000.000000E+00
 0.000000E+00  0.000000E+00  0.000000E+00  4.149000E+08  0.000000E+000.000000E+00
 0.000000E+00  0.000000E+00  0.000000E+00  0.000000E+00  1.088000E+100.000000E+00
 0.000000E+00  0.000000E+00  0.000000E+00  0.000000E+00  0.000000E+00  1.088000E+10

6.7.3.5 Exemple de “list-of-matrix-sections”

[Added_mass_Radiation_Damping]
[ADDED_MASS_LINE_1]
   4.00  8.770898E+04  0.000000E+00  2.122348E+05  0.000000E+00  1.844677E+07  0.000000E+00
  64.00  2.356385E+05  0.000000E+00  3.063730E+05  0.000000E+00  6.103379E+07  0.000000E+00
 125.00  2.334437E+05  0.000000E+00  3.058729E+05  0.000000E+00  6.044752E+07  0.000000E+00
[ADDED_MASS_LINE_2]
   4.00  0.000000E+00  2.406631E+06  0.000000E+00 -1.298266E+06  0.000000E+00  2.894931E+07
  64.00  0.000000E+00  7.130066E+06  0.000000E+00 -3.516670E+06  0.000000E+00  9.536960E+07
 125.00  0.000000E+00  7.100760E+06  0.000000E+00 -3.493418E+06  0.000000E+00  9.499156E+07

6.7.3.6 Exemple de “list-of-matrix-sections-with-id”

[FROUDE-KRYLOV_FORCES_AND_MOMENTS]
[INCIDENCE_EFM_MOD_001]   0.000000
  4.00  6.452085E+04  0.000000E+00  3.775068E+05  0.000000E+00  2.630894E+07  0.000000E+00
 64.00  8.296234E+04  0.000000E+00  2.098381E+07  0.000000E+00  1.280202E+08  0.000000E+00
125.00  2.179682E+04  0.000000E+00  2.105635E+07  0.000000E+00  1.258710E+08  0.000000E+00
[INCIDENCE_EFM_MOD_001]   30.00000
  4.00  5.988292E+04  3.453055E+04  5.220861E+05  6.375514E+05  2.533077E+07  1.149621E+06
 64.00  7.185571E+04  4.148640E+04  2.098683E+07  7.927488E+04  1.274491E+08  3.394846E+04
125.00  1.887675E+04  1.089865E+04  2.105658E+07  2.081770E+04  1.258309E+08  2.938749E+03
[INCIDENCE_EFM_PH_001]   0.000000
  4.00  5.079494E-01  1.570796E+00  1.171722E+00  1.570796E+00  1.426003E+00  1.570796E+00
 64.00 -3.141362E+00  1.570796E+00 -1.576680E+00  1.570796E+00 -1.768797E+00  1.570796E+00
125.00 -3.141476E+00  1.570796E+00 -1.572334E+00  1.570796E+00 -1.623406E+00  1.570796E+00

6.7.4 Sections utilisées par xdyn

Les sections du HDB actuellement utilisées par xdyn sont :

• Added_mass_Radiation_Damping pour les amortissements de radiation et les matrices de masses ajoutées,
• DIFFRACTION_FORCES_AND_MOMENTS pour les RAO des efforts de diffraction.

6.8 Repère de calcul hydrodynamique

Les efforts d’amortissement visqueux et de résistance à l’avancement sont calculés dans un repère appelé repère de calcul hydrodynamique, qui est un repère translaté par rapport au repère body. Le centre de ce repère est un point défini (dans le repère body) de la façon suivante :

• Son abscisse $x$ est celle du centre de la surface résultant de la projection du maillage sur le plan $(x,z)$
• Son ordonnée $y$ vaut zéro
• Son altitude $z$ est celle du centre de la surface résultant de la projection du maillage sur le plan $(x,y)$

Ce point est, en général, distinct du centre de gravité et du centre de volume. Il est défini dans la section dynamics/hydrodynamic forces calculation point in body frame du fichier YAML :

hydrodynamic forces calculation point in body frame:
    x: {value: 0.696, unit: m}
    y: {value: 0, unit: m}
    z: {value: 1.418, unit: m}

On note ${}^{\textrm{local}}T_{\textrm{body}}$ la transformation permettant de convertir des coordonnées dans le repère body en coordonnées du même point exprimées dans le repère de calcul hydrodynamique. ${}^{\textrm{local}}T_{\textrm{NED}}$ est celle permettant de convertir des coordonnées dans le repère NED en coordonnées du même point exprimées dans le repère de calcul hydrodynamique.

Il convient de distinguer ce repère de celui utilisé dans la base de données hydrodynamiques (fichiers HDB), utilisé pour l’expression des matrices d’amortissement de radiation, les RAO d’effort (pour le calcul des efforts de diffraction) et les masses ajoutées.

6.9 Rotations et constantes environnementales

6.9.1 Rotations

La convention de rotation permet de retrouver la matrice de rotation d’un repère par rapport à un autre, étant donné un triplet $(\phi, \theta, \psi)$. La convention utilisée est décrite dans la section rotations :

rotations convention: [psi,theta',phi'']

Cette ligne s’interprète de la façon suivante : étant donné un triplet $(\phi, \theta, \psi)$, on construit les matrices de rotation en effectuant d’abord une rotation d’angle $\psi$ autour de l’axe Z, ensuite une rotation d’angle $\theta$ autour de l’axe Y du repère précédemment transformé, puis une rotation d’angle $\phi$ autour de l’axe X du repère ainsi obtenu.

Si l’on avait noté :

rotations convention: [z,y',x'']

on aurait d’abord une rotation d’angle $\phi$ autour de l’axe Z, puis une rotation d’angle $\theta$ autour du nouvel axe Y, puis une rotation $\psi$ autour du nouvel axe X.

Des apostrophes sont utilisées pour indiquer des compositions de rotations par rapport au nouveau système d’axes, et donc une composition interne. Ainsi [x,y',z''] désignera une rotation autour X, suivie d’une rotation autour du nouvel axe Y, appelé Y’ et terminée par une rotation autour du nouvel axe Z, appelé Z’’. La double apostrophe fait référence au deuxième repère utilisée pour la composition de rotation.

La liste rotations convention comporte toujours trois éléments. Le deuxième élément est toujours différent du premier. Le troisième est soit différent des deux autres, soit égal au premier.

La convention des angles aéronautiques, fréquemment (et abusivement) dénotée “angles d’Euler” (lacet $\psi$, tangage $\theta$, roulis $\phi$), se définit de la façon suivante :

rotations convention: [psi, theta', phi'']

La convention d’orientation utilisée dans le logiciel ParaView est donnée par :

rotations convention: [psi,phi',theta'']

Pour plus de détails sur les conventions d’angles et d’axes, se référer à la documentation détaillée.

La définition d’un repère dans xdyn s’effectue à partir d’un repère connu (NED ou body). On définit le changement de coordonnées pour passer d’un repère connu à un nouveau repère de la manière suivante dans le fichier YAML :

• frame le nom du repère connu,
• x ,y ,z: le triplet de position où chaque position est définie par le dictionnaire des clés value et unit,
• phi ,theta ,psi, le triplet d’orientation dont l’interprétation en termes de matrices de rotations dépend de la convention d’orientation choisie

    frame: NED
    x: {value: 0, unit: m}
    y: {value: 0, unit: m}
    z: {value: 0, unit: m}
    phi: {value: 0, unit: rad}
    theta: {value: 0, unit: rad}
    psi: {value: 0, unit: rad}

Cette description est en particulier utilisée pour définir le changement de coordonnées entre le repère du maillage et le repère du corps (par la clef position of body frame relative to mesh) et la position initiale du corps à simuler (clef initial position of body frame relative to NED).

6.10 Convention de houle

xdyn utilise la convention $Z$ vers le bas pour l’amplitude. L’azimuth ou la direction de houle que l’on note $\gamma$ représente la direction de propagation de la houle.

Ainsi une direction nulle indique une propagation du sud vers le nord de la houle. Une direction de 90° représente une propagation de l’ouest vers l’est.

Les différentes positions du navire par rapport à la houle sont décrites comme suit :

• $\gamma - \psi = 0^\circ$ : navire dos à la houle, houle de l’arrière
• $0^\circ < \gamma - \psi < 90^\circ$ : Houle venant de bâbord arrière
• $\gamma - \psi = 90^\circ$ : Houle venant de bâbord
• $90^\circ < \gamma - \psi < 180^\circ$ : Houle venant de tribord arrière
• $\gamma - \psi = 180^\circ$: navire cap à la houle, houle de l’avant
• $180^\circ < \gamma - \psi < 270^\circ$ : Houle venant de tribord avant
• $\gamma - \psi = 270^\circ$ : Houle venant de tribord
• $270^\circ < \gamma - \psi < 360^\circ$ : Houle venant de bâbord avant

7 Modèles d’environnement

Les modèles d’environnement sont les modèles de houle (et, à terme, de vent, de courant…) utilisés par xdyn. Actuellement, seuls des modèles de houle sont implémentés. Leur paramétrisation figure dans la section environment du fichier YAML d’entrée. Elle peut être vide (par exemple, lors de la simulation simple du tutoriel 1).

Les modèles de houle interviennent pour le calcul des efforts hydrostatiques non-linéaires (par le truchement de l’élévation de la surface libre), les efforts de Froude-Krylov (par le biais de la pression dynamique), le calcul des efforts de diffraction (via le fréquentiel) et le modèle de safran (prise en compte des vitesses orbitales).

7.1 Constantes environnementales

L’accélération de la pesanteur (dénotée par g), la densité volumique de l’eau (rho) et sa viscosité nu sont des constantes qui interviennent dans plusieurs modèles physiques. Par conséquent, plutôt que d’être renseignées au niveau de chaque modèle et risquer ainsi des incohérences, elles figurent dans la section environmental constants qui a la forme suivante :

environmental constants:
    g: {value: 9.81, unit: m/s^2}
    rho: {value: 1025, unit: kg/m^3}
    nu: {value: 1.18e-6, unit: m^2/s}

Ces trois constantes sont l’ensemble de toutes les constantes environnementales actuellement utilisées par les modèles d’xdyn.

Comme expliqué dans une section précédente, les dimensions physiques ne sont pas vérifiées et simplement converties en unités du système international. Si xdyn rencontre une unité inconnue, il produit un message d’erreur du type :

unknow unit : hhm

7.2 Simulation sans houle

Pour simuler une surface libre parfaitement plane, on opère de la façon suivante :

    environment:
      - model: no waves
        constant sea elevation in NED frame: {value: 0, unit: m}

model: no waves indique que l’on souhaite une surface libre horizontale et constant sea elevation in NED frame représente l’élévation de la surface libre dans le repère NED.

Dans ce cas, les efforts d’excitation (Froude-Krylov et radiation) seront nuls.

7.3 Houle d’Airy

On peut définir une houle comme étant une somme de plusieurs spectres directionnels, c’est-à-dire un spectre de puissance et une dispersion spatiale. Pour dériver l’expression générale d’une houle composée de plusieurs spectres, on commence par le cas d’une houle monochromatique et monodirectionnelle.

7.3.1 Expression du potentiel de vitesse de la houle

Soit $V(x,y,z,t)=(u,v,w)$ la vitesse du fluide au point de coordonnées $(x,y,z)$ (dans le repère NED) et à l’instant $t$.

On suppose l’eau non visqueuse, incompressible, homogène et isotrope et l’on considère un écoulement irrotationnel. Supposer l’écoulement irrotationnel implique (d’après le lemme de Poincaré) que la vitesse dérive d’un potentiel que l’on appelle $\phi:(x,y,z,t)\mapsto\phi(x,y,z,t)$. Par définition, la vitesse en tout point de l’écoulement est donc donnée par :

$$V(x,y,z,t) = \textrm{grad}{\phi(x,y,z,t)}$$(7.1)

La pression $p$ vérifie l’équation de Bernoulli :

$$p + \rho g z -\rho\frac{\partial\phi}{\partial t} + \frac{\rho}{2} V\cdot V = C(t)$$(7.2)

où $C:t\mapsto C(t)$ est une fonction du temps arbitraire, donc l’équation est en particulier valable pour $C(t)=p_0$ (pression atmosphérique à la surface) :

$$p_0 + \rho g z -\rho\frac{\partial\phi}{\partial t} + \frac{\rho}{2} V\cdot V = p_0$$(7.3)

soit

$$g z -\frac{\partial\phi}{\partial t} + \frac{1}{2} V\cdot V = 0$$(7.4)

Le terme $\rho g z$ représente la pression hydrostatique et le terme $-\rho\frac{\partial\phi}{\partial t}$ est la pression dynamique.

Il s’agit de la première condition de surface libre.

On peut définir la fonction $F(x,y,z,t)=z-\eta(x,y,t)$

Pour une particule sur la surface libre, $F(x,y,z,t)=0$ ce qui implique que sa dérivée particulaire est nulle :

$$\frac{DF}{Dt} = \frac{\partial F}{\partial t} + V\cdot \nabla F = 0$$(7.5)

soit

$\frac{\partial\eta}{\partial t} + \frac{\phi}{\partial x}\frac{\partial \eta}{\partial x} + \frac{\phi}{\partial x}\frac{\partial\eta}{\partial y} - \frac{\partial\phi}{\partial z} = 0$ sur $z=\eta(x,y,t)$.

C’est la deuxième condition de surface libre.

En linéarisant ces deux conditions de surface libre, on obtient :

$$\frac{\partial \eta}{\partial t} = \frac{\partial\phi}{\partial z}$$(7.6)

$$g\eta - \frac{\partial\phi}{\partial t} = 0$$(7.7)

Par ailleurs, l’eau étant supposée incompressible, $\nabla\cdot V= \frac{\partial^2\phi}{\partial x^2} + \frac{\partial^2\phi}{\partial y^2} + \frac{\partial^2\phi}{\partial z^2} = 0$

Il s’agit d’une équation de Laplace dont la solution s’obtient par la méthode de séparation des variables. Plusieurs potentiels peuvent être solution. Par exemple :

• $\phi(x,y,z,t) = -\frac{g\eta_a}{\omega}\frac{\cosh(k\cdot(h-z))} {\cosh(k\cdot h)}\cos(k\cdot(x\cdot \cos(\gamma)+ y\cdot \sin(\gamma))-\omega\cdot t+\phi)$
• $\phi(x,y,z,t) = \frac{g\eta_a}{\omega}\frac{\cosh(k\cdot(h-z))} {\cosh(k\cdot h)}\sin(k\cdot(x\cdot \cos(\gamma)+ y\cdot \sin(\gamma))-\omega\cdot t+\phi)$

Ici, nous choisissons :

$$\phi(x,y,z,t) = -\frac{g\eta_a}{\omega}\frac{\cosh(k\cdot(h-z))} {\cosh(k\cdot h)}\cos(k\cdot(x\cdot \cos(\gamma)+ y\cdot \sin(\gamma))-\omega\cdot t+\phi)$$(7.8)

qui est le potentiel utilisé par le logiciel AQUA+.

• $h$ est la profondeur du fluide (hauteur du sol à la surface libre),
• $\eta_a$ est l’amplitude de la houle (en m),
• $x,y,z$ sont les coordonnées du point considéré, exprimées dans le repère NED,
• $k$ est le nombre d’onde, traduisant la périodicité spatiale.

7.3.2 Relation entre le nombre d’onde et la pulsation

Il existe une relation entre $k$ et $\omega$, appelée relation de dispersion, et qui s’écrit :

$$\omega^2 = g\cdot k \cdot \tanh(k\cdot h)$$(7.9)

où $h$ désigne la profondeur d’eau et $g$ l’accélération de la pesanteur.

En profondeur infinie ($k\cdot h > 3$), cette relation tend vers :

$$\omega^2 \sim g\cdot k$$(7.10)

7.3.3 Élévation de la houle

L’élévation de la houle découle de la deuxième condition de surface libre :

$\eta(x,y,t) = +\frac{1}{g} \frac{\partial\phi(x,y,z=0,t)}{\partial t} =-\sum_{i=1}^{nfreq} A(\omega_i) \sin(k_i\cdot(x\cdot \cos(\gamma)+ y\cdot \sin(\gamma))-\omega_i\cdot t+\phi_{i})$

où $\gamma$ désigne la direction de provenance de la houle, définie à la section convention de houle.

L’élévation $\eta$ de la houle en un point $(x,y)$ est un signal temporel $\eta(t)$. Par définition du plan $z=0$ (élévation moyenne de la houle), ce signal $\eta$ est centré. On le suppose également stationnaire, ce qui implique que sa fonction d’auto-corrélation $R$ ne dépend que de $\tau$:

$$R(\tau)=\mathbf{E} (\eta(t)\eta(t+\tau))$$(7.11)

La densité spectrale de puissance $G$ de $\eta$ est égale à la transformée de Fourier de sa fonction d’auto-corrélation $R$. Comme $R$ est paire et réelle, $G$ aussi et on ne considère usuellement que la partie positive (one-sided) en définissant la densité spectrale $S$ telle que :

$$S(\omega)=\left\{\begin{array}{cc}2G(\omega),&\omega\geq 0\\0,&\omega<0\end{array}\right.$$(7.12)

On a :

$$R(\tau) = \frac{1}{2} \sum_{i=1}^{nfreq} A(\omega_i)^2 \cos(\omega_i \tau)$$(7.13)

or on a également :

$$R(\tau) = \int_0^\infty S(\omega) \cos(\omega\tau)d\omega$$(7.14)

d’où, par identification :

$$A(\omega_i)^2=2 S(\omega_i)d\omega_i$$(7.15)

On peut généraliser cette formulation en faisant intervenir l’étalement directionnel $D(\gamma)$ de la houle :

$$A(\omega,\gamma)^2 = 2 S(\omega)d\omega D(\gamma) d\gamma$$(7.16)

On obtient, en définitive :

$$\eta(x,y,t) = -\sum_{i=1}^{nfreq} \sqrt{2 S(\omega_i)d\omega D(\gamma) d\gamma} \sin(k_i\cdot(x\cdot \cos(\gamma)+ y\cdot \sin(\gamma))-\omega_i\cdot t+\phi_{i})$$(7.17)

7.3.4 Pression dynamique

L’expression de la pression dynamique (champs de pression de la houle incidente), utilisée par le modèle de Froude-Krylov, est définie comme la pression totale moins la pression hydrostatique $\rho g z$ et son expression se déduit de la première condition de surface libre linéarisée :

$$p_{\textrm{dyn}} = -\rho \frac{\partial \Phi(x,y,z,t)}{\partial t}$$(7.18)

soit

$$p_{\textrm{dyn}} = \rho\cdot g \sum_{i=1}^{nfreq}A(\omega_i,\gamma) \frac{\cosh(k_i\cdot(h-z))}{\cosh(k_i\cdot h)}\sin(k_i\cdot(x\cdot \cos(\gamma)+ y\cdot \sin(\gamma))-\omega_i\cdot t+\phi_{i})$$(7.19)

• $g$ désigne l’accélération de la pesanteur (9.81 $m/s^2$)
• $\rho$ est la densité volumique du fluide (en $kg/m^3$)