3. Theoretical formulation of the problem#
The cumbersome structure-fluid interaction problem amounts to solving three problems simultaneously:
the structure is subject to a pressure field \(P\) imposed by the fluid on the wall \(\Sigma\);
the fluid is subject to a displacement field \({X}_{s}\) imposed by the structure on \(\Sigma\);
gravity acts on the free surface where \(p=\rho gz\).
We will first consider that the fluid is not heavy before introducing gravity in paragraph [§3.2].
3.1. Reminders on fluid-structure coupling#
In order to fully account for fluid-structure interaction, we will analyze separately the equations governing the behavior of the fluid and those that govern that of the structure, without considering in this chapter the boundary conditions concerning the free surface.
3.1.1. Fluid description#
It is considered that the system under study is subject to small disturbances around its state of equilibrium where the fluid and the structure are at rest: thus, \(P\mathrm{=}{p}_{0}+p\) and \({\mathrm{X}}_{\mathrm{s}}\mathrm{=}{\mathrm{x}}_{\mathrm{s}}({x}_{0}\mathrm{=}0)\). This allows [bib2] to be written:
\(\rho \mathrm{=}\mathrm{-}{\rho }_{f}\text{div}({\mathrm{x}}_{\mathrm{f}})\) from where \(p\mathrm{=}\mathrm{-}{\rho }_{f}{c}^{2}\text{div}({\mathrm{x}}_{\mathrm{f}})\).
With:
\(p\) the fluctuating pressure of the fluid,
\(\mathrm{\rho }\) the density disturbance of the fluid, \({\mathrm{\rho }}_{f}\) the density of the fluid at rest,
\({x}_{f}(r,t)\) the field of movement of a fluid particle.
The fluid is:
perfect (i.e. not viscous)
barotropic:
\(p=\rho {c}^{2}\) eq 3.1.1-1
and irrotational: there is a potential for \(\varphi\) trips, such as \(p={\rho }_{f}\frac{{\partial }^{2}\varphi }{\partial {t}^{2}}\)
The behavior of the volume of Eulerian fluid is therefore described by the following equations:
law of behavior:
\({T}_{\text{ij}}\text{=}\text{-}p{\delta }_{\text{ij}}\)
equation for the conservation of momentum in the fluid in the absence of a source:
\(\text{div}(\overline{\overline{T}})={\rho }_{f}\frac{{\partial }^{2}{x}_{f}}{\partial {t}^{2}}\) eq 3.1.1-2
mass conservation equation:
\(\frac{\partial \rho }{\partial t}+{\rho }_{f}\text{div}(\frac{\partial {x}_{s}}{\partial t})=0\) eq 3.1.1-3
By combining the equations for the conservation of momentum [éq 3.1.1-2] and mass [éq 3.1.1-3] written in harmonic regime at pulsation \(w\), we obtain, thanks to [éq 3.1.1-1], the Helmholtz equation:
\(\Delta p+\frac{{\omega }^{2}}{{c}^{2}}p=0\)
3.1.2. Structure description#
It is considered that the structure is linear elastic and that we remain in the domain of small disturbances. Taking these assumptions into account, we write:
the law of behavior in linear elasticity:
\({\sigma }_{\text{ij}}\text{=}{C}_{\text{ijkl}}{\varepsilon }_{\text{kl}}\)
the equation for the conservation of momentum in the structure in the absence of volume forces other than inertia forces:
\(\text{div}(\overline{\overline{\sigma }})={\rho }_{s}\frac{{\partial }^{2}{x}_{s}}{\partial {t}^{2}}\)
the compatibility equation on the deformation tensor:
\({\varepsilon }_{\text{kl}}=\frac{1}{2}({u}_{k,l}+{u}_{l,k})\)
3.1.3. Description of the fluid-structure interface#
At the interface (\(\Sigma\)) between the fluid and the structure, as the fluid is not viscous, there is continuity of normal stresses and normal speeds at the wall, and there is no tangential stress (absence of viscous friction). These boundary conditions are written as:
\(\{\begin{array}{c}{\sigma }_{\text{ij}}{n}_{i}={T}_{\text{ij}}{n}_{i}\text{=}\text{-}p{\delta }_{\text{ij}}{n}_{i}\\ \frac{\partial {x}_{f}}{\partial t}\text{.}n=\frac{\partial {x}_{s}}{\partial t}\text{.}n\end{array}\)
3.1.3.1. Formulation of the coupled problem#
Finally, the equation of the coupled fluid-structure problem is written, taking \(p\) as a variable describing the pressure field in the fluid and \({x}_{s}\) the field of displacements in the structure:
\(\{\begin{array}{cc}{C}_{\text{ijkl}}{x}_{{s}_{k,\mathrm{lj}}}+{\omega }^{2}{\rho }_{s}{x}_{{s}_{i}}=0& \text{dans}{V}_{s}\\ \Delta p+\frac{{\omega }^{2}}{{c}^{2}}p=0& \text{dans}{V}_{f}\\ {\sigma }_{\text{ij}}{n}_{i}={C}_{\text{ijkl}}{x}_{{s}_{k,l}}{n}_{i}\text{=}\text{-}p{\delta }_{\text{ij}}{n}_{i}& \text{sur}\Sigma \\ \frac{\partial p}{\partial n}={\rho }_{f}{\omega }^{2}{x}_{{f}_{i}}{n}_{i}& \text{sur}\Sigma \end{array}\)
The displacement fields \({x}_{s}\) for the structure and pressure \(p\) for the desired fluid minimize the functional:
\(L({x}_{s},p,z)=\frac{1}{2}\underset{{V}_{s}}{\int }\left[{\sigma }_{\text{ij}}({x}_{s}){\varepsilon }_{\text{ij}}({x}_{s})-{\rho }_{s}{\omega }^{2}{x}_{s}^{2}\right]-\underset{\Sigma }{\int }p{x}_{s}nd\Sigma +\frac{1}{2{\rho }_{f}}\underset{{V}_{f}}{\int }\left[\frac{1}{{\omega }^{2}}(\text{grad}p{)}^{2}-\frac{{p}^{2}}{{c}^{2}}\right]\text{dV}\)
3.2. Action of gravity on the free surface#
3.2.1. Formulation of the problem#
We recall here the linearized dynamic equations describing the small movements of a perfect fluid. An Eulerian description of the fluid is chosen:
\(\text{grad}P\mathrm{=}{\rho }_{f}(\mathrm{g}\mathrm{-}\frac{{\mathrm{\partial }}^{2}{\mathrm{x}}_{\mathrm{s}}}{\mathrm{\partial }{t}^{2}})\) in \({V}_{f}\)
At equilibrium the fluid particle was in \({M}_{0}\) and therefore: \(\text{grad}{P}_{0}\mathrm{=}{\rho }_{f}\mathrm{g}\) in \({V}_{f}\).
We consider movements of low amplitude around the state of equilibrium (this is the hypothesis of small disturbances): so \(M\text{=}{M}_{0}\text{+}{x}_{f}({M}_{0},t)\)
Let \(p\) be the Eulerian pressure fluctuation and \({p}_{L}\) the Lagrangian pressure fluctuation, then:
\(\begin{array}{}p(M,t)=P({M}_{0},t)-{P}_{0}({M}_{0})\\ {p}_{L}=P(M,t)-{P}_{0}({M}_{0})\end{array}\)
Taking into account the hypothesis of small trips:
\(\begin{array}{}{p}_{L}-p=\text{gradP}({M}_{o},t){x}_{f}({M}_{o},t)\\ \phantom{{p}_{L}-p}\text{=}\text{-}{\rho }_{f}g{x}_{f}({M}_{o},t)\end{array}\) eq 3.2.1-1
If we consider the case of a heavy fluid having a free surface in contact with a medium at constant pressure \({P}_{\text{atm}}\), we can write, while neglecting the effects of surface tension:
\(P(M,t)={P}_{\text{atm}}\) on the free surface \(\text{SL}\) i.e.: \({p}_{L}=0\). That is, with [éq 3.2.1-1], \(p={\rho }_{f}g({x}_{f}\text{.}z)\)
Taking into account the hypothesis of small movements, the instantaneous inclination of the tangential plane is an infinitely small of the first order. So \({x}_{f}\text{.}z\) merges to the second order with the vertical elevation \(h\).
Figure 3.2.1-a: approximation on the free area
Thus, if we add to the boundary conditions the condition of gravity on the free surface, this is the same as considering the linearized condition in \(z=H\):
\(p={\rho }_{f}gz\) eq 3.2.1-2
The equations of the global problem become:
\(\{\begin{array}{cc}{C}_{\text{ijkl}}{x}_{{s}_{k,\text{lj}}}+{\omega }^{2}{\rho }_{s}{x}_{{s}_{i}}=0& \text{dans}{V}_{s}\\ \Delta p+\frac{{\omega }^{2}}{{c}^{2}}p=0& \text{dans}{V}_{f}\\ {\sigma }_{\text{ij}}{n}_{i}={C}_{\text{ijkl}}{x}_{{s}_{k,l}}{n}_{i}\text{=}\text{-}p{\delta }_{\text{ij}}{n}_{i}& \text{sur}\Sigma \\ \frac{\partial p}{\partial n}={\rho }_{f}{\omega }^{2}{x}_{{f}_{i}}{n}_{i}& \text{sur}\Sigma \text{et}\text{sur}\text{SL}\\ p={\rho }_{f}gz& \text{sur}\text{SL}\end{array}\)
To express the functional, we use the law of behavior on the free surface. By considering an admissible field of displacement \(\mathit{dz}\) we obtain [bib2]:
\(\underset{\text{SL}}{\int }\rho gz\delta z\text{ds}=\underset{\text{SL}}{\int }p\delta z\text{ds}\)
Finally, let’s say the functional of the total fluid structure system subjected to gravity:
\(\begin{array}{}L({x}_{s},p,z)=\frac{1}{2}\underset{{V}_{s}}{\int }\left[{\sigma }_{\text{ij}}({x}_{s}){\varepsilon }_{\text{ij}}({x}_{s})-{\rho }_{s}{\omega }^{2}{x}_{s}^{2}\right]-\underset{\Sigma }{\int }p{x}_{s}\text{.}nd\Sigma +\frac{1}{2}{\rho }_{f}\underset{{V}_{f}}{\int }\left[\frac{1}{{\omega }^{2}}(\text{grad}p{)}^{2}-\frac{{p}^{2}}{{c}^{2}}\right]\text{dV}\\ \phantom{L({x}_{s},p,z)}+\frac{1}{2}\underset{\text{SL}}{\int }\rho g{z}^{2}\text{dS}-\underset{\text{SL}}{\int }pz\text{dS}\end{array}\)
This consideration of gravity implies two additional terms in the functional description of the fluid:
a potential energy term related to the free surface: \(\frac{1}{2}\underset{\text{SL}}{\int }\rho g{z}^{2}\text{ds}\)
a term due to the work of hydrodynamic pressure in the movement of the free surface: \(\underset{\text{SL}}{\int }pz\text{ds}\)
However, it should be noted that this is not the only effect of gravity since at every point on the wall \(\Sigma\) a permanent pressure \(\mathrm{-}\rho gz\) is exerted (where \(z\) is the altitude of the point \(M\) in question: it is assumed that \(z=0\) at the level of the free surface at equilibrium). Point \(M\) is animated by an infinitesimal \({X}_{s}\) movement, so the surface element \(d\Sigma\) varies and so does the effort due to permanent pressure. This effort is responsible for an additional stiffness term in addition to the structural rigidity in the system. It could cause buckling of the structure by cancelling out structural rigidity. This effect is negligible on vibratory characteristics ([bib2], [bib1]), so it is not taken into account.
3.2.2. Finite element discretization#
To obtain the discretized form of the functional, we replace each integral with a sum of integrals on each element \(i\) of the discretized system, then we use a finite element approximation of the unknown displacement and pressure functions on each element \(i\) [bib18].
The unknowns are \({\mathrm{X}}_{\mathrm{s}}(u,v,w),p,z\), so by setting \(\text{Ni}\) we have the shape functions (or nodal interpolation functions on the element \(i\)):
\(\{\begin{array}{c}\sigma =D\varepsilon \\ \varepsilon =Bu\end{array}\) \(\{\begin{array}{c}{x}_{s}={N}_{i}u\\ p(x,y,z)={N}_{i}p\\ \nabla p=\overline{{N}_{i}}p\end{array}\) \(\{\begin{array}{c}{x}_{s}\text{.}n={N}_{\Sigma i}u\\ z={N}_{\mathrm{Si}}z\end{array}\)
where \(\delta ,p\), are the unknowns at structural nodes and fluid nodes, and \(z\) are the unknowns at the free surface.
Hence the discretized expression of the functional associated with the problem:
\(L={u}^{t}(K-{\omega }^{2}M)u+{p}^{t}(\frac{H}{{\rho }_{f}{\omega }^{2}}-\frac{Q}{{\rho }_{f}{c}^{2}})p+{z}^{t}{\rho }_{f}g{K}_{z}z-2{p}^{t}{M}_{z}z-2{p}^{t}Cu\)
with
\(K=\sum _{i}\underset{\text{Vi}}{\int }{N}_{i}^{t}{B}_{i}^{t}{D}_{i}{B}_{i}{N}_{i}{\text{dV}}_{i}\) |
structural stiffness matrix |
\(M=\sum _{i}\underset{\text{Si}}{\int }{N}_{i}^{t}{\rho }_{f}{N}_{i}{\text{dV}}_{i}\) |
mass matrix of the structure |
and
\(Q=\sum _{i}\underset{{\text{Vi}}_{\text{fl}}}{\int }{N}_{i}^{t}{N}_{i}{\text{dV}}_{i}\) |
|
\({K}_{z}=\sum _{i}\underset{\text{Si}}{\int }{N}_{\text{Si}}^{t}{N}_{\text{Si}}{\text{dS}}_{i}\) |
|
\(H=\sum _{i}\underset{{\text{Vi}}_{\text{fl}}}{\int }\overline{{N}_{i}^{t}}\overline{{N}_{i}}{\text{dV}}_{i}\) |
where \(c\) is the speed of sound in the fluid, \({\mathrm{\rho }}_{f}\) the density of the fluid, and where \({K}_{f}\) corresponds to the potential energy of the fluid, \({K}_{z}\) to the potential energy of the free surface, \(H\) to the kinetic energy of the fluid, \(M\) to the fluid-solid coupling, and \({M}_{z}\) to the \(p\mathrm{-}z\) coupling on the free surface.
The finite element approximation of the complete problem then leads to the following matrix system:
\(\left[\begin{array}{ccc}K& -C& 0\\ 0& H& 0\\ 0& -{M}_{z}& {K}_{z}\end{array}\right]\left\{\begin{array}{c}u\\ p\\ z\end{array}\right\}-{\omega }^{2}\left[\begin{array}{ccc}M& 0& 0\\ {\rho }_{f}C& \frac{Q}{{c}^{2}}& {\rho }_{f}{M}_{z}\\ 0& 0& 0\end{array}\right]\left\{\begin{array}{c}u\\ p\\ z\end{array}\right\}=0\)
The first equation corresponds to the movement of the structure under pressure forces, the second to that of the movement of the fluid coupled to the structure and to the free surface, the third is the free surface equation.
However, the problem written in this way has non-symmetric mass and stiffness matrices, which prevents the use of conventional Code_Aster resolution algorithms.
3.2.3. Introduction of an additional variable#
To make the problem symmetric and to be able to use classical resolution methods, an additional variable is introduced: the potential of movements in the fluid \(j\) [bib2].
\({\mathrm{X}}_{\mathrm{f}}\mathrm{=}\text{grad}\varphi\) i.e. \(\rho {\omega }^{2}\varphi \mathrm{=}p\)
This additional unknown is linked to the unknowns of the problem, which leads to a singular stiffness matrix.
The complicated structure-fluid problem is reformulated:
\(\{\begin{array}{cc}{C}_{\text{ijkl}}{x}_{{s}_{k,\text{lj}}}+{\omega }^{2}{\rho }_{s}{x}_{{s}_{i}}=0& \text{dans}{V}_{s}\\ \Delta \varphi +\frac{{\omega }^{2}}{{\rho }_{f}{c}^{2}}p=0& \text{dans}{V}_{f}\\ p={\rho }_{f}{\omega }^{2}\varphi & \text{dans}{V}_{f}\\ {\sigma }_{\text{ij}}{n}_{i}={C}_{\text{ijkl}}{x}_{{s}_{k,l}}{n}_{i}\text{=}\text{-}{\rho }_{f}{\omega }^{2}\varphi {\delta }_{\text{ij}}{n}_{i}& \text{sur}S\\ \frac{\partial \varphi }{\partial n}={x}_{{f}_{i}}{n}_{i}& \text{sur}S\\ p={\rho }_{f}gz& \text{sur}\text{SL}\end{array}\)
This leads to the functionality of the coupled system:
\(\begin{array}{}L({x}_{s},p,j,z)=\frac{1}{2}\underset{{V}_{s}}{\int }\left[{\sigma }_{\text{ij}}({x}_{s}){\varepsilon }_{\text{ij}}({x}_{s})-{\rho }_{s}{\omega }^{2}{x}_{s}^{2}\right]+\frac{1}{2}\underset{{V}_{f}}{\int }\frac{{p}^{2}}{{c}^{2}}\text{dV}+\frac{1}{2}\underset{\text{SL}}{\int }\rho g{z}^{2}\text{ds}\\ \phantom{L({x}_{s},p,j,z)}-{\omega }^{2}\left[\underset{\Sigma }{\int }{\rho }_{f}\varphi {x}_{s}nd\Sigma +\underset{\text{SL}}{\int }{\rho }_{f}\varphi z\text{ds}+\underset{{V}_{f}}{\int }\left[\frac{{\rho }_{f}}{2}(\text{grad}\varphi {)}^{2}+\frac{p\varphi }{{c}^{2}}\right]\text{dV}\right]\end{array}\)
Or by discretizing:
\(\begin{array}{}L={\delta }^{t}(K-{\omega }^{2}M)\delta +\frac{1}{{\rho }_{f}{c}^{2}}{p}^{t}Qp+{\rho }_{f}g{z}^{t}{K}_{z}z\\ \phantom{\text{L=}}-2{\omega }^{2}\left[\frac{{\rho }_{f}}{2}{\varphi }^{t}H\varphi +{\rho }_{f}{\varphi }^{t}C\delta +{\rho }_{f}{\varphi }^{t}{M}_{z}z+\frac{1}{{c}^{2}}{p}^{t}Q\varphi \right]\end{array}\)
What is written, in matrix form:
\(\left[\begin{array}{cccc}K& 0& 0& 0\\ 0& \frac{Q}{{\rho }_{f}{c}^{2}}& 0& 0\\ 0& 0& 0& 0\\ 0& 0& 0& {\rho }_{f}g{K}_{z}\end{array}\right]\left[\begin{array}{c}\delta \\ p\\ \varphi \\ z\end{array}\right]-{\omega }^{2}\left[\begin{array}{cccc}M& 0& {\rho }_{f}C& 0\\ 0& 0& \frac{Q}{{c}^{2}}& 0\\ {\rho }_{f}{C}^{t}& \frac{{Q}^{t}}{{c}^{2}}& {\rho }_{f}H& {\rho }_{f}{M}_{z}\\ 0& 0& {\rho }_{f}{M}_{z}^{t}& 0\end{array}\right]\left[\begin{array}{c}\delta \\ p\\ \varphi \\ z\end{array}\right]=0\)