Numerical Simulation of Liquid Sloshing with Different Filling Levels Using OpenFOAM and Experimental Validation

A series of numerical simulations were performed to explore the influences of filling level, excitation frequency and amplitude on liquid sloshing by using the open source Computational Fluid Dynamics toolbox OpenFOAM (Open Field Operation and Manipulation), which was fully validated by the experimental data. The results show that the dynamic impact pressure is proportional to the external excitation amplitude only in non-resonance frequency ranges. Pressure-frequency response curves demonstrate a transition process from a ‘soft-spring’ response to a ‘hard-spring’ response following the changes of the filling level. Such a transition process is found to be dominated by the ratio of the filling level to tank length and the critical value can be obtained. It is also found that wave breaking influences the period of sloshing wave in tanks and ultimately alters the resonance frequency from the linear theory.


Introduction
Sloshing is a phenomenon commonly found in partially-filled liquid storage vessels, such as liquid cargo tanks and fuel tanks, in motion. Violent sloshing may result in serious structural damages to the liquid-tank or even overturn the liquid cargo ship. Thus, a reliable prediction of the sloshing is crucial for the design and deployment of such structures. For this purpose, theoretical analyses, physical experiments, and numerical simulations are commonly utilized.
As a pioneer, Moiseev [1] developed a nonlinear analytical solution for the sloshing problem by using the approximations and modal methods combined with the potential flow theory. Based on Moiseev theory, Faltinsen [2] developed a third-order steady-state solution for the liquid sloshing in a 2D rectangular tank under the swaying and rolling excitations. Faltinsen [3,4] also established nonlinear analytical solutions for liquid sloshing in a rectangular tank by using the multimodal approach. Ikeda and Nakagawa [5] studied liquid sloshing in a rectangular tank under a horizontal excitation using the potential flow theory and analyzed the influence of fluid sloshing on the nonlinear vibration of the structure. However, the theoretical analysis approach may lead to unreliable predictions when complex physical phenomena such as wave breaking and slamming occur due to their fundamental assumption of the potential flow.
Physical experiments are considered an accurate approach and have been widely used to investigate liquid sloshing particularly with a focus on the prediction of the impact pressure. Pistani and Thiagarajan [6] conducted a series of sloshing experiments to accurately measure the high pressure ∂ρ ∂t + ∇ · (ρU) = 0 (1) where ρ is the density, U is the fluid velocity vector, τ is the shear stress, C is the surface tension coefficient which is set to 0 in the current investigation, k is the interface curvature, α is the volume fraction, g is the acceleration of gravity, h is the position vector of the mesh centre measured from the coordinates origin, p rgh is the dynamic pressure. The fluid density ρ and the viscosity coefficient µ are respectively calculated by densities (ρ 1 , ρ 2 ) and viscosity (µ 1 , µ 2 ) of two fluids using the volume fraction α, OpenFOAM uses the Finite Volume Method to discretize its governing equations. A first-order implicit Euler discretization scheme is used for dealing with the time derivative terms, e.g., ∂U/∂t and ∂α/∂t. The Gauss linear discretization scheme is selected for gradient estimation, e.g., ∇·U. Gauss linear corrected is considered for laplacian schemes such as ∇p rgh , ∇ρ. With regard to the divergence terms such as ∇·(UU) and ∇·(Uα), the van Leer scheme is used. The PIMPLE algorithm, which is a combination of PISO (Pressure Implicit with Splitting of Operator) and SIMPLE (Semi-Implicit Method for Pressure-Linked Equations), is used for the velocity-pressure decoupling. More details of the InterDyMFoam can be found in the OpenFOAM website or other references. These will not be repeated here, only the solution process of the InterDyMFoam is illustrated in Figure 1 for completeness. module is the volume of fluid (VOF) based on the phase-fraction. The continuity equation, momentum equation, and phase equation are, respectively, as follows.
where ρ is the density, U is the fluid velocity vector, τ is the shear stress, C is the surface tension coefficient which is set to 0 in the current investigation, к is the interface curvature, α is the volume fraction, g is the acceleration of gravity, h is the position vector of the mesh centre measured from the coordinates origin, prgh is the dynamic pressure. The fluid density ρ and the viscosity coefficient μ are respectively calculated by densities (ρ1, ρ2) and viscosity (μ1, μ2) of two fluids using the volume fraction α, OpenFOAM uses the Finite Volume Method to discretize its governing equations. A first-order implicit Euler discretization scheme is used for dealing with the time derivative terms, e.g., ∂U/∂t and ∂α/∂t. The Gauss linear discretization scheme is selected for gradient estimation, e.g., ∇U. Gauss linear corrected is considered for laplacian schemes such as ∇prgh, ∇ρ. With regard to the divergence terms such as ∇·(UU) and ∇(Uα), the van Leer scheme is used. The PIMPLE algorithm, which is a combination of PISO (Pressure Implicit with Splitting of Operator) and SIMPLE (Semi-Implicit Method for Pressure-Linked Equations), is used for the velocity-pressure decoupling. More details of the InterDyMFoam can be found in the OpenFOAM website or other references. These will not be repeated here, only the solution process of the InterDyMFoam is illustrated in Figure 1 for completeness.

