#  About penalty-duality methods in fluid-structure interactions

Philippe Destuynder∗ and Erwan Liberge∗

###### Abstract

In fluid-structure interaction problems, some people use a penalty method for
positioning the structure inside the fluid. This is usually performed by
considering that the fluid is very stiff or/and very heavy at the place
occupied by the structure. These methods are very convenient for the
programming point of view but lead to ill conditioned operators. This is a
drawback in the numerical solution methods. In particular the forces applied
to the structure by the surrounding flow are not accurately estimated because
the penalty parameter -which is a very large number- appears in their
expressions. We suggest in this paper a mathematical analysis of the
difficulties encountered and we discuss how the penalty-duality method of D.
Bertsekas can be an interesting alternative to overcome them.

* Laboratoire LaSIE

Université de La Rochelle

Avenue Michel Crépeau, 17042 La Rochelle France

<EMAIL_ADDRESS><EMAIL_ADDRESS>
##  1 Introduction

Let us consider -for sake of brevity- a two dimensional bounded and connected
open set denoted by $\Omega$ with boundary $\Gamma$. The unit normal to
$\Gamma$ outward $\Omega$ is denoted by $\nu$. Inside $\Omega$ there is an
open subset $S$ with boundary $\partial S$. The complementary open set of
$\overline{S}$ in $\Omega$ is $\Omega_{S}=\Omega\setminus{\overline{S}}$. The
set $\overline{S}$ is occupied by a structure that we assume to be rigid
or/and heavy for sake of simplicity. But there is no difficulty for extending
our discussion to a flexible structure. The open set $\Omega_{S}$ is occupied
by a viscous fluid that we assume to be incompressible in order to close the
set of equations modeling the flow. The open set $S$ is moving with the
displacement of the structure. Hence a new discretization of $\Omega_{S}$
would be necessary in case where one chooses to model separately the two
media. Furthermore the transport of the physical field requires some special
transformation as a local (or global) Euler-Lagrange parametrization. In
addition the approximation and the transport of the forces (or stresses)
interacting between the two media is a tough problem in terms of precision.

$\nu$$\Omega_{s}$$S$$\Gamma=\Gamma_{0}\cup\Gamma_{1}$$\partial S$$\nu$ Figure
1: Schematic description of fluid structure interaction problem

An alternative which has been introduced by several authors [1] [12], consists
in considering that the two media obey to the same equilibrium laws but differ
by different constitutive relationships. More precisely, the coefficients of
the constitutive equation are very large in order to prescribe a rigid body
motion in the open set occupied by the structure. Another possibility consists
in forcing the velocity field of the fluid to be equal to a rigid body one
(the one of the structure) by adding a penalty term using the $L^{2}(S)$-norm
[4] [27] [5]. The advantage of these approaches is to work on a single mesh
(but requires some interpolations as far as the structure is moving on this
mesh and its boundary has no reason to respect the vertices of the mesh. The
second one is a very promising method for particles methods and lattice
Boltzmann modeling [4] [27] [5]. Vortex methods are also well adapted to this
strategy [26] [13] [14].

Several authors [31] [32] [24] have suggested to take into account the penalty
term only on the boundary of the structure. But the question of the ill
conditioning and the one of the mapping between the global mesh and this
boundary should be handled in order to improve the strategy, mainly in a
multiprocessor programming. Our purpose in this paper is only to discuss
mathematically the penalty method with these different possibilities
concerning the penalty term. We choose to remain at the theoretical level in
order to point out the mathematical aspects of the problems. Furthermore, our
analysis is performed for a given position of the structure and with a given
rigid body velocity of it. In a dynamical model one could consider that it is
the case when the classical algorithm (transport prediction-diffusion
correction) is used (see P.A. Raviart and V. Girault [34]. Concerning the
movement of the structure it is assumed to be sufficiently slow in order to
adopt an explicit scheme for its displacements (small reduced frequencies
approximation [22]-[16]).

Let us first state few notations used in the following. The velocity of the
fluid particles is denoted by ${\mathbb{u}}=\\{u_{i}\\},\;i=1,2$ (the
extension to 3D is not a difficulty for the theoretical approach but a real
one for the implementation). The strain ratio is
$\gamma({\mathbb{u}})=\\{\gamma_{ij}({\mathbb{u}})=\frac{1}{2}(\partial_{i}u_{j}+\partial_{j}u_{i})\\}$
and the stress ratio is ($\hbox{div}({\mathbb{u}})=0$):

$\sigma_{ij}=2\mu\gamma_{ij}({\mathbb{u}}),$

$\mu$ being the viscosity of the fluid. The unit normal to $\partial S$ inward
$S$ is also denoted by $\nu$ and the movement of the solid is represented by
its rigid body velocity (${\mathbb{k}}$ is the unit vector normal to the plan
containing $\Omega$; $r$ is the velocity rotation around ${\mathbb{k}}$ and
${\mathbb{a}}$ the velocity of point $O$). Hence the the velocity of the
structure at point ${\mathbb{x}}$ is:

${\mathbb{d}}=\\{d_{i}\\}={\mathbb{a}}+r{\mathbb{k}}\wedge\mathbb{Ox}=\left\\{\begin{array}[]{l}a_{1}-rx_{2}\\\
\\\ a_{2}+rx_{1}\end{array}\right\\},\;\;{\mathbb{x}}=(x_{1},x_{2})\;\hbox{
and }{\mathbb{a}}=(a_{1},a_{2}).$ (1)

The kinematical continuity between the fluid and the structure is traduced by
the following relation:

${\mathbb{u}}={\mathbb{a}}+r{\mathbb{k}}\wedge\mathbb{Ox}\hbox{ on }\partial
S.$ (2)

The flow model inside $\Omega_{S}$ consists in finding ${\mathbb{u}}$ such
that ($\Gamma_{0}$ and $\Gamma_{1}$ being the two distinct components of
$\Gamma$ where respectively the flow velocity -respectively the normal stress-
is given and $p$ is the static pressure):

$\left\\{\begin{array}[]{l}\sigma({\mathbb{u}})=\mu\gamma({\mathbb{u}}),\;\mathbb{div}(\sigma)-\mathbb{grad}(p)=0\;\hbox{
in }\;\Omega_{S},\\\ \\\ {\mathbb{div}}({\mathbb{u}})=0\;\hbox{ in
}\;\Omega_{S},\\\ \\\ {\mathbb{u}}={\mathbb{u}}_{0}\hbox{ on
}\;\Gamma_{0}\subset\Gamma,\;\sigma.\nu-p\nu={\mathbb{g}}\hbox{ on
}\Gamma_{1}\subset\Gamma,\\\ \\\ \hbox{ and finally the condition stated in
equation (\ref{eq2}) on $\partial S$.}\end{array}\right.$ (3)

This a quasi-static flow model (following the terminology of aeroelasticity
and hydroelasticity [22]-[16]) because the inertia forces aren’t taken into
account in the flow.

The equations of the model for the structure are the following ones (let us
notice that: $J_{O}=J_{G}+M||\mathbb{OG}||^{2}$ where $J_{G}$ is the inertia
around the center of mass $G$ and $J_{O}$ the one around the point $O$ and
finally the dot stands for the scalar product or more generally for the matrix
product):

$\left\\{\begin{array}[]{l}M\dot{\mathbb{a}}+M\dot{r}{\mathbb{k}}\wedge{\mathbb{OG}}+Mr{\mathbb{k}}\wedge{\mathbb{a}}-Mr^{2}\mathbb{OG}=\displaystyle\int_{\partial
S}p\nu-\sigma.\nu,\\\ \\\
J_{O}\dot{r}+M(\dot{\mathbb{a}},{\mathbb{k}},\mathbb{OG})+Mr({\mathbb{a}}.\mathbb{OG})+2Mr({\mathbb{k}}\wedge\mathbb{OG})(\mathbb{k}\wedge\mathbb{a})\\\
\\\ \hskip 113.81102pt=\displaystyle\int_{\partial
S}(p\nu-\sigma.\nu)\wedge\mathbb{Ox},\\\ \\\ +\hbox{initial conditions on
${\mathbb{a}}$ and $r$.}\end{array}\right.$ (4)

######  Remark 1.

In fact, we have considered in the previous modeling, that the movement of the
structure is quite slow compared to the velocity necessary for establishing
the permanent flow. In other words, it is the quasic-static approximation
corresponding to low reduce frequencies [22]-[16]. Therefore the problem to be
solved here is a simple first order differential equation in ${\mathbb{a}}$
and $r$ coupled with a static linear partial differential equation modeling
the flow around the moving structure. Therefore we focus on the solution
method for this last problem explicited at equation (3) assuming that
$\partial S$, ${\mathbb{u}}_{0}$ and ${\mathbb{g}}$ are given. $\Box$

##  2 The initial formulation for the flow

First of all let us introduce the problem (3) as an optimization one. We
introduce the functional space for the velocity field of the fluid:

$V_{0}=\\{{\mathbb{v}}=\\{v_{i}\\},\;i=1,2,\;v_{i}\in
H^{1}(\Omega_{S}),\;\mathbb{div}({\mathbb{v}})=0\;\hbox{ in
}\Omega_{S},\;{\mathbb{v}}=0\;\hbox{on}\;\Gamma_{0}\\}.$ (5)

Then we set ($\sigma.\gamma=\sum_{i,j\in\\{1,2\\}}\sigma_{ij}\gamma_{ij}$):

$\forall{\mathbb{v}}\in
V_{0},\;J({\mathbb{v}})=\displaystyle\frac{1}{2}\displaystyle\int_{\Omega_{S}}\sigma({\mathbb{v}}).\gamma({\mathbb{v}})-\displaystyle\int_{\Gamma_{1}}{\mathbb{g}}.{\mathbb{v}}$
(6)

It is classical that equations (3) are equivalent to:

$\min_{{\mathbb{v}}\in
V_{0}\;{\mathbb{v}}={\mathbb{a}}+r{\mathbb{k}}\wedge{\mathbb{O}x}\;\hbox{ on
}\partial S}J({\mathbb{v}}).$ (7)

Existence and uniqueness of the solution -say ${\mathbb{u}}$\- are also well
known (see for instance [34]).

######  Remark 2.

The most popular method for solving (3) is based on the mixed formulation
where both the velocity and the pressure are unknowns. The results are very
good and an extension with the acceleration term:

$\displaystyle\varrho\left[\displaystyle\frac{\partial{\mathbb{u}}}{\partial
t}+{\mathbb{u}}\nabla{\mathbb{u}}\right],$

works perfectly. But, this formulation is restricted to incompressible (or
slightly compressible) flows. The details are given for instance in [34].
$\Box$

##  3 The extension to the full domain using penalty

###  3.1 The definitions of the penalty models

The velocity ${\mathbb{u}}_{0}$ on the structure is defined by a rigid body
description. Such velocity fields satisfy the relation:

$\gamma_{ij}({\mathbb{u}}_{0})=0\;\;\forall i,j\in\\{1,2\\}.$ (8)

The idea of the penalty method in this case, consists in replacing problem (7)
by the next one ($\alpha\geq 0,\;\beta\geq 0,\;\gamma\geq 0\;\varepsilon>0$):

$\hskip-8.5359pt\left\\{\begin{array}[]{l}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\displaystyle\min_{{\mathbb{v}}\in
W_{0}}J^{\varepsilon}({\mathbb{v}})\\\ \\\ \hbox{where:}\\\ \\\
\hskip-5.69046ptJ^{\varepsilon}({\mathbb{v}})\hskip-1.42271pt=\displaystyle\hskip-1.42271pt\displaystyle\frac{1}{2}\left[\hskip-1.42271pt\displaystyle\int_{\Omega}\hskip-2.84544pt\sigma({\mathbb{v}}).\gamma({\mathbb{v}})\hskip-1.42271pt+\hskip-1.42271pt\displaystyle\frac{\alpha}{\varepsilon}\hskip-1.42271pt\displaystyle\int_{\partial
S}\hskip-4.26773pt||{\mathbb{v}}\hskip-1.42271pt-\hskip-1.42271pt{\mathbb{u}}_{0}||^{2}\hskip-1.42271pt+\hskip-1.42271pt\displaystyle\frac{\beta}{\varepsilon}\hskip-1.42271pt\displaystyle\int_{S}\hskip-1.42271pt\sigma({\mathbb{v}}){.\gamma({\mathbb{v}})}\hskip-1.42271pt+\hskip-1.42271pt\displaystyle\frac{\gamma}{\varepsilon}\hskip-1.42271pt\displaystyle\int_{S}\hskip-1.42271pt||{\mathbb{v}}\hskip-1.42271pt-\hskip-1.42271pt{\mathbb{u}}_{0}||^{2}\right]\\\
\\\ \hbox{and:}\\\ \\\
W_{0}=\left\\{{\mathbb{v}}=\\{v_{i}\\},\;i=1,2,\;v_{i}\in
H^{1}(\Omega),\;\mathbb{div}({\mathbb{v}})=0\hbox{ on
}\Omega,\;{\mathbb{v}}=0\hbox{ on }\Gamma_{0}\right\\}.\end{array}\right.$ (9)

Here again the existence and uniqueness of a solution are standard. The
solution of (9) is a function of $\alpha,\;\beta,\;\gamma$ and mainly of
$\varepsilon$. Our purpose in this section, is to analyze the behaviour of
this solution denoted by $\mathbb{u}^{\varepsilon}$ for several choices of
$\alpha,\beta$ and $\gamma$, when $\varepsilon$ tends to zero. Let us make
explicit the variational formulation of (9). The solution
${\mathbb{u}}^{\varepsilon}$ is characterized by:

$\left\\{\begin{array}[]{l}{\mathbb{u}}^{\varepsilon}\in W_{0}\;\hbox{such
that }\forall{\mathbb{v}}\in W_{0}:\\\ \\\
\displaystyle\int_{\Omega}\sigma({\mathbb{u}}^{\varepsilon}).\gamma({\mathbb{v}})+\displaystyle\frac{\alpha}{\varepsilon}\displaystyle\int_{\partial
S}({\mathbb{u}}^{\varepsilon}-{\mathbb{u}}_{0}).\mathbb{v}+\displaystyle\frac{\beta}{\varepsilon}\displaystyle\int_{S}\sigma({\mathbb{u}}^{\varepsilon}).\gamma({\mathbb{v}})\\\
\\\ \hskip
142.26378pt+\displaystyle\frac{\gamma}{\varepsilon}\displaystyle\int_{S}({\mathbb{u}}^{\varepsilon}-{\mathbb{u}}_{0}).\mathbb{v}=\displaystyle\int_{\Gamma_{1}}{\mathbb{g}}.{\mathbb{v}}\end{array}\right.$
(10)

######  Remark 3.

In the papers by Benamour-Liberge-Béghein [4] [27] [5], the authors set
$\alpha=0$ and also $\beta=0$ but $\gamma\neq 0$. This is a very convenient
strategy because one can use the mass matrix which is required for the dynamic
analysis (the mass density for this new term is artificially equal to
$\gamma/\varepsilon$ in $S$). The solution ${\mathbb{u}}^{\varepsilon}$ and
also the normal stress associated should be continuous across the boundary
$\partial S$ and this last condition (on the normal stresses) is not well
satisfied localy for $\varepsilon\simeq 0$. The reason is that the bilinear
form involved in the limit model for $\alpha=\beta=0$, which couples
${\mathbb{u}}^{0}$ on $\Omega_{S}$ and ${\mathbb{u}}^{1}$ on $S$, is only
coercive on the space:

$W_{\gamma}=\\{{\mathbb{v}}\hbox{ s.t.}\;\hbox{div}({\mathbb{v}})=0\hbox{ in
}\Omega\hbox{ and
}{\mathbb{v}}_{|\Omega}\in{\Huge[}H^{1}(\Omega_{S}){\Huge]}^{2},\;{\mathbb{v}}_{|S}\in{\Huge[}L^{2}(S){\Huge]}^{2}\\}.$
(11)

It is a larger space than $W_{0}$ and such that only the continuity of the
normal velocity ${\mathbb{v}}.\nu$ can be verified in the space
$H^{-1/2}(\partial S)$ using the trace Theorem (see for instance [29]) which
is a weak result, mainly concerning the evaluation of the stress due to the
fluid and applied to the structure. Furthermore, the tangential component of
${\mathbb{v}}_{|S}$ on $\partial S$ doesn’t exist in such a space.
Nevertheless this strategy is very tricky for the numerical implementation and
can lead to nice and cheap evaluation of the velocity field. Therefore we
discuss it in the following in order to check which kind of information can be
lost concerning the tangential component of the normal stress along $\partial
S$.

If $\beta\neq 0$ whatever would be $\alpha$ and $\gamma$, one can can ensure
the continuity of the two components of the velocity field
${\mathbb{u}}^{\varepsilon}$ solution of (9) when $\varepsilon\rightarrow 0$.
Hence only the rigid body motion of the open set $S$ is taken into account.
Let us point out that if $\alpha=\gamma=0$ nothing guarantees that it will be
the one prescribed by the equations of the mechanics. In fact
${\mathbb{u}}_{0}$ has disappeared from the penalty model and this statement
is propped by the numerical tests of section 4.

A discussion is carried out at subsection 3.2 in order to explain the
drawbacks of this nice strategy which is nevertheless used for convenience in
the implementation in existing softwares by many of authors. $\Box$

######  Remark 4.

The implementation of the penalty term on $S$ for $\beta>0$ is easy to do as
far it only consists in choosing a new expression for the viscosity
coefficient of the fluid. We set $\mu$ in $\Omega_{S}$ and
$\mu(1+\frac{\beta}{\varepsilon})$ in $S$. $\Box$

######  Remark 5.

In a numerical implementations, one needs to use two mappings because the
movement of the structure does respect the sides of the mesh used in the
approximation of the fluid movement. One -say $R^{h}$\- maps the position of
$\partial S$ on the mesh and the other -say $P^{h}$\- maps the mesh of
$\Omega$ on $\partial S$. But this embedding problem is fully and accurately
solved in a practical formulation by C. Farhat and his coworkers [10]-[25]. A
mathematical formulation of this arbitrary Euler-Lagrange parametrization is
briefly discussed in [15] using domain derivative tools for improving the
normal stresses approximation between the fluid and the structure. $\Box$

###  3.2 The asymptotic analysis

Let set a priori:

${\mathbb{u}}^{\varepsilon}={\mathbb{u}}^{0}+\varepsilon{\mathbb{u}}^{1}+\ldots$
(12)

By introducing this expression into the variational formulation (10) and by
equating the terms of same power in $\varepsilon$ one obtains the following
necessary conditions:

$\left\\{\begin{array}[]{l}\hbox{ i) terms of order
}\varepsilon^{-1}:\;\forall{\mathbb{v}}\in W_{0},\;(\hbox{ let us notice that
$\sigma({\mathbb{u}}_{0})=0$})\\\ \\\ \alpha\displaystyle\int_{\partial
S}({\mathbb{u}}^{0}-{\mathbb{u}}_{0}).{\mathbb{v}}+\beta\displaystyle\int_{S}\sigma({\mathbb{u}^{0}}).\gamma({\mathbb{v}})+\gamma\displaystyle\int_{S}({\mathbb{u}}^{0}-{\mathbb{u}}_{0}).{\mathbb{v}}=0,\\\
\\\ \hbox{ ii) terms of order zero: }\forall{\mathbb{v}}\in W_{0},\\\ \\\
\displaystyle\int_{\Omega}\sigma({\mathbb{u}^{0}}).\gamma({\mathbb{v}})+\alpha\displaystyle\int_{\partial
S}{\mathbb{u}}^{1}.{\mathbb{v}}+\beta\displaystyle\int_{S}\sigma({\mathbb{u}^{1}}).\gamma({\mathbb{v}})+\gamma\displaystyle\int_{S}{\mathbb{u}}^{1}.{\mathbb{v}}={\displaystyle\int_{\Gamma_{1}}{\mathbb{g}}.{\mathbb{v}}},\\\
\\\ \\\ \\\ \hbox{ iii) terms of order }\varepsilon:\;\forall{\mathbb{v}}\in
W_{0},\\\ \\\
\displaystyle\int_{\Omega}\sigma({\mathbb{u}^{1}}).\gamma({\mathbb{v}})+\alpha\displaystyle\int_{\partial
S}{\mathbb{u}}^{2}.{\mathbb{v}}+\beta\displaystyle\int_{S}\sigma({\mathbb{u}^{2}}).\gamma({\mathbb{v}})+\gamma\displaystyle\int_{S}{\mathbb{u}}^{2}.{\mathbb{v}}=0,\\\
\\\ \ldots\end{array}\right.$ (13)

Setting on $S$:

${\mathbb{u}}^{0}={\mathbb{u}}_{0}+{\mathbb{u}}^{00},\;{\mathbb{u}}^{00}\in
W_{0},$

with the condition ${\mathbb{u}}^{00}=0$ on $\partial S$, equation i) of (13)
implies that:

$\left\\{\begin{array}[]{l}\hbox{ if }\alpha,\beta>0\hbox{ but }\gamma\geq
0:\\\ \\\ \hskip 28.45274pt=>{\mathbb{u}}^{0}={\mathbb{u}}_{0}\hbox{ on
}\partial S\hbox{ and }\gamma_{ij}({\mathbb{u}}^{00})=0\hbox{ in
}S=>{\mathbb{u}}^{00}=0\hbox{ on }S,\\\ \\\ \hbox{ if }\beta=0\hbox{ but
}\gamma>0\hbox{ one has for any $\alpha\geq 0$}:\\\ \\\ \hskip
28.45274pt=>{\mathbb{u}}^{00}=0.\end{array}\right.$ (14)

Therefore ${\mathbb{u}}^{0}$ is a rigid body velocity on $S$ and equal to
${\mathbb{u}}_{0}$ as far as $\alpha$ or $\gamma$ is (are) strictly positive.
But it should be noticed that this last statement is valid only because
$\alpha>0$ or $\gamma>0$. For $\alpha=\gamma=0$ one could only claim that
${\mathbb{u}}^{0}$ is a rigid body velocity on $S$ but not necessarily equal
to ${\mathbb{u}}_{0}$. At this step, the case $\alpha=\beta=0$ seems to be the
most convenient, but the following contradicts this first conclusion.

We now introduce a new functional space by:

$K_{0}=\\{{\mathbb{v}}\in W_{0},\;{\mathbb{v}}=0\hbox{ in }S\\}.$ (15)

From equation ii) of (13), one deduces that:

$\forall{\mathbb{v}}\in
K_{0},\;\displaystyle\int_{\Omega_{S}}\sigma.({\mathbb{u}}^{0}).\gamma({\mathbb{v}})=\displaystyle\int_{\Gamma_{1}}{\mathbb{g}}.{\mathbb{v}},\;\hbox{
and if $\alpha^{2}+\gamma^{2}>0$ }=>{\mathbb{u}}^{0}={\mathbb{u}}_{0}\hbox{ on
}\partial S.$ (16)

The element ${\mathbb{u}}^{0}\in K_{0}$ is now perfectly defined and is
identical on $\Omega_{S}$ to the initial solution ${\mathbb{u}}$ characterized
at (3).

######  Remark 6.

If $\alpha=\gamma=0$ and $\beta>0$ one has only:

$\gamma_{ij}({\mathbb{u}}^{0})=0\hbox{ in
}S=>{\mathbb{u}}^{0}=a+r{\mathbb{k}}\wedge\mathbb{Ox}\hbox{ in }S.$ (17)

The term ${\mathbb{u}}^{0}$ belongs to the following closed subspace of
$W_{0}$ defined by:

$RBM=\\{{\mathbb{v}}\in W_{0},\;\gamma({\mathbb{v}})=0\hbox{ in }S\\}.$ (18)

Hence in this case ($\alpha=\gamma=0$) the term ${\mathbb{u}}^{0}\in RBM$ is
the unique solution of:

$\forall{\mathbb{v}}\in
RBM,\;\displaystyle\int_{\Omega}\sigma({\mathbb{u}}^{0}).\gamma({\mathbb{v}})=\displaystyle\int_{\Gamma_{1}}{\mathbb{g}}.{\mathbb{v}}\;\;\;\;.$
(19)

But the rigid body motion on $S$ is derived from the solution
${\mathbb{u}}^{0}$ and is no more the one that has been prescribed from the
movement of the structure. This is the fundamental difference that suggests to
avoid the choice $\alpha=0$ or $\gamma=0$. In fact the choice $\beta=0$ is
another possibility which will appear as a smart possibility in the discussion
contained in the next remark 8, even if it has its own drawbacks concerning
the regularity if the boundary $\partial S$ is not smooth enough as we mention
in the following. $\Box$

######  Remark 7.

Let us now consider the case where $\alpha>0$ and $\beta=\gamma=0$. This is
very interesting situation from the practical point of view, because we only
have to manage the position of the boundary $\partial S$ inside the mesh which
can be unchanged during the displacement of the structure (full eulerian
representation of the fluid). Let us also underline that this method is just
for the quasi-static case. If inertia terms are introduced in the fluid then
it doesn’t work in this case. But even for quasi-static cases (no inertia
terms), this assertion could be false for a multiprocessor programming where
each processor could be idling when waiting for the informations concerning
this boundary $\partial S$. At the order minus one in $\varepsilon$ in the
asymptotic expansion, one only gets the condition
${\mathbb{u}}^{0}={\mathbb{u}}_{0}$ on $\partial S$. But at the order zero in
$\varepsilon$ one has:

$\forall{\mathbb{v}}\in
W_{0},\;\displaystyle\int_{\Omega}\sigma({\mathbb{u}}^{0}).\gamma({\mathbb{v}})+\alpha\displaystyle\int_{\partial
S}{\mathbb{u}}^{1}.{\mathbb{v}}=\displaystyle\int_{\Gamma_{1}}{\mathbb{g}}.{\mathbb{v}}.$
(20)

Let us define three closed subspaces of $W_{0}$ by:

$\left\\{\begin{array}[]{l}W_{00}=\\{v\in W_{0},\;{\mathbb{v}}=0\hbox{ on
}\partial S\\},\\\ \\\ V_{00\Omega_{S}}=\\{{\mathbb{v}}\in V_{0}\hbox{ and
}{\mathbb{v}}=0\hbox{ on }\partial S,\;({\mathbb{v}}=0\hbox{ on }S)\\};\\\ \\\
V_{00S}={\Huge[}H^{1}_{0}(S){\Huge]}^{2}.\end{array}\right.$ (21)

It is worth to notice that one can identify $W_{00}$ with the space
$V_{00\Omega}\times V_{00S}$.

By restricting the virtual velocity ${\mathbb{v}}$ to the space $W_{00}$
(${\mathbb{u}}^{0}={\mathbb{u}}_{0}$ (the rigid body velocity) on $\partial
S$) one obtains two distinct sub-problems; one is set on $\Omega_{S}$ and the
other on $S$. They are given hereafter:

$\left\\{\begin{array}[]{l}{\mathbb{u}}^{0}_{|\Omega_{S}}\in
V_{0},\forall{\mathbb{v}}\in
V_{00\Omega_{S}}:\;\displaystyle\int_{\Omega_{S}}\sigma({\mathbb{u}}^{0}).\gamma({\mathbb{v}})=\displaystyle\int_{\Gamma_{1}}{\mathbb{g}}.{\mathbb{v}},\;\;{\mathbb{u}}^{0}={\mathbb{u}}_{0}\hbox{
on }\partial S,\\\ \\\ {\mathbb{u}}^{0}_{|S}\in
V_{00S},\;\forall{\mathbb{v}}\in
V_{00S},\;\displaystyle\int_{S}\sigma({\mathbb{u}}^{0}).\gamma({\mathbb{v}})=0,\;{\mathbb{u}}^{0}={\mathbb{u}}_{0}\hbox{
on }\partial S.\end{array}\right.$ (22)

The system of equations (22) leads to two independent problems for which the
solutions are uniquely defined respectively on $\Omega_{S}$ and on $S$. The
solution ${\mathbb{u}}^{0}$ on $\Omega_{S}$ is precisely the one of the
initial model we started from (3). The solution ${\mathbb{u}}^{0}$ on $S$ can
easily be characterized. Because the solution is unique and because the rigid
body ${\mathbb{u}}_{0}$ is a solution, it is the right one. Hence
${\mathbb{u}}^{0}$ is exactly the solution in the full domain $\Omega$ of the
initial problem.

Nevertheless, it should be underlined that in this case ($\beta=\gamma=0$ and
$\alpha>0$) the penalty term is reduced to a boundary term implying the
$L^{2}(\partial S)$ norm and therefore the continuous model is perfectly
defined. The numerical aspects are easy to handle if we use a separate
modeling of $\partial S$ (we do not need to divide the elements at the
boundary between $S$ and $\Omega_{S}$). But the condition number of the
discretized model remains an important question in the solution method
compared to the case where $\alpha>0$ and $\beta>0$. $\Box$

######  Remark 8.

In this remark we discuss the case where $\alpha=\beta=0$ and $\gamma>0$. The
first equation (13) leads directly to ${\mathbb{u}}^{0}={\mathbb{u}}_{0}$ in
$S$. Therefore, from the second equation (13) ${\mathbb{u}}^{0}$ should
connected to an element ${\mathbb{u}}^{1}$ such that:

$\forall{\mathbb{v}}\in
W_{0},\;\displaystyle\int_{\Omega}\sigma({\mathbb{u}}^{0}).\gamma({\mathbb{v}})+\gamma\displaystyle\int_{S}{\mathbb{u}}^{1}.{\mathbb{v}}=\displaystyle\int_{\Gamma_{1}}{\mathbb{g}}.{\mathbb{v}}.$
(23)

Because we proved that ${\mathbb{u}}^{0}={\mathbb{u}}_{0}$ in $S$, the
characterization of ${\mathbb{u}}^{0}$ on $\Omega_{S}$ is obtained by solving
the equation (which has a unique solution):

${\mathbb{u}}^{0}_{|\Omega_{S}}\in V_{0},\forall{\mathbb{v}}\in
V_{00\Omega_{S}},\;\displaystyle\int_{\Omega_{S}}\sigma({\mathbb{u}}^{0}).\gamma({\mathbb{v}})=\displaystyle\int_{\Gamma_{1}}{\mathbb{g}}.{\mathbb{v}},\;\;{\mathbb{u}}^{0}={\mathbb{u}}_{0}\hbox{
on }\partial S.$ (24)

Consequently, the term ${\mathbb{u}}^{1}$ should satisfy in $S$, the relation:

${\mathbb{u}}^{1}=0,$ (25)

which is not compatible with equation (23) which leads to ($\nu$ is here the
unit normal along $\partial S$ inward $S$):

$\forall{\mathbb{v}}\in W_{0},\;\displaystyle\int_{\partial
S}(\sigma({\mathbb{u}}^{0}).\nu).{\mathbb{v}}+\gamma\displaystyle\int_{S}{\mathbb{u}}^{1}.{\mathbb{v}}=0.$
(26)

In fact it is impossible to compute ${\mathbb{u}}^{1}$ which would appear as a
first order derivative at the origin of ${\mathbb{u}}^{\varepsilon}$ with
respect to $\varepsilon$. This phenomenon is due to a stiff boundary layer
near $\partial S$ which is fully analyzed in the linear case in the book of
J.L. Lions on singular perturbations [30]. Using interpolation technics
between Hilbert spaces, he proved that one could only hope a fractional
derivative order and furthermore in a larger space than the energy one. This
result is extended to Navier-Stokes equations in the paper by P. Angot & all
[1]. Nevertheless the convergence of ${\mathbb{u}}^{\varepsilon}$ to
${\mathbb{u}}^{0}$ is true in energy but for the weak topology for
$\alpha=\beta=0$ and $\gamma>0$ as we proved at Theorem 3.1 at the end of this
remark. But the speed of convergence in $\varepsilon$ is not obvious and
furthermore in a larger space than the energy one. This justifies to use this
penalty method in this framework when the programming advantages are clear
(multiprocessor programming for instance). A numerical discussion concerning
these theoretical results on a very simple model is suggested in section 4.
$\Box$

######  Theorem 3.1.

Let ${\mathbb{u}}^{\varepsilon}\in W_{0}$ be solution of (10) and
${\mathbb{u}}^{0}$ be the solution of (24) in $\Omega_{S}$ such that
${\mathbb{u}}^{0}={\mathbb{u}}_{0}$ in $S$. One has:

$\displaystyle\lim_{\varepsilon\rightarrow
0}||{\mathbb{u}}^{\varepsilon}-{\mathbb{u}}^{0}||_{0,\Omega}=0\hbox{ and
}||{\mathbb{u}}^{\varepsilon}-{\mathbb{u}}_{0}||_{0,S}\leq
c\sqrt{\varepsilon}.$ (27)

$\Box$

###### Proof.

Let set ${\mathbb{v}}={\mathbb{u}}^{\varepsilon}$ in equation (10). Using Korn
inequality [18], one can claim that there exists a constant independent on
$\varepsilon$ -say $c_{0}>0$\- such that:

$c_{0}||{\mathbb{u}}^{\varepsilon}||^{2}_{1,\Omega}+\displaystyle\frac{\gamma}{\varepsilon}||{\mathbb{u}}^{\varepsilon}-{\mathbb{u}}_{0}||^{2}_{0,S}\leq\displaystyle\frac{\gamma}{\varepsilon}\displaystyle\int_{S}{\mathbb{u}}_{0}{\mathbb{u}}^{\varepsilon}+\displaystyle\int_{\Gamma_{1}}{\mathbb{g}}.{\mathbb{u}}^{\varepsilon}.$
(28)

Using the triangular Cauchy-Schwarz111$\forall a,b,\zeta\geq 0,\;2ab\leq\zeta
a^{2}+\frac{1}{\zeta}b^{2}$ inequality and the continuity of the trace
operator from $H^{1}(\Omega)$ into $L^{2}(\Gamma_{1})$, on deduces the
existence of another constant independent on $\varepsilon$ -say $c_{1}$\- such
that:

$c_{0}||{\mathbb{u}}^{\varepsilon}||_{1,\Omega}^{2}+\displaystyle\frac{\gamma}{\varepsilon}||{\mathbb{u}}^{\varepsilon}-{\mathbb{u}}_{0}||_{0,S}^{2}\leq
c_{1}.$ (29)

Therefore one can extract from the sequence ${\mathbb{u}}^{\varepsilon}$ a
subsequence denoted by ${\mathbb{u}}^{\varepsilon^{\prime}}$ and such that:

$\left\\{\begin{array}[]{l}i)\hbox{ in
}S:\;\displaystyle\lim_{\varepsilon\rightarrow
0}||{\mathbb{u}}^{\varepsilon}-{\mathbb{u}}_{0}||_{0,S}=0\hbox{ (strong
convergence with the rate }\sqrt{\varepsilon}),\\\ \\\ i)\;\hbox{ in
}\Omega:\;\displaystyle\lim_{\varepsilon^{\prime}\rightarrow
0}{\mathbb{u}}^{\varepsilon^{\prime}}={\mathbb{u}}^{*}\hbox{ in
${\Huge[}H^{1}(\Omega){\Huge]}^{2}$-weakly; hence
${\mathbb{u}}^{*}={\mathbb{u}}_{0}$ in $S$}.\end{array}\right.$ (30)

Taking the limit in (10), one characterizes ${\mathbb{u}}^{*}$ as the unique
solution of:

${\mathbb{u}}^{*}={\mathbb{u}}_{0}\hbox{ on }\partial S\hbox{ and
}\forall{\mathbb{v}}\in
K_{0},\;\displaystyle\int_{\Omega_{S}}\sigma({\mathbb{u}}^{*}).\gamma({\mathbb{v}})=\displaystyle\int_{\Gamma_{1}}{\mathbb{g}}.{\mathbb{v}}.$
(31)

In other words ${\mathbb{u}}^{*}={\mathbb{u}}^{0}$ in all $\Omega$. Because of
the uniqueness of ${\mathbb{u}}^{0}$, all the sequence
${\mathbb{u}}^{\varepsilon}$ converges to ${\mathbb{u}}^{0}$ in
$H^{1}(\Omega)$-weakly and strongly in $L^{2}(S)$. But the embeding from
$H^{1}(\Omega)$ into $L^{2}(\Omega)$ is compact, hence the convergence of the
sequence ${\mathbb{u}}^{\varepsilon}$ to ${\mathbb{u}}^{0}$ is strong in
${\Huge[}L^{2}(\Omega){\Huge]}^{2}$. ∎

In the following of this section we assume that $\alpha>0$ and we start with
$\beta>0$ in order to reduce singular perturbations problems in the vicinity
of $\partial S$ which would damage the computation of the forces applied to
the structure by the fluid. Let us go on with the relation ii) of equation
(13) but now for arbitrary element ${\mathbb{v}}$ in $W_{0}$. Because
${\mathbb{u}}^{0}$ is perfectly determined, it remains to find an element (not
necessarily unique) ${\mathbb{u}}^{1}$ in $W_{0}$ such that (let us emphasize
on the fact that the right-hand side is a linear and continuous form on
$W_{0}$ without any additional regularity required on ${\mathbb{u}}^{0}$):

$\forall{\mathbb{v}}\in W_{0},\;\alpha\displaystyle\int_{\partial
S}{\mathbb{u}}^{1}.{\mathbb{v}}+\beta\displaystyle\int_{S}\sigma({\mathbb{u}}^{1}).\gamma({\mathbb{v}})=-\displaystyle\int_{\Omega_{S}}\sigma({\mathbb{u}}^{0}).\gamma({\mathbb{v}})+\displaystyle\int_{\Gamma_{1}}{\mathbb{g}}.{\mathbb{v}}\;\;.$
(32)

Or else from an integration by parts (the velocity fields are divergence
free):

$\forall{\mathbb{v}}\in W_{0},\;\alpha\displaystyle\int_{\partial
S}{\mathbb{u}}^{1}.{\mathbb{v}}+\beta\displaystyle\int_{S}\sigma({\mathbb{u}}^{1}).\gamma({\mathbb{v}})=-\displaystyle\int_{\partial
S}(\sigma({\mathbb{u}}^{0}).\nu).{\mathbb{v}}\;\;.$ (33)

Let us check that ${\mathbb{u}}^{1}$ is perfectly determined by the previous
relation on $S$ (but not yet on $\Omega_{S}$). We denote by
${\mathbb{u}}^{1p}$ this particular term defined on $S$. The bilinear form
$a^{1}$ defined by:

${\mathbb{u}},{\mathbb{v}}\in{\Huge[}H^{1}(S){\Huge]}^{2}\rightarrow
a^{1}({\mathbb{u}},{\mathbb{v}})=\alpha\displaystyle\int_{\partial
S}{\mathbb{u}}.{\mathbb{v}}+\beta\displaystyle\int_{S}\sigma({\mathbb{u}}).\gamma({\mathbb{v}}),$
(34)

is ${\Huge[}H^{1}(S){\Huge]}^{2}$-coercive (and symmetrical) and therefore,
from Lax-Milgram Theorem, one can claim that ${\mathbb{u}}^{1}$ is perfectly
defined on $S$ as a solution of (33) (using the expression given at (32) the
right-hand side is clearly a linear and continuous form on the functional
space ${\Huge[}H^{1}(S){\Huge]}^{2}$). It is denoted by ${\mathbb{u}}^{1p}$
and it also satisfies the incompressibility condition on $S$.

In order to complete the definition of ${\mathbb{u}}^{1}$ on the whole open
set $\Omega$, we go to the equations obtained at the order one. Choosing an
arbitrary element ${\mathbb{v}}\in K_{0}$ (see (39)), one obtains (let us
point out that the normal stress due to ${\mathbb{u}}^{1}$ is not continuous
across $\partial S$):

$\left\\{\begin{array}[]{l}\displaystyle\int_{\Omega_{S}}\sigma({\mathbb{u}}^{1}).\gamma({\mathbb{v}})=0,\\\
\\\ {\mathbb{u}}^{1}={\mathbb{u}}^{1p}\hbox{ on }\partial S\hbox{ and
}{\mathbb{u}}^{1}_{|\Omega_{S}}\in V_{0}.\end{array}\right.$ (35)

Finally, the term ${\mathbb{u}}^{1}$ is defined on the whole open set $\Omega$
and because of the continuity of the trace on $\partial S$, it belongs to the
space ${\Huge[}H^{1}(\Omega){\Huge]}^{2}$. Its restriction to $S$ is
${\mathbb{u}}^{1p}$ characterized at (33) and its restriction to $\Omega_{S}$
is solution of (35). Globally ${\mathbb{u}}^{1}$ belongs to the space $W_{0}$
because it is divergence free on each open subset of $\Omega$. It should be
underlined that no additional regularity is required in the computation of
${\mathbb{u}}^{1}$. But this will be different when $\beta=0$.

Let us finish this partial analysis of an assumed asymptotic expansion of
${\mathbb{u}}^{\varepsilon}$, by a partial characterization of
${\mathbb{u}}^{2}$ on $S$ denoted by ${\mathbb{u}}^{2p}$. Once
${\mathbb{u}}^{1}$ is known one should have (order 1 in $\varepsilon$ and
$[|.|]$ is the jump function across $\partial S$ in the direction $\nu$):

$\left\\{\begin{array}[]{l}\forall{\mathbb{v}}\in W_{0},\\\ \\\
\alpha\displaystyle\int_{\partial
S}{\mathbb{u}}^{2}.{\mathbb{v}}+\beta\displaystyle\int_{S}\sigma({\mathbb{u}}^{2}).\gamma({\mathbb{v}})=-\displaystyle\int_{\Omega}\sigma({\mathbb{u}}^{1}).\gamma({\mathbb{v}})=\displaystyle\int_{\partial
S}{[|}(\sigma({\mathbb{u}}^{1}){|]}.\nu.{\mathbb{v}}.\end{array}\right.$ (36)

The restriction of ${\mathbb{u}}^{2}$ to $\Omega_{S}$ is denoted by
${\mathbb{u}}^{2p}$. The existence and uniqueness of this term
${\mathbb{u}}^{2p}\in{\Huge[}H^{1}(S){\Huge]}^{2},\,\hbox{div(${\mathbb{u}}^{2p}$)$=0$}$
on $S$, are obtained exactly as we did for ${\mathbb{u}}^{1p}$. Then the
computation of ${\mathbb{u}}^{2}$ on $\Omega_{S}$ is also handled as we did
for ${\mathbb{u}}^{1}$. Therefore ${\mathbb{u}}^{2}$ is now defined on the
whole domain $\Omega$ and is an element of the space $W_{0}$. Here again, no
additional regularity is required for the computation of ${\mathbb{u}}^{2}$
but we shall have a different conclusion for $\beta=0$.

######  Remark 9.

If $\beta=\gamma=0$ and $\alpha>0$ the term ${\mathbb{u}}^{1}$ should be
solution of (terms of order zero in the assumed asymptotic expansion in
$\varepsilon$):

$\forall{\mathbb{v}}\in
W_{0},\;\displaystyle\int_{\Omega}\sigma({\mathbb{u}}^{0}).\gamma({\mathbb{v}})+\alpha\displaystyle\int_{\partial
S}{\mathbb{u}}^{1}.{\mathbb{v}}=\displaystyle\int_{\Gamma_{1}}{\mathbb{g}}.{\mathbb{v}}.$
(37)

This implies that necessarily one should have (here
$\sigma({\mathbb{u}}^{0}).\nu$ is the value from $\Omega_{S}$ and it is zero
for the contribution from $S$ because ${\mathbb{u}}^{0}$ is a rigid body
velocity on $S$):

${\mathbb{u}}^{1}=-\sigma({\mathbb{u}}^{0}).\nu\;\hbox{ on }\partial
S\;(\nu\hbox{ inward }S).$ (38)

It is necessary to prove that the normal component of the stress field
$\sigma({\mathbb{u}}^{0}).\nu$ is more regular than the classical result (see
for instance J.L. Lions-E. Magenes [29]) which states that it belongs to the
space ${\Huge[}H^{-1/2}(\partial S){\Huge]}^{2}$. In fact a hidden regularity
result that one can derive from the J. Hadamard domain derivation method which
has been widely used in mathematical control analysis [28], enables one to
prove that $\sigma({\mathbb{u}}^{0}).\nu\in{\Huge[}L^{2}(\partial
S){\Huge]}^{2}$. But this is insufficient because ${\mathbb{u}}^{1}$ should
belong to the space ${\Huge[}H^{1/2}(\partial S){\Huge]}^{2}$ in order to have
${\mathbb{u}}^{1}$ in the space ${\Huge[}H^{1}(S){\Huge]}^{2}$. Therefore one
can point out a weakness of this case which is the lost of some regularity
which will have some nasty consequences on the convergence of the penalty
method when $\varepsilon\rightarrow 0$. This is the case for instance if the
boundary of the structure has sharp corners (see P. Grisvard for details
[23]). But the penalty-duality method will overcome this difficulty and could
be seen as an improvement of this penalty method when $\alpha>0$ and
$\beta=\gamma=0$ as we show in the following.

Let us now turn to the characterization of ${\mathbb{u}}^{1}$ inside
$\Omega_{S}$ and $S$ (always with $\alpha>0$ and $\beta=\gamma=0$). Hence from
the equation at the order one in $\varepsilon$, one obtains:

$\forall{\mathbb{v}}\in
W_{0},\;\displaystyle\int_{\Omega}\sigma({\mathbb{u}}^{1}).\gamma({\mathbb{v}})+\alpha\displaystyle\int_{\partial
S}{\mathbb{u}}^{2}.{\mathbb{v}}=0\;\hbox{ and
}{\mathbb{u}}^{1p}=-\sigma({\mathbb{u}}^{0}).\nu\;\hbox{ on }\partial S.$ (39)

Here again, as in Remark 8, the system (37) can be split into two independent
problems. which characterizes ${\mathbb{u}}^{1}$ separately on $\Omega_{S}$
and on $S$.

$\left\\{\begin{array}[]{l}\hbox{ on }\Omega_{S}:\\\ \\\
\forall{\mathbb{v}}\in
V_{0},\;\displaystyle\int_{\Omega_{S}}\sigma({\mathbb{u}}^{1}).\gamma({\mathbb{v}})=0,\;\;{\mathbb{u}}^{1}={\mathbb{u}}^{1p}\hbox{
on }\partial S,\\\ \\\ \hbox{ on }S:\\\ \\\
\forall{\mathbb{v}}\in{\Huge[}H^{1}_{0}(S){\mathbb{]}}^{2},\;\displaystyle\int_{S}\sigma({\mathbb{u}}^{1}).\gamma({\mathbb{v}})=0\;\hbox{
and }{\mathbb{u}}^{1}={\mathbb{u}}^{1p}\hbox{ on }\partial
S.\end{array}\right.$ (40)

These two non-homogeneous boundary values problems (Dirichlet condition on
$\partial S$) have each one a unique solution because of the classical
properties of the linear elasticity operator. But in this case there are two
important features which should be emphasized.

  1. 1.

The regularity of the boundary of the immersed structure ($\partial S$) is
required in order to ensure that
$\sigma({\mathbb{u}}^{0}).\nu\in{\Huge[}H^{1/2}(\partial S){\Huge]}^{2}.$ If
it isn’t one could make use singular functions due to the corner on $\partial
S$ but this is beyond our scope in this paper.

  2. 2.

The solution ${\mathbb{u}}^{1}$ on $S$ is no more a rigid body velocity (in
general).

Nevertheless this method ($\alpha>0$ and $\beta=\gamma=0$) will be the most
adapted to the penalty-duality method that we discuss in section 5 as far as
one only considers quasi-static modeling. $\Box$

###  3.3 Convergence results (case $\alpha>0\>\beta>0$ and $\gamma=0$).

We introduce the error term:

$\underline{\mathbb{u}}^{\varepsilon}={\mathbb{u}}^{\varepsilon}-{\mathbb{u}}^{0}-\varepsilon{\mathbb{u}}^{1}-\varepsilon^{2}{\mathbb{u}}^{2}\in
W_{0}.$ (41)

From the definition of the different terms and from the triangular
inequality222$\forall a,b\in{\mathbb{R}},\;\xi\in{\mathbb{R}}^{+}\hbox{ one
has }ab\leq\frac{2\xi}{a}^{2}+\frac{1}{2\xi}b^{2}$ and $c$ denoting the
continuity constant of the bilinear form of the elasticity operator:

$\hskip-5.69046pt\left\\{\begin{array}[]{l}\forall\xi>0:\;\;\displaystyle\int_{\Omega}\hskip-1.42271pt\sigma(\underline{\mathbb{u}}^{\varepsilon}).\gamma(\underline{\mathbb{u}}^{\varepsilon})+\displaystyle\frac{\alpha}{\varepsilon}\displaystyle\int_{\partial
S}\hskip-1.42271pt||\underline{\mathbb{u}}^{\varepsilon}||^{2}+\displaystyle\frac{\beta}{\varepsilon}\displaystyle\int_{S}\hskip-1.42271pt\sigma(\underline{\mathbb{u}}^{\varepsilon}).\gamma(\underline{\mathbb{u}}^{\varepsilon})=\\\
\\\
-\varepsilon^{2}\displaystyle\int_{\Omega}\hskip-1.42271pt\sigma({\mathbb{u}}^{2}).\gamma(\underline{\mathbb{u}}^{\varepsilon})\leq
c\displaystyle\left\\{\displaystyle\frac{\xi}{2}\varepsilon^{4}||{\mathbb{u}}^{2}||_{1,\Omega}^{2}+\displaystyle\frac{1}{2\xi}||\underline{\mathbb{u}}^{\varepsilon}||_{1,\Omega}^{2}\right\\}.\end{array}\right.$
(42)

Finally, from classical coerciveness properties deduced from Korn inequality
[18], one has the following error estimates where $c_{0}$ and $c_{1}$ are two
positive constants independent on $\varepsilon$ (but dependent on the data of
the problem):

$c_{0}||\underline{\mathbb{u}}^{\varepsilon}||_{1,\Omega}^{2}+\displaystyle\frac{\alpha}{2\varepsilon}||\underline{\mathbb{u}}^{\varepsilon}||_{0,\partial
S}^{2}+\displaystyle\frac{\beta}{2\varepsilon}||\underline{\mathbb{u}}^{\varepsilon}||_{1,S}^{2}\leq
c_{1}\varepsilon^{4}.$ (43)

Let us summarize the previous results in the following statement.

######  Theorem 3.2.

Let ${\mathbb{u}}^{\varepsilon}$ be the solution of (10) and
${\mathbb{u}}^{0},{\mathbb{u}}^{1}$ the terms characterized in the asymptotic
expansion at (33) and (35). We also assume in this Theorem that $\alpha>0$ and
$\beta>0$. Then there exists a constant -say $c_{2}>0$\- independent on
$\varepsilon$ and such that:

$\left\\{\begin{array}[]{l}||{\mathbb{u}}^{\varepsilon}-{\mathbb{u}}^{0}-\varepsilon{\mathbb{u}}^{1}||_{1,\Omega_{S}}\leq
c_{2}\varepsilon^{2},\\\ \\\
||{\mathbb{u}}^{\varepsilon}-{\mathbb{u}}^{0}-\varepsilon{\mathbb{u}}^{1}||_{0,\partial
S}\leq
c_{2}\varepsilon^{5/2},\;\;||{\mathbb{u}}^{\varepsilon}-{\mathbb{u}}^{0}-\varepsilon{\mathbb{u}}^{1}||_{1,S}\leq
c_{2}\varepsilon^{5/2}.\end{array}\right.$ (44)

Furthermore the term ${\mathbb{u}}^{0}$ is exactly the solution of the model
(3) we started from where the displacement of the structure is the prescribed
to the boundary $\partial S$ of the fluid.

######  Corollary 1.

From Theorem 3.2 and making use of the triangular inequality, one obtains the
following error bounds between ${\mathbb{u}}^{\varepsilon}$ and
${\mathbb{u}}^{0}$ where $c_{3}$ is a constant independent on $\varepsilon$:

$\left\\{\begin{array}[]{l}||{\mathbb{u}}^{\varepsilon}-{\mathbb{u}}^{0}||_{1,\Omega_{S}}\leq
c_{3}\varepsilon,\\\ \\\
||{\mathbb{u}}^{\varepsilon}-{\mathbb{u}}^{0}||_{0,\partial S}\leq
c_{3}\varepsilon,\;\;||{\mathbb{u}}^{\varepsilon}-{\mathbb{u}}^{0}||_{1,S}\leq
c_{3}\varepsilon.\end{array}\right.$ (45)

######  Remark 10.

In the eventuality where $\alpha=0$ (no penalty term on the boundary of $S$)
the similar results as those of Theorem 3.2 are still true but
${\mathbb{u}}^{0}$ is the not the solution of the initial problem. Even if it
is a rigid body velocity on $S$, it is not the one of the structure deduced
from a correct writing of the principle of the mechanics. For instance the
mass of the structure is not appearing in the model. Furthermore we restrict
the analysis to quasi-static fluid-structure interaction (small reduce
frequencies). Hence this possibility should be forgotten. $\Box$

###  3.4 Convergence result for $\beta=\gamma=0$ and $\alpha>0$.

The results obtained in Theorem 3.2 are still true up to few modifications
implying an additional regularity assumption on the term ${\mathbb{u}}^{0}$.
We just give the statement without details because of the narrow similarity
with the one of Theorem 3.2. Similarly to the error bound derived for
$\alpha>0$ and $\beta>0$, and as far as we just upper bound the error between
${\mathbb{u}}^{\varepsilon}$ and ${\mathbb{u}}^{0}$, it is not necessary to
imply the term ${\mathbb{u}}^{2}$. At the opposite of the case $\beta>0$, it
would require also an additional regularity on $\sigma({\mathbb{u}}^{1}).\nu$
and therefore on ${\mathbb{u}}^{0}$ too.

######  Theorem 3.3.

The notations are those introduced in the previous subsection but
${\mathbb{u}}^{0}$ is now solution of (22). Furthermore, it is assumed that
the solution ${\mathbb{u}}^{0}$ of the equations (22) satisfies
$\sigma({\mathbb{u}}^{0}).\nu\in{\Huge[}H^{1/2}(\partial S){\Huge]}^{2}$. Then
there exists a constant $c_{4}>0$ independent on $\varepsilon$ such that:

$\left\\{\begin{array}[]{l}||{\mathbb{u}}^{\varepsilon}-{\mathbb{u}}^{0}||_{1,\Omega}\leq
c_{4}\varepsilon,\\\ \\\
||{\mathbb{u}}^{\varepsilon}-{\mathbb{u}}_{0}||_{0,\partial S}\leq
c_{4}\varepsilon.\end{array}\right.$ (46)

######  Remark 11.

One should notice that the convergence on $S$ to the rigid body velocity is in
$\varepsilon$ as in the case where $\beta>0$. The correctors
${\mathbb{u}}^{1}$ have in fact no real interest because the term that we need
to estimate is ${\mathbb{u}}^{0}$. Let us also emphasize on the interest of
this strategy which will clearly appear in the penalty-duality method. $\Box$

##  4 A simple numerical discussion between the various penalty methods

In order to compare the different penalty terms we discuss a very simple
example which shows the various behavior of the three penalty terms considered
in this paper.

###  4.1 The test case

The model is the following one:

$\hskip-14.22636pt\left\\{\begin{array}[]{l}V=\\{v\in
H^{1}(]0,L[),\;v(0)=0\\},\;\hbox{ find }u\in V\;\hbox{ that:
}\displaystyle\forall v\in V,\\\ \\\
\displaystyle\int_{0}^{L}\hskip-1.42271pt\displaystyle\frac{du}{dx}\displaystyle\frac{dv}{dx}\hskip-1.42271pt+\hskip-1.42271pt\displaystyle\frac{\alpha}{\varepsilon}\left(u\left(\frac{\tiny
L}{\tiny
2}\right)\hskip-1.42271pt-\hskip-1.42271ptu_{0}\right)v\left(\frac{L}{2}\right)\hskip-1.42271pt+\hskip-1.42271pt\displaystyle\int_{0}^{L}\hskip-1.42271pt\left\\{\displaystyle\frac{\beta}{\varepsilon}\displaystyle\frac{du}{dx}\displaystyle\frac{dv}{dx}\hskip-1.42271pt+\hskip-1.42271pt\displaystyle\frac{\gamma}{\varepsilon}\left(u-u_{0}\right)v\right\\}\hskip-1.42271pt=\hskip-1.42271pt0.\end{array}\right.$
(47)

The term $u_{0}$ is a given constant and $\varepsilon$ is the penalty
parameter in order to prescribe approximately $u=u_{0}$ on $]L/2,L[$. The
three parameters $\alpha,\beta$ and $\gamma$ are chosen equal to $0$ or $1$ in
order to compare their efficiency in the penalty strategy. The solution method
used is a first degree finite element method (1000 points). We have plotted
the results obtained for two choices of the penalty parameter $\varepsilon$ on
Figures 2-3. The space derivative of the solution is plotted on Figure 4 for
$\varepsilon=10^{-2}$

###  4.2 Comments on the results

First of all it is worth noting that the method with $\alpha=\gamma=0$ and
$\beta>0$ is not adapted to prescribe a given movement of the structure. As it
is explained in remark 10, the solution is constant on $]L/2,L[$ but it is not
the one that one hopes. Nevertheless, if we consider the two terms $\alpha>0$
and $\beta\geq 0$ or $\beta>0$ and $\gamma>0$ the results are are much better.
But in the second case the converge is not as good as than in the first one.
If one restricts the penalty term to $\alpha=\beta=0$ and $\gamma>0$ the
convergence occurs but is slower than in the case where $\alpha>0$, even if
this is the only term which is kept.

Concerning the derivatives of the solution with respect to the space
coordinate one can see on Figure 4 that the convergence is not satisfying if
$\alpha=0$. This the reason why we develop in the next section 5 a penalty-
duality method which implies only this boundary term. Furthermore, the dual
variable (Lagrange multiplier) is easier to handle when only this boundary
term is taken into account. In case of internal conditions the duality is a
little bit more complex. We refer to the so-called Arlequin method developed
by H. Ben Dhia [6]-[7] and in slightly different context by T. Belytschko &
all [3].

Figure 2: Penalty solution for $\varepsilon=0.1$ Figure 3: Penalty solution
for $\varepsilon=0.01$ Figure 4: Derivative in $x$ of the penalty solution for
$\varepsilon=0.01$

##  5 The penalty-duality method

In this section, we consider the case where $\alpha>0$ and $\beta=\gamma=0$.
Because the penalty model leads to ill conditioning, one can suggest to use a
variant of this strategy introduced fifty years ago by D. Bertsekas [9] and
known as the penalty-duality method. The point is to satisfy exactly the
constraints which are considered in the penalty term by a dual treatment.
Hence, for any $r>0$ ($r$ stands for $1/\varepsilon$ but will be moderate at
the opposite of what happens in the penalty strategy), we introduce the
Lagrangian defined by:

$L({\mathbb{v}},\lambda,\mu)=\displaystyle\frac{1}{2}\displaystyle\int_{\Omega}\sigma({\mathbb{v}}).\gamma({\mathbb{v}})+\displaystyle\frac{\alpha
r}{2}\displaystyle\int_{\partial
S}||{\mathbb{v}}-{\mathbb{u}}_{0}||^{2}+\displaystyle\int_{\partial
S}\mu.({\mathbb{v}}-{\mathbb{u}}_{0})-\displaystyle\int_{\Gamma_{1}}{\mathbb{g}}.{\mathbb{v}},$
(48)

where ${\mathbb{v}}\in W_{0},\;\mu\in{\Huge[}H^{-1/2}(\partial S){\Huge]}^{2}$
(we use the notation with the integral for the duality333the dual of
$H^{1/2}(\partial S)$ is $H^{-1/2}(\partial S)$ in $1D$ which is the case
here. but we know that it is an abuse; see remark 12). One can add the penalty
term inside $S$ (see [4]-[5]-[27]):

$\displaystyle\frac{\gamma}{\varepsilon}\displaystyle\int_{S}|{\mathbb{u}}-{\mathbb{u}}_{0}|^{2}.$

It is not really useful for the static case that we consider here but it is a
necessity for the dynamic case as those treated in the numerical tests of
subsections 8.2 and 8.3.

A saddle point of $L$ is an element $({\mathbb{u}}^{r},\lambda)\in
W_{0}\times{\Huge[}H^{-1/2}(\partial S){\Huge]}^{2}$ such that:

$\left\\{\begin{array}[]{l}\forall{\mathbb{v}}\in
W_{0},\;\displaystyle\int_{\Omega}\sigma({\mathbb{u}}^{r}).\gamma({\mathbb{v}})+{\alpha
r}\displaystyle\int_{\partial
S}({\mathbb{u}}^{r}-{\mathbb{u}}_{0}).{\mathbb{v}}+\displaystyle\int_{\partial
S}\lambda.{\mathbb{v}}=\displaystyle\int_{\Gamma_{1}}{\mathbb{g}}.{\mathbb{v}},\\\
\\\ \forall\mu\in{\Huge[}H^{-1/2}(\partial
S){\Huge]}^{2},\;\displaystyle\int_{\partial
S}\mu.({\mathbb{u}}^{r}-{\mathbb{u}}_{0})=0,\par\end{array}\right.$ (49)

######  Remark 12.

The writing with an integral of the duality between the spaces
${\Huge[}H^{-1/2}(\partial S){\Huge]}^{2}$ and ${\Huge[}H^{-1/2}(\partial
S){\Huge]}^{2}$ is a familiarity which is not correct because the functions
involved are not in the space $L^{2}(\partial S)$. It just a commodity (in our
mind). $\Box$

The existence and uniqueness are relevant of the general A.K. Aziz and I.
Babuska Theorem [2] that we recall hereafter for sake of convenience.

######  Theorem 5.1.

[Aziz-Babuska] Let $H_{1}$ and $H_{2}$ be two Hilbert spaces, $b$ a bilinear
and continuous form on $H_{1}\times H_{2}$ and $j$ a linear and continuous
form on $H_{2}$. It is asuumed that the following properties hold:

  1. #

1 if $q\in H_{2}$ is such that: $\forall\mu\in H_{1},\;b(\mu,q)=0$ then $q=0$;

  2. #

2 there exists a constant $\delta>0$ such that $\forall\mu\in
H_{1},\;\displaystyle\sup_{q\in
H_{2}}\displaystyle\frac{b(\mu,q)}{||q||_{H_{2}}}\geq\delta||\mu||_{H_{1}}.$

Then there exists a unique element $\lambda\in H_{1}$ such that:

$\forall q\in H_{2},\;b(\lambda,q)=j(q).$ (50)

Let us now apply this Theorem to our case.

######  Theorem 5.2.

Existence and uniqueness of a solution to (49). Let $r>0$ and $\alpha>0$. The
function ${\mathbb{g}}$ is assumed as previously to be in the space
${\Huge[}L^{2}(\Gamma_{1}){\Huge]}^{2}$. Hence the system (46) has a unique
solution $({\mathbb{u}}^{r},\lambda)$in the space
$W_{0}\times\left[H^{-1/2}(\partial S)\right]^{2}$. Furthermore
${\mathbb{u}}^{r}={\mathbb{u}}$ solution of the initial problem (3).

###### Proof.

If $({\mathbb{u}}^{r},\lambda)\in W_{0}\times\left[H^{-1/2}(\partial
S)\right]^{2}$ is a solution it should satisfy the following relation:

$\forall{\mathbb{v}}\in
W_{00},\;\displaystyle\int_{\Omega}\sigma({\mathbb{u}}^{r}).\gamma({\mathbb{v}})=\displaystyle\int_{\Gamma_{1}}{\mathbb{g}}.{\mathbb{v}},\;\hbox{
and necessarily:\;}{\mathbb{u}}^{r}={\mathbb{u}}_{0}\hbox{ on }\partial S.$
(51)

Therefore ${\mathbb{u}}^{r}$ is the initial solution of equation (3) which is
also the term ${\mathbb{u}}^{0}$ computed in the asymptotic method(s) applied
to the penalty method at section 3 (for both $\alpha>0,\;\beta=0$ or
$\alpha>0,\beta>0$). It also implies that it is unique and equal to the rigid
body velocity ${\mathbb{u}}_{0}$ on the whole domain $S$ and not only on its
boundary $\partial S$.

Finally, in order to complete the proof of Theorem 5.2 it is sufficient to
prove the existence and uniqueness of $\lambda$ using Aziz-Babuska Theorem
5.1.

Let us set:

$\left\\{\begin{array}[]{l}H_{1}=\left[H^{-1/2}(\partial
S)\right]^{2},\;H_{2}=\left[H^{-1/2}(\partial S)\right]^{2},\\\ \\\
\forall\mu\in H_{1},\;\forall{\mathbb{v}}\in
H_{2},\;b(\mu,{\mathbb{v}})=\alpha\displaystyle\int_{\partial
S}\mu.{\mathbb{v}},\\\ \\\ \forall{\mathbb{v}}\in
H_{2},\;j({\mathbb{v}})=\displaystyle\int_{\Omega_{S}}\sigma({\mathbb{u}}^{r}).\gamma({\mathbb{v}})-\displaystyle\int_{\Gamma_{1}}{\mathbb{g}}.{\mathbb{v}}.\end{array}\right.$
(52)

The first requirement of the Aziz-Babuska Theorem is obvious because of the
duality between $H_{1}$ and $H_{2}$. The second one is satisfied with
$\delta=1$ because of the definition of the norm in $H_{1}$.

The only point to be verified is the definition and the continuity of $j$ on
$H_{2}$. First of all $j$ is clearly linear. It depends on the value of
${\mathbb{v}}$ on the boundary $\partial S$ because of the definition of
${\mathbb{u}}^{r}={\mathbb{u}}^{0}$. Hence, for a given element ${\mathbb{v}}$
in $H_{2}$ one can associate any element in the space $V_{0}$ denoted also by
${\mathbb{v}}$ and which gives the same value of $j$. Let us recall that there
exists a linear and continuous operator -say $R$\- from $H_{2}$ into $V_{0}$
[29] such that:

$\exists c_{9}>0,\hbox{ independent on }{\mathbb{v}}\in H_{1}\hbox{ suh that:
}||R({\mathbb{v}})||_{1,\Omega_{S}}\leq c_{9}||{\mathbb{v}}||_{1/2,\partial
S}.$ (53)

Hence:

$|j({\mathbb{v}})|=\displaystyle|\displaystyle\int_{\Omega_{S}}\sigma(R{\mathbb{v}}).\gamma(R{\mathbb{v}})-\displaystyle\int_{\Gamma_{1}}{\mathbb{g}}.R({\mathbb{v}})|\leq
c_{10}||{\mathbb{v}}||_{1/2,\partial S},$ (54)

Hence the assumptions of Theorem 5.1 are satisfied and we can conclude the
proof of Theorem 5.2. ∎

######  Remark 13.

At the opposite of what occurs in the penalty method, the solution of the
penalty-duality model is exactly the one of the initial problem. $\Box$

######  Remark 14.

It is a basic point to notice that the solution ${\mathbb{u}}^{r}$ is
independent on $r$. Therefore it is not necessary to choose a large value as
in the penalty model ($r$ is the equivalent of $1/{\varepsilon}$) .
Furthermore one could point out that the strategy could also be applied
without the penalty term but in this case the algorithm studied in the next
section is more difficult to use as we underline in the following (difficulty
in adjusting the gradient step in the Uzawa solution method [9], [21]). If one
uses a conjugate gradient on the dual problem (in $\lambda$), the operator has
a better condition number if one add the penalty term. $\Box$

######  Remark 15.

A possibility that we do not recommend consists in considering the two penalty
terms ($\alpha>0$ and $\beta>0$). Nevertheless the dual space for the Lagrange
multiplier of the the constraint $\gamma_{ij}({\mathbb{u}})=0$ on $S$ is a
little bit complex and not very convenient in a practical application.
Furthermore the numerical implementation seems to be some more complicated.
$\Box$

##  6 The numerical algorithm and its solution method

In this section we focus on the Uzawa algorithm [21][11] for solving (49). Let
us start from a given value for $\lambda^{n_{0}}\in{\Huge[}H^{-1/2}(\partial
S){\Huge]}^{2}$ (the most classical choice is to start from
$\lambda^{n_{0}}=0$). For each $n\geq n_{0}$, we define ${\mathbb{u}}^{n}\in
W_{0}$ solution of:

$\forall{\mathbb{v}}\in
W_{0},\;\displaystyle\int_{\Omega}\sigma({\mathbb{u}}^{n}).\gamma({\mathbb{v}})+{r}\displaystyle\int_{\partial
S}({\mathbb{u}}^{n}-{\mathbb{u}}_{0}).{\mathbb{v}}=\displaystyle\int_{\Gamma_{1}}{\mathbb{g}}.{\mathbb{v}}-\displaystyle\int_{\partial
S}\lambda^{n}.{\mathbb{v}}.$ (55)

Then we upgrade $\lambda^{n}$ by setting:

$\lambda^{n+1}=\lambda^{n}+r({\mathbb{u}}^{n}-{\mathbb{u}}_{0})\;\hbox{ on
}\partial S.$ (56)

An important point is that if the first multiplier $\lambda^{n}$ belongs to
the space ${\Huge[}L^{2}(\partial S){\Huge]}^{2}$ all the sequence
$\lambda^{n}$ belongs to this space. But in general the solution $\lambda$ of
(49) doesn’t belongs to this space unless one has a regularity property. In
fact, one has (interpretation of (49)),
$\lambda=-\sigma({\mathbb{u}}^{r}).\nu$ on $\partial S$. This regularity
assumption is exactly the one used in the asymptotic analysis of the penalty
model for $\alpha>0$ and $\beta=0$. It is convenient (but not necessary) in
the following to adopt this hypothesis (the norm $L^{2}(\partial S)$ is more
convenient in the writings than the one of $H^{-1/2}(\partial S)$).

The convergence of the algorithm is classical [9]. But for sake of clarity let
us summarize the proof hereafter with few remarks concerning the regularity of
the multiplier $\lambda$ which are only valid for our case.

######  Theorem 6.1.

Let $({\mathbb{u}}^{n},\lambda^{n})\in W_{0}\times{\Huge[}H^{-1/2}(\partial
S){\Huge]}^{2}$ be the sequence defined by the algorithm (55)-(56). We assume,
just for sake of convenience in the writings, that:

$\lambda=-\sigma({\mathbb{u}}^{r}).\nu\in{\Huge[}L^{2}(\partial
S){\Huge]}^{2}.$

Then for any $r>0$:

$\left\\{\begin{array}[]{l}\displaystyle\lim_{n\rightarrow\infty}\;{\mathbb{u}}^{n}={\mathbb{u}}^{r}={\mathbb{u}}\hbox{
solution of (\ref{eq3})},\\\ \\\
\displaystyle\lim_{n\rightarrow\infty}\;\lambda^{n}=\lambda\hbox{ solution
with ${\mathbb{u}}^{r}$ of (\ref{eq27})}.\end{array}\right.$ (57)

###### Proof.

Let us introduce the gap variables:

$\underline{\mathbb{u}}^{n}={\mathbb{u}}^{n}-{\mathbb{u}}^{r},\;\;\underline{\lambda}^{n}=\lambda^{n}-\lambda.$

One has the following equalities:

$\left\\{\begin{array}[]{l}\forall{\mathbb{v}}\in
W_{0},\;\displaystyle\int_{\Omega}\sigma(\underline{\mathbb{u}}^{n}).\gamma(\underline{\mathbb{u}}^{n})+{r}\displaystyle\int_{\partial
S}|\underline{\mathbb{u}}^{n}|^{2}+\displaystyle\int_{\partial
S}\underline{\mathbb{u}}^{n}\underline{\lambda}^{n}=0,\\\ \\\
|\underline{\lambda}^{n+1}|^{2}=|\underline{\lambda}^{n}|^{2}+{r^{2}}|\underline{\mathbb{u}}^{n}|^{2}+{2r}\underline{\lambda}^{n}.\underline{\mathbb{u}}^{n}.\end{array}\right.$

We deduce that:

$\displaystyle\int_{\partial
S}|\underline{\lambda}^{n}|^{2}-\displaystyle\int_{\partial
S}|\underline{\lambda}^{n+1}|^{2}=r{\Huge[}-2\displaystyle\int_{\partial
S}\underline{\mathbb{u}}^{n}\underline{\lambda}^{n}-r\displaystyle\int_{\partial
S}|\underline{\mathbb{u}}^{n}|^{2}{\Huge]},$

or else:

$\begin{array}[]{l}\displaystyle\int_{\partial
S}|\underline{\lambda}^{n}|^{2}-\displaystyle\int_{\partial
S}|\underline{\lambda}^{n+1}|^{2}=r{\Huge[}2\displaystyle\int_{\Omega}\sigma(\underline{\mathbb{u}}).\gamma(\underline{\mathbb{u}}^{n})+2r\displaystyle\int_{\partial
S}|\underline{\mathbb{u}}^{n}|^{2}-r\displaystyle\int_{\partial
S}|\underline{\mathbb{u}}^{n}|^{2}{\Huge]}\\\ \\\
=r{\Huge[}2\displaystyle\int_{\Omega}\sigma_{ij}(\underline{\mathbb{u}}).\gamma(\underline{\mathbb{u}}^{n})+r\displaystyle\int_{\partial
S}|\underline{\mathbb{u}}^{n}|^{2}{\Huge]}\geq 0.\end{array}$

Because the real and positive sequence of numbers
$q^{n}=\displaystyle\int_{\partial S}|\underline{\lambda}^{n}|^{2}$ is
decreasing, it is convergent and therefore it is a Cauchy sequence. Hence
$\displaystyle\lim_{n\rightarrow\infty}q^{n}-q^{n+1}=0.$ As a consequence
(using the $H^{1}(\Omega)$-coerciveness of the bilinear form
$\displaystyle\int_{\Omega}\sigma({\mathbb{u}}).\gamma({\mathbb{v}})$):

$\left\\{\begin{array}[]{l}\displaystyle\lim_{n\rightarrow\infty}||{\mathbb{u}}^{n}-{\mathbb{u}}_{0}||_{0,\partial
S}=0,\\\ \\\
\displaystyle\lim_{n\rightarrow\infty}||{\mathbb{u}}^{n}-{\mathbb{u}}^{r}||_{1,\Omega}=0.\end{array}\right.$

The second convergence result reinforces the first one because of the
continuity of the trace mapping from $H^{1}(\Omega)$ into $H^{1/2}(\partial
S)$:

$\displaystyle\lim_{n\rightarrow\infty}||{\mathbb{u}}^{n}-{\mathbb{u}}_{0}||_{1/2,\partial
S}=0.$

The final step consists in proving the convergence of the multiplier
$\lambda^{n}$. From:

$\displaystyle\int_{\partial
S}\underline{\lambda}^{n}.{\mathbb{v}}=-\displaystyle\int_{\Omega}\sigma(\underline{\mathbb{u}}^{n}).\gamma({\mathbb{v}})-r\displaystyle\int_{\partial
S}\underline{\mathbb{u}}^{n}.{\mathbb{v}}-\displaystyle\int_{\Gamma_{1}}{\mathbb{g}}.{\mathbb{v}},$

and because of the property $\\#$ 2 mentioned in Theorem 5.1, one has:

$\displaystyle\lim_{n\rightarrow\infty}||\lambda^{n}-\lambda||_{-1/2,\partial
S}=0.$

Even if we assumed the regularity $\lambda\in{\huge[}L^{2}(\partial
S){\Huge]}^{2}$ the standard convergence is only in the space
${\huge[}H^{-1/2}(\partial S){\Huge]}^{2}$. ∎

######  Remark 16.

The main advantage in the choice of the penalty term reduced on the boundary
of the structure is that it enables to decouple completely the numerical
software and just requires a geometrical prolongation mapping from the
boundary $\partial S$ of the structure into the open set $\Omega$ and a
geometrical restriction mapping from $\Omega$ on to $\partial S$. This
advantage was underlined by many authors and a nice presentation by J.
Hovnanian is given in [24]. $\Box$

##  7 Few elementary tests for the penalty-duality strategy

Just in order to make a link with the elementary test introduced in section 4,
we apply the penalty-duality algorithm to the same model. The results are
plotted on Figures 5-6. They show the efficiency of this method regarding the
precision even for $\varepsilon$ not so small. The numerical implementation in
a complex model is certainly another discussion which will be carried out in
A. Falaise and E. Liberge [19]. In particular, the mapping connecting the
boundary $\partial S$ with the global mesh is a cornerstone problem regarding
the multiprocessor programming.

Figure 5: Solution obtained with the penalty-duality strategy with
$\varepsilon=.1$ Figure 6: Derivative of the solution obtained with the
penalty-duality strategy with $\varepsilon=.1$ during the four first
iterations

##  8 Application to convection-diffusion equation

We consider two examples: the first one is the linear advection-diffusion
model and the second one is the Burgers equation. Nevertheless it is worth to
point out that the boundary of the structure are not updated..

###  8.1 The linear advection diffusion model

We consider the following $1D$ advection-diffusion equation where
$u\in{\mathcal{C}}^{0}([0,T]\times[0,L])$ is the unknown (where $c>0$
$u_{0}(t)$ and $d(x)$ are given):

$\left\\{\begin{array}[]{l}x\in]0,a[\cup]b,L[,\;t\geq
0:\;\displaystyle\frac{\partial u}{\partial t}+c\displaystyle\frac{\partial
u}{\partial x}-\nu\displaystyle\frac{\partial^{2}u}{\partial x^{2}}=0,\\\ \\\
t\geq 0\;u(0,t)=u(L,t),\;\;x\in]0,a[\cup]b,L[:\;u(x,0)=u_{0}(x),\\\ \\\ t\geq
0,\;x\in]a,b[:\;u(x,t)=d_{0}(t)\;(\hbox{ sine function}).\end{array}\right.$
(58)

We denote by $\Delta x$ (respectively $\Delta t$) the space step (respectively
the time step). The lower index indicate the space dicretization and the upper
one the time discretization. The numerical scheme that we consider for this
equation in $(]0,a[\cup]b,L[)\times]0,T[$ leads to:

$\left\\{\begin{array}[]{l}\forall i\in\\{2,i_{a}\\}\cup\\{i_{b},N-1\\},\\\
\\\ \displaystyle\frac{u^{n+1}_{i}-u^{n}_{i}}{\Delta
t}+c(\displaystyle\frac{u_{i}^{n}-u^{n}_{i-1}}{\Delta
x})-\nu(\displaystyle\frac{u^{n+1}_{i+1}-2u^{n+1}_{i}+u^{n+1}_{i-1}}{\Delta
x^{2}})=0,\\\ \\\ u_{i}^{0}=u_{0}(i\Delta x)\hbox{ initial conditions,}\\\ \\\
\forall n\geq 0:\;u_{1}^{n}=u_{N}^{n},\hbox{ boundary conditions, }\\\ \\\
\;u_{i_{a}}^{n}=u_{i_{b}}^{n}=d_{0}^{n}\;\hbox{ boundary of the structure},\\\
\\\ \hbox{ where }i_{a}\hbox{ and }i_{b}\;\hbox{ are the indices corresponding
to $a$ and $b$}.\end{array}\right.$ (59)

The penalty-duality formulation is the following one ($\delta_{i,i_{c}}$ is
the Kronecker symbol $=1$ if $i=i_{c}$ and $0$ else and we introduce a penalty
term with the parameter $\gamma$ as before):

$\left\\{\begin{array}[]{l}\forall i\in\\{2,N-1\\},\\\ \\\
\displaystyle\frac{u^{n+1}_{i}-u^{n}_{i}}{\Delta
t}+\displaystyle\frac{\gamma}{\varepsilon}(u_{i}^{n+1}-d_{0}^{n+1})+c(\displaystyle\frac{u_{i}^{n}-u^{n}_{i-1}}{\Delta
x})+r[(u_{i}^{n+1}-d_{0}^{n+1})\delta_{i,i_{a}}+\\\ \\\
(u_{i}^{n+1}-d_{0}^{n+1})\delta_{i,i_{b}}]-\nu(\displaystyle\frac{u^{n+1}_{i+1}-2u^{n+1}_{i}+u^{n+1}_{i-1}}{\Delta
x^{2}})+\\\ \\\
\lambda^{n+1}(i_{a})\delta_{i,i_{a}}+\lambda^{n+1}(i_{b})\delta_{i,i_{b}}=0,\;\;u_{i}^{0}=u_{0}(i\Delta
x)\hbox{ initial conditions,}\\\ \\\ \forall n\geq
0:\;u_{1}^{n}=u_{N}^{n}\hbox{ and }u_{-1}^{n}=u_{N-1}^{n}\hbox{ periodic
boundary conditions, }\\\ \\\ \forall n\geq
0:\;u_{i_{a}}^{n}=u_{i_{b}}^{n}=d_{0}^{n}\;\hbox{ boundary of the
structure}.\end{array}\right.$ (60)

For each $n$ the duality algorithm (the iterations are indexed by $p$) is
defined as follows:

  * #

$p=0,\;\lambda_{i_{a}}^{0}=\lambda_{i_{b}}^{0}=0$;

  * #

compute $u^{n+1}_{i}$ solution of (60)
$\lambda^{p}=(\lambda^{p}_{i_{a}},\lambda^{p}_{i_{b}})$ being fixed;

  * #

set:
$\lambda^{p+1}=\lambda^{p}+r(u_{i_{a}}^{n}-d_{0}^{n},u_{i_{b}}^{n}-d_{0}^{n})$

  * #

convergence test => stop or $p=p+1$

###  8.2 The numerical tests for the advection-diffusion (linear)

We have tested the two possibilities: a) no duality and b) with duality. The
results are plotted on Figures 7(a) and 7(b) for the first case and on Figures
8(a) and 8(b) for the case with duality. We choose the following data set:

  * #

the kinematical viscosity is $\nu=0.001$;

  * #

the coefficient of the Bertsekas algorithm is $r=10$;

  * #

the length of the space interval is $L=1$ and the length of the structure is
$L/10$ and the coordinates of the two extremities are denoted by
$x_{a}=.45L,\;x_{b}=.55L$;

  * #

the time dependence of the prescribed velocity of the structure is
${u}_{st}(t)=\sin(2\pi t)$ and $c=1$ for the advection-diffusion case. The
space dependence in the structure is (just for an example):

$u_{sx}(x)=0.4+2(2x-x_{a}-x_{b})/L.$

Hence the velocity prescribed in the structure is $u_{st}(t)u_{sx}(x)$;

  * #

the initial condition outside the structure is (still for example):

$\begin{array}[]{l}u(x,0)=4.4(x_{a}-x-L/5)\hbox{ if }x\in[0,x_{a}-L/5],\\\
u(x,0)=4.4(x-x_{b}-L/5)\hbox{ if }x\in[x_{b}+L/5,L]\;\hbox{ and $u(x,0)=0$
elsewhere;}\end{array}$

  * #

the penalty parameter used for the periodic conditions at $x=0$ and $x=L$ but
also for the penalty term inside the structure is $\varepsilon=10^{-3}$;

  * #

the time delay is $T=2$ and there are $1000$ time steps and 500 space steps.

There is a meaningful difference between the two cases and mainly on the
derivatives with respect to $x$ if $\varepsilon>10^{-8}$. For
$\varepsilon<10^{-8}$ the penalized solution is quite the one obtained with
the penalty-duality method. In other words the duality algorithm is not
useful. Nevertheless, for such values of the penalty parameter the solver
requires a large enough viscosity which is not so drastic for larger values of
$\varepsilon$. Furthermore the one dimensional test aren’t meaningful
concerning the condition number of two or three dimensional models. Our goal
in this paper is only to discuss the advantages and the drawbacks of this
added duality algorithm for fluid-structure models from a theoretical point of
view. Figures 7 show the solution (ad) (Figure 7(a)) and the derivative versus
$x$ of the solution (ad) (Figure 7(b) obtained without neither penalty inside
the structure nor duality, but with a penalty term on the boundary between the
structure and the fluid. The structure is moving at a given frequency (2
periods have been plotted). This induces waves in the fluid. One should notice
on Figure 7(b) that the propagation of stress waves is very damped in the
fluid. The solutions obtained with the duality method are plotted on Figure 8.
The normal stress at the two boundaries versus time, with and without duality
are plotted Figure 9.

(a) Solution

(b) The derivative versus $x$ of the solution (ad)

Figure 7: Solution (ad) obtained without neither penalty inside the structure
nor duality. But there is a penalty term on the boundary between the structure
and the fluid.

(a) Solution (ad)

(b) The derivative versus $x$ of the solution (ad)

Figure 8:  Solution (ad) obtained with the duality method on the boundary of
the structure but no penalty term inside the structure.

(a) without duality

(b) with duality

Figure 9: The normal stress at the two boundaries of the structure versus the
time

Next, the penalty term inside the structure is added. Result with and without
duality are plotted on Figures 10 and 11. The duality algorithm does not
affect the normal stress at the fluid-solid interface as far as $\varepsilon$
is very small (Figure 12) but the use of the penalty term makes an important
difference.

(a) Solution

(b) Space derivative of the solution of the advection-convection equation

Figure 10: Solution of the advection-convection equation with two penalty
terms: one on the boundary between the fluid and the structure and the other
inside the structure. But no duality iteration have been used.

(a) Solution

(b) Space derivative of the solution of the advection-convection equation

Figure 11: solution of the advection-convection equation with the two penalty
terms and the duality algorithm.

(a) no duality

(b) duality

Figure 12: This Figure represents the normal stress at the f-s interfaces; the
two penalty terms have been used. The difference with Figure 9 (no penalty
inside the structure) can be noticed.

###  8.3 A non linear case (Burgers)

In this subsection we consider the following non linear model (Burgers
equation) but the boundary of the structure aren’t updated. The scheme used is
similar to the one used is the linear case (fully implicit for the diffusion
and Godounov scheme for the transport) and the data are also the same
concerning the dimensions and the initial values.

$\left\\{\begin{array}[]{l}x\in]0,a[\cup]b,L[,\;t\geq
0:\;\displaystyle\frac{\partial u}{\partial
t}+\displaystyle\frac{\partial}{\partial
x}(\displaystyle\frac{u^{2}}{2})-\nu\displaystyle\frac{\partial^{2}u}{\partial
x^{2}}=0;\\\ \\\ t\geq
0\;u(0,t)=u(L,t),\;\;x\in]0,a[\cup]b,L[:\;u(x,0)=u_{0}(x)\\\ \\\ t\geq
0,\;x\in]a,b[:\;u(x,t)=d_{0}(t)\;(\hbox{ sine function}).\end{array}\right.$
(61)

####  8.3.1 The numerical results

First of all let us point out that the results should be compared with the
quasi-exact solution obtained with a very small value of the penalty parameter
$\varepsilon$ but which requires some precautions in the adjustment of the
numerical tests (ill conditioning). The quasi-exact solution and its space
derivatives are plotted on Figures 21(a)-21(b)-22(a)-22(b).

The first results have been plotted on Figures 13(a)-13(b)-14 for the case
where there is no penalty term inside the structure and no duality. But there
is a moderate penalty term at the f-s interfaces. If we compare to the quasi-
exact solution plotted on Figure 21(a) one can see a big difference. The same
is true (and even amplified for the space derivatives (see Figures 13(b)-14
compared to Figures 21(a)-22(a)).

On the Figures 15(a)-15(b)-16(a)-16(b) we have plotted the solution of Burgers
equation obtained with the penalty-duality algorithm but without penalty term
($L^{2}$-norm between the velocity obtained by the scheme and the one
prescribed in the structure). The results compared with the quasi-exact
solution are much better but not satisfying.

The penalty term inside the structure without duality has been added in the
computations shown on Figures 17(a)-17(b)-18. The results are closer to the
quasi-exact one of Figures 21(a)-21(b)-22(a). The improvement with adding the
duality algorithm is small for small values of the penalty parameter as one
can see on Figures 19(a)-19(b)-20(a)-20(b) (compared to the quasi-exact
solution), but exists.

(a) Solution (Burgers)

(b) Derivative versus $x$ of the solution (Burgers)

Figure 13: Solution (Burgers) without duality and no penalty term inside the
structure (same data as for the linear case) but $r=.1$ Figure 14: Normal
stress at the f-s interface for the Burgers solution without duality and no
penalty term inside the structure (same data as for the linear case)

(a) Solution

(b) Derivative versus $x$

Figure 15: The solution with duality (Burgers) but no penalty term at the f-s
interface. The data are still the same as in the linear case. One can see a
meaningful difference with the case without duality.

(a) Normal stress at the f-s interface

(b) The Lagrange multipliers at the two f-s interfaces

Figure 16:  Burgers solution obtained with duality but no penalty term inside
the structure.

(a) solution

(b) Derivative versus $x$ of Burgers solution

Figure 17: Burgers solution obtained without duality but a penalty term inside
the structure ($L^{2}$).  Figure 18: Normal stress at the f-s interface (no
duality but a penalty term inside the structure.

(a) solution

(b) Derivative versus $x$ of Burgers solution

Figure 19: Burgers solution obtained with duality and a penalty term inside
the structure.

(a) Normal stress at the f-s interface

(b) Lagrange multipliers at the f-s interfaces

Figure 20: Normal stress and Lagrange multipliers at the f-s interface with
duality and a penalty term inside the structure.

(a) Reference solution

(b) Derivative versus $x$ of the reference solution

Figure 21: Reference solution of Burgers equation (duality + penalty term
inside the structure + refined mesh + time step)

(a) Normal stress at the f-s interfaces

(b) Lagrange multipliers at the f-s interfaces, (the sign of the Lagrange
multiplier is opposite to the one of the normal stress in this case)

Figure 22: Reference solution (Figure 21(a)) obtained with duality and a
penalty term inside the structure.

##  9 Conclusion

In this paper we have first studied several penalty methods in a simple fluid-
structure model assuming the case of a small reduce frequency in order to
partially decouple the inertia terms between the fluid and the structure. We
discussed the possibilities of several penalty terms inside the structure
($L^{2}$ or/and $H^{1}$ norms but also by penalizing the continuity of the
velocities at the boundary between the fluid and the structure. Our goal has
been to focus on the mathematical difficulties which can occur in the
numerical schemes due to the ill conditioning of the penalty models.We compare
the methods from the theoretical point of view and on a very simple 1D model.
We also introduced in this framework the Bertsekas penalty-duality algorithm
for the ensuring the velocity continuity at the fluid-structure interface. It
appears that the penalty-duality algorithm applied to this quasi-static model
without inertia terms is very efficient compared to other penalty methods.
Furthermore it leads to a much better condition number of the numerical scheme
as far as it not necessary to use a large value of the penalty coefficient.

In a second step we compared the various possibilities discussed in the first
step on two very simple 1D fluid-structure models. The first one is the
advection-diffusion equation with a prescribed movement of an immersed
structure and the second one is a similar model but with Burgers equation for
the fluid. In this case there a convection term (linear for the advection and
non linear for Burgers equation). The numerical tests show that the $L^{2}$
penalty term inside the structures is very efficient as far as the penalty
parameter is very small. The coupling with the penalty-duality algorithm
involving the velocity continuity at the interface enables to improve very
slightly the numerical tests as far as the inner penalty parameter is very
small. But it is much more efficient if this parameter is more moderate. This
can be therefore an interesting improvement in cases where the condition
numerber is an important point in the numerical scheme. One can forcast that
this is mainly the case for for 2D and more for 3D models.

Because this paper is mainly theoretically oriented, we focused on the
mathematical analysis of the penalty and penalty-duality methods for Stokes
equations coupled with a rigid structure in movement. The tests are just an
illustration in order to eillustrate our purpose and to point out the limits
of our conclusions. For a more physical study we refer to A. Falaise and E.
Liberge [19].

Nevertheless, an embedding in a global Euler-Lagrange representation is a
classical solution in order to be able to solve a fully time dependent model
where the structure (assumed to be rigid for sake of simplicity) is moving
inside the fluid with large displacements due to the forces applied on the
structure by the fluid. This aspect will be discussed in a forthcoming study
[15] where we use the so-called $\theta$-method [17] for transferring
informations between the fluid model and the one of the structure which can be
flexible and moving with large displacements.

## References

  * [1] P. Angot, C.H. Bruneau, and P. Fabrie, (1999), A penalization method to take into account obstacles in incompressible flows. Numer. Math. 81, 4, pp. 497-520.
  * [2] I. Babuska, and A. K. AzIz, (1972), Survey lectures on the mathematical foundations of the finite element method. The Mathematical Foundations of the Finite Element Method with Applications to Partial Differential Equations. A.K. AzIz ed., 5-359. Academic Press, New York, 1972.
  * [3] T. Belytschko, Y. Lu and L. Gu, (1994), Element-free Galerkin methods, International Journal for Numerical Methods in Engineering, vol. 37, 1994, pp. 229-256.
  * [4] M. Benamour , E. Liberge and C. Béghein, (2020), A volume penalization lattice Boltzmann method for simulating flows in the presence of obstacles. Journal of computational science, Elsevier, 2020, 39, pp.101050.
  * [5] M. Benamour, E. Liberge, C. Béghein, (2015), Lattice Boltzmann method for fluid flow around bodies using volume penalization. The international journal of multiphysics, Multi-Science Publishing, 2015, 9 (3), pp. 299-316.
  * [6] H. Ben Dhia, (1998), Problèmes mécaniques multi-échelles: la méthode Arlequin; Multiscale mechanical problems: the Arlequin method. Comptes Rendus de l’Académie des Sciences - Series IIB - Mechanics-Physics-Astronomy, Volume 326, Issue 12, 1998, pp. 899-904.
  * [7] H. Ben Dhia, (2005) Local-Global approaches. The Arlequin method. 5ème Colloque National de Calcul des Structures, Giens 2005. (archives ouvertes HAL: 00281019).
  * [8] M Bercovier, M Engelman, (1979), A finite element for the numerical solution of viscous incompressible flows. Journal of Computational Physics 30 (2), 181-201.
  * [9] D. Bertsekas (1975), Necessary and sufficient conditions for a penalty method to be exact. Mathematical programming, vol. $9$, no 1, pp.87-99., Springer-Verlag.
  * [10] R.Borker, D.Huang, S. Grimberg, C. Farhat, P. Avery and J. Rabinovitch, (2019), Mesh adaptation framework for embedded boundary methods for computational fluid dynamics and fluid-structure interaction, International Journal for Numerical Methods in Fluids, Vol. 90, no8, July 2019, pp. 389-424.
  * [11] J. Cea, (1968), Optimisation: Théorie et algorithmes. Dunod, Paris.
  * [12] M. Coquerelle and G.H. Cottet, (2008), A vortex level set method for the two-way coupling of an incompressible fluid with colliding rigid bodies. J. Comput. Phys. vol. 21, no 227, pp. 9121-9137.
  * [13] G.H. Cottet and E. Maitre, (2004), A Level-Set Formulation of Immersed Boundary Methods for Fluid-Structure Interaction Problems. Comptes Rendus de l’Académie des Sciences de Paris, 338, pp. 51-56.
  * [14] G. H. Cottet, F. Gallizio, A. Magni and I. Mortazavi, (2011), A Vortex Penalization Method for Flows with Moving Immersed Obstacles. In Proceedings of ASME-JSME-KSME Joint Fluids Engineering Conference, July, 2011, ASME, pp. 3703-3708.
  * [15] Ph. Destuynder and E. Liberge, (2021), Euler-Lagrange technics in fluid-structure modeling for the computation of forces. To appear.
  * [16] Ph. Destuynder Aeroelasticity (2020), An introduction to quasi-static aeroelasticity. Lectures notes, Springer-Verlag, Berlin-New-York.
  * [17] Ph. Destuynder, (2008) Aéroélasticitét Aéroacoustique, Hermès-Lavoisier, Paris-Londres.
  * [18] G.Duvaut et J.L. Lions, (1972), Inéquations variationnelles en mécanique et en physique. Dunod, Paris.
  * [19] A. Falaise and E. Liberge (2021), A penalty-duality strategy in a fluid structure problem. Internal research report LaSIE, La Rochelle University.
  * [20] P. Fedkiw, R., Aslam, T., Merriman, B., and Osher, S. A non-oscillatory eulerian approach to interfaces in multimaterial flows (the ghost fluid method). J. Comput. Phys. 152 (1999), 457-492.
  * [21] M. Fortin and R. Glowinski, (1982), Résolution numérique de problèmes aux limites par des méthodes de Lagrangien augmenté. Méthodes Mathématiques de l’informatique, no 9, Dunod, Paris.
  * [22] Y.C. Fung, (1969), An Introduction to the Theory of Aeroelasticity (Dover Books on Aeronautical Engineering).
  * [23] P. Grisvard (1992), Singularities in boundary value problems. Paris, Masson, 1992.
  * [24] J. Hovnanian, (2012), Méthode de frontières immergées pour la mécanique des fluides. Application à la simulation de la nage. https://tel.archives-ouvertes.fr/tel-00835013.
  * [25] D. Z. Huang, D. De Santis and C. Farhat, (2018), A family of position- and orientation-independent embedded boundary methods for viscous flow and fluid-structure interaction problems, Journal of Computational Physics, Volume 365, no15, July 2018, Pages 74-104.
  * [26] N. James, E. Maitre, I. Mortazavi. Immersed boundary methods for the numerical simulation of incompressible aerodynamic and fluid-structure interactions. Annales mathÃ©matiques Blaise Pascal, cedram, 2013, 20 (1), pp.139-173.
  * [27] E. Liberge et C. Béghein, (2019), Méthode de Lattice Boltzmann couplée avec la pénalisation volumique pour la simulation de problèmes d’interaction fluide structure. 24${}^{\hbox{\\`{e}me}}$ Congrès FranÃ§ais de Mécanique, Aug 2019, Brest, France.
  * [28] J. L. Lions, (1988), Contrôlabilité exacte perturbations et stabilisation de systèmes distribués. T.1, RMA n8, Masson, Paris.
  * [29] J.L. Lions and E. Magenes (1968) Problèmes aux limites non-homogènes et applications. vol. 1, Dunod, Paris. and (English version) Non-Homogeneous Boundary Value Problems and Applications Vol. 1, Springer, Berlin, (1972).
  * [30] J.L. Lions (1973), Perturbations singulièresdans les problèmes aux limites et en contròle optimal Springer-Verlag, Lectures notes in Mathematics, n0323, Springer-Verlag, Berlin-Heidelberg-New York.
  * [31] R. Mittal, and G. Iaccarino,(2005), Immersed boundary methods. Annu. Rev. Fluid. Mech (2005), 1-27.
  * [32] C. Peskin, (2002), The immersed boundary method. Acta Numer. 11 (2002), 1-39.
  * [33] P. Poncet, (2001), Méthodes particulaires pour la simulation des sillages tridimensionnels. PhD thesis, Université Joseph Fourier Grenoble, (2001).
  * [34] P.A. Raviart and V. Girault,(1986) Finite Element Methods for Navier-Stokes Equations. Springer Series in Computational Mathematics, Berlin-New-York.


# Nonstationary self-gravitating configurations of scalar and electromagnetic
fields

Ju V Tchemarina1 E G Alekseeva2 A N Tsirulev1 and N K Nuraliev1 1 Faculty of
Mathematics, Tver State University, 35 Sadovyi, Tver, Russia, 170002  2
Faculty of Robotics and Complex Automation, Bauman Moscow State Technical
University, 5/1 Baumanskaya 2-ya, Moscow, Russia, 105005
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Mathematical modeling of gravitating configurations of physical fields is one
of the priority directions of the modern theory of gravity. Most of the exact
solutions constructed within the framework of the general relativity are
static or stationary configurations. This is due to the objective complexity
of solving the Einstein equations under the assumption of nonstationarity. We
present an approach to constructing nonstationary configurations of a
spherically symmetric nonlinear real scalar field and the electromagnetic
field, which are assumed both to be minimally coupled to gravity. It is based
on the isolation of one invariant equation written in terms of the
characteristic function and scalar field potential. Using the proposed method,
an exact nonstationary solution with a nontrivial topology of space-time will
be constructed.

##  1 Introduction

Nowadays, the construction and study of nonstationary models of astrophysical
objects and objects of the microworld based on a self-gravitating scalar field
are becoming increasingly popular. One of the current trends in these studies
is the use of a scalar field to describe particle-like configurations with
nontrivial space-time topology, including spherically symmetric topological
geons [3, 5, 6, 8, 9, 11]. It is worth noting that almost all currently known
models of topological geons and wormholes are stationary and, accordingly, do
not take into account the dynamics of such objects. Also the problem of
stability is of great importance when studying configurations with scalar
fields. The stability properties of scalar wormholes under spherical
perturbations are the subject of research [13].

The purpose of this study is to demonstrate a method for constructing
nonstationary configurations of a spherically symmetric scalar field with a
charge. Using the proposed method, a new exact nonstationary solution will be
constructed. This solution is a wormhole whose throat size changes with time.

##  2 Action and basic equations

We start with the action for the gravitating system of nonlinear real scalar
field and electromagnetic field which assumed both to be minimally coupled to
gravity

$\Sigma=\int\\!\left(-\frac{1}{2}S+\varepsilon\langle
d\phi,d\phi\rangle-2V(\phi)-\frac{1}{2}{F}_{{i}{j}}{F}^{{i}{j}}\right)\sqrt[]{|g|}\,d^{\,4}x\,,$

where $F={F}_{{i}{j}}dx^{i}\wedge dx^{j}$ is the electromagnetic field tensor,
$S$ – the scalar curvature, $\phi$ – the scalar field, $V(\phi)$ – the self-
interaction potential, $\varepsilon=\pm 1$ – the sign of the scalar field
kinetic term, and the angle brackets denote the scalar product induced by a
given space-time metric $g$.

It is clear that the stress-energy tensor is completely specified by the
formulas

$T=T_{(\phi)}+T_{(em)}\,,$
$T_{{(\phi)}ij}\,=\,2\varepsilon\,\partial_{i}\phi\,\partial_{j}\phi\,-\,\left(\varepsilon\,g^{km}\,\partial_{k}\phi\,\partial_{m}\phi\,-\,2\,V\right)\,g_{ij}\,,$
$T_{{(em)}ij}=-2g_{ik}F_{jl}F^{kl}+\frac{1}{2}g_{ij}F_{kl}F^{kl}.$

The geometry of space-time is governing by the system of Einstein-Klein-
Gordon-Maxwell equations

$\mathcal{R}_{ij}\,-\frac{1}{2}\,S\,g_{ij}\,=\,T_{ij}\,,$ (1)
$\frac{1}{\sqrt{|g|}}\,\partial_{i}\,(\sqrt{|g|}\,g^{ij}\,\partial_{j}\phi)\,+\,\varepsilon\,V^{\prime}_{\phi}\,=\,0\,,\quad|g|\,=\,|det\,(g_{ij})|\,,$
(2) $\frac{1}{\sqrt{|g|}}\frac{\partial}{\partial
x^{i}}(\sqrt{|g|}F^{ij})=0\,.$ (3)

The spherical symmetry of space-time allows us to write the electromagnetic
field tensor in the following form

$F=F_{tr}dt\wedge dr,\quad\text{where}\quad F_{tr}=F_{tr}(t,r)\,.$

Solving equation (3) in the coordinate basis $(t,r,\theta,\varphi)$, in which
the metric takes the form

$g=A^{2}dt\otimes dt-B^{2}dr\otimes dr-C^{2}(d\theta\otimes
d\theta+\sin^{2}\theta\,d\varphi\otimes d\varphi)\,,$ (4)

we obtain

$F=\frac{ABq}{C^{2}}dt\wedge dr\,.$

To simplify calculations, we will make the transition to the orthonormal basis
of vector fields and dual basis of one-forms

$A\,dt=e^{0},\quad B\,dr=e^{1},\quad C\,d\theta=e^{3},\quad
C\sin\theta\,d\varphi=e^{4},\,\quad F=\frac{q}{C^{2}}e^{0}\wedge e^{1}\,.$

In orthonormal basis, we may write the system of Einstein equations

$-2\frac{C_{(1)(1)}}{C}+2\frac{B_{(0)}C_{(0)}}{BC}-\frac{C_{(1)}^{2}-C_{(0)}^{2}-1}{C^{2}}\,=\,\varepsilon\left(\phi_{(1)}^{2}+\phi_{(0)}^{2}\right)+2\,V+\frac{q^{2}}{C^{4}},$
(5)
$-2\frac{C_{(0)(0)}}{C}+2\frac{A_{(1)}C_{(1)}}{AC}+\frac{C_{(1)}^{2}-C_{(0)}^{2}-1}{C^{2}}\,=\,\varepsilon\left(\phi_{(1)}^{2}+\phi_{(0)}^{2}\right)-2\,V-\frac{q^{2}}{C^{4}},$
(6)
$\frac{A_{(1)(1)}}{A}-\frac{B_{(0)(0)}}{B}+\frac{C_{(1)(1)}}{C}-\frac{C_{(0)(0)}}{C}+\frac{A_{(1)}C_{(1)}}{AC}-\frac{B_{(0)}C_{(0)}}{BC}=\,\varepsilon\left(\phi_{(0)}^{2}-\phi_{(1)}^{2}\right)-2\,V+\frac{q^{2}}{C^{4}},$
(7)
$-2\frac{C_{(0)(1)}}{C}+2\frac{B_{(0)}C_{(1)}}{BC}\,\equiv\,-2\frac{C_{(1)(0)}}{C}+2\frac{A_{(1)}C_{(0)}}{AC}=\,2\,\varepsilon\phi_{(0)}\phi_{(1)}\,,\,\,\,$
(8)

and Klein-Gordon equation

$\phi_{(0)(0)}-\phi_{(1)(1)}+\phi_{(0)}\frac{(BC^{2})_{(0)}}{BC^{2}}-\phi_{(1)}\frac{(AC^{2})_{(1)}}{AC^{2}}+\,\varepsilon
V^{\prime}_{\phi}\,=\,0.$

Lower indexes in round parentheses denote the directional derivatives along
the corresponding vector field.

##  3 Characteristic function

Following [7, 9, 10, 11, 12], in an invariant manner, we will define the
function

$f=f(\phi,C)\,=\,-\,\left\langle
dC,dC\right\rangle\,=\,C_{(1)}^{2}\,-\,C_{(0)}^{2}\,,\,\,dC=\,C_{(0)}{e^{0}}\,-\,C_{(1)}{e^{1}}\,.$
(9)

The solutions of equation $f(\phi,C)=0$ define hypersurfaces on which the
1-form $dC$ becomes null. For stationary configurations it is true on event
horizons and hence the function $f(C,\phi)$ will be referred to as the
characteristic function. The asymptotic behavior of this function also makes
it possible to assign a solution to a certain type.

We call the scalar field configuration as nonstationary or stationary
according as $\phi\neq\phi(C)$ or $\phi=\phi(C)$. The latter case is well
studied. Therefore, we will focus on the consideration of non-stationary
configurations. Note that the radius $C$ is defined in an invariant manner by
the metric (4).

Making use of the formula (9), we find the derivatives of $f(\phi,C)$ in the
directions of basis vector fields

$f_{(1)}=2C_{(1)}C_{(1)(1)}-2C_{(0)}C_{(0)(1)}\,,\,\,f_{(0)}=2C_{(1)}C_{(1)(0)}-2C_{(0)}C_{(0)(0)}\,.$

We thus obtain the relation

$C_{(1)(1)}=\frac{f_{(1)}}{2C_{(1)}}+\frac{C_{(0)}C_{(0)(1)}}{C_{(1)}}\,.$

Substituting the last expression in equation (5), we obtain

$-\frac{2}{C}\left(\frac{f_{(1)}}{2C_{(1)}}+\frac{C_{(0)}C_{(0)(1)}}{C_{(1)}}\right)+\frac{2B_{(0)}C_{(0)}}{BC}-\frac{f-1}{C^{2}}=\varepsilon\left(\phi^{2}_{(1)}+\phi^{2}_{(0)}\right)+2V+\frac{q^{2}}{C^{4}}\,.$

After multiplication by $C_{(1)}\,C^{2}$, it gives

$[C(f-1)]_{(1)}=2\frac{B_{(0)}C_{(0)}CC_{(1)}}{B}-2CC_{(0)}C_{(0)(1)}-C_{(1)}\left(C^{2}\varepsilon\left(\phi^{2}_{(1)}+\phi^{2}_{(0)}\right)+2C^{2}V+\frac{q^{2}}{C^{2}}\right)\,.$
(10)

From equation (8), we find

$C_{(0)(1)}=\frac{B_{(0)}C_{(1)}}{B}-\varepsilon C\phi_{(0)}\phi_{(1)}\,.$
(11)

By combining equations (10) and (11), we obtain equation

$[C(f-1)]_{(1)}=C_{(1)}[\varepsilon
C^{2}(\phi^{2}_{(1)}-\phi^{2}_{(0)})-2C^{2}V-\frac{q^{2}}{C^{2}}]+\phi_{(1)}[2C^{2}\varepsilon(\phi_{(0)}C_{(0)}-C_{(1)}\phi_{(1)})]\,.$

By a similar sequence of conversion, we obtain from equation (6) the
appropriate equation

$[C(f-1)]_{(0)}=C_{(0)}[\varepsilon
C^{2}(\phi^{2}_{(1)}-\phi^{2}_{(0)})-2C^{2}V-\frac{q^{2}}{C^{2}}]+\phi_{(0)}[2C^{2}\varepsilon(\phi_{(0)}C_{(0)}-C_{(1)}\phi_{(1)})]\,.$

For a non-stationary configuration the scalar field $\phi$ and radius $C$ can
be considered as independent variables. Therefore, the last two equations are
equivalent to one invariant equation

$d[C(f-1)]=C^{2}(\varepsilon\left(\phi_{(1)}^{2}-\phi_{(0)}^{2}\right)-2V-\frac{q^{2}}{C^{4}})dC+2C^{2}\varepsilon(\phi_{(0)}C_{(0)}-\phi_{(1)}C_{(1)})d\phi\,,$
(12)

where

$\phi_{(0)}C_{(0)}-\phi_{(1)}C_{(1)}=<d\phi,dC>\,,\quad\phi_{(1)}^{2}-\phi_{(0)}^{2}=-<d\phi,d\phi>\,.$

According to equation (12),

$<d\phi,dC>=\frac{\varepsilon}{2C}\,f^{\prime}_{\phi}\,,\quad<d\phi,d\phi>=-\varepsilon\left(2V+\frac{f^{\prime}_{C}}{C}\,+\frac{f-1}{C^{2}}+\frac{q^{2}}{C^{4}}\right)\,.$
(13)

Thus, we isolate one invariant equation written in terms of the characteristic
function and scalar field potential.

##  4 Approach to constructing nonstationary configurations

A significant part of the proposed method is the use of coordinates
$(\phi,C,\theta,\varphi)$ which is possible, at least locally, for any
nonstationary configuration. Note that the expressions (9) and (13) define
metric coefficients

$g^{\phi\phi}=<d\phi,d\phi>,\quad g^{\phi
C}=<dC,d\phi>,\quad\,g^{CC}=<dC,dC>\,.$

Thus, finding the inverse matrix, we obtain the components of the covariant
metric tensor and write the metric in coordinates $(\phi,C,\theta,\varphi)$

$ds^{2}=-4\frac{\varepsilon C^{4}fd\phi^{2}+C^{3}f^{\prime}_{\phi}d\phi
dC+(2C^{4}V+C^{3}f^{\prime}_{C}+C^{2}(f-1)+q^{2})dC^{2}}{4fC^{2}(2C^{2}V+Cf^{\prime}_{C}+f-1)-\varepsilon(Cf^{\prime}_{\phi})^{2}+4fq^{2}}-$
$-C^{2}\left(d\theta^{2}+\sin^{2}\theta\,d\varphi^{2}\right)\,.$ (14)

Next we find the electromagnetic field tensor in the coordinates
$(\phi,C,\theta,\varphi)$ in the form

$F=\frac{2q\,d\phi\wedge dC}{\sqrt{-4\varepsilon
fC^{2}(Cf^{\prime}_{C}+f-1+2C^{2}V)-4\varepsilon
q^{2}f+C^{2}{f^{\prime}_{\phi}}^{2}}}\,.$

Then we sequentially find the components of the stress-energy tensor and write
down the system of Einstein equations (1). We do not give the corresponding
expressions here because of their bulkiness. It turns out that the resulting
system of equations is equivalent to the single equation (2)

$8V^{2}C^{8}+(6C^{7}f^{\prime}_{C}+8C^{6}f-2\varepsilon
C^{6}f^{\prime\prime}_{\phi\phi}+8C^{4}q^{2}-8C^{6})V+\varepsilon
C^{6}V^{\prime}_{\phi}f^{\prime}_{\phi}-f^{\prime\prime}_{CC}C^{6}f+$
$+\varepsilon C^{5}f^{\prime\prime}_{\phi
C}f^{\prime}_{\phi}+\varepsilon(C^{4}-fC^{4}-C^{2}q^{2}-f^{\prime}_{C}C^{5})f^{\prime\prime}_{\phi\phi}+(f^{\prime}_{C})^{2}C^{6}+(3C^{3}q^{2}+3C^{5}f-3C^{5})f^{\prime}_{C}+$
$+4f^{2}C^{4}+(8C^{2}q^{2}-6C^{4})f+2C^{4}+2q^{4}-4C^{2}q^{2}=0.$ (15)

Thus, the characteristic function and the scalar field potential are connected
by one single Klein-Gordon equation. The solution can be completed by
integration of equation (15). However, this equation is too complicated to
solve it explicitly, without imposing any additional restrictions. It seems
most convenient to find the unknown function ($f(\phi,C)$ or $V(\phi)$) as a
series or numerical solution.

##  5 Special case for metric function

Despite the complexity, equation (15) in some cases allows to obtain exact
analytical solutions. Choosing the characteristic function in special form, we
obtain an explicit solution for the scalar field potential. Now suppose that

$f\left(\phi,C\right)=1+C^{2}h(\phi),\quad\varepsilon=-1\,,\quad q=0.$ (16)

This choice of the characteristic function is not accidental. It is due to the
possibility of obtaining exact solutions with a nontrivial topology of space-
time. The use of a phantom scalar field with a negative kinetic term is
necessary for constructing solutions with nontrivial topology. The
impossibility of the existence of topological geons and wormholes without a
horizon with a positive kinetic term was proved in [1, 2, 4].

The substitution (16) reduces equation (15) to the form

$8\,V^{2}+\left(2\,h^{\prime\prime}_{\phi\phi}+20h\right)V-V^{\prime}_{\phi}\,h^{\prime}_{\phi}-2\,{h^{\prime}_{\phi}}^{2}+2\,h^{\prime\prime}_{\phi\phi}\,h+12\,{h}^{2}\,=0\,.$
(17)

The formal solution of equation (17) can be readily written down in terms of
integrals

$V(\phi)=-\frac{3\,h(\phi)}{2}-\frac{{h^{\prime}_{\phi}}^{2}}{8\,h(\phi)}-\frac{e^{\emph{F}}\,{h^{\prime}_{\phi}}^{2}}{8\,h^{2}(\phi)\,\int{\frac{e^{\emph{F}}\,h^{\prime}_{\phi}}{h^{2}(\phi)}d\,\phi}}\,,\,\,\emph{F}(\phi)=-4\int{\frac{h(\phi)}{h^{\prime}_{\phi}}d\,\phi}\,.$
(18)

These integral formulas allow us to recover the scalar field potential by a
given characteristic function.

##  6 An exact nonstationary solution with a nontrivial topology of space-time

Using the proposed method, we construct a model of nonstationary wormhole.
Substituting in the formulas (18) the function $h(\phi)=\,\phi^{2}\,-\,1$ , we
find,

$V(\phi)=\,1\,-\frac{3\,\phi^{2}}{2}\,+\frac{1\,+\,\phi^{2}\,-\,e^{\,\phi^{2}}}{2\left(1\,+\,e^{\,\phi^{2}}\left(\phi^{2}\,-\,1\right)\right)}\,.$

The integration constant was chosen so that the scalar field potential on the
wormhole throat $\left(\phi=0\right)$ is a regular function.

The metric (14) now takes the form

$ds^{2}=\frac{1}{\Delta}\left(\left(C^{2}\left(1\,-\,\phi^{2}\right)-\,1\right)d\phi^{2}\,+\,2\,\phi\,C\,d\phi\,dC\,+\,\frac{\phi^{2}\,\left(1\,-\,e^{\,\phi^{2}}\right)}{1\,+\,e^{\,\phi^{2}}\left(\phi^{2}\,-\,1\right)}dC^{2}\right)-$
$-C^{2}\left(d\theta^{2}+\sin^{2}\theta\,d\varphi^{2}\right)\,,$

where

$\Delta=\frac{\phi^{2}\,\left(e^{\,\phi^{2}}\,-\,\phi^{2}\,C^{2}\,-\,1\right)}{1\,+\,e^{\,\phi^{2}}\left(\phi^{2}\,-\,1\right)}\,.$

Space-time, in the general relativity, is considered as a manifold with a
Lorentzian signature. Therefore, the solution is defined in the area which is
determined by the inequality

$e^{\,\phi^{2}}\,-\,\phi^{2}\,C^{2}\,-\,1\,<\,0\,.$

Next, we revert to the usual coordinates $(t,r,\theta,\varphi)$.

By letting

$\phi\,=\,r\,,\quad
C\,=\,\cosh(t)\sqrt{\frac{e^{\,r^{2}}\,-\,1}{r^{2}}}\,=\,\cosh(t)\left(1\,+\,\frac{r^{2}}{4}\,+\,O(r^{4})\right)\,,\quad
r\rightarrow 0\,,$

the metric reduces to the standard form (4).

Therefore,

$A^{2}=\frac{e^{\,r^{2}}\,-\,1}{r^{2}}=1+\frac{r^{2}}{2}+O(r^{4})\,,\quad
B^{2}=\frac{e^{\,r^{2}}\left(r^{2}\,-\,1\right)\,+\,1}{r^{2}\left(e^{\,r^{2}}\,-\,1\right)}=\frac{1}{2}+\frac{r^{2}}{12}+O(r^{4})\,,\quad
r\rightarrow 0\,.$

The plots of the metric functions $A(r),\,B(r)$ and the surface $C(t,r)$ are
presented in Figure 1.

Figure 1: Left plot: metric functions $A(r)$ and $B(r)$. Right plot: surface
$C(t,r)$.

Metric functions are regular, positive and even with respect to $r\,$. The
coordinate $r$ is space-like everywhere. This allows the solution to be
interpreted as a wormhole, the size of which $C(t,0)\,=\,\cosh(t)$ varies over
time, taking the smallest value $C\,=\,1$ at $t\,=\,0\,$. The constructed
model of a spherically symmetric wormhole illustrates the use of new
mathematical methods for studying and solving the Einstein’s equations.

##  7 Conclusions

We present an approach to constructing nonstationary configurations of a
spherically symmetric nonlinear real self-gravitating scalar field and the
electromagnetic field. It is based on the isolation of one invariant equation
written in terms of the characteristic function and scalar field potential.
Despite the fact that the question of a general solution of this equation is
open, this approach allows us to test known numerical nonstationary solutions.
Under additional conditions, it makes it possible to construct nonstationary
solutions with specified properties.

## References

## References

  * [1] Morris M S and Thorne K S 1988 Am. J. Phys. 56, 395
  * [2] Hochberg D and Visser M 1997 Phys. Rev. D 56, 4745 (arXiv: gr-qc/9704082)
  * [3] Dowker H F and Sorkin R D 1998 Class. Quantum Grav. 15 1153 (arXiv: gr-qc/9609064)
  * [4] Bronnikov K A and Shikin G N 2002 Grav. Cosmol. 8 107 (arXiv: gr-qc/0109027)
  * [5] Bronnikov K A and Fabris J C 2006 Phys. Rev. Lett. 96, 251101 (arXiv: gr-qc/0511109)
  * [6] Bronnikov K A, Dehnen H and Melnikov V N 2007 Gen. Rel. Grav. 39 973 (arXiv: gr-qc/0611022)
  * [7] Tchemarina Ju V and Tsirulev A N 2009 Grav. Cosmol. 15 94
  * [8] Bolokhov S V, Bronnikov K A and Skvortsova M V 2012 Class. Quantum Grav. 29 245006
  * [9] Solovyev D A and Tsirulev A N 2012 Class. Quantum Grav. 29 055013
  * [10] Nikonov V V, Potashov I M and Tsirulev A N 2016 Math. Model. Geom. 4 No 2 http://mmg.tversu.ru
  * [11] Kratovitch P V, Potashov I M, Tchemarina Ju V and Tsirulev A N 2017 J. Phys.: Conf. Ser. 934 012047 (arXiv: gr-qc/1805.04447)
  * [12] Kratovich P V and Tchemarina Ju V 2017 Math. Model. Geom. 5 No 2 http://mmg.tversu.ru (arXiv: gr-qc/1805.02698)
  * [13] Bronnikov K A 2018 Particles 1, Issue 1, 5 (arXiv: 1802.00098)