# A numerical method for wave-structure interactions in the Boussinesq regime

###### Abstract.

The goal of this work is to study waves interacting with partially immersed
objects allowed to move freely in the vertical direction, and in a regime in
which the propagation of the waves is described by the one dimensional
Boussinesq-Abbott system. The problem can be reduced to a transmission problem
for this Boussinesq system, in which the transmission conditions between the
components of the domain at the left and at the right of the object are
determined through the resolution of coupled forced ODEs in time satisfied by
the vertical displacement of the object and the average discharge in the
portion of the fluid located under the object. We propose a new extended
formulation in which these ODEs are complemented by two other forced ODEs
satisfied by the trace of the surface elevation at the contact points. The
interest of this new extended formulation is that the forcing terms are easy
to compute numerically and that the surface elevation at the contact points is
furnished for free. Based on this formulation, we propose a second order
scheme that involves a generalization of the MacCormack scheme with nonlocal
flux and a source term, which is coupled to a second order Heun scheme for the
ODEs. In order to validate this scheme, several explicit solutions for this
wave-structure interaction problem are derived and can serve as benchmark for
future codes. As a byproduct, our method provides a second order scheme for
the generation of waves at the entrance of the numerical domain for the
Boussinesq-Abbott system.

Geoffrey Beck, David Lannes, Lisl Weynans

Mathematics Subject Classification: 35G61, 35Q35, 74F10, 65M08

Keywords: Wave-structure interactions, Initial boundary value problems,
Boussinesq system, Numerical analysis

##  1\. Introduction

###  1.1. Presentation of the problem

While the first studies of the interactions of waves with floating structures
go back at least to John’s paper [20], or to the phenomenological integro-
differential equation derived by Cummins to describe the linear motion of
floating structures [11], this research field became increasingly active in
recent years. A first reason for this renewed interest is related to the
development of renewable marine energies as one of the tools for energy
transition. Indeed, several devices of offshore wind-turbines and wave-energy
convertors involve partially immersed structures [3].

A second reason for the recent mathematical activity on wave-structure
interactions is that this has been made technically feasible thanks to the
recent progresses on the mathematical understanding of the propagation of
water waves. The initial value problem in domains without boundaries
(${\mathbb{R}}^{d}$ or ${\mathbb{T}}^{d}$) is now well understood for the full
water waves (also called free-surface Euler) equations, as well as for
asymptotic models in shallow water (such as the nonlinear shallow water
equations, the Boussinesq systems, the Serre-Green-Naghdi equations).
Recently, the initial value problem has also been studied in domains with a
boundary. When the fluid domain is delimited by vertical sidewalls, the water
waves equations have been studied in [2]; in the case of non-vertical
sidewalls, this problem has been considered in [34, 29] for the water waves
equations, and in [25] for the shallow water and Green-Naghdi equations. The
initial boundary value problem, in which one imposes initial and boundary
datas, has also been investigated for the Boussinesq equations [10, 26]. These
advances make it more realistic to address the issues raised by the presence
of a partially immersed object.

From the numerical point of view, efficient numerical codes based on shallow
water models have been developed recently and can be used to address realistic
submersion issues (see for instance [35, 14]); here also, it is now a
reasonable prospect to address the specific difficulties raised by wave-
structure interactions.

The present paper is a contribution to the theoretical and numerical
understanding of these interactions, inasmuch as it provides a precise
description of the motion of a partially immersed object allowed to move
freely in the vertical direction under the action of waves described by a
nonlinear dispersive model (the standard Boussinesq-Abbott system), see Figure
1. It has to be considered as a partial (affirmative) answer to the wider
question: can the efficient modelling of waves based on shallow water models
be extended to allow the presence of floating structures? If this happens to
be true, the gain in computational time would allow to investigate the
behavior of many floating structures (the so-called farms of wave-energy
convertors or offshore wind turbines), as well as their impact on the wave
fields, which can have significant consequences in coastal regions. Answering
such question is out of reach for CFD methods that can be used to describe the
behavior of one wave-energy convertor, and also, to a lower extent, for
potential methods (see for instance [13, 17]). On the other hand, the linear
methods based on Cummins’ equation used in commercial softwares such as Wamit
neglect the nonlinear effects that can be important [32], especially in
shallow water, and are unable to provide a precise description of the impact
of a wave-farm on the wave-field.

Figure 1. The floating object

The presence of a floating structure in a shallow water model can be taken
into account following the approach proposed in [24] where the horizontal
plane is decomposed into two regions: the interior region (below the floating
object), and the exterior region (below the free surface waves). In the
exterior regions, the standard (depth integrated) shallow water model is used,
while in the interior region, an additional pressure term is present. This
pressure term corresponds to the pressure exerted by the fluid on the object
(and which eventually makes it move through Newton’s equations), and can be
understood as the Lagrange multiplier associated with the constraint that,
under the object, the surface elevation of the waves is constrained as it must
by definition coincide with the bottom of the object. It is possible to relax
this constraint by approximating the pressure term by a pseudo-compressible
relaxation; one can then use the same kind of asymptotic preserving schemes as
for the low-Mach limit in compressible gases. This approach has been used in
the present context in [15, 16], and is also relevant for other instances of
partially congested flows [33, 12, 5]. In this paper, we rather consider the
original (non relaxed) problem, which requires to understand precisely the
coupling between the interior and exterior regions.

It turns out that this wave-structure interaction problem can be reduced to an
initial boundary value problem for the wave model in the exterior region, with
non standard boundary (or transmission) conditions. In the case where the
horizontal dimension is $d=1$, the object has vertical walls located at
$x=\pm\ell$ as in Figure 1 and is only allowed to move vertically, and if the
wave model is given by the nonlinear shallow water equations, it was shown in
[24] that this transmission problem takes the form (in dimensionless
variables, see Section 2 for details),

$\begin{cases}\partial_{t}\zeta+ \partial_{x}q=0,\\\
\partial_{t}q+\varepsilon\partial_{x}(\frac{1}{h}q^{2})+h\partial_{x}\zeta=0,\end{cases}\quad\mbox{
on }(-\infty,-\ell)\cup(\ell,+\infty),$

where $\zeta$ is the elevation of the surface and $q$ the horizontal
discharge. Using the notation $\llbracket q\rrbracket=q(\ell)-q(-\ell)$ and
$\langle q\rangle=\frac{1}{2}(q(-\ell)+q(\ell))$, the transmission conditions
are given by

$\llbracket q\rrbracket=-2\ell\dot{\delta}\quad\mbox{ and }\quad\langle
q\rangle=\langle q_{\rm i}\rangle$

where the function $\delta$ and $\langle q_{\rm i}\rangle$ (representing
respectively the vertical displacement of the object and the mean discharge
under the object) solve an ODE of the form

$\frac{d}{dt}\theta=F(\theta,\zeta_{|_{x=-\ell}},\zeta_{|_{x=\ell}}),$

with $\theta=(\delta,\dot{\delta},\langle q_{\rm i}\rangle)$ and $F$ a smooth
function of no importance at this stage of the discussion. The coupling acts
in two ways: it is necessary to known $\theta$ to solve the transmission
problem for $(\zeta,q)$, and it is necessary to know the solution $(\zeta,q)$
to determine the forcing term
$F(\theta,\zeta_{|_{x=-\ell}},\zeta_{|_{x=\ell}})$ in the ODE for $\theta$. A
key point in the mathematical analysis of this problem is the regularity of
the traces $\zeta_{|_{x=\pm\ell}}$. Such a control is furnished by the
construction of a Kreiss symmetrizer, as shown in [19] where general initial
boundary value problems, possibly with a free boundary, are considered for a
wide class of hyperbolic systems; they include the above transmission problem
as well as the more complex free boundary problem one has to deal with when
the lateral boundaries of the object are not vertical.

Numerically, the evaluation of the traces $\zeta_{|_{x=\pm\ell}}$ also
requires a careful treatment which relies on the Riemann invariants associated
with the nonlinear shallow water equations [24]; we also refer to [6] for a
higher order scheme, to [9] where a wave-energy device is simulated using this
approach (the oscillating water column), and to [36] where controlability
issues were also addressed. The more complex case of an object freely floating
and with non-vertical walls (and therefore nontrivial dynamics for the contact
points) has been solved theoretically in [19], and numerically in [18] using
ALE methods to treat the evolution of the contact points.

A variant of the above wave-structure interaction problem for the viscous
nonlinear shallow water equations was also considered in [28] and an extension
to the case of horizontal dimension $d=2$ with radial symmetry has been
considered theoretically in [7] and the so-called decay test (or return to
equilibrium) investigated in the same configuration under an additional
assumption of linearity in [8]. Let us also mention [31] where the dynamics of
trapped air pockets are studied.

We propose here an extension in another direction. The principal drawback of
the nonlinear shallow water equations is that they neglect the dispersive
effects that play an important role is some important situations (they allow
for instance the existence of solitary waves). The most simple models that
generalize the nonlinear shallow water equations by adding dispersive terms
are the Boussinesq equations (see [22] for a recent review on shallow water
models). Replacing the nonlinear shallow water equations by the so called
Boussinesq-Abbott system in the above example, one obtains the following
transmission problem (see Section 2 for more details),

$\begin{cases}\partial_{t}\zeta+ \partial_{x}q=0,\\\
(1-\kappa^{2}\partial_{x}^{2})\partial_{t}q+\varepsilon\partial_{x}(\frac{1}{h}q^{2})+h\partial_{x}\zeta=0,\end{cases}\quad\mbox{
on }(-\infty,-\ell)\cup(\ell,+\infty),$

and with transmission conditions

$\llbracket q\rrbracket=-2\ell\dot{\delta}\quad\mbox{ and }\quad\langle
q\rangle=\langle q_{\rm i}\rangle$

where the function $\delta$ and $\langle q_{\rm i}\rangle$ (representing
respectively the vertical displacement of the object and the mean discharge
under the object) solve an ODE of the form

$\frac{d}{dt}\theta=F_{\kappa}(\theta,\zeta_{|_{x=-\ell}},\zeta_{|_{x=\ell}},\kappa^{2}(\partial_{t}^{2}\zeta)_{|_{x=-\ell}},\kappa^{2}(\partial_{t}^{2}\zeta)_{|_{x=\ell}}),$

with $\theta=(\langle q_{\rm i}\rangle,\dot{\delta},\delta)^{\rm T}$ and
$F_{\kappa}$ a smooth function of no importance at this stage of the
discussion.

The differences with the nondispersive case considered above are the operator
$(1-\kappa^{2}\partial_{x}^{2})$ applied in front of $\partial_{t}q$ in the
evolution equations and a contribution of the trace of $\partial_{t}^{2}\zeta$
to the forcing term in the ODE for $\theta$. Formally, the nondispersive case
is obtained by setting $\kappa=0$, but the mathematical and numerical
differences are considerable. Contrary to the hyperbolic case mentioned above,
there is no general theory for initial boundary value problems associated with
nonlinear dispersive systems and this is why several approximations have been
used to bypass this issue. In [6], wave-structure interactions using a
Boussinesq model was used, but the issue at the boundary was avoided by using
the (dispersionless) nonlinear shallow water equations in a small region
around the object; in [30] the behavior at the boundary was approximated at
second order using Bessel expansions and matched asymptotics; in [21],
Boussinesq type equations where computed in the whole domain, neglecting the
singularities of the surface elevation and of the discharge at the contact
line, while the presence of the object is taken into account by adding an
additional pressure term in the interior region. There are also approximate
methods based on sponge layers and artificial source terms which are often
used to generate waves at the entrance of the numerical domain [37]. Such
methods are far too rough to be used in the present case, where a precise
description of the waves at the contact points is needed; indeed, as shown in
[10], the behavior at the contact points can be quite complex and exhibit
dispersive boundary layers. This is why a new method to handle non-homogeneous
initial boundary value problems for the Boussinesq equations was proposed in
[26] and numerically implemented with an order $1$ scheme.

The approach used in the present paper allows us to treat the issues related
to the initial boundary value problem for Boussinesq-type equations without
any approximation; a byproduct of independent interest of the present paper is
that it furnishes a second order method for the generation of waves at the
numerical boundary of the fluid domain for the Boussinesq equations, hereby
complementing the first order generation scheme of [26].

As for the hyperbolic case discussed previously, the control of the traces of
$\zeta$ (and a fortiori of $\partial_{t}^{2}\zeta$) is a key ingredient of the
analysis, both from the PDE and numerical perspectives; however, due to the
presence of dispersion, there are no such things as a Kreiss symmetrizer or
Riemann invariants to help us. To solve these issues, we propose in this paper
an extended formulation of the equations. After remarking that the traces
$\zeta_{|_{x=\pm\ell}}$ solve a second order forced ODE, we introduce new
unknowns $\underline{\zeta}_{\pm}$ defined as the solutions of this ODE. This
allows us to replace $\zeta_{|_{x=\pm\ell}}$ by $\underline{\zeta}_{\pm}$ in
the forcing term $F_{\kappa}$ in the ODE for $\theta$, hereby avoiding the
computation of the traces. The resulting extended formulation just consists in
replacing the above ODE for $\theta$ by the higher dimensional ODE

$\frac{d}{dt}\Theta={\mathcal{G}}(\Theta,(R_{1}{\mathfrak{f}}_{\rm
sw})_{|_{x=\pm\ell}})$

with $\Theta=(\langle q_{\rm
i}\rangle,\dot{\delta},\dot{\underline{\zeta}_{+}},\dot{\underline{\zeta}_{-}},\delta,\underline{\zeta}_{+},\underline{\zeta}_{-})^{\rm
T}$ and ${\mathcal{G}}$ is a smooth mapping. This ODE is forced by the terms
$(R_{1}{\mathfrak{f}}_{\rm sw})_{|_{x=\pm\ell}}$ which depend on the solution
$(\zeta,q)$ of the Boussinesq equations in the fluid domain; the precise
meaning of this term will be given in Section 2, the important thing being
that the control of the trace of the quantities $R_{1}{\mathfrak{f}}_{\rm sw}$
at $x=\pm\ell$ does not raise any theoretical nor numerical difficulty.

The second step of our approach consists in transforming this extended
transmission problem into an initial value problem coupled with forced ODEs:
this means that we do no longer have to bother with the boundary conditions
(which are automatically propagated by the flow). This new formulation can be
written as a system of conservation laws with nonlocal flux and an
exponentially localized source term,

$\partial_{t}U+\partial_{x}({\mathfrak{F}}_{\kappa}(U))={\mathcal{S}}_{\pm}(\Theta,(R_{1}{\mathfrak{f}}_{\rm
sw})_{|_{x=\pm\ell}}){\mathfrak{b}}(x\mp\ell)\quad\mbox{in
}\quad\pm(\ell,\infty),$

where ${\mathfrak{F}}_{\kappa}$ denotes the nonlocal flux, while
${\mathcal{S}}_{\pm}$ is a smooth function of its arguments and
${\mathfrak{b}}$ an exponentially localized function. The above ODE for
$\Theta$ allows one to compute the source term in this system of nonlocal
conservation laws; conversely, the resolution of this system allows one to
compute the forcing terms $(R_{1}{\mathfrak{f}}_{\rm sw})_{|_{x=\pm\ell}}$ in
the ODE for $\Theta$: the coupling acts therefore both ways.

One of the advantages of this new formulation is that we can implement on it a
second order scheme that couples a MacCormack predictor corrector scheme
(generalized to handle nonlocal fluxes and a source term) for the computation
of the waves, and a second order Heun scheme for the computation of the forced
ODEs. We also exhibit several exact explicit solutions that we use to study
the convergence of our code and that are of independent interest.

###  1.2. Organization of the paper

In Section 2, we derive the formulation of the problem our numerical scheme is
based on: we first recall in §2.1 the reduction of [4] to a transmission
problem for the Boussinesq equation on the two connected components of the
exterior region, and then show in §2.2 that the traces of the surface
elevation at the contact points satisfy a forced second order ODE that we use
to write the new augmented formulation of the transmission problem in §2.3;
this transmission problem is finally rewritten as an initial boundary value
problem in §2.4.

The numerical schemes are presented in Section 3. The initial value problem
obtained in the previous section is a set of two conservation equations with
nonlocal flux and an exponentially decaying source term whose coefficient is
found by solving a set of forced second order ODEs. We propose two numerical
schemes based on an abstract formulation of these equations. The first one,
described in §3.2, is of first order and is an adaptation of the Lax-
Friedrichs scheme to the present context. The second one, studied in §3.3, is
of second order. It is based on the MacCormack predictor-corrector scheme for
the two conservation PDEs (with adaptations to handle the nonlocal flux and
the source term), and on a Heun scheme for the ODE part.

Numerical simulation are then presented in Section 4. We investigate several
configurations exploring different aspects of the coupling between the
Boussinesq equations and the forced ODEs used in the transmission conditions.
Wave generation is considered in §4.1, and is of independent interest as it
provides a way to generate waves at the entrance of the numerical domain for
the Boussinesq equations. The return to equilibrium test in which an object
oscillates vertically after being released from an out of equilibrium position
is studied in §4.2; in the linear case, an explicit solution is exhibited and
computed via Laplace transforms, and this solution is used to assess the
precision of our scheme. Interactions of waves with a fixed object are then
investigated in §4.3; here also, an exact solution is derived in the linear
case and used for validation. The most general configuration of waves
interacting with an object allowed to move freely in the vertical direction is
then considered in §4.4.

Throughout this article, we work with an abstract and concise formulation of
the equations. The precise equations, with the expressions of the various
coefficients involved, is postponed to Appendix A.

###  1.3. Notation

\- The horizontal axis ${\mathbb{R}}$ is decomposed throughout this paper into
an interior region ${\mathcal{I}}=(-\ell,\ell)$ and an exterior region
${\mathcal{E}}={\mathcal{E}}_{+}\cup{\mathcal{E}}_{-}$ with
${\mathcal{E}}_{-}=(-\infty,-\ell)$ and ${\mathcal{E}}_{+}=(\ell,\infty)$, and
two contact points ${x=\pm\ell}$.
\- For any function $f\in C(\overline{{\mathcal{E}}})$, we denote

