A LATTICE BOLTZMANN APPROACH FOR THREE-DIMENSIONAL TSUNAMI SIMULATION BASED ON THE PLIC-VOF METHOD

Free surface flow problems occur in numerous disaster simulations, such as tsunamis inland penetration in urban area. Simulation models for these problems have to be non-hydrostatic and three-dimensional because of the strong non-linearity and higher-order physical phenomena. Despite all the progress in the modern computational fluid dynamics, such simulations still present formidable challenges both from numerical and computational cost point of view. The lattice Boltzmann method (LBM) has been attracting attention as an alternative fluid simulation tool to overcome the problems. In current study, LBM for three-dimensional tsunami simulations is developed which are coupled with the piecewise linear interface calculation with the Volume of Fluid (VOF) approach. This model is for an efficient three-dimensional tsunami simulation by a one-fluid formulation, where the lattice Boltzmann equation is assigned to solve for a single virtual fluid. Various benchmark problems are also carried out to validate the utility of the proposed models in term of coastal engineering.


INTRODUCTION
Three-dimensional free surface modelling is required to simulate tsunami inundation in urban areas with high accuracy.Despite offering a powerful tool to investigate tsunamis, such as gridbased (Choi et al., 2007;Hu et al., 2012;Wroniszewski et al., 2014;Wu et al., 2014) or meshfree methods (Koshizuka and Oka, 1996;Khayyer and Gotoh, 2009;Asai et al., 2012;Sun et al., 2015), three-dimensional tsunami simulations have been successfully applied only in limited computing environment.The overall implicit scheme of most Navier-Stokes solvers has to solve one or more linear system of equations in each time step (i.e., the pressure Poisson equation), this implies that the capillary effects could acts as a formidable computational barrier.
The lattice Boltzmann method (LBM) has been attracted attention as a novel computational fluid dynamics.LBM originates from the lattice gas automata for kinematic molecular dynamics (McNamara and Zanetti, 1988).In contrast to other solvers, LBM can be considered an alternative approach to solve macroscopic fluid equations.The key features of the LBM are as follows: 1. Fully explicit method in time integration, which means LBM does not have to solve the pressure Poisson equation, iteratively 2. Easy implementation with high-performance computing (Janßen and Krafczyk, 2011;Obrecht et al., 2011;Rinaldi et al., 2012;Obrecht et al., 2013;Safi and Turek, 2016) due to the algorithmic operations and data locality In this regard, we believe LBM has advantages as an alternative approach to perform three-dimensional tsunami simulations efficiently.In current study, we develop a three-dimensional free surface model by LBM based on the piecewise linear interface calculation (PLIC) method to enhance the model accuracy and overcome the traditional model problems (Thürey and Rüde, 2009).For the validation of the free surface model, we present simulations of dam breaking and breaking waves benchmarks.

THE LATTICE BOLTZMANN METHOD
LBM evolved from methods for the simulation of gases that computed the motion of each molecule in the gas purely integer operations.In Hardy et al. (1976), there was a first attempt to perform fluid simulations with this approach.Motivated by this improvement (McNamara and Zanetti, 1988), they developed the algorithm that was actually called LBM by performing simulations with averaged floating point values instead of single fluid molecules.The important contribution to the basic LBM was the simplified collision operator with a single relaxation parameter.This collision model is well-known as the Bhatnagar-Gross-Krook (BGK) approximation (Bhatnagar et al., 1954), and was derived independently by Chen et al. (1992) and Qian et al. (1992).

