Lagrangian particle path formulation of multilayer shallow-water flows dynamically coupled to vessel motion

The coupled motion—between multiple inviscid, incompressible, immiscible fluid layers in a rectangular vessel with a rigid lid and the vessel dynamics—is considered. The fluid layers are assumed to be thin and the shallow-water assumption is applied. The governing form of the Lagrangian functional in the Lagrangian particle path (LPP) framework is derived for an arbitrary number of layers, while the corresponding Hamiltonian is explicitly derived in the case of two- and three-layer fluids. The Hamiltonian formulation has nice properties for numerical simulations, and a fast, effective and symplectic numerical scheme is presented in the two- and three-layer cases, based upon the implicit-midpoint rule. Results of the simulations are compared with linear solutions and with the existing results of Alemi Ardakani et al. (J Fluid Struct 59:432–460, 2015) which were obtained using a finite volume approach in the Eulerian representation. The latter results are extended to non-Boussinesq regimes. The advantages and limitations of the LPP formulation and variational discretization are highlighted.


Introduction
The Lagrangian and Eulerian descriptions of fluid motion are two viewpoints for representing fluid motion, with the Eulerian description being the most widely used in theoretical fluid dynamics. However, there are some settings where the Lagrangian particle path (LPP) description has advantages, one of which is shallow-water hydrodynamics. In the Eulerian form, the classical non-conservative shallow-water equations (SWEs) are h t + uh x + hu x = 0, u t + uu x + gh x = 0, (1.1) where (a, τ ) are the label and time coordinate in the Lagrangian frame, fluid positions are represented by x(a, τ ) and h = h(a, τ ). The first equation in (1.2) can be integrated to hx a = χ(a) where χ(a) is determined by the initial data, and substitution into the second equation gives Hence, the pair of equations (1.1) has been reduced to a single equation for x(a, t). Moreover the Eq. (1.3) is the Euler-Lagrange equation deduced, with fixed endpoint variations, from the Lagrangian functional x a χ(a)da dτ, (1.4) where for definiteness 0 < a < L. The advantage of the transformation from (1.1) to (1.3) is that variational numerical schemes can be developed, by directly discretizing (1.4), which have excellent energy conservation properties. This energy conservation property is particularly important when the fluid motion is inside a vessel, and it is coupled to the vessel motion, as then it is of interest to accurately capture the energy partition between vessel and fluid. This strategy was used for simulating the dynamic coupling with a single-layer fluid in [1,2]. The aim of this paper is to derive the LPP formulation to shallow water flows, with multiple layers of differing density, in a vessel with dynamic coupling, and use it as a basis for a variational formulation and numerical scheme. Although this generalisation is straightforward in principle, the variational formulation has complex subtleties due to the integration over different label spaces. Stewart and Dellar [3] successfully developed a variational formulation for shallow-water multilayer hydrodynamics by starting with a variational formulation for the full three-dimensional problem and reducing. The resulting variational principle for shallow water involves integration over each layer with respect to the local labels. With an aim to discretize the variational formulation, we modify the Stewart-Dellar formulation by introducing an explicit mapping between label spaces. Then all the integrations are over a single reference label space. Another addition to the variational formulation is that the multilayer shallow-water flow is dynamically coupled to the vessel motion. The theory will be developed in detail first for two-layers in Sect. 2 and then generalised to an arbitrary but finite number of layers in Sect. 4.
A schematic of the problem of interest in the case of two layers is shown in Fig. 1. In Fig. 1, the fluid is coupled to a vessel undergoing horizontal motion only, as there is no vertical acceleration component.
This system is a model for the Offshore Wave Energy Ltd (OWEL) ocean wave energy converter [4]. The OWEL wave energy converter (WEC) is essentially a rectangular box, open at one end to allow waves to enter and, once they have entered the device, the interior sloshing causes the wave to grow pushing air through a turbine and generating electricity. This interior system is a two-layer flow of air and water confined between upper and lower surfaces. This paper considers a simplified model of the OWEL configuration by including two-layers of differing density, but in a closed vessel. In Fig. 1, the vessel displacement q(t) could be prescribed, i.e. the vessel is subjected to a prescribed horizontal time-dependent force, or it could be determined as part of the solution. In this paper, we Here the vessel is constrained to move in the horizontal direction with displacement q(t). Here the displacement is determined by attaching the vessel to a nonlinear spring consider the vessel to be attached to a nonlinear spring, and hence, the vessel motion is governed by a combination of the restoring force of the spring and the hydrodynamic force of the fluid on the side walls of the vessel. The moving vessel walls in turn create a force on the fluid causing it to move, thus the system undergoes complex coupled motions.
The literature on two-layer flows in open systems, with and without a rigid lid is vast ( [5][6][7][8] to name a few), but in closed sloshing systems the literature is much more limited. The theoretical and experimental works of [9,10] show excellent agreement for sloshing in a fixed rectangular tank with a rigid lid when a Lagrangian representation of the system is reduced to a system of ordinary differential equations with dissipative damping. Also, [11] examine two-layer sloshing in a forced vessel and derive a forced KdV equation when the layer thicknesses are comparable in size, an analysis of which shows that forcing induces chaos into the system through homoclinic tangles. The studies most relevant to the one in this paper examine the two-layer sloshing problem using a numerical scheme based upon a class of high resolution wave-propagating finite volume methods known as f-wave methods for both the forced [12] and the coupled problem [13]. This f-wave approach is very effective and can be readily be extended to multilayer systems [14] and systems with bottom topography [15], but [13] find the scheme is limited to layer density ratios of ρ 2 /ρ 1 0.7, where ρ 1 and ρ 2 are the fluid densities in the lower and upper layers, respectively, due to a linear growth in the system constraint error. Therefore this approach is not able to model the interior workings of the OWEL WEC, where the air/water density ratio is ρ 2 /ρ 1 ≈ 10 −3 . The current paper formulates a simple numerical approach based upon a discretization of the LPP scheme, generalizing the numerics of [1] to two layers with nonlinear vessel motion. The LPP approach allows two-layer simulations with ρ 2 /ρ 1 = 10 −3 to be produced.
The principal difficulty introduced by the rigid lid in the multilayer formulation is the Eulerian constraint where h i > 0 and x i are the thickness of and fluid position in the ith layer, respectively, and the sum is over all the layers. In the LPP description, it is necessary to synchronise the position of the Lagrangian particles, otherwise the Eulerian constraint (1.5) will no longer hold at every spatial position. Here we overcome this problem by introducing layer mappings φ i (a, τ ) : [0, L] → [0, L] such that the fluid position functions in layer i satisfy where a 1 is the Lagrangian label in layer 1, and τ is the Lagrangian time variable. The maps φ i (a 1 , τ ) are defined by these constraints. The maps φ i become part of the variational formulation, and the integrals in the Lagrangian functionals are over the single reference space with labels a 1 .
The paper is laid out as follows. In Sect. 2, the governing equations and variational principles for the two-layer rigid lid sloshing problem in the LPP description are presented. In Sect. 3 a variational discretization leading to a symplectic numerical integrator is introduced and simulations are presented. The results include validation of the scheme and extension of the numerical results into the non-Boussinesq regime. In Sect. 4, we demonstrate how the theory is extended to multilayer flows with a rigid upper lid and present simulations for the three-layer problem. Full details for the derivation of the governing three-layer equations is given in an appendix. Concluding remarks and discussion are presented in Sect. 5.

