Analysis of two-dimensional nonlinear sloshing in a rectangular tank by using a concentrated mass model

Abstract

Major problems can occur when liquid sloshes in a tank, such as in liquid storage tanks during an earthquake, and this is an important engineering problem to address. To analyze this phenomenon, the finite-element method is generally used but involves many degrees of freedom when the tank is large. In this paper, a nonlinear numerical model with relatively few degrees of freedom is established for vertical and horizontal two-dimensional nonlinear sloshing in a rectangular tank excited horizontally. In addition, a method is proposed for reducing the number of degrees of freedom in the two-dimensional model. The natural frequencies, modes, and frequency responses are then compared among the concentrated mass model, theoretical calculations, and experimental results. Good agreement was achieved among them, thus demonstrating the validity of the model.

Introduction

Major problems can occur when liquid sloshes in a tank. For instance, the 2003 Tokachi-oki earthquake in Japan damaged a large tank storing oil at an industrial facility, causing a fire (Hatayama, 2008). Numerical simulations have been attempted by many researchers to understand sloshing behavior.

In analytical research, Moiseev (1958) presented a generalized theory for determining nonlinear sloshing behavior in a tank of arbitrary shape by using a perturbation method, and Faltinsen (1974) applied Moiseev’s theory to sloshing phenomena in a rectangular tank. Hayama et al. (1983) derived a similar theory that considers nonlinear free-surface conditions not considered by Moiseev, and they compared numerical results from their theory with the results of experiments involving a rectangular tank. With regard to the multimodal method, Faltinsen et al. (2000) and Faltinsen and Timokha (2002) derived a multimodal method based on two-dimensional potential theory and used it to simulate sloshing phenomena in a rectangular tank at shallow and intermediate water depths. However, the aforementioned perturbation and multimodal methods, which involve closed-form solutions, cannot be applied to tanks with complex shapes because only rectangular and cylindrical tanks yield linear theoretical solutions (Komatsu, 1958). Komatsu (1958) used the multimodal method to derive a nonlinear ordinary differential equation for sloshing and then solved it using the perturbation method for application to tanks of arbitrary shape. Love and Tait (2011, 2015) used the multimodal method to derive a nonlinear ordinary differential equation, and they solved it by numerical integration for application to flat-bottom tanks with arbitrary three-dimensional geometries. However, in the methods used by Komatsu (1958) and Love and Tait (2011, 2015), the mode shapes must be obtained in advance by using an approach such as the finite-element method (Love and Tait, 2011).

Regarding numerical analysis methods based on potential flow theory, Ikegawa (2002) derived a functional considering nonlinear free-surface boundary conditions and applied it to a finite-element method to solve nonlinear sloshing problems in a two-dimensional rectangular tank excited horizontally. Nakayama and Washizu (1980) extended Ikegawa’s method and analyzed the problem of a rectangular tank subjected to forced pitching oscillation. Wu et al. (1998) developed a finite-element model by using the Galerkin method to analyze sloshing in three-dimensional rectangular tanks. Nagashima (2009) applied the level-set extended finite-element method to sloshing in a rigid tank; in that approach, a finite-element mesh covers the entire domain in the tank including the air domain, and a signed distance function is defined to represent the distance from the free surface. In another approach, Faltinsen (1978) and Nakatama and Washizu (1981) developed a boundary-element method based on potential flow theory to solve nonlinear sloshing problems in a two-dimensional rectangular tank. However, the finite-element method and the boundary-element method involves many degrees of freedom when the tank is large.

To avoid the problem of high computational costs for FEM, an equivalent mechanical model consisting of masses, springs (or pendulums), and dashpots is often used. The model is derived so that natural frequencies and force acting on the side wall correspond to other solution. Linear mechanical models have been introduced for sloshing in a rectangular tank (Graham and Rodriguez, 1952) and in a cylindrical tank (Abramson et al., 1961) by using theoretical solutions of potential theory. Bauer (1966) developed nonlinear mechanical models for a cylindrical tank by using theoretical equations. However, there is no theoretical solution for tanks other than rectangular and cylindrical ones. For other shapes of tanks, Housner (1957) proposed a method for generating an equivalent mechanical model of an arbitrarily shaped tank by dividing elements in the horizontal and vertical directions. However, this model is linear and cannot calculate variations in the water level accurately (Kaneko, 2016). Methods to determine the linear mechanical model parameters for arbitrary tank shapes were developed by using experimental data (Sumner, 1965) or FEM results (Dorosos et al., 2008). Godderidge et al. (2012) developed a nonlinear mechanical model by using computational fluid dynamics simulation. However, for mechanical models other than cylindrical and rectangular tanks, it is necessary to perform experiments and numerical calculations in advance.

The ultimate aim of the work reported here is to establish a practical analytical model comprising one-dimensional masses and springs for analyzing nonlinear sloshing phenomena in a tank of arbitrary shape without the need to perform experiments and numerical calculations in advance. In this paper, we propose a model of a rectangular tank as a first step. In the process of deriving the model, it is not necessary to perform numerical calculations in advance to obtain the mode shapes or any other information for the mechanical model; however, as reported herein, experiments and numerical calculations with the proposed model must be performed in advance to identify the damping parameter. In previous studies, we proposed a one-dimensional concentrated mass model to analyze large-amplitude nonlinear sloshing phenomena in a rectangular tank containing shallow liquid (Ishikawa et al., 2016), where shallow is defined as a water depth that is less than one-fifth of the wavelength.

Herein, we propose a two-dimensional concentrated mass model to analyze the nonlinear sloshing phenomena to enable analysis of not only nonlinear shallow-water wave conditions but also nonlinear deep-water wave conditions. In this model, a method is proposed to transform the two-dimensional model into a one-dimensional model. As a first step toward arbitrarily shaped models, we propose the model for a rectangular tank. To validate the two-dimensional and reduced-DOF models, we compare the linear natural frequencies of the concentrated mass model with the corresponding theoretical values for a rectangular tank. To validate the model in a nonlinear region, we report on an experiment involving a rectangular tank and compare its results with the numerical ones.

Concentrated mass model

We consider two-dimensional nonlinear sloshing of liquid in a rigid rectangular tank excited horizontally, as shown in Figure 1(a). We assume that the liquid is incompressible with no surface tension and oscillates at amplitudes that are insufficient to cause wave breaking. The coordinate system is as shown in Figure 1. The liquid has horizontal length L _x and equilibrium depth L _z and is divided into N _x and N _z uniform elements in the x and z directions, respectively (Figures 1(b) and 2(a)); the element lengths are l _x and l _z in the x and z directions, respectively. The width of the tank in the y direction is taken as the unit length. Mass is concentrated on the centers of the right and left sides of each element, and the mass points are indexed as (i, j) _x . Considering the mass of liquid in the shaded area in Figure 2(b), the mass m at each mass point is

$m=\rho l_x{l}_z$

(1)

where ρ is the density of the liquid. The mass points can move in only the x direction, and the displacement of mass point (i, j) _x is δx_{i, j}, as shown in Figure 2(a). Nodal points that can move vertically are set at the top and bottom sides of each element called z nodal points, and these are indexed as (i, j) _z and have displacement δz_{i, j}. At the liquid surface, we have i = 0.

Figure 1.

Sloshing phenomena in a rectangular tank. (a) Liquid in a rectangular tank. (b) elements and mass points.

Figure 2.