The governing equation
The governing equation of LBM consists of two steps, the streaming and the collision step.In this study, a simple lattice BGK model was implemented for collision step.The lattice BGK equation, using a single relaxation time τ is given as: in which, ρ and u are the macroscopic fluid density and velocity, given by the 0th and 1st order moments of the distribution function f i : and τ determines the macroscopic kinematic fluid viscosity ν as: where ∆x and ∆t are the spacing lattice size and time step interval, respectively.For each of the lattice velocities e i , a floating point number f 0,...,q−1 , representing the fraction of particles moving with this velocity, needs to be stored.The virtual particle movement is restricted into a limited number of directions.In this thesis, We used the three-dimensional 19-speed (D3Q19) square lattice model for three-dimensional fluid simulation.Alternatives are models with 15 or 27 velocities.However, the latter one has no apparent advantages over the 19 velocity model, while the model with 15 velocities has a decreased stability.The D3Q19 model is, thus, usually preferable as it requires less memory than the 27 velocity model.The D3Q19 model with its lattice velocity vectors e 0,...,18 is shown in Fig. 1 The velocity vectors take the following values: (5) In the D3Q19 model there are particles not moving at all ( f 0 ), moving with speed 1 ( f 1,...,6 ) and moving with speed √ 2 ( f 7,...,18 ).In the following, a subscript of I will denote the value from the inverse direction of a value with subscript i. f i and f I are, thus, opposite distribution functions with inverse velocity vectors e I = −e i .

The equilibrium distribution function
The continuity equation in LBM introduces a compressibility effect which become critical for the simulation of the incompressible flows (Chen et al., 1992).In order to reduce the compressibility effect for the time dependent flows, we need the following incompressible LBM (He and Luo, 1997), a slight modification of the compressible one with respect to the density.
It is well understood that in incompressible flows the density is a constant, say ρ 0 , and the density fluctuation δρ, should be of the order O Ma 2 in the limit of Ma → 0. If we explicitly substitute ρ = ρ 0 + δρ into the equilibrium distribution function (Qian et al., 1992) and neglect the terms proportional to δρ (u/c) and δρ (u/c) 2 , which are of the order O Ma 3 or higher, then the equilibrium distribution function becomes: where c s is the virtual sound speed, w i is the weight coefficient given as; for D3Q19 square lattice model.
The above equation is the equilibrium distribution function of the incompressible LBM.By the same Chapman-Enskog expansion procedure, we can obtain the incompressible Navier-Stokes equations accurate up to the second order of Mach number O Ma 2 for continuity equation and third order of Mach number O Ma 3 for momentum equations as: Notice that the above system is still biased from standard Navier-Stokes equations by the compressibility effect related to Mach number.In order to reduce this compressibility effect, small Mach number is preferable, which means that ∆t has to be small with respect to ∆x.

Solid wall boundary condition
In LBM, there is no explicit distribution function f i streaming from the boundary into the bulk, we have to construct the unknown distribution functions.For no-slip solid wall boundaries, we imposed a second-order bounce-back scheme (Bouzidi et al., 2001).It slightly violates the conservation of mass, but because it is second order in space, an overall scheme of second order can be achieved.The modified bounce-back scheme is given by: where f I and f i are anti-parallel incoming (missing) and outgoing distribution functions, respectively.

Interface boundary condition
When we apply the single-phase free surface model into LBM, another problem be occurred in the interface cells.The distribution functions of empty gas cells are never accessed.Interface cells, however, always have empty cell neighbours (Fig. 2).Therefore, during the streaming step only distribution functions from fluid cells or other interface cells are streamed normally, while the distribution functions that would be read from empty cells need to be reconstructed with corresponding free surface boundary conditions at the interface cell (Körner et al., 2005;Thürey and Rüde, 2009).This boundary condition does not require additional constructs, such as ghost layers around the interface.They, thus, can be treated locally for each cell.An atmospheric pressure of ρ A = 1.0 is used, as this is also the reference density and pressure of the fluid.Moreover, it is assumed that the viscosity of the fluid is significantly lower than that of the gas phase, while having a higher density.Hence, the gas follows the fluid motion at the interface.In terms of distribution functions, this means that if at (x + ∆te i ) there is an empty cell: where u is the velocity of the interface cell at position x and time t.The pressure of the atmosphere onto the fluid interface is introduced by using ρ A for the density of the equilibrium distribution functions.Applying Eq. 11 to all directions with empty neighbour cells would result in a full set of distribution functions for interface cells.This boundary condition, thus, can be used for a multiphase flow (Anderl et al., 2014) taking proper density instead of density for atmosphere.However, to balance the forces on each side of the interface, the distribution functions coming from the direction of the interface normal n are also reconstructed.Thus, as shown in Fig. 2, if the distribution function f i would be streamed from an empty gas cell, or if n • e I > 0 holds, f i is reconstructed using Eq.11.

