Edinburgh Research Explorer One-dimensional and two-dimensional Green-Naghdi equations for sloshing in shallow basins

This paper presents a verified model of weakly non-linear wave sloshing in shallow basins, based on level I Green – Naghdi (GN) mass and momentum equations derived for mild-sloped beds. The model is verified for sloshing of an initially sinusoidal free surface perturbation in a square tank with a horizontal bed. The model is also used to investigate free surface sloshing of an initial Gaussian hump in closed square basins, over horizontal and non-uniform bed topographies. Analysis of the free surface slosh motions demonstrates that the model gives predictions in satisfactory agreement with the analytical solution of linearised shallow water theory obtained by Lamb. Discrepancies between GN predictions and linear analytical solutions arise from the effect of wave non-linearities arising from the wave amplitude itself and wave – wave interactions.

This paper presents a verified model of weakly non-linear wave sloshing in shallow basins, based on level I Green-Naghdi (GN) mass and momentum equations derived for mild-sloped beds. The model is verified for sloshing of an initially sinusoidal free surface perturbation in a square tank with a horizontal bed. The model is also used to investigate free surface sloshing of an initial Gaussian hump in closed square basins, over horizontal and nonuniform bed topographies. Analysis of the free surface slosh motions demonstrates that the model gives predictions in satisfactory agreement with the analytical solution of linearised shallow water theory obtained by Lamb. Discrepancies between GN predictions and linear analytical solutions arise from the effect of wave non-linearities arising from the wave amplitude itself and wave-wave interactions.
Notation a amplitude a t ij ; b t ij ; c t ij tridiagonal matrix coefficients a z b bed amplitude F vector derived from Green-Naghdi equations f frequency h total water depth h 0 still water depth i, j grid points in x-and y-directions k wave number L x basin length L y basin width P pressure r, s wave modes T period t time u, v horizontal velocity components w vertical velocity component z b bed elevation Δt time interval Δx, Δy grid size in x-and y-directions ζ free surface elevation above still water ζ c crest-induced free surface elevation ζ t trough-induced free surface elevation η free surface elevation above horizontal datum λ n shape function ρ water density σ pq Kronecker's delta function ϕ phase ω wave angular frequency