Elements, mass points in x direction, and z nodal points. (a) Displacements of mass points and z nodal points. (b) Number of elements, mass points, and z nodal points.

Connecting nonlinear springs

The nonlinear restoring forces acting on the mass points are derived from the hydrostatic pressure caused by the vertical displacement of the liquid surface and the hydrodynamic pressure caused by the vertical acceleration of the liquid. First, we derive the relationship between the displacement δx_{i, j} of the mass points and the displacement δz_{i, j} of the z nodal points. Because the liquid is incompressible, we assume that the z nodal points move while keeping the volume of each element in an arbitrary column j constant when the masses move, as shown Figure 3. From these relations, we obtain

$\left(\delta z_{i-1,j}-\delta z_{i,j}\right)\left(l_x+\delta x_{i,j}-\delta x_{i,j-1}\right)=-l_z\left(\delta x_{i,j}-\delta x_{i,j-1}\right)$

(2)

the left-hand side of which is the extruded volume in the z direction and the right-hand side of which is the volume change by the displacements δx_{i, j} and δx_{i, j−1}. From equation (2), the relative displacement of an arbitrary z nodal point is

$\delta z_{i-1,j}-\delta z_{i,j}=-\frac{l_z\left(\delta x_{i,j}-\delta x_{i,j-1}\right)}{l_x+\delta x_{i,j}-\delta x_{i,j-1}}$

(3)

As the boundary condition at the bottom of the tank in Figure 4, the displacement in the z direction at the bottom is assumed as

$\delta z_{N_z,j}=0$

(4)

Assuming that the displacements in the z direction are stacked from the bottom, from equations (3) and (4), the displacement at z nodal point (i, j) _z is

$\delta z_{i,j}=-{\displaystyle \sum _{k=i+1}^{N_z}\frac{l_z\left(\delta x_{k,j}-\delta x_{k,j-1}\right)}{l_x+\delta x_{k,j}-\delta x_{k,j-1}}}\hspace{1em}\hspace{0.17em}\left(0\le i\le N_z-1\right)$

(5)

However, upon using equation (5) and considering the elements in column j + 1, overlaps or gaps occur as shown in Figure 5, and these change the volume of the whole liquid despite it being incompressible. To avoid this problem, we introduce $\delta {\overline{x}}_j$ as the average displacement of the mass points in column j, as shown in Figure 4. The average displacement is given by

$\delta {\overline{x}}_j=\frac{1}{N_z}{\displaystyle \sum _{i=1}^{N_z}\delta x_{i,j}}$

(6)

We use $\delta {\overline{x}}_j$ when considering the extruded volume in the z direction. Replacing $\delta x_{i,j}-\delta x_{i,j-1}$ on the left-hand side in equation 2) with $\delta {\overline{x}}_j-\delta {\overline{x}}_{j-1}$ , the relative displacement of an arbitrary z nodal point becomes

$\delta z_{i-1,j}-\delta z_{i,j}=-\frac{l_z\left(\delta x_{i,j}-\delta x_{i,j-1}\right)}{l_x+\delta {\overline{x}}_j-\delta {\overline{x}}_{j-1}}$

(7)

and the displacement at z nodal point (i, j) _z becomes

$\delta z_{i,j}=-\frac{l_z}{l_x+\delta {\overline{x}}_j-\delta {\overline{x}}_{j-1}}{\displaystyle \sum _{k=i+1}^{N_z}\left(\delta x_{k,j}-\delta x_{k,j-1}\right)}\hspace{0.17em}\left(0\le i\le N_z-1\right)$

(8)

Considering equation 6) and (8), the displacement δz_{0, j} of the z nodal points at the top matches the displacement η _j of the liquid surface

${\eta }_j=\delta z_{0,j}=-\frac{l_z{N}_z\left(\delta {\overline{x}}_j-\delta {\overline{x}}_{j-1}\right)}{l_x+\delta {\overline{x}}_j-\delta {\overline{x}}_{j-1}}=-\frac{L_z\left(\delta {\overline{x}}_j-\delta {\overline{x}}_{j-1}\right)}{l_x+\delta {\overline{x}}_j-\delta {\overline{x}}_{j-1}}$

(9)

Equations (8) and (9) are nonlinear relationships, which is caused by nonlinearity of the element deformation.

Figure 3.

Deformation of element [i, j]. (a) Before deformation. (b) After deformation.

Figure 4.

Deformation of element from the bottom.

Figure 5.

Overlap and gap among elements.

Second, we derive the pressure p_{i, j} in element [i, j]. In the area of column j in Figure 6(a), Euler’s equation in the vertical direction is written as

$\frac{\partial w_j\left(z,t\right)}{\partial t}=-\frac{1}{\rho }\frac{\partial p_j\left(z,t\right)}{\partial z}-g$

(10)

where w _j (z, t) is the z-directional velocity at the z coordinate in column j, p _j (z, t) is the pressure in column j, and g is the acceleration due to gravity. The z coordinate of the center position in element [i, j] is z _i = −(i – 1/2)l _z ( $i\ge 2$ ). Integrating equation (10) with respect to z from z _i to η _j , as shown in Figure 6(b), the pressure p_{i, j} at z _i becomes

$\begin{array}{ll}{p}_{1,j}=\rho {\displaystyle {\int }_{({\eta }_j-l_z)/2}^{{\eta }_j}\frac{\partial w_j}{\partial t}dz}+\frac{1}{2}\rho g\left({\eta }_j+l_z\right)& (i=1)\\ p_{i,j}=\rho {\displaystyle {\int }_{z_i}^{{\eta }_j}\frac{\partial w_j}{\partial t}dz}+\rho g\left\{{\eta }_j+\left(i-1\right)l_z+\frac{l_z}{2}\right\}& (i\ge 2)\end{array}\}$

(11)

Herein, the pressure at the liquid surface is zero. We assume that the pressure p_{i, j} of element [i, j] is represented by the pressure at the center position z _i . The top element [1, j] contains the region protruding above the equilibrium free surface when the liquid level changes, so the z coordinate of the center position of the top element is z₁ = (η _j − l _z )/2, as shown in Figure 7(a). On the right-hand sides in equation (11), the first terms are related to the hydrodynamic pressure and the second terms are related to the hydrostatic pressure.

Figure 6.

Control volume and interval of integration. (a) Control volume at column j. (b) Interval of integration.

Figure 7.

Acceleration of each element. (a) Around element [1, j]. (b) Around element [i, j].

Third, we derive the restoring force acting on mass point (i, j) _x from pressure p_{i, j}. Assuming that $\delta {\ddot{\stackrel{⌢}{z}}}_{i,j}$ is the average vertical acceleration of z nodal points ${(i-1,j)}_z$ and ${(i,j)}_z$ (Figure 7), $\delta {\ddot{\stackrel{⌢}{z}}}_{i,j}$ is

$\delta {\ddot{\stackrel{⌢}{z}}}_{i,j}=\frac{\delta {\ddot{z}}_{i-1,j}+\delta {\ddot{z}}_{i,j}}{2}$

(12)

We assume that the vertical acceleration in element [i, j] is represented by the acceleration at the center position $\delta {\ddot{\stackrel{⌢}{z}}}_{i,j}$ . In the integration in equation (11), $\partial w_j/\partial t$ is replaced by the vertical average of vertical accelerations $\delta {\ddot{\stackrel{⌢}{z}}}_{i,j}$ . The integration in equation (11) is replaced by the summation using ${\ddot{\stackrel{⌢}{z}}}_{i,j}$ from element [1, j] to element [i, j] as

