Numerical simulation of liquid sloshing using arbitrary Lagrangian-Eulerian level set method

Liquid sloshing in an oscillating tank is simulated numerically as a fluidstructure interaction problem using the arbitrary Lagrangian-Eulerian and the level set coupled method, in which the computational grid points are moved with the velocity of the tank. It is shown by comparing the simulation results with the existing experimental results that the sloshing behavior of the free surface is predicted well by the present method. The simulation results are also compared with the case using the body force, in which the body force term is included in the fluid equations and the grid points are not moved. The difference between the moving grid method and the body force method is made clear both theoretically and numerically, and the limitation of the body force method is discussed.


INTRODUCTION
Thermal-hydraulic phenomena with two-phase flows are seen widely in nuclear engineering fields, and predictions of complicated interfacial phenomena are of practical importance.Characteristics of two-phase flows have been intensively studied both experimentally and numerically under wide variety of thermal-hydraulic conditions concerning with nuclear reactor safety.Two-phase flow phenomena under seismic conditions are, however, not well known.In some reactor components, induced fluid motion may, in turn, result in large pressure impact on structures.Thermal conditions such as heat transfer between fluid and structure are also affected by variation of two-phase flow conditions.Free surface behaviors of liquid sodium have been studied for fast breeder reactors (FBRs) [1][2][3][4][5].Numerical simulations were performed in some studies to obtain the surface motion, where the motion of reactor tank was taken into account as the external acceleration term in fluid equations [1,2,4].Stability analyses of boiling water reactors (BWRs) under seismic conditions have been performed by modifying the safety analysis code TRAC-BF1 to take into account the effect of seismic oscillation on thermal hydraulics [6].The oscillating acceleration was added to the momentum equation as an external body force term, as was the case for FBRs, and the coupled effect of the thermal hydraulics and the point kinetics was discussed.In these studies, seismic effects on thermal-hydraulics were modeled through the additional body force term in the fluid equations, instead of taking into account the oscillation of reactor components.Although the method using the external body force could easily be applied for reactor safety analyses, validation of this methodology has not yet been discussed.
In this study, the stratified two-phase flow field and the motion of the free surface in an oscillating tank are simulated numerically as a sample problem of the thermal-hydraulic behavior under seismic conditions.The growth of the surface wave in a partially filled oscillating tank is known as sloshing, and is important for structural integrity of the tank.Sloshing has been studied experimentally and numerically, for instance in relation to sea transport of oil and liquefied natural gas [7,8], seismic response of liquefied petroleum gas tank in petrochemical industry [9], and so on.The effect of oscillating tank on fluid motion was, in some cases, taken into account by including a body force induced by the tank motion in the momentum equation of fluid [7,8], as was the case for nuclear reactors [1,2,4,6].The computational grid for fluid simulation was, in other cases, moved directory according to the tank motion [9].Although the method using the body force is easy from the view point of numerical simulation, the method using the moving grid seems apparently to be corresponding to the real phenomena.
A simple sloshing experiment [10] is simulated in the following first using the moving grid.A stratified two-phase flow is contained in a rectangular tank, and the tank is set in an oscillatory motion.Incompressible Navier-Stokes equations are solved using the level set method [11].In the level set method, the level set function, which is the distance function from the free surface, is calculated by solving the transport equation using the flow velocities.The motion of the tank is modeled by the Arbitrary Lagrangian-Eulerian (ALE) method [12], where the computational grid points are moved with the velocity of the tank.Both the liquid-phase and the gas-phase flow fields with the free surface motion induced by the oscillating tank are thus obtained in this study.It is shown by comparing the simulation results with the experimental results that the sloshing behavior of the free surface is predicted well using the moving grid.The simulation results are compared next with the case using the body force, and the numerical approach using the body force is discussed.The effects of amplitude and frequency on the simulation results are also shown.