$f_{\pm}=f_{|_{x=\pm\ell}},\qquad\llbracket f\rrbracket=f_{+}-f_{-}\quad\mbox{
and }\quad\langle f\rangle=\frac{1}{2}\big{(}f_{+}+f_{-}\big{)}.$

\- If $f\in C^{1}([0,T])$, we sometimes use the notation
$\dot{f}=\frac{d}{dt}f$.
\- We denote by ${\mathfrak{f}}_{\rm sw}$ the momentum flux associated with
the nonlinear shallow water equations,

(1) ${\mathfrak{f}}_{\rm
sw}=\left(\varepsilon\frac{q^{2}}{h}+\frac{h^{2}-1}{2\varepsilon}\right).$

##  2\. An augmented formulation of the wave-structure interaction equations

The goal of this section is to derive the augmented formulation of the wave-
structure equations that we shall use in Section 3 to propose numerical
schemes. We first sketch in §2.1 the main steps of the analysis of [4] that
led to a formulation of the problem as a transmission problem between the two
connected components of the fluid domain, and with transmission conditions
determined through the resolution of an ODE forced by a source term involving
the traces at the contact points of the surface elevation and of their second
order time derivative. We then remark in §2.2 that these traces solve
themselves a second order ODE, but which is forced by a source term which is
easier to compute. This observation is the key ingredient that allows us to
derive in §2.3 an augmented formulation. It has the same structure as the
formulation derived in §2.1, namely, it is a transmission problem coupled with
a forced ODE. The crucial difference is that this ODE does no longer require
the computation of the traces of the surface elevation at the contact points
and that it can easily be computed numerically. Finally, we show in §2.4 that
this augmented transmission problem can be rewritten as an initial value
problem, which is the structure the numerical schemes of Section 3 are based
on.

###  2.1. Reduction to a transmission problem coupled with scalar ODEs

We remind here the main steps of the derivation of the equations describing
the interactions of a partially immersed object with one-dimensional waves in
a regime where these waves can correctly be described by the Boussinesq-Abbott
equations (see [22]); the object is assumed to have vertical sidewalls and can
be either fixed, in forced vertical motion, or allowed to float freely in the
vertical direction under the action of the waves.

In dimensionless variables, the equations involve two coefficients
$\varepsilon$ and $\mu$, respectively called nonlinearity and shallowness
parameters, and that are defined as

$\varepsilon=\frac{\mbox{typical amplitude of the waves}}{\mbox{typical
depth}}\quad\mbox{ and }\quad\mu=\Big{(}\frac{\mbox{typical
depth}}{\mbox{typical horizontal scale}}\Big{)}^{2};$

in the weakly nonlinear shallow water regime in which the Boussinesq-Abbott
equations are known to provide a good approximation of the motion of the
waves, one has

$\mu\ll 1\quad\mbox{ and }\quad\varepsilon=O(\mu);$

these conditions are assumed throughout this article. For the sake of
conciseness, we also introduce the parameter $\kappa$ as

$\kappa=\big{(}\frac{\mu}{3}\big{)}^{1/2};$

this parameter plays an important role as it measures the size of the
dispersive boundary layers that appear in the analysis of mixed initial
boundary-value problems for the Boussinesq equations, which are a dispersive
perturbation of an hyperbolic system [10].

As displayed in Figure 1, in dimensionless coordinates, the surface of the
fluid is parametrized at time $t$ by the function
$x\in{\mathbb{R}}\mapsto\varepsilon\zeta(t,x)$, and the horizontal discharge
(the vertical integral of the horizontal component of the velocity field) at
time $t$ and position $x$ is denoted $q(t,x)$. We also sometimes denote by $h$
the water depth, $h=1+\varepsilon\zeta$. Finally, we denote by
$\underline{P}(t,x)$ the pressure at the surface of the fluid, namely,
$\underline{P}=P(t,x,\varepsilon\zeta(t,x))$ if $P$ denotes the pressure field
in the fluid.

Regarding the solid object, we denote by $\pm\ell$ the position of its
vertical sidewalls and by $\varepsilon\zeta_{\rm w}$ the parametrization on
$(-\ell,\ell)$ of its bottom (the subscript ”w” stands for ”wetted part”); we
also denote by $\varepsilon\delta(t)$ the vertical deviation of the object
from its equilibrium position, and by $h_{\rm eq}$ the water depth at rest.
These quantities are related through

$\zeta_{\rm w}(t,x)=\delta(t)+\frac{1}{\varepsilon}\big{(}h_{\rm
eq}(x)-1\big{)}.$

N.B. For the sake of simplicity, we assume throughout this article that the
center of mass is located at $\\{x=0\\}$ and that $h_{\rm eq}(x)$ is an even
function.

The Boussinesq-Abbott equations for the motion of the waves are given for
$t>0$, $x\in{\mathbb{R}}$ by

(2) $\begin{cases}\partial_{t}\zeta+\partial_{x}q=0,\\\
(1-\kappa^{2}\partial_{x}^{2})\partial_{t}q+\varepsilon\partial_{x}\big{(}\frac{1}{h}q^{2}\big{)}+h\partial_{x}\zeta=-\frac{1}{\varepsilon}h\partial_{x}\underline{P}\end{cases}$

(with $h=1+\varepsilon\zeta$). We now have to distinguish between the exterior
region ${\mathcal{E}}=(-\infty,-\ell)\cup(\ell,\infty)$ where the surface of
the water is in contact with the air and the interior region
${\mathcal{I}}=(-\ell,\ell)$ where it is in contact with the object,

  * •

In the exterior region ${\mathcal{E}}$, the surface elevation $\zeta$ is free,
but the surface pressure $\underline{P}$ is constrained, assumed to be equal
to the (constant) atmospheric pressure $P_{\rm atm}$,

$\underline{P}(t,x)=P_{\rm atm}\quad\mbox{ for }t>0,\quad x\in{\mathcal{E}};$

the right-hand side in the second equation of (2) therefore vanishes.

  * •

In the interior region ${\mathcal{E}}$, it is the reverse: the surface
elevation is constrained because it has to coincide with the bottom of the
object,

$\zeta(t,x)=\zeta_{\rm w}(t,x)\quad\mbox{ for }t>0,\quad x\in{\mathcal{I}},$

but there is no constraint on the surface pressure $\underline{P}$ which,
under the general approach of [24], can be understood as the Lagrange
multiplier associated with the constraint on the surface elevation. Plugging
the constraint equation in the first equation of (2) one directly gets that

$q(t,x)=-x\dot{\delta}+\langle q_{\rm i}\rangle(t),$

where $\langle q_{\rm i}\rangle$ is a time dependent function corresponding to
the average discharge over the interior region. Using this relation and
applying $\partial_{x}$ to the second equation in (2) provides an elliptic
equation for $\underline{P}$,

$-\partial_{x}\big{(}\frac{1}{\varepsilon}h_{\rm
w}\partial_{x}\underline{P}\big{)}=-\ddot{\delta}+\partial_{x}\big{[}h_{\rm
w}\partial_{x}\zeta_{\rm w}+\varepsilon\partial_{x}\big{(}\frac{1}{h_{\rm
w}}(-x\dot{\delta}+\langle q_{\rm i}\rangle\big{)}^{2})\big{]},$

for $x\in(-\ell,\ell)$ and with $h_{\rm w}=h_{\rm eq}+\varepsilon\delta$. If
we know the boundary values of $\underline{P}$ at $x=-\ell+0$ and $x=\ell-0$,
this elliptic equation can be solved and it provides an expression for
$\underline{P}$ in terms of $h_{\rm eq}$, $\delta$, $\langle q_{\rm i}\rangle$
and of this boundary data. Using this expression in the second equation of (2)
then provides an expression for $\frac{d}{dt}\langle q_{\rm i}\rangle$ in
terms of the same quantities.

We also need coupling conditions at the contact points $x=\mp\ell$ between the
exterior and interior region. There are two of them,

  * •

Continuity of the horizontal discharge. Taking into account the expression of
the discharge derived above in the interior region, this condition yields

(3) $q(t,-\ell-0)=\ell\dot{\delta}+\langle q_{\rm i}\rangle(t)\quad\mbox{ and
}\quad q(t,\ell+0)=-\ell\dot{\delta}+\langle q_{\rm i}\rangle(t).$

  * •

Conservation of the total energy. Imposing conservation of the total (i.e.,
fluid+solid) energy classically provides the boundary data needed to solve the
elliptic equation derived for the surface pressure in the interior region [28,
7, 10]. We refer to [4] for the derivation of these boundary data in the
present context, but do not provide it here explicitly for the sake of
conciseness.

To summarize, we have the standard Boussinesq-Abbott equation in the exterior
region, with boundary condition on the discharge $q$ at $\mp\ell$ that are
given in terms of two functions of time, namely, $\langle q_{\rm i}\rangle$
and $\delta$. As said above, the fact that the elliptic equation for the
pressure in the interior region has been solved provides an evolution equation
for $\langle q_{\rm i}\rangle$; the last thing to do is therefore to determine
$\delta$. If the object is fixed or in forced motion then $\delta$ is given;
otherwise, it is of course given by Newton’s equation. The three cases can be
considered simultaneously by allowing an external force to be applied to the
solid (if the solid is fixed or in forced motion, this external force $F_{\rm
ext}$ represents the vertical force exerted on the solid to maintain it fixed
or with the desired motion). The outcome of this analysis, as shown111The
presence of the external force is not taken into account in that reference. It
is however straightforward to add it in Newton’s equation; note that in the
present dimensionless setting, the force has been nondimensionalized by
$2\ell\rho g$. in Theorem 3.1 of [4] is that the wave-structure interaction
problem under consideration can be reduced to a transmission problem. Using
the notations

$\langle f\rangle=\frac{1}{2}\big{(}f(\ell)-f(-\ell)\big{)}\quad\mbox{ and
}\quad\llbracket f\rrbracket=f(\ell)-f(-\ell)$

for all $f\in C((-\infty,-\ell]\cup[\ell,\infty))$, this transmission problem
can be written

(4) $\begin{cases}\partial_{t}\zeta+\partial_{x}q=0,\\\
(1-\kappa^{2}\partial_{x}^{2})\partial_{t}q+\partial_{x}{\mathfrak{f}}_{\rm
sw}=0\end{cases}\quad\mbox{ for }\quad t>0,\quad x\in{\mathcal{E}}$

where ${\mathfrak{f}}_{\rm sw}$ is the shallow water momentum flux given by
(1), and with transmission conditions across the floating object given by

(5) $\langle q\rangle=\langle q_{\rm i}\rangle\quad\mbox{ and }\quad\llbracket
q\rrbracket=-2\ell\dot{\delta},$

where $\langle q_{\rm i}\rangle$ and $\delta$ are functions of time solving

(6) $\displaystyle\alpha(\varepsilon\delta)\frac{d}{dt}\langle q_{\rm
i}\rangle+\varepsilon\alpha^{\prime}(\varepsilon\delta)\dot{\delta}\langle
q_{\rm i}\rangle$
$\displaystyle=-\frac{1}{2\ell}\llbracket\zeta+{\mathfrak{G}}\rrbracket,$ (7)
$\displaystyle\tau_{\kappa}(\varepsilon\delta)^{2}\ddot{\delta}+\delta-\varepsilon\beta(\varepsilon\delta)\dot{\delta}^{2}-\varepsilon\frac{1}{2}\alpha^{\prime}(\varepsilon\delta)\langle
q_{\rm i}\rangle^{2}$ $\displaystyle=\langle\zeta+{\mathfrak{G}}\rangle+F_{\rm
ext},$

where ${\mathfrak{G}}$ is the function defined on ${\mathcal{E}}$ by

(8)
${\mathfrak{G}}=\varepsilon\frac{1}{2}\frac{q^{2}}{h^{2}}-\kappa^{2}\frac{1}{h}\partial_{x}\partial_{t}q,$

and where the explicit expression of the functions $\alpha$, $\tau_{\kappa}$
and $\beta$, of no importance at this point of the discussion, are provided in
§A.1 of Appendix A. We just want to emphasize that the coefficient
$\tau_{\kappa}(\varepsilon\delta)^{2}$ in front of $\ddot{\delta}$ in (7)
takes into account the contribution of the added mass effect (when a solid
moves in a fluid, not only must it accelerate its own mass but also the mass
of the fluid around it).

The initial value problem corresponding to (4)-(7) is studied and solved in
[4]. Its structure is that of a transmission problem coupled with a set of
ODEs on $\langle q_{\rm i}\rangle$ and $\delta$. This coupling acts in both
ways: on the one hand, it is necessary to know $\langle q_{\rm i}\rangle$ and
$\delta$ in order to solve the transmission problem (4)-(5) and on the other
hand, one needs to know the solution $(\zeta,q)$ of this transmission problem
to compute the source term in the right-hand side of (6)-(7). From the
numerical view point, this last step is not easy to treat since one has to
compute the numerical trace of $\zeta$ and $\partial_{t}\partial_{x}q$ at the
contact points $x=\pm\ell$. The key ingredient we propose here to overcome
this difficulty is to work with an augmented formulation of the problem, with
additional functions of time involved in the system of ODEs for $\delta$ and
$\langle q_{\rm i}\rangle$, but where the computation of such traces is no
longer needed.

###  2.2. The trace equations

The source terms in the right-hand sides of (6)-(7) involve the trace of
$\zeta+{\mathfrak{G}}$ at $x=\pm\ell$, with ${\mathfrak{G}}$ given by (8).
Since (5) implies that $q_{|_{x=\pm\ell}}=\mp\ell\dot{\delta}+\langle q_{\rm
i}\rangle$ and remarking that one deduces from the first equation of (4) that
$\partial_{x}\partial_{t}q=-\partial_{t}^{2}\zeta$, we have

${\mathfrak{G}}_{|_{x=\pm\ell}}=\varepsilon\frac{1}{2}\Big{(}\frac{\mp\ell\dot{\delta}+\langle
q_{\rm
i}\rangle}{1+\varepsilon\zeta_{\pm}}\Big{)}^{2}+\kappa^{2}\frac{1}{1+\varepsilon\zeta_{\pm}}\ddot{\zeta}_{\pm},$

with $\zeta_{\pm}:=\zeta_{|_{x=\pm\ell}}$. The difficulty therefore lies in
the computation of the trace of $\zeta$ at $x=\pm\ell$ and of their second
time derivative. The augmented formulation consists in treating $\zeta_{\pm}$
as a new unknown function of time instead of getting it by taking the traces
of $\zeta$ at the contact points. This is made possible by the following
proposition which provides a second order ODE satisfied by $\zeta_{+}$ and
$\zeta_{-}$. This requires first the introduction of the Dirichlet and Neumann
inverses of the operator $(1-\kappa^{2}\partial_{x}^{2})$ on ${\mathcal{E}}$,
respectively denoted by $R_{0}$ and $R_{1}$. They are defined for all $F\in
L^{2}({\mathcal{E}})$ by

$R_{0}F=u\quad\mbox{ with
}\begin{cases}(1-\kappa^{2}\partial_{x}^{2})u=F&\mbox{ on }{\mathcal{E}},\\\
u_{|_{x=\pm\ell}}=0,\end{cases}$

and

(9) $R_{1}F=v\quad\mbox{ with
}\begin{cases}(1-\kappa^{2}\partial_{x}^{2})v=F&\mbox{ on }{\mathcal{E}},\\\
\partial_{x}v_{|_{x=\pm\ell}}=0.\end{cases}$

We can now state the following proposition. Note that the ODEs satisfied by
$\zeta_{\pm}$ only make sense in the presence of dispersion ($\kappa>0$).

######  Proposition 2.1.

Let $f$ and $g$ be two continuous functions of time. If $(\zeta,q)$ is a
smooth solution to

(10) $\begin{cases}\partial_{t}\zeta+\partial_{x}q=0,\\\
(1-\kappa^{2}\partial_{x}^{2})\partial_{t}q+\partial_{x}{\mathfrak{f}}_{\rm
sw}=0,\end{cases}t>0,\quad x\in{\mathcal{E}}$

with ${\mathfrak{f}}_{\rm sw}$ as in (1) and with transmission conditions

(11) $\langle q\rangle(t)=f(t)\quad\mbox{ and }\quad\llbracket
q\rrbracket(t)=2g(t),\qquad t>0$

then $\zeta_{\pm}=\zeta_{|_{x=\pm\ell}}$ solve the ODEs

(12)
$\begin{cases}\displaystyle\partial_{t}^{2}\zeta_{+}+\frac{1}{\kappa^{2}}\zeta_{+}+\frac{\varepsilon}{\kappa^{2}}\big{(}\frac{1}{2}\zeta_{+}^{2}+\frac{(f+g)^{2}}{1+\varepsilon\zeta_{+}}\big{)}=\frac{1}{\kappa^{2}}({R}_{1}{\mathfrak{f}}_{\rm
sw})_{+}+\frac{1}{\kappa}(\dot{f}+\dot{g}),\\\
\displaystyle\partial_{t}^{2}\zeta_{-}+\frac{1}{\kappa^{2}}\zeta_{-}+\frac{\varepsilon}{\kappa^{2}}\big{(}\frac{1}{2}\zeta_{-}^{2}+\frac{(f-g)^{2}}{1+\varepsilon\zeta_{-}}\big{)}=\frac{1}{\kappa^{2}}({R}_{1}{\mathfrak{f}}_{\rm
sw})_{-}-\frac{1}{\kappa}(\dot{f}-\dot{g}),\end{cases}$

where we used the notation $({R}_{1}{\mathfrak{f}}_{\rm
sw})_{\pm}=({R}_{1}{\mathfrak{f}}_{\rm sw})_{|_{x=\pm\ell}}$.

###### Proof.

Applying $R_{0}$ to the second equation in (10) and using the boundary
condition (11), one gets

$\partial_{t}q+R_{0}\partial_{x}{\mathfrak{f}}_{\rm
sw}=(\dot{f}\pm\dot{g})\exp(-\frac{|x\mp\ell|}{\kappa})\quad\mbox{ on
}\quad{\mathcal{E}}^{\pm}.$

