Digital implementation
=======================




Solving the Laplace equation by volume finite elements
----------------------------------------------------------------

Let's go back to the fluid Laplace problem with von Neumann-type boundary conditions:

:math:`\{\begin{array}{c}\begin{array}{c}\mathrm{\Delta }p=0\phantom{\rule{6em}{0ex}}\text{dans}\phantom{\rule{2em}{0ex}}{\mathrm{\Omega }}_{f}\\ (\frac{\partial p}{\partial \mathrm{n}}{)}_{{\mathrm{\Gamma }}_{1}}\text{=}-{\rho }_{f}\ddot{{\mathrm{x}}_{s}}\mathrm{.}\mathrm{n}\phantom{\rule{6em}{0ex}}\text{sur}\phantom{\rule{2em}{0ex}}{\mathrm{\Gamma }}_{1},\phantom{\rule{4em}{0ex}}{\mathrm{\Gamma }}_{1}=\phantom{\rule{2em}{0ex}}\underset{l=1,K}{\cup }{\gamma }_{l}^{\mathit{fs}}\end{array}\\ (\frac{\partial p}{\partial \mathrm{n}}{)}_{{\mathrm{\Gamma }}_{2}}=0\phantom{\rule{4em}{0ex}}\text{sur}\phantom{\rule{2em}{0ex}}{\mathrm{\Gamma }}_{2},\phantom{\rule{4em}{0ex}}{\mathrm{\Gamma }}_{2}=\partial \mathrm{\Omega }-{\mathrm{\Gamma }}_{1}\end{array}`

Let's write a variational formulation of this problem:

:math:`{\int }_{{\mathrm{\Omega }}_{f}}v\mathrm{.}\mathrm{\Delta }pd\mathrm{\Omega }=0\phantom{\rule{2em}{0ex}}\forall v`

Using Green's formula with a normal that is assumed to be oriented from the structure to the fluid (therefore inside the fluid volume) and by asking :math:`\Gamma \mathrm{=}{\Gamma }_{1}\mathrm{\cup }{\Gamma }_{2}`:

:math:`{\int }_{{\mathrm{\Omega }}_{f}}\mathrm{grad}v\mathrm{.}\mathrm{grad}pd\mathrm{\Omega }+{\int }_{\mathrm{\Gamma }}v\frac{\partial p}{\partial \mathrm{n}}\mathrm{.}\mathrm{n}dS=0`

Either:


.. _RefEquation 4.1-1:

:math:`{\int }_{{\mathrm{\Omega }}_{f}}\mathrm{grad}v\mathrm{.}\mathrm{grad}pd\mathrm{\Omega }={\rho }_{f}{\int }_{\mathrm{\Gamma }}v{\ddot{x}}_{n}dS` eq 4.1-1

We consider a partition of fluid volume :math:`{\mathrm{\Omega }}_{f}` into a finite number of elements. On this discretization of the domain, we can write an approximate form of the hydrodynamic pressure field:

:math:`{p}_{}\mathrm{=}\mathrm{\sum }_{i\mathrm{=}1}^{N}{N}_{i}(\mathrm{r}){p}_{i}`

:math:`{N}_{i}` represents the nodal interpolation functions defined on the elements: they are equal to 1 at the node number :math:`i`, and :math:`0` on all the others for finite elements P1.

Then, taking the nodal interpolation functions as test functions :math:`v` in succession, we obtain a system of :math:`N` equations by transferring to [:ref:`éq 4.1-1 <éq 4.1-1>`]:

:math:`j=1,\mathrm{...},N;\phantom{\rule{2em}{0ex}}{\int }_{{\mathrm{\Omega }}_{f}}\sum _{i=1}^{N}{p}_{i}\mathrm{grad}{N}_{i}\left(\mathrm{r}\right)\mathrm{.}\mathrm{\Delta }{N}_{j}\left(\mathrm{r}\right)d\mathrm{\Omega }={\rho }_{f}{\int }_{{\gamma }_{l}^{\mathit{fs}}}{N}_{j}{\ddot{x}}_{n}\mathit{dS}`

which can be written in the form of the linear system, whose solutions are the degrees of freedom of pressure :math:`\mathrm{P}`:


.. csv-table::

    ":math:`\mathrm{H}\mathrm{.}\mathrm{P}=\mathrm{\Phi }` with :math:`\mathrm{\Phi }` component vector :math:`{\mathrm{\Phi }}_{j}={\rho }_{f}\underset{\mathrm{\Gamma }}{\int }{N}_{j}{\ddot{x}}_{n}\mathit{dS}`"


.. _RefEquation 4.1-2                       |:

| with:math: `\ mathrm {H}` component matrix:math: `{H} _ {\ text {ij}} =\ underset {{\ mathrm {\ Omega}} _ {f}} {f}} {\ f}}} {\ int} {\ int}\ mathrm {int}}\ mathrm {.} \ mathrm {grad} {N} _ {j} _ {j} d\ mathrm {\ Omega} `eq 4.1-2 |
+---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+

Strictly speaking, this system is unique. It admits an infinity of solutions that differ from one constant. It is therefore necessary to impose pressure (Dirichlet-type boundary condition) at a point in the fluid to remove the uncertainty about the solution.

These precautions taken, by inverting the [:ref:`éq 4.1-2 <éq 4.1-2>`] system, we obtain the pressure field throughout the :math:`{\mathrm{\Omega }}_{f}` volume of fluid, including at the fluid/structure interface, where we are obviously interested in it.




Calculation of the coefficients of the added mass matrix on a modal basis
------------------------------------------------------

The value of the integral must be estimated numerically:


.. _RefEquation 4.2-1:

:math:`{m}_{\mathit{il}\text{jk}}=\underset{{\gamma }_{l}^{\mathit{fs}}}{\int }{p}_{\text{jk}}{\mathrm{X}}_{\mathit{il}}\mathrm{.}\mathrm{n}\mathit{dS}` eq 4.2-1

from a field with pressure nodes represented by a column vector noted :math:`{\mathrm{P}}_{\mathit{jk}}` and from a field with movement nodes corresponding to a modal deformation of an air structure and represented by the column vector :math:`{\mathrm{X}}_{\mathit{il}}`.

However, on the fluid/structure interface, the approximate pressure field :math:`{p}_{\mathrm{jk}}` due to the discretization of the interface into :math:`N` edge elements can be written:

:math:`{p}_{\text{jk}}={\sum ^{N}}_{m=1}{N}_{m}\left(\mathrm{r}\right){p}_{{\text{jk}}_{m}}`

while the "modal" displacement field is written on this same discretization:

:math:`{\mathrm{X}}_{\mathit{il}}=\sum _{n=1}^{N}{N}_{n}\left(\mathrm{r}\right){\mathrm{X}}_{{\mathit{il}}_{n}}`

So, by bringing these two expressions into the integral [:ref:`éq 4.2-1 <éq 4.2-1>`], we get:

:math:`\begin{array}{c}{m}_{\mathit{il}\text{jk}}\simeq \underset{{\gamma }_{l}^{\mathit{fs}}}{\int }(\sum _{m=1}^{N}{N}_{m}\left(\mathrm{r}\right){p}_{{\text{jk}}_{m}})\left[\begin{array}{c}\sum _{n=1}^{N}{N}_{n}\left(\mathrm{r}\right){X}_{{\mathit{ilx}}_{n}}\mathrm{{\rm A}}{n}_{x}+\sum _{n=1}^{N}{N}_{n}\left(\mathrm{r}\right){X}_{{\mathit{ily}}_{n}}\mathrm{{\rm A}}{n}_{y}\end{array}\right]\mathit{dS}\\ {m}_{\mathit{il}\text{jk}}\simeq \sum _{m=1}^{N}\sum _{n=1}^{N}{p}_{{\text{jk}}_{m}}(\underset{{\gamma }_{l}^{\mathit{fs}}}{\int }{N}_{m}\left(\mathrm{r}\right){N}_{n}\left(\mathrm{r}\right){n}_{x}\mathit{dS}){X}_{{\mathit{ilx}}_{n}}+\sum _{m=1}^{N}\sum _{n=1}^{N}{p}_{{\text{jk}}_{m}}(\underset{{\gamma }_{l}^{\mathit{fs}}}{\int }{N}_{m}\left(\mathrm{r}\right){N}_{n}\left(\mathrm{r}\right){n}_{y}\mathit{dS}){X}_{{\mathit{ily}}_{n}}\end{array}`

It was assumed in the above demonstration that the problem is two-dimensional.

This can be in the form of a dot product, involving a vector matrix product:


+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|:math:`{m}_{\mathit{il}\text{jk}}={\mathrm{P}}_{\text{jk}}^{T}{\mathrm{A}}_{x}{\mathrm{X}}_{\mathit{ilx}}+{\mathrm{P}}_{\text{jk}}^{T}{\mathrm{A}}_{y}{\mathrm{X}}_{\mathit{ily}}` with :math:`{\mathrm{A}}_{x}` coefficient matrix :math:`\underset{{\gamma }_{l}^{\mathit{fs}}}{\int }{N}_{i}{N}_{j}{n}_{x}\mathit{dS}`|
+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
|and :math:`{\mathrm{A}}_{y}` coefficient matrix :math:`\underset{{\gamma }_{l}^{\mathit{fs}}}{\int }{N}_{i}{N}_{j}{n}_{y}\mathit{dS}`                                                                                                                                                                                    |
+-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+