Numerical Simulation of Sloshing in 2D Rectangular Tanks Based on the Prediction of Free Surface

A finite difference method for analyzing 2D nonlinear sloshing waves in a tank has been developed based on the potential flow theory. After σ-transformation, the free surface is predicted by the kinematic condition, and nonlinear terms are approximated; the governing equation and boundary conditions are discretized to linear equations in the iterative process of time. Simulations of standing waves and sloshing in horizontally excited tanks are presented. The results are compared with analytical and numerical solutions in other literatures, which demonstrate the effectiveness and accuracy of this numericalmethod.Thebeating phenomenon of sloshing in the tank with different aspect ratios is studied. The relationship between sloshing force and aspect ratio under the same external excitation is also discussed.


Introduction
Sloshing is a liquid vibration phenomenon caused by the movement of the tank.When the liquid cargo is in transit, the sloshing would affect stability of the system severely, leading to damage or fatigue of the structure.So it is necessary to diminish the impact of sloshing and avoid large amplitude resonance.Sloshing has been studied for many years by analytical, numerical, and experimental methods.In many early studies, the analytical method was dominant.Abramson [1] used potential theory to study the linear sloshing of small amplitude.His work was systematic and influential.Considering the importance of nonlinear effects in the sloshing response, Faltinsen [2] analyzed nonlinear sloshing by perturbation theory.Later, numerical simulation of large amplitude sloshing has been a hot topic in research.Complicated 2D or 3D models have been calculated, and fully nonlinear theory has been developed to get the accurate solutions of sloshing.Faltinsen [3] presented a nonlinear numerical method of 2D sloshing in tanks.Nakayama and Washizu [4] adopted the boundary element method to analyze the sloshing in a tank which is subjected to pitching oscillation.Wu et al. [5] calculated sloshing waves in a 3D tank by using a finite element method based on fully nonlinear potential theory.Chen and Nokes [6] simulated complete 2D sloshing motion by a finite difference method.The primitive 2D Navier-Stokes equations are solved, and both the nonlinear free surface condition and fluid viscosity are considered.Afterwards, Wu and Chen [7] simulated fluid sloshing in a 3D tank with the six degrees of freedom by this timeindependent finite difference method.Kim [8] employed SOLA scheme to solve Navier-Stokes equations and assumed the free surface profile to be a single-valued function.He used this method to simulate sloshing flows in 2D and 3D containers.Celebi and Akyildiz [9] calculated nonlinear viscous liquid sloshing in a partially filled rectangular tank.They used finite difference method to solve equations, and the free surface was captured by VOF technique.Wang and Khoo [10] considered 2D nonlinear sloshing under random excitation by the finite element method.Sriram et al. [11] studied random sloshing under both horizontal and vertical excitation.In the aspect of theory, Faltinsen et al. [12,13] developed the infinite dimensional modal analysis technique to describe nonlinear sloshing of incompressible fluid, where the free surface motion was expanded in generalized Fourier series; Wu [14] analyzed second-order resonance of sloshing.Their conclusions were discussed in many other research articles (Hill [15], Frandsen [16], Firouz-Abadi et al. [17], etc.).
It is difficult to solve the sloshing equations with nonlinear free surface boundary conditions because free surface is always moving in the sloshing process.Some special methods are used to treat the moving boundary of free surface, such as VOF and ALE.Recently, -transformation has been employed in sloshing problems.It is an efficient space transformation method, especially suitable for equations with curved boundary.The basic idea of -transformation is to introduce a new stretching variable in the vertical direction, which transforms the computational domain into a regular shape.Chern et al. [18] and Frandsen et al. [16,19] used transformation to map the liquid domain onto rectangular region and then calculated horizontally and vertically forced sloshing problems in tanks.Turnbull et al. [20] simulated several 2D free surface wave motions by -transformed finite element model.According to Frandsen's analysis, transformation has two major advantages in simulating free surface flows: remeshing due to the moving free surface is avoided and free surface smoothing is often not required.The disadvantage of it is that it is a unique transformation and restricts the wave shape to be nonoverturning.In this paper, a new numerical method is proposed based on transformation.In the iterative process of time, a forecast scheme is used to estimate the free surface boundary, and some nonlinear terms are approximated.Then the governing equation and boundary conditions are linearized and solved by finite difference method.In many numerical methods, free surface is often updated directly by explicit schemes, such as Adams-Bashforth or Runge Kutta methods, while in this method, a predictor-corrector scheme of free surface is used.The free surface is first predicted before the iteration, so the fluid boundary is more precise for the numerical simulation; after the governing equation and boundary conditions are solved, the free surface is updated by the solutions.Because of the predictive step of free surface, it is reasonable to consider high accuracy is one of the advantages of this method.Furthermore, semi-implicit Crank-Nicolsen scheme is employed to discretize linear terms in free surface boundary conditions, so this method has better numerical stability than explicit schemes.Nonlinear sloshing in fixed and horizontally excited rectangular tanks has been simulated by this numerical method.The results are compared with analytical solutions or other numerical results, which prove the efficiency of this method.The limitation of this numerical method is also described.
Many researches have focused on sloshing phenomenon in a fixed aspect ratio by different excitations.In fact, the aspect ratio is always changed in the process of storage or transportation of liquid cargo.So it is necessary to investigate the influence of aspect ratio on liquid sloshing.Chen and Chiang [21] analyzed the effect of fluid depth on nonlinear characteristics of sloshing.In this paper, beating phenomenon of sloshing in different aspect ratios is considered.The shape of the tank is unchangeable, but the depths of liquid are set to be many different values.Under the same external excitation, the vibrations of free surface and sloshing forces acting on the walls are calculated and analyzed.The natural frequency of sloshing is affected by different aspect ratios, which may change the whole sloshing situation.The effect of nonlinear factor on sloshing force is also discussed.