$\begin{array}{ll}{f}_{1,j}=\frac{\rho }{2}\delta {\ddot{\stackrel{⌢}{z}}}_{1,j}{\left({\eta }_j+l_z\right)}^2+\frac{\rho }{2}g\left({\eta }_j+l_z\right){\hspace{0.17em}}^2& (i=1)\\ f_{i,j}=\rho \delta {\ddot{\stackrel{⌢}{z}}}_{1,j}{l}_z\left({\eta }_j+l_z\right)+\rho {\displaystyle \sum _{k=2}^{i-1}\delta {\ddot{\stackrel{⌢}{z}}}_{k,j}{l}_z^2}+\rho \delta {\ddot{\stackrel{⌢}{z}}}_{i,j}\frac{l_z^2}{2}+\rho g{l}_z\left\{{\eta }_j+\left(i-1\right)l_z+\frac{l_z}{2}\right\}& (2\le i\le N_z)\end{array}\}$

(14)

Multiplying p_{i, j} in equation (11) by the vertical element length l _z and using equation (13), the force f_{i, j} acting on mass point (i, j) _x becomes

$\begin{array}{ll}{\beta }_{1,j}=\frac{1}{2}\rho \delta {\ddot{\stackrel{⌢}{z}}}_{1,j}\left({\eta }_j+l_z\right){\hspace{0.17em}}^2& (i=1)\\ {\beta }_{i,j}=\rho \left\{\delta {\ddot{\stackrel{⌢}{z}}}_{1,j}{l}_z\left({\eta }_j+l_z\right)+{\displaystyle \sum _{k=2}^{i-1}\delta {\ddot{\stackrel{⌢}{z}}}_{k,j}{l}_z^2}+\frac{1}{2}\delta {\ddot{\stackrel{⌢}{z}}}_{i,j}{l}_z^2\right\}& (2\le i\le N_z)\end{array}\}$

(18)

Herein, we assume that the vertical length of the top elements (i = 1) is (η _j + l _z ) when the liquid level changes. In the equilibrium state, the vertical displacement and the acceleration of the liquid surface are zero; therefore, the force $f_{i,j}^0$ acting on an arbitrary mass point in the equilibrium state is

$f_{i,j}^0=\rho g{l}_z^2\left(i-\frac{1}{2}\right)$

(15)

Therefore, the restoring force $\delta f_{i,j}^k$ from the pressure is written as

$\delta f_{i,j}^k=f_{i,j}-f_{i,j}^0={\alpha }_{i,j}+{\beta }_{i,j}$

(16)

$\begin{array}{ll}{\alpha }_{1,j}=\rho g{l}_z{\eta }_j+\frac{1}{2}\rho g{\eta }_j^2& (i=1)\\ {\alpha }_{i,j}=\rho g{l}_z{\eta }_j& (2\le i\le N_z)\end{array}\}$

(17)

$\begin{array}{ll}{\beta }_{1,j}^{\text{L}}=\frac{1}{2}\rho \delta {\ddot{\stackrel{⌢}{z}}}_{1,j}{l}_z^2& (i=1)\\ {\beta }_{i,j}^{\text{L}}=\rho \left(\delta {\ddot{\stackrel{⌢}{z}}}_{1,j}{l}_z^2+{\displaystyle \sum _{k=1}^{i-1}\delta {\ddot{\stackrel{⌢}{z}}}_{k,j}{l}_z^2}+\frac{1}{2}\delta {\ddot{\stackrel{⌢}{z}}}_{i,j}{l}_z^2\right)& (i\ge 2)\end{array}\}$

(22)

where α_{i, j} is the restoring force derived from the static pressure caused by the vertical displacement of the liquid surface, and β_{i, j} is the restoring force derived from the dynamic pressure that is caused by vertical acceleration of the liquid.

Linearizing equations (7), (9), (17), and (18), we obtain

$z_{i,j}^{\text{L}}=-\frac{l_z}{l_x}{\displaystyle \sum _{k=i+1}^{N_z}\left(\delta x_{k,j}-\delta x_{k,j-1}\right)}$

(19)

${\eta }_j^{\text{L}}=\delta z_{0,j}=-\frac{l_z}{l_x}{\displaystyle \sum _{k=1}^{N_z}\left(\delta x_{k,j}-\delta x_{k,j-1}\right)}$

(20)

${\alpha }_{i,j}^{\text{L}}=\rho g{l}_z{\eta }_j^{\text{L}}$

(21)

$\begin{array}{l}\mathbf{M}\ddot{\mathbf{x}}+\mathbf{C}\dot{\mathbf{x}}+\mathbf{K}\mathbf{x}=\mathbf{G}\\ \mathbf{G}=\frac{1}{2}\rho g{\left[\begin{array}{cccc}{\mathbf{\sigma }}^T& \mathbf{0}& \cdots & \mathbf{0}\end{array}\right]}^T\\ \mathbf{\sigma }={\left[\begin{array}{cccc}{\eta }_1^2-{\eta }_2^2& {\eta }_2^2-{\eta }_3^2& \cdots & {\eta }_{N-1}^2-{\eta }_N^2\end{array}\right]}^T\end{array}\}$

(24)

Note that the average displacement $\delta {\overline{x}}_j$ in equation (6) is not necessary in the linearized equations (19) and (20).

Equations of motion

The equation of motion of mass point (i, j) _x is expressed as follows by considering the restoring forces acting from elements [i, j] and [i, j+1] (Figure 2(b))

$m\hspace{0.17em}\delta {\ddot{x}}_{i,j}={\alpha }_{i,j}+{\beta }_{i,j}-{\alpha }_{i,j+1}-{\beta }_{i,j+1}-f_{i,j}^c+f_{i,j+1}^c+f_{i,j}^e\left(t\right)$

(23)

Here, $f_{i,j}^c$ is the linear viscous damping force, and $f_{i,j}^e\left(t\right)$ is an external force acting on mass point (i, j) _x . Equation (23) becomes an equation involving only displacements δx_{i, j} by substituting equation (9) into equations (17) and (18), and equations (8) into equations (18). The mass points at the wall of the container are assumed to be fixed, so the displacements δx_{i, 0} and $\delta x_{i,\hspace{0.17em}{N}_x}$ are zero. Writing the equations of motion for all the mass points and arranging them in matrix form, we obtain

${\mathbf{M}}_{i,j}^a=\left[\begin{array}{cccc}{m}_{j,1}^{z1}{m}_{i,1}^{z2}+m_{j,2}^{z1}{m}_{i,2}^{z2}& -m_{j,2}^{z1}{m}_{i,2}^{z2}& & \\ -m_{j,2}^{z1}{m}_{i,2}^{z2}& m_{j,2}^{z1}{m}_{i,2}^{z2}+m_{j,3}^{z1}{m}_{i,3}^{z2}& -m_{j,3}^{z1}{m}_{i,3}^{z2}& \\ & \ddots & \ddots & \ddots \\ & & -m_{j,N_x-1}^{z1}{m}_{i,N_x-1}^{z2}& m_{j,N_x-1}^{z1}{m}_{i,N_x-1}^{z2}+m_{j,N_x}^{z1}{m}_{i,N_x}^{z2}\end{array}\right]\hspace{0.17em}$

(29)