Coupled equations with a two-layer fluid
In this section, we develop the equations for two-layer sloshing in a vessel with rectangular cross-section with a rigid lid coupled to horizontal vessel motion. A schematic of the problem is shown in Fig. 1.
The special case of two-layer flow is of interest for two reasons: Firstly, to simplify the analysis and make the derivation of the governing equations and solution procedure tractable and readable, and secondly because the underlying motivation for this work comes from the two-layer air-water flow inside the OWEL WEC. In Sect. 4 we document how the method presented in this section can be extended to incorporate multilayer shallow-water flow, and present simulation results for the case of three layers.
The rectangular vessel is a rigid body of length L and height d and we consider it filled with two immiscible, inviscid fluids of constant density ρ 1 and ρ 2 with ρ 1 > ρ 2 . The problem is assumed to be two-dimensional with the effect of the front and back faces of the vessel neglected. In what follows, the subscripts 1 and 2 denote the lowerand upper-layer variables respectively. There are two frames of reference in this problem, the inertial frame with coordinates X = (X, Y ) and the body frame with coordinates x i = (x i , y i ) in each layer i = 1, 2. These coordinate systems are related via the time-dependent uniform translation q(t) in the x 1 −direction, and in particular In each layer, the thickness of the fluid h i (x i , t) is assumed to be small such that the layer can be modelled using the shallow-water equations with a corresponding shallow-water velocity field u i (x i , t). The rigid upper lid constrains the flow such that (2.1) As we are considering a vessel with vertical side walls, we could consider the case where x 1 = x 2 and thus only consider one spatial variable, but we leave our notation general for now, to highlight the interesting subtleties of the problem that arise when considering the Lagrangian form. The Eulerian form of the shallow-water mass and momentum equations in each layer in the body frame are where g > 0 is the gravitational acceleration constant, p s (x 2 , t) is the unknown pressure at the rigid lid, and the over dots denote a full derivative with respect to t [13,16]. The fluid in each layer must satisfy the no-penetration conditions on the vessel side walls, and hence, the boundary conditions are The time-dependent motion of the vessel is not known a priori and is determined by a combination of a restoring force, such as a spring or a pendulum [17] and a hydrodynamic force exerted on the vessel side walls by the sloshing fluid. We assume that the vessel is connected to the spatial origin by a nonlinear spring, and hence, the vessel motion is governed by where ν 1 and ν 2 are constant spring coefficients and m (i) is the fluid mass in the ith layer. Here the integrals on the LHS of (2.6) denote the hydrodynamic force contribution of each layer to the vessel motion.
Equations (2.1)-(2.6) can be derived from an Eulerian variational principle by considering variations to the Lagrangian functional where Here p s (x 2 , t) enters as a Lagrange multiplier for the constraint, and f 1 (x 1 , t) and f 2 (x 2 , t) are additional Lagrange multipliers for the mass conservation equations (2.2) and (2.4). The Lagrangian in (2.7) is comprised of three Lagrangian functionals, one for the dry vessel and one for each fluid layer, as discussed in Sect. 4.6 of [3], where the term −ρ 2 gh 1 h 2 in L 1 is identified as the work done on the upper surface of the lower layer by the layer above, and the terms proportional to (u i +q) 2 couple the fluid motion to the vessel motion. One feature of the Lagrangian (2.7) is that the additional work term in L 1 , −ρ 2 gh 1 h 2 , is a function of x 1 , x 2 and t, but the integral, as written above, is over x 1 , moreover, as discussed earlier, the Eulerian constraint h 1 (x 1 , t) + h 2 (x 2 , t) = d has to hold for Both of these issues are overcome in Sect. 2.1 by introducing the constraint that x 1 = x 2 into (2.7) and formulating the problem in terms of the lower layer coordinate only. The shallow-water equations (2.2)-(2.5) could be solved numerically via some implicit shallow-water numerical scheme, with the vessel equation (2.6) integrated via standard fourth-order Runge-Kutta integration. However, this scheme would not necessarily have good energy conservation properties. Hence, in order to model the longtime oscillatory behaviour of the system, we construct a Hamiltonian formulation of the system in order to utilise geometric integration schemes. We do this by transforming the above Eulerian variational formulation to an LPP Lagrangian variational formulation.