Mathematical Model
The fluid is assumed to be incompressible, irrotational, and inviscid water, so the potential flow theory is employed to describe the sloshing in tanks.Surface tension is neglected.The free surface is assumed to never become overturned or broken during the sloshing process. and  are the length and the width of the rectangular tank, respectively.ℎ is the still water depth.Two Cartesian coordinate systems are introduced in Figure 1.One is the inertial Cartesian coordinate system (, ) fixed in space; -axis is horizontal and -axis is vertical.The other is the moving coordinate system (, ) connected to the tank, with the origin at the intersection of undisturbed free surface and the left wall of the tank.The tank is expected to have a translational oscillation along the -axis, so there is no flow along the width direction of the tank, which means the sloshing of liquid is only 2D nonlinear motion in coordinate systems.The displacement of the tank is expressed as ().Being considered in moving coordinate system (, ) fixed to the tank, the sloshing equations can be gotten as Equation ( 1) is the continuity equation of ideal fluid;  denotes velocity potential.Equations ( 2) and (3) are rigid wall boundary conditions.They indicate that the components of the fluid relative velocity normal to the walls are equal to zero.Equations ( 4) and ( 5) are boundary conditions of the free surface.One is the kinematic condition; the other is the dynamic condition.(, ) is the free surface elevation deviating the undisturbed water level in the moving coordinate system (, ). denotes acceleration of gravity.If (1)-( 5) are solved, the pressure in the fluid can be calculated by the incompressible Bernoulli equation: where  is the pressure and  is the density of fluid.The sloshing force is determined by integrating the liquid pressure over the tank wall area.Because the tank is moving along the -axis and the sloshing models are 2D cases, only  components of sloshing force would have great influence on the dynamics of the tank.In this paper, the sloshing force acting on the left and right walls is considered.