Remarking further that $R_{0}\partial_{x}=\partial_{x}R_{1}$, the problem is
therefore reduced to

(13) $\begin{cases}\partial_{t}\zeta+\partial_{x}q=0,\\\
\partial_{t}q+\partial_{x}R_{1}{\mathfrak{f}}_{\rm
sw}=(\dot{f}\pm\dot{g})\exp(-\frac{|x\mp\ell|}{\kappa}).\end{cases}$

Differentiating with respect to $x$ the second equation of (13) and using the
fact that $\partial_{t}\partial_{x}q=-\partial_{t}^{2}\zeta$, one gets

$-\partial_{t}^{2}\zeta+\partial_{x}^{2}R_{1}{\mathfrak{f}}_{\rm
sw}=\mp\frac{1}{\kappa}(\dot{f}\pm\dot{g})\exp\big{(}-\frac{1}{\kappa}(|x\mp\ell|)\big{)}.$

Since moreover
$\partial_{x}^{2}=-\frac{1}{\kappa^{2}}(1-\kappa^{2}\partial_{x}^{2})+\frac{1}{\kappa^{2}}$,
we deduce that

$-\partial_{t}^{2}\zeta-\frac{1}{\kappa^{2}}{\mathfrak{f}}_{\rm
sw}+\frac{1}{\kappa^{2}}R_{1}{\mathfrak{f}}_{\rm
sw}=\mp\frac{1}{\kappa}(\dot{f}\pm\dot{g})\exp\big{(}-\frac{1}{\kappa}(|x\mp\ell|)\big{)}.$

Taking the trace at $x=\pm\ell$, and substituting ${{\mathfrak{f}}_{\rm
sw}}_{|_{x=\pm\ell}}=\zeta_{\pm}+\varepsilon\big{(}\frac{1}{2}\zeta_{\pm}^{2}+\frac{(f\pm
g)^{2}}{1+\varepsilon\zeta_{\pm}}\big{)}$, we obtain the equations stated in
the proposition. ∎

###  2.3. The augmented formulation

Proposition 2.1 can be applied to the wave-structure interaction system
(4)-(7) with $f=\langle q_{\rm i}\rangle$ and $g=-\ell\dot{\delta}$. Together
with (6)-(7), this shows that $\langle q_{\rm i}\rangle$, $\delta$,
$\zeta_{+}$ and $\zeta_{-}$ solve the second order differential system

(14)
${\mathcal{M}}[\varepsilon\delta,\varepsilon\zeta_{\pm}]\frac{d}{dt}\begin{pmatrix}\displaystyle\langle
q_{\rm i}\rangle\\\ \displaystyle\dot{\delta}\\\
\displaystyle\dot{{{\zeta}}}{}_{+}\\\
\displaystyle\dot{{\zeta}}{}_{-}\end{pmatrix}+\begin{pmatrix}\displaystyle\frac{1}{2\ell}\llbracket\zeta\rrbracket\\\
\displaystyle\delta-\langle\zeta\rangle\\\ \displaystyle{{{\zeta}}}{}_{+}\\\
\displaystyle{{\zeta}}{}_{-}\end{pmatrix}=\varepsilon{\boldsymbol{\mathfrak{Q}}}[\varepsilon\delta,\varepsilon\zeta_{\pm}](\langle
q_{\rm i}\rangle,\dot{\delta},\zeta_{\pm})+\begin{pmatrix}0\\\ F_{\rm ext}\\\
\displaystyle({R}_{1}{\mathfrak{f}}_{\rm sw})_{+}\\\
\displaystyle({R}_{1}{\mathfrak{f}}_{\rm sw})_{-}\end{pmatrix},$

where ${\mathcal{M}}[\varepsilon\delta,\varepsilon\zeta_{\pm}]$ is the
invertible matrix

(15)
${\mathcal{M}}[\varepsilon\delta,\varepsilon\zeta_{\pm}]:=\left(\begin{array}[]{cc|cc}\alpha(\varepsilon\delta)&0&\frac{\kappa^{2}}{2\ell}\frac{1}{1+\varepsilon\zeta_{+}}&-\frac{\kappa^{2}}{2\ell}\frac{1}{1+\varepsilon\zeta_{-}}\\\
0&\tau_{\kappa}(\varepsilon\delta)^{2}&-\frac{1}{2}\frac{\kappa^{2}}{1+\varepsilon\zeta_{+}}&-\frac{1}{2}\frac{\kappa^{2}}{1+\varepsilon\zeta_{-}}\\\
\hline\cr-{\kappa}&\ell\kappa&\kappa^{2}&0\\\
\kappa&\ell\kappa&0&\kappa^{2}\end{array}\right),$

while
${\boldsymbol{\mathfrak{Q}}}[\varepsilon\delta,\varepsilon\zeta_{\pm}](\langle
q_{\rm i}\rangle,\dot{\delta},\zeta_{\pm})$ is a four-dimensional vector whose
entries are quadratic forms in $(\langle q_{\rm
i}\rangle,\dot{\delta},\zeta_{\pm})$ with coefficients depending on
$\varepsilon\delta$, $\varepsilon\zeta_{+}$ and $\varepsilon\zeta_{-}$ (the
exact expression of these terms is of no importance at this point, and we
refer the reader to Appendix A). The second order differential system (14) can
classically be transformed into a first order ODE on $\langle q_{\rm
i}\rangle$, $\dot{\delta}$, $\dot{\zeta}_{+}$, $\dot{\zeta}_{-}$, $\delta$,
$\zeta_{+}$ and $\zeta_{-}$ with forcing terms $(R_{1}{\mathfrak{f}}_{\rm
sw})_{\pm}$ and $F_{\rm ext}$ (see Appendix A.2). The augmented formulation is
obtained by replacing $\zeta_{+}$ and $\zeta_{-}$ by two additional unknowns
${\underline{\zeta}{}_{+}}$ and $\underline{\zeta}{}_{-}$ in this first order
ODE. It reads therefore

(16) $\begin{cases}\partial_{t}\zeta+\partial_{x}q=0,\\\
(1-\kappa^{2}\partial_{x}^{2})\partial_{t}q+\partial_{x}{\mathfrak{f}}_{\rm
sw}=0\end{cases}\quad\mbox{ for }\quad t>0,\quad x\in{\mathcal{E}},$

where ${\mathfrak{f}}_{\rm sw}$ is as in (1), and with transmission conditions
across the floating object given by

(17) $\langle q\rangle=\langle q_{\rm i}\rangle\quad\mbox{ and
}\quad\llbracket q\rrbracket=-2\ell\dot{\delta},$

where $\langle q_{\rm i}\rangle$ and $\delta$ are functions of time determined
by the first order ODE

(18) $\frac{d}{dt}\Theta={\mathcal{G}}\big{(}\Theta,(R_{1}{\mathfrak{f}}_{\rm
sw})_{+},(R_{1}{\mathfrak{f}}_{\rm sw})_{-},F_{\rm ext}\big{)},$

with $\Theta:=\big{(}\langle q_{\rm
i}\rangle,\dot{\delta},\dot{\underline{\zeta}{}_{+}},\dot{}\underline{\zeta}{}_{-},\delta,{\underline{\zeta}{}_{+}},\underline{\zeta}{}_{-}\big{)}^{\rm
T}$ and where ${\mathcal{G}}$ is a smooth function of its arguments and whose
exact expression is given in Appendix A. It is a consequence of Proposition
2.2 below that if the initial data for $\underline{\zeta}{}_{\pm}$ and
$\dot{\underline{\zeta}}{}_{\pm}$ are chosen appropriately, then
$\zeta_{\pm}=\underline{\zeta}{}_{\pm}$ and
$\dot{\zeta}_{\pm}=\dot{\underline{\zeta}}{}_{\pm}$ for all times, as expected
(see Proposition 2.2 below).

######  Remark 2.1.

The difference between the augmented formulation (16)-(18) and the original
formulation (4)-(7) lies in the ODE used to determine the functions $\delta$
and $\langle q_{\rm i}\rangle$ involved in the transmission conditions. In the
original formulation, one has a first order $3$-dimensional ODE (on $\delta$,
$\dot{\delta}$ and $\langle q_{\rm i}\rangle$), which is forced by $F_{\rm
ext}$, $(R_{1}{\mathfrak{f}}_{\rm sw})_{|_{x=\pm\ell}}$,
$\zeta_{|_{x=\pm\ell}}$ and $\frac{d^{2}}{dt^{2}}\zeta_{|_{x=\pm\ell}}$. In
the augmented formulation, the ODE is of higher dimension, namely, it is a
first order $7$-dimensional ODE (on $\delta$, $\dot{\delta}$, $\langle q_{\rm
i}\rangle$, $\underline{\zeta}{}_{\pm}$ and
$\dot{\underline{\zeta}}{}_{\pm}$), but it is forced only by $F_{\rm ext}$ and
$(R_{1}{\mathfrak{f}}_{\rm sw})_{|_{x=\pm\ell}}$. These two quantities do not
raise any difficulty since $F_{\rm ext}$ is a given external force and
$(R_{1}{\mathfrak{f}}_{\rm sw})_{|_{x=\pm\ell}}$ can easily be computed
numerically (see §3.1.2 below), contrary to the traces of $\zeta$ and
$\partial_{t}^{2}\zeta$ at the contact points that appear in the original
formulation and that are very delicate to compute.

###  2.4. Transformation into an initial value problem

We reformulate in this section the wave-structure transmission problem
(16)-(18) in the form of an initial value problem that is easier to handle
from a numerical point of view. This formulation is the new augmented
formulation (with additional variables $\underline{\zeta}{}_{\pm}$) we shall
base our numerical schemes on. In [4] (see also the lecture notes [23]), the
well-posedness of the standard formulation (4)-(7) is proved, and it could
similarly be obtained for the augmented formulation; for the sake of
conciseness, we do not give here such a result and just prove that both
formulations have the same regular solutions, and that the additional
variables $\underline{\zeta}{}_{\pm}$ coincide with the traces
$\zeta_{|_{x=\pm\ell}}$ under certain compatibility conditions on the initial
data. We use the following notation for the source term in the reformulated
momentum equation,

(19) ${\mathcal{S}}_{\pm}\big{(}\Theta,(R_{1}{\mathfrak{f}}_{\rm
sw})_{\pm},F_{\rm
ext}\big{)}:={\mathcal{G}_{1}}\big{(}\Theta,(R_{1}{\mathfrak{f}}_{\rm
sw})_{\pm},F_{\rm
ext}\big{)}\mp\ell{\mathcal{G}_{2}}\big{(}\Theta,(R_{1}{\mathfrak{f}}_{\rm
sw})_{\pm},F_{\rm ext}\big{)},$

where ${\mathcal{G}}_{1}$ and ${\mathcal{G}}_{2}$ denote the first two
components of the mapping ${\mathcal{G}}$ in the right-hand side of the ODE
(18). We also recall that we denote ${\mathcal{E}}^{-}=(-\infty,-\ell)$ and
${\mathcal{E}}^{+}=(\ell,\infty)$ the two connected components of the fluid
domain ${\mathcal{E}}$.

######  Proposition 2.2.

Let $(\zeta,q)$ and $\Theta=\big{(}\langle q_{\rm
i}\rangle,\dot{\delta},\dot{\underline{\zeta}}{}_{+},\dot{\underline{\zeta}}{}_{-},\delta,{\underline{\zeta}{}_{+}},\underline{\zeta}{}_{-}\big{)}^{\rm
T}$ be a regular solution to the transmission problem (16)-(18) with initial
data $U=(\zeta^{\rm in},q^{\rm in})$ and $\Theta^{\rm in}$. Then $(\zeta,q)$
and $\Theta$ also solve the initial value problem

(20) $\begin{cases}\partial_{t}\zeta+\partial_{x}q=0,\\\
\partial_{t}q+\partial_{x}R_{1}{\mathfrak{f}}_{\rm
sw}={\mathcal{S}}_{\pm}\big{(}\Theta,(R_{1}{\mathfrak{f}}_{\rm
sw})_{\pm},F_{\rm
ext}\big{)}\exp(-\frac{1}{\kappa}|x\mp\ell|)\end{cases}\quad\mbox{ on
}{\mathcal{E}}^{\pm},$

and

(21) $\frac{d}{dt}\Theta={\mathcal{G}}\big{(}\Theta,(R_{1}{\mathfrak{f}}_{\rm
sw})_{+},(R_{1}{\mathfrak{f}}_{\rm sw})_{-},F_{\rm ext}\big{)}.$

The converse is true, provided that the initial data satisfy the compatibility
conditions

(22) $\langle q^{\rm in}\rangle=\Theta^{\rm in}_{1},\qquad\llbracket q^{\rm
in}\rrbracket=-2\ell\Theta^{\rm in}_{2}.$

If moreover the initial data also satisfy

(23) $\zeta^{\rm in}_{|_{x=\ell}}=\Theta^{\rm in}_{6},\qquad\zeta^{\rm
in}_{|_{x=-\ell}}=\Theta^{\rm in}_{7},\qquad-(\partial_{x}q^{\rm
in})_{|_{x=\ell}}=\Theta^{\rm in}_{3},\qquad-(\partial_{x}q^{\rm
in})_{|_{x=-\ell}}=\Theta^{\rm in}_{4},$

then for all times, one has $\zeta_{|_{x=\pm\ell}}=\underline{\zeta}{}_{\pm}$.

######  Remark 2.2.

The proposition deals with the most general situation to cover in a unified
way all the situations considered in this article. It can be simplified in
various cases, as shown in Appendix A. For instance,

  * •

When the object is freely floating, one takes $F_{\rm ext}\equiv 0$.

  * •

When the data and the object are symmetric with respect to the vertical axis
$x=0$, then $\langle q_{\rm i}\rangle=0$.

  * •

If the object is in forced motion, i.e. if $\delta\equiv\delta_{\rm forced}$
for some given function $\delta_{\rm forced}$, then the ODE (21) can be
reduced to a $5$-dimensional ODE on $\big{(}\langle q_{\rm
i}\rangle,\dot{\underline{\zeta}}_{+},\dot{\underline{\zeta}}_{+},{\underline{\zeta}}_{+},{\underline{\zeta}}_{-}\big{)}^{\rm
T}$. Note that in this situation, an external force is needed to maintain the
object fixed (the exact expression of this force is derived in Remark A.3).

###### Proof.

Let us first prove the direct implication. Proceeding as in the proof of
Proposition 2.1, we can rewrite the second equation of (16) on each component
${\mathcal{E}}^{\pm}$ of ${\mathcal{E}}$ under the form

$\partial_{t}q+\partial_{x}R_{1}{\mathfrak{f}}_{\rm
sw}=\frac{d}{dt}(q_{|_{x=\pm\ell}})\exp(-\frac{1}{\kappa}|x\mp\ell|)\quad\mbox{
on }\quad{\mathcal{E}}^{\pm}.$

From the transmission conditions (17), we have
$q_{|_{x=\pm\ell}}=\mp\ell\delta+\langle q_{\rm i}\rangle$ so that the result
follows from the observation that, owing to (18), one has

$\frac{d}{dt}\langle q_{\rm
i}\rangle={\mathcal{G}}_{1}(\Theta,(R_{1}{\mathfrak{f}}_{\rm sw})_{\pm},F_{\rm
ext})\quad\mbox{ and
}\quad\frac{d}{dt}\dot{\delta}={\mathcal{G}}_{2}(\Theta,(R_{1}{\mathfrak{f}}_{\rm
sw})_{\pm},F_{\rm ext}).$

Conversely, if $(\zeta,q)$ solves (20), it suffices to apply
$(1-\kappa^{2}\partial_{x}^{2})$ to the second equation to show that
$(\zeta,q)$ solves (16). The equation (18) on $\Theta$ is the same as (21), so
that the only thing we need to prove is that the transmission conditions (17)
hold. Taking the trace of the second equation of (20) at the contact points
and taking the average and the jump, we find that

$\frac{d}{dt}\langle
q\rangle={\mathcal{G}_{1}}(\Theta,(R_{1}{\mathfrak{f}}_{\rm sw})_{\pm},F_{\rm
ext})\quad\mbox{ and }\quad\frac{d}{dt}\llbracket
q\rrbracket=-2\ell{\mathcal{G}}_{2}(\Theta,(R_{1}{\mathfrak{f}}_{\rm
sw})_{\pm},F_{\rm ext}),$

or equivalently (from the definition of ${\mathcal{G}}$),

$\frac{d}{dt}\langle q\rangle=\frac{d}{dt}\langle q_{\rm i}\rangle\quad\mbox{
and }\quad\frac{d}{dt}\llbracket
q\rrbracket=\frac{d}{dt}\big{(}-2\ell\dot{\delta}\big{)}.$

This shows that the time derivative of the transmission conditions (17) are
satisfied; the compatibility conditions (22) show moreover that the
transmission condition condition is satisfied at $t=0$. It is therefore
satisfied for all times.
For the last assertion, we can use Proposition 2.1 to show that
$\zeta_{|_{x=\pm\ell}}$ and $\underline{\zeta}{}_{\pm}$ satisfy the same
second order ODE in time. The additional condition (23) ensures that these
initial data and the initial value of the first time derivative coincide (we
also used the first equation of (16) to substitute
$\frac{d}{dt}(\zeta_{|_{x=\pm\ell}})=-(\partial_{x}q)_{|_{x=\pm\ell}}$). They
are therefore identical for all times. ∎

##  3\. Numerical schemes

We present in this section one first order and one second order numerical
scheme for the resolution of the augmented formulations derived in this
article. We explain these schemes for the general formulation (20)-(21). We
recall that these equations are conservation laws with a nonlocal flux and and
an exponentially localized source term,

(24)
$\partial_{t}U+\partial_{x}\big{(}{\mathfrak{F}}_{\kappa}(U)\big{)}={\mathcal{S}}_{\pm}(\Theta,(R_{1}{\mathfrak{f}}_{\rm
sw})_{-},(R_{1}{\mathfrak{f}}_{\rm sw})_{+},F_{\rm
ext}){\mathfrak{b}}(x\mp\ell)\qquad\mbox{ in }{\mathcal{E}}_{\pm}$

