Frequency Characteristics of Sloshing Resonance in a Three-Dimensional Shallow-Water Rectangular Tank

: The frequency characteristics of free surface elevation time histories of shallow-water sloshing in a three-dimensional rectangular tank is presented. The numerical model for sloshing motion uses an accurate velocity potential Boussinesq-type equation. A particular solution is adopted to express the external excitations and a linear damping term is introduced to replace a friction effect produced by tank walls. The excitation amplitude of the forced movement of the tank is tiny, while the excitation frequency varies around the ﬁrst resonance frequency. The numerical results show that the wave energy in the rectangular tank under resonance is concentrated in the integer multiple of excitation frequencies, the sum excitation frequencies, and the difference excitation frequencies.


Introduction
Sloshing and its effects have become a crucial part of the structural and navigational safety of liquid-carrying ships such as LNG or FPSO. In addition to preventing the hazards due to sloshing, TLD and anti-rolling tanks that take advantage of the sloshing effect are also concerns of researchers [1][2][3][4]. Different relative filling levels in tanks can lead to significantly different regimes of liquid sloshing resonance. For intermediate water depth or more filling level cases, the nonlinear standing waves, local impact loads, and wave breaking phenomenon have been studied by the experimental methods [5][6][7], the numerical models based on Navier-Stokes equations [8][9][10], and inviscid incompressible flow theories [11][12][13].
For the shallow-water sloshing resonance, new shallow-water equations were derived by [14] under the assumption of inviscid vortical flow. The Boussinesq-type depth-averaged equations were used in [15,16] for a two-dimensional shallow-water sloshing. The Fourier decomposition was adopted for the reconstruction of the Boussinesq-type equations, and an artificial viscosity term was used for dissipation. The time histories of free surface elevations under the coupled sway-heave-roll excitation case were compared with the results calculated by the Smoothed Particle Hydrodynamics method.
For the strong nonlinear shallow water sloshing motion usually accompanied by wave breaking and hydraulic jump, the methods based on the Navier-Stokes equations and experimental methods were widely adopted. The experimental study of nonlinear sloshing characteristics in different shallow-water tanks was given in [17]. The results show that the excitation amplitude and dispersion parameters can greatly affect the resonance frequency. The hydraulic jump phenomenon in the shallow-water rectangular tank was analyzed in [18]. The frequency range of occurrence of hydraulic jump in tanks with different aspect ratios was discussed, and it has been found that the liquid depth plays a significant role.
Shallow-water sloshing in a rectangular tank was investigated in [19][20][21] by the Boussinesq-type equations with velocity potential as the variable to be solved. The Fourier series was applied for the particular solution of velocity potential, and a linear damping term was added in the dynamic boundary condition. The free surface profiles in tanks under six degrees of freedom excitations were shown and the single-bump, double-bump, and triple-bump traveling waves can be found when the excitation frequency changes around the first natural frequency. When the external excitations around the z-axis were added, the swirling waves appeared and were coupled with the nonlinear traveling waves.
The free surface profiles of shallow-water sloshing have an extremely complex wave compositions. In this paper, the power spectral densities of free surface elevations in tanks are used to analyze the wave compositions of the complex free surface profiles. Firstly, the numerical model based on the Boussinesq-type equations is built, and the linear damping term is given. Then, the numerical model is validated by comparing some calculations with experimental data in the literature. Finally, the power spectral densities of free surface elevations in the tank which suffered from different forced movements are shown. The results shown in this paper are useful for further understanding the regimes of nonlinear shallow-water sloshing motions.

Statement of the Problem
Sloshing of liquid in a rectangular tank is considered, symbol l is the length of the tank, b is the width, and h is the filling level. As shown in Figure 1, a Cartesian coordinate system oxyz attached to the tank is adopted. xoy plane coincides with the still water plane and the origin o with the center of it. The sides and bottom of the fluid are confined by the tank walls, and the upper boundary is the free surface z = η(x, y, t). Assuming potential flow, there is a velocity potential Φ(x, y, z, t) in the fluid region. Nonlinear free surface boundary conditions are adopted. With the Laplace equations in fluid domain and boundary conditions at tank walls, the sloshing problem can be described as: The symbol ∇ is the three-dimensional gradient operator defined as ∇ = (∂/∂x, ∂/∂y, ∂/∂z). r is the position vector, and n is the normal vector at the tank wall, and points tank outside. g is gravitational acceleration and v = v 0 + w × r = [v a , v b , v c ] represents the forced motion of the rectangular tank.
The translatory velocity along (x, y, z) coordinate axes is expressed However, for the convenience of computation, the rotation axes of the angular velocity w = [w 1 , w 2 , w 3 ] are different. The dotted lines in Figure 1 show the rotation axis of the three angular velocity components, which are (x, 0, −h), (0, y, −h), and z-axes, respectively.