Lagrangian variational formulation
To transform the Eulerian shallow-water equations into a LPP formulation, we need to consider mappings from the Lagrangian particle labels and Lagrangian time (a i , τ ) in each layer to the corresponding Cartesian coordinates and Eulerian time (x i , t). This again demonstrates another peculiar feature of the problem, because there is no guarantee that for all τ , x 1 (a 1 , τ ) = x 2 (a 2 , τ ) which we require so as to satisfy the Eulerian constraint (2.1). The approach to overcome this problem is to link the two LPP labels in each layer via a 2 = φ(a 1 , τ ) where φ(a 1 , τ ) is an unknown map to be determined. In the subsequent analysis, we shall drop the subscript 1 from the Lagrangian label a 1 with the understanding that this is the label in the lower layer, and our primary reference label.
To derive the LPP formulation of the problem, consider the mapping with the constraint that in the upper layer We assume that the mapping is non-degenerate (∂ x 1 /∂a = 0) and thus the derivatives in (2.2) and (2.3) map to where here the subscripts a and τ denote partial derivatives. Because we have assumed the constraint (2.9) the derivatives in (2.4) and (2.5) map in the same way as in (2.10), but we can show this formally. From the form of x 2 in (2.9), the derivatives in the LPP setting map on to But we note from (2.9) that and thus the derivatives in the upper layer also map as in (2.10) as noted above. Under this LPP transformation, the fluid equations in each layer, (2.2)-(2.5), transform to where the constraint (2.1) has been used to remove h 2 = d − h 1 from (2.12). Equation (2.11) implies that while adding (2.11) to (2.13) and using (2.1) lead to the mass flux condition after integrating and using the side wall boundary conditions to fix the time-dependent integration function. Eliminating p s between (2.12) and (2.14) and using (2.16) to eliminate x 2τ , (2.1) to eliminate h 2 and (2.15) to eliminate h 1 lead to a PDE in x 1 (a, τ ) and q(τ ) only, which is coupled to vessel equation (2.6) which in the LPP description is The pair of equations (2.17) and (2.18) can also be determined by a variational approach from the Lagrangian (2.7) converted into the LPP description. We directly impose the constraints, and use (2.15) and (2.16) to write the Lagrangian solely in terms of x 1 and q, which takes the form Taking variations, with fixed endpoints, with respect to x 1 and q (e.g. writing q = q +δq with δq(τ 1 ) = δq(τ 2 ) = 0) leads to (2.17) and (2.18) respectively. Note that in the case ρ 2 = 0 (with ν 2 = 0), (2.19) reduces to the one-layer coupled Lagrangian given in [1], i.e. in this case the fluid does not feel the effect of the rigid lid.