where M, C, and K are the mass, damping, and stiffness matrices, respectively, of size $N_z\left(N_x-1\right)\times N_z\left(N_x-1\right)$ . Pertaining to the nonlinearity of the displacement of the liquid surface, matrix G is derived from the second term of the first equation of equations (17) and is of size $N_z\left(N_x-1\right)\times 1$ . If x _i is a displacement vector containing the displacements of the mass points in line i and x is a vector containing each x _i , then x and x _i are expressed as

$\begin{array}{l}\mathbf{x}={\left[\begin{array}{cccc}{\mathbf{x}}_1^T& {\mathbf{x}}_2^T& \cdots & {\mathbf{x}}_{N_z}^T\end{array}\right]}^{\text{T}}\\ {\mathbf{x}}_i={\left[\begin{array}{cccc}\delta x_{i,1}& \delta x_{i,2}& \cdots & \delta x_{i,N_x-1}\end{array}\right]}^{\text{T}}\end{array}\}$

(25)

Because α_{i, j} in equation (17) is a function of the displacements of the mass points and β_{i, j} in equation (18) includes their accelerations, the former constitutes the stiffness matrix and the latter the mass matrix. The components of the mass matrix are written as follows

$\mathbf{M}={\mathbf{M}}^m+{\mathbf{M}}^d=\left[\begin{array}{cccc}{\mathbf{M}}_{11}& {\mathbf{M}}_{12}& \cdots & {\mathbf{M}}_{1N_z}\\ {\mathbf{M}}_{21}& {\mathbf{M}}_{22}& \cdots & {\mathbf{M}}_{2N_z}\\ ⋮& ⋮& \ddots & ⋮\\ {\mathbf{M}}_{N_{z}1}& {\mathbf{M}}_{N_{z}2}& \cdots & {\mathbf{M}}_{N_z{N}_z}\end{array}\right]$

(26)

where ${\mathbf{M}}^m$ and ${\mathbf{M}}^d$ are

${\mathbf{M}}^m=\left[\begin{array}{cccc}m& 0& \cdots & 0\\ 0& m& \ddots & ⋮\\ ⋮& \ddots & \ddots & 0\\ 0& \cdots & 0& m\end{array}\right]\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}{\mathbf{M}}^d=\left[\begin{array}{ccccc}{\mathbf{M}}_1^1/2& {\mathbf{M}}_2^1& {\mathbf{M}}_3^1& \cdots & {\mathbf{M}}_{N_z}^1\\ {\mathbf{M}}_{\mathrm{2,1}}^b& {\mathbf{M}}_{\mathrm{2,2}}^a& {\mathbf{M}}_{\mathrm{2,3}}^b& \cdots & {\mathbf{M}}_{2,N_z}^b\\ {\mathbf{M}}_{\mathrm{3,1}}^b& {\mathbf{M}}_{\mathrm{3,2}}^b& {\mathbf{M}}_{\mathrm{3,3}}^a& \cdots & {\mathbf{M}}_{3,N_z}^b\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ {\mathbf{M}}_{N_z,1}^b& {\mathbf{M}}_{N_z,2}^b& {\mathbf{M}}_{N_z,3}^b& \cdots & {\mathbf{M}}_{N_z,N_z}^a\end{array}\right]$

(27)

${\mathbf{M}}^m$ and ${\mathbf{M}}^d$ represent the mass of liquid and the dynamic pressure, respectively. ${\mathbf{M}}_i^1$ , which means the dynamic pressure acting on the masses in the first row, is written as

$\begin{array}{l}{\mathbf{M}}_i^1=\left[\begin{array}{cccc}{m}_{i,1}^{z1}+m_{i,2}^{z1}& -m_{i,2}^{z1}& & \\ m_{i,2}^{z1}& m_{i,2}^{z1}+m_{i,3}^{z1}& -m_{i,3}^{z1}& \\ & \mathrm{\ddots }& \mathrm{\ddots }& \mathrm{\ddots }\\ & & -m_{i,N_x-1}^{z1}& m_{i,N_x-1}^{z1}+m_{i,N_x}^{z1}\end{array}\right]\hspace{0.17em}\hspace{0.17em}\\ \hspace{0.17em}{m}_{i,j}^{z1}=\frac{\rho {\left(l_z+{\eta }_j\right)}^2{l}_x{l}_z}{2{\left(l_x+\delta x_{i,j}-\delta x_{i,j-1}\right)}^2}\end{array}\}.$

(28)

${\mathbf{M}}_{i,j}^a$ and ${\mathbf{M}}_{i,j}^b$ , which means the dynamic pressure acting on the masses in the rows other than the first one, are written as

$\begin{array}{ll}{\displaystyle {\int }_{({\eta }_j-l_z)/2}^{{\eta }_j}\frac{\partial w_j}{\partial t}dz}=\frac{\delta {\ddot{\stackrel{⌢}{z}}}_{1,j}\left({\eta }_j+l_z\right)}{2}\hspace{0.17em}& (i=1)\\ {\displaystyle {\int }_{z_i}^{{\eta }_j}\frac{\partial w_j}{\partial t}dz}=\delta {\ddot{\stackrel{⌢}{z}}}_{1,j}\left(l_z+{\eta }_j\right)+{\displaystyle \sum _{k=2}^{i-1}}\delta {\ddot{\stackrel{⌢}{z}}}_{k,j}{l}_z+\delta {\ddot{\stackrel{⌢}{z}}}_{i,j}\frac{l_z}{2}& (2\le i<N_z)\end{array}\}$

(13)

${\mathbf{M}}_{i,j}^b=\left[\begin{array}{cccc}{m}_{j,1}^{z1}{m}_{i,1}^{z3}+m_{j,2}^{z1}{m}_{i,2}^{z3}& -m_{j,2}^{z1}{m}_{i,2}^{z3}& & \\ -m_{j,2}^{z1}{m}_{i,2}^{z3}& m_{j,2}^{z1}{m}_{i,2}^{z3}+m_{j,3}^{z1}{m}_{i,3}^{z3}& -m_{j,3}^{z1}{m}_{i,3}^{z3}& \\ & \ddots & \ddots & \ddots \\ & & -m_{j,N_x-1}^{z1}{m}_{i,N_x-1}^{z3}& m_{j,N_x-1}^{z1}{m}_{i,N_x-1}^{z3}+m_{j,N_x}^{z1}{m}_{i,N_x}^{z3}\end{array}\right]$

(30)

where

$m_{i,k}^{z2}=\left(l_z+{\eta }_k\right)l_z+\left(i-\frac{7}{4}\right)l_z^2\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}{m}_{i,k}^{z3}=\left(l_z+{\eta }_k\right)l_z+\left(i-\frac{3}{2}\right)l_z^2$

(31)

From equations (26) and (27), M_{i, j} becomes

$\begin{array}{l}{\mathbf{M}}_{11}=\text{diag}[\begin{array}{cccc}m& m& \cdots & m\end{array}\hspace{0.17em}]+\frac{{\mathbf{M}}_1^1}{2}\\ {\mathbf{M}}_{1j}={\mathbf{M}}_j^1\hspace{1em}(j\ge 2)\\ {\mathbf{M}}_{ii}=\text{diag}[\begin{array}{cccc}m& m& \cdots & m\end{array}\hspace{0.17em}]+{\mathbf{M}}_{ii}^a\\ {\mathbf{M}}_{ij}={\mathbf{M}}_{ij}^b\hspace{1em}(i\ne j)\end{array}\}$