Numerical Process
According to Frandsen [16], -transformation method is used: A new variable  is introduced.The derivatives of the potential function (, , ) with respect to , , and  are transformed into derivatives of Φ(, , ).The first derivatives of  can be expressed as Similarly, it is not difficult to obtain the second derivatives of (, , ).Then ( 1)-( 5) can be transformed by using variable substitution: The fluid domain of ( 1) is 0 <  < , −ℎ <  < (, ).After transformation, it has changed into a fixed rectangular area: 0 <  < , 0 <  < 1.
A finite difference method has been developed to solve ( 9)- (13).Δ is set to be the step size of time.In the iterative process of time, assuming that the values of velocity potential function Φ and the free surface function  are known at time Δ, that is, on  time step, the aim is to obtain the values of the two functions on  + 1 time step.Some difficulties are encountered in the numerical algorithm.First, the governing equation cannot be solved directly on  + 1 time step because free surface , which is unknown at that time, is a boundary curve of the governing equation.Second, there are nonlinear terms existing in free surface boundary conditions.
For the first problem, a prediction method has been introduced.Considering ( 9)- (11) on  + 1 time step,  should have been assigned to the value on  + 1 time step, but it is unknown.On the other hand, (9) on  + 1 time step was transformed from does not appear in the governing equation; it is only a boundary curve.An explicit scheme of (4) can be adopted to make prediction of  +1 and then substitute it into (14).So the four boundary curves of ( 14) are determined, which means that there is only one unknown quantity  +1 in ( 14).This idea should be implemented in ( 9)- (11).Because the computational domain is transformed, an explicit scheme of ( 12) is used to make prediction of  +1 :

𝑘+1
denotes the prediction value.Substitute it into ( 9) and ( 10) on +1 time step; then only Φ +1 is the unknown quantity in the two equations.Equation (11) does not have , so there is no need to do so.These three equations are discretized as For the second problem, an approximate method has been introduced.Consider two boundary conditions ( 12) and ( 13) on  + 1/2 time step; they are complicated nonlinear equations.For example, is a nonlinear term; it should have been assigned to the value on  + 1/2 time step in the calculation process, but here it is assigned to the value on  time step as an approximation.
Other nonlinear terms are dealt with by the same method.Linear terms in equations are discretized by Crank-Nicolson scheme, where  +1 and Φ +1 ( = 1) are set to be unknown quantities.It is noteworthy that  +1 is set to be the unknown quantity in ( 12) and ( 13), so as to update the free surface when the equations are solved.Overall, these two boundary conditions are discretized as Semi-implicit scheme equations of  +1 and Φ +1 ( = 1) have been derived.Moreover, ( 16)-( 21) can be discretized in space and solved by finite difference method.The solution  +1 will be taken as the corrected value of  on  + 1 time step.
After -transformation, (6) has been transformed as When the values of  and Φ are solved, it is not difficult to discretize (22).In the numerical simulation, sloshing force on the left and right walls can be expressed as The resultant sloshing force acting on the tank is