Hamiltonian formulation
The coupled Lagrangian system (2.19) can also be written in terms of a Hamiltonian functional with canonical variables (x 1 , q, w, p). The momentum variables are which can be written in the more convenient form The Hamiltonian can then be formed by taking the Legendre transformation of the Lagrangian (2.19). The Hamiltonian functional is given by where I = L 0 Bw da, and the governing form of Hamilton's equations are Here, as in [1], the gradient of H is taken with respect to the weighted inner product such that The form of (2.23) is equivalent to that in (2.17), which was derived directly from the Eulerian form of the equations. This equivalence is shown in Appendix 1.

Linear solutions to LPP problem
The linear solution of the two-layer shallow-water sloshing problem with a rigid lid in the Eulerian framework has been discussed in detail in [13]. However, the form of this linear solution in the LPP framework would be desirable in order to use it as an initial condition when numerically integrating Hamilton's equations so to validate the scheme. Hence, we briefly outline the linear solution procedure here.
We seek a linear solution to (2.17) about a quiescent fluid with the lower layer of mean thickness h 0 1 of the form is the mean thickness of the fluid in the upper layer. Seeking a harmonic solution of these equations with frequency ω emits the separable variable forms The general solution to (2.29) satisfying where γ is an as yet undetermined constant, and when we satisfy x 1 (L , τ ) = L ( X 1 (L) = 0) we find the relation on γ and q that The linear form of the vessel equation (2.18) upon substitution of (2.28), leads to f . Hence, using the above form of X 1 , the vessel equation leads to a second equation linking γ and q Solving (2.31) and (2.32) for non-trivial solutions leads to a characteristic equation for the frequency ω of the form Once the value of s is found from (2.33) then ω is given by A full discussion of the properties of this characteristic equation can be found in [13]. Of interest to us in this paper are the linear forms which we will use to check the validity of the numerical scheme. We are interested in results away from the resonance case (i.e. D(s) = 0, with P(s) = 0). In this case, sin(s) = 0 in (2.31), and hence, γ = β q α 2 tan 1 2 αL .

