David Le Touzé

Wenjin Cao

Zhe Li
email: zhe.li@ec-nantes.fr

Xuhui Li

A regularized single-phase lattice Boltzmann method for free-surface flows

Keywords: lattice Boltzmann method, free-surface flows, regularization, pressure noises

A regularized single-phase lattice Boltzmann method for free-surface flows

Introduction

Free-surface flows can be widely observed in nature and industrial applications, such as the ocean waves, the sloshing problem in the water tank and the metal foaming process, etc. This special type of two-phase flows are often of large deformation and high density ratio, and sometimes with violent impacting phenomena.

Efficiently simulating the free-surface flows still remains a challenging task for the conventional computational fluid dynamics (CFD) methods. Over the past several decades, the lattice Boltzmann method (LBM) emerged as an effective CFD tool with some attractive features, such as the linear advection, local collision process, and being free from solving the Poisson equation, which makes it well-suite for the large scale parallel computation.

In this context, the LBM can be a powerful tool for studying the hydrodynamic problems involving free-surface flows resorting to appropriate free-surface boundary condition to model the atmospheric pressure.

Two branches of the existing free-surface LB models can be found in the literature, which are the multiphase LB models and the single-phase LB models. Most of the multi-phase LB models fall into four categories [START_REF] Li | Lattice Boltzmann methods for multiphase flow and phase-change heat transfer[END_REF][START_REF] Liu | Multiphase lattice Boltzmann simulations for porous media applications[END_REF], which are [START_REF] Bhatnagar | A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems[END_REF] the color-gradient model [START_REF] Gunstensen | Lattice Boltzmann model of immiscible fluids[END_REF] where two colored particle distribution functions are used for two immiscible fluids and a recoloring step is designed for the phase separation; [START_REF] Botella | Benchmark spectral results on the lid-driven cavity flow[END_REF] the free-energy model [START_REF] Swift | Lattice Boltzmann simulation of nonideal fluids[END_REF] where a non-ideal thermodynamic pressure tensor is introduced into the 2 nd -order moment of the distribution function, enabling the phase separation to be governed by the non-ideal equation of state; [START_REF] Chapman | The mathematical theory of non-uniform gases: an account of the kinetic theory of viscosity[END_REF] the pseudopotential model [START_REF] Shan | Lattice Boltzmann model for simulating flows with multiple phases and components[END_REF] which defines an inter-particle potential to mimic the repulsive force between particles from different phases; (4) the phase-field model [START_REF] He | A lattice Boltzmann scheme for incompressible multiphase flow and its application in simulation of Rayleigh-Taylor instability[END_REF] where an additional index distribution function is employed to recover the Cahn-Hillard-like equation that describes the phase separation.

As an alternative, the single-phase LB models reduce the multi-phase flow into a single-phase approach with a free-surface boundary, which are also effective in the assumption that the dense-phase flow is usually considered to be negligibly affected by the light-phase motion. In this spirit, several authors [START_REF] Ginzburg | Lattice Boltzmann model for free-surface flow and its application to filling process in casting[END_REF][START_REF] Körner | Lattice Boltzmann model for free surface flow for modeling foaming[END_REF] proposed the concept of liquid volume fraction from the volume-of-fluid (VOF) method to describe the free-surface. The free-surface movement is driven by the cell mass exchange, which is calculated by a built-in operation during the streaming step in LB methods. The distribution functions on the interface are reconstructed by a dynamic boundary condition that ensures the balanced hydrodynamic force against the gas pressure. Similarly, Janssen and Krafczyk [START_REF] Janssen | A lattice Boltzmann approach for free-surface-flow simulations on non-uniform block-structured grids[END_REF] introduced another VOF-based model, where the free-surface boundary condition is the same as [START_REF] Körner | Lattice Boltzmann model for free surface flow for modeling foaming[END_REF], while they adopted a piecewise linear reconstruction method (PLIC) [START_REF] Gueyffier | Volume-of-fluid interface tracking with smoothed surface stress methods for three-dimensional flows[END_REF][START_REF] Youngs | Time-dependent multi-material flow with large fluid distortion[END_REF] for better capturing the free-surface. In addition, Thürey and Rüde [START_REF] Rüde | Free surface lattice-Boltzmann fluid simulations with and without level sets[END_REF] proposed a level-set free-surface model based on the free-surface model in [START_REF] Körner | Lattice Boltzmann model for free surface flow for modeling foaming[END_REF], where a level-set front tracking method [START_REF] Sethian | Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science[END_REF] was used to determine the flag change.

In free-surface flows the dynamics of the light phase is neglected due to the high density ratio, such as O (1000) between water and air. In [START_REF] Grunau | A lattice Boltzmann model for multiphase fluid flows[END_REF][START_REF] Li | Lattice Boltzmann methods for multiphase flow and phase-change heat transfer[END_REF][START_REF] Liu | Multiphase lattice Boltzmann simulations for porous media applications[END_REF] it is summarized that the applicable density ratio for colorgradient model can be up to O (1000) for stationary tests [START_REF] Leclaire | Numerical evaluation of two recoloring operators for an immiscible two-phase flow lattice Boltzmann model[END_REF], whereas it is restricted to O (10) for dynamic problems due to numerical instability issues. Another comment is seen in [START_REF] Li | Lattice Boltzmann methods for multiphase flow and phase-change heat transfer[END_REF] that the cause of the instability of the free-energy model were assessed as the operation of adding density-gradient associated correction terms in order to remove the non-hydrodynamic terms in the macroscopic equations. As for the pseudo-potential models, there are some published works on applications with a density ratio up to O (1000), such as [START_REF] Huang | Forcing term in single-phase and Shan-Chen-type multiphase lattice Boltzmann models[END_REF][START_REF] Kupershtokh | On equations of state in a lattice Boltzmann method[END_REF][START_REF] Li | Forcing scheme in pseudopotential lattice Boltzmann model for multiphase flows[END_REF][START_REF] Sun | Evaluation of force implementation in pseudopotential-based multiphase lattice Boltzmann models[END_REF].

