Modelling of Nonlinear Wave-Buoy Dynamics Using Constrained Variational Methods

We consider a comprehensive mathematical and numerical strategy to couple water-wave motion with rigid ship dynamics using variational principles. We present a methodology that applies to three-dimensional potential flow water waves and ship dynamics. For simplicity, in this paper we demonstrate the method for shallow-water waves coupled to buoy motion in two dimensions, the latter being the symmetric motion of a crosssection of a ship. The novelty in the presented model is that it employs a Lagrange multiplier to impose a physical restriction on the water height under the buoy in the form of an inequality constraint. A system of evolution equations can be obtained from the model and consists of the classical shallow-water equations for shallow, incompressible and irrotational waves, and relevant equations for the dynamics of the wave-energy buoy. One of the advantages of the variational approach followed is that, when combined with symplectic integrators, it eliminates any numerical damping and preserves the discrete energy; this is confirmed in our numerical results.


INTRODUCTION
Fast ships deliver all kinds of services, for example transportation of wind farm maintenance personnel and provision of supplies for oil and gas platforms. Ships are especially vulnerable to extreme wave phenomena such as waves with anomalously high amplitudes relative to the ambient waves. The aim of this paper is to make progress towards investigating the effect of nonlinear water waves, as well as amplified rogue waves, on floating structures such as ships and wave-energy devices. With that ultimate goal in mind, we consider a reduced two-dimensional problem with waves impacting a floating wave-energy buoy with a simple linear structure, constrained to move upward only (as shown in Fig. 1). The simplification is performed in order to understand the mechanics behind the wave-buoy coupling and to test the accuracy of the method. Moreover, this allows our novel approach (described later) to be tested on a simpler problem, before generalising it to include three-dimensional ship dynamics coupled to potential flow water waves.
We present a mathematical and numerical model of the motion of a buoy in (non)linear waves. We derive a coupled model for the wave-buoy dynamics by following a variational methodology [1,2,3], in order to ensure zero numerical damping which is important for wave propagation. The final system of evolution equations comprises the classical shallow-water-wave equations for incompressible and irrotational waves, and a set of equations describing the dynamics of the wave-energy buoy. The novelty in our model is in the presence of a physical restriction on the water height under the buoy, which is enforced through an inequality constraint and by employing a Lagrange multiplier.
Our approach with the use of an inequality constraint avoids the need of splitting the domain in two subdomains: one with a free surface h(x,t) and a second subdomain where the free-water surface is constrained by the presence of the buoy and hence h(x,t) = h b (x, Z(t)) in that region. The latter condition physically means that the water under the buoy takes exactly the shape of the wetted part of the buoy. In Kalogirou and Bokhove [4], we imposed this constraint by using a Lagrange multiplier λ in the second subdomain only. The disadvantage of the approach presented in [4] is that the finite element discretisation leads to the emergence of local matrices valid in the region under the buoy only.
The use of inequality constraints is typical in molecular and polymer dynamics [5] and optimal control systems [6], but it was never seen used in the context of water waves. In our problem, the inequality constraint is imposed weakly in the variational principle and also includes the equality constraint presented in Kalogirou and Bokhove [4] as a sub-case. Alternative approaches of enforcing similar inequality constraints can be found in [7].
The model is solved numerically using a variational (dis)continuous Galerkin finite element method with a special, new and robust time integration method. The numerical imple-mentation employs second-order continuous Lagrange polynomial approximations in space and a constrained extension of the second-order (dis)continuous Störmer-Verlet symplectic scheme in time (i.e. the RATTLE algorithm). The model is implemented in Firedrake [8], an automated system for the solution of partial differential equations with the finite element method. We have shown the numerical results of the linearised coupled wave-buoy system before [4,9] and here we will confirm their validity using the new, novel approach described above. The numerical results confirm conservation of total mass and energy, as well as phasespace volume.