Numerical algorithm
To formulate the numerical scheme, we discretize the Lagrangian state space into N parcels by setting τ ) (note the dropping of the subscript '1' on the x here) and w i (t) := w(a i , τ ). The derivatives are discretized using forward differences, except when i = N + 1 where backward differences are used, and the integrals are approximated using the trapezoidal rule. Equations (2.24)-(2.26) can be discretized in a straightforward manner, as the variables for which variations are taken, do not appear differentiated with respect to a on the RHS of the equations. However, in (2.23), derivatives of x 1 with respect to a do appear in the RHS, and thus, it is not clear how to discretize these equations. To overcome this, we use a semi-discretization of the Hamiltonian, where the Hamiltonian is discretized and then variations with respect to x i are taken.
To form the semi-discretization, we note that the discretized form of Therefore the integrals which appear in H can be approximated using the trapezoidal rule where λ denotes either w 2 , w or 1, and therefore it can be shown that We also need to take variations with respect to x i which occur in A −1 , and again it is simple to show that Finally from [1] we note that Hence, the full discretized form of Hamilton's equations to leading order are This gives a set of 2N equations for the 2N + 4 unknowns. The remaining 4 equations come from the boundary conditions The discretized set of equations can be written as where we define p = ( p, w 1 , . . . , w N +1 ) and q = (q, x 1 , . . . , x N +1 ). This form of the equations is amenable to time integration by a geometric integration scheme, namely the implicit-midpoint rule approach. In this case, the system becomes the set of nonlinear algebraic equations where n denotes the time-step, such that p n = p(n τ ). This system of implicit equations are solved at each new time step via Newton iterations. In order to increase the speed of the iteration scheme, the method of [18] is employed to iteratively calculate the inverse Jacobian matrix after the first iteration of the first time step.
In this paper, we consider the following initial condition from the linear solution which is derived in Sect. 2.2 from (2.30). The value of q is the initial displacement of the vessel from its equilibrium point, while q 1 and q 2 are chosen as independent parameters. When q 1 = q 2 = q the initial condition gives the linear solution derived in Sect. 2.2 and we can verify our results against the exact solution, when q 1 = 0 with q 2 and q independent, then we have the same initial condition as in [13], and thus, we can verify against their nonlinear results, and finally when q 1 = q 2 = 0, we have an initial condition akin to those achievable in an experiment, namely a horizontally displaced vessel released from rest with a quiescent fluid.

Numerical results
In this section, we present results of the numerical scheme for several sets of parameter values. In order to validate the numerical scheme, we compare our results both with the linear solution, and the nonlinear f-wave numerical scheme results presented in [13]. Once validated we then present results in the non-Boussinesq limit, a limit which the f-wave scheme struggles to resolve due to difficulties satisfying the system constraints. For the results presented we set N = 200 and τ = 10 −3 (although N = 50 and τ = 10 −2 are sufficient for the linear results).
Results are presented for the vessel evolution q(t) and the surface interface evolution h 1 (x 1 , t) along with time evolutions of the total vessel energy E v (t) and the total fluid energy E f (t) defined by It is also possible to show via simple algebraic manipulation that and thus E v + E f = H . Therefore the Hamiltonian is the total energy of the system, and thus, the discretized form of H The parameter values for the simulations presented in this section are given in Table 1.
The linear results in Figs Table 1.
The dots in both figures signify the results of [13], and the agreement is excellent. There are some minor discrepancies in the fluid interface profiles in Fig. 7, mainly close to the side walls, but these differences do not  The density ratio ρ 2 /ρ 1 = 0.7 in Figs. 6 and 7 is on the borderline between the Boussinesq and non-Boussinesq regimes. The f-wave numerical scheme developed by [13] works most effectively in the Boussinesq regime, especially for weakly nonlinear simulations. The scheme encounters problems satisfying the system constraints for density ratios ρ 1 /ρ 1 0.7. The Hamiltonian scheme developed here has the rigid lid and mass-flux conditions (2.1) and (2.16) directly built in to the scheme and so is able to resolve simulations for for these density ratios. Figs. 8 and 9 show results for ρ 2 /ρ 1 = 10 −3 , which is the density ratio of air to water for an initial condition akin to those found in an experimental setup, q 1 = q 2 = 0. As the initial interface is flat, the initial condition consists of an infinite sum of all the sloshing modes in (2.33) at different amplitudes, and thus, the result is the lowest frequency mode superposed with higher frequency modes, giving the small oscillations on the results. The energy error H N (t) in Fig. 8b, although larger than the result in Fig. 6b, is still relatively small O(10 −5 ), and bounded for the time-scale of the simulations. The results in Fig. 9 depict the interface gently sloshing back and forth in the vessel, and as it does so it becomes increasingly more fine scaled. This is a well known characteristic when