Numerical Techniques
A linear damping term for replacing friction effect produced by tank walls is applied. According to [20], the linearized velocity potential Φ n is approximated by Φ n = −a n φ n (x, y, z) sin(ω n t + θ) (5) n represents the sloshing mode, and the parameter is estimated by Here, a n is the amplitudes and ω n the natural frequency. S C and S F represent the tank wall area and linear free surface area, respectively. ν is the kinematic viscosity. The velocity potential is divided into two parts Φ = φ + ϕ, where φ is a disturbance potential and ϕ is a particular solution satisfying the Laplace equation and boundary condition at the tank wall and bottom. According to [20], the particular solution considering the 6-DOF motion of the tank can be constructed as follows: β n cos(γ n y) sinh γ n (x − l/2) γ n cosh(γ n l/2) where λ n = nπ/l, α n = −4l/(n 2 π 2 ), γ n = nπ/b, and β n = −4b/(n 2 π 2 ). Based on the superposition principle, the decomposition of velocity potential leads to a simple boundary condition ∇φ · n = 0, which contributes to solving φ by the Boussinesq-type approach. More details can be found in [19]. With initial conditions known, the time stepping process of Φ and η can be advanced by Equations (3) and (4). The central difference scheme and fourth-order Runge-Kutta method are used for spatial derivative calculation and temporal iteration, respectively.

Validation
The numerical model including the linear damping term has been valid and applied for the shallow-water sloshing in [19][20][21]. Comparison is made between the present results with the mode n of the particular solution equal to 9 and the experimental data given in [22]. According to [22], the time histories of wave elevations were measured by the resistance wave gauges fixed on the tank made of lucite, and a wetting agent was added to the water for reducing the influence of surface tension.
The geometric parameters of the tank are l = 0.6095 m, b = 0.23 m, and h = 0.06 m. The harmonic movement of the tank occurs just in the x-direction, such that X tank = Asin(ωt). The amplitude A is set to 0.196 cm, and excitation frequencies ω are chosen to be around the first natural frequency ω r , which was obtained from the linear solution for sloshing in the tank. As shown in Figure 2, a uniform grid with 31 nodes in the length direction and 11 nodes in the width direction is employed for calculation. The time step is set to 0.06, and the simulation has been over 100 periods to ensure that the sloshing motion reached a stable state. The maximum transient dimensionless amplitudes at the tank wall versus normalized excitation frequencies are compared in Figure 3, where Exp represents the experimental data and Num the numerical results. The excitation frequencies divided by ω r = 3.8924 s −1 ranged from 0.9 to 1.12. The maximum value of the response curve appears around ω/ω r = 1.05 in experiment data while 1.06 in numerical data. The curves rise when approaching the maximum value from the left, but there is a sudden short pullback around ω/ω r = 0.97. Moving away from the maximum value from the right side, the curve drops sharply at first and then tends to decrease slowly. Overall, there are two discontinuities in the response curve. These discontinuities are produced by the effects of nonlinearities which result in a rapid change of sloshing motion form. More details about the sloshing motions can be found in [20].
In general, the results of the numerical model agree well with the experimental results. However, there were inconsistencies in the results from ω/ω r = 1.03 to 1.07. In addition to the occurrence of local breaking on free surface, the sensitivity of sloshing motion to the excitation frequency accounts for it. However, the correct prediction about location of the two discontinuities still shows that our numerical model has a good ability to capture the change of sloshing form caused by the variation of excitation frequency.

Results and Discussion
Based on the former numerical model, the free surface profiles in the rectangular tank are simulated. The sway, roll, coupled sway and roll, coupled surge and sway, coupled surge-sway-roll-pitch, and coupled surge (v 1 )-sway (v 2 )-roll (w 1 )-pitch (w 2 )-yaw (w 3 ) excitations are considered.