ONE-FLUID FREE SURFACE MODEL FOR TSUNAMI MODELLING
Free surface flow problems involve an air and a water phase, so that, technically, such flows can be classified as multi-phase flows.However, in tsunami simulation, the flow behaviour is dominated by the water phase, which means that it is not necessary to resolve the full aerodynamics inside the air phase.We, thus, developed a three-dimensional free surface model by one-fluid formulation.When the governing equations are solved on a fixed grid, using one set of equations for the whole flow field, the different fluids have to be identified in some way.This is generally done by using a colour function that takes different values in the different fluids.

Competing cells status
The Volume of Fluid (VOF) method (Hirt and Nichols, 1981;Youngs, 1982;Rudman, 1997) is the oldest approach to advect a colour function and -after many improvements and innovations -continues to be widely used for three-dimensional tsunami simulations due to the conservation of fluid mass.We also adopted VOF approach as a free surface model in LBM.The colour function c, which means fluid fraction value for each grid, is introduced as a new macroscopic value (Thürey and Rüde, 2009).This function is in the range [0, 1].
The piecewise linear interface reconstruction Hirt and Nichols (1981) proposed an interface reconstruction method, in which it is approximated by straight lines, parallel to the coordinate directions.To determine whether the interface should be horizontal or vertical, Hirt and Nichols (1981) have found the normal to the interface, using values of c in the neighboring cells, and selected the orientation of the interface depending on whether the normal was more closely aligned with the horizontal or the vertical axis.This approach is the simplest VOF model and commonly used in three-dimensional tsunami simulation as of today.Rudman (1997), however, suggested that the Hirt and Nichols (1981) method is not significantly more accurate.
Although the method of Hirt and Nichols (1981) perhaps did not improve significantly on the simple linear interface calculation (SLIC) method (Noh and Woodward, 1976), it nevertheless suggested that the key to improving the behaviour of the advection scheme was the reconstruction of the interface in each cell, using the value of the colour function in each cell, along with the value in the neighbouring cells.In the piecewise linear interface calculation (PLIC) method introduced by Youngs (1982), the interface is approximated by a straight-line segment in each cell.Despite this approach, thus, little more complicated in comparison with SLIC method, PLIC method is continue to be improved even now.For example, Multi-interface advection and reconstruction solver (MARS) has been proposed by Kunugi and Kino (2005) as an advanced PLIC algorithm for multiphase fluid analysis.
In this study, we implemented PLIC method as an interface reconstruction for more accurate free surface simulation by LBM.The orientation of the line is determined by the normal to the interface, which is found by considering the value c in both the cell under consideration and in the adjacent cells.Once the interface in each cell has been constructed, the fluxes from one cell to another are computed by geometric considerations.The result of the advection generally depends on the accuracy of the interface reconstruction; finding the normal accurately, therefore, becomes critical for PLIC methods.In three dimension, the interface shapes are approximated as a trapezoid as shown in Fig. 3, and the line segment function is given in each interface cell as follows: where n is the unit interface normal and α is the segment from the Cartesian origin.The calculation algorithms of the PLIC method, thus, are as follows: 1. Calculate the interface normal n 2. Calculate the segment parameter α

Interface normal
Although various algorithms have been proposed to determine the interface normal n, Parker and Youngs method (Pilliod and Puckett, 2004) is commonly used for its simplicity and accuracy.In this approach, the interface normal n is calculated by the gradient of the fluid fraction c direct, as follows: The gradient is discretized by the central difference scheme from the surrounding values: where c is the average fraction level of neighboring cells as: cy (x, y, z) = cz (x, y, z) = where w i, j is the weighting parameter.The following weights achieve the best accuracy.
Interface segment parameter Various algorithms can be used to determine the segment parameter α.All algorithms must produce the same result in this step when the interface normal n is known.We, thus, used the simplest analytical algorithms (Scardovelli and Zaleski, 2000), which does not need iterative algorithm.
First, we standardized the calculated interface normal n as n 1 > 0, n 2 > 0, n 3 > 0 and n 1 +n 2 +n 3 = 1.When the interface normals were negative values, we transformed the coordinates to the normal nonnegative.Second, we calculated the inverse problem α = f (c, n) by utilising the fact that the fluid fraction c is a function of the interface normal n and the segment parameter α as follows: We determined α by solving the cubic equations analytically, where α max is the summation of the standardized interface normal as: and F is the Heaviside function as:

Advection of the interface cell
In VOF method, the following advection equation has to be solved by an proper discretisation scheme: In this study, the Lagrangian-explicit method (Aulisa et al., 2003(Aulisa et al., , 2007) was adopted to discretise Eq. 22. Fig. 4 also shows a schematic illustration of the method.In this method, the surface, which is approximated by line ab, moves to line cd with the face velocity at the next time step.The blue area in Fig. 4 is the outgoing mass flux toward the right-hand cell in the next time step.As a result, the fluid fraction c t+1 can be determined by the incoming and outgoing mass flux as follows: where V L is the mass flux toward the left-hand cell, VR is the mass flux toward the right-hand cell, and VC is the remaining mass flux.

VALIDATION A classic dam breaking benchmark
We simulated some dam breaking benchmarks.Fig. 5 shows the initial condition of the dam breaking flow (Martin and Moyce, 1952).In this benchmark, we validated the convergence behaviour of the model accuracy depending on the competing grid size, which is an important factor in determining the applicability of the model to other problems from an engineering perspective.The model accuracy was evaluated based on the dimensionless position of the surge front in the time evolution.
Tab. 1 shows the calculation parameters.In this benchmark, the dimensionless time T and distance from the left side of the competing domain Z were defined using the realistic t (s), z (m), and a (m) shown in Fig. 5 and gravity g m/s 2 as follows (Martin and Moyce, 1952): The effect of gate opening was considered to shift the plots (Nichols and Hirt, 1971).Fig. 6 shows the dimensionless time series position of the surge front.Another experimental data from (Koshizuka and Oka, 1996) were also plotted.The results of our model were in good agreement with the experimental data from the Case2 and Case3 simulations.Our model has a convergence of the calculation capability that depends on the resolution through the simulations.We found that our model achieves results that match the experimental data well, and we can resolve at least 80 computing cells per the characteristic length a in the benchmark.
On the other hand, the underestimations in the Case1 simulation are due to the difference in physical properties caused by the coarse grid size.The Reynolds number in the flow field will be smaller in the LBM if we cannot ensure a sufficient spacing resolution.To overcome this problem, it can be considered that another parameterisation approach for the relaxation time τ, determined by not the viscosity ν but the Reynolds number in the flow.This method might fix the Reynolds number as turbulence.We, however, believe that a practical solution for the problem is to ensure high resolution in space because we could not determine the Reynolds number in unsteady flows, such a free surface flow.The solid wall boundary condition should also be treated carefully.Nishi and Doan (2013) have pointed out that three is room improvement for the modified bounce-back scheme especially in free surface model because the accuracy of the bounce-back scheme in the LBM depends on the relaxation time τ.

A dam breaking flow with a single obstacle
From the above section, the simple dam breaking benchmark was successfully simulated by our proposed free surface model.Here, this benchmark is extended with the addition of a solid obstacle in the centre of the flow domain by Kölke (2005) (Fig. 7).
Upon removal of gate, the collapsing of the water column is initiated and the front evolves and finally impacts the solid obstacle.At the moment, the flow pattern changes drastically, similar to what occurs during wave impact on structures.As the water height decreases, highly curved jet evolves, and the water is deflected to the top of the obstacle.
Fig. 8 compares between the calculated results with the experimental data by Kölke (2005).We fixed ∆x = 1.25 × 10 −3 and ∆t = 3.12 × 10 −6 to reduce the weak-compressibility. Very good agreement can be seen, at least for this qualitative comparison.The jet evolves slightly faster in computations tan in the experiment.Hence, we conclude that a proper grid resolution is essential, to capture at least the large-scale sprays and the resulting energy dissipation, and to accurately represent the solid body in conjunction with the above discussions.We also validated that, in LBM, a careful treatment of the corner boundaries is crucial to achieve good results.For example modified slip boundary conditions in LBM, an accurate normal vector is required in the corner cells, which is discontinuous points.If a smoothed geometry (Shigeto and Sakai, 2013;Sun and Sakai, 2015) is used, where the corner is rounded out and the corner node uses a linear interpolation of the neighbouring surface normal vectors, the jet shape is affected.For this reason, the normal vector at the obstacle corners is set to the normal vector of each respectively participating face.We should use a high-order scheme (Succi et al., 1989;Maier et al., 1996) for the boundary condition for a more precise simulation.In addition, some bubbles remained in the fluid domain and did not rise.Our model cannot address the bubble dynamics directly due to single-phase VOF approach.The remaining bubbles, however, negatively affected the stability of the model and the calculation accuracy of the interface normal n.Numerical treatments are required to address the bubbles.

Breaking wave in a rectangular tank
Finally, we simulated a three-dimensional breaking wave at an air-water interface.Numerical simulations based on fully non-linear potential flow theory (Longuet-Higgins and Cokelet, 1976;Grilli et al., 1989) showed that a periodic sinusoidal wave of large amplitude, with initial velocities specified from linear wave theory, is not stable in a fully non-linear model and rapidly overturns and breaks.
The initial velocity field and interface shape for a two-dimensional linear wave of height H in depth h, specified in a three-dimensional, are given by Lubin et al. (2006) with z denoting the vertical direction and z = 0 at the undisturbed air-water interface: where σ = 2π/T is the wave angular frequency, k = 2π/L is the wave number, and h * = h + , where (y) = 0.08y is a small linear perturbation of the seafloor geometry.In this thesis, we set L = 0.1 (m) and T = 0.308 (s).As a result, the initial velocity profiles is shown in Fig. 9. Fig. 10 compares the calculated interface shapes between our results and the reference data (Lubin et al., 2006).We fixed the spacing resolution as (x, y, z) = (1024, 128, 256) in this simulation.
As expected, the wave quickly overturns and develops a plunging jet.The agreement of our results with Lubin et al. (2006) results is good, although the shape after breaking wave was a little weakened.This is caused by the averaging extrapolation approach for empty cells, which become interface cells at the next time step.The extrapolation velocities for the empty cell crucial affects a numerical accuracy in the VOF approach.However, using high order schemes easily causes a numerical instability.We, thus, used averaging extrapolation approach in this study.

CONCLUSION
In this study, we presented the implementation the single-phase free surface model based on PLIC approach to LBM to enable three-dimensional tsunami modellings more efficient.Moreover, we performed the validation of the proposed model in classic dam breaking benchmark and breaking waves problems.The following conclusions were drawn from the discussions above: • We have developed a fully explicit three-dimensional free-surface model by LBM.Our model reduces interface oscillations and diffusions by PLIC approach.
• We found that the proposed model can overcome the defects of the conventional free surface model in LBM, and can simulate the flow field naturally.
• We determined the appropriate spacing resolution by the characteristic length of the dam breaking benchmarks.Moreover, we demonstrated the application to other problems.
LBM, however, might require finer resolution than that of other computational fluid dynamic.The adaptive mesh refinement (Zabelok et al., 2015;Onodera and Idomura, 2018) technique will be needed to execute more efficient simulations.In addition to model development, we will investigate wave forces acting on structures (Kleefsman et al., 2005) by our model, which allows the weak-compressibility to appear from the viewpoint of coastal engineering.Table 1: Calculation parameters of the dam breaking flow (Martin and Moyce, 1952).

Figure 2 :
Figure 2: The missing distribution functions f i on the interface cells, (red functions) functions from empty gas cells, (green functions) functions along the interface normal n determined by VOF algorithms (broken line).

Figure 3 :Figure 4 :
Figure 3: A cut area is the region inside the parallelepiped ABCDEFGH and below the plane IJK.

Figure 5 :
Figure 5: The calculation domain of the Martin and Moyce (1952) dam breaking flows.

Figure 6 :
Figure 6: Timeseries of the dimensionless position of the surge front.