Fig. 6 a The vessel displacement q(t) and b E f (t), E v (t) and H N (t)
for the initial condition q 1 = 0, q 2 = q and the initial parameter values given in column 3 of  symplectic schemes are applied to sloshing problems [19] and is due to the energy of the system cascading down to the high frequency modes, in what is essentially a spectral scheme. However, as the numerical time integrator is symplectic, it conserves this energy and so this energy remains in the high frequency modes as these high frequency oscillations. These could be removed using an artificial viscosity term or the filtering scheme used by [20], but the numerical scheme will then no longer be energy conserving.
The two-layer results presented here show the capabilities of the Hamiltonian approach for these multilayer sloshing problems. Note, however, despite the introduction of the mapping φ(a 1 , τ ) to ensure x 1 (a 1 , τ ) = x 2 (φ(a 1 , τ ), τ ), this mapping was never discussed or plotted. This is because the two-layer problem is in fact a special case of the multilayer sloshing problem, because equations (2.1) and (2.16) mean that the upper-layer variables can be eliminated and the problem can be formulated solely in terms of lower-layer variables. In the next section, we formulate the general M-layer shallow-water problem, and present results for three-layer sloshing, where the mappings φ i do need to be calculated.

Fig. 8 a The vessel displacement q(t) and b E f (t), E v (t) and H N (t)
for the initial condition q 1 = q 2 = 0, q = 0 and the initial parameter values given by column 4 of Table 1 Fig . 9 The interface position h 1 (x 1 , t) for the results in Fig. 8

Extension to multilayer shallow-water flows
The extension of the theory to the M-layer shallow-water problem is straightforward, with the biggest difference being the necessity to calculate the mapping functions φ i (a 1 , τ ). The derivation and analysis can get a bit lengthy so detail is recorded in Appendix 3 for the three-layer case. A schematic of the general M-layer problem is shown in Fig. 10. Here we will impose the constraint x 1 = x 2 = ... = x M from the outset in order to simplify the analysis.
The ith layer mass conservation and momentum equations (for i = 1, ..., M) in the M-layer shallow-water equations with a rigid lid are 2) Fig. 10 Schematic of M-layer shallow-water sloshing in a moving rectangular vessel, with the constraint with the Eulerian rigid lid constraint The derivation of this multilayer system is given in Appendix 2. The associated generalised vessel equation to (2.6) is This multilayer shallow-water system is the Euler-Lagrange equation associated with the following Lagrangian functional in the Eulerian setting and p s and f j are Lagrange multipliers. In order to construct a geometric integration scheme such as that used in Sect. 2, we must first transform the equations from the Eulerian to the Lagrangian description, secondly construct a Lagrangian functional in the LPP description, and then Legendre transform to obtain the Hamiltonian form. To do this, we first note that as in Sect. which can be derived in the same way as for the two-layer case. As in the two-layer system, these two equations are used to eliminate u i and h i for one layer, which w.l.o.g, we choose to be the upper layer, with suffix M. Now introducing the LPP mapping (2.8) into the layer 1 mass conservation equation leads again to (2.11) and hence (2.15). Thus, h 1 is now written in terms of x 1 only, with u 1 = x 1τ its associated momenta. However, unlike the two-layer case, we still have layer variables (h 2 , u 2 ), . . . , (h M−1 , u M−1 ) to eliminate from the Lagrangian and replace by some position variable and its associated momenta.
If we now consider the constraint that τ ) is a mapping variable, then we can show that in each layer the mass conservation equation in the LPP framework reduces to ∂ ∂τ Similarly can be constructed, and the geometric integration scheme of Sect. 3 applied to it.

Numerical implementation for three layers
To demonstrate that the numerical scheme of Sect. 3 generalises to the M-layer problem, we present a result for coupled three-layer sloshing in Figs. 11 and 12. Details of the derivation of the three-layer Hamiltonian and symplectic integration scheme, as well as validation of the scheme, are given in Appendix 3. The initial conditions for these simulations are with the simulation parameter values given in column 5 of Table 1