Introduction
Slosh motions occur in liquid ballast tanks of ships, liquefied natural gas tankers and containers subject to seismic excitation. The response of the liquid free surface depends on the amplitude of the initial disturbance and the natural frequencies of motion, and involves complex fluid-structure interactions (Ibrahim, 2005). In practice, civil engineers may use sloshing tests to assess the likelihood of such phenomena occurring in oil tanks, elevated water towers and in reservoirs. According to Sarpkaya and Isaacson (1981) the earliest theories of progressive and sloshing waves were developed by Airy and Stokes, which were based on potential theory idealisations. Airy's theory was linear, so that it was strictly derived for waves of zero amplitude. Stokes extended the theory to deal with waves of finite amplitude, and he obtained series solutions that were later computed to high order by many hydrodynamicists in the late twentieth century (for more details see Sarpkaya and Isaacson (1981)). To extend to more realistic domains, computational methods have become widely used. These include: (a) boundary-element and finite-element potential flow solvers; (b) computational fluid dynamics-Navier-Stokes solvers with volume of fluid treatment of the free surface, Navier-Stokes solvers with level-set treatment of the free surface, Navier-Stokes solvers with mappings of the free surface; and (c) smoothed-particle hydrodynamics. The above-mentioned three-dimensional (3D) computational methods undergo inherently high computational expense; therefore, considerable effort has gone into depth-averaged approaches, they being cheaper to compute and yet capturing much of the physics. Examples of such approaches include shallow water equations (SWEs) (e.g. Lamb, 1916), and more recently Green-Naghdi (GN) equations. Non-linear SWEs are the depth-averaged form of continuity and Navier-Stokes momentum equations. Since SWEs neglect vertical motions and the consequent hydrostatic pressure, these equations are usually restricted to long-wave behaviour. Thus, SWEs are limited to shallow depth (surf zone) of the ocean (Bonneton et al., 2011;Nadiga et al., 1996). Green and Naghdi (1976) pioneered the development of non-linear equations for two-dimensional (2D) incompressible inviscid fluid sheets. Green and Naghdi (1976) proposed a theory of fluid sheets known as GN theory to model the 2D continuum of unsteady inviscid 3D flows. The theory facilitated prediction of unsteady, non-periodic, free surface flows. GN theory utilises some aspects of perturbation analysis in building up first-, second-and higher-order approximations (called levels) to layer-averaged mass and momentum equations. According to Webster and Shields (1991), the GN approach assumes a particular flow kinematic structure in the vertical direction for shallow-water problems. The fluid velocity profile is a finite sum of coefficients depending on space and time multiplied by a weighting function. GN fluid sheet theory reduces the dimensions from three to two, yielding equations that can be solved efficiently so that no scale is introduced and no term is deleted (Webster and Shields, 1991). Nevertheless, the lowest level of GN theory permits the kinematic boundary conditions to be satisfied. There are two types of GN theory: restricted and unrestricted. The former successfully models irrotational shallow-water flow field. Restricted GN theory was derived from the first level of the direct theory by means of a constrained director (Shields and Webster, 1988). Later, this procedure was extended to the kth level theory (Demirbilek and Webster, 1992). In other words, in a restricted GN theory, the k components of the 2D velocity components are constrained. Demirbilek and Webster (1992) developed an unrestricted version of GN theory of shallow water by enforcing conservation of mass and momentum in the vertical direction and implementing exact boundary conditions. They demonstrated that GN theory can appropriately predict the behaviour of a non-linear numerical wave tank. According to Webster and Shields (1991), GN sheet theory is placed between classical perturbation methods and pure numerical schemes. Webster and Shields (1991) note that for classical perturbation methods, there is usually no evidence that the assumed series is convergent. However, in certain flow problems, such as 2D water waves in both shallow and deep water addressed by GN theory, there is ample evidence of convergence. There is another difference between classical perturbation methods and GN theory. The former exactly satisfies field equations but partially satisfies boundary conditions. Thus, the classical perturbation methods show inconsistencies, while the GN theory is self-consistent since the boundary conditions are exactly satisfied, when the field equations are partially approximated. With regard to coastal engineering, the present GN model is merely applicable to the shallow depth (before the surf zone) and the intermediate depth of ocean (Jalali, 2016). GN levels higher than I can be used for the deep water ocean (Webster and Shields, 1991). The present paper describes the application of 2D and one-dimensional (1D) level I GN equations (according to Jalali (2016)) to free surface sloshing in closed square basins, with horizontal and non-uniform bed topographies. Sections 2 and 3 present the derivation of governing equations. Sections 4, 5 and 6 outline the numerical implementation. Section 7 presents the results for initial sinusoidal and Gaussian free surface perturbations. Section 8 is the summary of the main findings.

Derivation of continuity equation of level I GN equation
In a 3D Cartesian system, the classical mass conservation equation is applied to drive the GN continuity equation In which, u, v and w are the velocity components in x, y and z directions. In the present derivation of GN equations, the total depth, h, is h = h 0 + ζ. Here h 0 is the still water depth and ζ is the free surface elevation above still water level. The elevation of free surface above the fixed horizontal datum, η, is η = h + z b .
Here, z b is the bed elevation above a fixed horizontal datum. The kinematic bed and free surface boundary condition are as follows 2: where d| b is the kinematic bed boundary condition and d| s is the kinematic free surface boundary condition. To derive the GN continuity equation, it is assumed that the velocity vector, V, can be written as follows W n x; y; t ð Þλ n z ð Þ where W n = (u n , v n , w n ) is a vector of velocity component approximations at level n; λ n are the assumed shape functions depending on the z-direction and e is the level of approximation of GN theory. Expansion of Equation 4 for level I gives the following velocity parameters 5: and u 1 (x, y, t) = v 1 (x, y, t) = 0 (for more details see Demirbilek and Webster (1992)).
By applying Equation 5 in Equations 1 and 2, the values of w 1 and w 0 are obtained 6: By using Equation 5 in Equation 3 and replacing Equations 6 and 7 into the resulted equation, the 2D level I GN continuity equation is derived in which h is the total depth and (u 0 , v 0 ) are the horizontal velocity components at a particular point (x,y) and time t.