In the phase-field category, Inamuro et al. [START_REF] Inamuro | An improved lattice Boltzmann method for incompressible two-phase flows with large density differences[END_REF] and Lee et al. [START_REF] Lee | A stable discretization of the lattice Boltzmann equation for simulation of incompressible two-phase flows at high density ratio[END_REF] proposed their models for large density ratio flows. However, the model of Inamuro et al. [START_REF] Inamuro | An improved lattice Boltzmann method for incompressible two-phase flows with large density differences[END_REF] involves pressure correction processes which might be time-consuming, and the model of Lee et al. [START_REF] Lee | A stable discretization of the lattice Boltzmann equation for simulation of incompressible two-phase flows at high density ratio[END_REF] is rather not suitable for large velocity flows [START_REF] Fakhari | Phase-field modeling by the method of lattice Boltzmann equations[END_REF][START_REF] Li | Lattice Boltzmann methods for multiphase flow and phase-change heat transfer[END_REF][START_REF] Zheng | A lattice Boltzmann model for multiphase flows with large density ratio[END_REF].

To summarize, it seems that the multi-phase models are less suitable for large density ratio free-surface flows.

Besides the single-phase models possess another attractive feature that the absence of calculation in the lightphase domain saves considerable computational time. Hence, in the present work, the single-phase LB model [START_REF] Körner | Lattice Boltzmann model for free surface flow for modeling foaming[END_REF] is adopted for simulating free-surface flows, such as the water-air flow in the presence of impacting problems.

It is well known that although the standard single-relaxation-time (SRT) LBM [START_REF] Bhatnagar | A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems[END_REF] has been widely applied because of its simplicity and efficiency, it may suffer from severe numerical instabilities at high Reynolds numbers (the relaxation time is close to 0.5). To overcome this imperfection, several improved LBM schemes can be found in the literature. The two-relaxation-time (TRT) LBM [START_REF] Ginzburg | Two-relaxation-time lattice Boltzmann scheme: About parametrization, velocity, pressure and mixed boundary conditions[END_REF] decomposes the population solution into its symmetric and anti-symmetric components, with a dependent relaxation time for each. The multi-relaxationtime (MRT) LBM [START_REF] Lallemand | Theory of the lattice Boltzmann method: Dispersion, dissipation, isotropy, Galilean invariance, and stability[END_REF] offers more freedom by projecting the distribution functions onto a moment space and relaxing each moment in an individual rate. The entropic LBM [START_REF] Karlin | Perfect entropy functions of the lattice Boltzmann method[END_REF] improves the numerical stability by restoring the second law of thermodynamics [START_REF] Chikatamarla | Entropic lattice Boltzmann models for hydrodynamics in three dimensions[END_REF]. Alternatively, the regularized LBM (RLBM) has also attracted significant attention. In 2006, Zhang et al. [START_REF] Zhang | Efficient kinetic method for fluid simulation beyond the Navier-Stokes equation[END_REF] provided a straightforward explanation based on the Hermite representation of the LB equation. The regularization is a pre-collision procedure prior to the standard SRT-LBM, thus it is expected to have better numerical stability while maintaining the simplicity of the SRT-LBM. Latt et al. [START_REF] Latt | Straight velocity boundaries in the lattice Boltzmann method[END_REF] introduced another regularization method based on the Chapman-Enskog analysis, where the non-equilibrium distribution function can be computed from the velocity gradient. Recently, a comprehensive comparison of collision models in the LB framework has been given by Corexias et al. [START_REF] Coreixas | Comprehensive comparison of collision models in the lattice Boltzmann framework: Theoretical investigations[END_REF]. Among these schemes, the RLBM is a successful and simple variation of LBM, however its application on free-surface problems has not been seen in the literature.

In the present work, we first introduce the regularized collision model [START_REF] Zhang | Efficient kinetic method for fluid simulation beyond the Navier-Stokes equation[END_REF] to the single-phase free-surface LB model [START_REF] Körner | Lattice Boltzmann model for free surface flow for modeling foaming[END_REF], and we demonstrate that it is effective to stabilize the violent free-surface flows and to filter out the spurious pressure noise. As indicated in [START_REF] Coreixas | Recursive regularization step for high-order lattice Boltzmann methods[END_REF], the adopted RLBM applies a projection-based regularization, which can help improve the stability in the numerical simulation of the isothermal weakly-compressible fluid flows. Moreover, inspired by the work of [START_REF] Latt | Straight velocity boundaries in the lattice Boltzmann method[END_REF], we propose a new reconstruction method for the distribution functions on the free-surface, which is expected to be more consistent with the regularized collision model.

The rest of the paper is organized as follows: Section 2 gives the basic equations of the LB method and presents the free-surface model. The theoretical explanation of the RLBM is illustrated in Section 3. A new reconstruction method is proposed in Section 4. Several numerical validations are discussed in Section 5. Finally, the conclusions are drawn in Section 6.

2. Lattice Boltzmann method with the single-phase free-surface model

Lattice Boltzmann equation with force term

Introducing the force term from [START_REF] Guo | Discrete lattice effects on the forcing term in the lattice Boltzmann method[END_REF], the standard lattice Boltzmann equation with the single-relaxation-time collision model [START_REF] Bhatnagar | A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems[END_REF] reads

f i (x + ξ i ∆t, t + ∆t) -f i (x, t) = - ∆t τ [f i (x, t) -f eq i (x, t)] + ∆tF i (x, t) , (1)
where f i (x, t) is the distribution function at position x and time t, in the i th of the q directions in a given d-dimensional DdQq lattice. Accordingly, ξ i represents the lattice speed in the i th direction. In the present work, the D2Q9 lattice is adopted, and its lattice speeds are given as

ξ i =              (0, 0) , i = 0 cos (i -1) π 2 , sin (i -1) π 2 , i = 1, ..., 4 √ 2 cos (2i -9) π 4 , sin (2i -9) π 4 , i = 5, ..., 8 (2)
On the right-hand side of Equation ( 1), the first term denotes a relaxation procedure from distribution function f i to its equilibrium f eq i with τ being the relaxation time. The equilibrium distribution function is determined by the macroscopic variables as