Numerical Results
Some examples of nonlinear sloshing are simulated by this numerical method.They have been investigated by other researchers using various numerical and analytical methods, and many discussions and conclusions have been derived.
A rectangular tank model is established.The natural frequencies of linear transverse sloshing along -axis are (Faltinsen [2]) 4.1.Standing Waves.The tank is fixed, so the moving coordinate system (, ) coincides with the inertial coordinate system (, ).The initial conditions are described as is the amplitude of the initial wave profile, which influences the nonlinearity of the sloshing. indicates the mode number.The convergence of this numerical simulation is investigated.
In the tank, , , and ℎ are set to be 2 m, 1 m, and 1 m, respectively. is set to be 0.2 m;  takes value of 1.The first natural frequency of the sloshing is  1 = 3.76 rad/s.are different.The time and the vibration displacement in the figures are dealt with by dimensionless method: It is found that the results from these meshes and time step sizes are in good agreement, so in this example, a mesh distribution of 40 × 40 and time step size of 0.005 s ensure a reasonable converged solution.
Nonlinear sloshing equations of standing waves can be solved by perturbation techniques.First, nondimensional variables are introduced to the equations.Then,  =  2  / is chosen as the characteristic small parameter.Approximate linear equations are derived by small parameter expansion, and approximate analytical solutions can be obtained.The second-order approximate analytical free surface elevation for the th sloshing mode is prescribed as (Frandsen [16]) where Figure 3(a) shows the time history of the free surface elevation for  = 1 at the left wall by the two methods: numerical solution by the finite difference method proposed in this paper and second-order approximate analytical solution by perturbation method.Initially, good agreement between the two results is achieved, but the phase difference between them increases with time.This phenomenon has been noticed and discussed by many researchers (Chern et al. [18], Frandsen [16]).The conclusion is that the second-order analytical solution is not very accurate for large amplitude sloshing.Figure 3(b) displays numerical wave profiles for nearly half a period from dimensionless time  1  = 31.8to  1  = 34.6.Stationary point does not exist, and nonlinearity is obvious.Figures 3(c) and 3(d) provide the sloshing forces on the left wall and right wall, respectively (the forces are also nondimensionalized).It shows that slight double peaks appear in the time history.This has also been observed by Wu et al. [5].
Figure 4(a) shows the comparison of numerical solutions and approximate analytical solutions.Figure 4(b) shows numerical wave profiles for nearly half a period.Phase difference and nonlinearity are also observed.

Sloshing Caused by Horizontal Excitations.
The tank is assumed to take the horizontal oscillation; the external excitation acting on the tank is defined as harmonic motion: ℎ and  ℎ are the amplitude and angular frequency of the excited motion.The fluid is assumed to be stationary at  = 0.Then, in the moving coordinate system, the initial conditions are described as    In order to verify the precision of this numerical method for the sloshing models in horizontally moving tanks, two examples are calculated and compared with the solutions of Frandsen [16]: (1)  ℎ = 0.7 1 ,  ℎ = 0.036; (2)  ℎ = 1.3 1 ,  ℎ = 0.072.The grid size of the numerical method is 40 × 40 and the time step size is 0.006 s. Figure 5(a) shows that the two simulations are in close agreement.In the second example, the nonlinear effect is strong, so there is a little difference between the two solutions in Figure 5 numerical methods, such as the higher peaks and the less deep troughs.

Resonance Phenomenon.
Two situations of resonance are considered: the first resonance mode and the third resonance mode.The external excitations of the tank are described as  1 () = 0.007 cos( 1 ),  2 () = 0.007 cos( 3 ).
The amplitudes are the same, but the excited frequencies are different.The time histories of free surface elevation at the left wall and the resultant sloshing forces are calculated.Figure 6(a) shows that the sloshing of the first resonance mode is violate; the amplitude of the vibration is increasing.Figure 6(b) shows that the amplitude of the third resonance mode is much smaller than that of the first resonance mode.
In Figure 7, it can be concluded that the resultant sloshing force of the third resonance mode is trivial compared with the force of the first resonance mode.In fact, the first natural frequency of sloshing is the most influential of all the natural frequencies.

The Limitation of the Numerical Method.
The present numerical method is employed to simulate nonlinear sloshing in moving tanks.If the nonlinearity is quite strong, this method may lose its effectiveness.Two excitation frequencies are considered in the simulations.One is  ℎ = 0.85 1 ; the other is  ℎ = 0.95 1 , which indicates near resonant sloshing situation.Different excitation amplitudes  ℎ and aspect ratios ℎ/ are tested, and the limitation of the convergence of this numerical method is found out.First, the aspect ratio is set to be ℎ/ = 1/2, just as before.Excitation amplitude is increased in the numerical simulation, and the sloshing becomes more and more violent.For  ℎ = 0.85 1 , if the excitation amplitude is  ℎ = 0.06 m, Mathematical Problems in Engineering the simulation is divergent after several times of iteration; if  ℎ = 0.05 m, the simulation lasts for a long time.So it is considered that the limitation of  ℎ for nonresonant sloshing is 0.05 m.For  ℎ = 0.95 1 , the excitation amplitude is much smaller.The maximum limit of  ℎ is estimated to be 0.012 m.The time histories of free surface elevation at the left wall in the limit cases are shown in Figure 8.
Then the aspect ratio is changed by decreasing the still water depth ℎ.For the aspect ratio ℎ/ = 1/8, the numerical method still works, with a mesh size of 80 × 20.If the aspect ratio is less than 1/8, the simulation is divergent.As a result, the minimum limit of aspect ratio ℎ/ is considered to be 1/8.In this critical aspect ratio situation, for  ℎ = 0.95 1 , the maximum limit of excitation amplitude  ℎ is still 0.012 m, while for  ℎ = 0.85 1 , the maximum limit of  ℎ is 0.03 m. Figure 9 shows the time histories of free surface elevation in the limit cases in this critical aspect ratio.The phenomena of beating and resonance are very obscure; it can be concluded that the sloshing in low aspect ratio is quite different from the one in normal aspect ratio.