Fig. 11 a The vessel displacement q(t) and b E f (t), E v (t) and H N (t)
for the initial condition (4.9) and the initial parameter values given in column 5 of Table 1 Fig . 12 The interface positions h 1 (x 1 , t) (lower curves in each panel) and x 1 x 1 whose amplitude is strongly modulated by the inclusion of the third layer. This modulated vessel displacement is due to the hydrodynamic force on the vessel walls slowly becoming out of phase with the restoring force of the spring, before slowly moving back in phase. This more complex behaviour is not a surprise as the characteristic equation for this system (8.40) has more solutions compared to the two-layer equation (2.33) due to the inclusion of the additional interface. The interface profiles again show fine scale structure at later times, but at t = 29 there exists fairly large oscillations at the lower interface. Also, the energy error H N (t) in Fig. 11b, while still small, O(10 −5 ), grows moderately over the time frame of the simulation. The reason for these two observations, we believe, is due to the Kelvin-Helmholtz instability on the interface [21], and we use a smaller time-step to stop the error growing more rapidly. This is more evident in the validation simulation in Appendix 3. Hence, one has to check the energy error H N (t) for a calculation to determine whether it is still within tolerable bounds. Again the introduction of artificial viscosity or filtering would help limit this instability by removing the fine-scale high-frequency modes from the system, which grow fastest in an inviscid system [22].

Conclusions and discussion
This paper documents the Lagrangian variational formulation of the LPP representation of multilayer shallowwater sloshing, coupled to horizontal vessel motion governed by a nonlinear spring. The Lagrangian variational formulation was transformed to a Hamiltonian formulation which has nice properties for numerical simulation. A symplectic numerical integration scheme was applied to the resulting set of Hamiltonian partial differential equations for the two-layer problem, and results of the simulations were found to be in excellent agreement with the linear solution and the nonlinear results of the f-wave scheme of [13]. Using this Hamiltonian formulation the results of [13] were extended into the non-Boussinesq regime, with a result presented for a density ratio ρ 2 /ρ 1 = 10 −3 , akin to that of air over water.
The Hamiltonian formulation was presented in detail for the two-fluid system, but the solution procedure was generalised in Sect. 4 to a system of M-fluid layers coupled to horizontal vessel motion where the vessel is attached to a nonlinear spring. Results were presented for a three-layer system, with the full derivation confined to Appendix 3. Results for the three-layer system showed a system energy error which grew slowly in time, due to the Kelvin-Helmholtz instability on the fluid interfaces. For the results presented in this paper, this error growth was small and thus tolerable for the time frame of the simulations. However, this error would need to be examined in fully nonlinear simulations or long-time simulations to make sure it was small compared to the fluid velocities and vessel displacement. Also, in thin layers, where fluid velocities tend to be larger to conserve the mass flux (4.6), this instability could be an issue. Surface tension or a filter could be added to mollify the instability.
As this work was motivated by studying the WEC of Offshore Wave Energy Ltd (OWEL), a direction of great interest is to extend the vessel motion to incorporate rotation (pitch) along with the translations considered here, and to incorporate influx-efflux boundary conditions at the side walls, which model the waves entering the device and leaving through the extraction route. In the OWEL WEC, the wetting and drying of the upper rigid lid is very important for the modelling of the power-take-off mechanism. The current two-layer approach considered in this paper cannot incorporate this phenomena. The reason for this comes from the mass-flux equation h 1 u 1 + h 2 u 2 = 0 which holds throughout the fluid. We find that as h 2 → 0 in this expression, despite the momentum h 2 u 2 reducing in size, the value of u 2 becomes large which causes numerical difficulties in the current scheme. Thus, another area of great interest is to incorporate this feature into the model.
where I = L 0 Bw da, and hence Making these substitutions into (6.1) leads to where λ represents w or w 2 . Also, Making these substitutions into (6.2) above, along with noting that gives the required form of (2.23).

Appendix 2: Derivation of multilayer shallow-water equations with a rigid lid
In this section, we summarise the derivation of the multilayer shallow-water equations given in (4.1) and (4.2). Consider the same M-layer system as considered in Sect. 4, so in the ith layer, the two-dimensional Euler equations are where i = 1, ..., M, x, y, t subscripts denote partial derivatives and the over dots represent ordinary derivatives with respect to time. For simplicity, we drop the i subscripts from the layer coordinates (x i , y i , t). In the shallowwater regime, the common assumptions on the flow field are that Dv i /Dt ≈ 0 and u iy = 0 [23]. Under these assumptions, the vertical momentum equation reduces to p iy = −gρ i , which can be integrated from a general y-value in the layer to the upper-surface, denoted by H i = i j=1 h j , to give the pressure in each layer as Here P i = p i (x, H i , t), P N = p s (x, t) and H 0 = 0. Back substitution from the rigid lid to eliminate the P i expressions gives the i th layer pressure as Substituting this into the horizontal momentum equation with u iy = 0 gives (4.2).
The conservation of mass equation is derived in the usual way by integrating the continuity equation across the fluid layer, noting that u iy = 0, so for i = 1, · · · , M. Finally using the kinematic boundary condition on each interface and noting that H i − H i−1 = h i leads to (4.1) for each layer i.

