On the return to equilibrium problem for axisymmetric floating structures in shallow water

TL;DR

This paper investigates the return to equilibrium problem for axisymmetric floating structures in shallow water, deriving a nonlinear integro-differential equation governing the solid motion and establishing global existence and uniqueness of solutions.

Contribution

It introduces a reduced delay differential equation involving an extension-trace operator and derives explicit forms under small amplitude wave assumptions.

Findings

The solid motion follows a nonlinear second order integro-differential equation.

The linearized model corresponds to the Cummins equation.

Global existence and uniqueness of solutions are proven using energy conservation.

Abstract

In this paper we address the return to equilibrium problem for an axisymmetric floating structure in shallow water. First we show that the equation for the solid motion can be reduced to a delay differential equation involving an extension-trace operator whose role is to describe the influence of the fluid equations on the solid motion. It turns out that the compatibility conditions on the initial data for the return to equilibrium configuration are not satisfied, so we cannot use the result from [3] for the nonlinear problem. Hence, assuming small amplitude waves, we linearize the equations in the exterior domain and we keep the nonlinear equations in the interior domain. For such configurations, the extension-trace operator can be computed explicitly and the delay term in the differential equation can be put in convolution form. The solid motion is therefore governed by a nonlinear second order integro-differential equation, whose linearization is the well-known Cummins equation. We show global in time existence and uniqueness of the solution using the conservation of the total fluid-structure energy.

Equations

$$\left(m+a_{\infty }\right){\ddot{\delta }}_G(t)=-c{\delta }_G(t)-{\int }_0^{t}K(s){\dot{\delta }}_G(t-s)ds,$$

$$(m+m_a({\delta }_G)){\ddot{\delta }}_G(t)=-\mathfrak{c}{\delta }_G(t)-\nu ({\dot{\delta }}_G)+\beta ({\delta }_G){\dot{\delta }}_G^2(t),$$

$$(m+m_a({\delta }_G)){\ddot{\delta }}_G=$$

$$\Omega (t)=\{(r,z)\in {\mathbb{R}}_{+}\times \mathbb{R}|-h_0<z<\zeta (t,r)\}.$$

$$f_i:=f_{{|}_{r\in (0,R)}}\text{mbox}and{f}_e:=f_{{|}_{r\in (R,+\infty )}}.$$

$${\zeta }_i(t,r)={\zeta }_w(t,r).$$

$$h_e(0,r)\ge h_m\text{mbox}forr\in (R,+\infty ).$$

$$h_w(0,r)\ge h_{min}\text{mbox}forr\in (0,R),$$