NUMERICAL SIMULATION 2.1. GOVERNING EQUATIONS
Governing equations and numerical schemes are similar to those used in [13], but the compressibility is not included and the gravitational acceleration is taken into account.The governing equations for the two-phase flow field are the equation of continuity and the incompressible Navier-Stokes equations: (1) and (2) where ρ, u, p and µ, respectively, are the density, the velocity, the pressure and the viscosity, D is the viscous stress tensor, F s is a body force due to the surface tension, and g is the gravitational acceleration.The surface tension force is given by (3) where σ, κ, δ and φ are the surface tension, the curvature of the interface, the Dirac delta function and the level set function, respectively.The level set function is a distance function defined as φ = 0 at the free surface, φ < 0 in the liquid region, and φ > 0 in the gas region.The curvature is expressed in terms of φ: The density and viscosity are given, respectively, by where the subscripts g and l denote gas and liquid phases, respectively, and H is the smeared Heaviside function defined by (7) where ε is a small positive constant for which ∇φ = 1 for | φ | ≤ ε.The time evolution of φ is given by (8) In this study, the ALE method is applied, and the computational grid is moving with the same velocity as the velocity of the oscillating tank.The substantial derivative terms in Eqs. ( 2) and ( 8) are thus defined by (9) where U is the velocity of the computational grid or the oscillating tank.
In order to maintain the level set function as a distance function, an additional equation is solved: (10) where τ and α are an artificial time and a small constant, respectively.The level set function becomes a distance function in the steady-state solution of the above equation.The following equation is also solved to preserve the total mass of liquid and gas phases in time [14]; where M denotes the mass corresponding to the level set function and M 0 denotes the mass for the initial condition.
The finite difference method is used to solve the governing equations.The staggered mesh is used for spatial discretization of velocities.The convection terms are discretized using the second order upwind scheme and other terms by the central difference scheme.Time integration is performed by the second order Adams-Bashforth method.The SMAC method [15] is used to obtain pressure and velocities.

SIMULATION CONDITIONS
The simulation conditions are almost the same as the conditions of the sloshing experiment [10].The size of the tank is 1.0 m × 1.2 m × 0.1 m, and the initial water level is 0.5 m as shown in Fig. 1.The tank is set in an oscillatory motion in one horizontal direction.The oscillation of the tank location in the horizontal direction is given by (12) where A = 0.0093 m and ω = 5.311 rad/s are, respectively, the amplitude and the angular frequency of the oscillation.The velocity of the oscillating tank is used in the present moving grid method and is given as the differential of the tank location, (13) In this study, the case with the oscillating body force is compared with the moving grid method.The oscillating body force is given as the differential of the tank velocity, (14)  The above body force is applied as the additional external force term in the momentum equation given by Eq. ( 2), and the tank is not moved and U = 0 in Eq. ( 9) for the body force method.The slip boundary conditions are applied at all walls for both the moving grid and the body force methods.

RESULTS AND DISCUSSION
3.1.GRID DEPENDENCY AND THREE-DIMENTIONAL EFFECT Two-dimensional calculations are performed first to evaluate the grid dependency.Time evolutions of the surface elevation at the left and right walls of the tank are shown in Fig. 2, where the number of calculation grid points is varied from 50 × 60 to 200 × 240.Equal spaced square grid is used, and thus the grid size is 0.02 m to 0.005 m.The Courant number is set at a small constant of about 0.35, and the time step size is from 0.002 s to 0.0005 s according to the grid size.The oscillation of surface elevation increases with time and the grid dependency is shown to be small for the grid size smaller than 0.01 m even for the largest elevation in Fig. 2. The grid size is thus set equal to 0.01 m, and the number of grid points is 100 × 120 for the two-dimensional case.Three-dimensional calculation is performed next to see the effect of the horizontal depth of the tank.The same grid size is used and the number of grid points in the direction of the horizontal depth is 10 for the three-dimensional calculation.Calculated surface elevations at the mid depth are compared with the two-dimensional results in Fig. 3.It is confirmed that the three-dimensional effect is negligibly small in this sloshing simulation, and the twodimensional calculations with 100 × 120 grid points are performed in the following.