Derivation of x-direction momentum of level I GN equation
The general momentum conservation equation expanded in the x-direction is as follows 9: @ρu @t þ @ρuu @x þ @ρuv @y þ @ρuw @z ¼ À @P @x in which ρ is the water density. Depth integrating Equation 9 and then using the chain rule for the fourth term and applying the Leibnitz rule for the right-hand side term result in 10: where P n ¼ @ Ð η zb Pλ n dz =@x;P is pressure at the free surface (hereP ¼ 0) andP is the pressure at the bottom. Applying Equation 4 in Equation 10 11:

12:
X e¼1 r¼0 @u r @x λ r þ X e¼1 r¼0 @v r @y λ r þ X e¼1 r¼0 w r λ n 0 ¼ 0 ρu m w r y m nr ¼ À @P n @x ÀP @z b @x λ n j zb Derivation of the GN equation in z-direction is similar to the derivation of GN equation in x-direction; therefore, the detailed derivation is not included here for brevity. The constrained z-direction momentum equation of level I GN equation is 15: ρ w m w r y m rn ¼Pλ n j z b þ P n 0 À ρgy n where P 0 n ¼ Ð η z b Pλ n 0 dz and g is the gravitational acceleration.
(To learn more about the derivation of GN momentum equations in this paper, refer to the PhD theses by Haniffah (2013) and Jalali (2016).) Equations 14 and 15 are expanded for level I {(u 0 , u 1 ), (v 0 , v 1 ), (w 0 , w 1 )}. In z-momentum GN Equation 15, w 0 and w 1 are replaced by Equations 6 and 7. Two sets of equations are derived for x-and z-momentum GN equations (n = 0, level 0, and n = 1, level I). In z-momentum GN Equation 15, P 0 0 ¼ 0 and P 0 1 ¼ P 0 . For eliminating the effect of P 0 differentiation is applied with respect to x for n = 1 set of z-momentum GN equation and the obtained equation is added to the n = 0 set of x-momentum GN equation. Moreover, the effect of the bottom pressure term,P, is eliminated by multiplying (∂z b /∂x) the n = 0 set of z-momentum GN equation and adding the resulted equation to the n = 1 set of z-momentum GN equation. Finally, the 2D level I x-momentum GN equation for the stationary bed, (∂z b /∂t = ∂ 2 z b /∂x∂t = ∂ 2 z b /∂y∂t = 0), and non-uniform bathymetry is derived By simplifying Equation 16, the 2D level I x-momentum GN equation for uniform bathymetry is 17: The 2D level I y-direction momentum GN equation and its derivation are not included here for brevity (for complete detailed derivations of level I GN equations, see the PhD theses by Haniffah (2013) and Jalali (2016)). The corresponding 1D level I x-momentum GN equation for uniform bathymetry is 18: The corresponding 1D level I GN continuity equation is Engineering and Computational Mechanics Volume 170 Issue EM2 One-dimensional and two-dimensional Green-Naghdi equations for sloshing in shallow basins Jalali and Borthwick

