Multi-scale semi-Lagrangian lattice Boltzmann method
MMulti-scale semi-Lagrangian lattice Boltzmann method
N. G. Kallikounis, B. Dorschner, and I. V. Karlin ∗ Department of Mechanical and Process Engineering, ETH Zurich, 8092 Zurich, Switzerland (Dated: February 26, 2021)We present a multi-scale lattice Boltzmann scheme, which adaptively refines particles’ velocityspace. Different velocity sets, i.e., higher- and lower-order lattices, are consistently and efficientlycoupled, allowing us to use the higher-order lattice only when and where needed. This includesregions of either high Mach number or high Knudsen number. The coupling procedure of differentlattices consists of either projection of the moments of the higher-order lattice onto the lower-orderlattice or lifting of the lower-order lattice to the higher-order velocity space. Both lifting andprojection are local operations, which enable a flexible adaptive velocity set. The proposed schemecan be formulated both in a static and an optimal, co-moving reference frame, in the spirit of therecently introduced Particles on Demand method. The multi-scale scheme is first validated througha convected athermal vortex and also studied in a jet flow setup. The performance of the proposedscheme is further investigated through the shock structure problem and a high Knudsen Couetteflow, typical examples of highly non-equilibrium flows in which the order of the velocity set plays adecisive role. The results demonstrate that the proposed multi-scale scheme can operate accurately,with flexibility in terms of the underlying models and with reduced computational requirements.
I. INTRODUCTION
With its roots in kinetic theory, the lattice Boltzmannmethod (LBM) describes the evolution of fluid flow viathe propagation and collision of discretized particle dis-tribution functions (populations) f i ( x , t ), which are asso-ciated with a set of discrete velocities and constructed torecover the macroscopic Navier-Stokes equations (NSE)in the hydrodynamic limit. LBM has matured to a com-petitive alternative to conventional numerical solvers,with a vast range of applications including compress-ible flows [1], complex moving geometries [2], multiphaseflows [3, 4] and rarefied gas dynamics [5], to mention afew.While LBM has indeed conquered a large range of fluiddynamics, most popular LBM models use so-called stan-dard lattices such as the D2Q9 or the D3Q27 in two orthree dimensions ( D = 2 ,
3) with Q=9 and Q=27 dis-crete velocities, respectively. While this is mainly dueto their simplicity and efficiency, the limited number ofspeeds puts severe restrictions on their range of validity.On the other hand, a systematic increase of the numberof velocities to so-called high-order or multispeed latticeshas been shown to extend the range of validity signifi-cantly. High-order lattices can be constructed systemat-ically either by discretizing the Boltzmann equation onthe roots of the Hermite polynomials [5, 6] or by entropyconsiderations, yielding a set of so-called admissible lat-tices [7, 8]. Note that the roots of the Hermite polyno-mials are irrational numbers and thus require off-latticepropagation schemes such as the semi-Lagrangian LBM[9, 10], whereas admissible lattices as in [7, 8] remainon-lattice with integer-valued velocities. The increase ofaccuracy of such lattices can be exploited in many appli- ∗ Corresponding author; [email protected] cations such high-speed flows, non-equilibrium gas flowsor relativistic fluids [1, 11–13] to name a few.In this paper, we will restrict our attention to com-pressible as well as non-equilibrium flows but the pro-posed concepts are generic and can be used for all multi-scale applications using high-order lattices. In particu-lar, it is well known that the lack of Galilean invarianceand insufficient isotropy of standard lattices, limits clas-sical LB models to isothermal, low-Mach number flows[14, 15] and the extension of LBM to high-speed com-pressible flows is still an active area of investigation. Forinstance, so-called augmented LB models have been de-veloped in [16, 17] to mitigate these shortcomings by in-troducing non-local corrections into the kinetic equationsto eliminate the error terms in the momentum and en-ergy equation, which arise due to the constraints of thestandard lattices. Promising results have been shownin recent contributions, featuring both variable Prandtlnumber and adiabatic exponent [17]. Moreover, the so-called Particles on Demand (PonD) method has recentlybeen proposed in [18, 19], which eliminates the Galileaninvariance errors of the standard lattices from the outsetby representing the populations in a co-moving referenceframe. Note that while both of these approaches get bywith minimal velocity sets, another alternative to lift theaforementioned constraints is the use of multispeed lat-tices. There, the increase of the number of speeds mod-erates the lattice constraints and pertinent moments torecover the full Navier-Stokes-Fourier (NSF) equations inthe hydrodynamic limit can be represented by the lattice[20–23]. It must be noted however that while multispeedlattices can extend the range of velocity significantly, theassociated temperature range typically decreases with anincrease of the number of particle’s velocities. Recently,an entropic LBM realization of multispeed lattices hasdemonstrated promising results for both trans- and su-personic flows[1].High-order lattices can also be used to increase ac- a r X i v : . [ phy s i c s . f l u - dyn ] F e b curacy in non-equilibrium flows. The degree of non-equilibrium or rarefaction is usually quantified by theKnudsen number, which is defined as the ratio of themolecular mean free path and a characteristic length. Ithas been shown both analytically and numerically [24, 25]that by increasing the number of discrete velocities (orderof Gauss-Hermite quadrature), LB models can capturenon-equilibrium effects of wall-bounded flows beyond theNSF level.With the examples from above in mind, it needs tobe mentioned that while high-order lattices can providea more accurate description of the flow, they come athigh computational costs, which can make these mod-els prohibitive for flows with realistic complexity in threedimensions. Fortunately, for most practical applications,the regions requiring high-order velocity sets are typicallyconfined to a small sub-region of the entire computationaldomain. Hence, significant computational resources canbe saved by using a multi-scale description, which useshigher-order lattices only when and where needed. Inthat spirit, a variety of different multi-scale frameworks,coupling different methods have been proposed in the lit-erature [26]. For example, a multi-scale model, couplingLBM for continuum regions to direct simulation MonteCarlo (DSMC) for the high Knudsen number regions, wasproposed in [27, 28] for steady-state simulations. Fur-thermore, so-called discrete velocity models (DVM) havebeen shown to be successful in simulating rarefied gases[29–31] and multi-scale schemes with adaptively refinedphase space meshes have been proposed to reduce bothcomputational time and memory [32, 33]. DVMs havealso been coupled with the more efficient LBM, whichwas used in continuum flow zones whereas the DVM wasrestricted to the rarefied regions only [34, 35]. Finally, inthe realm of LBM, a finite difference LB scheme for rar-efied gas dynamics was proposed in [36], where differentlattices are coupled in a static manner by a non-local ex-trapolation procedure. While interesting, this approachis limited to static phase space refinement and suffersfrom severe stability issues due to a non-local ad hoc cou-pling procedure of different lattices based on extrapola-tion procedures.In this work, we propose a multi-scale LBM scheme,which alleviates these issues and allows for adaptivephase space refinement using a consistent and local cou-pling procedure. For maximum lattice flexibility, we use asemi-Lagrangian advection procedure to naturally decou-ple the velocity space from the physical space [9, 18, 37–40]. Higher-order lattices are thus only used when andwhere necessary, while preserving second-order accuracy.We further shed light on the nature of multi-scale prob-lems and the range of validity of coupling procedures.The proposed scheme is then validated on examples inboth high-Mach and high-Knudsen number flows butthe scheme can be beneficial whenever large lattices areneeded in a confined region. For illustration of the cou-pling procedure, we use a dual population, multispeedLBM model with variable Prandtl number and adiabatic exponent as proposed in [1] for high-Mach regions. Fur-ther, while high-Mach number flows are intrinsically cap-tured in PonD through an adaptive reference frame, weextend PonD’s range of validity to non-equilibrium flowsusing multispeed lattices.The paper is organised as follows. In Sec. II the two-population model is presented and it’s semi-Lagrangianrealisation is explained. Subsequently, in Sec. III themulti-scale coupling scheme is introduced. Numerical re-sults are presented in Sec. IV, with simulations of anathermal convected vortex, a jet flow, the shock struc-ture problem and a high Knudsen Couette flow. Finally,concluding remarks are given in Sec. V. II. MODEL DESCRIPTIONA. Discrete velocities
Without a loss of generality, we consider discretespeeds in two dimensions, D = 2, formed by tensor prod-ucts of roots of Hermite polynomials c iα , c i = ( c ix , c iy ) . (1)The roots of the lowest three Hermite polynomials arecollected in Tab. I for the sake of completeness. Followingthe standard nomenclature, we refer to the correspondingdiscrete speeds (1) as D Q D Q
16 and D Q
25 models.Each model is characterized by the lattice temperature T L and the weights W i associated with the vectors (1), W i = w ix w iy , (2)where w iα are weights of the Gauss–Hermite quadrature,see Tab. I. Note that the lattice temperature is matchedat the outset, T L = 1, for all the models under consider-ation. TABLE I. Lattice temperature T L , roots of Hermite poly-nomials c iα and weights w iα of the D = 1 Gauss–Hermitequadrature, and nomenclature.Model T L c iα w iα D = 2 D Q , / D Q ±√ / D Q ± (cid:112) − √ √ / D Q ± (cid:112) √ − √ / D Q / D Q ± (cid:112) − √
10 (7 + 2 √ / ± (cid:112) √
10 (7 − √ / With the discrete speeds (1), the particles’ velocities v i are defined relative to a reference frame, specified by theframe velocity u ref and the reference temperature T ref , v i = (cid:114) RT ref T L c i + u ref , (3)where R is the gas constant. Two reference frames of in-terest shall be considered below. The local co-movingreference frame is specified by the local temperature T = T ( x , t ) and the local flow velocity u = u ( x , t ).The lattice reference frame is specified by u ref = and T ref = T L /R . B. Kinetic equations
For the sake of presentation, we consider a two-population kinetic model for ideal gas with a variableadiabatic exponent and Prandtl number [1], f i ( x , t ) − f i ( x − v i δt, t − δt ) = ω ( f eqi − f i ) + ( ω − ω )( f ∗ i − f eqi ) , (4) g i ( x , t ) − g i ( x − v i δt, t − δt ) = ω ( g eqi − g i ) + ( ω − ω )( g ∗ i − g eqi ) , (5)where f eqi and g eqi are local equilibria while f ∗ i and g ∗ i arequasi-equilibrium populations. Moreover, δt is the timestep and ω and ω are relaxation parameters relatedto the dynamic viscosity and thermal conductivity [1],respectively. The local conservation laws for the density ρ , momentum ρ u and the total energy ρE are, ρ = Q − (cid:88) i =0 f i = Q − (cid:88) i =0 f eqi , (6) ρ u = Q − (cid:88) i =0 v i f i = Q − (cid:88) i =0 v i f eqi , (7) ρE = Q − (cid:88) i =0 v i f i + Q − (cid:88) i =0 g i = Q − (cid:88) i =0 v i f eqi + Q − (cid:88) i =0 g eqi . (8)We consider ideal gas with the internal energy of theform U = C v T , where C v is the specific heat at constantvolume. The total energy is, ρE = C v ρT + ρu . (9)The quasi-equilibrium populations are, f ∗ i = f eqi + W i Q : ( θ / c i ⊗ c i ⊗ c i − RT θ / sym( c i ⊗ ))6( RT ) , (10) g ∗ i = g eqi + W i θ / q · c i RT , (11)where sym( . . . ) denotes symmetrization and θ is the re-duced temperature, θ = RTT L , (12) while the non-equilibrium third-order tensor Q and theheat flux vector q are, Q = Q − (cid:88) i =0 ( v i − u ) ⊗ ( v i − u ) ⊗ ( v i − u )( f i − f eqi ) , (13) q = Q − (cid:88) i =0 ( v i − u )( g i − g eqi ) . (14)When the local flow velocity u ( x , t ) and the local tem-perature T ( x , t ) are used to gauge particles’ velocities(3), u ref = u ( x , t ) , T ref = T ( x , t ) , (15)we say that kinetic equations (4) and (5) are formulatedin the co-moving reference frame, where the equilibriumpopulations depend only on the density and the temper-ature, f eqi = ρW i , (16) g eqi = (cid:18) C v − D R (cid:19) T ρW i . (17)With the co-moving reference frame and for any of themodels of sec. II A, the hydrodynamic limit of the ki-netic equations (4) and (5) are the standard equations ofcompressible gas dynamics, with the dynamic viscosity µ ,thermal conductivity κ and the bulk viscosity ξ relatedto the relaxation parameters ω and ω as follows, µ = (cid:18) ω − (cid:19) ρRT δt, (18) κ = C p (cid:18) ω − (cid:19) ρRT δt, (19) ξ = (cid:18) C v − DR (cid:19) µ. (20)Here C p is the specific heat of ideal gas at constant pres-sure, C p = C v + R . The Prandtl number is defined asPr = C p µ/κ , and the adiabatic exponent is γ = C p /C v .In the following, we set R = 1, without loss of generality. C. Semi-Lagrangian realization
1. Co-moving reference frame
The implementation of the propagation using the co-moving reference frame requires a transformation be-tween non-equal reference frames which we remind forthe sake of completeness [18]. Specification of a refer-ence frame shall be denoted λ , λ = { u , T } , (21) v λi = √ θ c i + u , (22)where θ is the reduced reference temperature (12). Ina given reference frame λ , the Q linearly independentmoments of the f - and of the g -populations are definedas, M λmn = Q − (cid:88) i =0 f λi (cid:16) √ θc ix + u x (cid:17) m (cid:16) √ θc iy + u y (cid:17) n , (23) N λmn = Q − (cid:88) i =0 g λi (cid:16) √ θc ix + u x (cid:17) m (cid:16) √ θc iy + u y (cid:17) n , (24)where m, n ∈ { , . . . , √ Q − } . Equations (23) and (24)establish a linear map of the Q -dimensional populationvectors f λ and g λ into the moment vectors M λ and N λ .Denoting M λ the Q × Q the matrix of this map, we write(23) and (24) as, M λ = M λ f λ , (25) N λ = M λ g λ . (26)Consider two reference frames, λ and λ (cid:48) . Following [18],populations are transformed based on the principle ofindependence of the moments on the reference frame, M λ = M λ (cid:48) , (27) N λ = N λ (cid:48) . (28)With (27) and (28), the populations are transformedfrom the reference frame λ to the reference frame λ (cid:48) , f λ (cid:48) = G λ (cid:48) λ f λ , (29) g λ (cid:48) = G λ (cid:48) λ g λ , (30)where the transfer matrix G λ (cid:48) λ reads, G λ (cid:48) λ = M − λ (cid:48) M λ . (31)Finally, populations are reconstructed at a point x andtime t using Lagrange interpolation,¯ f λi ( x , t ) = p (cid:88) s =1 a s ( x − x s ) G λλ s f λ s ( x s , t ) , (32)¯ g λi ( x , t ) = p (cid:88) s =1 a s ( x − x s ) G λλ s g λ s ( x s , t ) , (33) where the summation is carried out over the collocationpoints x s and a s are interpolation functions. Below, weuse the third-order Lagrange polynomials. Reconstruc-tion (32) and (33) takes into account the fact that thereference frames λ s at various collocation points x s maydiffer from one another. Thus, the corresponding popula-tions f λ s and g λ s are transformed into a target referenceframe λ using (29) and (30) prior to the interpolation.Evaluation of the populations at the monitoring point x at time t involves the propagation and the collisionsteps. In the propagation step, semi-Lagrangian advec-tion is performed using the reconstruction (32) and (33)at the departure points of characteristic lines, f i ( x , t ) = ¯ f i ( x − v i δt, t − δt ) , (34) g i ( x , t ) = ¯ g i ( x − v i δt, t − δt ) . (35)Here, particle’s velocities v i are defined relative to theco-moving local reference frame at x and t . In order tofind the co-moving reference frame, a predictor-correctorprocess is executed as follows: In the prediction step, thelocal reference frame λ = { u , T } , is initialized usingthe local flow velocity and the local temperature availablefrom the previous time step, u = u ( x , t − δt ) , T = T ( x , t − δt ). The density, momentum and temperatureare consequently computed, ρ = Q − (cid:88) i =0 f λ i , (36) ρ u = Q − (cid:88) i =0 v λ i f λ i , (37) ρ E = Q − (cid:88) i =0 ( v λ i ) f λ i + Q − (cid:88) i =0 g λ i . (38)The computed velocity (36) and temperature (37) definethe corrector reference frame λ = { u , T } at the moni-tor point and the propagation step (34) is repeated withthe updated reference frame. The predictor-correctorprocess is iterated until convergence with the limit values, ρ ( x , t ) , u ( x , t ) , T ( x , t ) , f λ ( x ,t ) i = lim n →∞ ρ n , u n , T n , f λ n i , defining the density, velocity, temperature and the pre-collision populations at the monitoring point x at time t . The predictor-corrector iteration loop ensures that thepropagation and the collision steps are performed at theco-moving reference frame, in which the local equilibriumpopulations (16, 17) are exact.
2. Lattice reference frame
If instead of an adaptive co-moving reference PonDframe, the fixed lattice reference frame λ L = { , T L } isused, we arrive at the semi-Lagrangian LBM [9, 10, 37].In this special case, the lattice equilibrium populations f Li , g Li are evaluated by transforming the exact equilib-rium populations (16) and (17) to the lattice referenceframe, f L = G { ,T L }{ u ,T } f eq , (39) g L = (cid:18) C v − D (cid:19) T f L . (40)Similarly, the quasi-equilibria are transformed to the lat-tice reference frame, f ∗ ,L = G { ,T L }{ u ,T } f ∗ , (41) g ∗ ,L = G { ,T L }{ u ,T } g ∗ . (42)A discussion is in order here. It is well known that thestandard D Q D Q Q λααα = (cid:88) i =0 (cid:0) v λiα (cid:1) f λi , (43) R λαααα = (cid:88) i =0 (cid:0) v λiα (cid:1) f λi . (44)While the off-diagonal third-order elements, Q xyy = M and Q yxx = M , as well as the off-diagonal fourth-orderelement R xy = M are the same in both, the co-movingand the lattice reference frame since they are among themoments transformed, the equilibrium values of the non-transformed moments are not equal in both references.Indeed, using the co-moving equilibrium (16) we find, Q eqααα = 3 ρT u α + ρu α , (45) R eqαααα = 3 ρT + 6 ρT u α + ρu α . (46)With the diagonal moments (45) and (46), the equilib-rium third-order moment tensor and the once-contractedequilibrium fourth-order moment tensor become, Q eqαβγ = ρT ( u α δ βγ + u β δ αγ + u γ δ αβ ) + ρu α u β u γ , (47) R eqαβ = ρT (( D + 2) T + u ) δ αβ + ρ (( D + 4) T + u ) u α u β . (48)The equilibrium Maxwell–Boltzmann relations (47) and(48) are precisely what is required to recover the com-pressible flow equations in the thermodynamic limit. On the contrary, the lattice equilibrium (41) returns, insteadof (45) and (46), Q Lααα = 3 ρT L u α , (49) R Lαααα = 3 ρ ( T L ) + 3 ρT L u α . (50)Finally, if the temperature is fixed to the lattice tem-perature T L , the athermal model is recovered with f L = G { ,T L }{ u ,T L } f eq . (51) III. MULTI-SCALE COUPLING SCHEMEA. The overlap and switching
Each discrete velocity set can be characterized by itsdomain of validity. Going by an example, the standard D Q u (cid:46) . D Q
25 of Tab.I provides a larger domain of validity [5]. It is thereforeexpected that the use of the D Q
25 velocity set instead ofthe D Q V q = { v qi , i = 0 , . . . , q − } , V Q = { v Qi , i = 0 , . . . , Q − } , where q < Q .At a monitoring point x , at time t , we can choose be-tween the two corresponding population vectors, f q and f Q . We further assume that at ( x , t ) the validity do-mains of the higher-order Q -model and the lower-order q -model overlap and both models are equally applicable.Switching between the models then requires one of thetwo operations: • Lifting: The lifting operation switches from thelower-order q -model to the higher-order Q -modeland is required to increase the accuracy if the avail-able state f q is close to the limit of its validity do-main and we need to proceed with the higher-order . . . . . . M a − − − − ǫ FIG. 1. Deviation of the xxx component of the D Q Q Lxxx (49) from it’s Maxwell–Boltzmann counterpart Q eq xxx (47) as a function of Mach num-ber Ma = u/c s . (cid:15) = ( Q Lxxx − Q eq xxx ) /Q eq xxx . The deviationvanishes for the case of the D Q
25 lattice. model. The lifting operation thus requires us tospecify a map, f q → f Q . (52) • Projection: The projection operation switches fromthe higher-order Q -model to the lower-order q -model and is required to maintain the accuracywith lesser velocities if the available state f Q iswithin the validity domain and we can to proceedwith the lower-order model. The projection opera-tion thus requires us to specify a map, f Q → f q . (53)Note that the lifting and projection must occur whenboth low- and high-order are equally valid and a propertransfer of information between the different models ispossible, i.e., the validity domains must overlap. This is ageneral feature that the models of any multi-scale methodmust respect, such that a coupling strategy can be phys-ically devised. We shall now specify the two switchingoperations separately. B. Lifting
We use notation m q for the q -dimensional vector ofmoments of the low-order model, available at ( x , t ). Withthe q × q matrix M q , we have m q = M q f q . (54)On the other hand, a moment vector of the Q -populations M Q can be considered as an element of a direct sum, M Q = M q ⊕ M Q − q , (55) where M q ∈ Im( M q ) is a vector from the image of thematrix M q while M Q − q is the rest of the higher-ordermoments. The lifting operation consists in specifying theindividual contributions as, M q = m q , (56) M Q − q = M eq Q − q . (57)This construction can be understood as a decomposi-tion between conserved, ”slow” and ”fast” components.The information of the low-order population vector arefully retained, both equilibrium and non-equilibrium con-tributions, and the moments M q are identified as the”slow” components. On the other hand, the miss-ing M Q − q are considered as ”fast” moments, which arestrongly enslaved by the dynamics of the slow moments.The final moment vector M q → Q can then be written as M q → Q = m q ⊕ M eq Q − q . (58)With the moments M q → Q completed, the populations f Q are found by moment inversion, f Q = M Q − M q → Q , (59)where M Q is the Q × Q moments-to-populations trans-form for the higher-order model. C. Projection
In the projection step, a high-order population vector f Q is mapped to a lower-order f q . Fortunately, the high-order lattice can represent Q linearly independent mo-ments, which contains the subset of the first q moments,which are required to construct the populations f q of thelow-order lattice. Hence, in contrast to the lifting pro-cedure, there is no missing information and all linearlyindependent moments m Q → q are operationally availablefrom f Q , m Q → q = M q . (60)The population vector f q is obtained by moment in-version, f q = M q − m Q → q . (61)The concepts of the two mappings are generic and ap-ply with either f (density, momentum conserving lattice)or g (energy conserving lattice) populations. Finally, westress that both lifting and projection are fully local oper-ations, which leads to high efficiency, numerical stabilityand flexibility for an adaptive velocity refinement. D. Semi-Lagrangian propagation
Now we discuss the semi-Lagrangian propagation inthe multi-scale setting, using the lifting and the projec-tion operators. We consider the population f i ( x , t ), witha discrete velocity v i , belonging to either the lower- orhigher-order velocity set. The semi-Lagrangian propa-gation starts with the calculation of the departure point x dep ,i = x − v i δt , and the identification of the surround-ing collocation points, (see, Eq. (32)). Here we must takeinto consideration that the collocation points use in gen-eral velocity sets of different order. Thus, the lifting andprojection operators are applied, such that the popula-tions involved in the reconstruction are all expressed inthe same velocity set. Figure 2 shows an example of thesemi-Lagrangian propagation in the multi-scale setting.The implementation can be summarised in the follow-ing steps:1. Calculation of departure point, x dep ,i = x − v i δt and identification of the collocation points.2. For all collocation points: • If velocity set order is lower than velocity setorder of ( x , t ): apply lifting. • If velocity set order is higher than velocity setorder of ( x , t ): apply projection.3. Calculation of f i ( x dep ,i , t − δt ), using the recon-struction formula Eq. (32).The first and third step are the same as for the caseof a single velocity set throughout the domain, whereasthe second step is active only near the interface of dif-ferent velocity sets. The collision step that follows thepropagation proceeds as in the case of a uniform velocityset. The modification of the semi-Lagrangian algorithmthat enables the multi-scale feature is independent of theunderlying model and can be used with single or dou-ble populations and with static or co-moving referenceframe (PonD). It must be noted that in a co-moving ref-erence frame, the semi-Lagrangian propagation step isperformed iteratively in a predictor-correct loop. Dur-ing the iterations the lifting/projection operations needto be performed once to each point that is close to theinterface region (step 2 above). IV. NUMERICAL RESULTS AND DISCUSSION
In this section we investigate the performance andaccuracy of the multi-scale scheme. The velocity setsthat are used in the simulations, generated by Gauss-Hermite quadrature, are described in section II A. Thetime step δt of the simulations that follow is such thatCFL = (max | c i | δt ) /δx = 1 .
0, where c i correspond to thelattice velocities. FIG. 2. Schematic of the multi-scale semi-Lagrangian propa-gation. For simplicity we use four collocating points for theinterpolation (inside dotted rectangle) and we consider thecase where the lifting operation is active. The departure point x dep ,i = x − v i δt is surrounded by the collocating points a to d . The black arrows of points a and b represent the low-order population vector and with the lifting operation thecorresponding population in the higher-order velocity set isapproximated (red arrows). The semi-Lagrangian step con-cludes with the interpolation of the populations at the collo-cating points. A. Multi-scale flows
Before going into to the details of our validation, it isinstructive to remind the nature of multi-scale problems.Multi-scale flows are characterized by large variations ofcharacteristic quantities which most commonly includesbut is not limited to spatial scales, time scales, Machnumber or Knudsen number. In what follows, we onlyconsider the Mach or Knudsen number, since in thosecases we will benefit from using high-order lattices. Forlarge spatio-temporal scales on the other hand, we referto existing grid-refinement techniques such as [41]. How-ever, when devising a numerical scheme, which bridgessuch large variations by coupling different methods fordifferent scales, it is important to realize that consis-tent coupling is only possible if there is a region, wherethe validity range of both schemes overlap. In our case,this corresponds to regions of the flow, where both lat-tices can provide an accurate description of the flow field.These regions will eventually become the interface regionbetween both lattices and thus correspond to the phasespace refinement criterion. Hence, if such overlap regionsdo not exist in the flow, the use of a multi-scale schemeis inappropriate.We will first test these ideas for variation in Mach us-ing an athermal convected vortex and an athermal jetflow using a fixed reference frame ( λ = { , T L } ) as exam-ples. The underlying model for these simulation consistof a single population, without quasi-equilibrium relax-ation. The motivation of using different velocity sets inthe case of a fixed reference frame stems from the errorsof the D2Q9 lattice due to non-Galilean invariance asthe Mach number grows and the refinement criterion be-tween high and low order velocity sets is based on a Machnumber threshold. We proceed with the shock structureproblem using a co-moving reference frame (PonD) andthe full compressible model (f,g populations with vari-able Prandtl number and adiabatic exponent), in whichthe velocity set is dictated by the variation of the Knud-sen number.
1. Athermal vortex convection
The reference case of an athermal, convected vortex isstudied with the proposed multi-scale scheme, using thestandard D2Q9 lattice and the D2Q25. The initial condi-tions of the velocity and density fields are the following, u x = U − (cid:15) (cid:18) y − y c R c (cid:19) exp (cid:20) − ( x − x c ) + ( y − y c ) R c (cid:21) , (62) u y = (cid:15) (cid:18) x − x c R c (cid:19) exp (cid:20) − ( x − x c ) + ( y − y c ) R c (cid:21) , (63) ρ = ρ exp (cid:20) − M a v (cid:18) − r R c (cid:19)(cid:21) , (64)where U is the advection velocity, (cid:15) is the strength ofthe vortex and R c is the characteristic radius [42]. Theadvection and vortex Mach numbers are set to 0.25 and0.6 respectively.For the simulations a [200 × Q eq αβγ . Figure 1 shows the deviation of the Q eq xxx compo-nent from it’s Maxwell-Boltzmann counterpart Q eq,MB xxx =3 ρc s + ρu as a function of Mach number, where c s is thespeed of sound. Here a threshold value M a thr = 0 . FIG. 3. Density (top) and velocity set (bottom) at twodifferent times (∆ t = 500 iterations). The high order D2Q25follows the vortex as it is advected.FIG. 4. Density contours of an athermal convected vortexsimulation with D2Q9, D2Q25 and coupled D2Q9/25.
2. Athermal jet flow
The next case to be analysed through the multi-scalescheme is an athermal jet flow. The inflow velocity profileconsists of a base low-speed flow with Mach number equalto
M a L = 0 .
15 and a high-speed jet region superimposedsymmetrically in the centre of the domain with diameter D jet and Mach number M a H = 0 .
6. At the inlet a zeropressure gradient and fixed velocity are imposed, whereasa non-reflecting boundary condition is prescribed at theoutlet. Free-stream boundary conditions are applied atthe top and bottom planes. The simulation is carried outon a [750 × D jet = 25. Theviscosity is adjusted such that the jet Reynolds numberis Re = D jet U jet /ν = 1000.The multi-scale model uses the D2Q9 and D2Q25 ve-locity sets, at a fixed reference frame at rest. The refine-ment criterion is based on the local Mach number andset to M a thr = 0 .
3. Figure 6 shows the instantaneous x-velocity field and the corresponding spatial distributionof the velocity set at the same time step. The time av-erage results for the velocity set and the x-velocity areshown in figure 7, where it is evident that the D2Q25 x . . . . . ρ D2Q9 D2Q25 D2Q9
D2Q25D2Q9D2Q9/25
FIG. 5. Density of an athermal convected vortex simulationwith D2Q9, D2Q25 and coupled D2Q9/25. The vertical dot-ted lines on the plot indicate the high-order lattice for thecase of the coupled scheme. is dominant across the center-line and near the jet inletwhereas the D2Q9 is used at low-speed regions of theflow. Note that the average ratio of the nodes using thehigh order velocity set is roughly 15% and the discus-sion regarding the computational speed-up is presentedin section IV C .In figure 8, the multi-scale simulation using the hy-brid D2Q9/25 lattice is further compared with simula-tions in which the D2Q9 and D2Q25 were applied uni-formly in the entire domain. The top plot shows theaverage x-velocity (normalised with the maximum veloc-ity of the jet) across the center-line, while the middleand the bottom plots show the average x-velocity acrossplanes normal to the flow, with distances from the in-let x c /D jet = 12 ,
16. It is apparent that the multi-scalesimulation is in good agreement with the simulations ob-tained with the D2Q25. In contrast the D2Q9 shows sig-nificant deviations with respect to D2Q25. This furthervalidates the proposed multi-scale scheme.
3. Shock structure
The shock structure problem is a classical problem inkinetic theory of gases, in which non-equilibrium phe-nomena dominate the flow. It is well known that since thethermodynamic variables vary on a scale of few mean freepaths, traditional continuum equation such as Navier-Stokes-Fourier equations fail to correctly describe theshock wave structure.The numerical setup consists a quasi one-dimensionalplane shock wave, with an initial step at the center of thecomputational domain. The underlying model consistsof the two-population realisation of the PonD method,with Prandtl number fixed to
P r = 2 / γ = 5 /
3. The upsteamand downstream flow values are connected through the
FIG. 6. Simulation of athermal jet flow with coupledD2Q9/D2Q25. (Top) instantaneous snapshot of the x-velocity field; (Bottom) velocity set throughout the domainat the same time.
Rankine-Hugoniot conditions [43]. The upstream meanfree path for hard sphere molecules is defined as, λ = 165 √ πγ (cid:20) µ α p (cid:21) , (65)where p , α , µ are the pressure, speed of sound and theviscosity of the gas upstream of the shock respectively.The viscosity varies with the temperature as, µ = µ (cid:18) TT (cid:19) s , (66)where for the case of hard spheres s = 0 . FIG. 7. Simulation of athermal jet flow with coupledD2Q9/D2Q25. (Top) Time average of the velocity set; (Bot-tom) Time average of the x-velocity. ture, normal stress and heat flux are defined as follows, ρ n = ρ − ρ ρ − ρ , T n = T − T T − T , ˆ σ xx = σ xx p , ˆ q x = q x p √ T , (67)where the subscripts 1 and 2 indicate the upstream anddownstream values w.r.t. the shock wave.The D2Q16 and D2Q25 velocity sets are coupled forthe computation of the shock structure. The high ordervelocity is used inside the region of the shock wave andthe low order in the rest of the computational domain.The refinement criterion is based on the local Knudsennumber, which can be computed as follows, Kn = λL , L = (cid:12)(cid:12)(cid:12)(cid:12) φdφ/dx (cid:12)(cid:12)(cid:12)(cid:12) , (68) x/D jet . . . . . u a v g / u m a x D2Q9D2Q25D2Q9/25 − − x/D jet . . . . . . u a v g / u m a x D2Q9D2Q25D2Q9/25 − − x/D jet . . . . . u a v g / u m a x D2Q9D2Q25D2Q9/25
FIG. 8. Comparison of athermal jet flow with D2Q9, D2Q25and coupled D2Q9/D2Q25. (Top) Mean streamwise veloc-ity profiles along the centerline; (Middle) Mean streamwisevelocity profiles at x/D jet = 12; (Bottom) Mean streamwisevelocity profiles at x/D jet = 16 where φ is density (in general other macroscopic fieldscan be used such as Mach number). The threshold valuewas set to Kn thr = 0 .
02 , a typical value suggested bythe literature [36, 44]. The simulations were carried outwith a quasi one-dimensional setup with the resolution ofthe lattice ∆ x corresponding to 0 . λ . The velocity setin the domain is initialized with the D2Q16 lattice andwherever the local Knudsen number exceeds the thresh-old value is refined to D2Q25. The time evolution of thevelocity set and the computed Knudsen number is shownin figure 9. The D2Q25 is initially concentrated in the1 − − − x n . . . . . K n t = 100 t = 500 t = 4000 − − − x n V e l o c i t y S e t t = 100 t = 500 t = 4000 FIG. 9. Evolution of Knudsen profile (top) and velocity set(bottom) at times t = 100 , , middle of the domain and gradually expands towards theentire shock transition layer.The steady-state results for M a = 1 . ρ n = 0 . . √ πλ .It can be seen that the density, temperature and veloc-ity profiles match very well with the reference data. Thenormal stress and heat flux profiles are also in good agree-ment with the reference data. B. Micro-Couette flow
The simulation of microscale flows is of both practi-cal and theoretical interest especially when the Knudsennumbers is not negligible. A classical case to test the be-haviour of numerical solvers for this regime is the sheardriven Couette flow. It has been established that thechoice of the discrete velocities set as well as the bound-ary conditions are of paramount importance for the accu-rate description of the fluid near the walls. In particularthe standard D2Q9 lattice does not suffice to capture theKnudsen layer and significant deviations occur for the ve-locity profile at finite Knudsen numbers. Furthermore it − − − x n . . . . . . ρ n , T n D2Q16 D2Q25 D2Q16 ρ n , D Q / T n , D Q / − − − x n . . . . . . u n D2Q16 D2Q25 D2Q16 u n , D Q / − − − x n − . . . σ xx , q x D2Q16 D2Q25 D2Q16 σ xx , D Q / q x , D Q / FIG. 10. (Multi-scale D2Q16/25) Density,temperature, ve-locity, normal stress and heat flux profiles for Ma = 1 . has been reported that even-orders velocity sets performsignificantly better than odd-orders [25].At this point, we must stress that the micro-Couetteflow is not a true multi-scale problem. Indeed the Knud-sen number is uniform in the entire domain and in con-trast to shock-structure problem, where the Knudsennumber varies within the shock transition layer. For aquantitative discussion, we use the exact solution of themicro-Couette flow, which was been found analytically2for D2Q9 and D2Q16 lattices in [11]. The authors con-cluded that while the D2Q9 can predict a slip-flow so-lution, it fails in the transient regime ( Kn (cid:39) . u w = ± . T w = T L . Thesimulations were performed with a quasi-one dimensional[300 ×
4] grid. We test the multi-scale scheme with theD2Q9 velocity set in the main flow and the D2Q16 nearthe wall boundaries. The underlying model consists ofthe two-population realisation of the PonD method, withPrandtl number fixed to
P r = 2 / γ = 5 /
3. Diffusive bound-ary conditions are implemented to efficiently capture thegas-wall interactions [25, 47].The results for the non dimensional velocity (nor-malised with the difference of the wall velocities) forKnudsen numbers equal to 0.1, 0.2 and 0.4 are pre-sented in figure 12, in which comparison is made withresults from linearised BGK [25]. The dotted verti-cal lines on the plots indicate the D2Q9/D2Q16 inter-face for the multi-scale simulations. The position of theD2Q9/D2Q16 interface is set according to the theoreticalerror analysis (1% above the main flow error).We observe that for all Knudsen numbers the multi-scale scheme matches well with D2Q9 solution in themain flow and also with the reference results near thewall boundaries, due to the local use of the D2Q16 lat-tice. However, a small discontinuity develops near theinterface region of the different velocity sets, which be-comes more prominent with increasing Knudsen number.This is an expected behaviour, caused by the discrepancybetween the solutions of D2Q9 and D2Q16, as shown inFig. 11.
C. Computational efficiency
In this section, we investigate the increase of computa-tional efficiency when using the multi-scale scheme. Tothat end, a computational domain of size ( N × N ) is de-composed in two rectangular zones which use the high-order D2Q25 and low-order D2Q9 velocity sets respec-tively. The ratio of the CPU times t D2Q25 /t hyb , where t D2Q25 refers to the time of a globally used D2Q25 and t hyb to the hybrid scheme, is measured as a function ofthe domain size N , keeping the ratio of the high-order . . . . . . x/L D Q / D Q % Kn = 0 . Kn = 0 . FIG. 11. Relative error between D2Q9 and D2Q16 for sheardriven Couette flow, based on analytical solution [11]. nodes to 10% of the computational domain fixed. Thesimulations were carried out on a eight-core PC (IntelCore i7-9700@3GHz) and the timings are shown in Fig-ure 13. The speedup that can be observed for the sim-ulation using the D2Q9 lattice in the entire domain is t D Q /t D Q = 2 .
9. Thus, the maximum speedup thatcan be achieved when the ratio of D2Q25 is fixed to10% is ( t D2Q25 /t hyb ) max = 2 .
43. For small domain sizesthe computational cost due to the lifting and projectiontransformations at the interface nodes limits the speedupbut for larger domains (
N >
V. CONCLUSION
In this work, we presented a novel multi-scale schemewith an adaptive velocity set, according to a refinementcriterion based on the local Mach or Knudsen number.Velocity sets of different order are coupled through thelifting and projection operators. Both operators involveonly local computations, which results in a robust andflexible adaptive velocity refinement. The multi-scalescheme can be implemented with either a static or co-moving reference frame and with different models (sin-gle/double population, quasi-equilibrium models). Thenumerical results with a variety of flows validated theaccuracy and efficiency of the proposed scheme. Thiswork focused on 2D applications, but the lifting andprojection, underlying the coupling scheme, hold alsoin three dimensions. Thorough investigation of the pro-posed scheme in three dimensions is left for future work.
ACKNOWLEDGMENTS
This work was supported by European Research Coun-cil (ERC) Advanced Grant No. 834763-PonD. Computa-3 . . . . . . x . . . . . . u D2Q9 D2Q16 Kn=0.1
Linearised BGKD2Q9D2Q9/160 . . . . . . x . . . . . . u D2Q9 D2Q16 Kn=0.2
Linearised BGKD2Q9D2Q9/160 . . . . . . x . . . . . . u D2Q9 D2Q16 Kn=0.4
Linearised BGKD2Q9D2Q9/16
FIG. 12. Shear driven Couette flow with D2Q9, D2Q16and coupled D2Q9/16. Profiles of the normalised velocityfor Knudsen numbers Kn = 0 . Kn = 0 . Kn = 0 . tional resources at the Swiss National Super ComputingCenter CSCS were provided under grant No. s897. Appendix A: Convergence order
The accuracy in space of the multi-scale scheme is stud-ied with the incompressible Taylor-Green vortex. Thehigh-order D2Q25 is used inside a circular disk centered in the domain and the D2Q9 is used outside. Figure 14shows the scaling of the L error, which indicates thatthe multi-scale scheme retains second-order accuracy. N . . . . . t D Q / t h y b D Q D Q FIG. 13. Speed up ratio t D Q /t hyb as a function of thedomain size N. Ratio of high-order nodes (D2Q25) is fixed to10%.
16 32 64 128 256 L − − L e rr o r D Q /
25 : L ∼ L − . D Q L ∼ L − . FIG. 14. Scaling of L error for D2Q9 and coupled D2Q9/25for Taylor-Green vortex. Appendix B: PonD transformation matrix
The implementation of the PonD method relies on thetransformation of a population vector in different refer-ence frames, λ = { u , T } and λ (cid:48) = { u (cid:48) , T (cid:48) } . The transfor-mation is a linear operation, based on moment matching, f λ (cid:48) = M − λ (cid:48) M λ f λ = G λ (cid:48) λ f λ . (B1)Below we provide a simple script to obtain analyticalexpression for the transformation matrix for the case ofD2Q9 velocity set, using the symbolic toolbox of Matlab, syms ux uy th; % syms ux_p uy_p th_p; %syms M M_p; %Q=9;cx=[0 1 0 -1 0 1 -1 -1 1];cy=[0 0 1 0 -1 1 1 -1 -1];pow_x=[0,1,0,1,2,0,2,1,2];pow_y=[0,0,1,1,0,2,1,2,2];for i=1:1:Qfor j=1:1:QM(i,j)=(th*cx(j)+ux)^(pow_x(i))*(th*cy(j)+uy)^(pow_y(i));M_p(i,j)=(th_p*cx(j)+ux_p)^(pow_x(i))*(th_p*cy(j)+uy_p)^(pow_y(i));endendTransfMatrix=inv(M_p)*M; The same script with can be applied for other velocitysets, with a redefinition of the lattice velocities and thelist of moments.
Appendix C: Lifting and projection operations
To illustrate the coupling scheme we consider the lift-ing and projection operations, for the case of D2Q9-D2Q16 coupling. The velocity set of the D2Q9 { v i } lattice is formed from the tensor product of the discretevelocities in 1D ( − c, , c ) and the D2Q16 lattice { v i } from ( − c , − c , c , c ).For the lifting operation, we need to con-struct the higher-order moment vector M → =[ M , M , · · · , M ], given the lower-order populationvector f . The 9 first moments are directly computedfrom f , M kl = i =8 (cid:88) i =0 ( v ix ) k ( v iy ) l f i , (C1)where k, l ∈ [0 , f eq, , M kl = i =15 (cid:88) i =0 ( v ix ) k ( v iy ) l f eq, i . (C2)Finally, the inverse of the moment matrix M − is ap-plied to M → , f → = M − M → . (C3)For the projection operation we compute the momentvector M → = [ M , M , · · · , M ],given the higher-order population vector f as, M kl = i =15 (cid:88) i =0 ( v ix ) k ( v iy ) l f i , (C4)where k, l ∈ [0 , M − is applied to M → , f → = M − M → . (C5)The analytical forms of the matrices M − , M − canbe obtained from the following Matlab script, syms M_min M_max;syms c;Q_min=9;cx_min=[0 c 0 -c 0 c -c -c c]; cy_min=[0 0 c 0 -c c c -c -c];syms c1 c2;Q_max=16;cx_max=[c1 c1 -c1 -c1 c2 c2 c1 c1 -c1 -c1 -c2 -c2 c2 c2 -c2 -c2];cy_max=[c1 -c1 c1 -c1 c1 -c1 c2 -c2 c2 -c2 c1 -c1 c2 -c2 c2 -c2];pow_x=[0,1,0,1,2,0,2,1,2,0,3,1,3,2,3,3];pow_y=[0,0,1,1,0,2,1,2,2,3,0,3,1,3,2,3];for i=1:1:Q_minfor j=1:1:Q_minM_min(i,j)=(cx_min(j))^(pow_x(i))*(cy_min(j))^(pow_y(i));endendfor i=1:1:Q_maxfor j=1:1:Q_maxM_max(i,j)=(cx_max(j))^(pow_x(i))*(cy_max(j))^(pow_y(i));endendM_min_inv=inv(M_min);M_max_inv=inv(M_max); [1] N. Frapolli, S. S. Chikatamarla, and I. V. Karlin. En-tropic lattice Boltzmann model for compressible flows. Phys. Rev. E , 92:061301, Dec 2015.[2] B. Dorschner, S. S. Chikatamarla, and I. V. Karlin. En-tropic multirelaxation-time lattice Boltzmann method formoving and deforming geometries in three dimensions.
Phys. Rev. E , 95:063306, Jun 2017.[3] A. Mazloomi M, S. S. Chikatamarla, and I. V. Karlin.Entropic lattice Boltzmann method for multiphase flows.
Phys. Rev. Lett. , 114:174502, May 2015.[4] M. W¨ohrwag, C. Semprebon, A. Mazloomi Moqaddam,I. Karlin, and H. Kusumaatmaja. Ternary free-energyentropic lattice Boltzmann model with a high densityratio.
Phys. Rev. Lett. , 120:234501, Jun 2018.[5] Xiaowen Shan, Xue-Feng Yuan, and Hudong Chen. Ki-netic theory representation of hydrodynamics: a way be-yond the Navier–Stokes equation.
Journal of Fluid Me-chanics , 550:413–441, 2006.[6] Xiaowen Shan and Xiaoyi He. Discretization of the ve-locity space in the solution of the Boltzmann equation.
Phys. Rev. Lett. , 80:65–68, Jan 1998.[7] Shyam S. Chikatamarla and Iliya V. Karlin. Entropy andGalilean invariance of lattice Boltzmann theories.
Phys.Rev. Lett. , 97:190601, Nov 2006.[8] Shyam S. Chikatamarla and Iliya V. Karlin. Lattices forthe lattice Boltzmann method.
Phys. Rev. E , 79:046701,Apr 2009. [9] Giovanni Di Ilio, Benedikt Dorschner, Gino Bella, SauroSucci, and Iliya V Karlin. Simulation of turbulent flowswith the entropic multirelaxation time lattice Boltzmannmethod on body-fitted meshes.
Journal of Fluid Mechan-ics , 849:35–56, 2018.[10] Andreas Kr¨amer, Knut K¨ullmer, Dirk Reith, WolfgangJoppich, and Holger Foysi. Semi-Lagrangian off-latticeBoltzmann method for weakly compressible flows.
Phys.Rev. E , 95:023305, Feb 2017.[11] S. Ansumali, I. V. Karlin, S. Arcidiacono, A. Abbas, andN. I. Prasianakis. Hydrodynamics beyond Navier-Stokes:Exact solution to the lattice Boltzmann hierarchy.
Phys.Rev. Lett. , 98:124502, Mar 2007.[12] Jianping Meng and Yonghao Zhang. Accuracy analysisof high-order lattice Boltzmann models for rarefied gasflows.
Journal of Computational Physics , 230(3):835 –849, 2011.[13] M. Mendoza, B. M. Boghosian, H. J. Herrmann, andS. Succi. Fast lattice Boltzmann solver for relativistichydrodynamics.
Phys. Rev. Lett. , 105:014502, Jun 2010.[14] Y. H Qian and S. A Orszag. Lattice BGK models forthe Navier-Stokes equation: Nonlinear deviation in com-pressible regimes.
Europhysics Letters (EPL) , 21(3):255–259, jan 1993.[15] Qian, Y.-H. and Zhou, Y. Complete Galilean-invariantlattice BGK models for the Navier-Stokes equation.
Eu-rophys. Lett. , 42(4):359–364, 1998. [16] Nikolaos I. Prasianakis and Iliya V. Karlin. Lattice Boltz-mann method for thermal flow simulation on standardlattices. Phys. Rev. E , 76:016702, Jul 2007.[17] Mohammad Hossein Saadat, Fabian B¨osch, and Ilya V.Karlin. Lattice Boltzmann model for compressible flowson standard lattices: Variable Prandtl number and adia-batic exponent.
Phys. Rev. E , 99:013306, Jan 2019.[18] B. Dorschner, F. B¨osch, and I. V. Karlin. Particles ondemand for kinetic theory.
Phys. Rev. Lett. , 121:130602,Sep 2018.[19] Ehsan Reyhanian, Benedikt Dorschner, and Iliya V Kar-lin. Thermokinetic lattice Boltzmann model of nonidealfluids.
Physical Review E , 102(2):020103, 2020.[20] F. J. Alexander, S. Chen, and J. D. Sterling. Lat-tice Boltzmann thermohydrodynamics.
Phys. Rev. E ,47:R2249–R2252, Apr 1993.[21] Takeshi Kataoka and Michihisa Tsutahara. LatticeBoltzmann model for the compressible Navier-Stokesequations with flexible specific-heat ratio.
Phys. Rev. E ,69:035701, Mar 2004.[22] Q. Li, Y. L. He, Y. Wang, and W. Q. Tao. Coupleddouble-distribution-function lattice Boltzmann methodfor the compressible Navier-Stokes equations.
Phys. Rev.E , 76:056705, Nov 2007.[23] Minoru Watari. Finite difference lattice Boltzmannmethod with arbitrary specific heat ratio applicable tosupersonic flow simulations.
Physica A: Statistical Me-chanics and its Applications , 382(2):502 – 522, 2007.[24] Seung Hyun Kim, Heinz Pitsch, and Iain D Boyd.Slip velocity and Knudsen layer in the lattice Boltz-mann method for microscale flows.
Physical review E ,77(2):026704, 2008.[25] Jianping Meng and Yonghao Zhang. Gauss-Hermitequadratures and accuracy of lattice Boltzmann modelsfor nonequilibrium gas flows.
Phys. Rev. E , 83:036704,Mar 2011.[26] A. N. Gorban, N. Kazantzis, Y. G. Kevrekidis, H. C.¨Ottinger, and K. Theodoropoulos.
Model Reduction andCoarse-Graining Approaches for Multiscale Phenomena .Springer-Verlag Berlin Heidelberg, 2006.[27] G. Di Staso, H.J.H. Clercx, S. Succi, and F. Toschi.DSMC–LBM mapping scheme for rarefied and non-rarefied gas flows.
Journal of Computational Science ,17:357–369, 2016.[28] G. Di Staso, H.J.H. Clercx, S. Succi, and F. Toschi.Lattice Boltzmann accelerated direct simulation MonteCarlo for dilute gas flow simulations.
Phil. Trans. R.Soc. A , 374:20160226, 2016.[29] James E. Broadwell. Study of rarefied shear flow by thediscrete velocity method.
Journal of Fluid Mechanics ,19(3):401–414, 1964.[30] Luc Mieussens. Discrete velocity model and implicitscheme for the BGK equation of rarefied gas dynamics.
Mathematical Models and Methods in Applied Sciences ,10(08):1121–1149, 2000.[31] V. A. Titarev. Efficient deterministic modelling of three- dimensional rarefied gas flows.
Communications in Com-putational Physics , 12(1):162–192, 2012.[32] C. Baranger, J. Claudel, N. H´erouard, and L. Mieussens.Locally refined discrete velocity grids for stationary rar-efied flow simulations.
Journal of Computational Physics ,257:572 – 593, 2014.[33] Robert R. Arslanbekov, Vladimir I. Kolobov, andAnna A. Frolova. Kinetic solvers with adaptive meshin phase space.
Phys. Rev. E , 88:063301, Dec 2013.[34] V.V. Aristov, O.V. Ilyin, and O.A. Rogozin. Kinetic mul-tiscale scheme based on the discrete-velocity and lattice-Boltzmann methods.
Journal of Computational Science ,40:101064, 2020.[35] O.V. Ilyin. A method for simulating the dynamics ofrarefied gas based on lattice Boltzmann equations andthe BGK equation.
Comput. Math. and Math. Phys. ,58:1817–1827, 2018.[36] Jianping Meng, Yonghao Zhang, and Xiaowen Shan.Multiscale lattice Boltzmann approach to modeling gasflows.
Phys. Rev. E , 83:046701, Apr 2011.[37] Dominik Wilde, Andreas Kr¨amer, Dirk Reith, and Hol-ger Foysi. Semi-lagrangian lattice Boltzmann method forcompressible flows.
Phys. Rev. E , 101:053306, May 2020.[38] C. Shu, X. D. Niu, and Y. T. Chew. Taylor-series expan-sion and least-squares-based lattice Boltzmann method:Two-dimensional formulation and its applications.
Phys.Rev. E , 65:036708, Mar 2002.[39] M. Cheng and K. C. Hung. Lattice Boltzmann methodon nonuniform mesh.
International Journal of Computa-tional Engineering Science , 05(02):291–302, 2004.[40] A Bardow, I. V Karlin, and A. A Gusev. Generalcharacteristic-based algorithm for off-lattice Boltzmannsimulations.
Europhysics Letters (EPL) , 75(3):434–440,aug 2006.[41] Benedikt Dorschner, Nils Frapolli, Shyam S Chikata-marla, and Ilya V Karlin. Grid refinement for en-tropic lattice Boltzmann models.
Physical Review E ,94(5):053311, 2016.[42] Gauthier Wissocq, Jean-Fran ¸cois Boussuge, and PierreSagaut. Consistent vortex initialization for the athermallattice Boltzmann method.
Phys. Rev. E , 101:043306,Apr 2020.[43] Anderson John D.
Modern compressible flow: with his-torical perspective . McGraw-Hill, 1990.[44] G.A. Bird.
Molecular Gas Dynamics and the Direct Sim-ulation of Gas Flows . The Oxford engineering scienceseries. Clarendon Press, 1994.[45] Taku Ohwada. Structure of normal shock waves: Directnumerical analysis of the Boltzmann equation for hard-sphere molecules.
Physics of Fluids A: Fluid Dynamics ,5(1):217–234, 1993.[46] C. Cercignani.
Theory and Application of the BoltzmannEquation . Scottish Academic Press, Edinburgh, 1975.[47] Santosh Ansumali and Iliya V. Karlin. Kinetic boundaryconditions in the lattice Boltzmann method.