f eq i = w i ρ 1 + ξ i • u c 2 s + uu : ξ i ξ i -c 2 s I 2c 4 s , (3)
where ρ is the fluid density, u is the macroscopic fluid velocity, and I is the identity tensor. The lattice sound speed c s and the weights w i are lattice dependent parameters, which are given as c s = 1 √ 3 ∆x ∆t and

w i =              4 9 , i = 0 1 9 , i = 1, ..., 4 1 36 , i = 5, ..., 8 (4)
for the adopted D2Q9 lattice.

The last term on the right-hand side of Equation ( 1) is the force term, which is given as [START_REF] Guo | Discrete lattice effects on the forcing term in the lattice Boltzmann method[END_REF]]

F i = 1 - ∆t 2τ w i ξ i -u c 2 s + (ξ i • u) ξ i c 4 s • F , ( 5
)
where ∆t is the time step and F is the external body force. One should note that when adopting the force model given in [START_REF] Guo | Discrete lattice effects on the forcing term in the lattice Boltzmann method[END_REF], the relation between the distribution function and the macroscopic variables are written as

         ρ = i f i , ρu = i ξ i f i + ∆t 2 F , (6)
2.2. Single-phase free-surface model

Description of the free-surface

In the present work, the single-phase free-surface model [START_REF] Körner | Lattice Boltzmann model for free surface flow for modeling foaming[END_REF] based on the liquid volume fraction concept is adopted, and it is noted as FS1 model hereafter. As shown in Figure 1, the free-surface is represented by a layer of lattice cells, which are called interface cells here, along the path of the free-surface profile on the Cartesian background grid. Depending on the volume fraction (noted as α) of the liquid component in each cell, three types of cell flag can be distinguished, which means α = 1 for liquid, 0 < α < 1 for interface and α = 0 for gas.

In the FS1 model, the interface cells must form a single-layered continuous line, preventing the liquid cells and the gas cells from direct contact.

Reconstruction of the distribution functions on the interface

The missing distribution functions on the free-surface boundary need to be reconstructed, since no calculation are carried out in gas cells. Korner et al. [START_REF] Körner | Lattice Boltzmann model for free surface flow for modeling foaming[END_REF] provided a reconstruction method for interface cells, as shown in Figure 1, which is expressed as

f pc i (x -ξ i ∆t, t) = f eq i (x, t) + f eq ī (x, t) -f pc ī (x, t) , ∀i such that n • ξ i < 0, ( 7
)
where n is the normal vector of the interface, ī denotes the opposite direction of i, and f pc i stands for the postcollision distribution function. f pc i (xξ i ∆t, t) is temporarily stored in the gas cell xξ i ∆t and transmitted to f i (x, t + ∆t) after streaming.

Interface evolution and mass conservation

The free-surface movement is realized with the adaption of the positions of interface cells. Numerically, this evolution is driven by the mass exchanges between cells. Letting M = αρV 0 be the cell mass with V 0 being the fixed cell volume which can be set as unit, the cell mass change in one time step can be computed as

M (x, t + ∆t) = M (x, t) + i ∆M i (x, t), (8)
where ∆M i is the net mass increment in direction i due to streaming, which is computed as

∆M i (x, t) = C i f pc ī (x + ξ i ∆t, t) -f pc i (x, t) . (9)
The value of the coefficient C i is related to the involved two cells and is given in Table 1.

Table 1: Value of the coefficient C i in the mass evolution equation [START_REF] Ginzburg | Lattice Boltzmann model for free-surface flow and its application to filling process in casting[END_REF].

f lag (x) f lag (x + ξ i ∆t) C i gas gas 0 interface 0 interface gas 0 interface 1 2 [α(x) + α(x + ξ i ∆t)] liquid 1 liquid interface 1 liquid 1
When an interface cell gains so much mass during the streaming process that its volume fraction temporarily passes the limit of α = 1, this filled cell is forced to transform into a liquid cell, as shown in reaction, its neighboring gas cells must become interface cells in order to ensure a continuous interface layer.

The transformation from the interface cell to liquid cell generates a sudden cell mass drop of the excess mass component above the limit, i.e. ∆M ex = (α -1) ρV 0 , which will be distributed to the neighboring interface cells. Such distribution is based on a weighting function associated with the interface normal vector, which leads the excess mass to be transferred in the interface motion direction [START_REF] Thürey | Interactive free surface fluids with the lattice Boltzmann method[END_REF]. It is important to note here that by saying in the interface motion direction, we mean the trend of distributing larger percentage of the excess mass along the normal vector direction, since the normal vector as calculated in the following formula is not exactly the real interface motion direction

n (i, j) = ∇α (i, j) ∇α (i, j) , (10)
where the two components of the volume fraction gradient ∇α (i, j) is approximately computed as

∇α ≈ α (i + 1, j) -α (i -1, j) 2∆x , α (i, j + 1) -α (i, j -1) 2∆x . ( 11
)
Marking the position of the filled interface cell as x, the interface cell at position x + ξ i ∆t and time t + ∆t will receive an amount of mass of

∆M receive (x + ξ i ∆t) = w * i i w * i ∆M ex (x) , (12)
where the weighting coefficient w * i is calculated as

w * i =      n • ξ i , if n • ξ i > 0, 0, if n • ξ i 0. (13)
Besides, the newly formed interface cells are initialized with the equilibrium distribution functions that are computed from its surrounding cells as

               ρ (x, t) = i ρ (x + ξ i ∆t, t) i 1 , ∀f lag (x + ξ i ∆t, t) = gas, u (x, t) = i u (x + ξ i ∆t, t) i 1 , ∀f lag (x + ξ i ∆t, t) = gas. (14)
Similarly, for an emptied interface cell, of which the volume fraction temporarily becomes α < 0, the excess mass is ∆M ex = -αρV 0 , and the weighting coefficient for distributing the excess mass is computed as

w * i =      -n • ξ i , if n • ξ i < 0, 0, if n • ξ i 0. ( 15
)
Regarding the mass conservation, by calculating the mass change during streaming as shown in Equation ( 9) and by adopting the aforementioned distribution rule for the excess cell mass due to the cell flag change, the mass loss in one cell is always balanced by the mass gained in its neighbors. Hence, the total mass of the system is conserved by construction.

Regularized lattice Boltzmann method