Numerical implementation
To develop a valid numerical solver of 2D level I GN equations, the researchers discretised Equations 8 and 17 using second-order central finite differences (numerical solver of 2D level I GN momentum equation for non-uniform bathymetry, Equation 16 is presented in chapter 3 of the PhD thesis by Jalali (2016)). In the present numerical solver, two different sections were developed to deal with the continuity equation and momentum equations. Continuity equations do not have any cross-derivative terms; therefore, an explicit second-order finite difference solves the equations properly. On the other hand, momentum equations contain cross-derivative terms (∂ 2 u 0 /∂x∂t and ∂ 3 u 0 /∂x 2 ∂t). Since an explicit predictor-corrector scheme is incapable of solving this kind of equation, an implicit finite difference scheme is used to solve the GN momentum equations. Herein, an implicit tridiagonal matrix inversion scheme is utilised to solve GN momentum equations (for more details see chapter 3 of the PhD thesis by Jalali (2016)). Nevertheless, other numerical schemes such as finite volume or finite element may be capable of solving the continuity and momentum equations simultaneously. The developed numerical solver is based on finite difference that cannot solve these equations simultaneously. The 2D GN continuity Equation 8 is discretised by applying second-order central difference 20: The 2D level I GN x-direction momentum Equation 17 is rearranged and solved using an implicit scheme. To handle Equation 17, let 21: In second-order central differences this becomes 22: The discretised equation is rewritten as where a t ij , b t ij and c t ij are the tridiagonal matrix coefficients for x-direction. Here, i refers to x-direction; j is y-direction and t is time. F is also equal to the remaining spatial derivative terms in Equation 17. These terms are also discretised in Equation 17 by using second-order central differences, giving 24: The above sets of discretised equations form the tridiagonal matrix system. The unknown valuesû t 0ij are obtained for j = 2,…,j max −1 and i = 1,…,i max using the Thomas algorithm (Press et al., 2007). It should be mentioned that i max and j max refer to the maximum number of grid points in x-and y-directions. A similar numerical approach was applied to 2D level I GN y-momentum equation (for a detailed development of GN numerical solver refer to chapter 3 of the PhD thesis by Jalali (2016)). Iteration is then used to centre correctly (in space and time) the cross-derivative terms that appear in both the x-and y-direction momentum equations. Runge-Kutta fourth-order time integration is used to update the total depth and horizontal velocity components at each time step.

Boundary conditions
To solve the 2D GN equations, it is necessary to impose flexible and compatible boundary conditions. For instance, solid wall boundaries are located at the ends of the domain when simulating sloshing of waves in a tank. The surface elevation at the boundary obtained by cubic Lagrange interpolation of interior values is assigned according to Haniffah (2013). The velocity is set to zero at solid wall boundaries. Additional ghost grid points are located outside the boundaries, with antisymmetry imposed for horizontal velocity (u) on y-direction and symmetry on x-direction.

53
Engineering and Computational Mechanics Volume 170 Issue EM2 One-dimensional and two-dimensional Green-Naghdi equations for sloshing in shallow basins Jalali and Borthwick The velocity boundary conditions on x-direction for numerical solver of 2D GN equations are Here, index 1, i max and j max present the location of the wall; index 2, 3, i max −1, i max −2, j max −1 and j max −2 indicate the grid points located inside of the computational boundary (before the wall); and index −1, 0, i max + 1, i max + 2, j max + 1 and j max + 2 show ghost points located outside of the computational boundary (after the wall). The velocity boundary conditions in y-direction are not included here for brevity (for more details see Jalali (2016)). The symmetry boundaries imposed for surface elevation are @h @t In 1D level I GN numerical solver, Equations 18 and 19 are the governing equations. The 1D GN continuity Equation 19 is discretised by applying the second-order central difference In Equation 18, cross-derivative terms of x and t are rearranged and solved by using an implicit scheme, similar to the procedure followed for Equation 21. The discretised equation is F is also equal to the remaining spatial derivative terms in Equation 18. Using second-order central differences discretisise the terms in Equation 18 and gives 27: The velocity boundary conditions for solver of 1D GN equation are In 1D GN numerical solver, the symmetry boundaries imposed for surface elevation are @h @t For more details, see chapter 3 of the PhD thesis by Jalali (2016) and chapter 4 of the PhD thesis by Haniffah (2013).