(32)

Regarding the linearized mass matrix M^L, $m_{i,j}^{z1}$ in equation (28) and $m_{i,k}^{z2}$ and $m_{i,k}^{z3}$ in equation (31) are treated as

$m_{i,j}^{z1,\text{L}}=\frac{\rho l_z^3}{l_x},\hspace{0.17em}\hspace{0.17em}{m}_{i,k}^{z2,\text{L}}=l_z^2+\left(i-\frac{7}{4}\right)l_z^2\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}{m}_{i,k}^{z3,\text{L}}=l_z^2+\left(i-\frac{3}{2}\right)l_z^2$

(33)

The components of the stiffness matrix are written as follows:

$\begin{array}{l}\mathbf{K}=\left[\begin{array}{cccc}\widehat{\mathbf{K}}& \widehat{\mathbf{K}}& \cdots & \widehat{\mathbf{K}}\\ \widehat{\mathbf{K}}& \widehat{\mathbf{K}}& \cdots & \widehat{\mathbf{K}}\\ ⋮& ⋮& \ddots & ⋮\\ \widehat{\mathbf{K}}& \widehat{\mathbf{K}}& \cdots & \widehat{\mathbf{K}}\end{array}\right]\\ \widehat{\mathbf{K}}=\left[\begin{array}{ccccc}{k}_1+k_2& -k_2& 0& \cdots & 0\\ -k_2& k_2+k_3& -k_3& \ddots & ⋮\\ 0& \ddots & \ddots & \ddots & 0\\ ⋮& \ddots & -k_{N-2}& k_{N-2}+k_{N-1}& -k_{N-1}\\ 0& \cdots & 0& -k_{N-1}& k_{N-1}+k_N\end{array}\right],\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}{k}_j=\frac{\rho g{l}_z^2}{l_x+\delta {\overline{x}}_j-\delta {\overline{x}}_{j-1}}\end{array}\}$

(34)

From equations (28), (31), and (34), the mass and stiffness matrices contain displacements.

where

In addition, the damping matrix C is represented as a proportion of the linearized stiffness matrix K^L

$\begin{array}{l}\mathbf{C}=\gamma {\mathbf{K}}^L\\ {\mathbf{K}}^L=\left[\begin{array}{cccc}\tilde{\mathbf{K}}& \tilde{\mathbf{K}}& \cdots & \tilde{\mathbf{K}}\\ \tilde{\mathbf{K}}& \tilde{\mathbf{K}}& \cdots & \tilde{\mathbf{K}}\\ ⋮& ⋮& \ddots & ⋮\\ \tilde{\mathbf{K}}& \tilde{\mathbf{K}}& \cdots & \tilde{\mathbf{K}}\end{array}\right],\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\tilde{\mathbf{K}}=\left[\begin{array}{ccccc}2\tilde{k}& -\tilde{k}& 0& \cdots & 0\\ -\tilde{k}& 2\tilde{k}& -\tilde{k}& \ddots & M\\ 0& \ddots & \ddots & \ddots & 0\\ M& \ddots & -\tilde{k}& 2\tilde{k}& -\tilde{k}\\ 0& \cdots & 0& -\tilde{k}& 2\tilde{k}\end{array}\right],\hspace{0.17em}\hspace{0.17em}\tilde{k}=\frac{\rho g{l}_z^2}{l_x}\end{array}\}$

(35)

where γ is a constant of the proportional viscous damping whose optimal value is determined experimentally. The rationale for the damping matrix is unclear. Herein, we assume proportional viscous damping and adjust its value so that the numerical results agree with the experimental results. Valid modeling of sloshing damping is a future task.

Reducing the number of degrees of freedom

In the model proposed in the section Concentrated Mass Model, the elements are divided vertically to consider the vertical velocity distribution, but, consequently, the model has an unnecessarily high number of DOFs, with zero-frequency eigenvalues being generated. In this section, we propose a way to reduce the number of DOFs of the concentrated mass model, as shown in Figure 8.

Figure 8.

Reducing the number of degrees of freedom (DOFs). (a) Two-dimensional model. (b) Reduced-DOF model.

Using equations (26), (34), and (35), the equations of motion in equations (24) are rewritten as

$\left[\begin{array}{cccc}{\mathbf{M}}_{11}& {\mathbf{M}}_{12}& \cdots & {\mathbf{M}}_{1N_z}\\ {\mathbf{M}}_{21}& {\mathbf{M}}_{22}& \cdots & {\mathbf{M}}_{2N_z}\\ ⋮& ⋮& \ddots & ⋮\\ {\mathbf{M}}_{N_{z}1}& {\mathbf{M}}_{N_{z}2}& \cdots & {\mathbf{M}}_{N_z{N}_z}\end{array}\right]\left[\begin{array}{c}{\dot{\mathbf{x}}}_1\\ {\dot{\mathbf{x}}}_2\\ ⋮\\ {\dot{\mathbf{x}}}_{N_z}\end{array}\right]+\gamma \left[\begin{array}{cccc}\tilde{\mathbf{K}}& \tilde{\mathbf{K}}& \cdots & \tilde{\mathbf{K}}\\ \tilde{\mathbf{K}}& \tilde{\mathbf{K}}& \cdots & \tilde{\mathbf{K}}\\ ⋮& ⋮& \ddots & ⋮\\ \tilde{\mathbf{K}}& \tilde{\mathbf{K}}& \cdots & \tilde{\mathbf{K}}\end{array}\right]\left[\begin{array}{c}{\dot{\mathbf{x}}}_1\\ {\dot{\mathbf{x}}}_2\\ ⋮\\ {\dot{\mathbf{x}}}_{N_z}\end{array}\right]+\left[\begin{array}{cccc}\widehat{\mathbf{K}}& \widehat{\mathbf{K}}& \cdots & \widehat{\mathbf{K}}\\ \widehat{\mathbf{K}}& \widehat{\mathbf{K}}& \cdots & \widehat{\mathbf{K}}\\ ⋮& ⋮& \ddots & ⋮\\ \widehat{\mathbf{K}}& \widehat{\mathbf{K}}& \cdots & \widehat{\mathbf{K}}\end{array}\right]\left[\begin{array}{c}{\dot{\mathbf{x}}}_1\\ {\dot{\mathbf{x}}}_2\\ ⋮\\ {\dot{\mathbf{x}}}_{N_z}\end{array}\right]=\frac{1}{2}\rho g\left[\begin{array}{c}\mathbf{\sigma }\\ \mathbf{0}\\ ⋮\\ \mathbf{0}\end{array}\right]$

(36)

Here, the matrices M_{i, j}, $\tilde{\mathbf{K}}$ , and $\widehat{\mathbf{K}}$ are of size $(N_x-1)\times (N_x-1)$ . The damping matrix C includes the block matrix $\tilde{\mathbf{K}}$ overall and the stiffness matrix K includes the block matrix $\widehat{\mathbf{K}}$ . Subtracting the equation of row 1 from those of rows 2 to N _z in equation (36), we obtain

