3. Transient anechoic fluid elements#

This part presents most of the general constraints for implementing anechoic fluid elements with absorbent boundary elements with the paraxial approximation of order 0 in code_aster. For reasons of simplicity related to the manipulation of scalar quantities such as pressure or displacement potential, as opposed to vector quantities such as displacement, we are first interested in fluid elements.

3.1. Standard formulation#

Here we repeat Modaressi’s reasoning by adapting it to a fluid acoustic domain. First, we are interested in the sole data of the pressure quantity in this fluid. We will then come back to this modeling to adapt to code_aster constraints, highlighting the adjustments to be made.

So let’s look at the following configuration, using the conventions of the previous part in the vicinity of the border:

_images/Shape1.gif

Defining a local coordinate system at the element level allows us to systematically bring us back into such a situation.

3.1.1. Finite element formulation#

Pressure \(p\) verifies the Helmholtz equation throughout domain \(\Omega\) modelled with finite elements, which gives, for any virtual pressure field \(q\):

\(-\underset{\Omega }{\int }\nabla p\text{.}\nabla q-\frac{1}{{c}^{2}}\frac{{\partial }^{2}}{\partial {t}^{2}}\underset{\Omega }{\int }\text{pq}+\underset{\Sigma }{\int }\frac{\partial p}{\partial n}q=0\)

\(\Sigma\) represents the border of the \(O\) domain.

The quantity to be estimated on \(\Sigma\) using the paraxial approximation is here \(\frac{\partial p}{\partial n}\).

3.1.2. Paraxial approximation#

In the proposed configuration, the term \(\frac{\partial p}{\partial n}\) corresponds to \(\frac{\partial p}{\partial {x}_{3}}\).

Let us therefore consider a harmonic plane wave propagating in the fluid:

\(p=A\text{exp}\left[i({k}_{1}{x}_{1}+{k}_{2}{x}_{2}+{k}_{3}{x}_{3}-\omega t)\right]\)

Substituting in the Helmholtz equation, we get:

\({k}_{3}=\frac{\omega }{c}\sqrt{1-\frac{{c}^{2}}{{\omega }^{2}}({k}_{1}^{2}+{k}_{2}^{2})}\)

The following development is then obtained, for high frequencies (\(\omega\) large) or in the vicinity of the border (\({k}_{1}\) and \({k}_{2}\) small):

\({k}_{3}=\frac{\omega }{c}(1-\frac{{c}^{2}}{2{\omega }^{2}}({k}_{1}^{2}+{k}_{2}^{2}))\)

Or, by multiplying by \(\omega\) to make this quantity disappear from the denominator and after an inverse Fourier transform in space and time:

\(\frac{{\partial }^{2}p}{\partial {x}_{3}\partial t}\text{=}\text{-}\frac{1}{c}\frac{{\partial }^{2}p}{\partial {t}^{2}}+\frac{1}{2}c(\frac{{\partial }^{2}p}{\partial {x}_{1}^{2}}+\frac{{\partial }^{2}p}{\partial {x}_{2}^{2}})\)

As presented by Modaressi, this equation involves the derivative with respect to time of the surface term. In the context of this part, we are only interested in the term of order 0, i.e., after an integration in time, this makes the annoying derivative disappear:

\(\frac{\partial p}{\partial {x}_{3}}\text{=}\text{-}\frac{1}{c}\frac{\partial p}{\partial t}\) or more generally: \(\frac{\mathrm{\partial }p}{\mathrm{\partial }n}\text{=}\text{-}\frac{1}{c}\frac{\mathrm{\partial }p}{\mathrm{\partial }t}\)

It is this impedance relationship that we will discretize on the border of the finite element domain.

Note:

Given the disappearance of the term of order 1 in the development of the square root, the minimum order of approximation for fluid paraxials is in fact 1 and not 0. We will keep the name elements of order 0 for consistency with the solid. However, we are talking about second-order fluid elements when considering strictly positive elements.

3.2. Impedance of vibro-acoustic elements in code_aster#