COMPARISON WITH EXPERIMENTAL RESULTS
The shapes of the free surface were recorded in the experiment [10].Simulated surface shapes at 3.54 s and 7.08 s are compared with the experimental observations in Fig. 4.These two timings are corresponding to the maximum surface deformation as shown in Figs. 2 and  3.At these two timings, the positive elevation at the right wall is larger than the negative elevation at the left wall, and the neutral elevation location is slightly shifted from the horizontal center to the right side.It is shown that the agreement between the simulation and the experiment is satisfactory, especially at 3.54 s.The sloshing experiment is thus found to be simulated well by the present numerical approach using the ALE and the level set couplet method.

COMPARISON WITH BODY FORCE METHOD
The time evolution of the surface elevation obtained by the present moving grid method is compared with that by the body force method in Fig. 5.The growth of the surface wave is shown to be almost the same for the two methods.The difference of the two methods is thus discussed in the following.The momentum equation for the moving grid method is given by Eqs. ( 2) and ( 9) as   The flow velocity is then assumed to be divided into two parts: the grid velocity U and the induced velocity u′, (16) The momentum equation then becomes (17) where D′ is the viscous stress tensor for the induced velocity u′, and (18) is assumed.Equation ( 18) indicates that the grid velocity is not varied spatially, since the induced velocity is not zero generally.The last term in the right hand side of Eq. ( 17) is the oscillating body force given by Eq. ( 14), and Eq. ( 17) is the momentum equation for the body force method.It is thus obvious that Eq. ( 15) for the moving grid method is equivalent to Eq. ( 17) for the body force method.The calculated surface elevation by the moving grid method thus becomes the same as the result by the body force method as shown in Fig. 5.
It should, however, be noted that the calculated velocity field by the moving grid method includes the grid velocity, while that by the body force method does not.In order to see this difference clearly, the flow fields are compared in Fig. 6, where the shapes of the free surface and the velocity fields at 3.54 s and 7.08 s are shown.The velocity field obtained by the moving grid method is shown on the left, the induced velocity field obtained by the body force method is in the center, and the sum of the induced velocity field and the grid velocity is on the right.It is noted that the vector scale at 3.54 s is 4 times larger than that at 7.08 s.The induced flow field calculated by Eq. ( 17) is much different from the velocity field by Eq. ( 15) at 3.54 s, though the difference becomes small as the velocity becomes large at 7.08 s.It is, however, shown in Fig. 6 that the velocity field by the moving grid method corresponds to the sum of the induced velocity field by the body force method and the grid velocity.In other words, the velocity field by the moving grid method is observed on the fixed coordinate, while that by the body force method is on the moving coordinate with the grid velocity.This difference of coordinate system should be reminded for estimation of the calculated results obtained by the body force method.

EFFECTS OF OSCILLATION AMPLITUDE AND FREQENCY
In order to see the effects of oscillation amplitude and frequency on the surface elevation, the amplitude A and the frequency ω of tank oscillation are, respectively, increased in Figs.7 and 8.The results with two times and four times larger amplitudes are shown in Fig. 7 (a) and (b), respectively.The grid velocity and the body force are simply increased in proportion to the oscillation amplitude of the tank according to Eqs. ( 13) and (14).It is shown in Fig. 7 that the surface elevation becomes larger with increase in the oscillation amplitude.The elevation is, however, not increased linearly with the amplitude, and the frequency of surface elevation becomes slightly smaller, though the oscillation frequency of the tank is the same.Besides the maximum and the minimum elevations at the wall, it is of interest as shown in Fig. 7(b) that the neutral elevation, for which the right and the left elevations become equal, increases from zero as the maximum elevation becomes large.This indicates that the neutral surface shape is not flat as the surface deformation becomes large during sloshing.These nonlinear effects are found to be simulated by both the moving grid and the body force methods, and the calculated results are almost the same as shown in Fig. 7.
The results with two times and four times larger frequency are shown in Fig. 8 (a) and (b), respectively.It is noted that not only the frequency but also the effective amplitude of the grid velocity become two or four times larger according to Eq. ( 13) as the oscillation frequency becomes two or four times larger.Furthermore, the effective amplitude of the body force becomes four or sixteen times larger according to Eq. ( 14).In any case, the sloshing phenomena with the growth of surface wave are not seen in Fig. 8, since the increased frequency is different from the resonant frequency used in the experiment.Large oscillation does not appear, but complicated surface fluctuations occur.It is seen that the surface elevation becomes larger as the oscillation frequency increases, since the effective amplitudes of grid velocity and the body force are simultaneously increased.These effects of oscillation frequency are found to be simulated by both the moving grid and the body force methods, and the calculated results are again almost the same as shown in Fig. 8.