Lagrangian formulation
From (4.1) and (4.2), the three-layer shallow-water equations are where we have assumed the constraint that x 1 = x 2 = x 3 . The interface thicknesses, h i (x 1 , t), are constrained by the rigid lid constraint The motion of the vessel is governed by f and m 1 is the fluid mass in the i th layer. Equation (4.5) gives the form for the M-layer shallow-water Lagrangian in the Eulerian framework. Thus, for the three-layer system this Lagrangian is Note that f j , j = 1, 2, 3 and p s act as Lagrange multipliers. We can eliminate the top-layer variables h 3 and u 3 by using the constraint on the layer thicknesses (8.5) and the conservation of mass flux (4.6) which give To write the Lagrangian (8.7) in terms of the LPP formulation, we again consider the mapping (2.8) and we drop the subscript 1 on a, for simplicity. The constraint in the middle and upper layers become x 1 (a, τ ) = x 2 (φ 2 (a, τ ), τ ), and x 1 (a, τ ) = x 3 (φ 3 (a, τ ), τ ), where φ 2 (a, τ ) and φ 3 (a, τ ) are mappings. As for the two-layer case we do not have to consider the mapping function φ 3 (a, τ ) because h 3 and u 3 are eliminated using (8.8). However, this time we do need to determine the mapping φ 2 (a, τ ) as part of the solution procedure.

Numerical algorithm
To discretize the equations (8.15)-(8.20), we use the same approach as for the two-layer problem. Thus, we discretize the Lagrangian state space into N parcels using (3.1) and let τ ). We again use the semi-discretization method to derive the resulting discretized form of Hamilton's equations and rather than write out the whole form of the discretized equations, we note that they can be derived from (8.15)-(8.20) by noting (μ) a := 1 a (μ i − μ i−1 ), in the semi-discretization. See (2.23) and (3.3).
The initial conditions used to validate of the scheme are the linear forms from Sects. (8.41)-(8.43) where U 1 (a) and U 2 (a) are given by (8.36) and (8.37) with x replaced by a, and Q 1 is an independent parameter. When Q 1 = 1, the initial condition is that given by the linear problem, while when Q 1 = 0 it is an initial condition achievable in an experiment, namely a horizontally displaced vessel released from rest with a quiescent fluid.   Table 2.
Note that for the three-fluid system, the vessel energy is given by while the fluid energy is ρ [1,3] χ 2 1 + ρ [2,3]  , where H N is the discretized form of the Hamiltonian (8.14). The slight increase in the error over the duration of the simulation is believed to be due to the Kelvin-Helmholtz instability, but this error growth is not large enough to affect the result and hence is tolerable.
The nonlinear results in Figs. 15 and 16 compare directly to the two-layer simulation in Figs. 8 and 9. This result is given by the red dots in Figs. 8a and 9. The comparison with the two-layer result is excellent, with the slight discrepancy in the two results at large times due to the two sets of simulation parameter values not being identical (ρ 2 = 990 kg m −3 not 1000 kg m −3 ). The energy error in Fig. 15b again grows slightly in time, due to the Kelvin-Helmholtz instability, but as it is O(10 −7 ), it is tolerable for the results presented.
t Fig. 13 a The vessel displacement q(t) and b E f (t), E v (t) and H N (t) for the linear initial condition Q 1 = 1 and the initial parameter values given in column 1 of Table 2. The dots represent the linear solution (8.26) for the lowest-frequency mode  As a final note, the growth in the energy error, which we believe is a consequence of the Kelvin-Helmholtz instability can also be observed in the two-layer simulations when ρ 1 ≈ ρ 2 through a growing energy error H N (t), but in this case, the growth is not as obvious as in the three-layer simulations presented here.