Numerical procedure
The GN program comprises four main subroutines: input, calculation, update and output. For each test case, the following initial values are put into the program: bed elevation, initial water depth, amplitude, length and width of the study basin, number of grid points, time step and duration of simulation time. The initial conditions supply for the bed elevation above fixed horizontal datum, local depth and local horizontal velocity components throughout the tank. Solving the discretised continuity equation provides new water depth values throughout the grid. Then, the discretised momentum equations are solved. Iteration is used to solve the cross-derivative velocity terms in the momentum equation. Next, boundary conditions are invoked. The values of u, v and h are updated after each time step. The calculation process is repeated until the simulation is complete. The two benchmark tests comprise: sloshing in a square tank and free surface sloshing of an initial Gaussian hump in a square basin. In sloshing in a square tank to test for grid independence, the time history of wave elevation at the corner of the tank (in positive x-direction) was obtained on grids of increasingly fine resolution (Δx = 25 m (coarse grids), Δx = 10 m (medium grids) and Δx = 1 m (fine grids)) and a fixed time step Δt = 1 s. The results demonstrated that Δx = 1 m was sufficient to achieve a converged solution. Three time steps are chosen (Δt = 0·25, 1·0 and 2·0 s) on the converged grid with Δx = 1 m. There was close agreement between the results; therefore, Δt = 1 s was selected as a fixed time step in the simulation of sloshing in the tank (see chapter 4 of the PhD thesis by Jalali (2016)). To determine the number of grid points required to produce an accurate simulation of free surface sloshing of an initial Gaussian hump, grid convergence test was performed.

Sloshing in a tank
First, the benchmark test of sinusoidal free surface sloshing in a square tank is considered. The wavelength, L, is 1000 m and the still water depth, h 0 , is 5 m. Sloshing motions may even occur by using a very small number of amplitude disturbance. For the present test case the amplitudes are a = 0·005 and 0·05 m. These small numbers of wave amplitude are applied in order to create minimum non-linear behaviour by the sloshing wave. The first-order analytical solution for the depth profile evolution in space and time of a standing wave in a tank (e.g. Dean and Dalrymple, 2004) is Here, ζ c refers to the crest-induced free surface elevation time series, a is the amplitude of the standing wave, k is the wave number, ω is the angular frequency of the wave, x is the distance along the tank, t is time and ϕ is the phase. Wave angular frequency ω is obtained by means of dispersion relation In this case: ω = 0·044 rad/s, T (period) = 2π/ω = 142·8 s and f = (frequency) = 0·007 Hz.
One-dimensional and two-dimensional Green-Naghdi equations for sloshing in shallow basins Jalali and Borthwick to deal with sloshing in the tank, so it is necessary to measure free surface elevation of time history at two different locations: (a) in the corner of the tank and (b) at the centre of the tank. Figure 1(a) presents crest-induced water level time histories in the corner of the tank (in positive x-direction) for a sloshing wave of a small-amplitude disturbance a = 0·005 m with phase ϕ = 0. Here, ζ c /a = 1 at time t = 0 s for free surface elevation of time history in the corner of tank. Excellent agreement is obtained between the first-order analytical solution and the numerical prediction in which cross symbols (numerical prediction) essentially overlay the solid line (analytical solution). The standing wave behaviour is periodic and of constant amplitude. This case verifies that the numerical scheme gives a correct representation of the underlying mathematical description, provided the waves are nearly linear. Figure 1(b) depicts the numerical prediction of the crest-induced free surface elevation time history in the corner of the tank for a = 0·05 m. The free surface elevation time history displayed in Figure 1(b) is shorter than that in Figure 1(a) because non-linear effects eventually cause shock-like steepening of the wave profiles (becoming visible at about t = 2000 s) in the larger amplitude case leading to the numerical model becoming unstable. A shock-capturing scheme would be needed to overcome this problem, and is recommended for future implementation. Figure 1(c) presents the numerically predicted trough-induced free surface elevation time history at the corner of the tank in which ζ t = −acos(kx)cos(ωt). Qualitatively, the results are almost the same as for the crest-induced case (i.e . Figure1(b)). Non-linearity can be presented by even harmonics (ζ c + ζ t /2a).
To separate even harmonics, harmonics are treated as orthogonal functions. Figure 1(d) shows the numerically predicted free surface elevation time history of the even harmonic components (obtained by taking the average of results obtained for ϕ = 0 and π) for the amplitude a = 0·05 m. Here, the amplitude of the second-order harmonics grows monotonically with time until the point at which numerical instability occurs. It is worth noting that the (linear) analytical solution is not able to show the non-linear behaviour of even harmonics. The developed numerical model is also applied to simulate other possible predictions for a sloshing wave in the centre of the tank. In this case, the wavelength is 1000 m and the still water depth h 0 is 1 m. Two values of wave amplitude, a = 0·005 and 0·015 m, are selected. In this case: ω = 0·0197 rad/s, T = 319·28 s and f = 0·003 Hz.  Figure 2(b) depicts the numerically predicted crest-induced free surface elevation time history for a = 0·015 m. It is clear that the numerical solver is unable to simulate the long-term sloshing behaviour of the wave, with high-order even harmonic oscillations appearing after t = 1735 s. The non-linear effects eventually caused shock-like steepening of the wave profiles. Therefore, a shock-capturing scheme is required to overcome this problem. Figure 2(c) shows the numerically predicted free surface elevation time history of the trough-induced sloshing at the centre of the tank for a = 0·015 m. The results are qualitatively almost the same as for the crest-induced case (i.e. Figure 2(b)). Figure 2(d) shows the numerically predicted free surface elevation time history of the even harmonic components (obtained by ζ c + ζ t /2a) for a = 0·015 m. The effect of non-linearity increases as the initial slosh amplitude increases, as would be expected (for more details see chapter 4 of the PhD thesis by Jalali (2016)).