Experimental Setup
For the purpose of validating the numerical model, as well as exploring detailed physics associated with the liquid sloshing, a series of experiments are carried out in the Laboratory of Vibration Test and Liquid Sloshing at the Hohai University, China. A six-degree of freedom (DoF) motion simulation platform, which is commonly known as a hexapod and is able to perform six-DoF motions regularly or randomly according to an appropriate input of time histories, is utilized to generate the forced motions of the liquid tank, as shown in Figure 2. The tank used in the experiments is a rectangular tank made of plexiglass with an 8-mm thickness. The internal dimensions of the rectangular tank are L = 600 mm in length, W = 300 mm in width and H = 650 mm in height. A horizontal sinusoidal motion t A x ω sin -= is assigned to the rectangular tank. The amplitude and frequency of the motion can be set up through a control computer. The displacement of the motion platform was recorded by a displacement sensor. Figure 3, which compares both the actual motion platform displacement and theoretical displacement, indicates that the experimental apparatus has high motion accuracy.

Experimental Setup
For the purpose of validating the numerical model, as well as exploring detailed physics associated with the liquid sloshing, a series of experiments are carried out in the Laboratory of Vibration Test and Liquid Sloshing at the Hohai University, China. A six-degree of freedom (DoF) motion simulation platform, which is commonly known as a hexapod and is able to perform six-DoF motions regularly or randomly according to an appropriate input of time histories, is utilized to generate the forced motions of the liquid tank, as shown in Figure 2. The tank used in the experiments is a rectangular tank made of plexiglass with an 8-mm thickness. The internal dimensions of the rectangular tank are L = 600 mm in length, W = 300 mm in width and H = 650 mm in height. A horizontal sinusoidal motion x = −A sin ωt is assigned to the rectangular tank. The amplitude and frequency of the motion can be set up through a control computer. The displacement of the motion platform was recorded by a displacement sensor. Figure 3, which compares both the actual motion platform displacement and theoretical displacement, indicates that the experimental apparatus has high motion accuracy.  In the experimental study, two filling levels, 13.8% and 30.8%, are considered. For measuring the dynamic impact pressure on the tank wall, five pressure sensors are embedded in the left wall, as shown in Figure 2. Their specific locations are illustrated in Figure 4. A camera was fixed in front of the rectangular tank to record the profile of the free surface.
With the assumption of potential flow theory, the natural frequencies can be analytically determined as [26] ) h tanh( L n L n g n π π ω = , n = 1, 2, 3, where n is the mode number and h is the filling level. According to Equation (5), the first-mode natural frequencies are ω1 = 4.749 rad/s and ω1 = 6.333 rad/s for the tank with the low (13.8%) and high (30.8%) filling level, respectively.  In the experimental study, two filling levels, 13.8% and 30.8%, are considered. For measuring the dynamic impact pressure on the tank wall, five pressure sensors are embedded in the left wall, as shown in Figure 2. Their specific locations are illustrated in Figure 4. A camera was fixed in front of the rectangular tank to record the profile of the free surface.
With the assumption of potential flow theory, the natural frequencies can be analytically determined as [26] ) h tanh( L n L n g n π π ω = , n = 1, 2, 3, where n is the mode number and h is the filling level. According to Equation (5), the first-mode natural frequencies are ω1 = 4.749 rad/s and ω1 = 6.333 rad/s for the tank with the low (13.8%) and high (30.8%) filling level, respectively. In the experimental study, two filling levels, 13.8% and 30.8%, are considered. For measuring the dynamic impact pressure on the tank wall, five pressure sensors are embedded in the left wall, as shown in Figure 2. Their specific locations are illustrated in Figure 4. A camera was fixed in front of the rectangular tank to record the profile of the free surface.
With the assumption of potential flow theory, the natural frequencies can be analytically determined as [26] where n is the mode number and h is the filling level. According to Equation (5), the first-mode natural frequencies are ω 1 = 4.749 rad/s and ω 1 = 6.333 rad/s for the tank with the low (13.8%) and high (30.8%) filling level, respectively.

Mesh Convergence Test
Mesh convergence studies are performed to determine the mesh size used in the following simulations. The computational mesh consists of a uniform square grid. Three mesh sizes are considered in the tests. They are 2 mm, 5 mm and 8 mm, respectively. The model tank size is exactly the same as the experimental tank and the water depth is 90 mm (i.e., the low filling level in the experiment). The test parameters are chosen as A = 7 mm and ω = 4.749 rad/s. The frequency is the same as the first mode natural frequency predicted using the linear theory (Equation (5)). Therefore, we expected to see a violent sloshing in the tank due to the resonance. If the mesh resolution leads to a convergent solution for this case, it shall be sufficient for other cases, especially non-resonance cases. The comparison of the pressure recorded at P1 in the cases with different mesh sizes is plotted in Figure 5. These results are all similar, although the result obtained by using the 8-mm mesh slightly differs from others. A mesh size of 5 mm is employed in the following studies. Furthermore, the time step size is automatically adjusted by the condition that the Courant number ≤0.5.

Experimental Validation
To validate the numerical model for the cases with different sloshing conditions, six sets of experimental data with two different filling levels and three different excitation frequencies are employed. The comparisons between the experimental data and the numerical results are illustrated

Mesh Convergence Test
Mesh convergence studies are performed to determine the mesh size used in the following simulations. The computational mesh consists of a uniform square grid. Three mesh sizes are considered in the tests. They are 2 mm, 5 mm and 8 mm, respectively. The model tank size is exactly the same as the experimental tank and the water depth is 90 mm (i.e., the low filling level in the experiment). The test parameters are chosen as A = 7 mm and ω = 4.749 rad/s. The frequency is the same as the first mode natural frequency predicted using the linear theory (Equation (5)). Therefore, we expected to see a violent sloshing in the tank due to the resonance. If the mesh resolution leads to a convergent solution for this case, it shall be sufficient for other cases, especially non-resonance cases. The comparison of the pressure recorded at P1 in the cases with different mesh sizes is plotted in Figure 5. These results are all similar, although the result obtained by using the 8-mm mesh slightly differs from others. A mesh size of 5 mm is employed in the following studies. Furthermore, the time step size is automatically adjusted by the condition that the Courant number ≤0.5.

Mesh Convergence Test
Mesh convergence studies are performed to determine the mesh size used in the following simulations. The computational mesh consists of a uniform square grid. Three mesh sizes are considered in the tests. They are 2 mm, 5 mm and 8 mm, respectively. The model tank size is exactly the same as the experimental tank and the water depth is 90 mm (i.e., the low filling level in the experiment). The test parameters are chosen as A = 7 mm and ω = 4.749 rad/s. The frequency is the same as the first mode natural frequency predicted using the linear theory (Equation (5)). Therefore, we expected to see a violent sloshing in the tank due to the resonance. If the mesh resolution leads to a convergent solution for this case, it shall be sufficient for other cases, especially non-resonance cases. The comparison of the pressure recorded at P1 in the cases with different mesh sizes is plotted in Figure 5. These results are all similar, although the result obtained by using the 8-mm mesh slightly differs from others. A mesh size of 5 mm is employed in the following studies. Furthermore, the time step size is automatically adjusted by the condition that the Courant number ≤0.5.

Experimental Validation
To validate the numerical model for the cases with different sloshing conditions, six sets of experimental data with two different filling levels and three different excitation frequencies are employed. The comparisons between the experimental data and the numerical results are illustrated

Experimental Validation
To validate the numerical model for the cases with different sloshing conditions, six sets of experimental data with two different filling levels and three different excitation frequencies are employed. The comparisons between the experimental data and the numerical results are illustrated in Figure 6, in which the red circles represent the experimental data and the black line represents the numerical results. Overall, the numerical results agree with the experimental data for the It is noticed that the discrepancies between the numerical results and experimental results in the cases shown in Figure 6a,b become relatively more significant after t = 15 s. Further analysis on the pressure spectra ( Figure 7) suggests that the difference in Figure 6a is mainly caused by the fundamental pressure component (i.e., ω 1 ) and the experimental data decays after t = 15 s. For the case shown in Figure 6b, the spectra shown in Figure 7b further confirm a satisfactory agreement. In addition, the free surface profiles are also compared. Some results in the case with a resonance excitation are shown in Figure 8, which compares the wave profiles at six instants between 8.55 s and 9.17 s. A good match has been observed.
in Figure 6, in which the red circles represent the experimental data and the black line represents the numerical results. Overall, the numerical results agree with the experimental data for the non-resonance (Figure 6a-d), near-resonance (Figure 6b-e), or fully-resonance (Figure 6c-f) conditions. It is noticed that the discrepancies between the numerical results and experimental results in the cases shown in Figure 6a,b become relatively more significant after t = 15 s. Further analysis on the pressure spectra ( Figure 7) suggests that the difference in Figure 6a is mainly caused by the fundamental pressure component (i.e., ω1) and the experimental data decays after t = 15 s. For the case shown in Figure 6b, the spectra shown in Figure 7b further confirm a satisfactory agreement. In addition, the free surface profiles are also compared. Some results in the case with a resonance excitation are shown in Figure 8, which compares the wave profiles at six instants between 8.55 s and 9.17 s. A good match has been observed.   in Figure 6, in which the red circles represent the experimental data and the black line represents the numerical results. Overall, the numerical results agree with the experimental data for the non-resonance (Figure 6a-d), near-resonance (Figure 6b-e), or fully-resonance (Figure 6c-f) conditions. It is noticed that the discrepancies between the numerical results and experimental results in the cases shown in Figure 6a,b become relatively more significant after t = 15 s. Further analysis on the pressure spectra ( Figure 7) suggests that the difference in Figure 6a is mainly caused by the fundamental pressure component (i.e., ω1) and the experimental data decays after t = 15 s. For the case shown in Figure 6b, the spectra shown in Figure 7b further confirm a satisfactory agreement. In addition, the free surface profiles are also compared. Some results in the case with a resonance excitation are shown in Figure 8, which compares the wave profiles at six instants between 8.55 s and 9.17 s. A good match has been observed.     For further validation, the experimental data in the literature are also used. The first case considered here is the 2D sloshing experiment carried out by Liu and Lin [12]. In this case, the tank dimensions are L = 570 mm and H = 300 mm. The water depth h is 150 mm. The tank undergoes a sinusoidal translation with an amplitude of 5 mm and the frequency of 6.0578 rad/s. The comparison of the surface elevation recorded near the left boundary (20 mm) between the numerical and experimental data is shown in Figure 9. Excellent agreement is observed. Another experiment was performed by Liu et al [27], in which the 2D tank had a length of L = 900 mm and a height of H = 508 mm. The filling level h/H is 18.3%. The tank undergoes a sinusoidal rolling motion with motion amplitude of 4 degrees. The pressure recorded on the wall and at the mean water surface is used for comparison. The results are shown in Figure 10. Through such comparisons, one may agree that the present OpenFOAM modelling results in a satisfactory accuracy and is suitable for numerically simulating the sloshing problems.  For further validation, the experimental data in the literature are also used. The first case considered here is the 2D sloshing experiment carried out by Liu and Lin [12]. In this case, the tank dimensions are L = 570 mm and H = 300 mm. The water depth h is 150 mm. The tank undergoes a sinusoidal translation with an amplitude of 5 mm and the frequency of 6.0578 rad/s. The comparison of the surface elevation recorded near the left boundary (20 mm) between the numerical and experimental data is shown in Figure 9. Excellent agreement is observed. Another experiment was performed by Liu et al. [27], in which the 2D tank had a length of L = 900 mm and a height of H = 508 mm. The filling level h/H is 18.3%. The tank undergoes a sinusoidal rolling motion with motion amplitude of 4 degrees. The pressure recorded on the wall and at the mean water surface is used for comparison. The results are shown in Figure 10. Through such comparisons, one may agree that the present OpenFOAM modelling results in a satisfactory accuracy and is suitable for numerically simulating the sloshing problems.  For further validation, the experimental data in the literature are also used. The first case considered here is the 2D sloshing experiment carried out by Liu and Lin [12]. In this case, the tank dimensions are L = 570 mm and H = 300 mm. The water depth h is 150 mm. The tank undergoes a sinusoidal translation with an amplitude of 5 mm and the frequency of 6.0578 rad/s. The comparison of the surface elevation recorded near the left boundary (20 mm) between the numerical and experimental data is shown in Figure 9. Excellent agreement is observed. Another experiment was performed by Liu et al [27], in which the 2D tank had a length of L = 900 mm and a height of H = 508 mm. The filling level h/H is 18.3%. The tank undergoes a sinusoidal rolling motion with motion amplitude of 4 degrees. The pressure recorded on the wall and at the mean water surface is used for comparison. The results are shown in Figure 10. Through such comparisons, one may agree that the present OpenFOAM modelling results in a satisfactory accuracy and is suitable for numerically simulating the sloshing problems.

Three-Dimensional Results
It is widely accepted that the single degree of freedom sloshing shown in Figures 7 and 8 can be simplified as 2D problems in the numerical simulation. However, in the cases under resonance conditions, the wave breaking occurs and the 3D effect may become important. To investigate the significance of the 3D effect, corresponding 3D numerical investigations are carried out. The numerical results by both the 2D model and the 3D model are compared with the experimental data to shed light on the 3D effect. In the 3D simulations, pressure probes are placed across the transverse direction (shown in Figure 11). Two cases subjected to resonant conditions are considered. These   Figure 12 compares the pressure--time histories at the probes located at the same height from the tank bottom in 3D numerical simulations but different transverse positions. As observed, the impact pressures at different transverse positions but the same height are largely the same, except the pressures at the corner of the tanks, e.g., P1, P4, and P7. Further comparison in Figure 13 shows that the 3D model leads to a slightly better prediction on the pressures on the corner points, compared to the 2D simulations. However, the superiority of the 3D model over the 2D model is not obvious in terms of the improvement of the computational accuracy. Considering the high

Three-Dimensional Results
It is widely accepted that the single degree of freedom sloshing shown in Figures 7 and 8 can be simplified as 2D problems in the numerical simulation. However, in the cases under resonance conditions, the wave breaking occurs and the 3D effect may become important. To investigate the significance of the 3D effect, corresponding 3D numerical investigations are carried out. The numerical results by both the 2D model and the 3D model are compared with the experimental data to shed light on the 3D effect. In the 3D simulations, pressure probes are placed across the transverse direction (shown in Figure 11). Two cases subjected to resonant conditions are considered. These include

Three-Dimensional Results
It is widely accepted that the single degree of freedom sloshing shown in Figures 7 and 8 can be simplified as 2D problems in the numerical simulation. However, in the cases under resonance conditions, the wave breaking occurs and the 3D effect may become important. To investigate the significance of the 3D effect, corresponding 3D numerical investigations are carried out. The numerical results by both the 2D model and the 3D model are compared with the experimental data to shed light on the 3D effect. In the 3D simulations, pressure probes are placed across the transverse direction (shown in Figure 11). Two cases subjected to resonant conditions are considered. These   Figure 12 compares the pressure--time histories at the probes located at the same height from the tank bottom in 3D numerical simulations but different transverse positions. As observed, the impact pressures at different transverse positions but the same height are largely the same, except the pressures at the corner of the tanks, e.g., P1, P4, and P7. Further comparison in Figure 13 shows that the 3D model leads to a slightly better prediction on the pressures on the corner points, compared to the 2D simulations. However, the superiority of the 3D model over the 2D model is not obvious in terms of the improvement of the computational accuracy. Considering the high Figure 11. The locations of pressure probes. Figure 12 compares the pressure-time histories at the probes located at the same height from the tank bottom in 3D numerical simulations but different transverse positions. As observed, the impact pressures at different transverse positions but the same height are largely the same, except the pressures at the corner of the tanks, e.g., P1, P4, and P7. Further comparison in Figure 13 shows that the 3D model leads to a slightly better prediction on the pressures on the corner points, compared to the 2D simulations. However, the superiority of the 3D model over the 2D model is not obvious in terms of the improvement of the computational accuracy. Considering the high computational demand of the 3D simulation, as demonstrated in Table 1, one may agree that the 2D model is preferable in terms of computational robustness for the problems concerned in this paper. computational demand of the 3D simulation, as demonstrated in Table 1, one may agree that the 2D model is preferable in terms of computational robustness for the problems concerned in this paper.  computational demand of the 3D simulation, as demonstrated in Table 1, one may agree that the 2D model is preferable in terms of computational robustness for the problems concerned in this paper.

Effects of External Excitation Amplitudes
The excitation amplitude is an important factor that influences the intensity of the sloshing in the forced sloshing problem, and it reflects the input energy of the system. In this study, the influence of amplitude on the sloshing under different conditions will be discussed after a series of numerical studies being conducted. These cases are summarized in Table 2.  Figure 14 displays the response curves of the maximum impact pressure vs. the amplitude of the forced motion at different combinations of the filling level and frequency. It is worth mentioning that the vertical coordinate represents the peak of impact pressure which, at P1, was non-dimensionalized based on water depth, density and gravity (P = P 0 /ρgh, h is set as 0.2 m). Data obtained through physical model experiments (the amplitudes are 3 mm, 5 mm and 7 mm respectively) are also added to Figure 14 for validation. It can be noted that the impact pressure increases as the amplitude increases and that increase is basically linear in the non-resonance range. In order to further illustrate this relationship, Figure 15 shows the pressure-time histories in the cases with different excitation amplitudes. For the convenience of the comparison, the pressure results in the cases with A = 1 mm multiplied by a scaling factor of 3 is able to be comparable with the corresponding results with A = 3 mm. If two curves in Figure 15 match with each other (e.g., Figure 15a-d), the pressure is then linearly dependent on the motion amplitude; otherwise, the nonlinearity becomes more important (e.g., Figure 15e-f). This condition often occurs following the occurrence of wave resonance and breaking.

Resonant Hysteresis and Resonance in Advance
In the sloshing problem, the natural frequency of the tank can be computed using the analytical solution based on the potential flow theory Equation (5). However, the frequency leading to the maximum impact pressure is generally not consistent with the natural frequency due to the nonlinearity. The theory [28] has indicated the change in the resonant response from a so-called 'hard-spring' to 'soft-spring' behaviour as the filling level increases. Kobine's [29] confirms the predicted change in the hysteresis behaviour with an increasing filling level. Based on these studies, this section mainly conducts studies to explore the mechanism of the resonant hysteresis.

Resonant Hysteresis and Resonance in Advance
In the sloshing problem, the natural frequency of the tank can be computed using the analytical solution based on the potential flow theory Equation (5). However, the frequency leading to the maximum impact pressure is generally not consistent with the natural frequency due to the nonlinearity. The theory [28] has indicated the change in the resonant response from a so-called 'hard-spring' to 'soft-spring' behaviour as the filling level increases. Kobine's [29] confirms the predicted change in the hysteresis behaviour with an increasing filling level. Based on these studies, this section mainly conducts studies to explore the mechanism of the resonant hysteresis.

Resonant Hysteresis and Resonance in Advance
In the sloshing problem, the natural frequency of the tank can be computed using the analytical solution based on the potential flow theory Equation (5). However, the frequency leading to the maximum impact pressure is generally not consistent with the natural frequency due to the nonlinearity. The theory [28] has indicated the change in the resonant response from a so-called 'hard-spring' to 'soft-spring' behaviour as the filling level increases. Kobine's [29] confirms the predicted change in the hysteresis behaviour with an increasing filling level. Based on these studies, this section mainly conducts studies to explore the mechanism of the resonant hysteresis.
A large number of numerical simulations are performed to explore the influence of the filling level on the response between the impact pressure of tank wall and external excitation frequency. The specific parameters of the cases are in Table 3. There are 10 filling levels, one amplitude of motion, 15 external excitation frequencies, 30 s of sloshing duration considered. In order to gather more stable and accurate pressure data, the pressure probe is set at the bottom of the tank wall which is away from the pressure impacts (h = 30 mm).  Figure 16 shows the pressure-frequency response curves at 10 different filling levels. In Figure 16c-h, the corresponding experimental are also included for comparison, which, once again, show good agreements and reinforce the credibility of experimental and numerical models. In all the cases, the pressure increases as the frequency increases and then decreases afterword. However, at different filling levels, the rate of increase/decrease and the maximum response frequency are not the same. It can be clearly seen that when the filling level changes from low to high, the ascending process before the maximum response frequency becomes shorter, and the descending process after that becomes longer, the maximum response frequency (Marked in the Figure 16) decreases from 1.18 ω 1 to 0.92 ω 1 instead of keeping the first-mode natural frequency computed by the potential theory.
As is shown in Figure 16a-c, there are 'soft-spring' responses in shallow water sloshing; it is characterized by a rapid decrease after a slow rise to the maximum response value, with the maximum response frequency higher than the natural frequency. Besides, Figure 16h-k show a 'hard-spring' response in deepwater sloshing. Its characteristic is a slow descent after a rapid rise to the maximum response value and maximum response frequency is less than natural frequency. There is a critical depth, as shown in Figure 16e, at which level the resonance curve becomes symmetric about the natural frequency, being consistent with the resonance in normal harmonic systems.
The occurrences of the 'soft-spring' response and the 'hard-spring' response may be mainly caused by the nonlinear effect of wave breaking. In order to study the effect of wave breaking, unbroken cases of shallow water and deepwater are set up. With all other conditions remaining the same, the amplitude was reduced to 0.5 mm to make sure that the free surface is not broken during the whole process of sloshing. Figures 17 and 18 show the pressure-frequency response curve, the mean squared error-frequency response curve, the pressure-time history curve, and the spectral analysis at the maximum response frequency in the case of breaking and non-breaking. It should be noted that when at the same depth but without wave breaking, the phenomenon of resonant hysteresis disappears, and the analysis of the max pressure and mean squared error show a similar trend. It is well known that the double peak phenomenon in the pressure-time history curve under resonance indicates two impacts on the tank wall after the sloshing wave is broken. It can be seen in Figure 17 that the resonant sloshing at the low filling level has broken in the third period of sloshing, far earlier than the resonant sloshing at the high filling level. The Fourier transform in Figures 17 and 18 shows that more multiplications of the natural frequency appear in the wave breaking case and the case at the low filling level, which also shows that the phenomenon of wave breaking, will produce nonlinear effects and the sloshing at the low filling level has stronger nonlinearity than the sloshing at a high filling level.
produce nonlinear effects and the sloshing at the low filling level has stronger nonlinearity than the sloshing at a high filling level.  Based on the previous section's conclusion that the external excitation amplitude has a linear effect on the sloshing pressure, the pressure of the A = 0.5 mm case has been expanded 14 times to compare with the pressure of the A = 7 mm case to investigate how the sloshing state will be if the wave cannot break. As shown in Figure 19, the two curves all coincide before wave breaking, but the difference happened after the waves break. A significant phase advance occurred in the pressure-time history curve of shallow water sloshing, and phase delay appears in the pressure-time history curve of shallow water sloshing. This result shows that wave breaking will change the period of liquid sloshing in the tank. In the condition of shallow water, the period of the liquid sloshing decrease with wave breaking, the maximum response can be achieved at an external excitation frequency greater than the natural frequency. As for deepwater sloshing, the period of the liquid sloshing goes up with the wave breaking, so the system resonance can be excited at an external excitation frequency lower than the natural frequency. These are the root causes of resonant hysteresis and resonance in advance.
Why does the wave breaking have different effects on shallow water and deepwater sloshing? The hydraulic jump is a common phenomenon in shallow water sloshing; Figure 20 shows the wave propagation process before and after the occurrence of the hydraulic jump in the resonant case. It should be noted that the process of a sloshing wave propagating from the highest point on the right wall to the highest point on the left wall takes 0.71 s before the occurrence of the hydraulic jump, but it becomes 0.67 s after the hydraulic jump happens. One reason for this phenomenon is that the liquid climbing up the wall falls vertically and hits the free surface as shown in Figure 20; this force accelerates the velocity at the free surface of the sloshing wave so that the period of the liquid sloshing decreases. The wavelength of deepwater sloshing is longer, a standing wave phenomenon occurs under the excitation of resonance. Figure 21 shows the complex nonlinear phenomenon in deepwater resonant sloshing. It can be seen that there are obvious roof slamming phenomena and aerification which may dissipate energy during each impact. Besides, a large amount of liquid splashing will reduce the depth of water, which will lead to lower wave speed and a longer period of deepwater sloshing.
a. u. Based on the previous section's conclusion that the external excitation amplitude has a linear effect on the sloshing pressure, the pressure of the A = 0.5 mm case has been expanded 14 times to compare with the pressure of the A = 7 mm case to investigate how the sloshing state will be if the wave cannot break. As shown in Figure 19, the two curves all coincide before wave breaking, but the difference happened after the waves break. A significant phase advance occurred in the pressure--time history curve of shallow water sloshing, and phase delay appears in the pressure--time history curve of shallow water sloshing. This result shows that wave breaking will change the period of liquid sloshing in the tank. In the condition of shallow water, the period of the liquid sloshing decrease with wave breaking, the maximum response can be achieved at an external excitation frequency greater than the natural frequency. As for deepwater sloshing, the period of the liquid sloshing goes up with the wave breaking, so the system resonance can be excited at an external excitation frequency lower than the natural frequency. These are the root causes of resonant hysteresis and resonance in advance.
Why does the wave breaking have different effects on shallow water and deepwater sloshing? The hydraulic jump is a common phenomenon in shallow water sloshing; Figure 20 shows the wave propagation process before and after the occurrence of the hydraulic jump in the resonant case. It should be noted that the process of a sloshing wave propagating from the highest point on the right wall to the highest point on the left wall takes 0.71 s before the occurrence of the hydraulic jump, but it becomes 0.67 s after the hydraulic jump happens. One reason for this phenomenon is that the liquid climbing up the wall falls vertically and hits the free surface as shown in Figure 20; this force accelerates the velocity at the free surface of the sloshing wave so that the period of the liquid sloshing decreases. The wavelength of deepwater sloshing is longer, a standing wave phenomenon occurs under the excitation of resonance. Figure 21 shows the complex nonlinear phenomenon in deepwater resonant sloshing. It can be seen that there are obvious roof slamming phenomena and aerification which may dissipate energy during each impact. Besides, a large amount of liquid

Conclusions
A liquid sloshing numerical model was developed based on OpenFOAM in this study, and a series of experiments were conducted to validate the accuracy of the model. The satisfactory

Conclusions
A liquid sloshing numerical model was developed based on OpenFOAM in this study, and a series of experiments were conducted to validate the accuracy of the model. The satisfactory agreement between experimental and numerical data gives support to the assumption that the model

Conclusions
A liquid sloshing numerical model was developed based on OpenFOAM in this study, and a series of experiments were conducted to validate the accuracy of the model. The satisfactory agreement between experimental and numerical data gives support to the assumption that the model is suitable for 2D tank sloshing problems. The numerical results of the 2D model and the 3D model were compared with the experimental results, suggesting an insignificant 3D effect associated with the cases concerned in this paper. Considering the extremely longer Central Processing Unit (CPU) time by the 3D model, the 2D model is preferable and is used in the systematic numerical investigations. A large number of numerical simulations were performed to explore the influence of frequency, amplitude, and the filling level on sloshing: (1) the research on the influence of amplitude indicated that the impact pressure increases as the amplitude grows, and the increase is basically linear in the non-resonance range. Moreover, the change in amplitude affects sloshing more greatly near the resonance frequency. However, under the excitation of resonance frequency, wave breaking may challenge such a principle because of nonlinear effect; (2) the influence of water depth on pressure-frequency response also has been studied, the process of the transition from a 'soft-spring' response to a 'hard-spring' response is observed following the change of the filling level, the ascending trend of the pressures before the maximum response frequency becomes shorter, and the descending process after that becomes longer. The maximum response frequency decreases from 1.18 ω 1 to 0.92 ω 1 moving away from the first-mode natural frequency ω 1 calculated by potential theory; (3) The main reason for this difference should be the nonlinear effect of wave breaking. The nonlinear effect of wave breaking was studied by a pressure-frequency response curve, the mean squared error-frequency response curve, the pressure-time history curve, and the spectral analysis at the maximum response frequency in the case of breaking and non-breaking. The result show that when at the same depth but without wave breaking, the phenomenon of resonant hysteresis disappears, and the resonance sloshing at a low filling level is more nonlinear than the resonance sloshing at a high filling level; (4) The pressure of the A = 0.5 mm case has been expanded 14 times to compare with the pressure of the A = 7 mm case to investigate how the sloshing state will be if wave cannot break. The result shows that the period of the shallow liquid sloshing decreases with wave breaking, therefore, the maximum response can be achieved at an external excitation frequency higher than the natural frequency. Combined with the fluid phenomenon in Figures 20 and 21, the reasons may be the liquid climbing up the wall falls vertically and hits the free surface; this force accelerates the velocity of the sloshing wave. In the condition of deep water, the period of the liquid sloshing goes up with wave breaking, so the system resonance can be excited at an external excitation frequency lower than the natural frequency. One reasonable reason is that there is an obvious roof slamming and aerification which may dissipate energy during each slamming. Besides, a large amount of liquid splashing will reduce the depth of the water, which will lead to lower wave speed and a longer period of deepwater sloshing.
Some reasonable reasons have been raised combined with the fluid phenomenon in Figures 20  and 21: (1) the drop of the fluid after the hydraulic jump accelerates the velocity of the sloshing wave; (2) the roof slamming and aerification may dissipate energy during each impact; (3) a large amount of liquid splashing will reduce the depth of water, which will lead to a lower wave speed and a longer period of deepwater sloshing.