The lattice Boltzmann equation with force term (1) can be regularized onto the Hilbert space H N spanned by the leading N Hermite polynomials, through a truncated Hermite expansion on each term [START_REF] Zhang | Efficient kinetic method for fluid simulation beyond the Navier-Stokes equation[END_REF]. Taking the distribution function f i as an example, it is approximated by an Hermite expansion truncated at N th -order as

f i (x, t) ≃ fi (x, t) = w i N n=0 1 n! a (n) (x, t) : H (n) (ξ i ), (16)
where the operator ":" denotes full contraction of two tensors.

a (n) (x, t) = i f i (x, t) H (n) (ξ i ) is the n th -order
expansion coefficient, and H (n) (ξ i ) is the n th -order Hermite basis, of which the leading several ones are given

as              H (0) (ξ i ) = 1,
H (1) (ξ i ) = 1 c s ξ i , H (2) (ξ i ) = 1 c 2 s ξ i ⊗ ξ i -I. (17)
Thanks to the orthogonal nature of the Hermite polynomials, the regularized distribution function fi (x, t) and the original one f i (x, t) have the same expansion coefficients up to order

N , i.e. â(n) (x, t) ≡ a (n) (x, t) , ∀n N .
The analysis in [START_REF] Shan | Kinetic theory representation of hydrodynamics: a way beyond the Navier-Stokes equation[END_REF] revealed that a truncated series of the Hermite expansion at N = 2 is adequate for retaining hydrodynamic macroscopic variables, hence we choose to regularize the lattice Boltzmann equation onto the H 2 space. In addition, as indicated in [START_REF] Li | On the coupling of a direct-forcing immersed boundary method and the regularized lattice Boltzmann method for fluid-structure interaction[END_REF], it is important to include the 1 st -order reconstruction term when adopting the force scheme of Guo et al. [START_REF] Guo | Discrete lattice effects on the forcing term in the lattice Boltzmann method[END_REF] in order to incorporate correctly the body force in the fluid flow. Noting that the equilibrium distribution in Equation ( 3) and the force term in Equation ( 5) are already in the H 2 space. One can rewrite the lattice Boltzmann equation (1) as

f i (x + ξ i ∆t, t + ∆t) = f eq i (x, t) + 1 - ∆t τ f neq i (x, t) + ∆tF i (x, t) , (18)
where f neq i is the regularized non-equilibrium component of the distribution function and can be computed as

f neq i (x, t) = w i 2 n=0 1 n! a (n) neq (x, t) : H (n) (ξ i ) (19)
with a

(n) neq (x, t) = i f neq i (x, t) H (n) (ξ i ).
The non-equilibrium distribution function is extracted as

f neq i (x, t) = f i (x, t) -f eq i (x, t) . ( 20
)
Following the idea of [START_REF] Shan | Kinetic theory representation of hydrodynamics: a way beyond the Navier-Stokes equation[END_REF], the regularized lattice Boltzmann equation ( 18) is able to provide equivalent macroscopic informations in the hydrodynamic Navier-Stokes level as the original lattice Botlzmann equation

(1) does for simulating isothermal weakly-compressible fluid flows. In other words, the regularization procedure can help filter out the undesired non-hydrodynamic components which are invisible in the Navier-Stokes level but can be a source of numerical instabilities.

4. An improved reconstruction method of the distribution functions on the interface

Discussion on the original reconstruction method

In [START_REF] Körner | Lattice Boltzmann model for free surface flow for modeling foaming[END_REF] it is suggested to employ the reconstruction operation [START_REF] Coreixas | Recursive regularization step for high-order lattice Boltzmann methods[END_REF] to half of the distribution functions on the interface (except f 0 ), i.e. ∀i such that n • ξ i < 0. An explanation was given in [START_REF] Körner | Lattice Boltzmann model for free surface flow for modeling foaming[END_REF] that the hydrodynamic stress in the normal direction σ normal on a small free-surface area within the x cell can be computed from the eight post-collision distribution functions which penetrate the free-surface during the streaming

σ normal = n • i,n•ξi>0 f pc i (x, t) (e i -u) (e i -u) + i,n•ξi 0 f pc i (x -ξ i ∆t, t) (e i -u) (e i -u) , = n • i,n•ξi 0 f eq i + f eq ī (e i -u) (e i -u) , = n • i f eq i (e i -u) (e i -u) , = n • p G I. (21)
By reconstructing half of the distribution functions, the gas pressure p G imposed on the interface is balanced by the hydrodynamic force from the liquid side. However, this strict requirement may encounter some dilemmas in certain situations. One of the possible dilemma can be found at the edge of a convex free-surface profile, as shown in Figure 3.

The interface cell at the corner of a water column has five missing distribution function, which are f 3 , f 4 , f 6 , f 7 and f 8 . The criterion of n • ξ i < 0 only involves f 3 , f 4 and f 7 . In order to include the other two, one has to modify the criterion as ∀i such that n • ξ i 0 or f lag (xξ i ∆t) is gas. This new criterion may lead to an over-reconstructed set of distribution functions on the interface. Another dilemma occurs at the three-phase point, as shown in Figure 4. At an interface cell close to a free-slip solid wall, the distribution functions f 2 , f 5 and f 6 need to be computed by the specular boundary condition treatment for free-slip walls [START_REF] Lim | Application of lattice Boltzmann method to simulate microchannel flows[END_REF]. However, the original free-surface model requires f 3 , f 4 , f 6 and f 7 to be reconstructed and f 1 , f 2 , f 5 and f 8 to be streamed. In this way, the two boundary conditions conflict, which means that if the free-surface boundary condition is satisfied by the reconstruction procedure, then the free-slip condition cannot be guaranteed, and vice versa. To make matters worse, sometimes the number of distribution functions that can be streamed from neighboring cells are very limited, due to the incomplete neighborhood near the wall, which may cause greater confusion in the priority sequence among streaming, the specular rule and the reconstruction.

A new reconstruction method based on regularization