Free surface sloshing of an initial Gaussian hump in a closed square, flat-bottomed basin
The numerical solver of the 2D GN equations is now verified for non-linear free surface sloshing motions arising One-dimensional and two-dimensional Green-Naghdi equations for sloshing in shallow basins Jalali and Borthwick from an initial Gaussian hump free surface profile in a closed basin. The well-established analytical solution (Lamb, 1916;Wei and Kirby, 1995;Yao, 2007) of the linearised SWEs for the evolution of an initial hump in a closed square basin is    Figure 5. Comparison between analytically predicted fast Fourier transform (FFT) spectrum for the free surface elevation time history of the initial Gaussian hump with the numerically predicted FFT spectrum for the free surface elevation time history of even harmonics component 57

Engineering and Computational Mechanics
One-dimensional and two-dimensional Green-Naghdi equations for sloshing in shallow basins Jalali and Borthwick Here, ζ is the free surface elevation above still water level, where ; L x and L y are the length and width of the basin; ω is angular frequency; p and q are the number of wave components; σ pq is the Kronecker delta  One-dimensional and two-dimensional Green-Naghdi equations for sloshing in shallow basins Jalali and Borthwick function and 32: where a is the wave amplitude and b is the spreading parameter.
Consider a basin of 7·5 m length and 7·5 m width. The constant water depth is h 0 = 0·45 m. The initial amplitude of the hump a is 0·045 m and the spreading parameter b = 2 m −2 . To obtain an accurate estimate of the analytical solution, different numbers of wave components ( p and q) and grid size (Δx and Δy) are selected to solve the double Fourier series in Equation 30. Table 1 presents the initial elevation of analytically predicted free surface perturbation ζ c /a at the centre of the basin at time t = 0 s for different numbers of wave components ( p, q) and grid sizes (Δx, Δy). This table shows that using p = q = 51 and Δx = Δy = 0·075 m is sufficient to obtain a converged analytical solution. In the present GN numerical solver, four iterations on the medium grid size, Δx = Δy = 0·0375 m, are sufficient for the numerical predictions to be in satisfactory agreement with the analytical results. A time step of Δt = 0·05 s is found to be sufficient to achieve accurate and stable results. Figure  One-dimensional and two-dimensional Green-Naghdi equations for sloshing in shallow basins Jalali and Borthwick the initial hump. Here, ζ c /a = 1 at time t = 0 s for free surface elevation of time history at the centre of the basin. There is close agreement between the numerical and analytical results for about 10 s after the initial hump is released, after which differences are discernible between the numerical predictions and analytical solution. The foregoing discrepancies are largely due to non-linear (second-and higher-order) wave interactions which are modelled by the 2D level I GN equations, but neglected in the analytical solution of the linearised SWEs. The even harmonics of the sloshing motions induced by the initial Gaussian hump can be determined by simulating the free surface time series resulting from releasing the initial hump, and the corresponding free surface time histories driven by an initial trough of equal but opposite shape to that of the hump (following the separation of harmonics method utilised by Borthwick et al. (2006), Hunt et al. (2004 and Johannessen and Swan (2001), among others). Here, the harmonics are treated as orthogonal functions, and the even harmonics obtained by addition as (ζ c + ζ t /2), where ζ c refers to the free surface elevation time series of the initial Gaussian hump and ζ t the equivalent time series for the initial Gaussian trough. Figure 4 shows that for a relatively large amplitude hump One-dimensional and two-dimensional Green-Naghdi equations for sloshing in shallow basins Jalali and Borthwick (a = h 0 /2 = 0·225 m), it is possible to see evidence of the nonlinear effect produced by even harmonics, which are nondimensionalised with respect to the amplitude of the initial hump. The even harmonics have amplitudes of up to about 20% of that of the initial hump, and are perhaps growing slightly over the duration of the simulation. Figure 5 presents a comparison of the analytically predicted fast Fourier transform (FFT) spectrum for the free surface elevation time history of the initial Gaussian hump and the numerically predicted FFT spectrum for the free surface elevation time history of even harmonic components. It can be observed that all five peaks of numerical even harmonics occur at the same frequency as that of the analytically predicted peaks of Gaussian hump. Table 2 lists the resonant frequencies associated with different modes for the basin (r and s): the first peak occurs at mode r = 2 and s = 0; the second peak at mode r = 2 and s = 2; the third peak at mode r = 3 and s = 2; the fourth peak at mode r = 4 and s = 0; and finally the fifth peak at mode r = 5 and s = 1. Figures 6 and 7 compare the numerical simulations with the analytical solutions for sloshing in the basin, using 3D visualisations of the water surface at times t = 1, 5, 10 and 20 s. One-dimensional and two-dimensional Green-Naghdi equations for sloshing in shallow basins Jalali and Borthwick The corresponding contour maps are given in Figures 8 and 9. Although satisfactory agreement is achieved between the numerical predictions and analytical solution at t = 1 and 5 s, discrepancies between the numerical and analytical simulations become evident at t = 10 s, and grow with simulation time, as can be seen at t = 20 s where phase differences are observable. Figure 10 shows the numerically predicted velocity vectors and magnitude contours for the water surface at times t = 1, 5, 10 and 20 s after releasing the Gaussian hump in the flat-bottomed basin. Here, the velocity vectors indicate the direction of water particles and magnitude contours show the value of velocity in different sections of the basin. Slosh motions evolve in the basin from the initial hump as it rapidly drops under its own weight, causing a deep trough at the centre of the basin with an associated circular wavefront. The initial free surface motions are remarkably similar to those generated by the collapse of a liquid column, as modelled by Toro (2001), among others, except that the central oscillations do not die away as quickly. The balance between potential and kinetic energy drives repeated up and down motions at the centre of the basin, generating circular waves that propagate radially away from the centre of the basin and reflect with the basin walls. The repeated reflections between the waves with each other and the walls promote increasingly complicated sloshing modes dominated by waves whose wavelength is half the length of the basin. One-dimensional and two-dimensional Green-Naghdi equations for sloshing in shallow basins Jalali and Borthwick simulations is now considered. The Gaussian hump is released at t = 0 s and the numerical solution then propagated forward in time until 20 s, after which the time step is made negative and the numerical scheme is forced to simulate the backward propagation of the water surface until time zero is again reached. The results should be in almost identical agreement, given that the problem is thermodynamically reversible. There is no viscosity present, no turbulence, no surface tension and no sources of friction (e.g. from the basin walls or bed). Figure 13 examines reversibility by plotting the free surface elevation time history at the centre of the basin. The forward part of the simulation is shown by the solid line and the backward part by means of cross symbols. For the vast majority of the simulation, the agreement between the forward and backward processes is excellent. Discrepancies only occur when recovering the last 3 s of simulation. It seems likely that by travelling forward and backward in time, the accumulated dissipative error of the numerical scheme is responsible for causing the recovered hump to lose amplitude. Increasing the number of grids, increasing the number of iterations and selecting very One-dimensional and two-dimensional Green-Naghdi equations for sloshing in shallow basins Jalali and Borthwick small time steps did not reduce the magnitude of this accumulated dissipative error at the last 3 s of simulation. Figure 14 indicates the effect of reversibility on free surface elevation time history simulation at the centre of the basin for very small amplitude disturbance, a = 0·001 m, and spreading parameter b = 0·2 m −2 . Excellent agreement is obtained between forward and backward propagations of the water surface since the effect of non-linearity is neglected.

7.3
Parameter tests for sloshing in a square basin with non-uniform bathymetry Simulations are now considered of the sloshing behaviour of an initial Gaussian hump of water released in a square basin with non-uniform bathymetry. Here, the bed elevation is given by where a z b is the bed amplitude and b z b ¼ 2 m À2 is a bed spreading parameter. The initial local free surface elevation is 34: where a is the amplitude of the initial Gaussian hump in free surface elevation and b = 2 m −2 is a measure of its spread. Figures 15 and 16, respectively, depict 3D visualisation and contour maps of the water surface at times t = 1, 5, 10 and 20 s for relatively large values of bed hump amplitude (a z b ¼ 0Á225 m) and Gaussian hump amplitude (a = 0·225 m). It is worth mentioning that the bed hump amplitude, a z , is obtained through trial and error method (for more details see chapter 5 of the PhD thesis by Jalali (2016)). Figure 17  One-dimensional and two-dimensional Green-Naghdi equations for sloshing in shallow basins Jalali and Borthwick presents velocity vectors and magnitude contours for the water surface at times t = 1, 5, 10 and 20 s, where the bed topography contains a central hump. The effect of the bed hump on the evolution of the water free surface is most obvious at the centre of the basin exactly where the bed hump has its peak. At first, the Gaussian free surface hump drops rapidly to form a trough at the centre of the basin, releasing a circular ring-like wave that propagates towards the basin walls, where reflections occur. The plunging free surface at the centre of the basin interacts with the bed hump, leading to the recovery of a second clapotis-like hump, which peaks and releases a second circular wave. After several cycles of central peaks and troughs, the water surface motions immediately above the hump degenerate into a patch of small waves that heave up and down over the hump (after t $ 10 s); elsewhere the sloshing behaviour is similar to that of the corresponding case without a bed hump, particularly the presence of sloshing components whose wavelength is half the length of the basin. Figures 18 and 19, respectively, show 3D visualisation and contour maps of the evolution of the water surface over a Gaussian trough in the bed (a zb ¼ À0Á225 m). Figure 20 shows the velocity vectors and magnitude contours for the water surface at (a) t = 1 s, (b) t = 5 s, (c) t = 10 s and (d) t = 20 s, where the bed topography contains a central trough. The bed trough has the greatest effect at the centre of the basin, coincident with the peak position of the initial Gaussian free surface One-dimensional and two-dimensional Green-Naghdi equations for sloshing in shallow basins Jalali and Borthwick hump. The water free surface at the centre of the basin is able to fall further than for the corresponding bed hump case, before interacting with the bed; localised sloshing of circular waves develops above the bed trough; the slosh behaviour away from the basin centre is similar to that in the corresponding basin with a flat bed, with modes at half basin wavelength dominating.

Conclusions
This study has presented level I GN equations for shallow flow over uniform and non-uniform bed topography in the context of slosh motions in a container. It has been demonstrated that level I GN equations can represent sloshing in a closed square basin resulting from initial sinusoidal and initial Gaussian free surface perturbations. Satisfactory agreement was obtained between the model predictions and the linear analytical solution for relatively small initial wave amplitude (a ≤ 0·005 m). At larger amplitudes of initial disturbance, the numerically predicted free surface elevation time history steepened up and eventually began to develop a saw-tooth profile. Non-linear effects were particularly noticeable in the even harmonic slosh components. For a = 0·045 m and b = 2 m −2 satisfactory agreement was also obtained between the numerical predictions and Engineering and Computational Mechanics Volume 170 Issue EM2 One-dimensional and two-dimensional Green-Naghdi equations for sloshing in shallow basins Jalali and Borthwick semi-analytical solution of the early stages of free surface motions in a square, flat-bottomed basin after the initial release of the Gaussian hump. Discrepancies later evolved partly due to non-linear wave interaction effects, which were not described by the analytical theory. It was found that the non-uniform bathymetry has a localised effect on slosh motions.