code_aster has vibro-acoustic elements. This paragraph recalls the formulation choices made at the time of their implementation. We are inspired to present the existing reference documentation r4.02.02.

3.2.1. Limits of the p formulation#

In the context of fluid-structure interaction in harmony, the pressure-only formulation of the acoustic fluid leads to non-symmetric matrices. In fact, the global system is expressed, in variational form, in the following way:

\(\underset{{\Omega }_{s}}{\int }{C}_{\text{ijkl}}\text{.}{u}_{k,l}{v}_{i,j}-{\omega }^{2}\underset{{\Omega }_{s}}{\int }{\rho }_{s}{u}_{i}{v}_{i}-\underset{\Sigma }{\int }p{v}_{i}\text{.}{n}_{i}=0\) for the structure

\(\frac{1}{{\rho }_{f}{\omega }^{2}}\underset{{\Omega }_{f}}{\int }\nabla p\text{.}\nabla q-{k}^{2}\underset{{\Omega }_{f}}{\int }\mathrm{pq}-\underset{\Sigma }{\int }{u}_{i}\text{.}{n}_{i}q=0\) for the fluid

with \(k=\frac{\omega }{c}\), wave number for the fluid, \(v\) and \(q\) two virtual fields in the structure and in the fluid respectively.

After discretization by finite elements, the following matrix system is obtained:

\(\left[\begin{array}{cc}K& -C\\ 0& H\end{array}\right]\left[\begin{array}{c}u\\ p\end{array}\right]-{\omega }^{2}\left[\begin{array}{cc}M& 0\\ {\rho }_{f}{C}^{T}& \frac{Q}{{c}^{2}}\end{array}\right]\left[\begin{array}{c}u\\ p\end{array}\right]=0\)

where:

\(K\) and \(M\) are the stiffness and mass matrices of the structure

\(H\) and \(Q\) are the fluid matrices obtained respectively from the bilinear forms: \(\underset{{\Omega }_{f}}{\int }\nabla p\text{.}\nabla q\) and \(\underset{{O}_{f}}{\int }\mathrm{pq}\)

\(C\) is the coupling matrix obtained from the bilinear form: \(\underset{\Sigma }{\int }p{u}_{i}\text{.}{n}_{i}\)

The non-symmetric nature of this system does not allow the use of the classic code_aster resolution algorithms. This motivates the introduction of an additional variable in the description of the fluid.

3.2.2. Symmetric formulation in p and phi#

The new magnitude introduced is the potential for \(\phi\) trips, such as \(u\text{=}\nabla \varphi\). From this, the new variational form of the fluid-structure coupled system is obtained:

\(\underset{{\Omega }_{s}}{\int }{C}_{\text{ijkl}}\text{.}{u}_{k,l}{v}_{i,j}-{\omega }^{2}\underset{{\Omega }_{s}}{\int }{\rho }_{s}{u}_{i}{v}_{i}-{\rho }_{f}{\omega }^{2}\underset{\Sigma }{\int }\phi p{v}_{i}\text{.}{n}_{i}=0\) for the structure

\(\frac{1}{{\rho }_{f}{c}^{2}}\underset{{\mathrm{\Omega }}_{f}}{\int }\text{pq}-{\rho }_{f}{\omega }^{2}\left[\frac{1}{{\rho }_{f}{c}^{2}}\underset{{\mathrm{\Omega }}_{f}}{\int }\left(\varphi q+p\chi \right)-\underset{{\mathrm{\Omega }}_{f}}{\int }\nabla \varphi \text{.}\nabla \chi +\underset{S}{\int }\chi {\mathrm{u}}_{\mathrm{i}}{\mathrm{n}}_{\mathrm{i}}\right]=0\) for the fluid

With: \(p={\rho }_{f}{\omega }^{2}f\) in the fluid and \(\chi\) a field of potential for virtual displacement

This leads us to the symmetric matrix system:

\(\left[\begin{array}{ccc}K& 0& 0\\ 0& \frac{{M}_{f}}{{\rho }_{f}{c}^{2}}& 0\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}u\\ p\\ \phi \end{array}\right]-{\omega }^{2}\left[\begin{array}{ccc}M& 0& {\rho }_{f}{M}_{\Sigma }\\ 0& 0& \frac{{M}_{\text{fl}}}{{c}^{2}}\\ {\rho }_{f}{M}_{\Sigma }^{T}& \frac{{M}_{\text{fl}}^{T}}{{c}^{2}}& {\rho }_{f}H\end{array}\right]\left[\begin{array}{c}u\\ p\\ \phi \end{array}\right]=0\)

where: \(K\) and \(\mathrm{M}\) are the stiffness and mass matrices of the structure

\({\mathrm{M}}_{\mathrm{\Sigma }}\) is the coupling matrix obtained from the bilinear form \(\underset{\mathrm{\Sigma }}{\int }\varphi {\mathrm{u}}_{\mathrm{i}}{\mathrm{n}}_{\mathrm{i}}\)

\({\mathrm{M}}_{\mathrm{f}},{\mathrm{M}}_{\text{fl}}\) and \(\mathrm{H}\) are the fluid matrices obtained from the bilinear forms: \(\underset{{\mathrm{\Omega }}_{f}}{\int }\text{pq}\), \(\underset{{\mathrm{\Omega }}_{f}}{\int }pq\) (or \(\underset{{\mathrm{\Omega }}_{f}}{\int }\varphi q\)) and \(\underset{{\mathrm{\Omega }}_{f}}{\int }\nabla \varphi \text{.}\nabla \chi\)

3.2.3. Formulation in psi#

There is a third formulation defined in \(\psi\), that is to say a speed potential. For more details, refer to the documentation. This formulation is also available in code_aster. Its negative point is that it does not have direct access to post-treatment under fluid pressure.

3.2.4. Imposition of an impedance of order 0#

In general, an impedance relationship at the fluid boundary is expressed as follows:

\(p=Z\mathrm{v}\text{.}\mathrm{n}\)

where:

  • \(Z=\rho c{q}_{\alpha }\) is the imposed impedance (\(\rho\) is the density of the fluid, \(c\) is the speed of the fluid waves and \({q}_{\alpha }\) is the reflection coefficient);

  • \(\mathrm{v}\text{.}\mathrm{n}\) is the normal outgoing speed of fluid particles.

We deduce from this, according to the law of fluid behavior, which relates pressure to the movement of fluid particles for an acoustic fluid \(\nabla p-{\rho }_{f}\frac{{\partial }^{2}u}{\partial {t}^{2}}=0\):

\(\frac{{\rho }_{f}}{Z}\frac{\partial p}{\partial t}=\frac{\partial p}{\partial n}\)

The discretization of such an equation leads to a non-symmetric term in a formulation in \(p\) and \(\varphi\). It is preferable to formulate the condition in relation to the displacement potential, i.e.:

\(\nabla \varphi +\frac{{\rho }_{f}}{Z}\frac{\partial \varphi }{\partial t}=0\)

As an expression for the edge term associated with the impedance relationship, we then obtain:

\({\rho }_{f}\frac{{\partial }^{2}}{\partial {t}^{2}}\underset{\mathrm{\Sigma }}{\int }\varphi \frac{\partial \psi }{\partial n}=\frac{{\partial }^{3}}{\partial {t}^{3}}\underset{\mathrm{\Sigma }}{\int }\frac{{\rho }_{f}^{2}}{Z}\varphi \psi\)

We then note the appearance (somewhat artificial) of a term in third derivative with respect to time. In harmonics, this is not a problem. You can easily treat a term in \({\omega }^{3}\). This third-order term will be calculated using option IMPE_MECA in CALC_MATR_ELEM (or ASSEMBLAGE) and used as an entry for MATR_IMPE_PHI in DYNA_VIBRA (TYPE_CALCUL = “HARM”). For the transient calculation, rather than introducing an approximation of a third derivative into the Newmark schema implemented in the dynamic direct integration operators in code_aster DYNA_VIBRA (or DYNA_LINE) and DYNA_NON_LINE, we can transform the third-order term into a first-order term and treat this as a classical damping term (for more details refer to). This term will be calculated using option AMOR_MECA in CALC_MATR_ELEM (or ASSEMBLAGE). This option is available for all three fluid formulations (for formulation details \(u-\psi\) refer to).