Beating Phenomenon in Different Aspect
Ratios.Sloshing in a fixed aspect ratio has been studied in many articles.However, little work has been done to analyze the response of liquid sloshing in different aspect ratios by the same external excitation.When the external excitation is close to the first natural frequency of sloshing, the beating phenomenon can be observed.The study here is intended to investigate the effect of different aspect ratios on beating phenomenon and sloshing force.
In the numerical examples calculated before, the length and width of the tank are set to be 2 m and 1 m; the still water depth is 1 m.The first natural frequency of sloshing is  1 = 3.76 rad/s.The aspect ratio is ℎ/ = 1/2.Now the length  and the width  of the tank are fixed, and the aspect ratio is changed by changing the still water depth ℎ.The former water depth and natural frequency are set to be the characteristic quantity: ℎ * = 1 m,  * 1 = 3.76 rad/s.The relationship between aspect ratio and the first natural frequency must be considered.From (25), it is easy to get the graph of the natural frequency as a function of aspect ratio.Figure 10 indicates that the slope of the function curve is steep when the aspect ratio is small; if the aspect ratio is greater than 0.6, the function curve is nearly horizontal.
Under the same external excitation, the sloshing forces of beating are affected by three factors: proximity of excitation frequency to the natural frequency, liquid mass, and nonlinearity of sloshing.When the aspect ratio of the liquid changes, these three factors would change along with it.In order to discuss the effect of aspect ratio, two horizontal oscillations of the tank are studied.First, the tank is subjected Mathematical Problems in Engineering  to a low-frequency excitation.The movement of the tank is described as  ℎ = 0.03 m,  ℎ = 3.4 rad/s.Figure 10 suggests that if the aspect ratio is 0.3, the corresponding natural frequency is 3.4 rad/s.Our aim is to analyze the influence of this low-frequency excitation on the high aspect ratio sloshing system.Figure 11 shows the time histories of free surface elevation at the left wall in four different aspect ratios, where the parameters are nondimensionalized: The amplitude of free surface elevation is quite large in Figure 11(a), which shows strong nonlinearity of sloshing, while in Figures 11(b), 11(c), and 11(d), the amplitudes are nearly the same, smaller than the amplitude in Figure 11(a), so the nonlinear effects are weaker in these three sloshing situations.
Figure 12 shows the time histories of resultant sloshing force acting on the tank (the forces are also nondimensionalized:   = /ℎ * ).The amplitudes of the sloshing force are extracted in Figure 12.In order to study the effect of nonlinear factor, linear sloshing models by Abramson [1] are calculated under the same condition to make comparisons.The relationship between the amplitude of sloshing force and aspect ratio is shown in Figure 13.It is obvious that both linear and nonlinear results have the same trend: decreasing first and then increasing.On the one hand, the rise of the aspect ratio increases the natural frequency, which makes the sloshing far from resonance and reduces sloshing forces.On the other hand, the rise of the aspect ratio increases the mass of the liquid, thereby increasing the sloshing forces.For the sloshing in aspect ratio ℎ/ = 0.5, the excitation frequency is close to the natural frequency.If aspect ratio is increased, the factor of reducing sloshing force is more dominant.However, for the sloshing in high aspect ratios ℎ/ > 0.75, the rise of aspect ratio brings about slight changes in the natural frequency, so the factor of increasing sloshing force is more dominant.Comparing nonlinear results with linear ones in Figure 13, some conclusions are obtained.First, the nonlinearity of sloshing always decreases sloshing force.Second, to a certain extent, the difference between linear and nonlinear amplitude reflects the effect of nonlinear factor on sloshing force.For the sloshing in aspect ratio ℎ/ = 0.5, the difference is the largest, while for the other sloshing situations, the differences are nearly the same.This result accords with the conclusion in Figure 11.Third, nonlinearity is an important factor in large amplitude sloshing, which should not be neglected.Then, the tank is subjected to a high-frequency excitation.The movement of the tank is described as  ℎ = 0.015 m,  ℎ = 4 rad/s.Our aim is to analyze the influence of this high-frequency excitation on the low aspect ratio sloshing system.Figure 14 shows the time histories of free surface elevation at the left wall in four different aspect ratios.As the aspect ratio decreases, the amplitude and period of free surface elevation are diminishing.For the sloshing in low aspect ratios ℎ/ < 0.5, the reduction of aspect ratio brings about significant changes in the natural frequency, so the beating phenomenon has changed greatly even if the aspect ratio is slightly decreased.The resultant sloshing forces acting on the tank are also calculated.The relationship between the amplitude of sloshing force and aspect ratio is showed in Figure 15.The reduction of the aspect ratio makes the natural frequency far from the excitation frequency; meanwhile, it decreases the mass of the liquid in the tank.So both of the factors reduce sloshing force.As the aspect ratio decreases, the difference between linear and nonlinear amplitude of sloshing force also decreases.The reason is that nonlinear effects are becoming weaker, and the sloshing force itself is also decreasing.