$\left[\begin{array}{cccc}{\mathbf{M}}_{11}& {\mathbf{M}}_{12}& \cdots & {\mathbf{M}}_{1N_z}\\ {\tilde{\mathbf{M}}}_{21}& {\tilde{\mathbf{M}}}_{22}& \cdots & {\tilde{\mathbf{M}}}_{2N_z}\\ ⋮& ⋮& \ddots & ⋮\\ {\tilde{\mathbf{M}}}_{N_{z}1}& {\tilde{\mathbf{M}}}_{N_{z}2}& \cdots & {\tilde{\mathbf{M}}}_{N_z{N}_z}\end{array}\right]\left[\begin{array}{c}{\dot{\mathbf{x}}}_1\\ {\dot{\mathbf{x}}}_2\\ ⋮\\ {\dot{\mathbf{x}}}_{N_z}\end{array}\right]+\gamma \left[\begin{array}{cccc}\tilde{\mathbf{K}}& \tilde{\mathbf{K}}& \cdots & \tilde{\mathbf{K}}\\ \mathbf{0}& \mathbf{0}& \cdots & \mathbf{0}\\ ⋮& ⋮& \ddots & ⋮\\ \mathbf{0}& \mathbf{0}& \cdots & \mathbf{0}\end{array}\right]\left[\begin{array}{c}{\dot{\mathbf{x}}}_1\\ {\dot{\mathbf{x}}}_2\\ ⋮\\ {\dot{\mathbf{x}}}_{N_z}\end{array}\right]+\left[\begin{array}{cccc}\widehat{\mathbf{K}}& \widehat{\mathbf{K}}& \cdots & \widehat{\mathbf{K}}\\ \mathbf{0}& \mathbf{0}& \cdots & \mathbf{0}\\ ⋮& ⋮& \ddots & ⋮\\ \mathbf{0}& \mathbf{0}& \cdots & \mathbf{0}\end{array}\right]\left[\begin{array}{c}{\dot{\mathbf{x}}}_1\\ {\dot{\mathbf{x}}}_2\\ ⋮\\ {\dot{\mathbf{x}}}_{N_z}\end{array}\right]=\frac{1}{2}\rho g\left[\begin{array}{c}\mathbf{\sigma }\\ -\mathbf{\sigma }\\ ⋮\\ -\mathbf{\sigma }\end{array}\right]$

(37)

Here, ${\tilde{\mathbf{M}}}_{i,j}$ is

${\tilde{\mathbf{M}}}_{i,j}={\mathbf{M}}_{i,j}-{\mathbf{M}}_{1,j}$

(38)

In addition, we introduce the average displacement vector $\overline{x}$ , whose elements are the average displacements in equation (6). The average displacement vector can be written as the summation of ${\mathbf{x}}_1$ to ${\mathbf{x}}_{N_z}$

$\overline{\mathbf{x}}=\left[\begin{array}{c}\delta {\overline{x}}_1\\ \delta {\overline{x}}_2\\ ⋮\\ \delta {\overline{x}}_{N_x-1}\end{array}\right]=\frac{1}{N_z}\left[\begin{array}{cccc}{\mathbf{I}}_{N_x-1}& {\mathbf{I}}_{N_x-1}& \cdots & {\mathbf{I}}_{N_x-1}\end{array}\right]\left[\begin{array}{c}{\mathbf{x}}_1\\ {\mathbf{x}}_2\\ \cdots \\ {\mathbf{x}}_{N_z}\end{array}\right]$

(39)

where I is the unit matrix. Taking the equations of rows 2 to N _z in equation (37) and inserting equation (39) into the top row, we obtain

$\left[\begin{array}{cccc}{\mathbf{I}}_{N_x-1}& {\mathbf{I}}_{N_x-1}& \cdots & {\mathbf{I}}_{N_x-1}\\ {\tilde{\mathbf{M}}}_{21}& {\tilde{\mathbf{M}}}_{22}& \cdots & {\tilde{\mathbf{M}}}_{2N_z}\\ ⋮& ⋮& \ddots & ⋮\\ {\tilde{\mathbf{M}}}_{N_{z}1}& {\tilde{\mathbf{M}}}_{N_{z}2}& \cdots & {\tilde{\mathbf{M}}}_{N_z{N}_z}\end{array}\right]\left[\begin{array}{c}{\dot{\mathbf{x}}}_1\\ {\dot{\mathbf{x}}}_2\\ ⋮\\ {\dot{\mathbf{x}}}_{N_z}\end{array}\right]=N_z\left[\begin{array}{c}{\mathbf{I}}_{N_x-1}\\ \mathbf{0}\\ ⋮\\ \mathbf{0}\end{array}\right]\ddot{\overline{\mathbf{x}}}-\frac{1}{2}\rho g\left[\begin{array}{c}\mathbf{0}\\ \mathbf{\sigma }\\ ⋮\\ \mathbf{\sigma }\end{array}\right]$

(40)

which can be written as

$\begin{array}{l}{\mathbf{A}}_1\ddot{\mathbf{x}}={\mathbf{A}}_2\ddot{\overline{\mathbf{x}}}-\tilde{\mathbf{S}}\\ {\mathbf{A}}_1=\left[\begin{array}{cccc}{\mathbf{I}}_{N_x-1}& {\mathbf{I}}_{N_x-1}& \cdots & {\mathbf{I}}_{N_x-1}\\ {\tilde{\mathbf{M}}}_{21}& {\tilde{\mathbf{M}}}_{22}& \cdots & {\tilde{\mathbf{M}}}_{2N_z}\\ ⋮& ⋮& \cdots & ⋮\\ {\tilde{\mathbf{M}}}_{N_{z}1}& {\tilde{\mathbf{M}}}_{N_{z}2}& \cdots & {\tilde{\mathbf{M}}}_{N_z{N}_z}\end{array}\right]\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}{\mathbf{A}}_2=N_z\left[\begin{array}{c}{\mathbf{I}}_{N_x-1}\\ \mathbf{0}\\ ⋮\\ \mathbf{0}\end{array}\right]\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\tilde{\mathbf{S}}=\frac{1}{2}\rho g\left[\begin{array}{c}\mathbf{0}\\ \mathbf{\sigma }\\ ⋮\\ \mathbf{\sigma }\end{array}\right]\hspace{0.17em}\hspace{0.17em}\hspace{0.17em}\end{array}\}$

(41)

Multiplying ${\mathbf{A}}_1^{-1}$ by equation (41) from the left-hand side, we obtain

$\begin{array}{l}\ddot{\mathbf{x}}=\mathbf{T}\ddot{\overline{\mathbf{x}}}-\mathbf{S}\hspace{0.17em}\\ \mathbf{T}={\mathbf{A}}_1^{-1}{\mathbf{A}}_2\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\mathbf{S}={\mathbf{A}}_1^{-1}\tilde{\mathbf{S}}\end{array}\}$

(42)

which represents the relationship between the two-dimensional acceleration $\ddot{\mathbf{x}}$ and the one-dimensional acceleration $\ddot{\overline{\mathbf{x}}}$ . T is the transformation matrix from the two-dimensional model to the one-dimensional model. We assume that the displacement vectors x and $\overline{\mathbf{x}}$ and the velocity vectors $\dot{\mathbf{x}}$ and $\dot{\overline{\mathbf{x}}}$ are also related by T

$\begin{array}{l}\mathbf{x}=\mathbf{T}\overline{\mathbf{x}}\\ \dot{\mathbf{x}}=\mathbf{T}\dot{\overline{\mathbf{x}}}\end{array}\}$

(43)

Substituting equation (42) and (43) into equation (24) and left-multiplying by ${\mathbf{T}}^{\text{T}}$ , we obtain