MATHEMATICAL MODEL Variational Principle For Shallow-Water Waves
A classical mathematical theory in the modelling of incompressible and irrotational water waves is the so-called potential flow theory. It's main advantage is that it reduces the degrees of freedom in the problem by describing the fluid velocity u u u in terms of the gradient of a scalar field, the velocity potential Φ. The dynamics of water waves are also Hamiltonian, namely they have the property of conserving the total (i.e. kinetic and potential) energy of the system. In this case, the dynamics can be described using variational methods, in particular the well-known Luke's variational principle for water waves [1].
Consider a two-dimensional domain with horizontal coordi- is the local wave height that evolves in time t ∈ [0, T ]. In the shallow-water approximation considered here, the dynamics are reduced to the dynamics on the free surface and are fully described by the total wave height h(x,t) and free-surface velocity potential φ (x,t) := Φ(x, h(x,t),t). The simplification of Luke's variational principle [1] for shallow-water waves is which essentially provides an expression for the fluid's pressure (as defined by Bernoulli's principle). The above variational principle states that the variations of the pressure are zero, which is equivalent to saying that the pressure remains constant throughout the fluid.

Linear Wave-Energy Buoy
Suppose we have a two-dimensional channel filled with water up to a uniform height H 0 , and a wave-energy buoy of mass M residing at the right end of the channel. The buoy's motion is restricted as it can only move in the vertical; the position of the buoy can be found by following the vertical displacement of the centre of mass Z(t) as a function of time. The buoy's vertical velocity is denoted by W (t). The buoy resides at z = h b (x, Z(t)) above the bottom in a certain horizontal interval from x p (t) ≤ x ≤ L, with x p the edge of the buoy at the water surface. As a clear example, a linear buoy has the form where Z − H is the position of the keel (at x = L), α the angle between the hull and the horizontal and L the length of the domain. A simple schematic diagram of the problem configuration is given in Fig. 1. We artificially extend this buoy function h b smoothly into the entire domain so that it is defined for all 0 ≤ x ≤ L. The waterline point at x p (t) is defined (implicitly) by the equality h(x p (t),t) = h b (x p (t), Z(t)), which essentially imposes that the height of the free surface should equal the buoy's surface height at that point x = x p (t).

Coupled Model For Nonlinear Wave-Buoy Dynamics
A coupled model describing the evolution of water waves and the motion of a buoy can be obtained by combining the Lagrangian densities of the two systems into one variational principle. We therefore extend Luke's variational principle (1) for shallow-water waves and add the additional degrees of freedom for the buoy, i.e. we write The above variational principle is still incomplete, as it doesn't include any condition on the water height under the buoy. In the region x ∈ [x p , L], the water height h(x,t) is not an independent variable but is restricted by the presence of the buoy. We impose this physical constraint in the form of an inequality constraint as seen next. The constrained variational principle is where we have first integrated by parts and changed the sign of the variational principle (3), and then imposed the constraint with global Lagrange multiplier λ . The use of the slack function µ 2 in the constraint is such that to enforce the inequality condition (h − h b ) ≤ 0 as will be seen later.
Taking the variations in (4) results in To go from step (6a) to (6b), we have integrated the second and last lines by parts and applied the end-point conditions δ φ | T 0 = 0 and δW | T Steady-state configuration with buoy and water being stationary. The water has a uniform depth H 0 and the waterline point is fixed at L p .
The resulting equations of motion are the following When looking at Eqs. (7c)-(7d) obtained by the variations of δ λ and δ µ, respectively, it becomes obvious that in different regions of the domain we can have one of the two following cases: (1) In the part of the domain where (2) In the part of the domain under the buoy where h = h b , i.e. for x p ≤ x ≤ L, we get from (7c)-(7d) that The above conditions (8) are called complementary slackness conditions [6,7]. Therefore, by introducing the constraint (5), we have imposed the non-negative nature of (h b − h) both as inequality as well as an equality.

Steady State
When the free-water surface is flat and the buoy is stationary, the problem admits a steady-state solution, given by as shown in Fig. 2. Here, the function H(x) denotes the water level at rest and is defined as the buoy's hull shape at the floating rest state. The rest water level is essentially equal to the flat state H 0 in the region where only water exists in the channel, whereas in the region with the buoy it is given by the linear shape of the buoy's hull H b (x).
The steady-state solution also satisfies the equations of motion (7), which become From the above system, we find On applying Archimedes' principle for the volume of displaced water (which should equal the volume of the buoy), we can find an exact expression for the rest position of the buoy's centre of mass, given bȳ Finally, Eqs. (11)-(13) also define the waterline point at rest x = L p , which is now independent of time and found to be

Linearised Dynamics
The nonlinear variational principle (4) can be linearised around the rest state by considering the steady-state solution (9)-(13) plus small perturbations In what follows, all perturbation variables are written without the tildes for simplicity.
In the linearised variational principle, the constant terms are discarded (as their variations are zero), terms linear in the perturbation variables cancel after using the steady-state solution, and only quadratic terms are kept. Hence, we find the following variational principle for the linearised equations of motion after using the end-point conditions δ φ | T 0 = 0 and δW | T 0 = 0. The above variations lead to the linear equations of motion Once again, we get the following two cases in different parts of the domain: (1) In region 0 ≤ x ≤ L p , we have thatλ = 0 andμ > 0 (see (12a)); therefore it follows from (17c)-(17d) that in the interval with a free surface λ = 0 and µ = 0.

NUMERICAL DISCRETISATION OF THE EQUATIONS Spatial Discretisation: Finite Element Method
We split the horizontal domain [0, L] in N finite elements 0 = x 0 < x 1 < · · · < x k < · · · < x N = L of equal sizes, and consider continuous piecewise-linear basis functions ϕ k (x) in each of the elements. The space-dependent variables η, φ , λ and µ are approximated by η h , φ h , λ h and µ h and the latter are expanded as follows where the Einstein summation convention is used for notational brevity. The above finite element expansions are substituted in the linearised variational principle in Eq. (16a) and the following matrices are introduced with indices k, used to denote the total number of nodes. The time-continuous but space-discrete linear variational principle is therefore given by

Time Discretisation: Symplectic Methods
After obtaining the variational principle (21) in a spacediscrete form, we then apply a modified RATTLE algorithm to discretise the variational principle in time. RATTLE [10] was developed alongside with SHAKE to provide a stable timediscretisation symplectic scheme for constrained Hamiltonian systems [11,12], and it can be seen as an extension of the Störmer-Verlet method. The final space-time discrete equations of motion for the linear system are the following Since the above system is linear, we can easily solve the two Eqs. (22e)-(22f) to determine solutions λ n+1/2 k and µ n+1/2 k first. In order to obtain these solutions, the unknowns η n+1 k and Z n+1 have to be eliminated from the discrete constraint Eq. (22e). This can be achieved by first using Eqs. (22c) and (22d), respectively, followed by use of Eqs. (22a) and (22b) to eliminate φ n+1/2 and W n+1/2 . Equations (22e)-(22f) therefore reduce to a linear system for λ n+1/2 k and µ n+1/2 k , given by Equivalently, we can defineλ ∆t B k , in which case the above system can be written in the following matrix form The above left-hand-side matrix is symmetric. Upon solving linear system (24) and knowing the solution forλ n+1/2 , the rest of the equations can be solved in the following order

Firedrake Implementation of semi-linear system
The linear system (24) and the final system of Eqs. (25) can be solved easily using linear algebra in a mathematical software such as MATLAB. However, the system becomes more complicated when the nonlinear problem is considered, in which case an iterative method needs to be used for the solution of a nonlinear system of equations (such as Newton's and Halley's methods). The numerical package Firedrake [8] is a user-friendly system which makes the implementation of nonlinear systems much simpler, by employing linear and nonlinear solvers provided by PETSc [13,14,15].
In this section, we present a semi-linear system that is obtained by considering linear waves and a nonlinear buoy (and therefore nonlinear motion of the waterline point x p ). Firedrake requires a time-discrete system to be provided by the user, while the finite element discretisation in space is performed on a level lower than the user interface. The user also specifies the following: 1. mesh: one-dimensional domain of length L, split into N elements; 2. function space: Lagrange polynomials, here chosen to be linear and continuous across the elements; 3. functions: numerical approximation of the unknowns η h , φ h , λ h , µ h ; 4. test functions: δ η h , δ φ h , δ λ h , δ µ h ; 5. constants: numerical approximation of the scalars Z h , W h ; 6. weak formulations: system of equations discretised in time.
At every time step ∆t, the current solutions will be denoted with superscript n, and the updated solutions at the next time step will have a superscript n + 1.
In our case we also introduce a new variablẽ and a new test function v, to overcome the difficulty of a term involving double integrals (represented by Q k Q λ n+1/2 k in the space-discrete case above -see Eq. (23a)). The variableĨ belongs in a function space consisting of functions on the real line (the function space R), which is typically employed to model global constraints. The presence of the unknown functionĨ from space R in our system (given next in (27)), leads to implementation difficulties due to the coupling with all of the other degrees of freedom.
The weak formulations implemented in Firedrake are the following ,Ĩ, forming a mixed system in Firedrake. This system of weak formulations is nontrivial as it involves a function over R which is constant over the domain. The remaining Eqs. (27f)-(27i) are then explicit and can be solved in the order provided.
Preconditioning. In order to be able to use iterative methods effectively, a suitable preconditioner P must be defined [16]. To find a suitable preconditioner for the mixed system (27a)-(27e), we look at the corresponding linear problem which allows us to define a left-hand-side matrix/operator A which has a natural block structure. The resulting preconditioned matrix M = P −1 A is then spectrally equivalent to A and has a smaller spectral condition number [17]. The goal is to reduce the matrix A into a block diagonal matrix with some operator for φ n+1/2 h and only unit matrices for the remaining variables η n+1 h ,λ n+1/2 h , µ n+1/2 h ,Ĩ. We therefore find the following (block diagonal) preconditioning matrix with I the identity matrix.

NUMERICAL RESULTS
We consider waves generated by a piston wavemaker located at the left wall of the channel, i.e. at x = R(t). Assuming that the wavemaker is linearised around x = 0, we take x = R(t) ≈ 0. The motion of the wavemaker is given by R(t) = A ω (1 − cos(ωt)), such that R(t) and dR/dt = A sin(ωt) are both zero initially (at t = 0). To accommodate the motion of the wavemaker, we add the term ρH 0 (dR/dt)φ x=0 in the continuum variational principle (4) or (16a). The corresponding term in the space-discrete variational principle (21) is ρH 0 (dR/dt)W k , with W k = ϕ k | x=0 . Note that due to the presence of the wavemaker at the left boundary, the W k vector has only one non-zero entry (the W 1 ) based on the contribution of the ϕ 1 basis (tent) function. In the final spacetime discrete system, the time-dependent term dR/dt is evaluated at time level t n + 1 2 = (n + 1 2 )∆t, n = 0, . . . , T . The initial value problem is solved with initial conditions φ = 0 and η = 0 (see Fig. 4, first panel). Waves are then generated from the left due to the sinusoidal motion of the wavemaker (see Fig. 4, second panel) with amplitude A = 0.0498 m and frequency ω = 8π L √ gH 0 (corresponding to a physical frequency of approximately 3.9618 Hz). The remaining of the parameters that appear in the system take the following values:  (14) and (13), respectively.
In what follows, we present numerical solutions of the linearised problem (24)-(25). The numerical calculation is performed on a uniform mesh with N = 100 elements and the time step is chosen such that to satisfy the CFL condition ∆t = L/(2πN √ gH 0 ) = 0.0016. The response of the buoy is depicted in Fig. 3, where the time-evolution of the vertical displacement Z(t) and velocity W (t) are shown. The numerically computed wave height h(x, y) and the vertical position of the buoy are demonstrated in Fig. 4 for times t = 0.0, 0.4, 0.8, 1.21, 1.65 seconds, while the corresponding velocity potential φ (x,t) at the same times is plotted in Fig. 5. The waterline is marked with black dots and is calculated using a second-order expansion of h(L p +x p ,t) = h b (L p +x p ,Z +Z). The results are identical to the ones obtained by the method presented in Kalogirou and Bokhove [4] with equality constraint.
We have also confirmed that the total energy (or Hamiltonian) of the full system is "constant" after the wavemaker is switched off, in the sense that it remains bounded and exhibits small-but-finite oscillations. This is illustrated in Fig. 6, where the numerically computed energy is shown in the time interval T w ≤ t ≤ T (for wavemaker switch-off time T w = 2 seconds and final time T = 3 seconds). The top-row panels depict the conservation of energy (top right panel) and the error from the "initial" energy E(T w ) (top left panel). The bottom-row panels clearly exhibit the anticipated exchange of energy between water waves (bottom left panel) and the buoy (bottom right panel).

DISCUSSION
A mathematical model for the motion of a wave-energy buoy in nonlinear water waves is presented. The model is developed using variational methods and is based on extending a shallowwater version of Luke's variational principle for water waves to include a floating buoy, allowed to move only vertically. The final system of evolution equations is nonlinear and couples the dynamics of the waves to the dynamics of the buoy through a hydrostatic pressure term. The novelty in our model lies in the use of an inequality constraint to impose the physical condition on the water height being equal to the shape of the buoy in the region under the buoy. The inequality constraint is imposed by the use of a Lagrange multiplier and includes the equality case presented before in Kalogirou and Bokhove [4].
The constraint model is discretised in space using the finite element method (with linear polynomial approximation in each element, and continuous across the elements) and in time using a variation of the RATTLE symplectic algorithm. The numerical implementation is performed using the automated system Firedrake. Numerical solutions are presented for the problem linearised around a steady state, in a two-dimensional channel and linear waves generated by the periodic motion of a piston wavemaker. The total energy of the system is conserved after the wavemaker is switched off, while the energies of water wave and buoy are exactly exchanged between each other.
The three-dimensional analogue of the problem considered here consists of a wave channel that allows the generation of nonlinear rogue-wave effects through a V-shaped contraction, and the impact of such waves on a floating tetrahedral wave-energy buoy constrained to move only vertically in the contraction. We aim to perform further numerical simulations for the dynamics of the coupled system in a hierarchy of increasing complexity: linear water-wave and nonlinear buoy dynamics, and fully-coupled nonlinear water-wave and buoy dynamics. As a validation of the model, laboratory experiments will be performed in a small-scale wave tank with a wave-energy buoy attached to an AC-induction motor, for which we have realised a working-scale model.
Finally, the methodology presented in this paper is generic and can there be readily extended to incorporate a ship in a straightforward manner. The additional complexity arises in the form of additional degrees of freedom: the position of centre of mass in the three-dimensional space X X X = (X,Y, Z), and the rotation of the ship around the three principal axis of rotation θ θ θ = (θ , ϕ, ψ) [18]. The variational principle (3) or (4) can be therefore easily extended to include the dynamics of the ship (such an extended principle can be found in [4]). The extension to three-dimensional potential flow or Benney-Like equations [19] is also possible and we can build on existing applications thereof in Firedrake.