Conclusions
A new finite difference method based on the prediction of free surface has been developed and used to simulate nonlinear sloshing phenomenon in tanks by the potential flow theory.After -transformation, the fluid domain is transformed to a rectangle, and a predictor-corrector scheme for free surface is used in the iteration of time; then the governing equation and boundary conditions are approximately linearized and solved in every time step.The convergence of this method is surveyed.The numerical examples of standing waves and horizontally forced sloshing models have been calculated, and good agreement between the numerical simulation and results in other literatures has been obtained.This numerical method is simple and easy in programming, with excellent accuracy and stability, and can be extended to other sloshing problems.The time histories of free surface elevation and sloshing forces in the first and third resonance situation are also presented and analyzed.The beating phenomenon of sloshing in different aspect ratios under the same excitation has been numerically simulated by this method.The results show that the effect of the aspect ratio on the sloshing force is complicated.For the lowfrequency excitation acting on the high aspect ratio models, the rise of aspect ratio makes the sloshing far from resonance and increases the liquid mass, so the sloshing force decreases first and then increases.For the high-frequency excitation acting on the low aspect ratio models, the reduction of aspect ratio makes the sloshing far from resonance and decreases the liquid mass, so the sloshing force decreases all the way.Since the first natural frequency of liquid sloshing in not a linear function of aspect ratio, the change of aspect ratio produces much stronger effect in low aspect ratio sloshing systems.Moreover, nonlinearity of sloshing always decreases sloshing force.The effect of nonlinear factor depends on the violence of sloshing.

Figure 2 (Figure 2 :Figure 3 :Figure 4 :
Figure 2: Time history of free surface elevation for different meshes and time step sizes.

Figure 5 :Figure 6 :
Figure 5: Comparison of free surface elevation at the left wall.

Figure 11 :
Figure 11: Free surface elevation at the left wall (low-frequency excitation).

Figure 14 :
Figure 14: Free surface elevation at the left wall (high-frequency excitation).