$\begin{array}{l}\overline{\mathbf{M}}\ddot{\overline{\mathbf{x}}}+\overline{\mathbf{C}}\dot{\overline{\mathbf{x}}}+\overline{\mathbf{K}}\overline{\mathbf{x}}=\overline{\mathbf{G}}+{\mathbf{T}}^T\mathbf{M}\mathbf{S}\\ \overline{\mathbf{M}}={\mathbf{T}}^T\mathbf{M}\mathbf{T}\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\overline{\mathbf{C}}={\mathbf{T}}^T\mathbf{C}\mathbf{T}\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\overline{\mathbf{K}}={\mathbf{T}}^T\mathbf{K}\mathbf{T}\hspace{0.17em}\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}\overline{\mathbf{G}}={\mathbf{T}}^T\mathbf{G}\end{array}\}$

(44)

Equation (44) is the equation of motion of the reduced-DOF model using the average displacement. In equation (44), there are $N_x-1$ DOFs, which is fewer than the $N_z\times (N_x-1)$ DOFs in equation (24). Having solved for $\overline{x}$ in equation (44), the displacements η _j of the liquid surface can be calculated by using equation (9).

However, the transformation matrix T in equation (42) and (43) contains the displacements δx_{i, j} because it contains the mass matrix ${\tilde{\mathbf{M}}}_{i,j}$ whose elements are a function of the displacements, as shown in equations (28) and (31). Therefore, T changes at each time step in the numerical calculation, as does the equation of motion equations (44). In the case of linear analysis, the transformation matrix T^L becomes

$\begin{array}{l}{\mathbf{T}}^L={\left({\mathbf{A}}_1^L\right)}^{-1}{\mathbf{A}}_2\\ {\mathbf{A}}_1^L=\left[\begin{array}{cccc}{\mathbf{I}}_{N_x-1}& {\mathbf{I}}_{N_x-1}& \cdots & {\mathbf{I}}_{N_x-1}\\ {\tilde{\mathbf{M}}}_{21}^L& {\tilde{\mathbf{M}}}_{22}^L& \cdots & {\tilde{\mathbf{M}}}_{2N_z}^L\\ ⋮& ⋮& \ddots & ⋮\\ {\tilde{\mathbf{M}}}_{N_{z}1}^L& {\tilde{\mathbf{M}}}_{N_{z}2}^{\cdots }& \cdots & {\tilde{\mathbf{M}}}_{N_z,N_z}^L\end{array}\right]\hspace{0.17em},\hspace{0.17em}\hspace{0.17em}{\tilde{\mathbf{M}}}_{i,j}^L={\mathbf{M}}_{i,j}^L-{\mathbf{M}}_{1,j}^L\hspace{0.17em}\end{array}\}$

(45)

In the case of shallow water waves, the elements do not need to be divided vertically in a two-dimensional model because the velocity distribution is constant vertically (Ishikawa et al., 2016). The shallow water model due to Lepelletier and Raichlen (1988) involves one-dimensional equations using the average flow velocity, which can be used because the flow velocity distribution is constant. On the other hand, deep water waves have a velocity distribution; therefore, a two-dimensional model is necessary. However, in the proposed approach, the two-dimensional model can be transformed into the reduced-DOF model keeping the information about the velocity distribution by assuming the average velocity. The velocity distribution can be obtained by the inverse transformation of equation (42) from the results of the reduced-DOF model. As mentioned above, the meaning of average velocity in the shallow-water wave model differs from that in the proposed model.

Comparison of natural modes in the rectangular tank

To validate the proposed model in the linear region, we compare the numerical results calculated using the concentrated mass model (two-dimensional model and reduced-DOF model) with the corresponding theoretical values.

Analytical method using the concentrated mass model

Using the linearized equation (33) and (35) and assuming $\mathbf{C}=\mathbf{G}=\mathbf{0}$ in equation (24), we obtain the linearized equation of motion of the two-dimensional model

${\mathbf{M}}^L\ddot{\mathbf{x}}+{\mathbf{K}}^L\mathbf{x}=\mathbf{0}$

(46)

The components in the linearized matrices M^L and K^L are constants.

By solving the eigenvalue problem of equation (46), we obtain the order-r natural angular frequency ω _r and modal vector ${\mathbf{X}}_{\mathrm{r}}$ of size $N_z(N_x-1)\times 1$ . Meanwhile, using the linearized transfer matrix in equation (45), we obtain the linearized equation of motion of the reduced-DOF model as follows

$\begin{array}{l}\overline{\mathbf{M}}\ddot{\overline{\mathbf{x}}}+\overline{\mathbf{K}}\overline{\mathbf{x}}=\mathbf{0}\\ \overline{\mathbf{M}}={\left({\mathbf{T}}^L\right)}^T{\mathbf{M}}^L{\mathbf{T}}^L,\hspace{0.17em}\hspace{0.17em}\overline{\mathbf{K}}={\left({\mathbf{T}}^L\right)}^T{\mathbf{K}}^L{\mathbf{T}}^L\hspace{0.17em}\end{array}\}$

(47)

from which we obtain the order-r natural angular frequency ω _r and modal vector ${\overline{\mathbf{X}}}_{\mathrm{r}}$ of size $N_x-1$ . From equation (43), the modal vector ${\mathbf{X}}_{\mathrm{r}}$ of the two-dimensional model and of size $N_z\times (N_x-1)$ is reproduced from ${\overline{\mathbf{X}}}_{\mathrm{r}}$ as

${\mathbf{X}}_r={\mathbf{T}}^L{\overline{\mathbf{X}}}_r$

(48)

Validating the concentrated mass model in the linear region

From small-amplitude wave theory, the order-r natural frequency f _r of the rectangular tank and the order-r natural modes u _r and w _r for the horizontal and vertical velocities, respectively, are given by Faltinsen and Timokha (2009)

$\begin{array}{l}{f}_r=\frac{1}{2\pi }\sqrt{g{\kappa }_r\tanh \left({\kappa }_r{L}_z\right)}\hspace{0.17em}\\ {\kappa }_r=\frac{r\pi }{L_x}\end{array}\}$

(49)

$\begin{array}{l}{u}_r\left(x,z\right)=A\sin {\kappa }_{r}x\cosh {\kappa }_r\left(z+L_z\right)\\ w_r\left(x,z\right)=-A\cos {\kappa }_{r}x\sinh {\kappa }_r\left(z+L_z\right)\end{array}\}$

(50)

Table 1 lists the parameter values. The division numbers for the concentrated mass model are N _x = 100 and N _z =50. The natural frequencies and natural modes of the two-dimensional model are obtained by solving the eigenvalue problem of equation (46), and those of the reduced-DOF model are obtained by solving the eigenvalue problem of 47).

Table 1.

Parameter values.

L x [m]	0.5	L z [m]	0.25
ρ [kg/m3]	1,000	g [m/s2]	9.81

In Table 2, the natural frequencies given by the concentrated mass model (two-dimensional and reduced-DOF models) are compared with the theoretical values given by equation (49). In addition, the theoretical natural modes given numerically and by equation (50) are shown in Figure 9. The numerical results for the natural frequencies and modal shapes agree well with the theoretical ones. Therefore, the proposed concentrated mass model is valid for sloshing in the linear region.

Table 2.

Comparison of natural frequencies.