with $U=(\zeta,q)^{T}$ and where ${\mathfrak{F}}_{\kappa}(U)$ is the nonlocal
flux given by

(25) ${\mathfrak{F}}_{\kappa}(U)=\big{(}q,R_{1}{\mathfrak{f}}_{\rm
sw}\big{)}^{T},$

while the source terms ${\mathcal{S}}_{\pm}$ are as in (19) and
${\mathfrak{b}}$ is the shape of source term

(26) ${\mathfrak{b}}(x)=\big{(}0,\exp(-\frac{|x|}{\kappa})\big{)}^{\rm T}.$

The quantity $\Theta$ is defined as $\Theta=\big{(}\langle q_{\rm
i}\rangle,\dot{\delta},\dot{\underline{\zeta}}{}_{+},\dot{\underline{\zeta}}{}_{-},\delta,{\underline{\zeta}{}_{+}},\underline{\zeta}{}_{-}\big{)}^{\rm
T}$ and solves a system of $7$ first order ODEs forced by
$(R_{1}{\mathfrak{f}}_{\rm sw})_{+}$, $(R_{1}{\mathfrak{f}}_{\rm sw})_{-}$ and
$F_{\rm ext}$,

(27) $\frac{d}{dt}\Theta={\mathcal{G}}\big{(}\Theta,(R_{1}{\mathfrak{f}}_{\rm
sw})_{+},(R_{1}{\mathfrak{f}}_{\rm sw})_{-},F_{\rm ext}\big{)}.$

######  Remark 3.1.

As explained in Remark 2.2, in some of the examples considered in this paper,
the ODE (27) can be reduced to a possibly lower dimensional ODE; we refer to
Appendix A where such simplifications are derived.

###  3.1. Notations

We gather here the main notations used to write our numerical schemes. We
first set our notations for the discretized quantities, and then explain how
we define the discrete version of the nonlocal operator $R_{1}$ defined in
(9).

####  3.1.1. Discretization

Figure 2. Space discretzation

We denote by $\Delta x$ the mesh size and decompose the two components
${\mathcal{E}}^{-}$ and ${\mathcal{E}}^{+}$ of the exterior domain into a
disjoint union of cells (see figure 2),

${\mathcal{E}}^{-}=\Big{(}\bigcup_{i=-\infty}^{-1}{\mathcal{C}}_{i}\Big{)}\cup{\mathcal{C}}_{0^{-}}\quad\mbox{
and
}\quad{\mathcal{E}}^{+}={\mathcal{C}}_{0^{+}}\cup\Big{(}\bigcup_{i=1}^{\infty}{\mathcal{C}}_{i}\Big{)}$

with

${\mathcal{C}}_{i}=(x_{i-1/2},x_{i+1/2})\mbox{ if }i\neq 0\quad\mbox{ and
}\quad{\mathcal{C}}_{0^{-}}=(x_{-1/2},-\ell),\qquad{\mathcal{C}}_{0^{+}}=(\ell,x_{1/2})$

and where

$x_{i+1/2}=-\ell+(i+1/2)\Delta x\quad\mbox{ if }i<0\quad\mbox{ and }\quad
x_{i-1/2}=\ell+(i-1/2)\Delta x\quad\mbox{ if }i>0.$

######  Remark 3.2.

Of course, the numerical domain is of finite size but we work with large
enough domains so that the influence of the left and right boundaries of the
numerical domain are not seen in the computations near the solid object. For
the sake of clarity, we do not mention these boundaries in the presentation of
the numerical scheme.

We also write $\Delta t>0$ the time stepping and denote by

$U^{n}=U(n\Delta t),\qquad\Theta^{n}=\Theta(n\Delta t)\quad\mbox{ and }\quad
F_{\rm ext}^{n}=F_{\rm ext}(n\Delta t)$

the values of $U=(\zeta,q)^{\rm T}$, of the ${\mathbb{R}}^{7}$-valued vector
$\Theta$ involved in the ODE (27), and of the external force $F_{\rm ext}$ at
each time step. We further denote by $U^{n}_{i}$
($i\in\big{(}{\mathbb{Z}}\backslash\\{0\\}\big{)}\cup\\{0^{-},0^{+}\\}$) the
approximation of ${\mathcal{U}}^{n}$ in the middle of the cell
${\mathcal{C}}_{i}$ furnished by the numerical scheme.

####  3.1.2. About the nonlocal operator $R_{1}$

The equations (24)-(27) involve the quantities $R_{1}{\mathfrak{f}}_{\rm sw}$
and $(R_{1}{\mathfrak{f}}_{\rm sw})_{\pm}$, where we recall that $R_{1}$ is
the inverse of $(1-\kappa^{2}\partial_{x}^{2})$ on
${\mathcal{E}}^{-}\cup{\mathcal{E}}^{+}$ with Neumann boundary condition at
$\pm\ell$, as defined in (9), and that $(R_{1}{\mathfrak{f}}_{\rm sw})_{\pm}$
stands for the trace of $R_{1}{\mathfrak{f}}_{\rm sw}$ at $\pm\ell$.

We keep the same notation $R_{1}$ for the discrete inverse of the operator
$(1-\kappa\partial_{x}^{2})$ with homogeneous Neumann condition at the
boundary. We use here a standard centered second order finite difference
approximation for the discretization of $\partial_{x}^{2}$. More precisely, if
$F=(f_{i})_{|i|\geq 1}$, we denote by $R_{1}F$ the vector $R_{1}F=V$ where
$V=(v_{i})_{|i|\geq 1}$ is given by the resolution of the equations

$v_{i}-\kappa^{2}\frac{v_{i+1}-2v_{i}+v_{i-1}}{\Delta
x^{2}}=f_{i},\qquad|i|\geq 2$

while, for $i=\pm 1$ a second order discretization of the Neumann boundary
condition leads to

$v_{-1}-\kappa^{2}\frac{2}{3}\frac{v_{-2}-v_{-1}}{\Delta
x^{2}}=f_{-1}\quad\mbox{ and }\quad
v_{1}-\kappa^{2}\frac{2}{3}\frac{v_{2}-v_{1}}{\Delta x^{2}}=f_{1}.$

Similarly, we still denote by $({R}_{1}F)_{\pm}$ the discrete version of the
traces ${R}_{1}F$ at the boundaries; they are naturally defined by the second
order approximation

(28)
$\big{(}{R}_{1}F\big{)}_{-}=\frac{4}{3}v_{-1}-\frac{1}{3}v_{-2}\quad\mbox{ and
}\quad\big{(}{R}_{1}F\big{)}_{+}=\frac{4}{3}v_{1}-\frac{1}{3}v_{2}.$

###  3.2. A first order scheme

We propose here an adaptation of the Lax-Friedrichs scheme for the
conservation laws with nonlocal flux (24). This scheme is an extension of the
scheme used in [26] for the numerical simulation of the Boussinesq equations
with generating boundary condition (i.e. with data on $\zeta$ at the entrance
of the numerical domain). It reads

(29) $\frac{U_{i}^{n+1}-U_{i}^{n}}{\Delta t}+\frac{1}{\Delta
x}\big{(}{\mathfrak{F}}^{n}_{\kappa,i+1/2}-{\mathfrak{F}}^{n}_{\kappa,i-1/2}\big{)}={\mathcal{S}}_{\pm}^{n}{\mathfrak{b}}_{i},\qquad\pm
i\geq 1,\quad n\geq 0,$

with

(30)
${\mathcal{S}}_{\pm}^{n}={\mathcal{S}}_{\pm}(\Theta^{n},(R_{1}{\mathfrak{f}}_{\rm
sw}^{n})_{+},(R_{1}{\mathfrak{f}}_{\rm sw}^{n})_{-},F_{\rm
ext}^{n})\quad\mbox{ and }\quad{\mathfrak{b}}_{i}={\mathfrak{b}}(i\Delta
x)\quad\mbox{ if }\pm i>0;$

the discrete flux correspond to the Lax-Friedrichs scheme,

(31)
$\begin{cases}{\mathfrak{F}}_{\kappa,i+1/2}^{n}=\frac{1}{2}\big{(}{\mathfrak{F}}_{\kappa,i+1}^{n}+{\mathfrak{F}}_{\kappa,i}^{n}\big{)}-\frac{\Delta
x}{2\Delta t}\big{(}U_{i+1}^{n}-U_{i}^{n}\big{)}&\mbox{if }i\leq-2,\\\
{\mathfrak{F}}_{\kappa,i-1/2}^{n}=\frac{1}{2}\big{(}{\mathfrak{F}}_{\kappa,i}^{n}+{\mathfrak{F}}_{\kappa,i-1}^{n}\big{)}-\frac{\Delta
x}{2\Delta t}\big{(}U_{i}^{n}-U_{i-1}^{n}\big{)}&\mbox{if }i\geq
2,\end{cases}$

with the notations

${\mathfrak{F}}_{\kappa,i}^{n}=\big{(}q_{i}^{n},(R_{1}{\mathfrak{f}}_{\rm
sw}^{n})_{i}\big{)}^{\rm T}\quad\mbox{ and }\quad{\mathfrak{f}}_{\rm
sw}^{n}={\mathfrak{f}}_{\rm sw}(U^{n});$

finally, for $i=\pm 1$, we must adapt (31) in the following way,

(32)
$\begin{cases}{\mathfrak{F}}_{\kappa,-1/2}^{n}=\frac{1}{2}\big{(}{\mathfrak{F}}_{\kappa,0^{-}}^{n}+{\mathfrak{F}}_{\kappa,-1}^{n}\big{)}-\frac{\Delta
x}{2\Delta t}\big{(}U_{0^{-}}^{n}-U_{-1}^{n}\big{)}\\\
{\mathfrak{F}}_{\kappa,1/2}^{n}=\frac{1}{2}\big{(}{\mathfrak{F}}_{\kappa,1}^{n}+{\mathfrak{F}}_{\kappa,0^{+}}^{n}\big{)}-\frac{\Delta
x}{2\Delta t}\big{(}U_{1}^{n}-U_{0^{+}}^{n}\big{)}\end{cases}$

with

(33)
${\mathfrak{F}}_{\kappa,0^{\pm}}^{n}=\big{(}q^{n}_{0^{\pm}},(R_{1}{\mathfrak{f}}_{\rm
sw}^{n})_{\pm}\big{)}^{\rm T};$

the component $(R_{1}{\mathfrak{f}}_{\rm sw}^{n})_{\pm}$ is computed according
to (28), but we still need to define $q_{0^{\pm}}^{n}$. By definition
$q^{n}_{0^{\pm}}$ is the approximation at time $n\Delta t$ of the trace of the
discharge $q$ at $\pm\ell$. From the transmission conditions (17) of the
continuous problem, we have $q_{|_{x=\pm\ell}}=\langle q_{\rm
i}\rangle\mp\ell\dot{\delta}$. Recalling also that $\langle q_{\rm i}\rangle$
and $\dot{\delta}$ are respectively the first and second components of
$\Theta$, this relation can be rewritten
$q_{|_{x=\pm\ell}}=\Theta_{1}\mp\ell\Theta_{2}$. At the discrete level, this
leads to the following definition for $q^{n}_{0^{\pm}}$,

(34) $q^{n}_{0^{\pm}}=\Theta^{n}_{1}\mp\ell\Theta^{n}_{2}.$

The equation (27) is discretized with a first-order explicit Euler scheme:

(35) $\frac{\Theta^{n+1}-\Theta^{n}}{\Delta
t}={\mathcal{G}}\big{(}\Theta^{n},(R_{1}{\mathfrak{f}}_{\rm
sw}^{n})_{+},(R_{1}{\mathfrak{f}}_{\rm sw}^{n})_{-},F_{\rm ext}^{n}\big{)}.$

The equations (29)-(35) furnish an induction relation that allows to compute
$U^{n+1}$ and $\Theta^{n+1}$ in terms of $U^{n}$ and $\Theta^{n}$. It need of
course to be initiated with initial data that are taken of the form

(36) $U_{i}^{0}=(\zeta^{\rm in}_{i},q^{\rm in})_{i}^{\rm
T}\qquad(i\in\big{(}{\mathbb{Z}}\backslash\\{0\\}\big{)}\cup\\{0^{-},0^{+}\\}),$

with $\zeta^{\rm in}$ and $q^{\rm in}$ describing the initial wave field in
the exterior domain, and

(37) $\Theta^{0}=\big{(}\langle q_{\rm i}\rangle^{\rm
in},\delta^{(1)},\underline{\zeta}_{+}^{(1)},\underline{\zeta}_{-}^{(1)},\delta^{(0)},\underline{\zeta}_{+}^{(0)},\underline{\zeta}_{-}^{(0)}\big{)}^{\rm
T}$

satisfies the discrete version of the compatibility conditions of Proposition
2.2, namely,

(38) $\langle q^{\rm in}\rangle=\langle q_{\rm i}\rangle^{\rm
in},\qquad\llbracket q^{\rm
in}\rrbracket=-2\ell\delta^{(1)},\qquad\underline{\zeta}_{\pm}^{(0)}=\zeta_{\pm}^{\rm
in},\qquad\underline{\zeta}^{(1)}_{\pm}=-(\partial_{x}q^{\rm in})_{\pm}.$

###  3.3. A second order scheme

We propose here an adaptation of the MacCormack scheme for the conservation
laws with nonlocal flux (24), coupled with a second-order Heun integration
scheme for the system of 7 first-order ODEs (27). Both are predictor-corrector
schemes. We use the same notations as in the previous subsection and can
decompose the scheme into four main steps:

\- Prediction step for the MacCormack scheme. This reads

(39) $\frac{U_{i}^{n,*}-U_{i}^{n}}{\Delta t}+\frac{1}{\Delta
x}\big{(}{\mathfrak{F}}^{n}_{\kappa,i}-{\mathfrak{F}}^{n}_{\kappa,i-1}\big{)}={\mathcal{S}}_{+}^{n}{\mathfrak{b}}_{i},\qquad
i>1,\quad n\geq 0,$

with $U_{i}^{n,*}=(\zeta_{i}^{n,*},q_{i}^{n,*})^{\rm T}$. We use a symmetric
scheme with respect to $x=0$ so that, for negative values of $i$, we use a
forward rather than backward derivative for the flux,

(40) $\frac{U_{i}^{n,*}-U_{i}^{n}}{\Delta t}+\frac{1}{\Delta
x}\big{(}{\mathfrak{F}}^{n}_{\kappa,i+1}-{\mathfrak{F}}^{n}_{\kappa,i}\big{)}={\mathcal{S}}_{-}^{n}{\mathfrak{b}}_{i},\qquad
i<-1,\quad n\geq 0,$

For $i=1$ and $i=-1$ it reads

(41) $\displaystyle\frac{U_{1}^{n,*}-U_{1}^{n}}{\Delta t}+\frac{1}{\Delta
x}\big{(}{\mathfrak{F}}^{n}_{\kappa,1}-{\mathfrak{F}}^{n}_{\kappa,0^{+}}\big{)}$
$\displaystyle={\mathcal{S}}_{+}^{n}{\mathfrak{b}}_{1},$ (42)
$\displaystyle\frac{U_{-1}^{n,*}-U_{-1}^{n}}{\Delta t}+\frac{1}{\Delta
x}\big{(}{\mathfrak{F}}^{n}_{\kappa,0^{-}}-{\mathfrak{F}}^{n}_{\kappa,-1}\big{)}$
$\displaystyle={\mathcal{S}}_{-}^{n}{\mathfrak{b}}_{-1},$

for $n\geq 0$ and with ${\mathfrak{F}}_{\kappa,0^{\pm}}^{n}$ as in (33).

\- Prediction step for the Heun scheme. This step is similar to a first-order
explicit Euler scheme,

(43) $\frac{\Theta^{n,*}-\Theta^{n}}{\Delta
t}={\mathcal{G}}\big{(}\Theta^{n},(R_{1}{\mathfrak{f}}_{\rm
sw}^{n})_{+},(R_{1}{\mathfrak{f}}_{\rm sw}^{n})_{-},F_{\rm ext}^{n}\big{)}.$

\- Corrector step for the MacCormack scheme. With the quantities computed in
the previous steps, we define

(44) ${\mathfrak{f}}_{\rm sw}^{n,*}={\mathfrak{f}}_{\rm
sw}(U^{n,*})\quad\mbox{ and }\quad
q^{n,*}_{0^{\pm}}=\Theta^{n,*}_{1}\mp\ell\Theta^{n,*}_{2}$

as well as an intermediate non-local flux and an intermediate source term,

(45) $\displaystyle{\mathfrak{F}}_{\kappa,i}^{n,*}$
$\displaystyle=\big{(}q_{i}^{n,*},(R_{1}{\mathfrak{f}}_{\rm
sw}^{n,*})_{i}\big{)}^{\rm T}\qquad|i|\geq 1,\quad n\geq 0,$ (46)
$\displaystyle{\mathcal{S}}_{\pm}^{n,*}$
$\displaystyle={\mathcal{S}}_{\pm}(\Theta^{n,*},(R_{1}{\mathfrak{f}}_{\rm
sw}^{n,*})_{+},(R_{1}{\mathfrak{f}}_{\rm sw}^{n,*})_{-},F_{\rm ext}^{n})\quad
n\geq 0.$

The correction step for the MacCormack scheme then reads

(47) $\frac{U_{i}^{n+1}-U_{i}^{n}}{\Delta
t}+\frac{{\mathfrak{F}}^{n}_{\kappa,i}-{\mathfrak{F}}^{n}_{\kappa,i-1}+{\mathfrak{F}}^{n,*}_{\kappa,i+1}-{\mathfrak{F}}^{n,*}_{\kappa,i}}{2\Delta
x}=\frac{{\mathcal{S}}_{+}^{n}+{\mathcal{S}}_{+}^{n,*}}{2}{\mathfrak{b}}_{i}\qquad
i\geq 1,$

for $n\geq 0$. Here again, we take a symmetric scheme so that for $i\leq-1$,
we take a forward difference of ${\mathfrak{F}}^{n}$ and a backward difference
of ${\mathfrak{F}}^{n,*}$,