$$\{\begin{array}{l}{\partial }_{t}h+{\partial }_{r}q+\frac{q}{r}=0,\\ {\partial }_{t}q+{\partial }_r\left(\frac{q^2}{h}\right)+\frac{q^2}{rh}+gh{\partial }_{r}h=-\frac{h}{\rho }{\partial }_r\underline{P}\end{array}$$

$$q_{e_{{|}_{r=R}}}=q_{i_{{|}_{r=R}}}.$$

$$\{\begin{array}{l}-\left({\partial }_r+\frac{1}{r}\right)\left(\frac{h_w}{\rho }{\partial }_r{\underline{P}}_i\right)=\left({\partial }_r+\frac{1}{r}\right)\left({\partial }_r\left(\frac{q_i^2}{h_w}\right)+\frac{q_i^2}{r{h}_w}+g{h}_w{\partial }_r{h}_w\right)-{\dot{w}}_G,\\ {{\underline{P}}_i}_{{|}_{r=R}}=P_{\mathrm{atm}}+\rho g({\zeta }_e-{\zeta }_i{)}_{{|}_{r=R}}+P_{\mathrm{cor}},\end{array}$$

$$m{\ddot{\delta }}_G(t)=-mg+{\int }_{r\le R}\left({\underline{P}}_i-P_{\mathrm{atm}}\right).$$

$$\{\begin{array}{l}{\partial }_t{h}_e+{\partial }_r{q}_e+\frac{q_e}{r}=0,\\ {\partial }_t{q}_e+{\partial }_r\left(\frac{q_e^2}{h_e}\right)+\frac{q_e^2}{r{h}_e}+g{h}_e{\partial }_r{h}_e=0,\\ {q_e}_{{|}_{r=R}}=-\frac{R}{2}{\dot{\delta }}_G,\end{array}$$

$$(m+m_a({\delta }_G)){\ddot{\delta }}_G(t)=-\mathfrak{c}{\delta }_G(t)+\mathfrak{c}{\zeta }_e(t,R)+\left(\mathfrak{b}(h_e)+\beta ({\delta }_G)\right){\dot{\delta }}_G^2(t)$$

$$\mathfrak{c}=\rho g\pi R^2,\mathfrak{b}(h_e)=\frac{b}{h_e^2(t,R)}\text{mbox}withb=\frac{\pi \rho R^4}{8},$$

$$m_a({\delta }_G)=\frac{\rho \pi }{2}{\int }_0^R\frac{r^3}{h_w({\delta }_G,r)}dr,$$

$$\beta ({\delta }_G)=\frac{b}{2h_w^2({\delta }_G,R)}+\frac{\pi \rho }{8}{\int }_0^R\frac{r^4}{h_w^3({\delta }_G,r)}{\partial }_r{h}_w({\delta }_G,r)dr,$$

$$q_i(t,r)=-\frac{r}{2}{\dot{\delta }}_G(t).$$

$$\{\begin{array}{l}{\partial }_{t}u+A(u){\partial }_{r}u+B(u,r)u=0,\\ {q_e}_{{|}_{r=R}}=-\frac{R}{2}{\dot{\delta }}_G,\\ u(0)=u_0,\end{array}$$

$$A(u)=\left(\begin{array}[]{cc}0&1\\ gh_{e}-\dfrac{q_{e}^{2}}{h_{e}^{2}}&\dfrac{2q_{e}}{h_{e}}\\ \end{array}\right),\ \ \ \ B(u,r)=\left(\begin{array}[]{cc}0&\dfrac{1}{r}\\[10.0pt] 0&\dfrac{q_{e}}{rh_{e}},\\ \end{array}\right)$$

$$u_0=({\zeta }_{e,0},q_{e,0}{)}^T.$$

$$X^k(T):=\bigcap _{j=0}^k{C}^j([0,T],H_r^{k-j}((R,+\infty )))$$

$$\Vert u{\Vert }_{X^k(T)}:=\underset{[0,T]}{\sup }\Vert u(t){\Vert }_{X^k},\Vert u(t){\Vert }_{X^k}=\sum _{j=0}^k\Vert {\partial }_t^{j}u(t){\Vert }_{H_r^{k-j}((R,+\infty ))},$$

$$\Vert u(t){\Vert }_{X^k}^2+\Vert u_{{|}_{r=R}}{\Vert }_{H^k((0,t))}^2\le C\left(T,\Vert u_0{\Vert }_{H_r^k((R,\infty ))}^2,\Vert {q_e}_{{|}_{r=R}}{\Vert }_{H^k((0,t))}^2\right)$$

$$\mathcal{B}:H^k((0,T))\times H_r^k((R,\infty ))$$

$$({\dot{\delta }}_G,u_0)$$

$$(m+m_a({\delta }_G)){\ddot{\delta }}_G(t)$$

$$=-\mathfrak{c}{\delta }_G(t)+\mathfrak{c}\mathcal{B}\left[{\dot{\delta }}_G,u_0\right](t)+\left(\frac{b}{(\mathcal{B}\left[{\dot{\delta }}_G,u_0\right](t)+h_0{)}^2}+\beta ({\delta }_G)\right){\dot{\delta }}_G^2(t).$$

$${\mathbf{e}}_2\cdot \text{mbox}\text{textquotedblleft}{\partial }_t^j{u}_{{|}_{t=0}}{\text{textquotedblright}}_{{|}_{r=R}}=-\frac{R}{2}\text{mbox}\text{textquotedblleft}\frac{d^{j+1}}{d{t}^{j+1}}{{\delta }_G}_{{|}_{t=0}}\text{textquotedblright},$$

$$\text{mbox}\text{textquotedblleft}{\partial }_t^1{u}_{{|}_{t=0}}\text{textquotedblright}=-A(u_0){\partial }_r{u}_0-B(u_0,r)u_0,$$

$$\mbox{\textquotedblleft}\frac{d^{2}}{dt^{2}}{\delta_{G}}_{|_{t=0}}\textquotedblright=\frac{-\mathfrak{c}\delta_{0}+\mathfrak{c}\mathcal{B}\left[\dot{\delta}_{G},u_{0}\right]_{|_{t=0}}+\bigg{(}\frac{b}{(\mathcal{B}\left[\dot{\delta}_{G},u_{0}\right]_{|_{t=0}}+h_{0})^{2}}+\beta(\delta_{0})\bigg{)}\delta_{1}^{2}}{m+m_{a}(\delta_{0})}.$$

$$\text{mbox}forr\in (R,+\infty )h_{e,0}(r)\ge h_{\mathrm{m}},\left(g{h}_{e,0}-\frac{q_{e,0}^2}{h_{e,0}^2}\right)(r)\ge c_{\mathrm{sub}},$$

Full text

1 Introduction

The return to equilibrium problem is a particular configuration of the floating structure problem. It consists in releasing a partially submerged solid body in a fluid initially at rest and letting it evolve towards its equilibrium position. The interest of this problem is that it can easily be done experimentally and it is used in engineering to determine several important characteristics of floating objects. More precisely, engineers assume that the solid satisfies a linear integro-differential equation, the Cummins equation (see [4]). The experimental data coming from the return to equilibrium problem (also called decay test) are then used to identify the coefficients of this linear equation. John in [8] studied the problem in shallow water in one horizontal dimension for an object with flat bottom: he considered the linearized fluid equations for small amplitude waves and he wrote an explicit expression for the solid motion under linear approximation. Ursell in [18] and Maskell and Ursell [12], using like John the linear approach, obtained an explicit solution in integral form for the vertical displacement of the object. Still under the linear approximation Cummins in [4] treated a general ship motion and reduced the free motion of the floating body to an integro-differential equation. From Wehausen and Laitone [21] we know that also Sretenskii, several years before Cummins, obtained an integro-differential equation for the vertical displacement which he solved numerically. The Cummins equation for the vertical displacement reads

[TABLE]

where $\delta_{G}(t)=z_{G}(t)-z_{G,eq}$ is the displacement from the equilibrium position of the vertical position of the centre of mass, $m$ is the mass of the structure, $a_{\infty}$ is the added mass at infinity frequency, $c$ is the hydrostatic coefficient and $K$ is the impulse response function (also known as retardation function and fluid memory). It appears in naval architecture and hydrodynamical engineering literature and it is used to study the motion of ships or wave energy converters. Recently Lannes in his paper [9] on the dynamics of floating structures modelled the return to equilibrium problem using a different formulation for the hydrodynamical model with the aim to take into account nonlinear effects. He wrote the explicit equations in the one-dimensional (horizontal) case and, considering the nonlinear shallow water model, he showed that the position of the solid is fully determined by the nonlinear second order damped ODE

[TABLE]

where $m_{a}(\delta_{G})$ is the nonlinear added mass and $\nu(\dot{\delta}_{G})$ is the nonlinear damping term. Numerical simulations for the one dimensional model proposed by Lannes are made in [19].

In our recent paper [3] we dealt with the two-dimensional (horizontal) case, we showed the local well-posedness for the axisymmetric floating structure problem in the shallow water regime for initial data regular enough, provided some compatibility conditions are satisfied. We considered a solid, with vertical side-walls and a cylindrical symmetry, forced to move only vertically. For such a configuration, the horizontal coordinates of the contact line between the air, the fluid and the solid, are time independent. For an object with no vertical walls, finding the horizontal coordinates of the contact line is a free boundary problem, recently solved in the one horizontal dimension case by Iguchi and Lannes in [7] where the contact line is replaced by two contact points. The floating structure problem for a viscous fluid in a one dimensional bounded domain is considered in [11].

The aim of this paper is to extend the work of Lannes on the return to equilibrium problem to the two-dimensional case taking into account nonlinear effects and using the same framework as in [3], which means that we consider here the axisymmetric setting, the shallow water approximation for the fluid and a solid with the properties we have stated before. An important change with respect to the one-dimensional case is the presence of delay terms in the equation governing the solid motion. The nonlinear coupled system can be treated in an abstract way but, as we show here, it requires compatibility conditions that are not satisfied in the return to equilibrium problem. For this reason, we linearize the equations in the exterior domain but we keep the nonlinear effects in the interior domain. This approach permits us to improve the classical linear model and we get a nonlinear second order delay differential equation on the vertical displacement of the structure. If we linearize around the equilibrium state we get the standard linear Cummins equation, hence we can see the result of this paper as a rigorous justification and an extension of Cummins’ work.

1.1 Outline of the paper

In Section 2 we first recall the notations for the floating structures that we have used in [3] and we write the equations for the coupled problem. Then we show that the differential equation for the solid motion can be written in a closed form by introducing an extension-trace operator, which takes the boundary value of the horizontal discharge, defined as the fluid horizontal velocity vertically integrated, in the exterior domain and gives the boundary value of the fluid height in the exterior domain. In Theorem 2.2 we solve the equation by a fixed point argument. Finally we consider the return to equilibrium configuration, giving the initial conditions on the fluid and solid unknowns. It turns out that the compatibility conditions, which are necessary in order to apply the existence and uniqueness theorem from [3], are not satisfied for these particular initial conditions.

In Section 3 we neglect the nonlinear effects in the exterior region, but we keep them under the object provided it does not touch the bottom of the fluid domain. We write a linear-nonlinear model for the floating structure problem: we linearize the equations in the exterior domain and we keep the nonlinearities in the interior domain. Hence the equations for the fluid in the exterior domain become the linear shallow water equations and the free surface elevation in the exterior domain satisfies a wave equation. Then, by applying a Fourier-Laplace transform argument, we can write the trace of the exterior free surface elevation at the boundary as a convolution product between the inverse Laplace transform of a Hankel function and the time derivative of the displacement $\delta_{G}$. Hence we have that the solid motion is described by the nonlinear integro-differential equation

[TABLE]

Differently from (2), the damping term is a linear function of $\dot{\delta}_{G}$ since the equations in the exterior domain are linear and it is given by a delay term due to dispersion which does not occur in one dimension. Its linearization around the equilibrium gives a reformulation of the Cummins equation for the vertical displacement (1). We show in Theorem 47 the global existence and uniqueness of its solution, provided an admissibility condition for the initial datum, using the conservation of the total fluid-structure energy. In Section 4 we explain the numerical method we use to plot the time evolution of the vertical displacement of the structure for the return to equilibrium problem. We compare the numerical solution to the nonlinear integro-differential equation with the solution to the linear Cummins equation and we note that for large initial data the nonlinear effects should not be neglected. In Appendix A we define the Hankel functions and we show some properties and results.

Acknowledgements

The author warmly thanks Kai Koike for his precious advises on the long time behaviour of the impulse response function and Pierre Magal for his remarks on functional differential equations with infinite delay.

The author was supported by the grants ANR-17-CE40-0025 NABUCO and ANR-18-CE40-0027 SingFlows, the Junior Chair BOLIDE of the IDEX of the University of Bordeaux, the Simone and Cino Del Duca Foundation and the Nouvelle Aquitaine Regional Council.

2 Nonlinear floating structure equations

Let us recall the following notations: $\zeta(t,r)$ is the elevation of the free surface, $h(t,r)=\zeta(t,r)+h_{0}$ is the fluid height, $q(t,r)$ is the horizontal discharge, i.e. the radial component of the fluid velocity vertically integrated, $\underline{P}$ is the trace of the pressure at the surface and $\zeta_{w}(t,r)$ is the parametrization of the bottom of the solid. The centre of mass of the solid is $G(t)=(0,0,z_{G}(t))$ and its velocity is $\mathbf{U}_{G}(t)=(0,0,w_{G}(t))$. We define $\delta_{G}(t)=z_{G}(t)-z_{G,eq}$ the displacement from the equilibrium of the vertical position of the centre of mass. We denote by $\rho_{m}$ the density of the floating body and $H$ its height. The fluid domain is

[TABLE]

Moreover, as shown in Figure 1, the presence of the solid permits us to divide the radial line in two regions, the interior domain $(0,R)$ and the exterior domain $(R,+\infty)$, whose boundary is the projection $r=R$ of the contact line between the fluid, the air and the body. Throughout all the paper we will note, for a function $f(r)$,

[TABLE]

We have the contact constraint in the interior domain

[TABLE]

As in the standard water waves theory we assume that the initial height of the fluid in the exterior domain does not vanish, i.e. there exists $h_{m}>0$ such that

[TABLE]

For the sake of the problem, we suppose also that the solid does not touch the bottom of the domain during its motion. Hence we assume that the initial height of the fluid under the solid does not vanish, i.e. there exists $h_{min}>0$ such that

[TABLE]

with $h_{w}(0,r)=h_{i}(0,r)$ in the interior domain due to (4).

We showed in [3] that the floating structure problem in the case of an axisymmetric flow without swirl is described by

[TABLE]

coupled with the transition condition

[TABLE]

We have $\underline{P}_{e}=P_{\mathrm{atm}}$, where $P_{\mathrm{atm}}$ is the constant atmospheric pressure, while $\underline{P}_{i}$ is given by the following elliptic problem in the interior domain $(0,R)$:

[TABLE]

with $P_{\mathrm{cor}}=\dfrac{\rho}{2}q^{2}_{i_{|_{r=R}}}\left(\dfrac{1}{h^{2}_{e}}_{|_{r=R}}-\dfrac{1}{h^{2}_{i}}_{|_{r=R}}\right)$. We replace $h_{i}=\zeta_{i}+h_{0}$ with $h_{w}=\zeta_{w}+h_{0}$ due to the contact constraint (4). The boundary condition on the pressure is chosen in order to have exact conservation of the energy for the fluid-solid system (see [3]).

The free motion of the solid is described by Newton’s law for the conservation of the linear momentum

[TABLE]

Using the elliptic equation (9) we can formulate the floating structure problem in the axisymmetric case as the following coupled problem (for details see [3]):

• •

the quasilinear hyperbolic boundary problem for the fluid motion in the exterior domain

[TABLE]

• •

Newton’s equation for the conservation of the linear momentum can be put under the form

[TABLE]

with

[TABLE]

with $h_{w}(\delta_{G},r)=h_{w,eq}(r)+\delta_{G}(t)$. Due to this decomposition of $h_{w}$ and the contact constraint (4), we get the boundary condition in (10) from (8) and the explicit resolution of the first equation in (7) in the interior domain, that is

[TABLE]

The term $h_{w,eq}(r)$ is the fluid height under the solid at the equilibrium position and $\zeta_{w,eq}(r)=h_{w,eq}(r)-h_{0}$ is the elevation of the bottom of the solid at the equilibrium position. They both depend on the density of the fluid $\rho$, the density of the solid $\rho_{m}$, the depth $h_{0}$ and the height of the solid $H$ (see Section 3 for the explicit expressions in the flat bottom case).

2.1 Extension-trace operator for the coupling with the exterior domain

In this section we want to show that, in the ODE for the solid part of the coupled system (10) - (11), we can write the coupling term $\zeta_{e}(t,R)$ (also $h_{e}^{2}(t,R)$), the trace of the free surface elevation in the exterior domain at the boundary $r=R$, as an extension-trace operator applied to the trace of the horizontal discharge in the interior domain at the boundary $r=R$, that is $-\frac{R}{2}\dot{\delta}_{G}$.

We consider the exterior quasilinear hyperbolic initial boundary value problem (10) and using $u=(\zeta_{e},q_{e})^{T}$ we can write it as

[TABLE]

with

[TABLE]

and

[TABLE]

We consider the functional space

[TABLE]

endowed with the norm

[TABLE]

where $H^{k}_{r}:=H^{k}(rdr)$ is the weighted Sobolev space. In Theorem 5.3 of [3] we showed that, for $k\geq 2$, there exists $T>0$ and a unique solution $u=(\zeta_{e},q_{e})^{T}\in X^{k}(T)$ to (13), provided the initial data $u_{0}\in H^{k}_{r}((R,+\infty))$, the boundary datum ${q_{e}}_{|_{r=R}}\in H^{k}((0,T))$ and compatibility conditions are satisfied up to order $k-1$. Moreover $u$ satisfies the following energy estimate:

[TABLE]

for all $t\in(0,T)$. Then we can define an operator $\mathcal{B}$ such that

[TABLE]

We call it an extension-trace operator since it takes the trace of $q_{e}$, that is $-\frac{R}{2}\dot{\delta}_{G}$, the initial data $u_{0}$ and it extends to the couple $(\zeta_{e},q_{e})$ by solving the initial boundary value problem (13) and then it takes the trace of $\zeta_{e}$. One can easily note that $\mathcal{B}$ is nonlinear. Then, using the fact that $h_{e}=\zeta_{e}+h_{0}$ and assuming $u_{0}$ to be given, we can write the equation (11) for the solid motion as a second order delay differential equation only in terms of $\delta_{G}$, namely

[TABLE]

It is a delay differential equation since we need to know $\dot{\delta}_{G}$ for all $t^{\prime}\in[0,t]$ in order to know the value of $\mathcal{B}\left[\dot{\delta}_{G},u_{0}\right]$ at time $t$. This equation can be solved by a standard fixed point argument. Let us recall the compatibility conditions on the initial data (see [16]):

Definition 2.1**.**

The data $u_{0}\in H_{r}^{k}((R,+\infty))$, $\delta_{0}\in\mathbb{R}$ and $\delta_{1}\in\mathbb{R}$ of the floating structure coupled system (13) - (17) satisfy the compatibility conditions up to order $k-1$ if, for $0\leq j\leq k-1$, the following holds:

[TABLE]

where $\mbox{\textquotedblleft}u_{|_{t=0}}\textquotedblright=u_{0}$, $\mbox{\textquotedblleft}\frac{d}{dt}{\delta_{G}}_{|_{t=0}}\textquotedblright=\delta_{1}$, and for $j\geq 1$ $\mbox{\textquotedblleft}\partial_{t}^{j}u_{|_{t=0}}\textquotedblright$ and $\mbox{\textquotedblleft}\frac{d^{j+1}}{dt^{j+1}}{\delta_{G}}_{|_{t=0}}\textquotedblright$ are inductively defined by formally taking $j-1$ time derivatives of system (13) and of (17) respectively, and evaluating at $t=0$. For instance,

[TABLE]

[TABLE]

Then, we can state the following existence result whose proof (in the case of a solid with a flat bottom) is detailed by the author in [3] (Theorem 5.3):

Theorem 2.2**.**

For $k\geq 2$, let $u_{0}=(\zeta_{e,0},q_{e,0})$, $\delta_{0}$ and $\delta_{1}$ satisfy the compatibility conditions in Definition 2.1 up to order $k-1$. Assume that there exist some constants $h_{\mathrm{m}},c_{\mathrm{sub}}>0$ such that

[TABLE]

with $h_{e,0}=h_{0}+\zeta_{e,0},$ and that

[TABLE]

Then, there exists $T>0$ such that the Cauchy problem for (17) with initial data

[TABLE]

admits a unique solution $\delta_{G}\in H^{k+1}((0,T))$.

2.2 The return to equilibrium configuration

We want to focus now on a particular configuration of the floating structure problem, the return to equilibrium problem. It consists in dropping the solid, with no initial velocity, into a fluid initially at rest from a non-equilibrium position. By the definition of this particular configuration, we have specific initial conditions for the coupled problem (10) - (11).

The initial conditions for the solid equation are

[TABLE]

and for the fluid equations are

[TABLE]

for all $r\in(R,+\infty)$. In order to apply the theory of the initial boundary value problem we need these specific initial data to satisfy the compatibility conditions defined in [3]. The compatibility conditions of order 0 and 1 are respectively:

[TABLE]

Due to the nature of the return to equilibrium configuration, we have

[TABLE]

Therefore the compatibility condition of order 0 is satisfied but not the one of order 1. Then Theorem 5.3 of [3] cannot be applied since one hypothesis required is that the initial and boundary data must satisfy the compatibility conditions at least up to order $1$. When the compatibility conditions at order 1 are not satisfied, sonic waves propagate (we refer to Métivier [13] for the existence of such waves).

Remark 2.3**.**

One can choose a different value for $\delta_{1}$ in order to satisfy the compatibility conditions and be able to apply the results of Theorem 5.3 in [3].

3 Linear-nonlinear model for floating structures

The impossibility to apply the theory of initial boundary value problems to the particular configuration of the return to equilibrium brings us to consider a linearization of the equations (7) in the exterior domain, which describes the case of small amplitude waves. We generalize however the works by Cummins and other authors in the literature by keeping the nonlinear effects in the interior domain. We only assume that the solid does not touch the bottom of the fluid domain. In this section we introduce the linear-nonlinear model for the floating structure problem, we prove the conservation of the total energy for this model and then we show that with this linear approximation we can write the extension-trace operator $\mathcal{B}[\dot{\delta}_{G},u_{0}]$ (simply written $\mathcal{B}[\dot{\delta}_{G}]$ from now on) as a linear convolution operator. Then the delay differential equation (17) for the solid motion becomes a nonlinear second order integro-differential equation.

3.1 An energy conserving linear-nonlinear model

We consider the following linear-nonlinear model for the floating structure problem:

• •

in the exterior domain $(R,+\infty)$

[TABLE]

• •

in the interior domain $(0,R)$

[TABLE]

and the boundary conditions

[TABLE]

[TABLE]

with $P_{\mathrm{cor}}=-\dfrac{\rho}{2}{\dfrac{q_{i}^{2}}{h_{i}^{2}}}_{|_{r=R}}$. As in the full nonlinear case the condition (21) can be written in terms of the solid vertical displacement $\delta_{G}$ and it becomes

[TABLE]

Furthermore we have the conservation of the energy for the new linear-nonlinear model (see [3] for the conservation of the energy in the full nonlinear model):

Proposition 3.1**.**

Let us define the shallow water fluid energy for the linear-nonlinear shallow water equations (19) - (20)

[TABLE]

and the solid energy (only with vertical motion)

[TABLE]

Then the total fluid-structure energy $E_{tot}=E_{SW}+E_{sol}$ is conserved, i.e.

[TABLE]

Proof.

Multiplying the first equation of (19) by $\rho g\zeta_{e}r$, the second equation by $\dfrac{\rho q_{e}}{h_{0}}r$ and summing them together, we have local conservation of the energy, i.e.

[TABLE]

where $e_{ext}$ is the local fluid energy in the exterior domain

[TABLE]

and $F_{ext}$ is the flux in the exterior domain

[TABLE]

We consider the equations (20) in the interior domain. Analogously, multiplying the first equation by $\rho g\zeta_{i}r$, the second equation by $\dfrac{\rho q_{i}}{h_{i}}r$ and summing them together, we obtain

[TABLE]

where $e_{int}$ is the local fluid energy in the interior domain

[TABLE]

and $F_{int}$ is the flux in the interior domain

[TABLE]

For the sake of clarity we remark that, in order to derive $\partial_{t}e_{int}$ in (26), we have used the following identity for the convective term:

[TABLE]

where the last equality is due to the first equation in (20) and to the fact that $\partial_{t}h_{i}=\partial_{t}\zeta_{i}$. We integrate (25) on $[R,+\infty)$ and (26) on $[0,R]$ and multiplying them by $2\pi$ we obtain

[TABLE]

where \left\llbracket f\right\rrbracket is the jump of a function $f$ at the boundary $r=R$ defined as

[TABLE]

By integration by parts we get

[TABLE]

On the other hand, from the definition of $E_{sol}$, we have

[TABLE]

where we used Newton’s law for the conservation of the linear momentum and, since the structure moves only vertically,

[TABLE]

coming from standard solid mechanics. From the contact constraint (4) and the mass conservation equation in (20) we get

[TABLE]

Therefore

[TABLE]

Using the expression of the interior pressure $\underline{P}_{i}$ on the boundary $r=R$ in (22) and the transition condition (21) we get the conservation of the total energy. ∎

3.2 Linear equations in the exterior domain

In this subsection we focus on the linear shallow water equations in the exterior domain

[TABLE]

with $v_{0}=\sqrt{gh_{0}}$, coupled with the transition condition

[TABLE]

Taking the derivative of the first equation in (30) with respect to time and replacing the value of $\partial_{t}q_{e}$ with the expression in the second equation we find the linear wave equation

[TABLE]

with $\Delta_{r}:=\partial_{rr}+\dfrac{1}{r}\partial_{r}$.

We consider only positive time $t$ (we can treat $\zeta_{e}$ as a causal function, i.e. $\zeta_{e}=0$ for $t<0$). In the same way as John did in [8], we apply the Laplace transform

[TABLE]

to the wave equation and we get the following Helmholtz equation with complex coefficients:

[TABLE]

We have $\mathcal{L}\left(\partial_{tt}\zeta_{e}\right)=s^{2}\mathcal{L}\left(\zeta_{e}\right)+\partial_{t}\zeta_{e}(0)+s\zeta_{e}(0)$ but in this configuration we have in addition $\partial_{t}\zeta_{e}(0)=0$ and $\zeta_{e}(0)=0$ from (18). The general solution of (32) is

[TABLE]

where $H_{0}^{(1)}$ and $H_{0}^{(2)}$ are the Hankel functions of first order and second order respectively with index 0.

Remark 3.2**.**

Let us consider the Bessel functions of the first kind and of the second kind, respectively $J_{n}$ and $Y_{n}$, solutions to

[TABLE]

The Hankel functions of first order with index $n$ are defined as

[TABLE]

and the Hankel functions of second order with index $n$ as

[TABLE]

From the asymptotic behaviour of the Hankel functions (see Appendix A) we know that $H_{0}^{(1)}(z)\sim\sqrt{\dfrac{2}{\pi z}}e^{iz}$ and $H_{0}^{(2)}(z)\sim\sqrt{\dfrac{2}{\pi z}}e^{-iz}$ for large $|z|$ and $0<\arg z<\pi.$ Therefore for large $|s|r$ and $-\dfrac{\pi}{2}<\arg s<\dfrac{\pi}{2}$

[TABLE]

[TABLE]

Thus for $\Re(s)>0$ and large $r$

[TABLE]

[TABLE]

These terms represent respectively an outgoing progressive wave and an incoming progressive wave. Since in this problem we consider only outgoing waves, we impose $a_{2}(s)=0$.

Applying the Laplace transform to the second equation of (30), we get the following boundary condition for the exterior Helmholtz problem:

[TABLE]

using the transition condition (31). Therefore we finally have

[TABLE]

using the relation $(H_{0}^{(1)})^{\prime}=-H_{1}^{(1)}$ between the derivative of $H_{0}^{(1)}$ and the Hankel function of first order with index 1. From Appendix A we have $\dfrac{H_{0}^{(1)}(s)}{H_{1}^{(1)}(s)}\rightarrow i$ for large $|s|$. Adding and subtracting this limit we have

[TABLE]

with

[TABLE]

with $f(s)\rightarrow 0$ as $|s|\rightarrow+\infty$. It turns out that we can write $f$ as a Laplace transform of some function:

Lemma 3.3**.**

There exists a unique function $F\in L^{2}\left(\mathbb{R}_{+}\right)\cap C\left([0,+\infty)\right)$ such that $f(s)=\mathcal{L}\left(F\right)(s)$, with

[TABLE]

independent of $c>0$, in the sense of $L^{2}$-Fourier transform and

[TABLE]

with $\lambda=\frac{1}{4}$, in the sense of Lebesgue integral.

Proof.

We know that both $H_{0}^{(1)}(is)$, $H_{1}^{(1)}(is)$ are holomorphic functions on $\mathbb{C}_{+}$, and $H_{1}^{(1)}(is)\neq 0$ in $\mathbb{C}_{+}$ (see [1],[5]), then $f(s)$ is holomorphic on $\mathbb{C}_{+}$. Moreover $f$ is bounded in $\mathbb{C}_{+}$ since $f\rightarrow 0$ at infinity and $f$ is bounded around the boundary $i\mathbb{R}$ (from Appendix A we have $\frac{H_{0}^{(1)}(is)}{H_{1}^{(1)}(is)}\sim-is\log(is)$ for $s\rightarrow 0$). Hence $f\in H^{\infty}(\mathbb{C}_{+})$. Now we want to show that $f\in L^{2}(i\mathbb{R})$: $f$ is defined also in $\overline{\mathbb{C}_{+}}$ if we consider the one-valued functions $H_{0}^{(1)}$ and $H_{1}^{(1)}$ (considering the one-valued logarithm in the definition of the Hankel functions in Appendix A). Moreover we have that

[TABLE]

as $|s|\rightarrow+\infty$, hence

[TABLE]

Therefore by the Smirnov theorem (see [14]) $f\in H^{2}(\mathbb{C}_{+})$, where $H^{2}\left(\mathbb{C}_{+}\right)$ is the so-called Hardy space, and by the Paley-Wiener theorem (see [6, 22]) there exists a unique function $F\in L^{2}\left(\mathbb{R}_{+}\right)$ such that $\mathcal{L}\left(F\right)(s)=f(s)$ with

[TABLE]

is to be understood in the sense of $L^{2}$ Fourier transforms for any $c>0$. On the other hand, from (35) we have $g(s)=f(s)-\frac{1}{4s}$ is Lebesgue integrable on the line $\mathrm{Re}s=c$ for any $c>0$. From Lemma 3.9. of [15] there exists a function $\widetilde{F}\in C([0,+\infty))$ such that $\mathcal{L}\left(\widetilde{F}\right)(s)=g(s)$, with

[TABLE]

independent of $c>0.$ Hence, writing $f(s)=g(s)+\frac{1}{4s}$ and using the fact that $\mathcal{L}\left(\lambda\right)=\frac{\lambda}{s}$ for all complex constant $\lambda$, we have that $\mathcal{L}\left(F\right)(s)=f(s)$ with

[TABLE]

and $\lambda=\frac{1}{4}.$ ∎

Then we can write the coupling term with the fluid motion $\zeta_{e}(t,R)$ as an explicit function of the solid velocity $\dot{\delta}_{G}$ under convolution form:

Proposition 3.4**.**

Considering the linearized shallow water equations in the exterior domain, the following holds:

[TABLE]

with $F(t)$ as in Lemma 3.3.

Remark 3.5**.**

It is easy to see that there exists a function $F_{0}$ such that

[TABLE]

with

[TABLE]

Now we show the long time behaviour of $F$, which due to the previous remark is is independent on the choice of the parameters $v_{0}$ and $R$.

Proposition 3.6**.**

$F$* has a long time polynomially decaying behaviour. In particular there exist $M>0$ such that*

[TABLE]

for all $t\geq 0.$

Proof.

Without loss of generality we prove the lemma for $F_{0}$, which is a function independent of the parameters $v_{0}$ and $R$. As said before, $F(t)=F_{0}(\frac{v_{0}}{R}t)$ and the asymptotic behaviour of $F$ is the same as the one of $F_{0}$.

First, let us recall that $F_{0}$ is defined as the inverse Laplace transform of $f_{0}$ that has a branch cut on the real negative semi-axis by choosing the value of the complex logarithm that has a branch cut in the lower imaginary semi-axis in the series expansion of the Bessel function $Y_{n}$ for $n=0,1$ in Appendix A. Let us consider the closed curve (see Figure 2)

[TABLE]

where

[TABLE]

with $\varepsilon,\,\delta_{1},\,\delta_{2},L>0$ and $c>0$ as in Lemma 3.3. We choose $\varepsilon$ small enough such that the zeros of $H_{1}^{(1)}(is)$ lay in the left of the curve $\mathcal{C}$. This can be done since from [5] the zeros of $H_{1}^{(1)}(is)$ all have real part less than $-\alpha$ for some $\alpha>0$. Hence $f_{0}$ is holomorphic in the interior of $\mathcal{C}$ and by Cauchy’s integral theorem

[TABLE]

or equivalently

[TABLE]

Taking the limit $L\rightarrow+\infty$, the limit $\delta_{2}\rightarrow 0^{+}$ and the limit $\delta_{1}\rightarrow 0^{+}$, by definition of the inverse Laplace transform as a Bromwich integral we have

[TABLE]

where we have used the fact that

[TABLE]

by uniform convergence. Let us deal with the first two terms in (37). As we did in Lemma 3.3, we write $f_{0}(s)=g_{0}(s)+\frac{1}{4s}$ with $g_{0}(s)$ integrable on the line $\mathrm{Re}(s)=-\varepsilon$. Hence we have that

[TABLE]

is finite and is exponentially decaying in $t$. We focus now on the contour integration of $\frac{1}{4s}$. We compute that

[TABLE]

Since the integrand is even and it is bounded close to zero, we have that

[TABLE]

The integral in the right-hand side can be computed by considering the contour integration of the associated complex-valued function along the semicircle of radius $L$ centered at the origin and closed on the real axis. Since the function is holomorphic, by Cauchy’s integral theorem

[TABLE]

where $C_{L}$ is the arch of the semicircle. The second term vanishes as $L\rightarrow+\infty$, hence it yields

[TABLE]

For the last three terms in (37) we have that

[TABLE]

by uniform convergence. Using the analytic continuation formulas in Appendix A, we compute that

[TABLE]

and using the series expansion of $J_{n}$ and $Y_{n}$ for $n=0,1$ in Appendix A we get

[TABLE]

Therefore in (38) we have

[TABLE]

which implies the statement of the lemma. ∎

The polynomial decay of the impulse response function $F$ given by Proposition 3.6 is numerically showed in Figure 3 for a particular set of parameters. Moreover, the kernel $F$ satisfies the following equality, which will be used in the proof of the Theorem 47:

Lemma 3.7**.**

The convolution kernel $F$ is such that

[TABLE]

Proof.

By the definition of the Laplace transform and by Lemma 3.3,

[TABLE]

From Appendix A, we have that, as $s\rightarrow 0$,

[TABLE]

Hence, taking the limit $s\rightarrow 0^{+}$ in (41) we get

[TABLE]

where we have used Lebesgue’s dominated convergence theorem due to Lemma 3.3 and Proposition 3.6. ∎

3.3 Integro-differential equation for the solid motion

From now on we suppose for simplicity that the bottom of the structure is flat, then $\zeta_{w}$ (as well as $h_{w}$) does not depend on the space variable $r$, but Proposition 3.8 holds for a structure with non-flat bottom as well. We know from Proposition 3.4 that, considering the linear shallow water equations (19) in the exterior domain, we can write the trace of the surface elevation $\zeta_{e}$ at the boudary $r=R$ as a function of the time derivative of the displacement $\delta_{G}.$ Then the nonlinear differential equation (11) describing the solid motion can be written as a nonlinear delay differential equation.

Proposition 3.8**.**

Considering the linear shallow water equations (30) for the fluid motion in the exterior domain, the solid motion is described by the following second order nonlinear integro-differential equation:

[TABLE]

with $\mathfrak{c}$ as in (12), $\nu=\dfrac{\mathfrak{c}R}{2v_{0}}$, $m_{a}(\delta_{G})=\dfrac{b}{h_{w}(\delta_{G})},$ $\beta(\delta_{G})=\dfrac{b}{2h^{2}_{w}(\delta_{G})}.$ The convolution kernel is defined by

[TABLE]

for any $c>0$ and

[TABLE]

Remark 3.9**.**

In the integro-differential equation (43) $m_{a}(\delta_{G})$ is the time dependent added mass, $\mathfrak{c}$ is the hydrostatic coefficient and $\nu$ is the damping coefficient. The convolution term, whose kernel $F$ is the so-called impulse response (or retardation) function, accounts for fluid-memory effects that incorporate the energy dissipation due to the radiated waves coming from the motion of the structure. Moreover, linearizing (43) around the equilibrium state, we get

[TABLE]

This linear equation is nothing but the well-known Cummins equation for the vertical displacement (1). Proposition 3.8 therefore provides a rigorous justification of the Cummins equation and generalizes it to take into account the nonlinear effects in the interior domain.

Remark 3.10**.**

Recall that in Proposition 3.4 we show that

[TABLE]

Therefore, considering the linear equations (19), the extension-trace operator (16) becomes a linear convolution operator, that is

[TABLE]

with the convolution kernel $K(s)=F(s)-\dfrac{R}{2v_{0}}\mathfrak{d}_{s=0},$ where $\mathfrak{d}_{s=0}$ is the Dirac delta distribution.

We state now the following global existence and uniqueness result of the solution to the solid motion equation in the case of linear shallow water equations for the fluid motion. Let us recall that $\rho_{m}$ and $H$ denote the density and the height of the solid, respectively.

Theorem 3.11**.**

The Cauchy problem for the nonlinear second order integro-differential equation (43) with initial data

[TABLE]

admits a unique solution $\delta_{G}\in C^{2}([0,+\infty),\mathbb{R})$ provided

[TABLE]

Remark 3.12**.**

One needs to consider the parameters of the problem such that

[TABLE]

which means that the fluid height under the solid at the equilibrium position is positive.

Proof.

Let us write (43) as the following first order nonlinear integro-differential equation on $x(t)=(\delta_{G}(t),\dot{\delta}_{G}(t))^{T}$:

[TABLE]

with $K(t,s)=F(t-s).$ The local existence and uniqueness of the solution to (48) can be obtained using an iterative method in the Banach space $C([0,T),\mathbb{R}^{2})$. This integro-differential equation can also be seen as a functional differential equation with infinite delay by extending $\delta_{G}$ and $\dot{\delta}_{G}$ respectively by $\delta_{0}$ and by [math] for $t<0$. We refer to Liu and Magal [10] for the analysis of this type of delayed equations, which would give asymptotic stability of the equilibrium position in the case of an exponentially decaying convolution kernel.

We show here the extension of the existence interval $[0,T)$ to $[0,+\infty)$ using the conservation of total fluid-structure energy. First, we prove in the following lemma that the solution is bounded.

Lemma 3.13**.**

The displacement $\delta_{G}$ and its derivative $\dot{\delta}_{G}$ are both bounded.

Proof.

From Proposition 3.1 we know that the energy of the coupled floating structure system considering the linear shallow water equations for the fluid motion

[TABLE]

is conserved. Moreover, $E_{SW}(t)$ can be written as the sum of the fluid energy in the interior domain,

[TABLE]

and the fluid energy in the exterior domain,

[TABLE]

To get the expression of the fluid energy in the interior domain we use the constraint (4) and we add the constant term $\frac{1}{2}\rho g\,2\pi\int_{0}^{R}\zeta^{2}_{w,eq}rdr$ in order to have zero energy at the equilibrium. From Archimedes’ principle we have

[TABLE]

and, since the bottom of the solid is flat, we have

[TABLE]

Then

[TABLE]

and

[TABLE]

Using (51) and the fact that $q_{i}(t,r)=-\dfrac{r}{2}\dot{\delta}_{G}(t)$ (see Section 2), the fluid energy in the interior domain $E_{int}(t)$ becomes

[TABLE]

In particular the total energy at instant $t=0$ is

[TABLE]

using $\delta_{G}(0)=\delta_{0}$ and $\dot{\delta}_{G}(0)=0.$ By the conservation of total energy and using the identity $m=\rho_{m}\pi R^{2}H$ we have

[TABLE]

Consider $t^{*}=\sup\{t\in[0,T)\ |\ h_{w}(s)>0\mbox{ for}\ s\in(0,t)\}$. From condition (47) we have $h_{w}(0)=h_{w,eq}+\delta_{0}>0$ , hence $t^{*}>0$. Suppose $t^{*}<T$. Then, for $t\in(0,t^{*})$ the right-hand side of (52) has to be non-negative. By solving the inequality with respect to $\delta_{G}(t)$, we have

[TABLE]

By the non-negativity of $E_{ext}(t)$ we get the bound

[TABLE]

Moreover, from (52) we have

[TABLE]

which yields the bound

[TABLE]

Using condition (47), by continuity we have $h_{w}(t^{*})\geq h_{w,eq}-|\delta_{0}|>0$ and there exists $\epsilon>0$ small enough such that $h_{w}(t^{*}+\epsilon)>0$, where $t^{*}$ is the maximal time such that $h_{w}(t)>0$ for $t\in(0,T)$. Then necessarily $t^{*}=T,$ which implies that the bounds (53) - (54) hold in the existence interval $(0,T)$. ∎

The bounds (53) - (54) give

[TABLE]

[TABLE]

We remark that $F\in L^{1}(\mathbb{R}_{+})$ due to Lemma 3.3 and Proposition 3.6. The admissibility condition (47) on $\delta_{0}$ guarantees that for all $t\geq 0$

[TABLE]

[TABLE]

Therefore, the existence interval can be extended to $[0,+\infty)$ by iterating the previous argument. ∎

Remark 3.14**.**

The conditions (55) - (56) express the physical fact that, during all the motion, both the solid and the fluid trace at the solid walls do not touch the bottom of domain.

4 Numerical method

In order to solve numerically the delay differential equation (43) we write it under the form

[TABLE]

with $d_{1}(t),...,d_{k}(t)$ the components of the non-constant delays vector $d(t)$. In our case we have chosen $d_{k}(t)=t-kdt$ with $dt=0.1$ and $k=1,...,N$ for $N=100$. Then, we implement in our code the MATLAB solver ddesd, which integrates with the explicit Runge-Kutta (2,3) pair and interpolant of ode23. For more details on the solver we refer to Shampine [17]. Moreover, we compute the convolution integral applying the trapezoidal integration method following Armesto et al. [2]. In an analogous way, we compute the convolution kernel $F$ for a given set of time steps $n\Delta t$ with $n=1,...,N$ since the influence of the Kernel is negligible after some time $t^{*}=N\Delta t$. Then, we compare the numerical result given by the nonlinear integro-differential equation (43) with the one obtained from its linear approximation. In Figure 4 we consider $h_{0}=15\ \text{m}$, $R=10\ \text{m}$, $H=10\ \text{m}$, $\rho=1000\ \text{kg}/\text{m}^{3}$ and the volume density of the solid $\rho_{m}=0.5\ \rho.$ We choose two different initial data: $\delta_{0}=1\ \text{m}$ and $\delta_{0}=5\ \text{m}$. One can see that for large amplitudes the nonlinear effects should not be neglected in order to better describe the solid motion. This difference justifies the approach to keep nonlinearities in the equation of the floating body problem in the interior domain. Moreover,one can note that the displacement goes to zero but the structure definitely does not reach its equilibrium position: this is due to the motion of the fluid which makes the solid constantly move.

Appendix A Hankel functions

In this appendix we show some results and properties for the Hankel functions. Let us consider the following differential equation:

[TABLE]

This differential equation is called Bessel equation of index $\nu$. Solutions to this equation are called Bessel functions. Let us consider the case when $\nu=n$, with $n\in\mathbb{Z}$. Bessel functions of the first kind, denoted by $J_{n}(z)$,

[TABLE]

are entire in z.

Bessel functions of the second kind, denoted by $Y_{n}(z)$

[TABLE]

where $\psi=\dfrac{\Gamma^{\prime}}{\Gamma}$, with $\Gamma$ the Gamma function, have a branch point in $z=0$. Both $J_{n}$ and $Y_{n}$ are real valued if $z$ is real. Let us define

[TABLE]

[TABLE]

We call them respectively $Hankel\mbox{ }functions$ of first order and second order with index $n$, and they are solutions to the Bessel equation. Each solution has a branch point at $z=0$ for all $n$. The principal branches of $H_{n}^{(1)}(z)$ and $H_{n}^{(2)}(z)$ are two-valued and discontinuous on the cut along the negative real axis. They are holomorphic functions of $z$ throughout the complex plane cut (see Chapter 9 of [1]).

Now let us show some representations of these functions useful for our problem. From [20] we have an integral representation for $z=x>0$:

[TABLE]

and

[TABLE]

and a series representation for large $|z|$ and $0<\arg z<\pi$:

[TABLE]

[TABLE]

with

[TABLE]

[TABLE]

Last we recall analytic continuation formulas for $m\in\mathbb{Z}$ (see [5]):

[TABLE]

[TABLE]

[TABLE]