The Sway Excitations
A three-dimensional square tank with l = b = 0.5875 m and h = 0.03 m is considered. A uniform grid with 31 nodes in the length direction and 11 nodes in the width direction is employed for calculation. The time step is set to 0.06, and the simulation has been over 100 periods to ensure that the sloshing motion reached a stable state.
The translatory external excitations along the y-axis are adopted here. The excitation amplitude is 0.001 m, and the excitation frequencies vary between 0.92 ω 1 to 1.08 ω 1 . Here, ω 1 = 2.8881 s −1 is the lowest natural frequency. According to the previous research, the violent sloshing motions appear in the vicinity of the lowest natural frequency and the transformation of free surface wave patterns happens. The time-history curves of the relative free surface elevation of the tank wall under different cases are displayed in As shown in Figures 6 and 7, the power spectral densities of time histories of free surface elevations at an observation point within 600 s are calculated. The Matlab function pwelch is used, and the frequencies divided by the excitation frequencies give the values of abscissa. As shown in the figures, the wave energy is distributed in the integer multiple of excitation frequencies. It should be noted that the rule can be applied to all the excitation frequencies from 0.92 ω 1 to 1.08 ω 1 . For the excitation frequencies 0.92 ω 1 and 1.08 ω 1 , the linear standing waves lead to the single frequency energy distribution. For the nonlinear double-bump traveling wave (Figure 5a) and triple-bump traveling wave (Figure 4c), the wave energy can be found on the double, triple, and higher multiple of excitation frequencies, but the main wave energy still focuses on the excitation frequency. In Figure  7b, the wave energy located on the double excitation frequency is higher than that on the excitation frequency, which leads to the nonlinear single-bump traveling waves ( Figure  5b). Different free surface wave patterns in Figures 4 and 5 are due to the superposition of different wave components (Figures 6 and 7). The frequency distribution rule of the wave components has been found, but the amplitude distribution rule should be studied further.

The Roll Excitations
The three-dimensional rectangular tank used in Section 4.1 is adopted here and the following sections. The excitation amplitude is 0.001 rad and frequencies vary between 0.92 ω 1 and 1.08 ω 1 . According to the numerical results, the free surface wave patterns of different excitation frequencies are similar to that in the tank excited by the sway motions. In Figures 8 and 9, the time histories of relative free surface elevations and the power spectral densities with different excitation frequencies 0.96 ω 1 , 0.98 ω 1 and 1.02 ω 1 are given. The power spectral densities of free surface elevations are likewise distributed in the integer multiple of excitation frequencies. In addition, the amplitude distribution rule is similar to that in Section 4.1.

Coupled Sway and Roll Excitations
The amplitudes of sway excitations are 0.001 m and of roll excitations are 0.001 rad. The sway and roll motions are excited at the same frequencies, which ranged from 0.92 ω 1 to 1.08 ω 1 . According to the numerical results, the wave amplitudes excited by the coupled sway and roll motions are much lower than that excited by the single sway motions or single roll motions. Meanwhile, the resonance frequency interval is narrower and closer to the natural frequency ω 1 . The conclusion and details can be found in [20]. For cases excited by the combined motion of sway and roll, the distribution positions of wave energy are also on the integer multiple of excitation frequencies ( Figure 10). As shown in Figure 10c, the wave energy distribution indicates the linear standing waves present in the tank when the coupled excitation frequency equals 1.02 ω 1 . However, the strong nonlinear traveling waves can be observed under cases excited by the single sway motions or single roll motions at the same excitation frequency.

Coupled Surge and Sway Excitations with Different Amplitudes
The surge and sway motions have the same excitation frequency 1.02 ω 1 but different excitation amplitudes. The direction of external combined excitation is not along the diagonal of oxy plane in the tank since the surge motion has a 0.001 m excitation amplitude, while the sway motion has a 0.002 m excitation amplitude. The power spectral densities of free surface elevations at four tank corners (l/2, b/2), (−l/2, b/2), (l/2, −b/2) and (−l/2, −b/2) marked by A, B, C and D are shown in Figures 13 and 14. As shown in the figures, the power spectral density distributions at the corners A and D are the same, and the wave energy locates on the integer multiple of excitation frequency 1.02 ω 1 . The power spectral densities of the other two tank corners B and C are also the same. There is an angle between the direction of combined excitations and the diagonal of oxy plane. The diagonal component of combined excitations is dominant. As shown in Figures 13 and 14, the spectral densities located on the excitation frequency at the corners A and D are much higher than that at the corners B and C.