(48) $\frac{U_{i}^{n+1}-U_{i}^{n}}{\Delta
t}+\frac{{\mathfrak{F}}^{n}_{\kappa,i+1}-{\mathfrak{F}}^{n}_{\kappa,i}+{\mathfrak{F}}^{n,*}_{\kappa,i}-{\mathfrak{F}}^{n,*}_{\kappa,i-1}}{2\Delta
x}=\frac{{\mathcal{S}}_{-}^{n}+{\mathcal{S}}_{-}^{n,*}}{2}{\mathfrak{b}}_{i}\qquad
i\leq-1;$

in particular, there is no need to defined boundary values
${\mathfrak{F}}_{\kappa,0^{\pm}}^{n,*}$, of the intermediate flux.

-Correction step for the Heun scheme. This reads, for $n\geq 0$,

(49) $\frac{\Theta^{n+1}-\Theta^{n}}{\Delta
t}=\frac{{\mathcal{G}}\big{(}\Theta^{n},(R_{1}{\mathfrak{f}}_{\rm
sw}^{n})_{\pm},F_{\rm
ext}^{n}\big{)}+{\mathcal{G}}\big{(}\Theta^{n,*},(R_{1}{\mathfrak{f}}_{\rm
sw}^{n,*})_{\pm},F_{\rm ext}^{n+1}\big{)}}{2}.$

The initial data have the same form as for the first order scheme described in
the previous subsection.

##  4\. Numerical simulations

We have seen in §2.1 that the wave-structure interaction problem under
consideration in this paper can be reduced to a transmission problem
potentially coupled to two forced ODEs for the vertical displacement $\delta$
of the object and the mean discharge $\langle q_{\rm i}\rangle$ under the
object.

We first consider in §4.1 a situation where this coupling is absent. This
corresponds to the case where a wave is generated in a wave tank by moving the
object vertically with a prescribed motion. This example is of particular
interest since it provides an efficient way to generate waves for the
Boussinesq equations at the entrance of a numerical domain if we have at our
disposal time series of the horizontal discharge at the boundary, hereby
extending the result of [26] where data on the surface elevation were used.

We then consider in §4.2 the return to equilibrium problem (also called decay
test or drop test by engineers) which consists in releasing an object from an
out of equilibrium position and to observe its oscillations. These examples
involve the coupling of the transmission problem with the ODE on $\delta$. In
the linear case, we are able to derive exact explicit solutions that we
compute to check the numerical convergence of our scheme; the nonlinear case
is then investigated and the importance of the dispersive effects pointed out
by comparing with simulations based on the nonlinear shallow water equations
instead of the Boussinesq system.

We then investigate in §4.3 a configuration where the transmission problem is
coupled to the interior discharge $\langle q_{\rm i}\rangle$, namely, the
interaction of waves with a fixed partially immersed object. Here again, we
derive an explicit exact solution in the linear case that we use to validate
that this coupling is also of second order. The nonlinear case is then
considered.

