Lattice Boltzmann Method for wave propagation in elastic solids with a regular lattice: Theoretical analysis and validation
Maxime Escande, Praveen Kumar Kolluru, Louis Marie Cléon, Pierre Sagaut
LLattice Boltzmann Method for wave propagation in elastic solids with aregular lattice: Theoretical analysis and validation
Maxime Escande a,b , Praveen Kumar Kolluru c , Louis Marie Cl´eon d , Pierre Sagaut b, ∗ a Ecole Polytechnique, Route de Saclay, 91120 Palaiseau, France. b Aix Marseille Univ, CNRS, Centrale Marseille, M2P2 UMR 7340, 13451 Marseille, France. c Engineering Mechanics Unit, Jawaharlal Nehru Centre for Advanced Scientific Research, Jakkur, Bangalore 560064, India. d Sorbonne Universit´e, CNRS, Institut Jean Le Rond d’Alembert, F-75005 Paris, France.
Abstract
The von Neumann stability analysis along with a Chapman-Enskog analysis is proposed for a single-relaxation-time lattice Boltzmann Method (LBM) for wave propagation in isotropic linear elastic solids,using a regular D2Q9 lattice. Different boundary conditions are considered: periodic, free surface, rigidinterface. An original absorbing layer model is proposed to prevent spurious wave reflection at domainboundaries. The present method is assessed considering several test cases. First, a spatial Gaussian forcemodulated in time by a Ricker wavelet is used as a source. Comparisons are made with results obtainedusing a classical Fourier spectral method. Both P and S waves are shown to be very accurately predicted.The case of Rayleigh surface waves is then addressed to check the accuracy of the method.
Keywords:
Lattice Boltzmann Method, regular lattices, elastic solids, Navier equation, Poisson ratio,surface waves
1. Introduction
Lattice Boltzman Methods (LBM) are now a popular approach in Computational Fluid Dynamics (e.g.see [1, 2]), with a wide range of applications, ranging from nearly-incompressible flows to multiphase flows,flows in porous media, combustion and recently to supersonic flows, e.g. [3, 4, 5]. They exhibit a very highefficiency when compared to many classical Naviers-Stokes-based solvers, leading to a rapidly growing usefor both academic and industrial purposes. Their main advantage comes from their explicit and compactcharacter, along with the use of Cartesian grids. The common feature of Lattice-Boltzmann-like methods isthat they rely on the following set of m advection-collision equations: ∂f i ∂t + c i · ∇ f i = − τ ( M [ f ] − f i ) , i = 0 , m − c i is a set of arbitrary solution-independent velocities, M [ f ] is a function ∗ Corresponding author
Email address: [email protected] (Pierre Sagaut)
Preprint submitted to Elsevier September 15, 2020 a r X i v : . [ c s . C E ] S e p f f designed to recover the targeted macroscopic equations (the Maxwellian in classical hydrodynamics)and τ a relaxation time. The computational efficiency of LBM-like approaches comes from the Strangsplitting used to decouple the advection and the collision step [6]: the advection step simplifies as a swapin memory without additional floating point operations, while the collision step is strictly local in eachcell. The macroscopic quantities such as density, velocity and pressure in fluid mechanics are recovered bypost-processing the functions f i . A key element in LBM is the choice of the lattice, i.e. the set of m nodes(one per couple ( f i , c i )) needed to solve Eq. (1) at a given node. The dimension N comes from a trade-offbetween the accuracy, stability and computational efficiency of the method. A lattice is generally denoted DnQm , where n and m are the space dimension and the number of nodes in the lattice, respectively.While LBM originally appeared as a way to solve numerically the Boltzmann equation that statisticallydescribes fluid flow physics, it can also be interpreted and extended as a smart way to solve sets of PDEs,especially for wave propagation [7]. In that case, they must be seen as a purely numerical way to solve anequation, starting from an ad hoc change of variables to recover a set of coupled advection-reaction equations.In order to preserve the efficiency of the method, the advection operator must be linear and the reactionterm should be written as a relaxation term. Using this approach, LBM-like methods have been proposedfor many physical models, e.g. the heat conduction equation (e.g. [8]), the Schr¨odinger equation (e.g. [9])and Maxwell’s equations (e.g. [10, 11]).The present paper addresses the issue of developing a LBM-like method for the propagation of waves inelastic solids. Up to now, extension of LBM for solid mechanics has been addressed by a very few authorsonly, and a general fully satisfactory method is still to be developed.As a matter of fact, a LBM approach has been proposed to solve the Lam´e equation, which governsthe static deformation for linear elastic solids [12]. An unsteady model, that was able to simulate wavepropagation in Poisson solids has been proposed in [13], using a finite difference propagation operator witha limiter to enhance the stability of the method. This approach is restricted to the case of Poisson solids,and necessitates the use of very small time steps to prevent numerical instability. A more general methodwas recently proposed by [14], which is not restricted to Poisson solids and it is based on the collide-streamimplementation that characterizes the fluid LBM approaches. This method is based on the use of complexcrystallographic lattices. The case of propagation of shock waves in solids was addressed in the 1D case in[15], in which a finite element method is used for the propagation step.Previous unsteady methods use what is referred to as ”multi-speed” lattices, since the lattices involvea large number of nodes resulting in wide computational stencils. In some of them, lattice have eitherhalf-speeds: velocities in between the neighbour lattice site (so a double staggered grid is needed), or double-speeds (or more), i.e. velocities above the neighbour lattice site. O’Brien et al. [13] uses a D2Q13 (D2Q9 +( ± ,
0) form vectors) in 2D or a D3Q25 (D3Q19 + ( ± , , f i . Second, ELM is based on Newton’s second law, each nodes in the lattice being tied to the node underconsideration by a spring, the characteristic of the later being tuned to recover the desired effect. In LBM,one solves coupled advection-relaxation equations, leading to different numerical features.Therefore, defining a LBM method based on a regular lattice is a first step towards simulation of complexsolid media. It is proposed here to adapt Murthy’s method to regular lattices (see Section 2), and to analyzeits features by carrying a Linear Stability Analysis (see Section 4) and deriving rigorously the associatedmacroscopic equations via a multiscale analysis that mimics the Chapman-Enskog expansion used in thefluid case (see Section 3). The method is then validated considering propagation of S and P waves in aninfinite medium (see Section 5).Another extension of the model proposed in [14] is then performed by addressing finite computationaldomains and boundary conditions (see Section 6) , along with the definition of an original sponge layertechnique to avoid spurious wave reflection. These developments are validated considering the propagationof Rayleigh surface waves. 3 . Lattice Boltzmann Method for wave propagation in elastic solids The LB method presented here is based on the one proposed by Murthy et al. [14] which employs higher-order Crystallographic Lattices [16, 17]. These lattices were shown to result in a gain in terms of compactnessof the computational stencil and number of unknowns per node. The proposed LB model for elastic solidsshares some features with the well-known basic LBM for fluids, more precisely with weakly compressibleisothermal LBM with BGK collision operator. In this work, we extend this LB model for regular lattices,particularly to the widely used D2Q9 lattice instead of a crystallographic lattice.The goal here is to recover the Navier equation for homogeneous linear isotropic elastic solids, whichyields the following wave propagation equation: ρ ∂ t u α = ( λ + µ ) ∂ α ∂ β u β + µ∂ β u α + F α (2)where ρ is the density (assumed to be homogeneous), u α is the displacement, λ and µ are the Lam´ecoefficients and F α is a bulk external force. Introducing the mass flux j α = ρ∂ t u α , we have ∂ t j α = λ + µρ ∂ α ∂ β j β + µρ ∂ β j α + 1 ρ ∂ t F α (3)By separating a longitudinal and transverse parts in this equation, one can recover both P waves (com-pressive) and S waves (shear) with their respective wave speeds being v P = (cid:115) λ + 2 µρ , v S = (cid:114) µρ . (4)The velocity ratio of the primary and secondary waves is v P v S = (cid:114) − ν − ν , (5)with Poisson ratio ν defined as ν = λ λ + µ ) . (6)Equation (3) can be recast with a moment chain with a source term as follows [31] ∂ t ρ + ∂ α j α = 0 ∂ t j α + ∂ β P αβ = µ − λρ ∂ α ρ + F α ∂ t P αβ + ∂ γ Q eqαβγ = 0 (7)with Q eqαβγ = µρ ( j α δ βγ + j β δ αγ + j γ δ αβ ) (8)and P αβ behaving as the stress tensor but with the opposite sign, i.e. P αβ = − σ αβ , leading to a set ofconservation equations whose mathematical structure is close to the one used for non-Newtonian fluids, withthe important difference that the mass flux appears in place of the momentum.4he force term in the mass flux equation is necessary for non-Poisson solids (Poisson solids are solids with λ = µ , i.e. ν = 0 . ρ is still used as a constant. From now on, we merge theelastic force µ − λρ ∂ α ρ and external force F α in a generic source term S α = µ − λρ ∂ α ρ + F α (9)In practice, the spatial derivative ∂ α ρ is computed via a second-order centered finite difference scheme ∂ α ρ ( xxx ) = ρ ( xxx + ∆ x α ∆ x α ∆ x α ) − ρ ( xxx − ∆ x α ∆ x α ∆ x α )2∆ x α . (10)We define the forcing term S i as S i = w i c iα S α b (11)The macroscopic quantities are recovered computing the moments of f eqi according to the followingrelations (cid:88) i f eqi = ρ (cid:88) i f eqi c iα = j α − ∆ t S α (cid:88) i f eqi c iα c iβ = P αβ (cid:88) i f eqi c iα c iβ c iγ = Q eqαβγ = b ( j α δ βγ + j β δ αγ + j γ δ αβ ) (12)with b being the “sound speed” of the set, matching the S wave speed v S = µ/ρ in Eq. (8). So one can seethat the classical concept of sound speed in LBM for fluids c s matches with the S wave speed here. So forregular lattices, one obtains v S = b = 1 / √ (cid:88) i f eq,fluidsi = ρ (cid:88) i f eq,fluidsi c iα = ρv α (cid:88) i f eq,fluidsi c iα c iβ c iγ = ρc s ( v α δ βγ + v β δ αγ + v γ δ αβ ) (13)where v α = ∂ t u α = j α /ρ is the fluid velocity.The equilibrium distribution functions for the elastic solid LBM are given by f eqi = w i (cid:32) ρ + j α c iα b + P nαβ (cid:0) c iα c iβ − b δ αβ (cid:1) b (cid:33) (14)5ith P nαβ = P αβ − ρb δ αβ (15)This is the same form as the equilibrium distribution for fluids, replacing u α u β by P nαβ : f eq,fluidsi = w i ρ (cid:32) u α c iα b + u α u β (cid:0) c iα c iβ − b δ αβ (cid:1) b (cid:33) (16)This is consistent with the second order moment in equation (12) since for fluids (cid:88) i f eq,fluidsi c iα c iβ = ρu α u β + ρb δ αβ ≡ P αβ (17)The resulting discrete Lattice Boltzmann equation with the source terms is: f i ( xxx + c i c i c i ∆ t, t + ∆ t ) − f i ( xxx, t ) = − ∆ tτ ( f i ( xxx, t ) − f eqi ( xxx, t )) + ∆ t (cid:18) − ∆ t τ (cid:19) S i (18)A second order accurate explicit time integration can be obtained using the same change of variables asin fluid LBM. In the rest of the article we consider a non-dimensionalized system, with the simplest latticeunits ∆ t = 1, ∆ x α = 1.In the next section, we present a Chapman-Enskog analysis for this model.
3. Associated macroscopic conservation equation: multiscale Chapman-Enskog-type analysis
As in the fluid LBM case, it is assumed that the solution f i corresponds to a small perturbation aboutthe equilibrium state f eqi . To recover the macroscopic evolution equations associated to Eqs. (18), (14)and (9), it is proposed to mimic the Chapman-Enskog multiscale expansion classically used in fluid LBM tothis end. Here, this procedure is purely mathematical, since there is no link with the physical Boltzmannequation for statistical fluid physics in the case of solid LBM. Introducing the small parameter (cid:15) (which isthe Knudsen number in the dilute gas case), one considers the following expansion f i = f eqi + (cid:15)f (1) i + (cid:15) f (2) i + ... (19)As in the classical Chapman-Enskog expansion, it is assumed that both S i and S α are of order O ( (cid:15) ) S α = (cid:15)S (1) α (20)and time and spatial derivatives are expanded as ∂ t = (cid:15)∂ (1) t + (cid:15) ∂ (2) t + ...∂ α = (cid:15)∂ (1) α . (21)6rom the mass, momentum and stress conservation equations Eq. (7) we obtain the following solvabilityconditions. It should be noted that the forcing is responsible for the non-zero terms: (cid:88) i f neqi = − ∆ t (cid:88) i S (1) i = 0 , (cid:88) i f neqi c iα = − ∆ t (cid:88) i S (1) i c iα = − ∆ t S (1) α , (cid:88) i f neqi c iα c iβ = − ∆ t (cid:88) i S (1) i c iα c iβ = 0 (22)Since S (1) α ∼ O ( (cid:15) ) it only affects f (1) i , hence, the order-by-order solvability conditions are (cid:88) i f ( n ≥ i = 0 , (cid:88) i f (1) i c iα = − ∆ t S (1) α , and (cid:88) i f ( n ≥ i c iα = 0 , (cid:88) i f ( n ≥ i c iα c iβ = 0 (23)Incompressible hydrodynamic equations (for Newtonian fluids) lack the last condition as only mass andmomentum are conserved quantities.By using Eq. (19) in a Taylor expansion of the Lattice Boltzmann equation Eq. (18), and by identifyingthe terms of order O ( (cid:15) ) and O ( (cid:15) ), we have: (cid:16) ∂ (1) t + c iα ∂ (1) α (cid:17) f eqi − (cid:18) − ∆ t τ (cid:19) w i c iβ S (1) β b = − τ f (1) i ∂ (2) t f eqi + (cid:18) − ∆ t τ (cid:19) (cid:16) ∂ (1) t + c iα ∂ (1) α (cid:17) (cid:32) f (1) i + ∆ t w i c iβ S (1) β b (cid:33) = − τ f (2) i (24)Taking the zeroth through third moment at O ( (cid:15) ) from Eq. (24), one obtains ∂ (1) t ρ + ∂ (1) α j α = 0 ∂ (1) t j α + ∂ (1) β P αβ = S (1) α ∂ (1) t P αβ + ∂ (1) γ Q eqαβγ = 0 ∂ (1) t Q eqαβγ + ∂ (1) κ R eqαβγκ = − τ Q (1) αβγ + b (cid:18) − ∆ t τ (cid:19) (cid:16) S (1) α δ βγ + S (1) β δ αγ + S (1) γ δ αβ (cid:17) , (25)at O ( (cid:15) ) we have ∂ (2) t ρ = 0 ∂ (2) t j α = 0 ∂ (2) t P αβ + (cid:18) − ∆ t τ (cid:19) ∂ (1) γ Q (1) αβγ = − b ∆ t (cid:18) − ∆ t τ (cid:19) (cid:16) ∂ (1) γ S (1) γ δ αβ + ∂ (1) α S (1) β + ∂ (1) β S (1) α (cid:17) . (26)7s can be seen from the above set of equations the additional term Q (1) αβγ appearing in Eq. (26) dependson R eqαβγκ via Eq. (25) which in turn depends on the choice of the discrete velocity model. For regularlattices, D2Q9, D3Q15, D3Q19 and D3Q27 , satisfying the following isotropy conditions (cid:88) i w i = 1 , (cid:88) i w i c iα c iβ = b δ αβ , (cid:88) i w i c iα c iβ c iγ c iκ = b ∆ (4) αβγκ = b ( δ αβ δ γκ + δ αγ δ βκ + δ ακ δ βγ ) , (cid:88) i w i c iα c iβ c iγ c iκ c iθ c iψ = b (cid:16) ∆ (6) αβγκθψ − δ αβγκθψ (cid:17) = b (cid:16) δ αβ ∆ (4) γκθψ + δ αγ ∆ (4) βκθψ + δ ακ ∆ (4) βγθψ + δ αθ ∆ (4) βγκψ + δ αψ ∆ (4) βγκθ − δ αβγκθψ (cid:17) , (27)the modified chain (cid:15) (25) + (cid:15) (26) reduces to ∂ t ρ + ∂ α j α = 0 ∂ t j α + ∂ β P αβ = S α ∂ t P αβ + ∂ γ Q eqαβγ = (cid:18) τ − ∆ t (cid:19) b ∂ γ (cid:0) ∂ α P nβγ + ∂ β P nαγ + ∂ γ P nαβ − ∂ γ P nκκ δ αβγ (cid:1) (28)The details of the calculations are given in Appendix A for the sake of brevity.It should be noted that the last relation in Eq. (27) is specific to regular lattices and takes into accountthe lack of isotropy of such sets. Ideally, the correct 6th order moment should be ∆ (6) αβγκθψ , but is known tobe satisfied only by very high-order on-lattice models [32].The multiscale analysis reveals the existence of the source term in the associated macroscopic evolutionequation for the stresses. This source term is of diffusive/dissipative nature, and therefore will tend to dampand smooth the stresses P nαβ . It is proportional to b ( τ − ∆ t/ t .
4. Linear Stability Analysis
To study the stability of the scheme, the modified moment chain (28) could be used to find the modifiedNavier equation and study its modes with plane analysis. However, finding the modified Navier equationdoesn’t carry all the useful information on the LBM method.In order to supplement the modified equation analysis, the von Neumann analysis of the distributionfunctions f i is now carried out (see [33] for a discussion of the correct spectral analysis and [34] for theinterpretation within the LBM framework). We consider disturbances of the plane wave form: f i ( xxx, t ) = f i exp( − i ( ωt + k · xk · xk · x )) (29)8 Real −1.0−0.50.00.51.0 I m ag i na r y −1.0 −0.5 0.0 0.5 1.0 Real −1.0−0.50.00.51.0 I m ag i na r y −1.0 −0.5 0.0 0.5 1.0 Real −1.0−0.50.00.51.0 I m ag i na r y −1.5 −1.0 −0.5 0.0 0.5 1.0 Real −1.0−0.50.00.51.0 I m ag i na r y Figure 1: Eigenvalues and stability disks for different Poisson ratio. Top-left: 0; Top-right: 0.25; Bottom-left:0.3; Bottom-right:0.4. The stability disk is plotted in green, and the eigenvalues in blue (refer to the web version of the article for the colors).
For the vector F = ( f , f , ..., f q − ) T , a matrix equation is obtained by using the moments (12) of Eq.(18) F ( t + ∆ t ) = M · F ( t ) (30)where the coefficients m ij of M are the following m lj = exp( − ik γ c lγ ∆ t ) (cid:20)(cid:18) − ∆ tτ (cid:19) δ lj + w l ∆ tτ (cid:18) c lα c jα b + i ∆ t Λ c lα sin( k α ∆ x α )2∆ x α + 12 b (cid:0) c lα c lβ − b δ αβ (cid:1) (cid:0) c jα c jβ − b δ αβ (cid:1)(cid:19) + (cid:18) − ∆ t τ (cid:19) ∆ t Λ b c lα sin( k α ∆ x α )∆ x α (cid:21) (31)with Λ = µ − λρ b = 1 − ν − ν (32)using µ/ρ = b and Eq. (6). In Eq. (31), the matrix index i is replaced by l to avoid any confusion withthe imaginary unit i = − M is then diagonalized to obtain the complex eigenvalues. These eigenvalues are plotted inFig. 1 for different Poisson ratio, with wave vectors of norm in [0 , π ], directed by the unit vector (1 , T / √ ° directed wave vector, for a Poisson ratio different than 0 .
25, there are someeigenvalues outside the stability disk, and particularly when the ν gets closer to the incompressible limit 0 . π/2 π k x k y k x k y k x k y k x k y Figure 2: Unstable wave vectors for different Poisson ratio. Top-left: ν = 0; Top- right: ν = 0 .
3; Bottom-left: ν = 0 . ν = 0 .
4. The case ν = 0 .
25 is not shown since all wave vectors are stable. Only one quadrant is displayed dueto the elementary symmetries of the set.
To investigate the other direction causing instability, we plotted on the ( k x , k y ) plane, the points where someeigenvalues were outside the stability disk (i.e. with a modulus more than one). The results are shown inFigure 2.We see that for low Poisson ratio, the instability occurs only near the 45 ° direction and for high wavevectors only (near the π -radius circle). This results in 45 ° -direction instability pattern (see figure 3), aftera long time (since only high wave-numbers are amplified, and since a Gaussian shape has very few highwave-numbers).For high Poisson ratio, the domain of instability increases until all directions are unstable at high enoughwave-numbers. This happens at ν = 0 .
4. Since we observe that instability occurs at lower wave-numbersfor ν > .
25 than for ν < .
25, the high Poisson ratio simulations diverge more quickly. But is is useful10 igure 3: Example of instability patterns. Left: 45 ° ; Right: all- directions. The color-map is saturated to blue-red colors dueto great positive-negative values displayed: a sign of instability. In the 45 ° instability pattern, the diagonal stripes are spacesby 2 . ° patterns can be observed in low ν simulation after a long time, all-directions patterns can be observedin high ν simulation (here ν = 0 . noticing that only very high wave vectors are unstable, typically | k ∆ x | > π/ ν = 0 .
4. Thereforestable simulations can be recovered by refining the grid, or by filtering high frequencies at the end of a timestep.For a Poisson ratio ν > .
25 the velocity of P waves is greater than the lattice speed: v P ( ν > . > √ v S = 1 = ∆ x ∆ t (33)It is a classic source of instability to have an information supposed to propagate faster than the lattice speed.This is another way of understanding the instability of Poisson ratio greater than 0 . ∂ α ρ in (10), the informationcan in reality propagate at speed 2 instead of 1. So, the absolute limitation of ν beyond which v P wouldfaster than 2 is ν lim = 5 / ≈ . ν = 0 . b = v S .Nevertheless, the present method can be used for practical simulations. It is unconditionnaly stable forPoisson solids, i.e. ν = 0 .
25, and can be used for finite-time simulations for ν (cid:54) = 0 .
25 at classical LBM CFLnumber. Since the instability arises only at high wave numbers.
5. Validation: Bulk wave propagation
As a first validation, one addresses the propagation of both P and S waves in the bulk of a 2D linearisotropic elastic medium. To this end, a x -force Gaussian source is placed in the middle of a domain periodicin x and y directions. The source is spatially Gaussian (instead of point source for example) to avoid high11patial frequencies (high wave numbers) which behave less precisely and can lead to instability. The typicalradius of the source is 4 ∆ x . This source is modulated through time with a classic Ricker wavelet with aperiod of 20 ∆ t .The D2Q9 simulation is compared with reference results obtained using a Fourier spectral method (withCrank-Nicolson time integration to allow for CFL=1, i.e. ∆ x = ∆ t to correspond with the LBM iterations)for several reasons: •
2D analytical solutions of the wave equation are not trivial (they make use of Bessel’s functions), • The source is not point source (which would involve convolutions for an analytical solution), • The simulation is periodic in both x and y directions.The detail of this spectral method is explained in Appendix B.The simulations are performed on 64x64, 128x128, 256x256 and 512x512 grids using τ = 0 . t alongwith ν = 0 . , . , . , .
25 and 0 .
3. Instantaneous results obtained at time 70∆ t on the 128x128 grid for ν = 0 . x = 2 N x / y = 2 N y / L norm of the error committed on the longitudinal mass flux j x is defined as p = (cid:118)(cid:117)(cid:117)(cid:117)(cid:117)(cid:116) (cid:80) N x ,N y i,j =1 , (cid:16) j LBM x ( i, j ) − j spectral x ( i, j ) (cid:17) (cid:80) N x ,N y i,j =1 , (cid:16) j spectral x ( i, j ) (cid:17) (34)Values computed are given in Table 1 and plotted in Fig. 6. It is observed that the error exhibits a weaksensitivity to the Poisson ratio, which is a good property, and that a first order of convergence is obtainedon the mass flux.To check the velocity of the bulk waves, some measurement of the travel time of P waves (time to reachthe left-right borders from the center) and S waves (time to reach the top-bottom borders from the center)have been performed for several values of the Poisson ratio(see figure 7). The maximum relative error is1.1% on the 128x128 grid, corresponding to a very accurate capture of the wave physics.The present results show that the D2Q9 Lattice Boltzmann method yields very accurate prediction ofbulk P and S waves, with an accuracy similar to those reported with other methods, e.g. [27, 25, 13, 35].12
50 100x050100 y −1.0−0.50.00.51.0 M a ss f l u x A m p li t ude y −1.0−0.50.00.51.0 M a ss f l u x y M a ss f l u x Figure 4: Comparison of j x for LBM D2Q9 with the Fourier spectral method, for Poisson ratio ν = 0 .
1. Top-left: LBM D2Q9;Bottom-left: Fourier spectral; Top-right: Ricker wavelet modulation of the source; Bottom-right: absolute difference of j x between the two models. t −0.02−0.010.000.010.02 j x SpectralLBM D2Q9
Figure 5: Comparison of the seismogram of the longitudinal mass flux j x for LBM D2Q9 with the spectral method, for a Poissonratio ν = 0. able 1: L norm of the error computed on the longitudinal mass flux j x for several combinations of the grid resolution Nx × N y and the Poisson ratio ν . ν = 0 ν = 0 . ν = 0 . ν = 0 . ν = 0 . N x = N y = 64 0.2508 0.2351 0.2227 0.2185 0.2141 N x = N y = 128 0.1112 0.1111 0.1137 0.1155 0.1164 N x = N y = 256 0.0633 0.0644 0.0667 0.0678 0.0681 N x = N y = 512 0.0349 0.0355 0.0366 0.0371 0.0371 Figure 6: Grid convergence of the L norm of the error on the longitudinal mass flux j x . ν v P / v S Figure 7: Comparison of the theoretical (dashed line) and LBM D2Q9 results (diamond markers) relation between ν and v P /v S . Boundary conditions The developments given in [14] and in the above section were focused on the bulk wave propagation ininfinite media, and only periodic boundary conditions were used. We now extend the method by derivingboundary conditions for more realistic cases: rigid interfaces with null displacement, free surface and non-reflecting outlet boundary condition.
The simplest BC is probably the ”Bounce-Back” boundary condition (BB). For fluids it represents the no-slip condition , typically used to implement a steady solid wall in a viscous fluid, leading to (for instance fora bottom wall): v x,wall = 0 with v the velocity. For solids, a rigid interface interface would mean j x,wall = 0,i.e. the displacement and the mass flux are null on the boundary. This would exist with the interface of oursolid of interest with a rigid wall (for instance another solid with a much greater density or stiffness). Theusual BB condition imposes the wall location to lie between two nodes : y wall = y b + ∆ y (cid:126)n (35)with y b the last layer near the interface and (cid:126)n the normal interface unit vector pointing outwards the solids.Therefore, the present Bounce Back condition is implemented as f ¯ i ( x b , t + ∆ t ) = f i ( x b , t ) (36)with c ¯ i being the opposite of c i : c ¯ i = − c i We consider now an interface with a medium with a much lower density and stiffness than the solid, e.g.gas or void. In the following, the case of the interface between the solid and air is used for the sake of clarity.The solid/air interface location is ∆ x/ x b located in the solid medium: x air = x b + ∆ x (cid:126)n (37)with (cid:126)n the normal interface unit vector pointing outwards the solids. The bouncing condition is now : f ¯ i ( x b , t + ∆ x ) = − f i ( x b , t ) + 2 w i (cid:18) ρ air + 12 b P nαβ,air (cid:0) c iα c iβ − b δ αβ (cid:1)(cid:19) (38)where ρ air is the density estimated at the border with a linear interpolation ρ air = 32 ρ b − ρ in (39)15
20 40 60 80 100 120 x y −1.00−0.75−0.50−0.250.000.250.500.751.00 M a ss f l u x x y −1.00−0.75−0.50−0.250.000.250.500.751.00 M a ss f l u x Figure 8: Test of the BCs: reflection of a vertical P wave on different interfaces. Left: bottom rigid wall (BB); Right: top freesurface (ABB). The star marker indicates the position of the source where ρ in is the density one node further inside the solid x in = x b − ∆ x(cid:126)n (40)The P nαβ,air tensor at the interface is : P nαβ,air = P αβ,air − b ρ air δ αβ (41)To explain how to find P αβ,air , we consider the case of a vertical interface: P αβ,air = P xx,air
00 0 (42)because P αβ,air n β = F α = 0, F α is the force applied by air on the solid at the interface, which is assumedto be negligible in usual pressure and temperature conditions. P xx,air is obtained with linear interpolationsimilar to ρ air . We now validate the two boundary conditions discussed above by considering wave reflection on thesetwo kinds of boundaries.To this end, a domain with periodic lateral conditions is defined, the BCs of interest being implementedon the top (free surface case) or bottom boundary (rigid wall case), as illustrated in Fig. 9. A directionalsource pointing vertically is added near the interface of interest. We observe that the waves are bouncing offthese interfaces. The larger circle arc is the first P-wave with a direct path and the second concentric circleis the P-wave reflected. 16 igure 9: Configuration of the first BC test
Two criteria indicate that the boundary conditions are accurately implemented.The first one deals with the sign of the velocity. For a free surface, there is no change in the sign of thespeed (i.e. phase change) after the reflection, since the condition is a Neumann condition for velocity. Forthe rigid interface, there is a change in the sign of the speed (i.e. phase inversion) because it is a Dirichletcondition on velocity. This can be understood considering the following simplified mass flux equation: ∂ t j α + ∂ β P αβ = 0 (43)Considering an horizontal boundary at y = 0, it leads to ∂ t j = − ∂ y P , yielding P ( c s t ∓ y ) = ± c s j ( c s t ∓ y ),so that P = 0 at the border is equivalent to P (0 , t ) = P i ( c s t − P r ( c s t +0) = 0, where the i and r subscriptsare related to incident and reflected waves, respectively. For rigid walls with j = 0, the Dirichlet conditiongives the result in a straightforward way.The second criterion is the speed magnitude at the border. For a free surface the incident and reflectedwave interfere constructively while for the rigid interface, the interference is destructive. This is visible witha darker (respectively lighter) area near the border. To simulate a semi-infinite medium, we need to introduce non-reflecting boundary conditions in order todefine outlet boundary conditions. Different approaches have been tested for fluids, based either on the use ofimproved BCs on the outlet plane or the implementation of sponge layers in which fluctuations are damped,or a combination of these two approaches. Numerical experiments performed in fluid LBM framework haveshown that the use of a sponge layer is mandatory to get very accurate results for internal flows or on smallcomputational domain.Contrary to fluids where the viscosity can be changed with τ to introduce a progressive damping effect,the Chapman-Enskog-like analysis given above shows that τ doesn’t induce a viscous behaviour in the presentcase. Therefore, in the present solid LBM case, it is chosen to design absorbing layers. The principle is to17moothly increase a damping factor near the borders to completely remove the reflected waves. This is doneby adding a linear penalty term in the mass flux equation: ∂ t ρ + ∂ α j α = 0 ∂ t j α + ∂ β P αβ = µ − λρ ∂ α ρ − Aj α ∂ t P αβ + ∂ γ Q eqαβγ = 0 (44)where A is the damping factor and the term − Aj α is treated as a forcing (exactly like an external forcing orthe elastic artificial forcing).The associated modified Navier equation is : ∂ t j α + A∂ t j α = λ + µρ ∂ α ∂ β j β + µρ ∂ β j α + 1 ρ ∂ t F α (45)A plane-wave analysis shows that all modes are damped in the absorbing layer.
7. Validation: Rayleigh surface waves
The last validation case is related to the simulation of surface waves, that are very important in seismology.The motion being restricted in the plane in the present paper (2D simulations), the only surface wavespossible are the Rayleigh waves (R waves). They stem from the combination of both P and S wavesand their reflections at the surface. Their relative speed with P and S waves are also varying with thePoisson ratio but they are slightly slower than S waves. In the particular case of Poisson solids ( ν = 0 . v R = v S √ . ≈ . v S [36]. Rayleigh waves have an amplitude that decays exponentially with depth.The motion of the particle is ellipsoidal, like the the motion of particle of fluids in the swell. Rayleigh waveshave been reported to be accurately captured using ELM (see [28, 26] and references given herein) but, tothe knowledge of the authors, they have not been computed via LBM up to now.To isolate and observe surface waves, the following numerical experiment has be conducted on a 100x300computational grid. A free surface condition is implemented on the top boundary and the three remainingborders are treated with absorbing layers with a thickness of 30 nodes each.As a first step, absorbing layers are used to damp bulk S and P waves in order to isolate Rayleigh surfacewaves. In the second step, absorbing layers are removed and periodic conditions are imposed on lateralboundaries in order to let the Rayleigh waves propagate freely. The source is located on the upper leftcorner so that all waves going to the left are quickly removed (absorbed by the left sponge layer). Thefastest wave (P waves) going to the right are also removed by the right absorbing layer, then (after 400 timesteps) the left-right absorbing layers are replaced by a simple left-right periodic condition.18
50 100 150 200 250 x y A Figure 10: Absorbing layers
A 2D circular bulk wave with constant energy has a typical amplitude decaying as 1 / √ r while a surfacewave in the present case will be a 1D wave, whose amplitude will be constant in the energy-preserving case.Thanks to that difference, surface waves should last longer than the S waves that are travelling to the rightjust slightly faster).The results are displayed in Fig. 11. The simulation qualitatively captures the behaviour of Rayleighwaves. At time steps 600 and 1200 there remains some S waves also going to the right because they are onlyslightly faster than Rayleigh waves (at ν = 0 . v R ≈ v S ). The speed of R waves have been measuredwith a 0.14% relative error, by measuring the number of time steps needed to achieve for two complete turnsaround the periodic surface: v measured R = 0 . , v theo R = (cid:114) . . x = 83. In fact this time was chosen, despite the weak amplitude, because the remainingS waves and other artefacts are clearly separated. The results are displayed in the Figure 12.The exponential model is simply j y = A exp( − y/d ) (47)With the exponential regression (fitting with the amplitude A and the typical depth d ), we see that thedamping is approximately of the exponential form, yet, the amplitude near the surface has a bump. In fact,19
50 100 150 200 250 300 x y −0.10.00.10 50 100 150 200 250 300 x y −0.10.00.10 50 100 150 200 250 300 x y −0.10.00.10 50 100 150 200 250 300 x y −0.10.00.10 50 100 150 200 250 300 x y −0.10.00.1 Figure 11: Vertical mass flux profile at different time, for ν = 0 .
25. The star marker indicates the position of the source. Fromtop to bottom: iteration 30, 100, 600, 1200, 1800
20 40 60 80 100
Depth A m p li t ude LBM D2Q9ExponentialAnalytical
Figure 12: Profile of the amplitude of the mass flux in a vertical section x = 83 at the iteration 1800. The regression were donewith the least square method. For the exponential model the amplitude was A = 0 . d = 21 . − . For the analytical model, the amplitude is 0.012, the wavelength is λ R = 37 . − a rigorous study of the analytical solution for the Rayleigh equation gives a more complex solution [37]: j y = − α S (cid:20) exp( − kα S y ) −
21 + α P exp( − kα P y ) (cid:21) (48)with α S ≡ (cid:112) − v R /v S and α P ≡ (cid:112) − v R /v P .This analytical model was integrated and fitted with the amplitude and the wavelength λ R of the wave(giving k = 2 π/λ R ). It can be observed that the simulation results show a better match with these analyticalresults and the regression even recovers the correct wavelength λ R (within 1.3% error) that can be alsomeasured between two maxima of the mass flux.These results clearly show that the present solid LBM is able to accurately capture surface waves, andthat the proposed boundary conditions for rigid surfaces, free surfaces and sponge layers are efficient.
8. Conclusions
The Lattice Boltzmann method proposed in [14] has been extended in several ways, by i) using a regularlattice in place of the more complex Crystallographic Lattice used in the original version, ii) defining boundaryconditions for free surface and rigid surfaces and iii) introducing efficient sponge layer technique to preventspurious wave reflection. These new elements have been supplemented by two theoretical analyses, i.e. theChapman-Enskog expansion to recover the associated macroscopic equations and the von-Neumann-typelinearized stability analysis. 21t is shown that the method the first-order accurate considering the L norm of the error on the mass flux,and that the method is unconditionnally stable for Poisson solids and stable for low-wave-number solutionfor other values of the Poisson ratio ν .Numerical experiments have demonstrated the accuracy of the present method for both bulk and surfacewaves and, to the knowledge of the authors, it is the first time that surface waves are successfully computedusing a Lattice Boltzmann method for elastic solids. Acknowledgments
The authors warmly acknowledged Dr. G. S. O’Brien for sharing his C codes and useful discussions.
Appendix A. Chapman-Enskog Analysis
The moment R eqαβγκ for regular lattices is R eqαβγκ = ρb ∆ (4) αβγκ + 12 b P nθψ (cid:32)(cid:88) i w i c iα c iβ c iγ c iκ c iθ c iψ − b δ θψ (cid:88) i w i c iα c iβ c iγ c iκ (cid:33) = ρb ∆ (4) αβγκ + b P nθψ (cid:16) ∆ (6) αβγκθψ − δ αβγκθψ − δ θψ ∆ (4) αβγκ (cid:17) = ρb ∆ (4) αβγκ + b (cid:16) P nαψ ∆ (4) βγκψ + P nβψ ∆ (4) αγκψ + P nγψ ∆ (4) αβκψ + P nκψ ∆ (4) αβγψ − P nψψ δ αβγκψ (cid:17) = ρb ∆ (4) αβγκ + b (cid:0) P nαβ δ γκ + P nαγ δ βκ + P nακ δ βγ + P nβγ δ ακ + P nβκ δ αγ + P nλκ δ αβ − P nψψ δ αβγκψ (cid:1) (49)Taking a spatial derivative ∂ (1) κ of the above moment gives ∂ (1) κ R eqαβγκ = ρb (cid:16) δ βγ ∂ (1) α ρ + δ αγ ∂ (1) β ρ + δ αβ ∂ (1) γ ρ (cid:17) + b (cid:16) ∂ (1) α P nβγ + ∂ (1) β P nαγ + ∂ (1) γ P nαβ (cid:17) + b ∂ (1) κ (cid:0) P nακ δ βγ + P nβκ δ αγ + P nγκ δ αβ − P nψψ δ αβγκψ (cid:1) = b (cid:16) ∂ (1) α P βγ + ∂ (1) β P αγ + ∂ (1) γ P αβ + ∂ (1) κ (cid:0) P nακ δ βγ + P nβκ δ αγ + P nγκ δ αβ − P nκκ δ αβγκ (cid:1)(cid:17) . (50)The time derivative of Q eqαβγ can be computed using Eq. (12) and replacing the ∂ (1) t j α terms with − ∂ (1) κ P ακ + S (1) α via Eq. (25) to be ∂ (1) t Q eqαβγ = b ∂ (1) t ( j α δ βγ + j β δ αγ + j γ δ αβ )= − b ∂ (1) κ ( P ακ δ βγ + P βκ δ αγ + P ακ δ αβ )+ b (cid:16) S (1) α δ βγ + S (1) β δ αγ + S (1) γ δ αβ (cid:17) (51)22sing the above relations Q (1) αβγ can be calculated from the final equation of Eq. (25) to be Q (1) αβγ = − τ (cid:20) ∂ (1) t Q eqαβγ + ∂ (1) κ R eqαβγκ − b (cid:18) − ∆ t τ (cid:19) (cid:16) S (1) α δ βγ + S (1) β δ αγ + S (1) γ δ αβ (cid:17)(cid:21) = − τ b (cid:104) ∂ (1) α P nβγ + ∂ (1) β P nαγ + ∂ (1) γ P nαβ − ∂ (1) κ P nκκ δ αβγκ (cid:105) − b (cid:20) ∆ t (cid:16) S (1) α δ βγ + S (1) β δ αγ + S (1) γ δ αβ (cid:17)(cid:21) (52) Appendix B. Custom Fourier spectral method
A Fourier spectral method with Crank-Nicholson time integration is used to compare bulk waves. Wedefine ( ˆ j x , ˆ j y ) T as the spatial Fourier transform of the mass flux vector ( j x , j y ) T . Let the ˆ J vector beˆ J = ˆ j x ˆ j y ∂ t ˆ j x ∂ t ˆ j y . (53)To simplify the notation, we denote a = ( λ + 2 µ ) /ρ , b = µ/ρ and d = ( λ + µ ) /ρ , such that the Navierequation for j x and j y in the Fourier domain can be written as ∂ tt ˆ j x = − a k x ˆ j x − d k x k y ˆ j y − b k y ˆ j x + 1 ρ ∂ t ˆ F x ∂ tt ˆ j y = − b k x ˆ j y − d k x k y ˆ j x − a k y ˆ j y + 1 ρ ∂ t ˆ F y (54)Hence, ∂ t ˆ J = ∂ t ˆ j x ∂ t ˆ j y − a k x ˆ j x − d k x k y ˆ j y − b k y ˆ j x + ρ ∂ t ˆ F x − b k x ˆ j y − d k x k y ˆ j x − a k y ˆ j y + ρ ∂ t ˆ F y (55)This equation is of the form ∂ t ˆ J = f (cid:16) ˆ J, ˆ F (cid:17) . The Crank-Nicholson scheme is A-Stable, which allows usto take ∆ t = ∆ x , furthermore, it is a symplectic integrator, which is beneficial for equations of the form¨ X + ω X = 0 like (55).The Crank-Nicholson scheme is ˆ J n +1 = ˆ J n + ∆ t (cid:104) f (cid:16) ˆ J n (cid:17) + f (cid:16) ˆ J n +1 (cid:17)(cid:105) (56)On simplification, the above equations can be evaluated in matrix form to be M ˆ J n +1 = N ˆ J n ˆ G x ˆ G y , (57)23ence, ˆ J n +1 = M − N ˆ J n ˆ G x ˆ G y (58)with ˆ G x = 1 ρ (cid:16) ∂ t ˆ F x,n + ∂ t ˆ F x,n +1 (cid:17) ˆ G y = 1 ρ (cid:16) ∂ t ˆ F y,n + ∂ t ˆ F y,n +1 (cid:17) , (59) M = − ∆ t
00 1 0 − ∆ t A ∆ t B ∆ t B ∆ t C ∆ t (60) N = ∆ t ∆ t − A ∆ t − B ∆ t ∆ t − B ∆ t − C ∆ t ∆ t (61)and, A = a k x + b k y B = d k x k y C = a k y + b k x (62)Then the equation (58) is just a simple matrix relation that can be iterated along with the LBM. Whenneeded, the mass flux can be obtained by an inverse Fourier transform. References [1] T. Kru¨uger, H. Kusumaatmaja, A. Kuzmin, O. Shardt, G. Silva, E. Viggen, The Lattice BoltzmannMethod. Principles and Practice, Springer, 2017.[2] Z. Guo, C. Shu, The Lattice Boltzmann Method and its applications in engineering, World Scientific,2013.[3] J. Jacob, O. Malaspinas, P. Sagaut, A new hybrid recursive regularised bhatnagar–gross–krook collisionmodel for lattice boltzmann method-based large eddy simulation, Journal of Turbulence 19 (11-12)(2018) 1051–1076. 244] Y. Feng, P. Boivin, J. Jacob, P. Sagaut, Hybrid recursive regularized thermal lattice boltzmann modelfor high subsonic compressible flows, Journal of Computational Physics 394 (2019) 82–99.[5] S. Guo, Y. Feng, J. Jacob, P. Sagaut, An efficient lattice Boltzmann method for industrial aerodynamicsflows, I: basic model on D3Q19 lattice, Journal of Computational Physics 418 (2020) 109570.[6] P. Dellar, Lattice boltzmann formulation for linear viscoelastic fluids using an abstract second stress,SIAM Journal of Scientific Computing 36 (2014).[7] Y. Guangwu, A lattice boltzmann equation for waves, Journal of Computational Physics 161 (2000)61–69.[8] W. Jiaung, J. Ho, C. Kuo, Lattice boltzmann method for the heat conduction problem with phasechange, Numerical Heat Transfer, Part B 39 (2001).[9] L. Zhiong, S. Feng, P. Dong, S. Gao, Lattice boltzmann schemes for the nonlinear schr¨odinger equation,Physical Review E 74 (2006).[10] S. Hasanoge, S. Succi, S. Orszag, Lattice boltzmann method for electromagnetic wave propagation,EuroPhysics Letters 96 (2011).[11] Y. Liu, G. Yan, A lattice Boltzmann Model for maxwell’s equations, Applied Mathematical Modelling38 (2014).[12] X. Yin, G. Yan, T. Li, Direct simulations of the linear elastic displacements field based on a lattice Boltz-mann model: Direct simulations of the linear elastic displacements field based on a lattice Boltzmannmodel, Int. J. Numer. Meth. Engng 107 (3) (2016) 234–251. doi:10.1002/nme.5167 .[13] G. S. O’Brien, T. Nissen-Meyer, C. J. Bean, A Lattice Boltzmann Method for Elastic Wave Propagationin a Poisson Solid, Bulletin of the Seismological Society of America 102 (3) (2012) 1224–1234. doi:10.1785/0120110191 .[14] J. S. N. J.Murthy, P. K. Kolluru, V. Kumaran, S. Ansumali, Lattice Boltzmann Method for WavePropagation in Elastic Solids, CiCP 23 (4) (2018). doi:10.4208/cicp.OA-2016-0259 .[15] S. Xiao, A lattice Boltzmann method for shock wave propagation in solids, Communications in Numer-ical Methods in Engineering 23 (2007).[16] M. Namburi, S. Krithivasan, S. Ansumali, Crystallographic Lattice Boltzmann Method, Sci Rep 6 (1)(2016) 27172. doi:10.1038/srep27172 . 2517] P. K. Kolluru, M. Atif, M. Namburi, S. Ansumali, Lattice boltzmann model for weakly compressibleflows, Physical Review E 101 (1) (2020) 013309.[18] I. Ispolatov, M. Grant, Lattice boltzmann method for viscoelastic fluids, Physical Review E 65 (2002).[19] P. Lallemand, D. d’Humieres, L. Luo, R. Rubinstein, Theory of the lattice boltzmann method: Three-dimensional model for the linear viscoelastic fluids, Physical Review E 67 (2003).[20] O. Malaspinas, N. F´etier, M. Deville, Lattice boltzmann method for the simulation of viscoelastic fluidflows, Journal of Non-Newtonian Fluid Mechanics 165 (2010).[21] T. Phillips, G. Roberts, Lattice boltzmann models for non-newtonian flows, Journal of Applied Mathe-matics 76 (2011).[22] G. N. Frantziskonis, Lattice Boltzmann method for multimode wave propagation in viscoelastic mediaand in elastic solids, Phys. Rev. E 83 (6) (2011) 066703. doi:10.1103/PhysRevE.83.066703 .[23] A. Gupta, M. Sbragaglia, A. Scagliarini, Hybrid lattice boltzmann/finite difference simulations of vis-coelastic multicomponent flows in confined gepmetries, Journal of Computational Physics 291 (2015).[24] G. O’Brien, C. Bean, A 3d discrete numerical elastic lattice method for seismic wave propagation inheterogeneous media with topography, Geophysical Research Letters 31 (2004).[25] G. O’Brien, C. Bean, An irregular lattice method for elastic wave propagation, Geophysical JournalInternational 187 (2011).[26] G. O’Brien, Elastic lattice modelling of seismic waves including a free surface, Computers and Geo-sciences 67 (2014).[27] G. O’Brien, Dispersion analysis and computational efficiency of elastic lattice methods for seismic wavepropagation, Computers and Geosciences 35 (2009).[28] R. del Valle-Garcia, F. Sanchez-Sesma, Rayleigh waves modeling using an elastic lattice model, Geo-physical Research Letters 30 (2003).[29] M. Xia, H. Zhou, Q. Li, H. Chen, Y. Wang, S. Wang, A general 3d lattice spring model for modelingelastic waves, Bulletin of the Seismological Society of America 107 (2017).[30] D. Polyzos, D. Fotiadis, Derivation of mindlin’s first and second strain gradient elastic theory via simplelattice and continuum models, International Journal of Solids and Structures 49 (2012).[31] S. Godunov, E. Romenskii, Elements of continuum mechanics and conservation laws, Springer, 2003.2632] H. Chen, X. Shan, Fundamental conditions for N-th-order accurate lattice Boltzmann models, PhysicaD: Nonlinear Phenomena 237 (14-17) (2008) 2003–2008. doi:10.1016/j.physd.2007.11.010 .[33] T. Sengupta, A. Dipankar, P. Sagaut, Error dynamics: Beyond von Neumann analysis, Journal ofComputational Physics 226 (2007).[34] G. Wissocq, P. Sagaut, J.-F. Boussuge, An extended spectral analysis of the lattice Boltzmann method:modal interactions and stability issues, Journal of Computational Physics 380 (2019) 311–333. doi:10.1016/j.jcp.2018.12.015 .[35] M. Xia, S. Wang, H. Zhou, X. Shan, H. Chen, Q. Li, Q. Zhang, Modelling viscoacoustic wave propagationwith the lattice Boltzmann method, Sci Rep 7 (1) (2017) 10169. doi:10.1038/s41598-017-10833-wdoi:10.1038/s41598-017-10833-w