3.2.5. Imposition of an impedance of order 1#

The radiation condition is exact for a plane wave. However, to asymptotically approach the infinity behavior of a spherical wave, for a given radius \(R\), a first-order approximation is used. This condition is recalled in, as well as the complete description of the terms imposed by this condition. Depending on the formulation chosen, we will impact the mass or the stiffness of the model:

  • If we choose the formulation \(u-p-\varphi\), we add a mass term calculated using the MASS_MECA option in CALC_MATR_ELEM (or ASSEMBLAGE);

  • If we choose the formulations \(u-p\) or \(u-\psi\), we add a stiffness term calculated using the option RIGI_MECA in CALC_MATR_ELEM (or ASSEMBLAGE).

Radius \(R\) will be given via option LONG_CARA in DEFI_MATERIAU (FLUIDE…). This parameter is zero by default (in this case the first-order impedance is not taken into account).

3.2.6. The reflection rate#

If we want to represent a reflection rate \(\alpha\), for example in the case of the presence of alluvium at the bottom of the reservoir, we must then correct the speed value \(c\) by a new value \(c\text{'}\) such as \(c\text{'}=c\frac{1+\alpha }{1-\alpha }\). Parameter \(\alpha\) will be set via option COEF_AMOR in DEFI_MATERIAU (FLUIDE…). This parameter is zero by default (in this case the border is perfectly absorbent). This parameter influences both the order 0 and the 1 order impedance.

3.2.7. Detailed formulation#

Here we propose the precise formulation for an acoustic fluid modelled on a domain \(\mathrm{\Omega }\) with an anechoic condition on a part \({\mathrm{\Sigma }}_{a}\) of the domain border \(\mathrm{\Sigma }\). Apart from that, the border is broken down into a free surface and a part in contact with a rigid solid. The introduction of external stresses or the presence of an elastic structure is easily modelled by current methods. Volume and surface elements are formulated in \(p\) and \(\varphi\).

The equations in the fluid are:

\({\rho }_{f}\mathrm{\Delta }\varphi +\frac{1}{{c}^{2}}p=0\) in volume \(\mathrm{\Omega }\) eq 3.2.4-1

\(p={\rho }_{f}\frac{{\partial }^{2}\varphi }{\partial {t}^{2}}\) in volume \(\mathrm{\Omega }\) eq 3.2.4-2

\(p=0\) on the free surface eq 3.2.4-3

\(\frac{\partial \varphi }{\partial n}=0\) on the rigid wall eq 3.2.4-4

\(\frac{\partial p}{\partial n}\text{=}\text{-}\frac{1}{c}\frac{\partial p}{\partial t}\) on the part of the border with anechoic condition eq 3.2.4-5

We multiply the equation [éq 3.2.4-1] by a virtual potential field \(y\) and we integrate into \(\mathrm{\Omega }\):

\(\underset{{\mathrm{\Omega }}_{f}}{\int }\left[\frac{1}{{c}^{2}}p\psi +{\rho }_{f}\frac{{\partial }^{2}}{\partial {t}^{2}}\left(\nabla \varphi \text{.}\nabla \psi \right)\right]+\underset{\mathrm{\Sigma }}{\int }\psi {\rho }_{f}\frac{{\partial }^{2}}{\partial {t}^{2}}\left(\frac{\partial \varphi }{\partial n}\right)=0\) according to Green’s formula

Or, with the boundary conditions on \(\mathrm{\Sigma }\) and equation [éq 3.2.4-2]:

\(\underset{{\mathrm{\Omega }}_{f}}{\int }\left[\frac{1}{{c}^{2}}p\psi +{\rho }_{f}\frac{{\partial }^{2}}{\partial {t}^{2}}\left(\nabla \varphi \text{.}\nabla \psi \right)\right]+\underset{\mathrm{\Sigma }}{\int }\psi {\rho }_{f}\frac{\partial p}{\partial n}=0\)

It is therefore possible to apply the impedance condition formulated under pressure:

\(\underset{{\mathrm{\Sigma }}_{a}}{\int }\psi {\rho }_{f}\frac{\partial p}{\partial n}\text{=}-\frac{1}{c}\underset{{\mathrm{\Sigma }}_{a}}{\int }\psi {\rho }_{f}\frac{\partial p}{\partial t}\)

In addition, to arrive at a symmetric formulation of the volume terms, we multiply the equation [éq3.2.4-2] by a virtual pressure field \(q\) and we integrate into \(\mathrm{\Omega }\):

\(\underset{{\mathrm{\Omega }}_{f}}{\int }\frac{\text{pq}}{{\rho }_{f}{c}^{2}}-\frac{{\partial }^{2}}{\partial {t}^{2}}\underset{{\mathrm{\Omega }}_{f}}{\int }\frac{\varphi q}{{c}^{2}}=0\)

Summarizing the two variational equations, we get:

\(\frac{1}{{\rho }_{f}{c}^{2}}\underset{{\mathrm{\Omega }}_{f}}{\int }\text{pq}+{\rho }_{f}\frac{{\partial }^{2}}{\partial {t}^{2}}\left[\frac{1}{{\rho }_{f}{c}^{2}}\underset{{\mathrm{\Omega }}_{f}}{\int }\left(\varphi q+p\psi \right)-\underset{{\mathrm{\Omega }}_{f}}{\int }\nabla \varphi \text{.}\nabla \psi \right]-\frac{1}{c}\underset{{\mathrm{\Sigma }}_{a}}{\int }\psi {\rho }_{f}\frac{\partial p}{\partial t}=0\)

Matricially:

\(\left[\begin{array}{cc}{\mathrm{M}}_{\mathrm{f}}& 0\\ 0& 0\end{array}\right]\left[\begin{array}{c}p\\ \mathrm{\varphi }\end{array}\right]-\frac{1}{c}\left[\begin{array}{cc}0& 0\\ \mathrm{A}& 0\end{array}\right]\left[\begin{array}{c}\dot{p}\\ \dot{\mathrm{\varphi }}\end{array}\right]+\left[\begin{array}{cc}0& \frac{{\mathrm{M}}_{\text{fl}}}{{c}^{2}}\\ \frac{{\mathrm{M}}_{\text{fl}}^{\mathrm{T}}}{{c}^{2}}& {\mathrm{\rho }}_{f}\mathrm{H}\end{array}\right]\left[\begin{array}{c}\ddot{p}\\ \ddot{\mathrm{\varphi }}\end{array}\right]=0\)

where sub-matrices \({\mathrm{M}}_{\mathrm{f}},{\mathrm{M}}_{\text{fl}}\) and \(\mathrm{H}\) discretize the same previous bilinear forms.

The \(\mathrm{A}\) sub-matrix discredits the term \(\underset{{\mathrm{\Sigma }}_{a}}{\int }\psi {\rho }_{f}\frac{\partial p}{\partial t}\). It should be noted that the damping matrix obtained is not symmetric.

3.3. Usage in code_aster#

Taking into account anechoic fluid elements and calculating their impedance requires specific modeling on absorbent boundaries:

  • in 2D with “2D_ FLUI_ABSO “or “AXIS_FLUI_ABSO” modeling on the absorbing edges at \(n\) knots;

  • in 3D with the “3D_ FLUI_ABSO “modeling on the absorbent faces with \(n\) knots.

In harmonic analysis for the formulation \(u-p-\varphi\) with the operator DYNA_VIBRA, a mechanical impedance is calculated beforehand by the option IMPE_MECA of the operator CALC_MATR_ELEM and it is entered in DYNA_LINE_HARM (keyword MATR_IMPE_PHI). In all other cases, the 0th-order impedance will be taken into account via option AMOR_MECA of the CALC_MATR_ELEM operator. The first-order impedance will be taken into account via stiffness or mass terms (depending on the formulation) as explained above. The material characteristics of fluid absorbent borders will be defined in DEFI_MATERIAU (FLUIDE…).