Coupled Surge-Sway-Roll-Pitch Excitations
In this section, the coupled surge-swing-roll-pitch excitation is researched. The excitation amplitude of the translational motion excitation (surge and sway) is 0.001 m and 0.001 rad of the rotational excitation (roll and pitch). The excitation frequency of the coupled motions is selected to be 0.98 ω 1 . The power spectral density of free surface elevations at the tank corner (l/2, b/2) is shown in Figure 15. The wave energy distribution in the tank can be found on the integer multiple of excitation frequency 0.98 ω 1 . Using the same excitation frequency and amplitude, the superposition of results under coupled sway-roll excitations and coupled surge-pitch excitations equal results of the coupled surge-sway-roll-pitch case. However, the free surface height at the corner of the tank cannot be obtained by the same superposition due to the combined effect of multiple excitations. Compared to the results in Figure 10a, the spectral density due to coupled surge-sway-roll-pitch excitations cannot be calculated by the double spectral density due to coupled sway-roll excitations.

Coupled Surge-Sway-Roll-Pitch-Yaw Excitations
The excitation frequencies of surge, sway, roll, and pitch motions are equal to ω 1 . For the yaw excitations, the excitation frequency is calculated by ω 2 r = gν tanh νh and ν = (π/l) 2 + (π/b) 2 . The excitation amplitude of the translational motion(surge and sway) excitation is 0.001 m and 0.001 rad of the rotational excitation (roll, pitch and yaw). The power spectral densities of free surface elevations at the tank corners (l/2, b/2; −l/2, b/2) are calculated. Because the lowest natural frequencies of yaw excitation ω r = 4.0672 s −1 and others are different, the abscissa is valued by the frequency.
As shown in Figure 16a, the integer multiple of excitation frequency (ω 1 ) is marked by red * , and the integer multiple of excitation frequency (ω r ) is denoted by red o. The sum excitation frequencies ω 1 + ω r and 2ω 1 + ω r are expressed by red +. The wave energy distribution of the free surface elevations at the tank corner A (l/2, b/2) can be found clearly. At the other corner B (−l/2, b/2), the wave energy is focused on the integer multiple of excitation frequency ω r and excitation frequency ω 1 doubling. The power spectral density at the excitation frequency ω 1 can not be found because the corner is not located on the direction of coupled surge-sway-roll-pitch excitations.

Conclusions
The shallow-water sloshing resonance in the three-dimensional tank was researched by Boussinesq-type equations with the linear damping term. The power spectral densities of free surface elevations in tank under resonance were analyzed for different external excitations which include single degree of freedom and multiple degrees of freedom excitations. Research on the sloshing phenomenon of external excitation frequency near the first natural frequency shows that: (1) According to the Figures 7, 9, and 10, whether it is a single degree of freedom excitation or two degrees of freedom excitations with the same excitation frequency, the power spectral density distributions are mainly at the integer multiple of excitation frequencies.
(2) Combining the curves in Figures 11 and 12, when the excitation frequencies of two directions are different, the power spectral density is distributed in the difference excitation frequencies and sum excitation frequencies in addition to the integral multiple excitation frequency. (3) The coupling effect of external excitation on sloshing motion is highly nonlinear. By comparing Figures 15 and 10a, it can be found that the spectral density due to coupled surge-sway-roll-pitch excitations cannot be calculated by the double spectral density due to coupled sway-roll excitations. (4) At present, the power spectral densities of the free surface elevations has been given, and the power spectral densities of the forces should be consistent with it. The wave loads on the tank walls should be further studied.
Understanding the energy spectrum characteristics of sloshing motion can serve to establish a reduced-order model in the future, in which we can focus on these frequencies of the main distribution of sloshing energy. The amplitude distribution rules of power spectral densities should be studied further.
Author Contributions: Conceptualization, Y.S.; validation, X.Y.; formal analysis, X.Y. and P.X.; writing and editing, X.Y. and Y.S. All authors have read and agreed to the published version of the manuscript.