CONCLUSION
The stratified two-phase flow field and the motion of the free surface in the oscillating tank have been simulated numerically to assess the methodology for simulating reactor thermal hydraulics under seismic conditions.The moving grid method, where the oscillating tank was modeled directly using the moving grid of the ALE method with the oscillating velocity,  was compared with the body force method, where the effect of tank motion was simulated using an oscillating body force acting on the fluid in the stationary tank.The stratified twophase flow field in the tank was simulated using the level set method in both cases.
It was found from the comparison with the existing experimental reuslts that the sloshing behavior of the free surface was predicted well by the moving grid method.The momentum equation of the body force method was shown to be equivalent to that of the moving grid method, and the calculated results by the body force method coinsided with those by the moving grid method.The calculated velocity field by the body force method was, however, the induced velocity field, which did not include the oscillating velocity of the tank.It was found that the sum of the induced velocity field and the oscillating velocity was corresponding to the velocity field calculated by the moving grid method.It should be noted that the oscillating velocity was assumed to be spatially constant for the body force method.If whole the simulation region does not move simultaneously, for instance in case of deformation of components under seismic conditions, the body force method is not equivalent to the moving grid method.The governing equations would be complicated for simulating cases with spatially different body forces.
The effects of oscillation amplitude and frequency were also discussed.The oscillation of surface elevation was increased with increase in the oscillation amplitude of the tank.The frequency of the surface oscillation was found to be affected by the amplitude.It was shown that the effect of oscillation amplitude of the tank could be discussed simply by changing the amplitude of the grid velocity or of the body force.The large sloshing phenomena were not seen and the oscillation of surface elevation became complicated as the oscillation frequency increased.It was shown that the increase in the oscillation frequency of the tank corresponded not only to the increase in the frequencies of the grid velocity and the body force but also to the increase in effective amplitudes of the grid velocity and the body force.This point should be reminded for discussing the effect of frequency on the results simulated both by the moving grid and by the body force methods.
In this study, the simple sloshing experiment was simulated numerically as a sample problem, and the governing equations had no empirical correlations.For simulating reactor thermal hydraulics, two-fluid model codes such as TRAC would be used generally.Such safety analysis codes include large number of empirical correlations, which are obtained under the static conditions with constant gravitational acceleration.The empirical correlations and related terms should thus be treated carefully for simulating thermal hydraulics under oscillating conditions not only by the body force method but also by the moving grid method.

Figure 1
Figure 1 Schematic of sloshing experiment.

Figure 2
Figure 2 Effect of grid size on time evolution of surface elevation.

344Figure 3
Figure 3 3D effect on time evolution of surface elevation.

Figure 4
Figure 4  Comparison of free surface shape between simulation and experiment.

Figure 5
Figure 5  Comparison of surface elevation between moving grid and body force methods.

Figure 6
Figure 6 Comparison of flow field between moving grid and body force methods.

348Figure 7
Figure 7 Effects of oscillation amplitude on the time evolution of surface elevation.

Figure 8
Figure 8 Effects of oscillation frequency on the time evolution of surface elevation.