To provide a countermeasure for the aforementioned dilemmas of the original model, a new way of reconstructing the distribution functions at the free-surface is introduced here, which is labeled as FS2 model in the rest part of the paper. The idea is inspired by the work of Latt et al. [START_REF] Latt | Straight velocity boundaries in the lattice Boltzmann method[END_REF], which is originally proposed for straight velocity boundaries. The new reconstruction method requires to compute the density ρ and the fluid velocity u at first. Once they are known, the equilibrium distribution function can be computed from its definition (3), and the non-equilibrium part can be approximated from [START_REF] Latt | Straight velocity boundaries in the lattice Boltzmann method[END_REF]

f neq i ≈ τ ρ c 2 s ξ i ξ i -c 2 s I : ∇u. (22)
Afterwards, the distribution function can be reconstructed as

f i = f eq i + f neq i . (23)
In the present work, the fluid density at the interface is chosen to be ρ = p G /c 2 s , and the fluid velocity is predicted through an inverse distance extrapolation as

u(x) = i w † i u(x i ) i w † i , (24)
where x i stands for a liquid cell within the neighborhood area of a searching radius of r around the interface cell at x, and w † i is the corresponding weight as

w † i = 1 d i
, in which d i is the distance between x and x i .

Consequently, the velocity gradient in Equation ( 22) can be obtained using a classic finite difference scheme.

It is important to note here that the reconstruction ( 22) of the non-equilibrium distribution function can be explained through a Hermite regularization [START_REF] Li | On the coupling of a direct-forcing immersed boundary method and the regularized lattice Boltzmann method for fluid-structure interaction[END_REF]. Following the Chapman-Enskog analysis [START_REF] Chapman | The mathematical theory of non-uniform gases: an account of the kinetic theory of viscosity[END_REF], the distribution function and the derivative operators in the LBE can be expanded into a series with respect to a small parameter

ε as          f i = ε 0 f (0) i + ε 1 f (1) i + ε 2 f (2) i + O ε 3 , ∂ t = ε 1 ∂ t1 + ε 2 ∂ t2 + O ε 3 , ∇ = ε 1 ∇ 1 + O ε 2 , ( 25
)
where t 1 = ε 1 t and t 2 = ε 2 t. One can obtain the 2 nd -order moment Π (1) of f

Π

(1) = -τ ρc 2 s ∇ 1 u + (∇ 1 u) T , (26)
By assuming

f neq i ≃ εf (1)
i , the 2 nd -order expansion coefficient of f neq i can be obtained as

a (2) neq = i 1 c 2 s ξ i ⊗ ξ i -I f neq i ≃ ε c 2 s Π (1) = -τ ρ ∇u + (∇u) T . ( 27
)
Hence f neq i can be regularized through

f neq i = w i 2 n=0 1 n! a (n) neq (x, t) : H (n) (ξ i ), ≃ -w i τ ρ c 2 s ξ i ξ i -c 2 s I : ∇u. (28)
The analysis above shows that Equation ( 22) approximately restores the regularization procedure [START_REF] Zhang | Efficient kinetic method for fluid simulation beyond the Navier-Stokes equation[END_REF] of projecting f neq i onto the H 2 space. In this sense the present reconstruction method is expected to be more consistent with the RLBM scheme.

Algorithm

The proposed reconstruction method is to substitute the original one, while the rest ingredients of the free-surface model stay unchanged, such as cell mass evolution and interface update. The present FS2-RLBM algorithm is summarized in the following table, where the cell flag in the box at the beginning of each item indicates the target cell type of the corresponding command:

Algorithm: FS2-RLBM Require: f i (x, t n ), ρ (x, t n ), u (x, t n ), F (x, t n ), α (x, t n ), m (x,

Numerical results and discussions

Lid driven cavity flow

The lid-driven cavity flow with a Reynolds number of Re = 1000 is firstly studied to assess the accuracy and efficiency of the RLBM. Simulations are carried out by the SRT-LBM and RLBM, with 65 × 65, 129 × 129, 257 × 257 and 513 × 513 lattices, and the relaxation time is kept at τ /∆t = 0.506 for all cases. The treatment for boundary condition in [START_REF] Latt | Straight velocity boundaries in the lattice Boltzmann method[END_REF] is adopted.

Table 3 gives the consumed CPU times and the L 2 -error of the velocity components on the centerlines,

Error = (u -u Ref ) 2 u 2 Ref , (29)
where the reference data is chosen as the benchmark results of spectral methods [START_REF] Botella | Benchmark spectral results on the lid-driven cavity flow[END_REF].

One can observe from the L 2 -errors that

• with the same lattice resolution, RLBM results show an overall slightly better accuracy than the SRT-LBM results,

• RLBM consumes about 1.3 times CPU time than SRT-LBM. However with the same CPU time, i.e.

RLBM with 209 × 209 lattice and SRT-LBM with 257 × 257 lattice, RLBM can provide as accurate results

as SRT-LBM yet with a smaller lattice,

• the L 2 -error slopes in Figure 5 show an overall 2 nd -order convergence for both schemes, which means the regularization does not affect the convergence of LBM.

Furthermore, from the pressure contours in Figure 6, it can be clearly observed that the SRT-LBM results have some severe chessboard-liked pressure (density) fluctuations introduced from the singularity at the upper corners, while RLBM succeeds in filtering out these undesired numerical noises.

Viscous standing wave

In this part, a periodic standing wave of wavelength λ = 2.0 m and wave steepness ε = 0.05 is considered.

The wave steepness is the ratio of the wave height 2A to the wavelength λ, i.e. ε = 2A/λ where A is the wave amplitude. It is pointed out in [START_REF] Colagrossi | Theoretical analysis and numerical verification of the consistency of viscous smoothed-particle-hydrodynamics formulations in simulating free-surface flows[END_REF][START_REF] Lighthill | Waves in fluids[END_REF] that the potential theory gives an approximate solution to this problem, and the velocity potential reads

ϕ (x, y, t) = - Ag ω cosh [k (y + H)] cosh (kH) cos (kx) cos (ωt) , ( 30
)
where g is the gravity acceleration, k = 2π/λ is the wave number, H is the water depth of calm water, and consequently ω = gk tanh (kH) is the angular frequency, and T = 2π/ω is the wave period. One can observe that Equation ( 30) not only is periodic in the x-direction, but also satisfies symmetry at x = iλ 2 , ∀i ∈ Z. Hence the computation domain is chosen to be of width L = λ/2 with symmetric boundary conditions used on the two vertical walls, as shown in Figure 7. At time t = 0, the free surface profile is a horizontal line at y = 0, and the pressure field is assumed to be hydrostatic. Besides, the velocity field is initialized as the gradient of the velocity potential ∇ϕ 0 , where ϕ 0 = ϕ (x, y, 0). After simple computation, one can know that the maximum velocity appears at t = 0, which is u max = Agk ω . Consequently, the Reynolds number is defined as Re = λu max ν .

If the fluid is inviscid, the total kinetic energy of the system is conserved in time, and the Reynolds number is infinite. However, for a viscous case as discussed here, the Reynolds number is of finite value, and the kinetic energy dissipates due to the viscous effect. For the viscous standing wave at Reynolds number Re = 100, the snapshots of the flow field at some typical time instants are shown in Figure 8, where the results are obtained the present FS1-RLBM scheme with a lattice resolution of 80 in one wave height, i.e. 2A = 80∆x.

Generally speaking, the RLBM works well with the FS1 free-surface model and provides reasonable results.

However, small horizontal velocities are observed in the area close to the free-surface at t = 0.25T and t = 0.75T when the velocity in the fluid domain is supposed to be zero, and the free-surface is not perfectly horizontal at t = 0.50T and t = 1.00T . One of the reasons is that a slight difference of the time period exists between the numerical and the analytical solutions, another reason may lie in the discrete nature of the interface evolution mechanism of the adopted free-surface model which might introduce small perturbations from the staircase-liked interface layer.

For a more precise validation, we extract the evolution of the total kinetic energy of the system in Figure 9, where the analytical solution of the kinetic energy damping is given in [START_REF] Colagrossi | Theoretical analysis and numerical verification of the consistency of viscous smoothed-particle-hydrodynamics formulations in simulating free-surface flows[END_REF][START_REF] Lighthill | Waves in fluids[END_REF] as

E t = λA 2 g 8 e -4νk 2 t [1 + cos (2ωt)] . (31)
As one can observe, for the adopted four lattice resolutions (2A/∆x = 20, 40, 60 and 80), the results converge as the grid is refined. The numerical results have a time period very close to the analytical solution, whereas the damping rate of the kinetic energy is a little higher than in the analytical prediction. It is necessary to note that the analytical solution, which is obtained from the linear theory, has its own error at the wave steepness ε = 0.05. This error would be reduced for smaller amplitude waves. However, in order to get converged numerical results

t / with an adequate lattice resolution, the calculation would become extremely costly, therefore we stop at the present wave steepness. In addition, the relative kinetic energy value is close to 0 at t = 0.25T and t = 0.75T , which suggests that the aforementioned velocity errors near the surface are acceptable. Moreover, the numerical results of the SRT-LBM and RLBM schemes are compared to each other, where the same lattice resolution 2A/∆x = 80 is adopted.

Figure 10a and 10b show the time evolution of the total kinetic energy and the water depth at the left (x = 0) and the right (x = L) boundaries of the computation domain. These two results are nearly coincident, which indicates that the regularization procedure brings no extra numerical dissipation for free-surface flows at small Reynolds numbers.

Dambreak flow

At last, we consider a more violent dambreak flow test-case. Figure 11 gives the initial configuration of the simulation, which is of the same size as in the experiment conducted by Lobovsky et al. [START_REF] Lobovskỳ | Experimental investigation of dynamic pressure loads during dam break[END_REF]. A 600 mm (L) × 300 mm (H) water column is reserved to the right side in a 1610 mm × 900 mm water tank. The water column collapses due to gravity and impacts on the left vertical wall. Four pressure sensors, marked as P1-P4, are distributed vertically on the left wall to record the impact pressure signals (at heights 3 mm, 15 mm, 30 mm and 80 mm). Besides, the water front and the water level at H1-H4 positions are examined (at 300 mm, 865 mm, 1114 mm and 1362.5 mm downstream from the right wall).

In the present work, the water column is discretized by a 600 × 300 lattice and initialized as hydrostatic.

The discrete time step is ∆t = 1 × 10 -5 s. The gravity acceleration is g = -9.81 m/s 2 . The effective Reynolds number of this first simulation is Re = Av/ν = 2.0 × 10 4 , where the characteristic velocity v is the predicted water front speed √ 2gH, and the characteristic length A is chosen to be the same as in the experiment which is the distance between the water column and the left wall. Since the viscous effect is not dominant, a half-way specular reflecting boundary condition [START_REF] Lim | Application of lattice Boltzmann method to simulate microchannel flows[END_REF] is adopted to mimic free-slip solid walls.

Reynolds number dependence

One may notice that the Reynolds number in the numerical simulation is smaller than that in the experiment which is as high as Re = 3.8 × 10 6 . Hence we will firstly demonstrate that the present Reynolds number is high enough to obtain a converged result, and that higher Reynolds number only affect the numerical stability rather than accuracy. In this context, simulations by FS1-RLBM with Re = 1.0 × 10 4 , 4.0 × 10 4 and 8.0 × 10 4 are carried out. Their water front position and water level evolutions are shown in Figure 12.

The RLBM result does not change much when the Reynolds number is increased by one order of magnitude, which indicates that even Reynolds number 1.0 × 10 4 is large enough that the viscous effect is relatively too t ( g / H ) 0.5 small to be observed. With a Reynolds number starting from Re = 4.0 × 10 4 , small oscillations can be seen on the RLBM curves, which indicates an increasing numerical instability. We shall come to this during the comparison to SRT-LBM results later. With Re = 2.0 × 10 4 , two sets of simulations are carried out at first, i.e. FS1-SRT-LBM and FS1-RLBM. Some snapshots are illustrated in Figure 13. We would like to note here that the pressure field is badly disturbed when the falling water jet hits on the free-surface. As Marrone et al. [START_REF] Marrone | Analysis of free-surface flows through energy considerations: Single-phase versus two-phase modeling[END_REF] pointed out, the single-phase free-surface approach is valid for the liquid flow up to when gas is entrapped within the liquid. When the latter occurs, the assumption of neglecting the effect of gas motion is questionable. In this context, we shall focus on the time period before such hitting behavior occurs.

Comparison of FS1-SRT-LBM and FS1-RLBM with

It is observed that the simulated water front moves a little faster than the experimental one. Two reasons may contribute to this time difference. The first one is the effect of the gate motion, namely in the experiment, a boundary layer is developed near the lifted gate, which will hold the water front for a while and thus make it slower than the numerical results. The other reason may be the use of free-slip boundary condition in the simulation. In this aspect, the study of Marrone et al. [START_REF] Marrone | δ-SPH model for simulating violent impact flows[END_REF] provides some evidences, where both free-slip boundary and non-slip boundary conditions are tested in terms of the arriving time when the water front hits the wall. Their results confirmed that the water impact did occur earlier in the free-slip case than in the non-slip cases, even when the Reynolds number for the non-slip case is as large as Re = 10000. This indicates that adopting the non-slip boundary condition may be closer to the physics nature, but it requires accurately discretizing the boundary layer and thus needs a very fine mesh close to the boundary, which is considered to be expensive for a simulation. Besides, the SPH results in [START_REF] Marrone | δ-SPH model for simulating violent impact flows[END_REF] also showed that the time difference of the water impact is reduced when the Reynolds number is increasing for the non-slip case. Based on this consideration, we think that the present practice of applying the free-slip boundary condition is acceptable.

The pressure signals are shown in Figure 14, where the numerical curves are based on raw pressure data without any filtering operation. As previously discussed, the numerical water front has a larger speed than the experimental one, hence the first peaks of the numerical pressure signals appear earlier than the experimental ones for the same reason. It is also observed from the snapshots in Figure 13 that a smooth free-surface profile and pressure field are obtained by the RLBM while the smoothness cannot be maintained by the SRT-LBM. This is confirmed in Figure 14 with more details. The first peaks of the SRT-LBM pressure goes up to extremely high values at the four pressure sensors, whereas the RLBM pressure peaks are similar with the experimental ones. Furthermore, the SRT-LBM pressure is oscillating wildly, while this severe fluctuation is remarkably reduced by the RLBM scheme. From this comparison, the regularization in the RLBM scheme is verified to be able to largely improve the numerical stability of the LB method in high Reynolds number applications. Four combinations of FS1/FS2 and SRT-LBM/RLBM are then tested with Re = 3200. We chose to reduce the Reynolds number here in order to let the influence of the free-surface models show up when the instability of SRT-LBM does not contribute too much. The dambreak flow is reproduced successfully by the four schemes, and the numerical flows are very similar to the experimental one although the Reynolds number is largely reduced in the simulation (we only extract the snapshots at t = 159.9 and 862.3ms for demonstration as shown in Figure 15). From the pressure contours, one can observe that the present model (FS2) performs as well as the original one (FS1), if not better, and it is the regularization procedure in the RLBM which dominates in smoothing the pressure field.

However, if we zoom in a little by extracting the pressure signals at positions P1-P4, the advantage of the present model can be seen. Figure 16 gives the pressure signals recorded with the four schemes, where the results of both free-surface models with the standard SRT-LBM are shown in the left panels, and those of the two free-surface models with the RLBM are plotted in the right panels.

In association with the standard SRT-LBM, the present FS2 model provides results that are close to that of the FS1 model: the pressure oscillations are in the same range. As a contrast, in the RLBM framework, although the pressure oscillations have already been remarkably reduced by the RLBM, the present FS2 model is able to achieve even better results. Indeed, this improvement is not that prominent in the pressure curves, but this is due to the fact that the number of interface cells that are caught in the aforementioned dilemmas of the original model, is of a relatively small percentage during the whole calculation. Fairly speaking, one is still able to conclude that the present FS2 model is more consistent with the RLBM.

Conclusions

In this paper, the applicability of the lattice Boltzmann method for free-surface flows is investigated. The single-phase free-surface approach [START_REF] Körner | Lattice Boltzmann model for free surface flow for modeling foaming[END_REF] is adopted to avoid the difficulty of the large density ratio for multi-phase models, and that it can keep a sharp interface with only one cell. As an improvement to the original model, we adopted several techniques to enhance the robustness of the scheme, including the cell mass distribution rule during flag change and the countermeasure for adjacent filled and emptied interface cells.

In order to alleviate the numerical instability that the standard single-relaxation-time lattice Boltzmann method (SRT-LBM) suffers from at high Reynolds numbers, the projection-based regularization procedure [START_REF] Zhang | Efficient kinetic method for fluid simulation beyond the Navier-Stokes equation[END_REF] is introduced to the free-surface lattice Boltzmann solver. Inspired by the similar spirit of the Hermite regularization on the non-equilibrium distribution functions [START_REF] Latt | Straight velocity boundaries in the lattice Boltzmann method[END_REF], we proposed to calculate the equilibrium distribution function of the interface cell based on the extrapolated macroscopic variables and reconstruct the non-equilibrium one from the local velocity gradient. The new reconstruction method is expected to not only ensure the balance of the hydrodynamic force and the gas pressure in the normal direction on the interface, but also take into account the viscous effect. Moreover, it eliminates the dilemmas of over-reconstruction or boundary condition conflicts that the original model may encounter.

Through several numerical test-cases, the following conclusions are drawn:

• the RLBM is verified to be able to remarkably reduce the numerical noises generated by the singularity at the lid-cavity corners and the free-surface;

• the RLBM does not decrease the convergence order, and it brings no extra energy dissipation compared with the SRT-LBM;

• although RLBM consumes about 30% CPU time more than the SRT-LBM with the same lattice resolution, it can achieve an equivalent accuracy with a smaller lattice within the same CPU time;

• the proposed regularization-based reconstruction method of the distribution functions on the interface is more consistent with the RLBM, and it can further reduce the undesired numerical noises in the pressure field.

Generally speaking, the proposed RLBM free-surface solver with the regularization-based reconstruction method is an efficient and robust alternative for modeling free-surface flows, with the major benefit of a largely improved numerical stability and consequently a better applicability to large deformation free-surface flows.

Figure 1 :

 1 Figure 1: Lattice cells are classified into liquid, interface and gas cells based on its volume fraction. The distribution functions in the direction of ∀n • ξ i < 0 (i = 1, 2, • • • , 8) are reconstructed, where n is the normal unit vector at the free-surface.

Figure 2 .Figure 2 :

 22 Figure 2: When a filled interface cell at position x turns into a liquid cell, the excess cell mass is distributed to the neighboring interface cells at position x + ξ i ∆t that satisfy n • ξ i > 0. The distribution is based on a weighting coefficient related to the interface normal vector n.

Figure 3 :

 3 Figure 3: Over-reconstructed distribution functions for interface cells on a convex free-surface profile.

Figure 4 :

 4 Figure 4: Conflict of boundary condition treatments at the three-phase point.

Figure 5 :

 5 Figure 5: L 2 -error of the velocity components along the centerlines.

Figure 6 :

 6 Figure 6: Pressure fields of the SRT-LBM and RLBM results.

Figure 7 :

 7 Figure 7: Initial state of the viscous standing wave problem.

Figure 8 :

 8 Figure 8: Snapshots of the viscous standing wave flow at Re = 100, where the color indicates the relative velocity magnitude u/umax. The results are obtained from the FS1-RLBM scheme with a lattice resolution of 2A/∆x = 80.

Figure 9 :

 9 Figure 9: Decay of the total kinetic energy in the viscous standing wave problem.

Figure 10 :

 10 Figure 10: Time evolution of the total kinetic energy of the system and the water depth at both the left (x = 0) and the right (x = L) boundaries of the computation domain, based on the results from FS1-SRT-LBM and FS1-RLBM with a lattice resolution of 2A/∆x = 80.

Figure 11 :

 11 Figure 11: Initial configuration of the dambreak problem (the unit length is 1mm)

t ( 8 /

 8 H ) 0.5

Figure 12 :

 12 Figure 12: Evolution of water front and water level at H1-H4, obtained by the FS1-RLBM with Re = 1.0 × 10 4 , 2.0 × 10 4 , 4.0 × 10 4 and 8.0 × 10 4 .

  Re = 2.0 × 10 4

  Results from the experiment (left), the FS1-SRT-LBM (middle) and the FS1-RLBM (right) at t = 0.0ms. Results from the experiment (left), the FS1-SRT-LBM (middle) and the FS1-RLBM (right) at t = 159.9ms. Results from the experiment (left), the FS1-SRT-LBM (middle) and the FS1-RLBM (right) at t = 276.6ms. Results from the experiment (left), the FS1-SRT-LBM (middle) and the FS1-RLBM (right) at t = 373.3ms. Results from the experiment (left), the FS1-SRT-LBM (middle) and the FS1-RLBM (right) at t = 449.9ms. Results from the experiment (left), the FS1-SRT-LBM (middle) and the FS1-RLBM (right) t = 573.3ms. Results from the experiment (left), the FS1-SRT-LBM (middle) and the FS1-RLBM (right) t = 862.3ms. Results from the experiment (left), the FS1-SRT-LBM (middle) and the FS1-RLBM (right) t = 1023.3ms.

  Results from the experiment (left), the FS1-SRT-LBM (middle) and the FS1-RLBM (right) t = 1166.6ms.

Figure 13 :

 13 Figure 13: Snapshots of free-surface profiles and pressure fields at time 0.0, 159.9, 276.6, 373.3, 449.9, 573.3, 862.3, 1023.3, 1166.6 ms. The experimental results (left), the FS1-SRT-LBM (middle) and the FS1-RLBM (right) are shown (reprinted with permission from Elsevier).

Figure 14 :

 14 Figure 14: Evolution of the pressure signals at P1, P2, P3 and P4 positions.

5. 3 . 3 .

 33 Comparison of FS1/FS2 in combination with SRT-LBM/RLBM with Re = 3200

  t = 159.9ms: Experiment (left), FS1-SRT-LBM (middle left), FS2-SRT-LBM (middle), FS1-RLBM (middle right), and FS2-RLBM (right).

  t = 862.3ms: Experiment (left), FS1-SRT-LBM (middle left), FS2-SRT-LBM (middle), FS1-RLBM (middle right), and FS2-RLBM (right).

Figure 15 :Figure 16 :

 1516 Figure 15: Snapshots of the free-surface profiles and the pressure fields at time t = 159.9 and 862.3ms. The exhibited results are from the experiment (left), the FS1-SRT-LBM (middle left), the FS2-SRT-LBM (middle), the FS1-RLBM (middle right) and the FS2-RLBM (right).

  t n ) and f lag array from t n time-step 1. liquid interface Compute f eq i (x, t n ) with Equation (3) and F i (x, t n ) with Equation (5),

	2. liquid interface Extract f neq i	(x, t n ) by Equation (20),
	3. liquid interface Regularize f neq

