# A NUMERICAL METHOD FOR WAVE-STRUCTURE INTERACTIONS IN THE BOUSSINESQ REGIME

Geoffrey Beck, David Lannes, Lisl Weynans

<figure><img src="image_1.png" /></figure>

## To cite this version:

Geoffrey Beck, David Lannes, Lisl Weynans. A NUMERICAL METHOD FOR WAVE-STRUCTURE INTERACTIONS IN THE BOUSSINESQ REGIME. 2023. hal-04151128

HAL Id: hal-04151128

https://hal.science/hal-04151128

Preprint submitted on 4 Jul 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
A NUMERICAL METHOD FOR WAVE-STRUCTURE
INTERACTIONS IN THE BOUSSINESQ REGIME

Geoffrey Beck, David Lannes, Lisl Weynans

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.

**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><img src="image_3.png" /><figcaption>FIGURE 1. The floating object</figcaption></figure>

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 \text{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 \text{and} \quad \langle q \rangle = \langle q_i \rangle $$

where the function $\delta$ and $\langle q_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_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 know $\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 controllability 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 non-linear 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),

$$ \left\{ \begin{array}{ll} \partial_t \zeta + \partial_x q = 0, & \text{on } (-\infty, -\ell) \cup (\ell, +\infty), \\ (1 - \kappa^2 \partial_x^2) \partial_t q + \varepsilon \partial_x (\frac{1}{h} q^2) + h \partial_x \zeta = 0, & \text{on } (-\infty, -\ell) \cup (\ell, +\infty), \end{array} \right. $$

and with transmission conditions

$$ \|q\| = -2\ell\dot{\delta} \quad \text{and} \quad \langle q \rangle = \langle q_i \rangle $$

where the function $\delta$ and $\langle q_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_i \rangle, \dot{\delta}, \delta)^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 $\zeta_{\pm}$ defined as the solutions of this ODE. This allows us to replace $\zeta|_{x=\pm\ell}$ by $\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 f_{sw})|_{x=\pm\ell}) $$

with $\Theta = (\langle q_i \rangle, \dot{\delta}, \dot{\zeta}_+, \dot{\zeta}_-, \delta, \zeta_+, \zeta_-)^T$ and $\mathcal{G}$ is a smooth mapping. This ODE is forced by the terms $(R_1 f_{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 f_{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)) = S_\pm(\Theta, (R_1 f_{sw})|_{x=\pm\ell}) \mathfrak{b}(x \mp \ell) \quad \text{in} \quad \pm(\ell, \infty), $$

where $\mathfrak{F}_\kappa$ denotes the nonlocal flux, while $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 f_{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(\bar{\mathcal{E}})$, we denote

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

- If $f \in C^1([0, T])$, we sometimes use the notation $\dot{f} = \frac{d}{dt}f$.

- We denote by $f_{\text{sw}}$ the momentum flux associated with the nonlinear shallow water equations,

$$(1) \qquad f_{\text{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{\text{typical amplitude of the waves}}{\text{typical depth}} \quad \text{and} \quad \mu = \left( \frac{\text{typical depth}}{\text{typical horizontal scale}} \right)^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 \text{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 = \left( \frac{\mu}{3} \right)^{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_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_{\text{eq}}$ the water depth at rest. These quantities are related through

$$ \zeta_w(t, x) = \delta(t) + \frac{1}{\varepsilon}(h_{\text{eq}}(x) - 1). $$

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_{\text{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) \qquad \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 = -\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_{\text{atm}}$,

$$ \underline{P}(t, x) = P_{\text{atm}} \quad \text{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_w(t, x) \quad \text{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_i \rangle(t), $$

where $\langle q_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 \left( \frac{1}{\varepsilon} h_w \partial_x \underline{P} \right) = -\ddot{\delta} + \partial_x \left[ h_w \partial_x \zeta_w + \varepsilon \partial_x \left( \frac{1}{h_w} \left( -x\dot{\delta} + \langle q_i \rangle \right)^2 \right) \right], $$

for $x \in (-\ell, \ell)$ and with $h_w = h_{\text{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_{\text{eq}}$, $\delta$, $\langle q_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_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) \qquad q(t, -\ell - 0) = \ell \dot{\delta} + \langle q_i \rangle(t) \quad \text{and} \quad q(t, \ell + 0) = -\ell \dot{\delta} + \langle q_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_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_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_{\text{ext}}$ represents the vertical force exerted on the solid to maintain it fixed or with the desired motion). The outcome of this analysis, as shown<fnref n="1" /> 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} (f(\ell) - f(-\ell)) \quad \text{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) \qquad \begin{cases} \partial_t \zeta + \partial_x q = 0, \\ (1 - \kappa^2 \partial_x^2) \partial_t q + \partial_x f_{\text{sw}} = 0 \end{cases} \quad \text{for } t > 0, \quad x \in \mathcal{E} $$

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

$$ (5) \qquad \langle q \rangle = \langle q_i \rangle \quad \text{and} \quad \llbracket q \rrbracket = -2\ell\dot{\delta}, $$

where $\langle q_i \rangle$ and $\delta$ are functions of time solving

$$ (6) \qquad \alpha(\varepsilon\delta) \frac{d}{dt} \langle q_i \rangle + \varepsilon\alpha'(\varepsilon\delta)\dot{\delta}\langle q_i \rangle = -\frac{1}{2\ell} \llbracket \zeta + \mathfrak{G} \rrbracket, $$

$$ (7) \qquad \tau_{\kappa}(\varepsilon\delta)^2 \ddot{\delta} + \delta - \varepsilon\beta(\varepsilon\delta)\dot{\delta}^2 - \varepsilon\frac{1}{2}\alpha'(\varepsilon\delta)\langle q_i \rangle^2 = \langle \zeta + \mathfrak{G} \rangle + F_{\text{ext}}, $$

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

$$ (8) \qquad \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

<footnote n="1">¹The 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$.</footnote>
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_i \rangle$ and $\delta$. This coupling acts in both ways: on the one hand, it is necessary to know $\langle q_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_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_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} \left( \frac{\mp \ell \dot{\delta} + \langle q_i \rangle}{1 + \varepsilon \zeta_{\pm}} \right)^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 \text{with} \quad \begin{cases} (1 - \kappa^2 \partial_x^2) u = F & \text{on } \mathcal{E}, \\ u|_{x=\pm \ell} = 0, \end{cases} $$

and

$$ (9) \qquad R_1 F = v \quad \text{with} \quad \begin{cases} (1 - \kappa^2 \partial_x^2) v = F & \text{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) \qquad \begin{cases} \partial_t \zeta + \partial_x q = 0, \\ (1 - \kappa^2 \partial_x^2) \partial_t q + \partial_x f_{sw} = 0, \end{cases} \quad t > 0, \quad x \in \mathcal{E} $$

with $f_{sw}$ as in (1) and with transmission conditions

$$ (11) \qquad \langle q \rangle(t) = f(t) \quad \text{and} \quad \|q\|(t) = 2g(t), \quad t > 0 $$
then $\zeta_{\pm} = \zeta_{|x=\pm\ell}$ solve the ODEs

$$ (12) \qquad \begin{cases} \partial_t^2 \zeta_+ + \frac{1}{\kappa^2} \zeta_+ + \frac{\varepsilon}{\kappa^2} \left( \frac{1}{2} \zeta_+^2 + \frac{(f+g)^2}{1+\varepsilon\zeta_+} \right) = \frac{1}{\kappa^2} (R_1 f_{sw})_+ + \frac{1}{\kappa} (\dot{f} + \dot{g}), \\ \partial_t^2 \zeta_- + \frac{1}{\kappa^2} \zeta_- + \frac{\varepsilon}{\kappa^2} \left( \frac{1}{2} \zeta_-^2 + \frac{(f-g)^2}{1+\varepsilon\zeta_-} \right) = \frac{1}{\kappa^2} (R_1 f_{sw})_- - \frac{1}{\kappa} (\dot{f} - \dot{g}), \end{cases} $$

where we used the notation $(R_1 f_{sw})_{\pm} = (R_1 f_{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 f_{sw} = (\dot{f} \pm \dot{g}) \exp\left(-\frac{|x \mp \ell|}{\kappa}\right) \quad \text{on } \mathcal{E}^{\pm}. $$

Remarkng further that $R_0 \partial_x = \partial_x R_1$, the problem is therefore reduced to

$$ (13) \qquad \begin{cases} \partial_t \zeta + \partial_x q = 0, \\ \partial_t q + \partial_x R_1 f_{sw} = (\dot{f} \pm \dot{g}) \exp\left(-\frac{|x \mp \ell|}{\kappa}\right). \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 f_{sw} = \mp \frac{1}{\kappa} (\dot{f} \pm \dot{g}) \exp \left( -\frac{1}{\kappa} (|x \mp \ell|) \right). $$

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} f_{sw} + \frac{1}{\kappa^2} R_1 f_{sw} = \mp \frac{1}{\kappa} (\dot{f} \pm \dot{g}) \exp \left( -\frac{1}{\kappa} (|x \mp \ell|) \right). $$

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

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

$$ (14) \quad \mathcal{M}[\varepsilon\delta, \varepsilon\zeta_{\pm}] \frac{d}{dt} \begin{pmatrix} \langle q_i \rangle \\ \dot{\delta} \\ \dot{\zeta}_+ \\ \dot{\zeta}_- \end{pmatrix} + \begin{pmatrix} \frac{1}{2\ell} \llbracket \zeta \rrbracket \\ \delta - \langle \zeta \rangle \\ \zeta_+ \\ \zeta_- \end{pmatrix} = \varepsilon \mathfrak{Q} [\varepsilon\delta, \varepsilon\zeta_{\pm}] (\langle q_i \rangle, \dot{\delta}, \zeta_{\pm}) + \begin{pmatrix} 0 \\ F_{\text{ext}} \\ (R_1 f_{\text{sw}})_+ \\ (R_1 f_{\text{sw}})_- \end{pmatrix}, $$

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

$$ (15) \quad \mathcal{M}[\varepsilon\delta, \varepsilon\zeta_{\pm}] := \begin{pmatrix} \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_-} \\ -\kappa & \ell\kappa & \kappa^2 & 0 \\ \kappa & \ell\kappa & 0 & \kappa^2 \end{pmatrix}, $$

while $\mathfrak{Q}[\varepsilon\delta, \varepsilon\zeta_{\pm}] (\langle q_i \rangle, \dot{\delta}, \zeta_{\pm})$ is a four-dimensional vector whose entries are quadratic forms in $(\langle q_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_i \rangle$, $\dot{\delta}$, $\zeta_+$, $\zeta_-$, $\delta$, $\zeta_+$ and $\zeta_-$ with forcing terms $(R_1 f_{sw})_{\pm}$ and $F_{\text{ext}}$ (see Appendix A.2). The augmented formulation is obtained
by replacing $\zeta_+$ and $\zeta_-$ by two additional unknowns $\zeta_+$ and $\zeta_-$ in this first order ODE. It reads therefore

$$ (16) \qquad \begin{cases} \partial_t \zeta + \partial_x q = 0, \\ (1 - \kappa^2 \partial_x^2) \partial_t q + \partial_x f_{sw} = 0 \end{cases} \quad \text{for } t > 0, x \in \mathcal{E}, $$

where $f_{sw}$ is as in (1), and with transmission conditions across the floating object given by

$$ (17) \qquad \langle q \rangle = \langle q_i \rangle \quad \text{and} \quad \llbracket q \rrbracket = -2\ell\delta, $$

where $\langle q_i \rangle$ and $\delta$ are functions of time determined by the first order ODE

$$ (18) \qquad \frac{d}{dt} \Theta = \mathcal{G}(\Theta, (R_1 f_{sw})_+, (R_1 f_{sw})_-, F_{\text{ext}}), $$

with $\Theta := (\langle q_i \rangle, \dot{\delta}, \dot{\zeta}_+, \dot{\zeta}_-, \delta, \zeta_+, \zeta_-)^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 $\zeta_+$ and $\dot{\zeta}_+$ are chosen appropriately, then $\zeta_+ = \zeta_+$ and $\dot{\zeta}_+ = \dot{\zeta}_+$ 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_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_i \rangle$), which is forced by $F_{\text{ext}}$, $(R_1 f_{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_i \rangle, \zeta_+$ and $\dot{\zeta}_+$), but it is forced only by $F_{\text{ext}}$ and $(R_1 f_{sw})_{|_{x=\pm\ell}}$. These two quantities do not raise any difficulty since $F_{\text{ext}}$ is a given external force and $(R_1 f_{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 $\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 $\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) \qquad S_{\pm}(\Theta, (R_1 f_{sw})_{\pm}, F_{\text{ext}}) := G_1(\Theta, (R_1 f_{sw})_{\pm}, F_{\text{ext}}) \mp \ell G_2(\Theta, (R_1 f_{sw})_{\pm}, F_{\text{ext}}), $$

where $G_1$ and $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 = (\langle q_i \rangle, \dot{\delta}, \dot{\zeta}_+, \dot{\zeta}_-, \delta, \zeta_+, \zeta_-)^T$ be a regular solution to the transmission problem (16)-(18) with initial data $U = (\zeta^{\text{in}}, q^{\text{in}})$ and $\Theta^{\text{in}}$. Then $(\zeta, q)$ and $\Theta$ also solve the initial value problem

$$ (20) \qquad \left\{ \begin{array}{l} \partial_t \zeta + \partial_x q = 0, \\ \partial_t q + \partial_x R_1 f_{sw} = S_{\pm}(\Theta, (R_1 f_{sw})_{\pm}, F_{\text{ext}}) \exp(-\frac{1}{\kappa}|x \mp \ell|) \end{array} \right. \quad \text{on } \mathcal{E}^{\pm}, $$

and

$$ (21) \qquad \frac{d}{dt} \Theta = \mathcal{G}(\Theta, (R_1 f_{sw})_+, (R_1 f_{sw})_-, F_{\text{ext}}). $$

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

$$ (22) \qquad \langle q^{\text{in}} \rangle = \Theta_1^{\text{in}}, \qquad \llbracket q^{\text{in}} \rrbracket = -2\ell\Theta_2^{\text{in}}. $$

If moreover the initial data also satisfy

$$ (23) \qquad \zeta_{|x=\ell}^{\text{in}} = \Theta_6^{\text{in}}, \quad \zeta_{|x=-\ell}^{\text{in}} = \Theta_7^{\text{in}}, \quad -(\partial_x q^{\text{in}})|_{x=\ell} = \Theta_3^{\text{in}}, \quad -(\partial_x q^{\text{in}})|_{x=-\ell} = \Theta_4^{\text{in}}, $$

then for all times, one has $\zeta_{|x=\pm\ell} = \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_{\text{ext}} \equiv 0$.

* When the data and the object are symmetric with respect to the vertical axis $x=0$, then $\langle q_i \rangle = 0$.

* If the object is in forced motion, i.e. if $\delta \equiv \delta_{\text{forced}}$ for some given function $\delta_{\text{forced}}$, then the ODE (21) can be reduced to a 5-dimensional ODE on $(\langle q_i \rangle, \dot{\zeta}_+, \dot{\zeta}_+, \zeta_+, \zeta_-)^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 f_{sw} = \frac{d}{dt} (q_{|x=\pm\ell}) \exp(-\frac{1}{\kappa}|x \mp \ell|) \quad \text{on } \mathcal{E}^{\pm}. $$

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

$$ \frac{d}{dt} \langle q_i \rangle = \mathcal{G}_1(\Theta, (R_1 f_{sw})_{\pm}, F_{\text{ext}}) \quad \text{and} \quad \frac{d}{dt} \dot{\delta} = \mathcal{G}_2(\Theta, (R_1 f_{sw})_{\pm}, F_{\text{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 f_{sw})_{\pm}, F_{\text{ext}}) \quad \text{and} \quad \frac{d}{dt} \llbracket q \rrbracket = -2\ell \mathcal{G}_2(\Theta, (R_1 f_{sw})_{\pm}, F_{\text{ext}}), $$
or equivalently (from the definition of $\mathcal{G}$),

$$ \frac{d}{dt} \langle q \rangle = \frac{d}{dt} \langle q_i \rangle \quad \text{and} \quad \frac{d}{dt} \llbracket q \rrbracket = \frac{d}{dt} (-2\ell\dot{\delta}). $$

This shows that the time derivative of the transmission conditions (17) are satisfied; the compatibility conditions (22) show moreover that the transmission 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 $\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. $\square$

### 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 an exponentially localized source term,

$$ (24) \quad \partial_t U + \partial_x (\mathfrak{F}_\kappa(U)) = S_\pm(\Theta, (R_1 f_{sw})_-, (R_1 f_{sw})_+, F_{\text{ext}})b(x \mp \ell) \quad \text{in } \mathcal{E}_\pm $$

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

$$ (25) \quad \mathfrak{F}_\kappa(U) = (q, R_1 f_{sw})^T, $$

while the source terms $S_\pm$ are as in (19) and $b$ is the shape of source term

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

The quantity $\Theta$ is defined as $\Theta = (\langle q_i \rangle, \dot{\delta}, \dot{\zeta}_+, \dot{\zeta}_-, \delta, \zeta_+, \zeta_-)^T$ and solves a system of 7 first order ODEs forced by $(R_1 f_{sw})_+, (R_1 f_{sw})_-$ and $F_{\text{ext}}$,

$$ (27) \quad \frac{d}{dt} \Theta = \mathcal{G}(\Theta, (R_1 f_{sw})_+, (R_1 f_{sw})_-, F_{\text{ext}}). $$

*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).

<figure><img src="image_4.png" /><figcaption>FIGURE 2. Space discretization</figcaption></figure>
3.1.1. Discretization. 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}^- = \left( \bigcup_{i=-\infty}^{-1} C_i \right) \cup C_{0^-} \quad \text{and} \quad \mathcal{E}^+ = C_{0^+} \cup \left( \bigcup_{i=1}^{\infty} C_i \right) $$

with

$$ C_i = (x_{i-1/2}, x_{i+1/2}) \text{ if } i \neq 0 \quad \text{and} \quad C_{0^-} = (x_{-1/2}, -\ell), \quad C_{0^+} = (\ell, x_{1/2}) $$

and where

$$ x_{i+1/2} = -\ell + (i + 1/2)\Delta x \quad \text{if } i < 0 \quad \text{and} \quad x_{i-1/2} = \ell + (i - 1/2)\Delta x \quad \text{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), \quad \Theta^n = \Theta(n\Delta t) \quad \text{and} \quad F^n_{\text{ext}} = F_{\text{ext}}(n\Delta t) $$

the values of $U = (\zeta, q)^T$, of the $\mathbb{R}^7$-valued vector $\Theta$ involved in the ODE (27), and of the external force $F_{\text{ext}}$ at each time step. We further denote by $U^n$ ($i \in (\mathbb{Z}\setminus\{0\}) \cup \{0^-, 0^+\}$) the approximation of $U^n$ in the middle of the cell $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 f_{\text{sw}}$ and $(R_1 f_{\text{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 f_{\text{sw}})_\pm$ stands for the trace of $R_1 f_{\text{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| \ge 1}$, we denote by $R_1 F$ the vector $R_1 F = V$ where $V = (v_i)_{|i| \ge 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, \quad |i| \ge 2 $$

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

$$ v_{-1} - \kappa^2 \frac{2v_{-2} - v_{-1}}{3\Delta x^2} = f_{-1} \quad \text{and} \quad v_1 - \kappa^2 \frac{2v_2 - v_1}{3\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) \quad (R_1 F)_- = \frac{4}{3}v_{-1} - \frac{1}{3}v_{-2} \quad \text{and} \quad (R_1 F)_+ = \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) \qquad \frac{U_i^{n+1} - U_i^n}{\Delta t} + \frac{1}{\Delta x} (\tilde{\mathfrak{F}}_{\kappa,i+1/2}^n - \tilde{\mathfrak{F}}_{\kappa,i-1/2}^n) = S_{\pm}^n \mathfrak{b}_i, \quad \pm i \ge 1, \ n \ge 0, $$

with

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

the discrete flux correspond to the Lax-Friedrichs scheme,

$$ (31) \qquad \begin{cases} \tilde{\mathfrak{F}}_{\kappa,i+1/2}^n = \frac{1}{2}(\tilde{\mathfrak{F}}_{\kappa,i+1}^n + \tilde{\mathfrak{F}}_{\kappa,i}^n) - \frac{\Delta x}{2\Delta t}(U_{i+1}^n - U_i^n) & \text{if } i \le -2, \\ \tilde{\mathfrak{F}}_{\kappa,i-1/2}^n = \frac{1}{2}(\tilde{\mathfrak{F}}_{\kappa,i}^n + \tilde{\mathfrak{F}}_{\kappa,i-1}^n) - \frac{\Delta x}{2\Delta t}(U_i^n - U_{i-1}^n) & \text{if } i \ge 2, \end{cases} $$

with the notations

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

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

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

with

$$ (33) \qquad \tilde{\mathfrak{F}}_{\kappa,0\pm}^n = (q_{0\pm}^n, (R_1 \mathfrak{f}_{\text{sw}}^n)_{\pm})^\text{T}; $$

the component $(R_1 \mathfrak{f}_{\text{sw}}^n)_\pm$ is computed according to (28), but we still need to define $q_{0\pm}^n$. By definition $q_{0\pm}^n$ 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_i \rangle \mp \ell\delta$. Recalling also that $\langle q_i \rangle$ and $\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_{0\pm}^n$,

$$ (34) \qquad q_{0\pm}^n = \Theta_1^n \mp \ell \Theta_2^n. $$

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

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

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) \qquad U_i^0 = (\zeta_i^{\text{in}}, q_i^{\text{in}})_i^\text{T} \quad (i \in (\mathbb{Z}\setminus\{0\}) \cup \{0^-, 0^+\}), $$

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

$$ (37) \qquad \Theta^0 = (\langle q_i \rangle^{\text{in}}, \delta^{(1)}, \zeta_+^{(1)}, \zeta_-^{(1)}, \delta^{(0)}, \zeta_+^{(0)}, \zeta_-^{(0)})^\text{T} $$

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

$$ (38) \qquad \langle q_i^{\text{in}} \rangle = \langle q_i \rangle^{\text{in}}, \quad [\![q_i^{\text{in}}] \!]= -2\ell\delta^{(1)}, \quad \underline{\zeta}_{\pm}^{(0)} = \zeta_{\pm}^{\text{in}}, \quad \underline{\zeta}_{\pm}^{(1)} = -(\partial_x q_i^{\text{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) \qquad \frac{U_i^{n,*} - U_i^n}{\Delta t} + \frac{1}{\Delta x} (\mathfrak{F}_{\kappa,i}^n - \mathfrak{F}_{\kappa,i-1}^n) = S_+^n b_i, \quad i > 1, n \ge 0, $$

with $U_i^{n,*} = (\zeta_i^{n,*}, q_i^{n,*})^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) \qquad \frac{U_i^{n,*} - U_i^n}{\Delta t} + \frac{1}{\Delta x} (\mathfrak{F}_{\kappa,i+1}^n - \mathfrak{F}_{\kappa,i}^n) = S_-^n b_i, \quad i < -1, n \ge 0, $$

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

$$ (41) \qquad \frac{U_1^{n,*} - U_1^n}{\Delta t} + \frac{1}{\Delta x} (\mathfrak{F}_{\kappa,1}^n - \mathfrak{F}_{\kappa,0+}^n) = S_+^n b_1, $$

$$ (42) \qquad \frac{U_{-1}^{n,*} - U_{-1}^n}{\Delta t} + \frac{1}{\Delta x} (\mathfrak{F}_{\kappa,0-}^n - \mathfrak{F}_{\kappa,-1}^n) = S_-^n b_{-1}, $$

for $n \ge 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) \qquad \frac{\Theta^{n,*} - \Theta^n}{\Delta t} = \mathcal{G}(\Theta^n, (R_1 f_{sw}^n)_+, (R_1 f_{sw}^n)_-, F_{\text{ext}}^n). $$

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

$$ (44) \qquad f_{sw}^{n,*} = f_{sw}(U^{n,*}) \quad \text{and} \quad q_{0\pm}^{n,*} = \Theta_1^{n,*} \mp \ell \Theta_2^{n,*} $$

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

$$ (45) \qquad \mathfrak{F}_{\kappa,i}^{n,*} = (q_i^{n,*}, (R_1 f_{sw}^{n,*})_i)^T \quad |i| \ge 1, n \ge 0, $$

$$ (46) \qquad S_{\pm}^{n,*} = S_{\pm}(\Theta^{n,*}, (R_1 f_{sw}^{n,*})_+, (R_1 f_{sw}^{n,*})_-, F_{\text{ext}}^n) \quad n \ge 0. $$

The correction step for the MacCormack scheme then reads

$$ (47) \qquad \frac{U_i^{n+1} - U_i^n}{\Delta t} + \frac{\mathfrak{F}_{\kappa,i}^n - \mathfrak{F}_{\kappa,i-1}^n + \mathfrak{F}_{\kappa,i+1}^{n,*} - \mathfrak{F}_{\kappa,i}^{n,*}}{2\Delta x} = \frac{S_+^n + S_+^{n,*}}{2} b_i \quad i \ge 1, $$

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

$$ (48) \qquad \frac{U_i^{n+1} - U_i^n}{\Delta t} + \frac{\mathfrak{F}_{\kappa,i+1}^n - \mathfrak{F}_{\kappa,i}^n + \mathfrak{F}_{\kappa,i}^{n,*} - \mathfrak{F}_{\kappa,i-1}^{n,*}}{2\Delta x} = \frac{S_-^n + S_-^{n,*}}{2} b_i \quad i \le -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 \ge 0$,

$$ (49) \qquad \frac{\Theta^{n+1} - \Theta^n}{\Delta t} = \frac{\mathcal{G}(\Theta^n, (R_1 \mathfrak{f}_{\text{sw}}^n)_{\pm}, F_{\text{ext}}^n) + \mathcal{G}(\Theta^{n,*}, (R_1 \mathfrak{f}_{\text{sw}}^{n,*})_{\pm}, F_{\text{ext}}^{n+1})}{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_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_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_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) \qquad \left\{ \begin{array}{ll} \partial_t \zeta + \partial_x q = 0, & t > 0, \quad x > \ell \\ (1 - \kappa^2 \partial_x^2) \partial_t q + \partial_x f_{\text{sw}} = 0, & t > 0, \quad x < \ell \end{array} \right. $$

with $f_{\text{sw}}$ as in (1) and with boundary condition

$$ (51) \qquad q|_{x=\ell}(t) = g(t), \quad t > 0 $$
and initial condition

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

and where g, ζ<sup>in</sup> and q<sup>in</sup> are some given functions satisfying the compatibility
condition

$$
(53) \qquad q^{\text{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 ζ rather than q was considered and where a
first order scheme was proposed.

Remark 4.1. For the wave generation problem, one has $(\zeta^{\text{in}}, q^{\text{in}}) = (0,0)$ and $g(t) = -\ell \dot{\delta}_{\text{forced}}$, where $\delta_{\text{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) \qquad \begin{cases} \partial_t \zeta + \partial_x q = 0, \\ \partial_t q + \partial_x R_1 f_{sw} = S_+(t) \exp\left(-\frac{x-\ell}{\kappa}\right) \end{cases} \quad \text{with} \quad 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:

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

and for the second-order scheme:

$$
S^n = \frac{g^{n+1} - g^{n-1}}{2\Delta t} \quad \text{and} \quad 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) = (\zeta_c(x - x_0 - ct), c\zeta_c(x - x_0 - ct)),
$$

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) \qquad c^2 \kappa^2 \zeta_c'' - 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_{\max}^2 + \varepsilon\zeta_{\max}^3}{\zeta_{\max} - \frac{\ln(1+\varepsilon\zeta_{\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_{\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^{\text{in}}, q^{\text{in}}) = (\zeta_c(x - x_0 - 0), c\zeta_c(x - x_0 - 0))$.

**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} (\zeta_{i+1}^n - 2\zeta_i^n + \zeta_{i-1}^n) \quad 1 \le i \le 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_{\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.

## 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
<figure><img src="image_5.png" /><figcaption>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$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>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$.</figcaption></figure>

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\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_{\text{ext}} \equiv 0$, $\langle q_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 ($\epsilon = 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,

$$ (\tau_{\kappa}(0)^2 + \ell\kappa)\ddot{\delta} + \ell K_{\kappa}^1 * \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 $K_{\kappa}^1$ is given in terms of the first Bessel function $J_1$ by the relation

$$ K_{\kappa}^{1}(t) = \frac{1}{t} J_{1}\left(\frac{t}{\kappa}\right). $$

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

$$ (56) \qquad \hat{\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, \quad 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_{\text{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.

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
<figure><img src="image_7.png" /><figcaption>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).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>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).</figcaption></figure>

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) \qquad \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><img src="image_9.png" /><figcaption>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$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>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).</figcaption></figure>

**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_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_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_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
<figure><img src="image_11.png" /><figcaption>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).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>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.</figcaption></figure>

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_i \rangle$ we exhibit an explicit solution of the equations in the linear case ($\epsilon = 0$). In that case, the wave-structure equations (4)-(7) take the form

$$ (58) \qquad \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 \text{in } \mathcal{E}^{\pm} $$

with transmission conditions

$$ (59) \qquad \llbracket q \rrbracket = 0 \quad \text{and} \quad \langle q \rangle = \langle q_i \rangle $$
and where $\langle q_i \rangle$ solves the forced ODE

$$ (60) \qquad \alpha(0) \frac{d}{dt} \langle q_i \rangle = -\frac{1}{2\ell} [\zeta + \kappa^2 \partial_t^2 \zeta]. $$

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 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

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

solve (58)-(59) with initial data

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

and with

$$ \langle q_i \rangle(t) = \omega [q^c \cos(\omega t) - \zeta^s \sin(\omega t)]. $$

If moreover

$$ q^c = -\frac{1}{2\ell\alpha(0)k}(q_+^s - q_-^s), \quad \text{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_i \rangle^{\text{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.

### 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_{\text{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
<figure><img src="image_13.png" /><figcaption>FIGURE 11. Interaction with a fixed object, linear case: convergence results for $\langle q_i \rangle$ with the second-order scheme, $3\kappa^2 = 0.1$ (left) and $3\kappa^2 = 0.3$ (right).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>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).</figcaption></figure>

convergence for $\langle q_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.

**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,
<figure><img src="image_15.png" /><figcaption>FIGURE 13. Interaction with a fixed object, non-linear case: mesh-convergence results for $\langle q_i \rangle$ with the first-order scheme, $3\kappa^2 = 0.1$ (left) and $3\kappa^2 = 0.3$ (right).</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIGURE 14. Interaction with a fixed object, non-linear case: mesh-convergence results for $\langle q_i \rangle$ with the second-order scheme, $3\kappa^2 = 0.1$ (left) and $3\kappa^2 = 0.3$ (right).</figcaption></figure>

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 is fixed.

## 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_i \rangle$ under the object; we have also shown that, in the augmented formulation, $\delta$ and $\langle q_i \rangle$ are found through the resolution of the 7-dimensional ODE (18) on $\Theta :=$
<figure><img src="image_17.png" /><figcaption>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$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>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).</figcaption></figure>

($\langle q_i \rangle, \dot{\delta}, \dot{\zeta}_+, \dot{\zeta}_-, \delta, \zeta_+, \zeta_-$)T that we wrote in abstract form as

$$ (61) \qquad \frac{d}{dt} \Theta = \mathcal{G}(\Theta, (R_1 f_{\text{sw}})_+, (R_1 f_{\text{sw}})_-, F_{\text{ext}}), $$

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_{\text{eq}}(x)$ is an even function.
<figure><img src="image_19.png" /><figcaption>FIGURE 17. Interaction with a freely floating object, non-linear case: mesh-convergence results for $\langle q_i \rangle$ with the first-order scheme, $3\kappa^2 = 0.1$ (left) and $3\kappa^2 = 0.3$ (right).</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>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).</figcaption></figure>

A.1. Explicit expressions of some coefficients. We did not make 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_{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 the water depth under the object at time $t$ is $h_{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_{eq})$$
<figure><img src="image_21.png" /><figcaption>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).</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>FIGURE 20. Interaction with a freely floating object, non-linear case: mesh-convergence results for $\langle q_i \rangle$ with the second-order scheme, $3\kappa^2 = 0.1$ (left) and $3\kappa^2 = 0.3$ (right).</figcaption></figure>

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

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

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_{\text{eq}} + \varepsilon\delta} dx + \kappa^2 \langle \frac{1}{h_{\text{eq}} + \varepsilon\delta} \rangle. $$
<figure><img src="image_23.png" /><figcaption>FIGURE 21. Interaction with a freely floating object, non-linear case: mesh-convergence results for $ζ_+$ with the second-order scheme, $3κ^2 = 0.1$ (left) and $3κ^2 = 0.3$ (right).</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>FIGURE 22. Interaction with a freely floating object, non-linear case: profile of solutions, $3κ^2 = 0.1$ (left) and $3κ^2 = 0.3$ (right), computed with the second-order scheme and $N = 400$.</figcaption></figure>

### A.2. The general case.
As seen in (14), the first four components of $Θ$ satisfy

$$ (62) \quad \mathcal{M}[\varepsilon\delta, \varepsilon\underline{\zeta}_{\pm}] \frac{d}{dt} \begin{pmatrix} \langle q_i \rangle \\ \dot{\delta} \\ \dot{\zeta}_{+} \\ \underline{\zeta}_{-} \end{pmatrix} + \begin{pmatrix} \frac{1}{2\ell}[[\zeta]] \\ \delta - \langle \zeta \rangle \\ \frac{\zeta_{+}}{\zeta_{-}} \\ \end{pmatrix} = \varepsilon \mathcal{Q} [\varepsilon\delta, \varepsilon\underline{\zeta}_{\pm}] (\langle q_i \rangle, \dot{\delta}, \underline{\zeta}_{\pm}) + \begin{pmatrix} 0 \\ F_{\text{ext}} \\ (R_1 f_{\text{sw}})_+ \\ (R_1 f_{\text{sw}})_- \end{pmatrix} $$

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

$$ (63) \quad \mathcal{M}[\varepsilon\delta, \varepsilon\underline{\zeta}_{\pm}] := \begin{pmatrix} \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}_-} \\ -\kappa & \ell\kappa & \kappa^2 & 0 \\ \kappa & \ell\kappa & 0 & \kappa^2 \end{pmatrix}, $$
with $\underline{h}_{\pm} = 1 + \varepsilon\underline{\zeta}_{\pm}$; simple computations also show that the quadratic term $\mathfrak{Q}[\varepsilon\delta, \varepsilon\underline{\zeta}_{\pm}]$ is of the form

$$ \mathfrak{Q}[\varepsilon\delta, \varepsilon\underline{\zeta}_{\pm}] = (\mathfrak{Q}_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}])^{\mathrm{T}} $$

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

$$ \mathfrak{Q}_i(\langle q_i \rangle, \dot{\delta}, \underline{\zeta}_{\pm}) = -\alpha'(\varepsilon\delta)\langle q_i \rangle \dot{\delta} - \frac{1}{4\ell} \left[ \left( \frac{-\ell\dot{\delta} + \langle q_i \rangle}{\underline{h}_+} \right)^2 - \left( \frac{\ell\dot{\delta} + \langle q_i \rangle}{\underline{h}_-} \right)^2 \right], $$

$$ \mathfrak{Q}_{\delta}(\langle q_i \rangle, \dot{\delta}, \underline{\zeta}_{\pm}) = \beta(\varepsilon\delta)\dot{\delta}^2 + \frac{1}{2}\alpha'(\varepsilon\delta)\langle q_i \rangle^2 + \frac{1}{4}\left[\left(\frac{-\ell\dot{\delta} + \langle q_i \rangle}{\underline{h}_+}\right)^2 + \left(\frac{\ell\dot{\delta} + \langle q_i \rangle}{\underline{h}_-}\right)^2\right], $$

$$ \mathfrak{Q}_{+}(\langle q_i \rangle, \dot{\delta}, \underline{\zeta}_{\pm}) = \frac{1}{2}\underline{\zeta}_{+}^2 - \frac{1}{\underline{h}_+}(-\ell\dot{\delta} + \langle q_i \rangle)^2, $$

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

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

$$ \begin{pmatrix} 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}_-} \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{pmatrix}; $$

the resulting equation takes the form

$$ (64) \quad \tilde{\mathcal{M}}[\varepsilon\delta, \varepsilon\underline{\zeta}_{\pm}] \frac{d}{dt} \begin{pmatrix} \langle q_i \rangle \\ \dot{\delta} \\ \dot{\underline{\zeta}}_+ \\ \dot{\underline{\zeta}}_- \end{pmatrix} + \begin{pmatrix} 0 \\ \delta \\ \underline{\zeta}_+ \\ \underline{\zeta}_- \end{pmatrix} = \varepsilon \tilde{\mathfrak{Q}}(\langle q_i \rangle, \dot{\delta}, \underline{\zeta}_{\pm}) + \begin{pmatrix} -\frac{1}{2\ell} \llbracket \frac{1}{\underline{h}} R_1 f_{\mathrm{sw}} \rrbracket \\ \langle \frac{1}{\underline{h}} R_1 f_{\mathrm{sw}} \rangle + F_{\mathrm{ext}} \\ (R_1 f_{\mathrm{sw}})_+ \\ (R_1 f_{\mathrm{sw}})_- \end{pmatrix} $$

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

$$ (65) \quad \tilde{\mathcal{M}}[\varepsilon\delta, \varepsilon\underline{\zeta}_{\pm}] := \begin{pmatrix} \alpha(\varepsilon\delta) + \frac{\kappa}{\ell} \langle \frac{1}{\underline{h}} \rangle & -\frac{\kappa}{2} \llbracket \frac{1}{\underline{h}} \rrbracket \\ -\frac{\kappa}{2} \llbracket \frac{1}{\underline{h}} \rrbracket & \tau_{\kappa}(\varepsilon\delta)^2 + \kappa\ell \langle \frac{1}{\underline{h}} \rangle \\ -\kappa & \ell\kappa \\ \kappa & \ell\kappa \end{pmatrix} \begin{pmatrix} 0 & 0 \\ 0 & 0 \\ \kappa^2 & 0 \\ 0 & \kappa^2 \end{pmatrix}, $$

and the components $\tilde{\mathfrak{Q}}_i, \tilde{\mathfrak{Q}}_{\delta}$ and $\tilde{\mathfrak{Q}}_{\pm}$ are given by

$$ \tilde{\mathfrak{Q}}_i(\langle q_i \rangle, \dot{\delta}, \underline{\zeta}_{\pm}) = -\frac{1}{4\ell} \llbracket \frac{1}{\underline{h}} \underline{\zeta}^2 \rrbracket + \frac{1}{4\ell} \llbracket \frac{1}{\underline{h}} \underline{\delta}^2 \rrbracket + \frac{1}{4\ell} \llbracket \frac{1}{\underline{h}} \underline{\zeta} \rangle \langle q_i \rangle^2 - (\alpha'(\varepsilon\delta) + \langle \frac{1}{\underline{h}^2} \rangle) \dot{\delta} \langle q_i \rangle, $$

$$ \tilde{\mathfrak{Q}}_{\delta}(\langle q_i \rangle, \dot{\delta}, \underline{\zeta}_{\pm}) = \frac{1}{2} \langle \frac{1}{\underline{h}} \underline{\zeta}^2 \rangle + (\beta(\varepsilon\delta) - \frac{1}{2}\ell^2 \langle \frac{1}{\underline{h}^2} \rangle) \dot{\delta}^2 + \frac{1}{2}(\alpha'(\varepsilon\delta) - \langle \frac{1}{\underline{h}^2} \rangle) \langle q_i \rangle^2 + \frac{1}{2}\ell \llbracket \frac{1}{\underline{h}^2} \rrbracket \dot{\delta} \langle q_i \rangle, $$
and

$$
\begin{align*}
\tilde{\mathcal{Q}}_+(\langle q_i \rangle, \dot{\delta}, \underline{\zeta}_{\pm}) &= -\frac{1}{2}\underline{\zeta}_{\pm}^2 - \frac{1}{\underline{h}_+} \langle q_i \rangle^2 - \ell^2 \frac{1}{\underline{h}_+} \dot{\delta}^2 + 2\ell \frac{1}{\underline{h}_+} \langle q_i \rangle \dot{\delta}, \\
\tilde{\mathcal{Q}}_-(\langle q_i \rangle, \dot{\delta}, \underline{\zeta}_{\pm}) &= -\frac{1}{2}\underline{\zeta}_{\pm}^2 - \frac{1}{\underline{h}_-} \langle q_i \rangle^2 - \ell^2 \frac{1}{\underline{h}_-} \dot{\delta}^2 - 2\ell \frac{1}{\underline{h}_-} \langle q_i \rangle \dot{\delta}.
\end{align*}
$$

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

$$
(66) \quad \frac{d}{dt} (\langle q_i \rangle, \dot{\delta}, \underline{\zeta}_+, \underline{\zeta}_-)^{\mathrm{T}} = \mathcal{G}_I(\Theta, (R_1 f_{\mathrm{sw}})_{\pm}, F_{\mathrm{ext}})
$$

with the notation $\mathcal{G} := (\mathcal{G}_1, \mathcal{G}_2, \mathcal{G}_3, \mathcal{G}_4)^{\mathrm{T}}$; we deduce from (64) that

$$
\begin{align*}
& \mathcal{G}_I(\Theta, (R_1 f_{\text{sw}})_{\pm}, F_{\text{ext}}) \\
&= \tilde{\mathcal{M}}[\varepsilon\delta, \varepsilon\underline{\zeta}_{\pm}]^{-1} \left[ - \begin{pmatrix} 0 \\ \delta \\ \underline{\zeta}_{+} \\ \underline{\zeta}_{-} \end{pmatrix} + \varepsilon \tilde{\mathcal{Q}}(\langle q_i \rangle, \dot{\delta}, \underline{\zeta}_{\pm}) + \begin{pmatrix} -\frac{1}{2\ell} \left[ \frac{1}{\underline{h}} R_1 f_{\text{sw}} \right] \\ \left\langle \frac{1}{\underline{h}} R_1 f_{\text{sw}} \right\rangle + F_{\text{ext}} \\ (R_1 f_{\text{sw}})_+ \\ (R_1 f_{\text{sw}})_- \end{pmatrix} \right],
\end{align*}
$$

while

$$
\mathcal{G}_5(\Theta) = \dot{\delta}, \quad \mathcal{G}_6(\Theta) = \dot{\underline{\zeta}}_+, \quad \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 $\tilde{\mathcal{M}}[\varepsilon\delta, \varepsilon\underline{\zeta}_{\pm}]$, namely,

$$
(67) \quad \tilde{\mathcal{M}}[\varepsilon\delta, \varepsilon\underline{\zeta}_{\pm}]^{-1} = \begin{pmatrix} -\frac{4}{D}(\tau_{\kappa}(\varepsilon\delta)^2 + \kappa\ell\langle\frac{1}{\underline{h}}\rangle) & -\frac{2\kappa}{D}\left[\frac{1}{\underline{h}}\right] \\ -\frac{2\kappa}{D}\left[\frac{1}{\underline{h}}\right] & -\frac{4}{D}(\alpha(\varepsilon\delta) + \frac{\kappa}{\ell}\langle\frac{1}{\underline{h}}\rangle) \\ -\frac{4}{\kappa D}(\tau_{\kappa}(\varepsilon\delta)^2 + \kappa\ell\langle\frac{1}{\underline{h}_-}\rangle) & \frac{4}{\kappa D}(\kappa\langle\frac{1}{\underline{h}_-}\rangle + \ell\alpha(\varepsilon\delta)) \\ \frac{4}{\kappa D}(\tau_{\kappa}(\varepsilon\delta)^2 + \kappa\ell\langle\frac{1}{\underline{h}_+}\rangle) & \frac{4}{\kappa D}(\kappa\langle\frac{1}{\underline{h}_+}\rangle + \ell\alpha(\varepsilon\delta)) \end{pmatrix} \begin{pmatrix} 0 & 0 \\ 0 & 0 \\ \frac{1}{\kappa^2} & 0 \\ 0 & \frac{1}{\kappa^2} \end{pmatrix}
$$

with

$$
(68) \quad D = -4 \left( \alpha(\varepsilon\delta) + \frac{\kappa}{\ell} \langle \frac{1}{\underline{h}} \rangle \right) \times \left( \tau_{\kappa}(\varepsilon\delta)^2 + \kappa \ell \langle \frac{1}{\underline{h}} \rangle \right) + \kappa^2 \left[ \frac{1}{\underline{h}} \right]^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_{\text{forced}}$ ($\delta_{\text{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 := (\langle q_i \rangle, \dot{\zeta}_+, \dot{\zeta}_-, \zeta_+, \zeta_-)^{\mathrm{T}}$ and (64) can be simplified into

$$
\begin{align}
\tilde{\mathcal{M}}_{\text{forced}}[t, \varepsilon \underline{\zeta}_{\pm}] \frac{d}{dt} \begin{pmatrix} \langle q_i \rangle \\ \dot{\zeta}_+ \\ \dot{\zeta}_- \end{pmatrix} + \begin{pmatrix} 0 \\ \dot{\zeta}_+ \\ \dot{\zeta}_- \end{pmatrix} &= \varepsilon \tilde{\mathcal{Q}}_{\text{forced}}[t, \varepsilon \underline{\zeta}_{\pm}] (\langle q_i \rangle, \underline{\zeta}_{\pm}) \nonumber \\
(69) \qquad & \qquad + \begin{pmatrix} -\frac{1}{2\ell} \left[ \frac{1}{\underline{h}} R_1 f_{\text{sw}} \right] \\ (R_1 f_{\text{sw}})_+ \\ (R_1 f_{\text{sw}})_- \end{pmatrix} + \begin{pmatrix} \frac{1}{2}\kappa \left[ \frac{1}{\underline{h}} \right] \\ -\ell\kappa \\ -\ell\kappa \end{pmatrix} \ddot{\delta}_{\text{forced}},
\end{align}
$$
where the matrix $\tilde{\mathcal{M}}_{\text{forced}}[t, \varepsilon\underline{\zeta}_{\pm}]$ is block-triangular,

$$ (70) \qquad \tilde{\mathcal{M}}_{\text{forced}}[t, \varepsilon\underline{\zeta}_{\pm}] := \begin{pmatrix} \alpha(\varepsilon\delta_{\text{forced}}) + \frac{\kappa}{\ell}\langle\frac{1}{h}\rangle & 0 \\ -\kappa & \kappa^2 \\ \kappa & 0 \end{pmatrix}, $$

and

$$ \tilde{\mathfrak{Q}}_{\text{forced}}[t, \varepsilon\underline{\zeta}_{\pm}] (\langle q_i \rangle, \underline{\zeta}_{\pm}) := \begin{pmatrix} \tilde{\mathfrak{Q}}_i[\varepsilon\delta_{\text{forced}}, \varepsilon\underline{\zeta}_{\pm}] (\langle q_i \rangle, \dot{\delta}_{\text{forced}}, \underline{\zeta}_{\pm}) \\ \tilde{\mathfrak{Q}}_+ [\varepsilon\delta_{\text{forced}}, \varepsilon\underline{\zeta}_{\pm}] (\langle q_i \rangle, \dot{\delta}_{\text{forced}}, \underline{\zeta}_{\pm}) \\ \tilde{\mathfrak{Q}}_- [\varepsilon\delta_{\text{forced}}, \varepsilon\underline{\zeta}_{\pm}] (\langle q_i \rangle, \dot{\delta}_{\text{forced}}, \underline{\zeta}_{\pm}) \end{pmatrix}, $$

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

**Remark A.2.** We have made explicit the dependence of $\tilde{\mathcal{M}}_{\text{forced}}$ and $\tilde{\mathfrak{Q}}_{\text{forced}}$ on the time variable $t$ because $\delta_{\text{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_{\text{forced}} \equiv 0$).

With $\mathcal{G}_1$ now being three dimensional (but with an extra dependence on $t$), $\mathcal{G}_1 := (\mathcal{G}_1, \mathcal{G}_2, \mathcal{G}_3)^T$, we have therefore

$$ (71) \qquad \frac{d}{dt} (\langle q_i \rangle, \dot{\underline{\zeta}}_+, \dot{\underline{\zeta}}_-)^T = \mathcal{G}_1(t, \Theta, (R_1 f_{\text{sw}})_\pm) $$

and

$$ \mathcal{G}_1(t, \Theta, (R_1 f_{\text{sw}})_\pm) = \tilde{\mathcal{M}}_{\text{forced}}[t, \varepsilon\underline{\zeta}_\pm]^{-1} \times \left[ - \begin{pmatrix} 0 \\ \underline{\zeta}_+ \\ \underline{\zeta}_- \end{pmatrix} + \varepsilon \tilde{\mathfrak{Q}}_{\text{forced}}[t, \varepsilon\underline{\zeta}_\pm] (\langle q_i \rangle, \underline{\zeta}_\pm) + \begin{pmatrix} -\frac{1}{2\ell} \llbracket \frac{1}{h} R_1 f_{\text{sw}} \rrbracket \\ (R_1 f_{\text{sw}})_+ \\ (R_1 f_{\text{sw}})_- \end{pmatrix} + \begin{pmatrix} \frac{1}{2} \kappa \llbracket \frac{1}{h} \rrbracket \\ -\ell\kappa \\ -\ell\kappa \end{pmatrix} \ddot{\delta}_{\text{forced}} \right], $$

while $\mathcal{G}_4(\Theta) = \underline{\zeta}_+$, $\mathcal{G}_5(\Theta) = \underline{\zeta}_-$ and

$$ \tilde{\mathcal{M}}_{\text{forced}}[t, \varepsilon\underline{\zeta}_{\pm}]^{-1} = \begin{pmatrix} \frac{1}{\alpha + \frac{\kappa}{\ell} \langle \frac{1}{h} \rangle} & 0 & 0 \\ \frac{1}{\kappa \alpha + \frac{\kappa}{\ell} \langle \frac{1}{h} \rangle} & \frac{1}{\kappa^2} & 0 \\ -\frac{1}{\kappa \alpha + \frac{\kappa}{\ell} \langle \frac{1}{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_{\text{forced}}$? The answer is explicitly given (writing $\mathfrak{Q}_\delta = \mathfrak{Q}_\delta[\varepsilon\delta_{\text{forced}}, \varepsilon\underline{\zeta}_\pm]$, etc.)

$$ F_{\text{ext}} = \delta - \left\langle \frac{1}{h} R_1 f_{\text{sw}} \right\rangle - \varepsilon \mathfrak{Q}_{\delta} (\langle q_i \rangle, \dot{\delta}_{\text{forced}}, \zeta_{\pm}) \\ - \frac{1}{4} \frac{D_{\text{forced}}}{\alpha(\varepsilon\delta_{\text{forced}}) + \frac{\kappa}{\ell} \langle \frac{1}{h} \rangle} \ddot{\delta}_{\text{forced}} \\ - \frac{1}{2} \frac{\kappa}{\alpha(\varepsilon\delta_{\text{forced}}) + \frac{\kappa}{\ell} \langle \frac{1}{h} \rangle} \llbracket \frac{1}{h} \rrbracket [\varepsilon \mathfrak{Q}_i (\langle q_i \rangle, \dot{\delta}_{\text{forced}}, \zeta_{\pm}) - \frac{1}{2\ell} \llbracket \frac{1}{h} R_1 f_{\text{sw}} \rrbracket], $$
where $D_{\text{forced}}$ is deduced from the expression given for $D$ in (68) by substituting $\delta_{\text{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_{\text{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_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) \qquad \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 \text{for } t > 0, x \in (\ell, \infty) $$

with boundary condition

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

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

$$ (74) \qquad \frac{d}{dt} \Theta = \mathcal{G}(\Theta, (R_1 f_{\text{sw}})_+, F_{\text{ext}}), $$

with $\Theta := (\dot{\delta}, \dot{\zeta}_+, \delta, \zeta_+)^{\text{T}}$ and where $\mathcal{G} = (\mathcal{G}_1, \mathcal{G}_2, \mathcal{G}_3, \mathcal{G}_4)^{\text{T}}$ is given by

$$ \left( \begin{array}{c} \mathcal{G}_1 \\ \mathcal{G}_2 \end{array} \right) = \mathcal{M}_{\text{sym}}[\varepsilon\delta, \varepsilon\zeta_+]^{-1} \left[ \begin{array}{l} - \left( \begin{array}{c} \delta \\ \zeta_+ \end{array} \right) + \varepsilon \left( \begin{array}{c} \mathcal{Q}_\delta^{\text{sym}}[\varepsilon\delta, \varepsilon\zeta_+](\dot{\delta}, \zeta_+) \\ \mathcal{Q}_+^{\text{sym}}[\varepsilon\delta, \varepsilon\zeta_+](\dot{\delta}, \zeta_+) \end{array} \right) \\ + \left( \begin{array}{c} \frac{1}{h_+} (R_1 f_{\text{sw}})_+ + F_{\text{ext}} \\ (R_1 f_{\text{sw}})_+ \end{array} \right) \end{array} \right] $$

and $\mathcal{G}_3 = \dot{\delta}$, $\mathcal{G}_4 = \dot{\zeta}_+$, and with

$$ \mathcal{M}_{\text{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 $\mathcal{Q}_\delta^{\text{sym}}$ and $\mathcal{Q}_+^{\text{sym}}$ are obtained by replacing $\langle q_i \rangle = 0$ and $\zeta_- = \zeta_+$ in the formula derived above for $\mathcal{Q}_\delta$ and $\mathcal{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_{\text{ext}} = 0$ and initial conditions $(\zeta, q)|_{t=0} = (0, 0)$ and $(\dot{\delta}, \dot{\zeta}_+, \delta, \zeta_+)|_{t=0} = (0, 0, \delta^{\text{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_{\text{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}_{\text{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. Tucnak, 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 equa- tions, *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.