Order	Theoretical value [Hz]	Two-dimensional model [Hz]	Reduced-DOF model [Hz]
––	––	0 (4851 eigenvalues)	—
1	1.1966	1.1966	1.1966
2	1.7638	1.7637	1.7637
3	2.1641	2.1637	2.1637
4	2.4990	2.4982	2.4982
5	2.7940	2.7926	2.7926
6	3.0607	3.0584	3.0584

Figure 9.

Comparison of shapes of natural modes. (a) Theoretical modes (equation (50)). (b) Two-dimensional model. (c) Reduced-DOF model.

However, zero-frequency natural frequencies appear in the two-dimensional model (Table 2) because the stiffness matrix is singular. There are 4851 zero-frequency eigenvalues, and that number increases with the number of vertical separations. With a zero-frequency natural mode, the water surface does not move. By contrast, the reduced-DOF model has no zero-frequency natural frequencies, and its other natural frequencies agree exactly with those of the two-dimensional model. In addition, the reduced-DOF model reproduces the natural modal shapes of the two-dimensional model (Figure 9(c)). Therefore, our method for reducing the number of DOFs is valid.

This linear problem has the theoretical solution of equation (49) and (50). However, actual sloshing is a nonlinear phenomenon and differs from the theoretical solution when the amplitude is large. We confirm the validity of the nonlinear model in the section Comparison With Experimental Results.

Comparison with experimental results

In this section, numerical results from the nonlinear concentrated mass model are compared with experimental results to assess the validity of the proposed model in the nonlinear region.

Experimental apparatus and analysis method

To validate the proposed models, we performed sloshing experiments in a rectangular tank, as shown in Figure 10. The rectangular tank of width 100 mm was formed from 10-mm-thick acrylic boards. It was set on a linear guide and excited horizontally by a vibration exciter. A capacitance-type water-level gauge was installed at a distance of L_s from the left wall. We performed experiments in four configurations, which we refer to as cases 1–4. The values of the horizontal length L _x , water depth L _z , and excitation displacement Δ in each case are given in Table 3. In cases 1 and 2, the water was deeper than in cases 3 and 4; the latter two cases satisfy the shallow-water wave condition (i.e., the water depth is less than one-fifth of the wavelength).We assume that the inertial force $f^{\text{I}}\left(t\right)$ expressed by

$f^{\text{I}}\left(t\right)=-m\mathrm{\Delta }{\mathit{\Omega }}^2\sin \mathit{\Omega }t$

(51)

acts on each mass point by horizontal excitation, where Ω is the excitation angular frequency. The frequency response is calculated using the harmonic balance method (Urabe, 1965; Kondou et al., 1986). The harmonic balance method is a frequency-domain numerical method based on Galerkin approximation using Fourier series to find a steady periodic solution; the stability of the periodic solution can also be determined using the harmonic balance method. The proposed model also can be calculated by numerical integration. The order of the Fourier series in the harmonic balance method is eight, and the analysis considers nonlinearities; the calculation parameters are given in Table 3. The value of γ is determined so that the amplitudes of the numerical results agree with those of the experimental results in Figure 11. The values of the density ρ and the gravity acceleration g are listed in Table 1. The division numbers for the concentrated mass model are N _x = 30 and N _z = 10. As explained in the section Reducing the Number of Degrees of Freedom, the transformation matrix T in equation (42) contains the displacements δx_{i, j} and so must be calculated at each time step. In this section, the linearized transfer matrix T^L in equation (45) is used for the simulations.

Figure 10.

Experimental apparatus.

Table 3.

Parameter values for each case.

	L x [mm]	L z [mm]	L s [mm]	Δ [mm]	γ
Case 1	500	250	30	0.25	0.3
Case 2	500	250	30	0.3	0.25
Case 3	1,000	50	50	2.25	0.05
Case 4	1,000	50	50	3.5	0.087

Figure 11.

Frequency responses (reduced-DOF model). (a) Case 1. (b) Case 2. (c) Case 3. (d) Case 4.

Experimental and numerical results

For cases 1 and 2 (deep water), Figure 11(a) and (b) show the frequency responses around the first-order resonance. The linear first-order natural frequency is 1.20 Hz. The resonance curves tilt slightly to the left in these cases. In general, the resonance curve tilts to the left when L _z /L _x > 0.337 (Flanagan and Timokha, 2009), and in cases 1 and 2 we have L _z /L _x = 0.5. The numerical results agree well with the experimental results, and the numerical results calculated using the reduced-DOF model agree well with those calculated using the two-dimensional model.

Figures 11(c) and (d) show the frequency responses around the first-order resonance for cases 3 and 4 (shallow water) with the reduced-DOF model. The first-order natural frequency is 0.349 Hz. The resonance curves tilt to the right in the case of shallow water. In general, the resonance curve tilts to the right when L _z /L _x < 0.337, and in cases 3 and 4 we have L _z /L _x = 0.05. In cases 3 and 4, there are some peaks around the first-order resonance, and at these peak frequencies, solitary waves move back and forth in the tank (Ishikawa et al., 2016). From Figure 11, the numerical results from the proposed model agree well with the experimental results. Therefore, the proposed model is valid for numerical analysis of nonlinear sloshing phenomena.

Table 4 gives the natural frequencies in cases 1–4. For cases 3 and 4, half of the second-order natural frequency is 0.3446 Hz (= 0.6892/2), which is close to the frequency of (iii) in Figure 11(c); therefore, second-order superharmonic resonance occurs at (iii). A third of the third-order natural frequency is 0.3380 Hz, which is close to the frequency of (ii), and a quarter of the fourth-order natural frequency is 0.3297 Hz, which is close to the frequency of (i); therefore, these peaks are superharmonic resonances. Meanwhile, for cases 1 and 2, half of the second-order natural frequency is 0.8819 Hz, and a third of the third-order natural frequency is 0.7214 Hz, but these frequencies are far from the first-order natural frequency (1.1966 Hz); therefore, there is only one peak around the first-order natural frequency in the case of deep water, as shown in Figure 11(a) and (b).

Table 4.

Natural frequencies in each case.

Order	Cases 1 and 2 [Hz]	Cases 3 and 4 [Hz]
1	1.1966	0.3487
2	1.7638	0.6892
3	2.1641	1.0142
4	2.4990	1.3187
5	2.7940	1.5999
6	3.0607	1.8572

Conclusions

We proposed a concentrated mass model for analyzing nonlinear two-dimensional sloshing in a rectangular tank, and our conclusions are as follows.

1. A concentrated mass model was proposed for two-dimensional sloshing in a rectangular tank excited horizontally. The restoring forces were derived from static and dynamic pressures, the latter of which considers vertical acceleration of the z nodal points.

2. A method for reducing the number of DOFs of the proposed two-dimensional model was proposed, by which the number of DOFs in the vertical direction can be reduced to one.

3. The advantage of the proposed method is that DOFs of the model are less than that of the two-dimensional finite-element method, thereby reducing the computation time dramatically.

4. The natural frequencies in a rectangular tank as obtained from the concentrated mass model (two-dimensional and reduced-DOF models) agreed well with the theoretical natural frequencies.

5. The numerical results agreed well with the experimental results both qualitatively and quantitatively.

6. In a tank with complex shape, the water depth varies from location to location. In this case, the water area is divided into rectangular elements, and the DOF reduction is performed at each location to apply the proposed model to the tank with complex shape. In the case of a three-dimensional problem, the water region is divided into rectangular parallel-piped elements, and mass points are located at the faces of the elements. However, these models have not been validated, and doing so is a future task, along with modeling the damping.