Finally, a configuration involving the most general coupling (with both
$\delta$ and $\langle q_{\rm i}\rangle$ is considered in §4.4; it consists in
the interaction of a solitary wave with an object freely floating in the
vertical direction.

###  4.1. Wave generation

The first physical configuration we consider consists in creating waves in a
fluid initially at rest by moving up and down a partially immersed object. By
symmetry, it is enough to consider the waves in the right component
${\mathcal{E}}^{+}=(\ell,\infty)$ of the fluid domain. As shown in §A.3 of
Appendix A, the mathematical formulation of this problem is a particular case
of the following initial boundary value problem with boundary condition on the
discharge $q$, namely,

(50) $\begin{cases}\partial_{t}\zeta+\partial_{x}q=0,\\\
(1-\kappa^{2}\partial_{x}^{2})\partial_{t}q+\partial_{x}{\mathfrak{f}}_{\rm
sw}=0,\end{cases}t>0,\quad x>\ell$

with ${\mathfrak{f}}_{\rm sw}$ as in (1) and with boundary condition

(51) $q_{|_{x=\ell}}(t)=g(t),\qquad t>0$

and initial condition

(52) $(\zeta,q)_{|_{t=0}}(x)=(\zeta^{\rm in},q^{\rm in})(x),\qquad x>\ell,$

and where $g$, $\zeta^{\rm in}$ and $q^{\rm in}$ are some given functions
satisfying the compatibility condition

(53) $q^{\rm in}(x=0)=g(t=0),$

which is obviously necessary to obtain solutions that are continuous at the
origin in time and space. This problem is somehow symmetric to the one
considered in [26] where a boundary condition on $\zeta$ rather than $q$ was
considered and where a first order scheme was proposed.

######  Remark 4.1.

For the wave generation problem, one has $(\zeta^{\rm in},q^{\rm in})=(0,0)$
and $g(t)=-\ell\dot{\delta}_{\rm forced}$, where $\delta_{\rm forced}$ is the
prescribed vertical displacement of the center of mass of the object.

Contrary to the other physical configurations we consider in this article, the
wave generation problem (or more generally, the initial boundary value problem
(50)-(53)) does not require the resolution of an ODE to determine the boundary
data on the discharge. The formulation as an initial value problem given in
Proposition 2.2 then reduces to

(54) $\begin{cases}\partial_{t}\zeta+\partial_{x}q=0,\\\
\partial_{t}q+\partial_{x}R_{1}{\mathfrak{f}}_{\rm
sw}={\mathcal{S}}_{+}(t)\exp(-\frac{x-\ell}{\kappa})\end{cases}\quad\mbox{
with }\quad{\mathcal{S}}_{+}(t)=\dot{g}(t),$

for $x>\ell$ and with initial condition (52) satisfying (53).
The numerical scheme presented in §3.2 and §3.3 can be simplified by skipping
the second and fourth step related to the Heun scheme, and by taking simply
for the first-order scheme:

${\mathcal{S}}^{n}=\frac{g^{n+1}-g^{n}}{\Delta t},$

and for the second-order scheme:

${\mathcal{S}}^{n}=\frac{g^{n+1}-g^{n-1}}{2\Delta t}\quad\mbox{ and
}\quad{\mathcal{S}}^{n,*}=\frac{g^{n+2}-g^{n}}{2\Delta t}.$

The wave generation problem gives us the opportunity to validate our numerical
code with a nonlinear case. The Boussinesq-Abbott equations admits solitary
waves solutions of the form

$(\zeta,q)(t,x)=\big{(}\zeta_{c}(x-x_{0}-ct),c\zeta_{c}(x-x_{0}-ct)\big{)},$

where $c>0$, $x_{0}\in{\mathbb{R}}$ and $\zeta_{c}$ is a smooth, even and
fastly decaying function. These solutions can be used to test the precision of
the code. For the Boussinesq-Abbott equations, there is no explicit formula
for $\zeta_{c}$ and it is determined by the resolution of a nonlinear second
order ODE, namely,

(55)
$c^{2}\kappa^{2}\zeta_{c}^{\prime\prime}-c^{2}\frac{\zeta}{1+\varepsilon\zeta}+\frac{\varepsilon^{2}\zeta^{2}+2\varepsilon\zeta}{2\varepsilon}=0,$

with

$c^{2}=\frac{\varepsilon}{6}\frac{3\zeta_{\rm max}^{2}+\varepsilon\zeta_{\rm
max}^{3}}{\zeta_{\rm max}-\frac{\ln(1+\varepsilon\zeta_{\rm
max})}{\varepsilon}}$

(see for instance [26] for more details on the computations); these formula
furnish a family of solitary waves parametrized by their maximal amplitude
$\zeta_{\rm max}$. Solving the above ODE with a standard high precision ODE
solver provides us a solution to (50) that we use to assess the precision of
the numerical solution obtained with our numerical scheme for (54) with
discharge boundary data $g(t)=c\zeta_{c}(0-x_{0}-ct)$ and initial data
$(\zeta^{\rm in},q^{\rm
in})=(\zeta_{c}(x-x_{0}-0),c\zeta_{c}(x-x_{0}-0)\big{)}$.

######  Remark 4.2.

For very fine meshes, spurious oscillations may appear. These oscillations are
reminiscent of the oscillations that appear when using dispersive schemes
(such as the Lax-Wendroff or MacCormack schemes) to simulate shock waves.
Flux-limiters methods are typically used to control this phenomenon [27].
Here, these oscillations are created at the boundary, whose position is fixed,
and we use a very simple efficient method consisting in adding an artificial
viscosity on a finite number $n_{0}$ of cells near the boundary. More
precisely, in the right-component of the fluid domain (the left component is
treated symmetrically) we add the following term in the right-hand side of the
first component of (48),

$\nu\frac{\Delta x}{\Delta
t}\big{(}\zeta^{n}_{i+1}-2\zeta^{n}_{i}+\zeta^{n}_{i-1}\big{)}\qquad 1\leq
i\leq n_{0}$

with $\nu>0$ a fixed coefficient, that we take equal to $2.136$. This
corresponds to an artificial viscosity $\nu\frac{(\Delta x)^{3}}{\Delta
t}\partial_{x}^{2}\zeta$; for a fixed ration $\Delta x/\Delta t$, this
viscosity is of order $2$ and therefore does not alter the overall second
order of the MacCormack scheme.

We choose for this test $\zeta_{{\rm max}}=1$ and
$3\kappa^{2}=\varepsilon=0.3$. Once $c$ is computed, we solve the differential
equation (55) with a high order numerical method in order to obtain our
reference solution. The size of the computational domain is $L=6$. The space
step is computed as $\Delta_{x}=L/N$, with $N=200,240,300,400$. We take a
constant time step $\Delta_{t}=0.8\,\Delta_{x}$. The maximum of the soliton is
initially located on the left of the computational domain, at $x=-L/2$, so
that the initial datum in the small domain is almost zero, and then the
soliton propagates inside it. The numerical results at final time $T_{f}=20$
are presented for both schemes on Figures 3 and 4, showing respectively a
first-order and a second-order convergence.

These results correspond to a final time where the soliton has completely
entered the computational domain, so that the influence of the dispersive
boundary layer due to the generating condition at the left side of the domain
is nearly zero. However, if one would choose a final time where the soliton is
still entering the computational domain, then one would notice that the error
for the variable $\zeta$ is only first-order in the vicinity of the left
boundary of the computational domain, while it is still second-order for the
variable $q$. This first-order error is probably due to a lack of accuracy in
the numerical evaluation of the spatial derivative of $q$ near the left
boundary.

|
---|---
Figure 3. Solitary wave with generating condition on $q$, $L^{\infty}$ error for $3\kappa^{2}=0.3$ with first-order scheme, left: $\zeta$, right: $q$. |
---|---
Figure 4. Solitary wave with generating condition on $q$, $L^{\infty}$ error
for $3\kappa^{2}=0.3$ with second-order scheme, left: $\zeta$, right: $q$.

###  4.2. Return to equilibrium

We consider here the return to equilibrium problem (also called decay test),
which consists in dropping the floating object from an out of equilibrium
position and to let it oscillate vertically and stabilize towards its
equilibrium position. This is a problem of practical importance because it is
used by engineers to characterize some buoyancy properties of the solid, and
theoretically because it leads to simpler equations than the general wave-
structure equations. For instance, in the nonlinear non dispersive case
($\varepsilon\neq 0$, $\kappa=0$), it is possible to show that the
dimensionless vertical displacement $\varepsilon\delta$ of the solid with
respect to its equilibrium position is fully described by a second order
nonlinear scalar ODE [24, 4] and that in the linear dispersive case
($\varepsilon=0$, $\kappa\neq 0$) it is governed by a second order linear
integro-differential equation [4]. In the nondispersive case, similar
equations have also been derived in the presence of viscosity [28] in the
linear case, as well as in the $2D$ radial and partially linear case [8]. In
the presence of nonlinear and dispersive effects ($\varepsilon\neq 0$,
$\kappa\neq 0$), it does not seem possible to derive such a simple equation
for the motion of the solid and the wave-structure equations must therefore be
solved.

As for the wave generation problem, there is a symmetry in this problem which
allows to consider only the right part ${\mathcal{E}}^{+}$ of the fluid domain
and the governing wave-structure interaction equations reduce to an initial
boundary value problem of the form (50) with $g=-\ell\dot{\delta}$. The
difference is that the vertical displacement $\delta$ is no longer a given
function $\delta$ but is found through the resolution of Newton’s equation
(see (7) for its general expression). Since this equation involves the trace
of $\zeta$ at the contact point, we have to work with the augmented
formulation provided by Proposition 2.2. Since in this particular case, one
has $F_{\rm ext}\equiv 0$, $\langle q_{\rm i}\rangle\equiv 0$ and
$\zeta_{+}=\zeta_{-}$, the $7$-dimensional ODE on $\Theta$ can be simplified
into a simpler $4$-dimensional ODE (see §A.3 in Appendix A). The interest of
this test case is that, since the interior discharge identically vanishes, it
allows us to investigate specifically the coupling between the waves and the
vertical displacement of the object. We first consider in §4.2.1 the linear
case for which explicit solutions exist and can be used to investigate the
precision of the code, and then show in §4.2.2 some simulations in the
nonlinear case.

####  4.2.1. Convergence error in the linear case

We first consider the linear case ($\varepsilon=0$) since in the case, it was
shown in [4] that the evolution of $\delta$ can be found by solving a linear
second order integro-differential equation, namely,

$\big{(}\tau_{\kappa}(0)^{2}+\ell\kappa\big{)}\ddot{\delta}+\ell{\mathcal{K}}^{1}_{\kappa}*\dot{\delta}+\delta=0,$

with initial conditions $\delta(0)=\delta_{0}$ and $\dot{\delta}(0)=0$ and
where $\tau_{\kappa}(\cdot)$ is defined in Appendix A and the kernel
${\mathcal{K}}_{\kappa}^{1}$ is given in terms of the first Bessel function
$J_{1}$ by the relation

${\mathcal{K}}^{1}_{\kappa}(t)=\frac{1}{t}J_{1}\big{(}\frac{t}{\kappa}\big{)}.$

The solution of this integro-differential equation is given explicitly by
taking the Laplace transform (denoted with a hat),

(56)
$\widehat{\delta}(s)=\frac{\tau_{\kappa}(0)^{2}s+\ell\sqrt{1+\kappa^{2}s^{2}}}{\tau_{\kappa}(0)^{2}s^{2}+s\ell\sqrt{1+\kappa^{2}s^{2}}+1}\delta_{0},\qquad
s\in{\mathbb{C}}_{0},$

where ${\mathbb{C}}_{0}$ is the half-plane of complex numbers $s$ such that
$\Re s>0$. The vertical displacement deduced from the exact formula (56), and
denoted $\delta_{\rm exact}$ is compared with the surface elevation $\delta$
found by solving the wave-structure equations using the numerical schemes
presented in §3. In order to discard possible numerical errors in the
computation of the inverse Laplace transform one has to apply to (56) two
different inversion methods (the Euler and Talbot methods [1]); we impose that
they match up to $10^{-4}$ terms to consider the solution provided as relevant
to be considered as an exact solution for our convergence studies.

In our numerical tests we chose $3\kappa^{2}=0.3$ or $3\kappa^{2}=0.1$,
$h_{eq}=1-0.3$, $l=4$ and the size of the computational domain $L=30$. The
space steps $\Delta_{x}$ were computed as $\Delta_{x}=(L-l)/(N+1)$, with
$N=300,400,500,600$ for the first-order scheme and $N=60,120,240,320$ for the
second-order scheme. The time step was computed as $\Delta_{t}=0.9\Delta_{x}$.
The numerical results at final time $T_{f}=15$ computed with the first-order
and the second-order schemes show respectively a first-order convergence, see
Figure 5 and a second-order convergence, see Figure 6. On Figure 7 we compare
the numerical results for the two schemes for $3\kappa^{2}=0.3$ and $N=100$,
showing evidence that it is advantageous to use the second-order scheme.

|
---|---
Figure 5. Return to equilibrium, linear case: convergence results for $\delta$ with the first-order scheme, $3\kappa^{2}=0.1$ (left) and $3\kappa^{2}=0.3$ (right). |
---|---
Figure 6. Return to equilibrium, linear case: convergence results for $\delta$ with the second-order-order scheme, $3\kappa^{2}=0.1$ (left) and $3\kappa^{2}=0.3$ (right). |
---|---
Figure 7. Return to equilibrium, linear case: temporal evolution of $\delta$,
comparison between the first-order scheme and the second-order scheme,
$3\kappa^{2}=0.3$ and $N=100$.

####  4.2.2. The nonlinear case

We do not have any exact solution to compare with in the nonlinear case, and
we therefore use mesh-convergence to study the order of our schemes. The
reference solution is computed with a very refined mesh: $N=2400$. We chose
$\epsilon=0.3$, $3\kappa^{2}=0.3$ or $3\kappa^{2}=0.1$, $h_{eq}=1-0.3$, $l=4$
and the size of the computational domain $L=30$. The space steps $\Delta_{x}$
were computed as $\Delta_{x}=(L-l)/(N+1)$, with $N=160,200,240,300,400$ for
the first-order scheme and $N=120,160,200$ for the second-order scheme. The
meshes are defined so that the points of the coarse meshes always coincide
with the points of the very refined mesh of the reference solution. The time
step was computed as $\Delta_{t}=0.7\Delta_{x}$. The numerical results at
final time $T_{f}=15$ computed with the first-order and the second-order
schemes show respectively a first-order convergence, see Figure 8 and a
second-order convergence, see Figure 9.

We perform another test to study qualitatively the influence of the dispersion
on the nonlinear decay test. We compare the trajectories obtained for
different values of $\kappa$ with the trajectory obtained in the non
dispersive case ($\kappa=0$). In the latter case, it was shown in [4] that the
evolution of $\delta$ can be found, under some smallness assumptions, by
solving a second order differential equation of the form

(57)
${{\tau_{0}(\varepsilon\delta)^{2}}}\ddot{\delta}+\ell\dot{\delta}+\delta+\varepsilon
B(\varepsilon\delta)\dot{\delta}^{2}=0,$

where $B(\cdot)$ is a smooth function whose exact expression can be found in
Corollary 4.3 of [4]. This nondispersive solution is used as reference to
illustrate the contribution of the dispersive terms. On Figure 10 we compare
this solution with the numerical results for $\varepsilon=0.3$,
$h_{eq}=1-0.3$, $l=0.25$ and various values of $3\kappa^{2}$ ranging from
$0.05$ to $1$.

|
---|---
Figure 8. Return to equilibrium, non-linear case: mesh-convergence results for $\delta$ with the first-order scheme, $3\kappa^{2}=0.1$ (left) and $3\kappa^{2}=0.3$ (right). |
---|---
Figure 9. Return to equilibrium, non-linear case: mesh-convergence results for
$\delta$ with the second-order-order scheme, $3\kappa^{2}=0.1$ (left) and
$3\kappa^{2}=0.3$ (right).
---
Figure 10. Return to equilibrium, non-linear case: temporal evolution of
$\delta$ for $\epsilon=0.3$, comparison between results with the second-order
scheme obtained with $N=400$, different values of $\mu=3\kappa^{2}$, and the
exact solution in the non dispersive case.

###  4.3. Waves interacting with a fixed object

We consider here waves that are interacting with a fixed partially immersed
object. This is a particular example of prescribed motion ($\delta\equiv 0$),
but contrary to the wave generation problem, the waves are not supposed to be
symmetric with respect to the central axis $\\{x=0\\}$. It follows that the
interior discharge $\langle q_{\rm i}\rangle$ does not vanish identically and
that it must be found by solving the ODE (6). From a mathematical point of
view, this physical configuration is somehow symmetric to the return to
equilibrium problem in the sense that it allows one to focus on the coupling
of the fluid equation with the dynamic of the interior discharge $\langle
q_{\rm i}\rangle$ (since $\delta\equiv 0$), while for the return to
equilibrium problem the coupling was only with the vertical displacement
$\delta$ (since in that case $\langle q_{\rm i}\rangle\equiv 0$). In this
case, the $7$-dimensional ODE on $\Theta$ of the augmented formulation given
in Proposition 2.2 can be reduced to a $5$-dimensional ODE, as explained in
§A.3 of Appendix A. We first study in §4.3.1 the linear case for which we
exhibit a family of explicit solutions that we use to validate our code; the
nonlinear case is then considered in §4.3.2

####  4.3.1. Convergence error in the linear case

In order to investigate the ability of our scheme to correctly describe the
coupling of the Boussinesq-Abbott equation with the average interior discharge
$\langle q_{\rm i}\rangle$ we exhibit an explicit solution of the equations in
the linear case ($\varepsilon=0$). In that case, the wave-structure equations
(4)-(7) take the form

(58) $\begin{cases}\partial_{t}\zeta+\partial_{x}q=0,\\\
(1-\kappa^{2}\partial_{x}^{2})\partial_{t}q+\partial_{x}\zeta=0\end{cases}\quad\mbox{
in }\quad{\mathcal{E}}^{\pm}$

with transmission conditions

(59) $\llbracket q\rrbracket=0\quad\mbox{ and }\quad\langle q\rangle=\langle
q_{\rm i}\rangle$

and where $\langle q_{\rm i}\rangle$ solves the forced ODE

(60) $\alpha(0)\frac{d}{dt}\langle q_{\rm
i}\rangle=-\frac{1}{2\ell}\llbracket\zeta+\kappa^{2}\partial_{t}^{2}\zeta\rrbracket.$

A family of exact solutions which are periodic in time is given in the
following proposition (which can be checked with basic computations omitted
here).

######  Proposition 4.1.

Let $k\neq 0$ and $\omega\neq 0$ satisfy the the dispersion relation

$\omega^{2}=\frac{k^{2}}{1+\kappa^{2}k^{2}}.$

For all
$(\zeta_{+}^{c},\zeta_{-}^{c},\zeta^{s},q_{+}^{s},q_{-}^{s},q^{c})\in{\mathbb{R}}^{6}$,
the functions $(\zeta,q)$ defined in ${\mathcal{E}}^{\pm}$ by

$\begin{cases}\zeta(x\mp\ell,t)=&\frac{k}{2}\big{[}(\zeta_{\pm}^{c}+q^{c})\cos(kx-\omega
t)+(\zeta_{\pm}^{c}-q^{c})\cos(kx+\omega t)\\\
&+(\zeta^{s}+q_{\pm}^{s})\sin(kx-\omega
t)+(\zeta^{s}-q_{\pm}^{s})\sin(kx+\omega t)\big{]},\\\
q(x\mp\ell,t)=&\frac{\omega}{2}\big{[}(\zeta_{\pm}^{c}+q^{c})\cos(kx-\omega
t)-(\zeta_{\pm}^{c}-q^{c})\cos(kx+\omega t)\\\
&+(\zeta^{s}+q_{\pm}^{s})\sin(kx-\omega
t)-(\zeta^{s}-q_{\pm}^{s})\sin(kx+\omega t)\big{]},\end{cases}$

solve (58)-(59) with initial data

$\begin{cases}\zeta^{\text{\rm
in}}(x\pm\ell)=k\zeta_{\pm}^{c}\cos(kx)+k\zeta^{s}\sin(kx),&x\in\mathcal{E}_{\pm}\\\
q^{\text{\rm in}}(x\pm\ell)=\omega q^{c}\cos(kx)+\omega
q_{\pm}^{s}\sin(kx),&x\in\mathcal{E}_{\pm}\end{cases}$

and with

$\langle q_{\rm i}\rangle(t)=\omega\big{[}q^{c}\cos(\omega
t)-\zeta^{s}\sin(\omega t)\big{]}.$

If moreover

$q^{c}=-\frac{1}{2\ell\alpha(0)k}(q_{+}^{s}-q_{-}^{s}),\quad\mbox{ and
}\quad\zeta^{s}=\frac{1}{2\ell\alpha(0)k}(\zeta_{+}^{c}-\zeta_{-}^{c})$

then (60) is also satisfied, with initial data $\langle q_{\rm i}\rangle^{\rm
in}=\omega q^{c}$.

In our numerical tests we chose $3\kappa^{2}=0.3$ or $3\kappa^{2}=0.1$,
$h_{eq}=1-0.2$, $l=1$, $k=2$ and the size of the computational domain $L=10$.
The space steps $\Delta_{x}$ were computed as $\Delta_{x}=(L-l)/(N+1)$, with
$N=200,240,300,360,400$ for the second-order scheme. The time step was
computed as $\Delta_{t}=0.9\Delta_{x}$. To impose the exact solution on both
left and right outer boundaries we use the wave generation method described in
§4.1. The numerical results at final time $T_{f}=1$ show a second-order
convergence, see Figure 11. On Figure 12 one can observe the shape of this
exact solution, computed with $N=400$ points.

|
---|---
Figure 11. Interaction with a fixed object, linear case: convergence results for $\langle q_{\rm i}\rangle$ with the second-order scheme, $3\kappa^{2}=0.1$ (left) and $3\kappa^{2}=0.3$ (right). |
---|---
Figure 12. Interaction with a fixed object, linear case: profile of solutions,
$3\kappa^{2}=0.1$ (left) and $3\kappa^{2}=0.3$ (right).

####  4.3.2. The nonlinear case

In the absence of explicit solution in the nonlinear case, we use mesh-
convergence to study the precision of our schemes. In this test the initial
condition is the solitary wave described in §4.1, with $\zeta_{\rm max}=0.2$,
centered at $x=-15$, at the left side of the fixed object. The size of the
computational domain is $L=30$. For this test the reference solution is
computed with a very refined mesh: $N=2400$. We chose $\epsilon=0.3$,
$3\kappa^{2}=0.3$ or $3\kappa^{2}=0.1$, $h_{eq}=1-0.3$, $l=4$. The space steps
$\Delta_{x}$ were computed as $\Delta_{x}=(L-l)/(N+1)$, with
$N=160,200,240,300,400$ for the first-order scheme and $N=100,120,160$ for the
second-order scheme. The meshes are defined so that the points of the coarse
meshes always coincide with the points of the very refined mesh of the
reference solution. The time step was computed as $\Delta_{t}=0.7\Delta_{x}$.
The numerical results at final time $T_{f}=20$ computed with the first-order
and the second-order schemes show respectively a first-order convergence for
$\langle q_{\rm i}\rangle$, see Figure 13 and a second-order convergence, see
Figure 14. On Figure 15 one can observe the shape of the numerical solution,
computed with $N=400$ points.

|
---|---
Figure 13. Interaction with a fixed object, non-linear case: mesh-convergence results for $\langle q_{\rm i}\rangle$ with the first-order scheme, $3\kappa^{2}=0.1$ (left) and $3\kappa^{2}=0.3$ (right). |
---|---
Figure 14. Interaction with a fixed object, non-linear case: mesh-convergence results for $\langle q_{\rm i}\rangle$ with the second-order scheme, $3\kappa^{2}=0.1$ (left) and $3\kappa^{2}=0.3$ (right). |
---|---
Figure 15. Interaction with a fixed object, non-linear case: profile of
solutions, $3\kappa^{2}=0.1$ (left) and $3\kappa^{2}=0.3$ (right), computed
with the second-order scheme and $N=400$.

###  4.4. Waves interacting with a freely floating object

In the general case, when the object moves under the influence of the waves
and the waves are in return modified by the presence of the object, one has to
consider the full augmented system presented in Proposition 2.2.

We use exactly the same test as in §4.3.2, the only difference being that the
object is now allowed to move vertically under the action of the waves, so
that we also have to study the convergence of $\delta$. The space steps
$\Delta_{x}$ were computed as $\Delta_{x}=(L-l)/(N+1)$, with $N=240,300,400$
for the first-order scheme, $N=120,160,200,240$ for the second-order scheme
and $3\kappa^{2}=0.1$, and $N=80,100,120,160$ for the second-order scheme and
$3\kappa^{2}=0.3$. The numerical results at final time $T_{f}=20$ computed
with the first-order and the second-order schemes show respectively a first-
order convergence, see Figures 16, 17 and 18, and a second-order convergence,
see Figures 19, 20 and 21. On Figure 22 one can observe the shape of the
numerical solution, computed with $N=400$ points. A comparison between Figure
22 and Figure 15 shows that the profiles of the reflected and transmitted
waves differ. In particular, when the object is allowed to move, the reflected
and transmitted wave are preceded by a depression wave that is not present
when the object if fixed.

|
---|---
Figure 16. Interaction with a freely floating object, non-linear case: mesh-convergence results for $\delta$ with the first-order scheme, $3\kappa^{2}=0.1$ (left) and $3\kappa^{2}=0.3$ (right). |
---|---
Figure 17. Interaction with a freely floating object, non-linear case: mesh-convergence results for $\langle q_{\rm i}\rangle$ with the first-order scheme, $3\kappa^{2}=0.1$ (left) and $3\kappa^{2}=0.3$ (right). |
---|---
Figure 18. Interaction with a freely floating object, non-linear case: mesh-convergence results for $\zeta_{+}$ with the first-order scheme, $3\kappa^{2}=0.1$ (left) and $3\kappa^{2}=0.3$ (right). |
---|---
Figure 19. Interaction with a freely floating object, non-linear case: mesh-convergence results for $\delta$ with the second-order scheme, $3\kappa^{2}=0.1$ (left) and $3\kappa^{2}=0.3$ (right). |
---|---
Figure 20. Interaction with a freely floating object, non-linear case: mesh-convergence results for $\langle q_{\rm i}\rangle$ with the second-order scheme, $3\kappa^{2}=0.1$ (left) and $3\kappa^{2}=0.3$ (right). |
---|---
Figure 21. Interaction with a freely floating object, non-linear case: mesh-convergence results for $\zeta_{+}$ with the second-order scheme, $3\kappa^{2}=0.1$ (left) and $3\kappa^{2}=0.3$ (right). |
---|---
Figure 22. Interaction with a freely floating object, non-linear case: profile
of solutions, $3\kappa^{2}=0.1$ (left) and $3\kappa^{2}=0.3$ (right), computed
with the second-order scheme and $N=400$.

##  Appendix A Exact expressions of the quantities involved in the wave-
structure models

We have shown in §2.3 that the wave-structure interaction equations can be
formulated as a transmission problem between the two components
${\mathcal{E}}^{-}=(-\infty,-\ell)$ and ${\mathcal{E}}^{+}=(\ell,\infty)$ of
the fluid domain, with transmission conditions involving the vertical
displacement $\delta$ of the object and the average horizontal discharge
$\langle q_{\rm i}\rangle$ under the object; we have also shown that, in the
augmented formulation, $\delta$ and $\langle q_{\rm i}\rangle$ are found
through the resolution of the $7$-dimensional ODE (18) on
$\Theta:=\big{(}\langle q_{\rm
i}\rangle,\dot{\delta},\dot{\underline{\zeta}}{}_{+},\dot{\underline{\zeta}}{}_{-},\delta,{\underline{\zeta}{}_{+}},\underline{\zeta}{}_{-}\big{)}^{\rm
T}$ that we wrote in abstract form as

(61) $\frac{d}{dt}\Theta={\mathcal{G}}\big{(}\Theta,(R_{1}{\mathfrak{f}}_{\rm
sw})_{+},(R_{1}{\mathfrak{f}}_{\rm sw})_{-},F_{\rm ext}\big{)},$

where ${\mathcal{G}}$ is a smooth function of its arguments. The goal of this
section is to derive the explicit expression of the mapping ${\mathcal{G}}$
which is used for our numerical computations. We first provide in §A.1 the
explicit expression of various coefficients that appear in the wave-structure
equations and then derive in §A.2 the explicit expression of the mapping
${\mathcal{G}}$ in the most general case. We then point out the
simplifications that can be performed when the object is in fixed or forced
motion (in §A.3) and when the system has a symmetry with respect to the
vertical axis $\\{x=0\\}$ (in §A.4 ).

N.B. We recall that for the sake of simplicity, we assume throughout this
article that $h_{\rm eq}(x)$ is an even function.

###  A.1. Explicit expressions of some coefficients

We did not made explicit in the main text of the article most of the constants
that appear in the wave-structure interaction equations studied in this paper
and derived in [4] because they were not relevant for the mathematical and
numerical analysis of these equations. Of course, they are necessary for
realistic simulations of wave-structure interactions and we provide them here.
Let us first remind that the configuration under consideration is a floating
object with vertical sidewalls located at $x=\pm\ell$ and that can only move
in the vertical direction. In dimensionless variables, we denote by $h_{\rm
eq}(x)$ the water depth below the object at equilibrium and by
$\varepsilon\delta(t)$ the displacement of the object at time $t$ from its
equilibrium position, so that that the water depth under the object at time
$t$ is $h_{\rm eq}(x)+\varepsilon\delta(t)$. The dimensionless mass $m$ of the
object can be defined through Archimedes’ principle,

$m=\frac{1}{2\ell}\int_{-\ell}^{\ell}(1-h_{\rm eq})$

and the formulas below will also involve two scalar functions
$\alpha(\varepsilon\delta)$ and $\beta(\varepsilon\delta)$ defined as

$\displaystyle\alpha(\varepsilon\delta)=\frac{1}{2\ell}\int_{-\ell}^{\ell}\frac{1}{h_{\rm
eq}(x)+\varepsilon\delta}{\rm d}x\quad\mbox{ and
}\quad\beta(\varepsilon\delta)=\frac{1}{2}\frac{1}{2\ell}\int_{-\ell}^{\ell}\frac{x^{2}}{(h_{\rm
eq}(x)+\varepsilon\delta)^{2}}{\rm d}x;$

the quantities $\tau_{\kappa}(\varepsilon\delta)^{2}$ that appears in Newton’s
equation is given by

$\tau_{\kappa}(\varepsilon\delta)^{2}=3\kappa^{2}m+\frac{1}{2\ell}\int_{-\ell}^{\ell}\frac{x^{2}}{h_{\rm
eq}+\varepsilon\delta}{\rm d}x+\kappa^{2}\langle\frac{1}{h_{\rm
eq}+\varepsilon\delta}\rangle.$

###  A.2. The general case

As seen in (14), the first four components of $\Theta$ satisfy

(62)
${\mathcal{M}}[\varepsilon\delta,\varepsilon\underline{\zeta}{}_{\pm}]\frac{d}{dt}\begin{pmatrix}\displaystyle\langle
q_{\rm i}\rangle\\\ \displaystyle\dot{\delta}\\\
\displaystyle\dot{\underline{\zeta}}{}_{+}\\\
\displaystyle\dot{\underline{\zeta}}{}_{-}\end{pmatrix}+\begin{pmatrix}\displaystyle\frac{1}{2\ell}\llbracket\underline{\zeta}\rrbracket\\\
\displaystyle\delta-\langle\underline{\zeta}\rangle\\\
\displaystyle{\underline{\zeta}{}_{+}}\\\
\displaystyle\underline{\zeta}{}_{-}\end{pmatrix}=\varepsilon{\boldsymbol{\mathfrak{Q}}}[\varepsilon\delta,\varepsilon\underline{\zeta}{}_{\pm}](\langle
q_{\rm i}\rangle,\dot{\delta},\underline{\zeta}{}_{\pm})+\begin{pmatrix}0\\\
F_{\rm ext}\\\ \displaystyle({R}_{1}{\mathfrak{f}}_{\rm sw})_{+}\\\
\displaystyle({R}_{1}{\mathfrak{f}}_{\rm sw})_{-}\end{pmatrix}$

where
$\langle\underline{\zeta}\rangle:=\frac{1}{2}({\underline{\zeta}{}_{+}}+\underline{\zeta}{}_{-})$,
$\llbracket\underline{\zeta}\rrbracket:=\zeta_{+}-\underline{\zeta}{}_{-}$,
and ${\mathcal{M}}[\varepsilon\delta,\varepsilon\underline{\zeta}{}_{\pm}]$ is
the invertible matrix

(63)
${\mathcal{M}}[\varepsilon\delta,\varepsilon\underline{\zeta}{}_{\pm}]:=\left(\begin{array}[]{cc|cc}\alpha(\varepsilon\delta)&0&\frac{\kappa^{2}}{2\ell}\frac{1}{\underline{h}_{+}}&-\frac{\kappa^{2}}{2\ell}\frac{1}{\underline{h}_{-}}\\\
0&\tau_{\kappa}(\varepsilon\delta)^{2}&-\frac{1}{2}\frac{\kappa^{2}}{\underline{h}_{+}}&-\frac{1}{2}\frac{\kappa^{2}}{\underline{h}_{-}}\\\
\hline\cr-{\kappa}&\ell\kappa&\kappa^{2}&0\\\
\kappa&\ell\kappa&0&\kappa^{2}\end{array}\right),$

with $\underline{h}_{\pm}=1+\varepsilon\underline{\zeta}_{\pm}$; simple
computations also show that the quadratic term
${\boldsymbol{\mathfrak{Q}}}[\varepsilon\delta,\varepsilon\underline{\zeta}{}_{\pm}]$
is of the form

${\boldsymbol{\mathfrak{Q}}}[\varepsilon\delta,\varepsilon\underline{\zeta}{}_{\pm}]=\big{(}{\mathfrak{Q}}_{\rm
i}[\varepsilon\delta,\varepsilon\underline{\zeta}{}_{\pm}],{\mathfrak{Q}}_{\delta}[\varepsilon\delta,\varepsilon\underline{\zeta}{}_{\pm}],{\mathfrak{Q}}_{+}[\varepsilon\delta,\varepsilon\underline{\zeta}{}_{\pm}],{\mathfrak{Q}}_{-}[\varepsilon\delta,\varepsilon\underline{\zeta}{}_{\pm}]\big{)}^{\rm
T}$

with (writing simply ${\mathfrak{Q}}_{\rm i}={\mathfrak{Q}}_{\rm
i}[\varepsilon\delta,\varepsilon\underline{\zeta}{}_{\pm}]$, etc.)

$\displaystyle{\mathfrak{Q}}_{\rm i}(\langle q_{\rm
i}\rangle,\dot{\delta},\underline{\zeta}{}_{\pm})$
$\displaystyle=-\alpha^{\prime}(\varepsilon\delta)\langle q_{\rm
i}\rangle\dot{\delta}-\frac{1}{4\ell}\big{[}\Big{(}\frac{-\ell\dot{\delta}+\langle
q_{\rm
i}\rangle}{\underline{h}_{+}}\Big{)}^{2}-\Big{(}\frac{\ell\dot{\delta}+\langle
q_{\rm i}\rangle}{\underline{h}_{-}}\Big{)}^{2}\big{]},$
$\displaystyle{\mathfrak{Q}}_{\delta}(\langle q_{\rm
i}\rangle,\dot{\delta},\underline{\zeta}{}_{\pm})$
$\displaystyle=\beta(\varepsilon\delta)\dot{\delta}^{2}+\frac{1}{2}\alpha^{\prime}(\varepsilon\delta)\langle
q_{\rm i}\rangle^{2}+\frac{1}{4}\big{[}\Big{(}\frac{-\ell\dot{\delta}+\langle
q_{\rm
i}\rangle}{\underline{h}_{+}}\Big{)}^{2}+\Big{(}\frac{\ell\dot{\delta}+\langle
q_{\rm i}\rangle}{\underline{h}_{-}}\Big{)}^{2}\big{]},$
$\displaystyle{\mathfrak{Q}}_{+}(\langle q_{\rm
i}\rangle,\dot{\delta},\underline{\zeta}{}_{\pm})$
$\displaystyle=\frac{1}{2}{\underline{\zeta}{}_{+}}^{2}-\frac{1}{\underline{h}_{+}}\big{(}-\ell\dot{\delta}+\langle
q_{\rm i}\rangle\big{)}^{2},$ $\displaystyle{\mathfrak{Q}}_{-}(\langle q_{\rm
i}\rangle,\dot{\delta},\underline{\zeta}{}_{\pm})$
$\displaystyle=\frac{1}{2}\underline{\zeta}{}_{-}^{2}-\frac{1}{\underline{h}_{-}}\big{(}\ell\dot{\delta}+\langle
q_{\rm i}\rangle\big{)}^{2}.$

The matrix ${\mathcal{M}}$ is a $4\times 4$ matrix whose inverse is quite
complicated; we therefore transform it into a block-triangular matrix by
multiplying (62) by the matrix

$\left(\begin{array}[]{cc|cc}1&0&-\frac{1}{2\ell}\frac{1}{\underline{h}_{+}}&\frac{1}{2\ell}\frac{1}{\underline{h}_{-}}\\\
0&1&\frac{1}{2}\frac{1}{\underline{h}_{+}}&\frac{1}{2}\frac{1}{\underline{h}_{-}}\\\
\hline\cr 0&0&1&0\\\ 0&0&0&1\end{array}\right);$

the resulting equation takes the form

(64)
$\widetilde{\mathcal{M}}[\varepsilon\delta,\varepsilon\underline{\zeta}{}_{\pm}]\frac{d}{dt}\begin{pmatrix}\displaystyle\langle
q_{\rm i}\rangle\\\ \displaystyle\dot{\delta}\\\
\displaystyle\dot{\underline{\zeta}}{}_{+}\\\
\displaystyle\dot{\underline{\zeta}}{}_{-}\end{pmatrix}+\begin{pmatrix}\displaystyle
0\\\ \displaystyle\delta\\\ \displaystyle{\underline{\zeta}{}_{+}}\\\
\displaystyle\underline{\zeta}{}_{-}\end{pmatrix}=\varepsilon\widetilde{\boldsymbol{\mathfrak{Q}}}(\langle
q_{\rm
i}\rangle,\dot{\delta},\underline{\zeta}{}_{\pm})+\begin{pmatrix}\displaystyle-\frac{1}{2\ell}\llbracket\frac{1}{\underline{h}}R_{1}{\mathfrak{f}}_{\rm
sw}\rrbracket\\\
\displaystyle\langle\frac{1}{\underline{h}}R_{1}{\mathfrak{f}}_{\rm
sw}\rangle+F_{\rm ext}\\\ \displaystyle({R}_{1}{\mathfrak{f}}_{\rm sw})_{+}\\\
\displaystyle({R}_{1}{\mathfrak{f}}_{\rm sw})_{-}\end{pmatrix}$

where the matrix
$\widetilde{\mathcal{M}}[\varepsilon\delta,\varepsilon\underline{\zeta}{}_{\pm}]$
is block-triangular,

(65)
$\widetilde{\mathcal{M}}[\varepsilon\delta,\varepsilon\underline{\zeta}{}_{\pm}]:=\left(\begin{array}[]{cc|cc}\alpha(\varepsilon\delta)+\frac{\kappa}{\ell}\langle\frac{1}{\underline{h}}\rangle&-\frac{\kappa}{2}\llbracket\frac{1}{\underline{h}}\rrbracket&0&0\\\
-\frac{\kappa}{2}\llbracket\frac{1}{\underline{h}}\rrbracket&\tau_{\kappa}(\varepsilon\delta)^{2}+\kappa\ell\langle\frac{1}{\underline{h}}\rangle&0&0\\\
\hline\cr-{\kappa}&\ell\kappa&\kappa^{2}&0\\\
\kappa&\ell\kappa&0&\kappa^{2}\end{array}\right),$

and the components $\widetilde{\mathfrak{Q}}_{\rm i}$,
$\widetilde{\mathfrak{Q}}_{\delta}$ and $\widetilde{\mathfrak{Q}}_{\pm}$ are
given by

$\displaystyle\widetilde{\mathfrak{Q}}_{\rm i}(\langle q_{\rm i}\rangle,$
$\displaystyle\dot{\delta},\underline{\zeta}{}_{\pm})=-\frac{1}{4\ell}\llbracket\frac{1}{\underline{h}}\underline{\zeta}^{2}\rrbracket$
$\displaystyle+\frac{1}{4}\ell\llbracket\frac{1}{\underline{h}^{2}}\rrbracket\dot{\delta}^{2}+\frac{1}{4\ell}\llbracket\frac{1}{\underline{h}^{2}}\rrbracket\langle
q_{\rm
i}\rangle^{2}-\big{(}\alpha^{\prime}(\varepsilon\delta)+\langle\frac{1}{\underline{h}^{2}}\rangle\big{)}\dot{\delta}\langle
q_{\rm i}\rangle,$ $\displaystyle\widetilde{\mathfrak{Q}}_{\delta}(\langle
q_{\rm i}\rangle,$
$\displaystyle\dot{\delta},\underline{\zeta}{}_{\pm})=\frac{1}{2}\langle\frac{1}{\underline{h}}\underline{\zeta}^{2}\rangle$
$\displaystyle+\big{(}\beta(\varepsilon\delta)-\frac{1}{2}\ell^{2}\langle\frac{1}{\underline{h}^{2}}\rangle\big{)}\dot{\delta}^{2}+\frac{1}{2}\big{(}\alpha^{\prime}(\varepsilon\delta)-\langle\frac{1}{\underline{h}^{2}}\rangle\big{)}\langle
q_{\rm
i}\rangle^{2}+\frac{1}{2}\ell\llbracket\frac{1}{\underline{h}^{2}}\rrbracket\dot{\delta}\langle
q_{\rm i}\rangle,$

and

$\displaystyle\widetilde{\mathfrak{Q}}_{+}(\langle q_{\rm
i}\rangle,\dot{\delta},\underline{\zeta}{}_{\pm})$
$\displaystyle=-\frac{1}{2}{\underline{\zeta}{}_{+}}^{2}-\frac{1}{\underline{h}_{+}}\langle
q_{\rm
i}\rangle^{2}-\ell^{2}\frac{1}{\underline{h}_{+}}\dot{\delta}^{2}+2\ell\frac{1}{\underline{h}_{+}}\langle
q_{\rm i}\rangle\dot{\delta},$
$\displaystyle\widetilde{\mathfrak{Q}}_{-}(\langle q_{\rm
i}\rangle,\dot{\delta},\underline{\zeta}{}_{\pm})$
$\displaystyle=-\frac{1}{2}\underline{\zeta}{}_{-}^{2}-\frac{1}{\underline{h}_{-}}\langle
q_{\rm
i}\rangle^{2}-\ell^{2}\frac{1}{\underline{h}_{-}}\dot{\delta}^{2}-2\ell\frac{1}{\underline{h}_{-}}\langle
q_{\rm i}\rangle\dot{\delta}.$

Since by definition of ${\mathcal{G}}$ one has

(66) $\frac{d}{dt}\big{(}\langle q_{\rm
i}\rangle,\dot{\delta},\dot{\underline{\zeta}}{}_{+},\dot{\underline{\zeta}}{}_{-}\big{)}^{\rm
T}={\mathcal{G}}_{\rm I}(\Theta,(R_{1}{\mathfrak{f}}_{\rm sw})_{\pm},F_{\rm
ext})$

with the notation ${\mathcal{G}}_{\rm
I}:=\big{(}{\mathcal{G}}_{1},{\mathcal{G}}_{2},{\mathcal{G}}_{3},{\mathcal{G}}_{4}\big{)}^{\rm
T}$; we deduce from (64) that

$\displaystyle{\mathcal{G}}_{\rm I}$
$\displaystyle(\Theta,(R_{1}{\mathfrak{f}}_{\rm sw})_{\pm},F_{\rm ext})$
$\displaystyle=\widetilde{\mathcal{M}}[\varepsilon\delta,\varepsilon\underline{\zeta}{}_{\pm}]^{-1}\Big{[}-\begin{pmatrix}\displaystyle
0\\\ \displaystyle\delta\\\ \displaystyle{\underline{\zeta}{}_{+}}\\\
\displaystyle\underline{\zeta}{}_{-}\end{pmatrix}+\varepsilon\widetilde{\boldsymbol{\mathfrak{Q}}}(\langle
q_{\rm
i}\rangle,\dot{\delta},\underline{\zeta}{}_{\pm})+\begin{pmatrix}\displaystyle-\frac{1}{2\ell}\llbracket\frac{1}{\underline{h}}R_{1}{\mathfrak{f}}_{\rm
sw}\rrbracket\\\
\displaystyle\langle\frac{1}{\underline{h}}R_{1}{\mathfrak{f}}_{\rm
sw}\rangle+F_{\rm ext}\\\ \displaystyle({R}_{1}{\mathfrak{f}}_{\rm sw})_{+}\\\
\displaystyle({R}_{1}{\mathfrak{f}}_{\rm sw})_{-}\end{pmatrix}\Big{]},$

while

${\mathcal{G}}_{5}(\Theta)=\dot{\delta},\qquad{\mathcal{G}}_{6}(\Theta)=\dot{\underline{\zeta}{}_{+}},\qquad{\mathcal{G}}_{7}(\Theta)=\dot{}\underline{\zeta}{}_{-}.$

######  Remark A.1.

For the numerical computations, we use the explicit expression for the inverse
of the matrix
$\widetilde{\mathcal{M}}[\varepsilon\delta,\varepsilon\zeta_{\pm}]$, namely,

(67)
$\widetilde{\mathcal{M}}[\varepsilon\delta,\varepsilon\underline{\zeta}{}_{\pm}]^{\rm-1}=\left(\begin{array}[]{cc|cc}-\frac{4}{D}\big{(}\tau_{\kappa}(\varepsilon\delta)^{2}+\kappa\ell\langle\frac{1}{\underline{h}}\rangle\big{)}&-\frac{2\kappa}{D}\llbracket\frac{1}{\underline{h}}\rrbracket&0&0\\\
-\frac{2\kappa}{D}\llbracket\frac{1}{\underline{h}}\rrbracket&-\frac{4}{D}\big{(}\alpha(\varepsilon\delta)+\frac{\kappa}{\ell}\langle\frac{1}{\underline{h}}\rangle\big{)}&0&0\\\
\hline\cr-\frac{4}{\kappa
D}\big{(}\tau_{\kappa}(\varepsilon\delta)^{2}+\kappa\ell{\frac{1}{\underline{h}_{-}}}\big{)}&\frac{4}{\kappa
D}\big{(}\kappa{\frac{1}{\underline{h}_{-}}}+\ell\alpha(\varepsilon\delta)\big{)}&\frac{1}{\kappa^{2}}&0\\\
\frac{4}{\kappa
D}\big{(}\tau_{\kappa}(\varepsilon\delta)^{2}+\kappa\ell{\frac{1}{\underline{h}_{+}}}\big{)}&\frac{4}{\kappa
D}\big{(}\kappa{\frac{1}{\underline{h}_{+}}}+\ell\alpha(\varepsilon\delta)\big{)}&0&\frac{1}{\kappa^{2}}\end{array}\right)$

with

(68)
$D=-4\big{(}\alpha(\varepsilon\delta)+\frac{\kappa}{\ell}\langle\frac{1}{\underline{h}}\rangle\big{)}\times\big{(}\tau_{\kappa}(\varepsilon\delta)^{2}+\kappa\ell\langle\frac{1}{\underline{h}}\rangle\big{)}+\kappa^{2}\llbracket\frac{1}{\underline{h}}\rrbracket^{2}.$

###  A.3. The case of an object fixed or in forced motion

When the object is fixed or in forced motion, the position of the center of
mass is known and $\delta$ is therefore equal to some given function
$\delta_{\rm forced}$ ($\delta_{\rm forced}\equiv 0$ if the solid is fixed).
The ODE (61) can be reduced to an ODE on ${\mathbb{R}}^{5}$ instead of
${\mathbb{R}}^{7}$. The variable $\Theta$ now stands for
$\Theta:=\big{(}\langle q_{\rm
i}\rangle,\dot{\underline{\zeta}{}_{+}},\dot{}\underline{\zeta}{}_{-},{\underline{\zeta}{}_{+}},\underline{\zeta}{}_{-}\big{)}^{\rm
T}$ and (64) can be simplified into

$\displaystyle\widetilde{\mathcal{M}}_{\rm
forced}[t,\varepsilon\underline{\zeta}{}_{\pm}]\frac{d}{dt}\begin{pmatrix}\displaystyle\langle
q_{\rm i}\rangle\\\ \displaystyle\dot{\underline{\zeta}}{}_{+}\\\
\displaystyle\dot{\underline{\zeta}}{}_{-}\end{pmatrix}+\begin{pmatrix}\displaystyle
0\\\ \displaystyle{\underline{\zeta}{}_{+}}\\\
\displaystyle\underline{\zeta}{}_{-}\end{pmatrix}=$
$\displaystyle\varepsilon\widetilde{\boldsymbol{\mathfrak{Q}}}_{\rm
forced}[t,\varepsilon\underline{\zeta}{}_{\pm}](\langle q_{\rm
i}\rangle,\zeta_{\pm})$ (69)
$\displaystyle+\begin{pmatrix}\displaystyle-\frac{1}{2\ell}\llbracket\frac{1}{\underline{h}}R_{1}{\mathfrak{f}}_{\rm
sw}\rrbracket\\\ \displaystyle({R}_{1}{\mathfrak{f}}_{\rm sw})_{+}\\\
\displaystyle({R}_{1}{\mathfrak{f}}_{\rm
sw})_{-}\end{pmatrix}+\begin{pmatrix}\frac{1}{2}\kappa\llbracket\frac{1}{\underline{h}}\rrbracket\\\
-\ell\kappa\\\ -\ell\kappa\end{pmatrix}\ddot{\delta}_{\rm forced},$

where the matrix $\widetilde{\mathcal{M}}_{\rm
forced}[t,\varepsilon\underline{\zeta}{}_{\pm}]$ is block-triangular,

(70) $\widetilde{\mathcal{M}}_{\rm
forced}[t,\varepsilon\underline{\zeta}{}_{\pm}]:=\left(\begin{array}[]{cc|cc}\alpha(\varepsilon\delta_{\rm
forced})+\frac{\kappa}{\ell}\langle\frac{1}{\underline{h}}\rangle&0&0\\\
-{\kappa}&\kappa^{2}&0\\\ \kappa&0&\kappa^{2}\end{array}\right),$

and

$\widetilde{\boldsymbol{\mathfrak{Q}}}_{\rm
forced}[t,\varepsilon\zeta_{\pm}](\langle q_{\rm
i}\rangle,\underline{\zeta}{}_{\pm}):=\begin{pmatrix}\widetilde{\mathfrak{Q}}_{\rm
i}[\varepsilon\delta_{\rm forced},\varepsilon\zeta_{\pm}](\langle q_{\rm
i}\rangle,\dot{\delta}_{\rm forced},\underline{\zeta}{}_{\pm})\\\
\widetilde{\mathfrak{Q}}_{+}[\varepsilon\delta_{\rm
forced},\varepsilon\zeta_{\pm}](\langle q_{\rm i}\rangle,\dot{\delta}_{\rm
forced},\underline{\zeta}{}_{\pm})\\\
\widetilde{\mathfrak{Q}}_{-}[\varepsilon\delta_{\rm
forced},\varepsilon\zeta_{\pm}](\langle q_{\rm i}\rangle,\dot{\delta}_{\rm
forced},\underline{\zeta}{}_{\pm})\end{pmatrix},$

and ${\mathfrak{Q}}_{\rm i}$ and ${\mathfrak{Q}}_{\pm}$ as in the previous
section.

######  Remark A.2.

We have made explicit the dependence of $\widetilde{\mathcal{M}}_{\rm forced}$
and $\widetilde{\boldsymbol{\mathfrak{Q}}}_{\rm forced}$ on the time variable
$t$ because $\delta_{\rm forced}$ is now an explicit function of time, and is
now an non autonomous contribution to the ODE for $\Theta$ (except of course
if the object is fixed, in which case $\delta_{\rm forced}\equiv 0$).

With ${\mathcal{G}}_{\rm I}$ now being three dimensional (but with an extra
dependence on $t$), ${\mathcal{G}}_{\rm
I}:=\big{(}{\mathcal{G}}_{1},{\mathcal{G}}_{2},{\mathcal{G}}_{3}\big{)}^{\rm
T}$, we have therefore

(71) $\frac{d}{dt}\big{(}\langle q_{\rm
i}\rangle,\dot{}\underline{\zeta}_{+},\dot{}\underline{\zeta}_{-}\big{)}^{\rm
T}={\mathcal{G}}_{\rm I}(t,\Theta,(R_{1}{\mathfrak{f}}_{\rm sw})_{\pm})$

and

$\displaystyle{\mathcal{G}}_{\rm I}$
$\displaystyle(t,\Theta,(R_{1}{\mathfrak{f}}_{\rm
sw})_{\pm})=\widetilde{\mathcal{M}}_{\rm
forced}[t,\varepsilon\underline{\zeta}{}_{\pm}]^{-1}$ $\displaystyle\times$
$\displaystyle\Big{[}-\begin{pmatrix}\displaystyle 0\\\
\displaystyle{\underline{\zeta}{}_{+}}\\\
\displaystyle\underline{\zeta}{}_{-}\end{pmatrix}+\varepsilon\widetilde{\boldsymbol{\mathfrak{Q}}}_{\rm
forced}[t,\varepsilon\underline{\zeta}{}_{\pm}](\langle q_{\rm
i}\rangle,\underline{\zeta}{}_{\pm})+\begin{pmatrix}\displaystyle-\frac{1}{2\ell}\llbracket\frac{1}{\underline{h}}R_{1}{\mathfrak{f}}_{\rm
sw}\rrbracket\\\ \displaystyle({R}_{1}{\mathfrak{f}}_{\rm sw})_{+}\\\
\displaystyle({R}_{1}{\mathfrak{f}}_{\rm
sw})_{-}\end{pmatrix}+\begin{pmatrix}\frac{1}{2}\kappa\llbracket\frac{1}{\underline{h}}\rrbracket\\\
-\ell\kappa\\\ -\ell\kappa\end{pmatrix}\ddot{\delta}_{\rm forced}\Big{]},$

while ${\mathcal{G}}_{4}(\Theta)=\dot{\underline{\zeta}{}_{+}}$,
${\mathcal{G}}_{5}(\Theta)=\dot{}\underline{\zeta}{}_{-}$ and

$\widetilde{M}_{\rm
forced}[t,\varepsilon\underline{\zeta}{}_{\pm}]^{-1}=\begin{pmatrix}\displaystyle\frac{1}{\alpha+\frac{\kappa}{\ell}\langle\frac{1}{\underline{h}}\rangle}&0&0\\\
\displaystyle\frac{1}{\kappa}\frac{1}{\alpha+\frac{\kappa}{\ell}\langle\frac{1}{\underline{h}}\rangle}&\frac{1}{\kappa^{2}}&0\\\
\displaystyle-\frac{1}{\kappa}\frac{1}{\alpha+\frac{\kappa}{\ell}\langle\frac{1}{\underline{h}}\rangle}&0&\frac{1}{\kappa^{2}}\end{pmatrix}.$

######  Remark A.3.

The second component of (66) (evolution equation on $\delta$) does not appear
any longer in (71) but remains of course valid. It can be used to answer the
following control problem: what external force should we apply to the object
so that the vertical displacement of its center of gravity coincides with
$\delta_{\rm forced}$? The answer is explicitly given (writing
${\mathfrak{Q}}_{\delta}={\mathfrak{Q}}_{\delta}[\varepsilon\delta_{\rm
forced},\varepsilon\zeta_{\pm}]$, etc.)

$\displaystyle F_{\rm ext}=$
$\displaystyle\delta-\langle\frac{1}{h}R_{1}{\mathfrak{f}}_{\rm
sw}\rangle-\varepsilon{\mathfrak{Q}}_{\delta}(\langle q_{\rm
i}\rangle,\dot{\delta}_{\rm forced},\zeta_{\pm})$
$\displaystyle-\frac{1}{4}\frac{D_{\rm forced}}{\alpha(\varepsilon\delta_{\rm
forced})+\frac{\kappa}{\ell}\langle\frac{1}{h}\rangle}\ddot{\delta}_{\rm
forced}$ $\displaystyle-\frac{1}{2}\frac{\kappa}{\alpha(\varepsilon\delta_{\rm
forced})+\frac{\kappa}{\ell}\langle\frac{1}{h}\rangle}\llbracket\frac{1}{h}\rrbracket\Big{[}\varepsilon{\mathfrak{Q}}_{\rm
i}(\langle q_{\rm i}\rangle,\dot{\delta}_{\rm
forced},\zeta_{\pm})-\frac{1}{2\ell}\llbracket\frac{1}{h}R_{1}{\mathfrak{f}}_{\rm
sw}\rrbracket\Big{]},$

where $D_{\rm forced}$ is deduced from the expression given for $D$ in (68) by
substituting $\delta_{\rm forced}$ to $\delta$.

###  A.4. Simplifications in the symmetric case

When the object is symmetric with respect to the vertical axis $\\{x=0\\}$
(i.e. if $h_{\rm eq}$ is an even function), as assumed throughout this
article, it is possible to consider symmetric flows for which $\zeta$ is an
even function, $q$ is odd, and $\langle q_{\rm i}\rangle\equiv 0$ (such
conditions are propagated by the equations from the initial data). This is for
instance the case for waves generated by a floating object in a fluid
initially at rest. By symmetry, the augmented transmission problem (16)-(18)
reduces to an augmented initial boundary value problem on the half-line
${\mathcal{E}}^{+}=(\ell,\infty)$,

(72) $\begin{cases}\partial_{t}\zeta+\partial_{x}q=0,\\\
(1-\kappa^{2}\partial_{x}^{2})\partial_{t}q+\varepsilon\partial_{x}\big{(}\frac{1}{h}q^{2}\big{)}+h\partial_{x}\zeta=0\end{cases}\quad\mbox{
for }\quad t>0,\quad x\in(\ell,\infty)$

with boundary condition

(73) $q_{|_{x=\ell}}=-\ell\dot{\delta},$

where $\delta$ is a function of time determined by the first order ODE

(74) $\frac{d}{dt}\Theta={\mathcal{G}}\big{(}\Theta,(R_{1}{\mathfrak{f}}_{\rm
sw})_{+},F_{\rm ext}\big{)},$

with
$\Theta:=\big{(}\dot{\delta},\dot{\underline{\zeta}{}_{+}},\delta,{\underline{\zeta}{}_{+}}\big{)}^{\rm
T}$ and where
${\mathcal{G}}=({\mathcal{G}}_{1},{\mathcal{G}}_{2},{\mathcal{G}}_{3},{\mathcal{G}}_{4})^{\rm
T}$ is given by

$\displaystyle\begin{pmatrix}{\mathcal{G}}_{1}\\\
{\mathcal{G}}_{2}\end{pmatrix}={\mathcal{M}}_{\rm
sym}[\varepsilon\delta,\varepsilon\zeta_{+}]^{-1}\Big{[}$
$\displaystyle-\begin{pmatrix}\delta\\\
\zeta_{+}\end{pmatrix}+\varepsilon\begin{pmatrix}{\mathfrak{Q}}_{\delta}^{\rm
sym}[\varepsilon\delta,\varepsilon\zeta_{+}](\dot{\delta},\zeta_{+})\\\
{\mathfrak{Q}}_{+}^{\rm
sym}[\varepsilon\delta,\varepsilon\zeta_{+}](\dot{\delta},\zeta_{+})\end{pmatrix}$
$\displaystyle+\begin{pmatrix}\frac{1}{h_{+}}(R_{1}{\mathfrak{f}}_{\rm
sw})_{+}+F_{\rm ext}\\\ (R_{1}{\mathfrak{f}}_{\rm
sw})_{+}\end{pmatrix}\Big{]}$

and ${\mathcal{G}}_{3}=\dot{\delta}$, ${\mathcal{G}}_{4}=\dot{\zeta}_{+}$, and
with

${\mathcal{M}}_{\rm
sym}[\varepsilon\delta,\varepsilon\zeta_{+}]=\begin{pmatrix}\tau_{\kappa}(\varepsilon\delta)^{2}+\kappa\ell{\frac{1}{h_{+}}}&0\\\
\ell\kappa&\kappa^{2}\end{pmatrix},$

while ${\mathfrak{Q}}_{\delta}^{\rm sym}$ and ${\mathfrak{Q}}_{+}^{\rm sym}$
are obtained by replacing $\langle q_{\rm i}\rangle=0$ and
$\zeta_{-}=\zeta_{+}$ in the formula derived above for
${\mathfrak{Q}}_{\delta}$ and ${\mathfrak{Q}}_{\delta}$.

Two particular physical situations of particular interest fit into the
symmetric framework and are investigated in this paper:

  * •

The return to equilibrium. An object is released from an out of equilibrium
position in a fluid initially at rest. This situation is described by
(72)-(74) with $F_{\rm ext}=0$ and initial conditions
$(\zeta,q)_{|_{t=0}}=(0,0)$ and
$(\dot{\delta},\dot{\zeta}_{+},\delta,\zeta_{+})_{|_{t=0}}=(0,0,\delta^{\rm
in},0)$.

  * •

Wave generation. Waves are generated in a fluid initially at rest by moving
the object up and down with a prescribed motion $\delta_{\rm forced}$. The
problem then reduces to an initial boundary value problem with boundary
condition on $q$, namely, $q_{|_{x=\ell}}=g$, with $g=-\ell\dot{\delta}_{\rm
forced}$. This boundary data is explicitly given and does not require the
resolution of a first order ODE as the other problems considered here.

## Acknowledgment

D. L. was partially supported by the grant ANR-18-CE40-0027-01 Singflows.

## References

  * [1] J. Abate, W. Ward, A Unified Framework for Numerically Inverting Laplace Transforms, INFORMS Journal of Computing 18(4)  (2006) 408–421
  * [2] T. Alazard, N. Burq, C. Zuily, Low regularity Cauchy theory for the water-waves problem: canals and swimming pools, Journées équations aux dérivées partielles (2011) 1–20.
  * [3] A. Babarit, Ocean wave energy conversion: resource, technologies and performance, Elsevier (2017).
  * [4] G. Beck, D. Lannes, Freely floating objects on a fluid governed by the Boussinesq equations, Ann. Inst. H. Poincaré Anal. Non Linéaire 39 (2022) 575–646.
  * [5] R. Bianchini, C. Perrin, Soft congestion approximation to the one-dimensional constrained Euler equations, Nonlinearity 34(10) (2021) 6901.
  * [6] U. Bosi, A. P. Engsig-Karup, C. Eskilsson, M. Ricchiuto, A spectral/hp element depth-integrated model for nonlinear wave–body interaction, Comput. Methods Appl. Mech. Eng. 348 (2019) 222–249.
  * [7] E. Bocchi, Floating structures in shallow water: local well-posedness in the axisymmetric case, SIAM J. Math. Anal. 52 (2020) 306–339.
  * [8] E. Bocchi, On the return to equilibrium problem for axisymmetric floating structures in shallow water, Nonlinearity 33 (2020) 3594.
  * [9] E. Bocchi, J. He, G. Vergara-Hermosilla, Modelling and simulation of a wave energy converter, ESAIM Proc. 70 (2021) 68–83
  * [10] D. Bresch, D. Lannes, G. Métivier, Waves interacting with a partially immersed obstacle in the Boussinesq regime, Anal. PDE 14 (2021) 1085–1124.
  * [11] W. Cummins, The Impulse Response Function and Ship Motions, Report (David W. Taylor Model Basin), Navy Department, David Taylor Model Basin, 1962.
  * [12] A.-L. Dalibard, C. Perrin, Partially congested propagation fronts in one-dimensional Navier-Stokes equations, J. Elliptic Parabol Equ. 7 (2021) 491–507.
  * [13] A. P. Engsig‐Karup, W. L. Laskowski, An efficient p‐multigrid spectral element model for fully nonlinear water waves and fixed bodies, Int. J. Numer. Methods Fluids 93 (2021) 2823–2841.
  * [14] A. Filippini, S. de Brye, V. Perrier, F. Marche, M. Ricchiuto, D. Lannes, P. Bonneton, UHAINA : A parallel high performance unstructured near-shore wave model, Journées Nationales Génie Côtier – Génie Civil, May 2018, La Rochelle, France. 47–56,
  * [15] E. Godlewski, M. Parisot, J. Sainte-Marie, F. Wahl, Congested shallow water model: roof modeling in free surface flow, ESAIM: Math. Model. Numer. Anal. 52 (2018) 1679–1707.
  * [16] E. Godlewski, M. Parisot, J. Sainte-Marie, F. Wahl, Congested shallow water model: on floating body, SMAI J. Comput. Math. 6 (2020) 227–251.
  * [17] O. I. Gusev, G. S. Khakimzyanov, L. B. Chubarov, Numerical investigation of the wave force on a partially immersed rectangular structure: Long waves over a flat bottom, Ocean Eng. 221 (2021) 108540.
  * [18] A. Haidar, F. Marche, F. Vilar, A robust DG-ALE formulation for nonlinear shallow water interactions with a partially immersed object, preprint hal-03764650 (2022)
  * [19] T. Iguchi, D. Lannes, Hyperbolic free boundary problems and applications to wave-structure interactions, Indiana Univ. Math. J. 70 (2021) 353–464.
  * [20] F. John, On the motion of floating bodies. I, Commun. Pure Appl. Math. 2 (1949) 13–57.
  * [21] T. Karambas, E. Loukogeorgaki, A Boussinesq-type model for nonlinear wave-heaving cylinder interaction, Energies 15 (2022) 469.
  * [22] D. Lannes, Modeling shallow water waves, Nonlinearity 33 (2020) R1
  * [23] D. Lannes, Initial boundary value problems for hyperbolic systems, and dispersive perturbations, Lecture notes of the Bressanone Winter School, to appear
  * [24] D. Lannes, On the dynamics of floating structures, Ann. PDE 3 (2017).
  * [25] D. Lannes, G. Métivier, The shoreline problem for the one-dimensional shallow water and Green-Naghdi equations, J. Ec. Polytech. Math. 5 (2018) 455–518.
  * [26] D. Lannes, L. Weynans, Generating boundary conditions for a Boussinesq system, Nonlinearity 33 (2020) 6868.
  * [27] R.J. LeVeque, Numerical Methods for Conservation Laws, Birkhauser-Verlag (1990).
  * [28] D Maity, J San Martín, T Takahashi, M Tucsnak, Analysis of a simplified model of rigid structure floating in a viscous fluid, J. Nonlinear Sci. 29 (2019) 1975–2020.
  * [29] M. Ming, C. Wang, Water‐Waves Problem with Surface Tension in a Corner Domain II: The Local Well‐Posedness, Commun. Pure Appl. Math. 74 (2021) 225–285.
  * [30] S. C. Mohapatra, H. Islam, C. Guedes Soares, Boussinesq model and CFD simulations of non-linear wave diffraction by a floating vertical cylinder, J. Mar. Sci. Eng. 8 (2020) 575
  * [31] M. Parisot, Congested shallow water model: trapped air pockets modeling, preprint hal-03748169 (2022).
  * [32] M. Penalba, G. Giorgi, J. Ringwood, Mathematical modelling of wave energy converters: A review of nonlinear approaches, Renew. Sustain. Ener. Rev. 78 (2017) 1188–1207.
  * [33] C. Perrin, K. Saleh, Numerical staggered schemes for the free-congested Navier-Stokes equations, SIAM J. Numer. Anal. 60 (2022) 1824–1852.
  * [34] T. de Poyferré, A priori estimates for water waves with emerging bottom, Arch. Ration. Mech. Anal. 232 (2019) 763–812.
  * [35] F. Shi, J. T. Kirby, J. C. Harris, J. D. Geiman, S. T. Grilli, A high-order adaptive time-stepping tvd solver for Boussinesq modeling of breaking waves and coastal inundation, Ocean Model. 43-44 (2012) 36–51.
  * [36] P. Su, M. Tucsnak, Shallow water waves generated by a floating object: a control theoretical perspective, Math. Control Relat. Fields 16 (2021).
  * [37] G. Wei G, J. Kirby, A. Sinha, Generation of waves in Boussinesq models using a source function method, Coast. Eng. 36 (1999) 271–99.