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.

Keywords

Sloshing modeling storage tank nonlinear dynamics concentrated mass 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 = ρ 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

(δ z_{i - 1, j} - δ z_{i, j}) (l_{x} + δ x_{i, j} - δ x_{i, j - 1}) = - l_{z} (δ x_{i, j} - δ x_{i, j - 1})

(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

δ z_{i - 1, j} - δ z_{i, j} = - \frac{l_{z} (δ x_{i, j} - δ x_{i, j - 1})}{l_{x} + δ x_{i, j} - δ 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

δ 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

δ z_{i, j} = - \sum_{k = i + 1}^{N_{z}} \frac{l_{z} (δ x_{k, j} - δ x_{k, j - 1})}{l_{x} + δ x_{k, j} - δ x_{k, j - 1}} (0 \leq i \leq N_{z} - 1)

(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

δ {\bar{x}}_{j}

as the average displacement of the mass points in column j, as shown in Figure 4. The average displacement is given by

δ {\bar{x}}_{j} = \frac{1}{N_{z}} \sum_{i = 1}^{N_{z}} δ x_{i, j}

(6)

We use

δ {\bar{x}}_{j}

when considering the extruded volume in the z direction. Replacing

δ x_{i, j} - δ x_{i, j - 1}

on the left-hand side in equation 2) with

δ {\bar{x}}_{j} - δ {\bar{x}}_{j - 1}

, the relative displacement of an arbitrary z nodal point becomes

δ z_{i - 1, j} - δ z_{i, j} = - \frac{l_{z} (δ x_{i, j} - δ x_{i, j - 1})}{l_{x} + δ {\bar{x}}_{j} - δ {\bar{x}}_{j - 1}}

(7)

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

δ z_{i, j} = - \frac{l_{z}}{l_{x} + δ {\bar{x}}_{j} - δ {\bar{x}}_{j - 1}} \sum_{k = i + 1}^{N_{z}} (δ x_{k, j} - δ x_{k, j - 1}) (0 \leq i \leq N_{z} - 1)

(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

η_{j} = δ z_{0, j} = - \frac{l_{z} N_{z} (δ {\bar{x}}_{j} - δ {\bar{x}}_{j - 1})}{l_{x} + δ {\bar{x}}_{j} - δ {\bar{x}}_{j - 1}} = - \frac{L_{z} (δ {\bar{x}}_{j} - δ {\bar{x}}_{j - 1})}{l_{x} + δ {\bar{x}}_{j} - δ {\bar{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} (z, t)}{\partial t} = - \frac{1}{ρ} \frac{\partial p_{j} (z, t)}{\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 \geq 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}{l} p_{1, j} = ρ \int_{(η_{j} - l_{z}) / 2}^{η_{j}} \frac{\partial w_{j}}{\partial t} d z + \frac{1}{2} ρ g (η_{j} + l_{z}) & (i = 1) \\ p_{i, j} = ρ \int_{z_{i}}^{η_{j}} \frac{\partial w_{j}}{\partial t} d z + ρ g {η_{j} + (i - 1) l_{z} + \frac{l_{z}}{2}} & (i \geq 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 $δ {\ddot{\overset{⌢}{z}}}_{i, j}$ is the average vertical acceleration of z nodal points ${(i - 1, j)}_{z}$ and ${(i, j)}_{z}$ (Figure 7), $δ {\ddot{\overset{⌢}{z}}}_{i, j}$ is

δ {\ddot{\overset{⌢}{z}}}_{i, j} = \frac{δ {\ddot{z}}_{i - 1, j} + δ {\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

δ {\ddot{\overset{⌢}{z}}}_{i, j}

. In the integration in equation (11),

\partial w_{j} / \partial t

is replaced by the vertical average of vertical accelerations

δ {\ddot{\overset{⌢}{z}}}_{i, j}

. The integration in equation (11) is replaced by the summation using

{\ddot{\overset{⌢}{z}}}_{i, j}

from element [1, j] to element [i, j] as

\begin{array}{l} f_{1, j} = \frac{ρ}{2} δ {\ddot{\overset{⌢}{z}}}_{1, j} {(η_{j} + l_{z})}^{2} + \frac{ρ}{2} g (η_{j} + l_{z})^{2} & (i = 1) \\ f_{i, j} = ρ δ {\ddot{\overset{⌢}{z}}}_{1, j} l_{z} (η_{j} + l_{z}) + ρ \sum_{k = 2}^{i - 1} δ {\ddot{\overset{⌢}{z}}}_{k, j} l_{z}^{2} + ρ δ {\ddot{\overset{⌢}{z}}}_{i, j} \frac{l_{z}^{2}}{2} + ρ g l_{z} {η_{j} + (i - 1) l_{z} + \frac{l_{z}}{2}} & (2 \leq i \leq 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}{l} β_{1, j} = \frac{1}{2} ρ δ {\ddot{\overset{⌢}{z}}}_{1, j} (η_{j} + l_{z})^{2} & (i = 1) \\ β_{i, j} = ρ {δ {\ddot{\overset{⌢}{z}}}_{1, j} l_{z} (η_{j} + l_{z}) + \sum_{k = 2}^{i - 1} δ {\ddot{\overset{⌢}{z}}}_{k, j} l_{z}^{2} + \frac{1}{2} δ {\ddot{\overset{⌢}{z}}}_{i, j} l_{z}^{2}} & (2 \leq i \leq 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} = ρ g l_{z}^{2} (i - \frac{1}{2})

(15)

Therefore, the restoring force

δ f_{i, j}^{k}

from the pressure is written as

δ f_{i, j}^{k} = f_{i, j} - f_{i, j}^{0} = α_{i, j} + β_{i, j}

(16)

\begin{array}{l} α_{1, j} = ρ g l_{z} η_{j} + \frac{1}{2} ρ g η_{j}^{2} & (i = 1) \\ α_{i, j} = ρ g l_{z} η_{j} & (2 \leq i \leq N_{z}) \end{array}}

(17)

\begin{array}{l} β_{1, j}^{L} = \frac{1}{2} ρ δ {\ddot{\overset{⌢}{z}}}_{1, j} l_{z}^{2} & (i = 1) \\ β_{i, j}^{L} = ρ (δ {\ddot{\overset{⌢}{z}}}_{1, j} l_{z}^{2} + \sum_{k = 1}^{i - 1} δ {\ddot{\overset{⌢}{z}}}_{k, j} l_{z}^{2} + \frac{1}{2} δ {\ddot{\overset{⌢}{z}}}_{i, j} l_{z}^{2}) & (i \geq 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}^{L} = - \frac{l_{z}}{l_{x}} \sum_{k = i + 1}^{N_{z}} (δ x_{k, j} - δ x_{k, j - 1})

(19)

η_{j}^{L} = δ z_{0, j} = - \frac{l_{z}}{l_{x}} \sum_{k = 1}^{N_{z}} (δ x_{k, j} - δ x_{k, j - 1})

(20)

α_{i, j}^{L} = ρ g l_{z} η_{j}^{L}

(21)

\begin{array}{l} M \ddot{x} + C \dot{x} + K x = G \\ G = \frac{1}{2} ρ g {[\begin{matrix} σ^{T} & 0 & \dots & 0 \end{matrix}]}^{T} \\ σ = {[\begin{matrix} η_{1}^{2} - η_{2}^{2} & η_{2}^{2} - η_{3}^{2} & \dots & η_{N - 1}^{2} - η_{N}^{2} \end{matrix}]}^{T} \end{array}}

(24)

Note that the average displacement

δ {\bar{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 δ {\ddot{x}}_{i, j} = α_{i, j} + β_{i, j} - α_{i, j + 1} - β_{i, j + 1} - f_{i, j}^{c} + f_{i, j + 1}^{c} + f_{i, j}^{e} (t)

(23)

Here,

f_{i, j}^{c}

is the linear viscous damping force, and

f_{i, j}^{e} (t)

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

δ x_{i, N_{x}}

are zero. Writing the equations of motion for all the mass points and arranging them in matrix form, we obtain

M_{i, j}^{a} = [\begin{matrix} m_{j, 1}^{z 1} m_{i, 1}^{z 2} + m_{j, 2}^{z 1} m_{i, 2}^{z 2} & - m_{j, 2}^{z 1} m_{i, 2}^{z 2} \\ - m_{j, 2}^{z 1} m_{i, 2}^{z 2} & m_{j, 2}^{z 1} m_{i, 2}^{z 2} + m_{j, 3}^{z 1} m_{i, 3}^{z 2} & - m_{j, 3}^{z 1} m_{i, 3}^{z 2} \\ ⋱ & ⋱ & ⋱ \\ - m_{j, N_{x} - 1}^{z 1} m_{i, N_{x} - 1}^{z 2} & m_{j, N_{x} - 1}^{z 1} m_{i, N_{x} - 1}^{z 2} + m_{j, N_{x}}^{z 1} m_{i, N_{x}}^{z 2} \end{matrix}]

(29)

where M, C, and K are the mass, damping, and stiffness matrices, respectively, of size

N_{z} (N_{x} - 1) \times N_{z} (N_{x} - 1)

. 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} (N_{x} - 1) \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} x = {[\begin{matrix} x_{1}^{T} & x_{2}^{T} & \dots & x_{N_{z}}^{T} \end{matrix}]}^{T} \\ x_{i} = {[\begin{matrix} δ x_{i, 1} & δ x_{i, 2} & \dots & δ x_{i, N_{x} - 1} \end{matrix}]}^{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

M = M^{m} + M^{d} = [\begin{matrix} M_{11} & M_{12} & \dots & M_{1 N_{z}} \\ M_{21} & M_{22} & \dots & M_{2 N_{z}} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ M_{N_{z} 1} & M_{N_{z} 2} & \dots & M_{N_{z} N_{z}} \end{matrix}]

(26)

where

M^{m}

and

M^{d}

are

M^{m} = [\begin{matrix} m & 0 & \dots & 0 \\ 0 & m & ⋱ & ⋮ \\ ⋮ & ⋱ & ⋱ & 0 \\ 0 & \dots & 0 & m \end{matrix}], M^{d} = [\begin{matrix} M_{1}^{1} / 2 & M_{2}^{1} & M_{3}^{1} & \dots & M_{N_{z}}^{1} \\ M_{2,1}^{b} & M_{2,2}^{a} & M_{2,3}^{b} & \dots & M_{2, N_{z}}^{b} \\ M_{3,1}^{b} & M_{3,2}^{b} & M_{3,3}^{a} & \dots & M_{3, N_{z}}^{b} \\ ⋮ & ⋮ & ⋮ & ⋱ & ⋮ \\ M_{N_{z}, 1}^{b} & M_{N_{z}, 2}^{b} & M_{N_{z}, 3}^{b} & \dots & M_{N_{z}, N_{z}}^{a} \end{matrix}]

(27)

M^{m}

and

M^{d}

represent the mass of liquid and the dynamic pressure, respectively.

M_{i}^{1}

, which means the dynamic pressure acting on the masses in the first row, is written as

\begin{array}{l} M_{i}^{1} = [\begin{matrix} m_{i, 1}^{z 1} + m_{i, 2}^{z 1} & - m_{i, 2}^{z 1} \\ m_{i, 2}^{z 1} & m_{i, 2}^{z 1} + m_{i, 3}^{z 1} & - m_{i, 3}^{z 1} \\ ⋱ & ⋱ & ⋱ \\ - m_{i, N_{x} - 1}^{z 1} & m_{i, N_{x} - 1}^{z 1} + m_{i, N_{x}}^{z 1} \end{matrix}] \\ m_{i, j}^{z 1} = \frac{ρ {(l_{z} + η_{j})}^{2} l_{x} l_{z}}{2 {(l_{x} + δ x_{i, j} - δ x_{i, j - 1})}^{2}} \end{array}} .

(28)

M_{i, j}^{a}

and

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}{l} \int_{(η_{j} - l_{z}) / 2}^{η_{j}} \frac{\partial w_{j}}{\partial t} d z = \frac{δ {\ddot{\overset{⌢}{z}}}_{1, j} (η_{j} + l_{z})}{2} & (i = 1) \\ \int_{z_{i}}^{η_{j}} \frac{\partial w_{j}}{\partial t} d z = δ {\ddot{\overset{⌢}{z}}}_{1, j} (l_{z} + η_{j}) + \sum_{k = 2}^{i - 1} δ {\ddot{\overset{⌢}{z}}}_{k, j} l_{z} + δ {\ddot{\overset{⌢}{z}}}_{i, j} \frac{l_{z}}{2} & (2 \leq i < N_{z}) \end{array}}

(13)

M_{i, j}^{b} = [\begin{matrix} m_{j, 1}^{z 1} m_{i, 1}^{z 3} + m_{j, 2}^{z 1} m_{i, 2}^{z 3} & - m_{j, 2}^{z 1} m_{i, 2}^{z 3} \\ - m_{j, 2}^{z 1} m_{i, 2}^{z 3} & m_{j, 2}^{z 1} m_{i, 2}^{z 3} + m_{j, 3}^{z 1} m_{i, 3}^{z 3} & - m_{j, 3}^{z 1} m_{i, 3}^{z 3} \\ ⋱ & ⋱ & ⋱ \\ - m_{j, N_{x} - 1}^{z 1} m_{i, N_{x} - 1}^{z 3} & m_{j, N_{x} - 1}^{z 1} m_{i, N_{x} - 1}^{z 3} + m_{j, N_{x}}^{z 1} m_{i, N_{x}}^{z 3} \end{matrix}]

(30)

where

m_{i, k}^{z 2} = (l_{z} + η_{k}) l_{z} + (i - \frac{7}{4}) l_{z}^{2}, m_{i, k}^{z 3} = (l_{z} + η_{k}) l_{z} + (i - \frac{3}{2}) l_{z}^{2}

(31)

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

\begin{array}{l} M_{11} = diag [\begin{matrix} m & m & \dots & m \end{matrix}] + \frac{M_{1}^{1}}{2} \\ M_{1 j} = M_{j}^{1} (j \geq 2) \\ M_{i i} = diag [\begin{matrix} m & m & \dots & m \end{matrix}] + M_{i i}^{a} \\ M_{i j} = M_{i j}^{b} (i \neq j) \end{array}}

(32)

Regarding the linearized mass matrix M^L,

m_{i, j}^{z 1}

in equation (28) and

m_{i, k}^{z 2}

and

m_{i, k}^{z 3}

in equation (31) are treated as

m_{i, j}^{z 1, L} = \frac{ρ l_{z}^{3}}{l_{x}}, m_{i, k}^{z 2, L} = l_{z}^{2} + (i - \frac{7}{4}) l_{z}^{2}, m_{i, k}^{z 3, L} = l_{z}^{2} + (i - \frac{3}{2}) l_{z}^{2}

(33)

The components of the stiffness matrix are written as follows:

\begin{array}{l} K = [\begin{matrix} \hat{K} & \hat{K} & \dots & \hat{K} \\ \hat{K} & \hat{K} & \dots & \hat{K} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ \hat{K} & \hat{K} & \dots & \hat{K} \end{matrix}] \\ \hat{K} = [\begin{matrix} k_{1} + k_{2} & - k_{2} & 0 & \dots & 0 \\ - k_{2} & k_{2} + k_{3} & - k_{3} & ⋱ & ⋮ \\ 0 & ⋱ & ⋱ & ⋱ & 0 \\ ⋮ & ⋱ & - k_{N - 2} & k_{N - 2} + k_{N - 1} & - k_{N - 1} \\ 0 & \dots & 0 & - k_{N - 1} & k_{N - 1} + k_{N} \end{matrix}], k_{j} = \frac{ρ g l_{z}^{2}}{l_{x} + δ {\bar{x}}_{j} - δ {\bar{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} C = γ K^{L} \\ K^{L} = [\begin{matrix} \tilde{K} & \tilde{K} & \dots & \tilde{K} \\ \tilde{K} & \tilde{K} & \dots & \tilde{K} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ \tilde{K} & \tilde{K} & \dots & \tilde{K} \end{matrix}], \tilde{K} = [\begin{matrix} 2 \tilde{k} & - \tilde{k} & 0 & \dots & 0 \\ - \tilde{k} & 2 \tilde{k} & - \tilde{k} & ⋱ & M \\ 0 & ⋱ & ⋱ & ⋱ & 0 \\ M & ⋱ & - \tilde{k} & 2 \tilde{k} & - \tilde{k} \\ 0 & \dots & 0 & - \tilde{k} & 2 \tilde{k} \end{matrix}], \tilde{k} = \frac{ρ 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

[\begin{matrix} M_{11} & M_{12} & \dots & M_{1 N_{z}} \\ M_{21} & M_{22} & \dots & M_{2 N_{z}} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ M_{N_{z} 1} & M_{N_{z} 2} & \dots & M_{N_{z} N_{z}} \end{matrix}] [\begin{matrix} {\dot{x}}_{1} \\ {\dot{x}}_{2} \\ ⋮ \\ {\dot{x}}_{N_{z}} \end{matrix}] + γ [\begin{matrix} \tilde{K} & \tilde{K} & \dots & \tilde{K} \\ \tilde{K} & \tilde{K} & \dots & \tilde{K} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ \tilde{K} & \tilde{K} & \dots & \tilde{K} \end{matrix}] [\begin{matrix} {\dot{x}}_{1} \\ {\dot{x}}_{2} \\ ⋮ \\ {\dot{x}}_{N_{z}} \end{matrix}] + [\begin{matrix} \hat{K} & \hat{K} & \dots & \hat{K} \\ \hat{K} & \hat{K} & \dots & \hat{K} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ \hat{K} & \hat{K} & \dots & \hat{K} \end{matrix}] [\begin{matrix} {\dot{x}}_{1} \\ {\dot{x}}_{2} \\ ⋮ \\ {\dot{x}}_{N_{z}} \end{matrix}] = \frac{1}{2} ρ g [\begin{matrix} σ \\ 0 \\ ⋮ \\ 0 \end{matrix}]

(36)

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

[\begin{matrix} M_{11} & M_{12} & \dots & M_{1 N_{z}} \\ {\tilde{M}}_{21} & {\tilde{M}}_{22} & \dots & {\tilde{M}}_{2 N_{z}} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ {\tilde{M}}_{N_{z} 1} & {\tilde{M}}_{N_{z} 2} & \dots & {\tilde{M}}_{N_{z} N_{z}} \end{matrix}] [\begin{matrix} {\dot{x}}_{1} \\ {\dot{x}}_{2} \\ ⋮ \\ {\dot{x}}_{N_{z}} \end{matrix}] + γ [\begin{matrix} \tilde{K} & \tilde{K} & \dots & \tilde{K} \\ 0 & 0 & \dots & 0 \\ ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & 0 & \dots & 0 \end{matrix}] [\begin{matrix} {\dot{x}}_{1} \\ {\dot{x}}_{2} \\ ⋮ \\ {\dot{x}}_{N_{z}} \end{matrix}] + [\begin{matrix} \hat{K} & \hat{K} & \dots & \hat{K} \\ 0 & 0 & \dots & 0 \\ ⋮ & ⋮ & ⋱ & ⋮ \\ 0 & 0 & \dots & 0 \end{matrix}] [\begin{matrix} {\dot{x}}_{1} \\ {\dot{x}}_{2} \\ ⋮ \\ {\dot{x}}_{N_{z}} \end{matrix}] = \frac{1}{2} ρ g [\begin{matrix} σ \\ - σ \\ ⋮ \\ - σ \end{matrix}]

(37)

Here,

{\tilde{M}}_{i, j}

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

(38)

In addition, we introduce the average displacement vector

\bar{x}

, whose elements are the average displacements in equation (6). The average displacement vector can be written as the summation of

x_{1}

x_{N_{z}}

\bar{x} = [\begin{matrix} δ {\bar{x}}_{1} \\ δ {\bar{x}}_{2} \\ ⋮ \\ δ {\bar{x}}_{N_{x} - 1} \end{matrix}] = \frac{1}{N_{z}} [\begin{matrix} I_{N_{x} - 1} & I_{N_{x} - 1} & \dots & I_{N_{x} - 1} \end{matrix}] [\begin{matrix} x_{1} \\ x_{2} \\ \dots \\ x_{N_{z}} \end{matrix}]

(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

[\begin{matrix} I_{N_{x} - 1} & I_{N_{x} - 1} & \dots & I_{N_{x} - 1} \\ {\tilde{M}}_{21} & {\tilde{M}}_{22} & \dots & {\tilde{M}}_{2 N_{z}} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ {\tilde{M}}_{N_{z} 1} & {\tilde{M}}_{N_{z} 2} & \dots & {\tilde{M}}_{N_{z} N_{z}} \end{matrix}] [\begin{matrix} {\dot{x}}_{1} \\ {\dot{x}}_{2} \\ ⋮ \\ {\dot{x}}_{N_{z}} \end{matrix}] = N_{z} [\begin{matrix} I_{N_{x} - 1} \\ 0 \\ ⋮ \\ 0 \end{matrix}] \ddot{\bar{x}} - \frac{1}{2} ρ g [\begin{matrix} 0 \\ σ \\ ⋮ \\ σ \end{matrix}]

(40)

which can be written as

\begin{array}{l} A_{1} \ddot{x} = A_{2} \ddot{\bar{x}} - \tilde{S} \\ A_{1} = [\begin{matrix} I_{N_{x} - 1} & I_{N_{x} - 1} & \dots & I_{N_{x} - 1} \\ {\tilde{M}}_{21} & {\tilde{M}}_{22} & \dots & {\tilde{M}}_{2 N_{z}} \\ ⋮ & ⋮ & \dots & ⋮ \\ {\tilde{M}}_{N_{z} 1} & {\tilde{M}}_{N_{z} 2} & \dots & {\tilde{M}}_{N_{z} N_{z}} \end{matrix}], A_{2} = N_{z} [\begin{matrix} I_{N_{x} - 1} \\ 0 \\ ⋮ \\ 0 \end{matrix}], \tilde{S} = \frac{1}{2} ρ g [\begin{matrix} 0 \\ σ \\ ⋮ \\ σ \end{matrix}] \end{array}}

(41)

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

\begin{array}{l} \ddot{x} = T \ddot{\bar{x}} - S \\ T = A_{1}^{- 1} A_{2}, S = A_{1}^{- 1} \tilde{S} \end{array}}

(42)

which represents the relationship between the two-dimensional acceleration

\ddot{x}

and the one-dimensional acceleration

\ddot{\bar{x}}

. T is the transformation matrix from the two-dimensional model to the one-dimensional model. We assume that the displacement vectors x and

\bar{x}

and the velocity vectors

\dot{x}

and

\dot{\bar{x}}

are also related by T

\begin{array}{l} x = T \bar{x} \\ \dot{x} = T \dot{\bar{x}} \end{array}}

(43)

Substituting equation (42) and (43) into equation (24) and left-multiplying by

T^{T}

, we obtain

\begin{array}{l} \bar{M} \ddot{\bar{x}} + \bar{C} \dot{\bar{x}} + \bar{K} \bar{x} = \bar{G} + T^{T} M S \\ \bar{M} = T^{T} M T, \bar{C} = T^{T} C T, \bar{K} = T^{T} K T, \bar{G} = T^{T} 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

\bar{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{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} T^{L} = {(A_{1}^{L})}^{- 1} A_{2} \\ A_{1}^{L} = [\begin{matrix} I_{N_{x} - 1} & I_{N_{x} - 1} & \dots & I_{N_{x} - 1} \\ {\tilde{M}}_{21}^{L} & {\tilde{M}}_{22}^{L} & \dots & {\tilde{M}}_{2 N_{z}}^{L} \\ ⋮ & ⋮ & ⋱ & ⋮ \\ {\tilde{M}}_{N_{z} 1}^{L} & {\tilde{M}}_{N_{z} 2}^{\dots} & \dots & {\tilde{M}}_{N_{z}, N_{z}}^{L} \end{matrix}], {\tilde{M}}_{i, j}^{L} = M_{i, j}^{L} - M_{1, j}^{L} \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 $C = G = 0$ in equation (24), we obtain the linearized equation of motion of the two-dimensional model

M^{L} \ddot{x} + K^{L} x = 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 $X_{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} \bar{M} \ddot{\bar{x}} + \bar{K} \bar{x} = 0 \\ \bar{M} = {(T^{L})}^{T} M^{L} T^{L}, \bar{K} = {(T^{L})}^{T} K^{L} T^{L} \end{array}}

(47)

from which we obtain the order-r natural angular frequency ω_r and modal vector

{\bar{X}}_{r}

of size

N_{x} - 1

. From equation (43), the modal vector

X_{r}

of the two-dimensional model and of size

N_{z} \times (N_{x} - 1)

is reproduced from

{\bar{X}}_{r}

X_{r} = T^{L} {\bar{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 π} \sqrt{g κ_{r} \tanh (κ_{r} L_{z})} \\ κ_{r} = \frac{r π}{L_{x}} \end{array}}

(49)

\begin{array}{l} u_{r} (x, z) = A \sin κ_{r} x \cosh κ_{r} (z + L_{z}) \\ w_{r} (x, z) = - A \cos κ_{r} x \sinh κ_{r} (z + L_{z}) \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/m³]	1,000	g [m/s²]	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^{I} (t)$ expressed by

f^{I} (t) = - m Δ Ω^{2} \sin Ω 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.

Footnotes

Acknowledgments

The authors would like to thank Mr Yuki Amano at the Railway Technical Research Institute and Prof. Kenichiro Matsuzaki at Kagoshima University for discussions about this research. The experimental equipment was made by Mr Yosuke Koba at Kyushu University.

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

ORCID iD

Satoshi Ishikawa

References

Abrason

Chu

Ransleben

(1961) Representation of fuel slosing in cylindrical tanks by an equivalent mechanical model. ARS Jouranl 31(12): 1697–1705.

Bauer

(1966) Nonlinear mechanical model for the description of propellant sloshing. AIAA Journal 4(9): 1662–1668.

Drosos

Dimas

Karabalis

(2008) Discrete models for seismic analysis of liquid storage tanks of arbitrary shape and fill height. Journal of Pressure Vessel Technology 130: 041801.

Faltinsen

(1974) A nonlinear theory of sloshing in rectangular tank. Journal of Ship Research 18(4): 224–241.

Faltinsen

(1978) A numerical nonlinear method of sloshing in tanks with two-dimensional flow. Journal of Ship Research 22(3): 193–202.

Faltinsen

Rognebakke

Lukovsky

, et al. (2000) Multidimensional modal analysis of nonlinear sloshing in a rectangular tank with finite water depth. Journal of Fluid Mechanics 407: 201–234.

Faltinsen

Timokha

(2002) Asymptotic modal approximation of nonlinear resonant sloshing in a rectangular tank with small fluid depth. Journal of Fluid Mechanics 470: 319–357.

Faltinsen

Timokha

(2009) Sloshing. New York, NY: Cambridge University Press.

Graham

Rodriguez

(1952) The characteristics of fuel motion which affect airplane dynamics. Journal of Applied Mechanics 19(3): 381–388.

10.

Godderidge

Turnock

Tan

(2012) A rapid method for the simulation of sloshing using a mathematical model based on the pendulum equation. Computers & Fluids 57: 163–171.

11.

Hatayama

(2008) Lessons from the 2003 Tokachi-oki, Japan, earthquake for prediction of long-period strong ground motions and sloshing damage to oil storage tanks. Journal of Seismology 12(2): 255–263.

12.

Hayama

Aruga

Watanabe

(1983) Nonlinear responses of sloshing in rectangular tanks: 1st report, nonlinear responses of surface elevation. Bulletin of JSME 26(219): 1641–1648.

13.

Housner

(1957) Dynamic pressures on accelerated fluid containers. Bulletin of the Seismological Society of America 47: 15–35.

14.

Ikegawa

(2002) Finite element analysis of fluid motion in a container. In: Oden

(ed) Finite Element Methods in Flow ProblemsHuntsville, AL: UAH Press.

15.

Ishikawa

Kondou

Matsuzaki

, et al. (2016) Analysis of nonlinear shallow water waves in a tank by concentrated mass model. Journal of Sound and Vibration 371: 171–182.

16.

Kaneko

(ed) (2016) Flow-Induced Vibrations: Classifications and Lessons From Practical Experiences. 2nd edition. London: Academic Press.

17.

Komatsu

(1958) Non-linear sloshing analysis of liquid in tanks with arbitrary geometries. International Journal of Non-linear Mechanics 22(3): 193–207.

18.

Kondou

Tamura

Sueoka

(1986) On a method of higher approximation and determination of stability criterion for steady oscillations in nonlinear systems. Bulletin of JSME 29(248): 525–532.

19.

Lepelletier

Raichlen

(1988) Nonlinear oscillations in rectangular tanks. Journal of Engineering Mechanics 114(1): 1–23.

20.

Love

Tait

(2011) Non-linear multimodal model for tuned liquid dampers of arbitrary tank geometry. International Journal of Non-Linear Mechanics 46: 1065–1075.

21.

Love

Tait

(2015) The response of structures equipped with tuned liquid dampers of complex geometry. Journal of Vibration and Control 21(6): 1171–1187.

22.

Moiseev

(1958) On the theory of nonlinear vibrations of a liquid of finite volume. Journal of Applied Mathematics and Mechanics 22(5): 860–872.

23.

Nagashima

(2009) Sloshing analysis of a liquid storage container using level set X-FEM. Communications in Numerical Methods in Engineering 25: 357–379.

24.

Nakayama

Washizu

(1980) Nonlinear analysis of liquid motion in a container subjected to forced pitching oscillation. International Journal for Numerical Methods in Engineering 15: 1207–1220.

25.

Nakayama

Washizu

(1981) The boundary element method applied to the analysis of two-dimensional nonlinear sloshing problems. International Journal for Numerical Methods in Engineering 17: 1631–1646.

26.

Sumner

(1965) Experimentally determined pendulum analogy of liquid sloshing in speherical and oblate-spheroidal tanks. NASA Technical Note. NASA TN D-2737.

27.

Urabe

(1965) Galerkin’s procedure for nonlinear periodic systems. Archive for Rational Mechanics and Analysis 20: 120–152.

28.

Eatock Taylor

(1998) Numerical simulation of sloshing waves in a 3D tank based on a finite element method. Applied Ocean Research 20: 337–355.