One-dimensional and two-dimensional Green–Naghdi equations for sloshing in shallow basins

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.

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₀ + ζ. Here h₀ 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

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

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

In which

and

u₁(x, y, t) = v₁(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₁ and w₀ are obtained

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₀, v₀) are the horizontal velocity components at a particular point (x,y) and time t.

The general momentum conservation equation expanded in the x-direction is as follows

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

where $Pn=∂∫zbη⁡Pλndz/∂x$⁠; $P^$ is pressure at the free surface (here $P^=0$⁠) and $P¯$ is the pressure at the bottom. Applying Equation 4 in Equation 10

where $ymn=∫zbηλmλndz,ymrn=∫zbηλmλrλndz$ and $ymrn=$$∫zbηλmλrλ′ndz$ (for more details see chapter 2 of the PhD thesis by Jalali (2016) and chapter 4 of the PhD thesis by Haniffah (2013)).

Applying Equation 4 in Equation 1 and using the Krylov–Kantorovich method for the third term of the obtained equation result in

Here, the index n changes to r in Equation 4.

By summing over m and then depth integrating, Equation 12 becomes

By inserting Equation 13 into Equation 11 and implementing the chain rule for the second and third terms of Equation 11, the constrained x-direction momentum equation of level I GN equation is derived

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

where $P′n=∫zbηPλn′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₀, u₁), (v₀, v₁), (w₀, w₁)}. In z-momentum GN Equation 15, w₀ and w₁ 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$ and $P′1=P0$⁠. For eliminating the effect of P₀ 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 = ∂²z_b/∂x∂t = ∂²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

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

The corresponding 1D level I GN continuity equation is

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 (∂²u₀/∂x∂t and ∂³u₀/∂x²∂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

The 2D level I GN x-direction momentum Equation 17 is rearranged and solved using an implicit scheme. To handle Equation 17, let

In second-order central differences this becomes

where $u^0ijt=∂u0/∂t|ijt$⁠. 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

The above sets of discretised equations form the tridiagonal matrix system. The unknown values $u^0ijt$ 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.

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 anti-symmetry imposed for horizontal velocity (u) on y-direction and symmetry on x-direction.

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

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

The velocity boundary conditions for solver of 1D GN equation are

In 1D GN numerical solver, the symmetry boundaries imposed for surface elevation are

For more details, see chapter 3 of the PhD thesis by Jalali (2016) and chapter 4 of the PhD thesis by Haniffah (2013).

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. To this end, 3D visualisations and contour maps of the free surface elevation patterns in the basin were obtained on increasing grid size with Δx = Δy = 0·15 m (coarse grids), Δx = Δy = 0·0375 m (medium grids) and Δx = Δy = 7·5 × 10⁻³ m (fine grids). The medium grid size Δx = Δy = 0·0375 m is sufficient for convergence. Therefore, Δx = Δy = 0·0375 m were chosen for numerical simulations. The numerical predictions of free surface sloshing of an initial Gaussian hump for different time steps (Δt = 0·05, 0·1 and 0·2 s) have shown that Δt = 0·05 s was sufficient for producing an accurate simulation (see chapter 5 of the PhD thesis by Jalali (2016)).

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₀, 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.

Figures 1 and 2 are obtained by applying the numerical solver of 1D GN equations. This program was specifically developed 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₀ 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(a) shows satisfactory agreement between the numerical prediction (cross symbols) and the analytical solution (solid line) of the crest-induced free surface elevation time histories at the centre of the tank by applying a = 0·005 m. Here, ζ_c/a = −1 at time t = 0 s for the crest-induced free surface elevation of time history at the centre of the tank. 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)).

The numerical solver of the 2D GN equations is now verified for non-linear free surface sloshing motions arising 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

Here, ζ is the free surface elevation above still water level, where

In Equation 30, $i=−1$⁠; 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 function and