i by Equation (19), 4. liquid interface Carry out the collision and streaming processes by Equation (18) 5. liquid Adapt ρ x, t n+1 , u x, t n+1 by Equation (6), 6. interface Set ρ x, t n+1 , and extrapolate u x, t n+1 by Equation (24), 7. interface Reconstruct the distribution functions by Equation (23), 8. interface Adapt m x, t n+1 by Equation (8) and modify α x, t n+1 , 9. interface Check α x, t n+1 and determine f lag x, t n+1 array, 10. interface Initilize the newly formed interface cells by Equation (14), 11. interface For the cells involved in the flag change, modify m x, t n+1 by Equation (12) and adapt α x, t n+1 , 12. Goto Step.1 for the next time-step.

Table 3 :

 3 L 2 -error of the velocity components on the centerlines and the consumed CPU times. The listed CPU times are relative values with respect to that of the 65 × 65 SRT-LBM case.

	Numerical scheme Lattice resolution Error u Error v	t CPU
		65 × 65	5.1765	9.9370	1.00
	SRT-LBM	129 × 129 257 × 257	0.5637 0.1441	0.9638 0.2263	14.08 199.49
		513 × 513	0.0505	0.0663 3008.21
		65 × 65	5.2084	9.8242	1.27
		129 × 129	0.5753	1.0117	18.93
	RLBM	209 × 209	0.1843	0.1979	199.46
		257 × 257	0.0734	0.1070	266.12
		513 × 513	0.0389	0.0487 4077.48

Declaration of interests

☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

☐The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: