We propose a new direct coupling scheme based on the overset technique to tackle moving boundary problems within the lattice Boltzmann framework. The scheme is based on the interpolation of distribution functions rather than moments, that is, macroscopic variables, and includes an additional hypothesis ensuring mass and momentum conservation at the interface nodes between fixed and moving grids. The method is assessed considering four test cases and considering both the vortical and the acoustic fields. It is shown that the direct coupling method results are in very good agreement with reference results on a configuration without any moving subdomain. Moreover, it is demonstrated that the direct coupling method provides an improvement of the accuracy of the lattice Boltzmann overset algorithm for aeroacoustics. In particular, a convected vortex test case is studied and reveals that the direct coupling approach leads to a better ability to conserve the vortex structure over time, as well as a reduction in spurious acoustic distorsions at the fixed/moving interface.

Fluid flows around rotating systems are nowadays found in numerous industrial applications, for example, fans, pumps, propellers, turbomachinery components, wind turbines, and blades. To better understand the physical phenomena involved in this type of flows, experiments may fail at providing enough data for a precise analysis of the induced complex flow field because of technological limitations of the measurements techniques. Consequently, the use of computational methods to tackle the modeling of flows around moving and/or deformable geometries has increasingly grown into a popular research topic.

In the past few decades, the ever-increasing advancements in computational technology have enabled the regular use of computational fluid dynamics for numerous types of flows. Most commonly, Navier–Stokes solvers using classical finite differences, finite volume, or finite element methods are used. However, in recent years, the lattice Boltzmann (LB) method has also emerged as a powerful alternative to solve the weakly compressible Navier–Stokes equations. Rather than solving the macroscopic conservation equations, the LB method is a statistical approach that describes a given fluid flow through the projection of the Boltzmann equation over a Hermite polynomial family, which provides a discrete equation based on particle distribution functions. The LB method is easily parallelizable and can deal with complex geometries and, on that account, is now widely used to model and simulate a large variety of phenomena [multiphase flows, combustion, turbulence, compressible, and thermal flows, see, e.g., Shan and Chen (1993), Filippova and Hänel (2000), Succi (2001), Sagaut (2010), Shan et al. (2006), Feng et al. (2019), Ma et al. (2020), Guo et al. (2020), Zhao et al. (2020), and Zhan et al. (2021)].

However, while many problems are now solved using classical fixed grids, the case of fluid flows involving moving or deformable bodies brings out additional complexity. Such flows, as opposed to the ones involving stationary bodies only, indeed naturally generate significant computational costs when fixed grids are used, due to the constant re-meshing needed every time the body displaces in space.

Some strategies were therefore proposed over the past to tackle this issue. One of them is the well-known immersed boundary method (IBM), introduced by Peskin (1972) and (2002), that consists in representing the fluid–structure interactions by the construction of a series of Lagrangian marker points that exactly follow the fluid and for each of which a judiciously chosen force is imposed. Originally developed within the Navier–Stokes framework, the IBM has been extended to LB methods in several papers, for example, Feng and Michaelides (2004), Zhu et al. (2011), and Favier et al. (2014). Another popular method is the Arbitrary Lagrangian–Eulerian (ALE) approach, which combines the Eulerian approach (far from the body), and the Lagrangian viewpoint (in the vicinity of the body). Thus, the equations are written on a mesh that moves in an arbitrary fashion in the vicinity of the body. The ALE approach has also been extended in several works to LB methods, see, for example, Meldi et al. (2013) for the simulation of immersed moving solids in low-speed incompressible flows, or more recently Saadat and Karlin (2020) for LB compressible flow simulations on unstructured moving meshes. In addition to this non-exhaustive list of methods, one can mention the so-called overset technique, on which the present paper is based. The overset method consists in decomposing the domain into one fixed grid and one moving grid, the rotating body being tied to the moving one. Both grids share an overlapping region and communicate through their common interfaces thanks to some interpolations. The overset method has recently been used in several LB studies of flows around moving bodies, among which Li (2011), Zhang et al. (2011), Far et al. (2020), and Lallemand and Luo (2020).

In particular, in the present work, we propose two different ways of communicating the information at the interface between moving and fixed grids, within the same computational framework. The first one consists in a classical interpolation of macroscopic values, which further leads to the reconstruction of the distribution functions at the interface. This method is referred to in the present paper as the “original overset” method. The second approach is based on the interpolation of distribution functions rather than moments, as opposed to the previous overset methods reported in the literature. Additionally, hypotheses ensuring mass and momentum conservation at the interface nodes between fixed and moving grids are drawn. In the present paper, this second approach is referred to as the “direct coupling (DC) overset” method. Even though the present work does not entirely address this issue in the sense that it still makes use of an overset framework, it is worth noting that one of the final aims of the direct coupling approach is to completely eliminate the existent overlapping area, thereby helping minimize the effort for parallelization as duplication of points would then be prevented. In this regard, the present work can be seen as a first step toward this objective by providing a tighter link between moving and fixed grids, as will be further detailed over Secs. II and III.

The direct coupling method is inspired by the work of Astoul et al. (2020), which deals with the modeling of abrupt resolution transitions of non-uniform grids and was shown to outstandingly improve the accuracy of the grid interface for aeroacoustic applications. The objective of the present work is to propose an extension of this method to the moving subdomain boundary problem and to assess its capabilities both for aerodynamic and aeroacoustic applications, in comparison to the original overset method developed within the same framework. In particular, one of the main motivations is to improve the ability of lattice Boltzmann methods (LBM) to capture acoustic properties in the context of moving boundaries, given the above-mentioned outstanding improvements observed in Astoul et al. (2020)'s work on grid refinement. Aeroacoustic applications are of particular interest because of the existing significant gap in research on this topic when moving boundaries are involved.

More precisely, the extension to overset rotating/fixed grid transitions puts forward different physical and mathematical elements from the ones exposed in Astoul et al. (2020), for example,

  • There is no re-scaling of quantities since the mesh resolution is constant.

  • The transition from fixed to moving and moving to fixed grids must properly take into account both the geometrical rotation of the grid and the non-inertial/inertial grid transition.

In this context, the present work describes a new algorithm capable of handling this type of transition, thus allowing us to simulate fluid flows involving rotating subdomains with a direct coupling approach.

This paper is organized as follows. In Sec. II, we briefly describe the concepts of the lattice Boltzmann method. The hybrid recursive regularized (HRR) collision model and the overset scheme as well as their coupling in the present methodology are presented. Section III discusses in detail the theoretical and numerical aspects of the direct coupling overset method developed in our work. Finally, Sec. IV presents a numerical validation of the direct coupling overset method with three aerodynamic test cases (an uniform flow, a Poiseuille flow, and an uniform flow past a rotating cylinder) and two aeroacoustic test cases (an acoustic pulse and a convected vortex).

Only the basic concepts of the lattice Boltzmann method will be presented here. For further details, the reader is referred to Krüger et al. (2017), Succi (2001), and Aidun and Clausen (2010).

The lattice Boltzmann equation can be used to describe the evolution in time and space of the discrete particle distribution function fi(x,t), which stands for the probability to observe a particle at position x and at time t with discrete velocity ci. With no external force, the lattice Boltzmann equation reads

fi(x+ci,t+1)fi(x,t)=Ωi(x,t)Δt,
(1)

where Ωi(x,t) is the collision operator and Δt is the time step, which is taken equal to one hereafter.

The collision model that has been adopted in this work is the hybrid recursive regularized (HRR) model (Latt and Chopard, 2006; Malaspinas, 2015; Coreixas, 2018; Jacob et al., 2018), for its overall increased stability and accuracy in comparison with the classical Bhatnagar–Gross–Krook (BGK) collision model (Bhatnagar et al., 1954). In particular, the HRR model has been validated across our investigations in different configurations involving moving subdomains and a clear enhancement in stability has been observed.

The LB method scheme for the HRR model is given by

fi=fi(0)+fi(1),
(2)

where fi(0) and fi(1) are respectively the equilibrium and off-equilibrium distribution functions, which are expressed as

fi(0)=g(0)(ρ,u)=wi(ρ+ci·(ρu)cs2+12cs4Hi(2)·a0(2)+12cs6(Hi,xxy(3)+Hi,yyz(3))(a0,xxy(3)+a0,yyz(3))+12cs6(Hi,xzz(3)+Hi,xyy(3))(a0,xzz(3)+a0,xyy(3))+12cs6(Hi,yyz(3)+Hi,xxz(3))(a0,yyz(3)+a0,xxz(3))+16cs6(Hi,xxy(3)Hi,yyz(3))(a0,xxy(3)a0,yyz(3))+16cs6(Hi,xzz(3)Hi,xyy(3))(a0,xzz(3)a0,xyy(3))+16cs6(Hi,yyz(3)Hi,xxz(3))(a0,yyz(3)a0,xxz(3)))
(3)

and

fi(1)=g(1)(ρ,u,Π(1))=wi(12cs4Hi(2)·a1(2)+12cs6(Hi,xxy(3)+Hi,yyz(3))(a1,xxy(3)+a1,yyz(3))+12cs6(Hi,xzz(3)+Hi,xyy(3))(a1,xzz(3)+a1,xyy(3))+12cs6(Hi,yyz(3)+Hi,xxz(3))(a1,yyz(3)+a1,xxz(3))+16cs6(Hi,xxy(3)Hi,yyz(3))(a1,xxy(3)a1,yyz(3))+16cs6(Hi,xzz(3)Hi,xyy(3))(a1,xzz(3)a1,xyy(3))+16cs6(Hi,yyz(3)Hi,xxz(3))(a1,yyz(3)a1,xxz(3))),
(4)

where wi and cs are, respectively, the lattice weights and the speed of sound of the lattice, the Hi(n) are the discrete Hermite polynomials of order n, and a0(n) and a1(n) are the Hermite coefficients of the equilibrium and off-equilibrium distribution, respectively,

a0(n)=i=0q1Hi(n)fi(0),
(5)
a1(n)=i=0q1Hi(n)fi(1),
(6)

where q is the number of discrete lattice velocities.

In addition, Π(1) is the deviatoric stress tensor, defined as

Π(1)=a1(2)=i=0q1Hi(n)fi(1).
(7)

Equations (3) and (4) involve macroscopic quantities, namely the density ρ and the velocity u, which are defined as the zeroth and first moments of the distribution functions, respectively,

ρ=ifi,
(8)
ρu=icifi.
(9)

In the framework of the HRR collision model, the equilibrium Hermite coefficients write (Jacob et al., 2018)

a0(n)=a0(n1)u,
(10)

with

a0(0)=ρ,
(11)

and the off-equilibrium Hermite coefficients write

a1,αβ(2)=iHi,αβ(2)(fifi(0)),
(12)
a1,αβγ(3)=uαa1,βγ(2)+uβa1,γα(2)+uγa1,αβ(2).
(13)

In order to model the motion of moving or deformable bodies within a given flow field, the method adopted in this work relies on the overset technique.

To illustrate of this technique, let us take the example of a flow past a rotating (square) cylinder (see Fig. 1).

FIG. 1.

Flow past a rotating square cylinder: representation of the interaction between fixed and moving grids.

FIG. 1.

Flow past a rotating square cylinder: representation of the interaction between fixed and moving grids.

Close modal

Two grids along with their two respective coordinate systems are used: a fixed one and a moving one. The solid body is tied to the moving grid, and the red area in Fig. 1 is an overlapping region shared by the two grids. The lattice Boltzmann equation is solved on each grid: on the fixed grid, it is done through Eq. (1), but as the moving grid is non-inertial, inertial forces have to be taken into account through an external forcing term Fi(x,t) as

fi(x+ci,t+1)fi(x,t)=Ωi(x,t)+(112τ)Fi(x,t).
(14)

The inertial forces included in the forcing term Fi(x,t) are, namely, the Coriolis and centrifugal forces. They are expressed as

F=ρ(ω×(ω×r)2ω×ucdωdt×r),
(15)

where ω is the angular velocity of the moving grid, r is the distance vector from the origin point of the moving grid to the fluid position, and uc is the velocity containing the so-called half-force correction (Krüger et al., 2017)

uc=u+F2ρ.
(16)

Thus, uc is equal to the velocity u shifted by F/(2ρ), which allows to obtain second-order accuracy (Zhang et al., 2011; Far et al., 2016; Krüger et al., 2017). This velocity is used at the collision step to calculate the equilibrium distribution functions.

More precisely, the external forcing term Fi(x,t) uses the scheme introduced by Guo et al. (2002), leading to

Fi(x,t)=wi((ciuc)·Fcs2+(ci·uc)(ci·F)cs4),
(17)

where the discrete lattice velocities ci are the ones related to the moving grid as the source term is always computed in the moving grid in its own reference frame.

In this section, we present the border condition problem in the lattice Boltzmann overset framework. As previously mentioned, there are two grids along with their two respective coordinate systems: a fixed one and a moving one. The solid body, if there is one, is tied to the moving grid. The two grids share an overlapping region, which is represented in red in Fig. 1. Information between the two grids is exchanged at the borders as shown in Fig. 1 and act as somewhat boundary conditions for the two grids. We refer to this issue as the “border condition problem.”

Let us consider the D2Q9 lattice for the sake of simplicity and zoom in at the fixed to moving border (Fig. 2). In Fig. 2, the dashed arrows represent the unknown distribution functions fi after a collision step, while the continuous arrows stand for the known distribution functions. The border condition problem can be rephrased as: How to reconstruct the unknown distribution functions?

FIG. 2.

Representation of the fixed to moving border. The black dashed arrows represent the unknown distribution functions after a streaming step, while the black continuous arrows stand for the known distribution functions. On the right part of the figure, a graphical representation of the interpolation process from the fixed to the moving grid is also shown, where the red circle stands for the interface moving node to be interpolated and the four blue squares stand for the four surrounding fixed interpolating nodes.

FIG. 2.

Representation of the fixed to moving border. The black dashed arrows represent the unknown distribution functions after a streaming step, while the black continuous arrows stand for the known distribution functions. On the right part of the figure, a graphical representation of the interpolation process from the fixed to the moving grid is also shown, where the red circle stands for the interface moving node to be interpolated and the four blue squares stand for the four surrounding fixed interpolating nodes.

Close modal

As mentioned in Sec. I, we propose two ways of addressing this issue. The first one, referred to as the original overset method, is widely used in the literature and consists in a classical interpolation of macroscopic values further leading to the reconstruction of the distribution functions at the mesh interface. On the other hand, the second approach, referred to as the direct coupling method, is based on the interpolation of distribution functions rather than moments and additionally draws complying hypotheses on mass and momentum conservation.

The detailed theoretical and numerical steps involved in both these methods are now exposed in the following.

In the following, any quantity related to the fixed grid or moving one will be referred to with a f or m superscript, respectively.

1. Fixed to moving border

Let xm be the position of a border node of the moving grid.

In the original overset-HRR algorithm, the macroscopic variables (ρ, u, Π(1)) on the fixed grid are interpolated onto the border nodes of the moving grid through a bilinear interpolation using k =8 interpolating nodes as

ρIf(xm)=ϕ(ρf(xkf)),k[1..8],
(18a)
uIf(xm)=ϕ(uf(xkf)),k[1..8],
(18b)
ΠIf,(1)(xm)=ϕ(Πf,(1)(xkf)),k[1..8],
(18c)

where ϕ is the bilinear interpolation operator and the I index means that it is an interpolated value.

The newly interpolated velocities have then to be modified to take into account (i) the geometrical rotation of the grid and (ii) the acceleration of the non-inertial moving frame of reference.

Point (i) involves a classical rotation matrix:

uI,Rf(xm)=RuIf,
(19)
ΠI,Rf,(1)(xm)=RΠIf,(1),
(20)

where the R index means that it is a geometrically rotated value.

Point (ii) implies the removal of the relative velocity of the moving grid with respect to the fixed grid vr=r×ω, where ω is the angular velocity of the moving grid and r is the distance vector from the origin point of the moving grid to the fluid position

uI,R,Pf=uI,Rfvr,
(21)

where the P index means that it is a value taking into account the transition to a non-inertial moving reference frame.

In conclusion, the macroscopic variables on the moving grid interface nodes are reconstructed as

ρm(xm)=ρIf(xm),
(22)
um(xm)=uI,R,Pf(xm),
(23)
Πm,(1)(xm)=ΠI,Rf,(1)(xm).
(24)

The distribution functions at location x=xm can be reconstructed using Eqs. (3) and (4) and the newly interpolated and transformed ρm and um as follows:

fi,regm=g(0)(ρm,um)+g(1)(ρm,um,Πm,(1)),
(25)

where g(0)(ρm,um)=fi(0) and g(1)(ρm,um,Πm,(1))=fi(1) of Eqs. (3) and (4).

Finally, at this stage, the transition to a non-inertial moving frame must also be taken into account in the distribution functions as

fim=fi,regm12Si(ρm,um)=fîf,
(26)

where Si is the source term standing for the Coriolis and centrifugal forces related to the non-inertial rotating frame, namely,

Si=wi((ciuc)·Fcs2+(ci·uc)(ci·F)cs4).
(27)

Also, the ̂ symbol means that f̂if is the final value of the distribution function using the macroscopic data coming from the fixed grid at the very end of the algorithm, meaning after all transformations have been performed (interpolation, geometrical rotation, change of reference frame).

2. Moving to fixed border

The same rationale may be used for the moving to fixed border, except that:

  • instead of removing vr=r×ω from the velocity, it is added. This is done at the very beginning of the algorithm, when extracting the eight interpolating velocities from the moving grid, since this is where vr is exactly known (i.e., in the moving region)
    k[1..8],uPm(xkm)=um(xkm)+vr(xkm),
    (28)
  • the fixed frame being the inertial reference frame, there is no source term Si and
    fif=fi,regf=g(0)(ρf,uf)+g(1)(ρf,uf,Πf,(1))=fîm.
    (29)

As mentioned above, the direct coupling method sets its starting point on the interpolation of distribution functions rather than moments, as opposed to the previous overset methods reported in the literature. This raises a number of issues, which are now detailed in the following paragraph.

1. Fixed to moving border

In order to interpolate the distribution functions from one grid to the other, the first step is to geometrically rotate the set of distribution functions through the rotation of their Hermite polynomial basis, yielding

fi;Rf=n=0Nj=0q11n!fjfH(n)(cjf)·H(n)(cim),
(30)

where the R index means that it is a geometrically rotated value and cim is obtained through the rotation matrix R

cim=Rcif.
(31)

Then, fi;Rf can be interpolated onto the moving border nodes with a bilinear interpolation method using k =8 interpolating nodes as

fi;R,Im(xm)=ϕ(fi;Rf(xkf)),k[1..8],
(32)

where ϕ is the bilinear interpolation operator and the I index means that it is an interpolated value.

The macroscopic variables can then be reconstructed as

ρm=i=0q1fi;R,Im,
(33)
um=i=0q1cifi;R,Im.
(34)

It should be noted that the deviatoric stress tensor is reconstructed by interpolation and rotation, exactly as in the original method, that is, Πm,(1)(xm)=ΠI,Rf,(1)(xm) following Eqs. (18c), (20), and (24).

The ultimate step is to take into account the transition from a fixed frame to a non-inertial moving frame. This is done by removing the relative velocity of the moving grid with respect to the fixed grid vr from the newly calculated velocity um. Then, the equilibrium and non-equilibrium functions are reconstructed as

fi;R,I,Pf,(0)(xm)=g(0)(xm;ρm,umvr),
(35)
fi;R,I,Pf,(1)(xm)=g(1)(xm;ρm,umvr,Πm,(1)),
(36)

where the P index means that it is a value taking into account the transition to a non-inertial moving reference frame.

Finally, the source term Si standing for the Coriolis and centrifugal forces related to the non-inertial rotating frame can be calculated and removed from the distribution functions, so that the transition to a non-inertial moving frame is also taken into account in the distribution functions

fim=fi;R,I,Pf,(0)+fi;R,I,Pf,(1)12Si(ρm,um)=f̂if.
(37)

Also, the ̂ symbol means that f̂if is the final value of the distribution function initially coming from the fixed grid at the very end of the algorithm, meaning after it has gone through all transformations (geometrical rotation, interpolation, change of reference frame).

2. Moving to fixed border

The same rationale may be used for the moving to fixed border, except that:

  • instead of removing vr=r×ω from the interpolating velocities ukm, it is added. This is done at the very beginning of the algorithm, prior to the geometrical rotation and interpolation of the distribution functions, since this is where vr is exactly known (i.e., in the moving region)
    fi;Pm,(0)=g(0)(xkm;ρkm,ukm+vr),
    (38)
    fi;Pm,(1)=g(1)(xkm;ρkm,ukm),
    (39)
  • the fixed frame being the inertial reference frame, there is no source term Si and
    fif=fi;P,R,If,(0)+fi;P,R,If,(1)=f̂im.
    (40)

For the sake of synthesis, in the following, the whole approach is detailed for the fixed to moving border but is easily transferrable to the moving to fixed border by interverting the f and m superscripts.

Let us first recall what was done in Astoul et al. (2020) in terms of direct coupling for grid refinement, in order to compare and identify the developments to perform in our case to elaborate our extension of the method to moving subdomains.

Firstly, as explained by Astoul et al. (2020), mesh refinement implies rescaling of physical quantities. In particular, viscosity and, consequently, relaxation time were modified in their work. This provided the following rescaling relation for the non-equilibrium part of the distribution functions on the co-located fine/coarse nodes of the interface between fine and coarse grids

fifine,(1)=Rficoarse,(1),
(41)

where fifine,(1) and ficoarse,(1) are the distribution functions at the interface related, respectively, to the fine and coarse grids, and R=0.5τfine/τcoarse is the rescaling factor.

In our case, there is no mesh refinement between fixed and moving frames; therefore, R =1. However, additional complexity is introduced by the presence of two different reference frames. Geometrical rotation and transition from a fixed frame to a non-inertial moving frame of reference (and vice versa) must therefore be addressed properly.

Second, in the work of Astoul et al. (2020), the interface shows some co-located coarse and fine nodes, which are used to reconstruct the missing populations on both coarse and fine grids. In the case of moving subdomains, we do not necessarily have co-located fixed and moving nodes as the moving grid moves freely independently from the fixed grid (see Fig. 2). Consequently, considering the case of the fixed to moving border, we assume at the interface the following relation:

fim,(1)(xm)=g(1)(ρm,um,ΠIm,(1)),
(42)

where ΠIm,(1) is the deviatoric stress tensor that has been interpolated from the fixed region (with the same bilinear interpolation method with k =8 interpolating nodes as for the previously exposed interpolation of density or velocity)

ΠIm,(1)(xm)=ϕ(Πf,(1)(xkf)),k[1..8],
(43)

where ϕ is the bilinear interpolation operator and the I index means that it is an interpolated value. In Eq. (43), Πf,(1)(xkf) is defined as

Πf,(1)(xkf)=i=0q1Hi(n)(fif(xkf)fif,(0)(xkf)).
(44)

The idea behind Astoul et al. (2020)'s work on grid refinement was to reconstruct the missing populations based on the following hypothesis ensuring mass and momentum conservation at the co-located fine/coarse nodes

iΦificoarse,(1)=iΦififine,(1)=0,
(45)

where Φi=(1,cx,i,cy,i,cz,i)T.

As previously mentioned, in our case, there are no co-located fixed and moving nodes. Therefore, we propose to re-write Eq. (45) at each of the interface moving nodes xm (red circles in Fig. 2)

iΦifim,(1)=0.
(46)

Using the relation fi(1)=fifi(0)(X), where X=(ρ,ux,uy,uz) is the vector of relevant macroscopic variables, it yields

iΦi(fimfim,(0)(X))=0.
(47)

The populations fim are computed following the whole procedure of (i) geometrical rotation (ii) interpolation (iii) change of reference frame, exposed in Subsection III A, meaning that fim of Eq. (47) is equal to final f̂if of Subsection III A. This reads

iΦi(f̂iffim,(0)(X))=0.
(48)

The equilibrium function being non-linear on X, Eq. (48), is a four-dimensional non-linear system on X.

In this system, all distribution functions f̂if are known. The only unknown is X. Thus, solving Eq. (47) allows to obtain all macroscopic variables X=(ρ,ux,uy,uz) at the interface moving nodes xm. Consequently, the equilibrium fim,(0) and off-equilibrium fim,(1) populations can be reconstructed, leading to the reconstruction of the total moving populations at the interface through the relation

fim=fim,(0)+fim,(1),
(49)

where fim,(1)(xm)=g(1)(ρ¯m,u¯m,ΠIm,(1)) as assumed in Eq. (42) and ρ¯m,u¯m are, respectively, the newly obtained density and velocity after solving the non-linear system Eq. (47).

The four-dimensional non-linear system on X of Eq. (48) can be re-written as follows:

G(X)=0,
(50)

where G(X)=iΦi(f̂iffim,(0)(X)).

Equation (50) is solved using an iterative Newton–Raphson method. This method is detailed in Astoul et al. (2020) and is recalled here for the sake of completeness.

The Newton–Raphson method is a root-finding method, which produces successively better approximations to the roots of a real-valued set of equations, in this case Eq. (50), G(X)=0.

Let X0 be an initial guess for a root of G, and δX a small variation around the root X0. Then, the Newton–Raphson method provides a good approximation for the root as

X0+δX=X0JG1(X0)G(X0),
(51)

where JG(X0)=dG(X0)/dX is the Jacobian matrix of G evaluated at X0.

This algorithm can then be applied iteratively to obtain

n0,Xn+1=XnJG1(Xn)G(Xn).
(52)

The convergence criteria are set at the user's convenience; here, we consider the method as converged when ||δX||<1010.

In Eq. (52), the Jacobian matrix is obtained thanks to the formal computation software Maxima, which provides an analytical expression for it. Then, the inversion of the Jacobian matrix as well as the Newton–Raphson algorithm is implemented in the LB method code, namely the ProLB solver (www.prolb-cfd.com), which has been used to carry out this research.

Let us summarize the present direct coupling overset methodology. The method is composed of two parts.

The first part is the interpolation of distribution functions rather than macroscopic variables, which is believed to be more versatile and adequate to replace the current method with the conventional macroscopic value-based interpolations, as it is easier to couple it to other distribution function-based techniques based on the computation of high-order moments of the distribution functions. Additionally, it can lead to better conservation properties as all macroscopic variables are reconstructed based on the same distribution functions.

The second part is the solving of the non-linear system through the Newton–Raphson method to filter out the undesired zeroth- and first-order moments of the non-equilibrium distribution functions. This second part is responsible for the major improvements observed on the convected vortex test case as will be seen later on. On that account, the DC method would be especially recommended for problems where, at the fixed/moving grid interface, non-nullity of the zeroth- and first-order moments of the non-equilibrium populations is observed.

As for any new model, we have made the choice to validate the direct coupling overset-HRR algorithm first with academic test cases. The objective in this section is to show the numerical validations and the associated methodology, covering and separating different physical effects, from the simplest case to more complex configurations.

In this context, the three following aerodynamic test cases are first considered: an uniform flow (free shear-flow), a Poiseuille flow (wall-bounded flow), and an uniform flow past a rotating cylinder (complex flow including a moving object). Then, we address two aeroacoustic test cases: an acoustic pulse and a convected vortex. It should be noted here that the aerodynamic test cases serve a validation purpose, while the major improvements brought out by the present research are obtained for aeroacoustic applications (see the convected vortex test case below).

In this entire section, the results provided by the direct coupling approach are compared to ones obtained with the original method and to a reference simulation using a single fixed grid with no moving subdomain.

1. Uniform flow

A pseudo-2D uniform flow including an empty moving subdomain is considered to investigate the capability of the method to preserve uniform flows, with density ρref=1 and velocity uref=(0.1,0.0), leading to a Mach number Ma =0.17. A Dirichlet velocity boundary condition is imposed at the inlet (uinlet=uref) and a Dirichlet density boundary condition at the outlet (ρoutlet=ρref). The Reynolds number is Re =90. The computational domain is a box of dimensions Lx=100×Ly=100. It contains a fixed grid of the same dimensions and a moving grid of dimensions Lxm=14×Lym=14. The origin and center of rotation are, respectively, located at points (0, 0) and (50, 50). The moving grid rotates at the angular velocity ωz=0.01. Both fixed and moving grids within the computational domain along with a representation of the uniform velocity field are sketched on Fig. 3.

FIG. 3.

Sketch of the computational domain for the uniform flow, including both the fixed and the moving grids and a representation of the velocity magnitude field (uref=(0.1,0.0)).

FIG. 3.

Sketch of the computational domain for the uniform flow, including both the fixed and the moving grids and a representation of the velocity magnitude field (uref=(0.1,0.0)).

Close modal

The L1 norm and the L1 norm error of a given field X are defined, respectively, as

||X||1=1Ni=1N|X(xi)|
(53)

and

||EX||1=1Ni=1N|Xref(xi)X(xi)||Xref(xi)|,
(54)

where N is the number of nodes in the fixed grid.

The temporal evolution of the spatial L1 norm over the whole computational domain of the density field ||ρ||1 and the longitudinal velocity field ||ux||1 is plotted in Fig. 4. This L1 norm is computed in the fixed grid relative to the overset computation so that comparisons with the reference computation, which only includes a fixed grid, can be performed.

FIG. 4.

Time evolution of the shifted spatial L1 norm of the density field ||ρ||1 (left) and the velocity field ||ux||1 (right), for the uniform flow test case. The direct coupling overset method is compared to the original method and to a reference computation not containing any moving subdomain.

FIG. 4.

Time evolution of the shifted spatial L1 norm of the density field ||ρ||1 (left) and the velocity field ||ux||1 (right), for the uniform flow test case. The direct coupling overset method is compared to the original method and to a reference computation not containing any moving subdomain.

Close modal

It can be observed that the direct coupling method results agree very closely with the reference results. More precisely, with regard to the reference computation, the errors write ||Eρ||1=2.86·107 % for the density and ||Eux||1=6.54·106 % for the velocity at the last time step.

Additionally, leaving out the machine precision error, no difference is observed between the direct coupling and the original spatial L1 norm computed values for both density and velocity (|||ρ||1Original||ρ||1DC|tfinal=|||ux||1Original||ux||1DC|tfinal=0).

This result was expected and can be explained by the following two points:

  • First, it is worth noting that in the special case of the uniform flow (mathematically, this case is equivalent to the macroscopic vector X being a constant in 4), the interpolation operator is commutable with the distribution function. This means that interpolating the distribution functions then computing the macroscopic values provides the very same result as interpolating the macroscopic values then reconstructing the distribution functions. In other words, the reconstructed distribution function from the original model [Eq. (26)] is equal to the distribution function obtained after the rotation/interpolation/change of reference frame process [Eq. (37)].

  • Second, during our investigations, we have computed the zeroth- and first-order moments of the non-equilibrium function at the moving/fixed interfaces, and for the uniform flow test case, they are found to be of the order of 1018 using the original model regardless. Since the purpose of the direct coupling approach is to impose the condition of nullity of both these zeroth- and first-order moments, and given that for the uniform flow both distributions from Eqs. (26) and (37) are equal [see point (i) above], it can be deduced that the effect of the Newton–Raphson iteration will be negligible. This is why, for the special case of the uniform flow, both original and direct coupling methods end up providing the same results.

Figure 5 shows the density and velocity vertical profiles upstream (x =15) and downstream (x =85) of the rotating region, at the last time step. It can be seen that both overset methods (original and direct coupling) provide results that are in a very good agreement with the ones obtained by the reference computation. The mean and maximum L1 norm errors on these profiles with respect to the reference simulation are drawn in Table I. Full data with maximum precision are shown here as the differences between the computations are tiny in all test cases. It can be seen that the error overall does not exceed 5.960×106 % for density and 1.490×104 % for velocity.

FIG. 5.

Density (top) and longitudinal velocity (bottom) profiles upstream (x =15) and downstream (x =85) of the rotating region, for the uniform flow test case. The direct coupling overset method is compared to the original method and to a reference computation not containing any moving subdomain.

FIG. 5.

Density (top) and longitudinal velocity (bottom) profiles upstream (x =15) and downstream (x =85) of the rotating region, for the uniform flow test case. The direct coupling overset method is compared to the original method and to a reference computation not containing any moving subdomain.

Close modal
TABLE I.

Mean and maximum L1 norm error (%) on the density and longitudinal velocity vertical profiles of Fig. 5, for the original and DC overset methods and for the uniform flow test case.

x =15Original methodDC method
Mean(||Eρ||1) (%) 1.668800001297654×106 1.668800001297654×106 
Max(||Eρ||1) (%) 5.960000004634480×106 5.960000004634480×106 
Mean(||Eux||1) (%) 1.362759979612352×105 1.362759979612352×105 
Max(||Eux||1) (%) 2.235999966892892×105 2.235999966892892×105 
x =15Original methodDC method
Mean(||Eρ||1) (%) 1.668800001297654×106 1.668800001297654×106 
Max(||Eρ||1) (%) 5.960000004634480×106 5.960000004634480×106 
Mean(||Eux||1) (%) 1.362759979612352×105 1.362759979612352×105 
Max(||Eux||1) (%) 2.235999966892892×105 2.235999966892892×105 
x =85Original methodDC method
Mean(||Eρ||1) (%) 1.311200001019586×106 1.311200001019586×106 
Max(||Eρ||1) (%) 5.960000004634480×106 5.960000004634480×106 
Mean(||Eux||1) (%) 4.357639934587714×105 4.357639934587714×105 
Max(||Eux||1) (%) 1.489999977559619×104 1.489999977559619×104 
x =85Original methodDC method
Mean(||Eρ||1) (%) 1.311200001019586×106 1.311200001019586×106 
Max(||Eρ||1) (%) 5.960000004634480×106 5.960000004634480×106 
Mean(||Eux||1) (%) 4.357639934587714×105 4.357639934587714×105 
Max(||Eux||1) (%) 1.489999977559619×104 1.489999977559619×104 

2. Poiseuille flow

We now consider a pseudo-2D Poiseuille flow, which adds more complexity with respect to the uniform flow test case, due to the presence of velocity gradients. Both fixed and moving grids within the computational domain have the same size as for the uniform flow test case, that is, respectively, Lx=100×Ly=100 and Lxm=14×Lym=14. The origin and center of rotation are, respectively, located at points (0, 0) and (50, 50). The moving grid rotates at the angular velocity ωz=0.01. No-slip boundary conditions are imposed on the upper and lower walls, while at the inlet a velocity Poiseuille profile is imposed as follows (Re =90, Ma =0.17):

ux,inlet=ux,max(4yLy4y2Ly2),
(55)
uy,inlet=0,
(56)
ρinlet=1,
(57)

where ux,max=0.1.

The corresponding computational domain velocity magnitude field is sketched in Fig. 6.

FIG. 6.

Sketch of the computational domain for the Poiseuille flow, including both the fixed and the moving grids and a representation of the velocity magnitude field as imposed in Eq. (55).

FIG. 6.

Sketch of the computational domain for the Poiseuille flow, including both the fixed and the moving grids and a representation of the velocity magnitude field as imposed in Eq. (55).

Close modal

The temporal evolutions of the spatial L1 norm of the density field ||ρ||1 and the longitudinal velocity field ||ux||1 are plotted in Fig. 7. Similarly to the uniform flow test case, the DC method's results agree very closely with the ones obtained by the original method simulation: at the last time step, the absolute difference between the two L1 norm computations for density and velocity is, respectively, |||ρ||1Original||ρ||1DC|tfinal=3.216×1010 and |||ux||1Original||ux||1DC|tfinal=4.894×1011. In the end, the following error values are found when compared to the reference computation: ||Eρ||1=7.89×105 % for the density and ||Eux||1=7.10×105 % for the velocity at the last time step. These low levels of error support the fact that both original and DC overset algorithms are thoroughly able to pass through the moving subdomain without altering the quality of the channel flow solution.

FIG. 7.

Time evolution of the spatial L1 norm of the density field ||ρ||1 (left) and the velocity field ||ux||1 (right), for the Poiseuille flow test case. The direct coupling overset method is compared to the original method and to a reference computation not containing any moving subdomain.

FIG. 7.

Time evolution of the spatial L1 norm of the density field ||ρ||1 (left) and the velocity field ||ux||1 (right), for the Poiseuille flow test case. The direct coupling overset method is compared to the original method and to a reference computation not containing any moving subdomain.

Close modal

Figure 8 shows the density and velocity vertical profiles upstream (x =15) and downstream (x =85) of the rotating region, at the last time step. It can be seen that both methods provide very accurate results, which are in a satisfactory agreement with the reference simulation results. The mean and maximum L1 norm errors on all these profiles drawn in Table II indeed show that the error overall does not exceed 1.547×104 % for density and 2.430×102 % for velocity.

FIG. 8.

Density (top) and longitudinal velocity (bottom) vertical profiles upstream (x =15) and downstream (x =85) of the rotating region, for the Poiseuille flow test case. The direct coupling overset method is compared to the original method and to a reference computation not containing any moving subdomain.

FIG. 8.

Density (top) and longitudinal velocity (bottom) vertical profiles upstream (x =15) and downstream (x =85) of the rotating region, for the Poiseuille flow test case. The direct coupling overset method is compared to the original method and to a reference computation not containing any moving subdomain.

Close modal
TABLE II.

Mean and maximum L1 norm error (%) on the density and longitudinal velocity vertical profiles of Fig. 8, for the original and DC overset methods and for the Poiseuille flow test case.

x =15Original methodDC method
Mean(||Eρ||1) (%) 1.176583906798957×104 1.174208598866126×104 
Max(||Eρ||1) (%) 1.546827095575228×104 1.428284107822683×104 
Mean(||Eux||1) (%) 7.881696528406885×104 7.880878689514660×104 
Max(||Eux||1) (%) 3.638914170708×103 3.638914170708×103 
x =15Original methodDC method
Mean(||Eρ||1) (%) 1.176583906798957×104 1.174208598866126×104 
Max(||Eρ||1) (%) 1.546827095575228×104 1.428284107822683×104 
Mean(||Eux||1) (%) 7.881696528406885×104 7.880878689514660×104 
Max(||Eux||1) (%) 3.638914170708×103 3.638914170708×103 
x =85Original methodDC method
Mean(||Eρ||1) (%) 2.955730755523385×105 2.955730755523385×105 
Max(||Eρ||1) (%) 7.157658657838396×105 7.157658657838396×105 
Mean(||Eux||1) (%) 1.3648512435990×102 1.3648230615681×102 
Max(||Eux||1) (%) 2.4302884137810×102 2.4302884137810×102 
x =85Original methodDC method
Mean(||Eρ||1) (%) 2.955730755523385×105 2.955730755523385×105 
Max(||Eρ||1) (%) 7.157658657838396×105 7.157658657838396×105 
Mean(||Eux||1) (%) 1.3648512435990×102 1.3648230615681×102 
Max(||Eux||1) (%) 2.4302884137810×102 2.4302884137810×102 

3. Uniform flow past a rotating cylinder

In this section, additional complexity is introduced again by considering the case of the flow past a rotating circular cylinder. The computational domain is a box of size Lx=798×Ly=798. It includes a circular moving grid of radius 30. The moving grid rotates at the angular velocity ωz=0.0017. The circular cylinder is of radius R =10. It is attached to the moving grid and therefore rotates in time at the same angular velocity rate as the moving grid. At the inlet, an uniform velocity profile is imposed: uref=(U=0.017,0.0), while at the outlet ρref=1. The viscosity is imposed as ν=1.70×103, which yields a Reynolds number of Re =200. Defining the rotation ratio a as the ratio between the local velocity at the rotating cylinder wall and the free stream velocity (a=Rωz/U), the latter physical values yield a =1.0. Both fixed and moving grids within the computational domain are sketched on Fig. 9, which also shows the velocity field. In addition, velocity streamlines and vorticity field are shown in Fig. 10. The velocity streamlines clearly show the disruption of wake symmetry due to the cylinder rotation, while the vorticity field puts forward the generation of the Kármán vortex structure.

FIG. 9.

Closeup at the fixed and the moving grids for the uniform flow past a rotating cylinder test case. The velocity magnitude field is also represented.

FIG. 9.

Closeup at the fixed and the moving grids for the uniform flow past a rotating cylinder test case. The velocity magnitude field is also represented.

Close modal
FIG. 10.

Velocity streamlines (left) and vorticity field (right) for the uniform flow past a rotating cylinder test case.

FIG. 10.

Velocity streamlines (left) and vorticity field (right) for the uniform flow past a rotating cylinder test case.

Close modal

Lift and drag coefficients have been computed at the end of the simulation using the far-field integral method (Toubin and Bailly, 2015; Gariepy et al., 2013; Wilhelm et al., 2018) and are shown in Table III. They are compared with Navier–Stokes direct numerical simulation (DNS) results given by Mittal and Kumar (2003). It can be seen that the aerodynamic coefficients obtained by both the original and DC methods are in a good agreement with the reference DNS results. In addition, the pressure coefficient Cp on the cylinder surface has been plotted in Fig. 11. Again, the results obtained by both overset methods compare very well with the reference DNS results of Mittal and Kumar (2003). The small observed differences originates in the post-processing step of reconstruction of the flow variables on the solid surface since a body-fitted grid is not used in the present LBM simulations (Cai et al., 2021).

TABLE III.

Aerodynamic coefficients (uniform flow past a rotating cylinder test case, Re =200, a =1).

Original methodDC methodMittal and Kumar (2003) 
CL 1.168 043 262 926 71 1.168 034 110 111 57 1.068 
CD 2.476 433 333 358 73 2.476 384 452 553 85 2.4 
Original methodDC methodMittal and Kumar (2003) 
CL 1.168 043 262 926 71 1.168 034 110 111 57 1.068 
CD 2.476 433 333 358 73 2.476 384 452 553 85 2.4 
FIG. 11.

Pressure coefficient Cp on the cylinder surface for the uniform flow past a rotating cylinder test case. Original and DC methods are in a very close agreement and compare very well with Mittal and Kumar (2003) DNS reference results.

FIG. 11.

Pressure coefficient Cp on the cylinder surface for the uniform flow past a rotating cylinder test case. Original and DC methods are in a very close agreement and compare very well with Mittal and Kumar (2003) DNS reference results.

Close modal

Figure 12 shows the density and velocity vertical profiles upstream (x/D=10) and downstream (x/D=10, where D =20 is the diameter of the cylinder) of the rotating cylinder. It can be seen that the original and DC methods agree very closely with each other, although slight differences are noticeable in comparison with the reference computation. These differences have been quantified in Table IV, which shows the mean and maximum L1 norm errors on the profiles with respect to the reference simulation. It can be seen that the error overall does not exceed 4.167×103 % for density and 6.813% for the velocity magnitude. Therefore, both original and DC methods still compare fairly well to the reference results considering the increased complexity along with a much higher Reynolds number brought out by the flow past a rotating cylinder test case.

FIG. 12.

Density (top), longitudinal velocity (center), and transverse velocity (bottom) vertical profiles upstream (x/D=10) and downstream (x/D=10) of the rotating region, for the uniform flow past a rotating cylinder test case. Original and DC methods are in a very close agreement.

FIG. 12.

Density (top), longitudinal velocity (center), and transverse velocity (bottom) vertical profiles upstream (x/D=10) and downstream (x/D=10) of the rotating region, for the uniform flow past a rotating cylinder test case. Original and DC methods are in a very close agreement.

Close modal
TABLE IV.

Mean and maximum L1 norm error (%) on the density and longitudinal velocity vertical profiles of Fig. 12, for the original and DC overset methods and for the uniform flow past a rotating cylinder test case.

x/D=10Original methodDC method
Mean(||Eρ||1) (%) 2.9803635257429236×104 2.9803635257429236×104 
Max(||Eρ||1) (%) 8.333402778411081×104 8.333402778411081×104 
Mean(||E||u||||1) (%) 3.799049107179938×102 3.7991061361194046×102 
Max(||E||u||||1) (%) 6.264432073432659×102 6.264 431 322 722 924 ×102 
x/D=10Original methodDC method
Mean(||Eρ||1) (%) 2.9803635257429236×104 2.9803635257429236×104 
Max(||Eρ||1) (%) 8.333402778411081×104 8.333402778411081×104 
Mean(||E||u||||1) (%) 3.799049107179938×102 3.7991061361194046×102 
Max(||E||u||||1) (%) 6.264432073432659×102 6.264 431 322 722 924 ×102 
x/D=10Original methodDC method
Mean(||Eρ||1) (%) 5.436273590404883×104 5.436273590404883×104 
Max(||Eρ||1) (%) 4.166805560179329×103 4.166805560179329×103 
Mean(||E||u||||1) (%) 5.67788198711541×101 5.679093365767347×101 
Max(||E||u||||1) (%) 6.812 550 348 059 371 6.813 554 108 265 357 
x/D=10Original methodDC method
Mean(||Eρ||1) (%) 5.436273590404883×104 5.436273590404883×104 
Max(||Eρ||1) (%) 4.166805560179329×103 4.166805560179329×103 
Mean(||E||u||||1) (%) 5.67788198711541×101 5.679093365767347×101 
Max(||E||u||||1) (%) 6.812 550 348 059 371 6.813 554 108 265 357 

Additionally, if one calculates the absolute difference values between the two spatial L1 norm computations for density and velocity over the whole domain and at the last time step, tiny difference is found: |||ρ||1Original||ρ||1DC|tfinal=1.141×109 and |||u||1Original||u||1DC|=4.852×106.

In conclusion, these low levels of error observed using the DC method (interpolating distribution functions along with a Newton–Raphson algorithm) allow to validate it as an alternative method to the original one (interpolating macroscopic values) for the computation and analysis of laminar flows past rotating solid bodies.

1. Acoustic pulse

A pseudo-2D acoustic pulse including an empty moving subdomain is now considered. The computational domain is a periodic box of dimensions Lx=100×Ly=100. It contains a fixed grid of the same dimensions and a moving grid of dimensions Lxm=39×Lym=39. The origin and center of rotation are, respectively, located at points (0, 0) and (50, 50). The moving grid rotates at the angular velocity ωz=0.00059, and the kinematic viscosity is ν=1×104. The acoustic pulse is initialized at the center of the box as follows:

ρ(x,y,z)=1+A((xxc)2+(yyc)22Rc2),
(58)
u(x,y,z)=0,
(59)

where A =0.001, xc=yc=50, and Rc = 10.

Both fixed and moving grids within the computational domain along with a representation of the density field at an arbitrary time step during the simulation are sketched on Fig. 13.

FIG. 13.

Sketch of the computational domain for the acoustic pulse test case, including both the fixed and the moving grids and a representation of the density field.

FIG. 13.

Sketch of the computational domain for the acoustic pulse test case, including both the fixed and the moving grids and a representation of the density field.

Close modal

Figure 14 shows the horizontal profiles located at the center of the domain (y =50) of density, both components of velocity ux and uy and velocity divergence, at the last time step. It can be seen that great agreement is achieved between the results of both overset methods and the reference computation. The velocity divergence profile is fairly smooth, and no parasitic distorsions are observed whatsoever. The presence of a rotating domain is thus not affecting the algorithm, and the acoustic wave is very accurately propagated from the moving to the fixed grid.

FIG. 14.

Density, velocity magnitude, and velocity divergence horizontal profiles (y =0.5), for the acoustic pulse test case. The direct coupling overset method is compared to the original method and to a reference computation not containing any moving subdomain.

FIG. 14.

Density, velocity magnitude, and velocity divergence horizontal profiles (y =0.5), for the acoustic pulse test case. The direct coupling overset method is compared to the original method and to a reference computation not containing any moving subdomain.

Close modal

For a more precise analysis, as with the previous test cases, quantified error values are drawn in Table V. It has been chosen here to show only absolute difference errors rather than errors in % as defined in Eq. (54). Indeed, given that the reference data set shows at some points very tiny values (that tend to zero), using the error in % of Eq. (54) type of definition leads to very high error values that have poor physical meaning. It can be seen that the error overall does not exceed 1.193·106 for density and 6.367·107 for velocity whatever the algorithm. These levels of error remain fairly low and very similar whatever the algorithm, which allows to validate both original and direct coupling overset methods as adequate for the computation of an acoustic pulse with a rotating subdomain.

TABLE V.

Mean and maximum absolute errors between the overset computations (original and DC) and the reference computation on the density and velocity horizontal profiles of Fig. 14, for the acoustic pulse test case.

Original methodDC method
Mean|ρ1ρ1REF| 8.126333333202674×107 8.126333333202674×107 
Max| ρ1ρ1REF| 1.192500000080088×106 1.192500000080088×106 
Mean(|u1||u||1REF|) 3.783709888991750×107 3.783664374708644×107 
Max(|u1||u||1REF|) 6.367026105555673×107 6.366852250907429×107 
Original methodDC method
Mean|ρ1ρ1REF| 8.126333333202674×107 8.126333333202674×107 
Max| ρ1ρ1REF| 1.192500000080088×106 1.192500000080088×106 
Mean(|u1||u||1REF|) 3.783709888991750×107 3.783664374708644×107 
Max(|u1||u||1REF|) 6.367026105555673×107 6.366852250907429×107 

2. Convected vortex

More complexity is now introduced by addressing the case of a convected vortex. This test case was previously shown to be likely to generate spurious acoustics when crossing a plane refinement interface in Gendre et al. (2017) and Astoul et al. (2020). Additionally, it has been observed across our investigations that the original overset algorithm may fail at conserving the vortex structure and consistency when the vortex crosses a rotating subdomain. For these reasons, the convected vortex is believed to be a challenging test case to assess the efficiency of our new algorithm.

The computational domain for the pseudo-2D convected vortex test case is a periodic box of dimensions Lx=200×Ly=200. It contains a fixed grid of the same dimensions and a moving grid of dimensions Lxm=79×Lym=79. The origin and center of rotation are, respectively, located at points (0, 0) and (100, 100). The moving grid rotates at the angular velocity ωz=0.006, and the kinematic viscosity is ν=1·1010. The convected vortex is initialized at the center of the box as follows (Re=4·1010, Ma =0.17):

ρ(x,y,z)=ρ0exp(ϵ22cs2exp((xxc)2+(yyc)2Rc2)),
(60)
ux(x,y,z)=Uxϵ(yycRc)exp((xxc)2+(yyc)22Rc2),
(61)
uy(x,y,z)=ϵ(xxcRc)exp((xxc)2+(yyc)22Rc2),
(62)
uz(x,y,z)=0,
(63)

where ρ0=1,Ux=0.1,ϵ=0.15,xc=yc=100, and Rc = 20.

Both fixed and moving grids within the computational domain along with a representation of the density field at an arbitrary time step during the simulation are sketched on Fig. 15.

FIG. 15.

Sketch of the computational domain for the convected vortex test case, including both the fixed and the moving grids and a representation of the density field.

FIG. 15.

Sketch of the computational domain for the convected vortex test case, including both the fixed and the moving grids and a representation of the density field.

Close modal

Figure 16 shows the horizontal profiles located at the center of the domain (y =50) of density, both components of velocity ux and uy and velocity divergence, at the last time step. It is worth noting that this corresponds to 20 times the advection time of the vortex, meaning that at the last time step, the vortex has crossed the rotating subdomain 20 times. Figure 16 reveals that while a poor agreement is found between the original overset method and the reference computation, excellent results are obtained with the direct coupling approach. The direct coupling method shows a good ability to conserve the vortex structure, as opposed to both the original approach. This phenomena can also been observed in Fig. 17, which shows the density contours generated by the original, the direct coupling, and the reference computations, at t =0, t=10tadv, and t=20tadv, where tadv is the advection time. It can clearly be seen that using the direct coupling method, the intensity of the vortex does not fade away as rapidly as with the original overset method.

FIG. 16.

Density, longitudinal velocity ux, transverse velocity uy, and velocity divergence horizontal profiles (y =50), for the acoustic pulse test case. The direct coupling overset method is compared to the original method and to a reference computation not containing any moving subdomain.

FIG. 16.

Density, longitudinal velocity ux, transverse velocity uy, and velocity divergence horizontal profiles (y =50), for the acoustic pulse test case. The direct coupling overset method is compared to the original method and to a reference computation not containing any moving subdomain.

Close modal
FIG. 17.

Closeup at density contours generated by the original (left), the direct coupling (center), and the reference (right) computations, at t =0 (top), t=10tadv (center), and t=20tadv (bottom), where tadv is the advection time. The direct coupling method allows a better conservation of the vortex structure and consistency than the original overset method. Color maps are tightened to enhance visualization of the vortex diffusion over time.

FIG. 17.

Closeup at density contours generated by the original (left), the direct coupling (center), and the reference (right) computations, at t =0 (top), t=10tadv (center), and t=20tadv (bottom), where tadv is the advection time. The direct coupling method allows a better conservation of the vortex structure and consistency than the original overset method. Color maps are tightened to enhance visualization of the vortex diffusion over time.

Close modal

Additionally, the velocity divergence horizontal profiles in Fig. 16 of reveals another interesting phenomena: high parasitic distorsions are created at the vicinity of the moving grid interface for both overset methods. This shows that spurious acoustics are generated by the presence of the rotating domain. It can be observed that the direct coupling approach allows to reduce the amplitude of these distorsions in comparison with the original method.

For a more precise analysis, quantified L1 norm error values are drawn in Table VI. For the same reasons as for the acoustic pulse test cases, it has been chosen here to show only absolute difference errors rather than errors in % as defined in Eq. (54). It can be seen that the direct coupling method allows a clear reduction of the error on all the profiles: the maximum absolute errors on density, longitudinal and transverse velocity are, respectively, about 1.76×,6.18×,2.04× lower using the direct coupling approach over the original method. Additionally, we have quantified the amplitude of the parasitic distorsions at the vicinity of the moving grid interface that are observed in Fig. 16. The distortion amplitudes are 1.16× and 2.11× lower using the direct coupling approach over the original method, for the left one (x =25) and the right one (x =76), respectively.

TABLE VI.

Mean and maximum absolute errors between the overset computations (original and DC) and the reference computation on the density and velocity horizontal profiles of Fig. 16, for the convected vortex test case.

Original methodDC method
Mean|ρ1ρ1REF| 9.333193049999910×104 4.074471959999959×104 
Max| ρ1ρ1REF| 4.724383400000×103 2.678275100000×103 
Mean(|ux1||ux||1REF|) 1.530553475055×103 3.368238869740786×104 
Max(|ux1||ux||1REF|) 4.819702046274×103 7.792670588818718×104 
Mean(|uy1||uy||1REF|) 4.765100393857×103 2.212539263029×103 
Max(|uy1||uy||1REF|) 1.9273253639242×102 9.425645768786×103 
Original methodDC method
Mean|ρ1ρ1REF| 9.333193049999910×104 4.074471959999959×104 
Max| ρ1ρ1REF| 4.724383400000×103 2.678275100000×103 
Mean(|ux1||ux||1REF|) 1.530553475055×103 3.368238869740786×104 
Max(|ux1||ux||1REF|) 4.819702046274×103 7.792670588818718×104 
Mean(|uy1||uy||1REF|) 4.765100393857×103 2.212539263029×103 
Max(|uy1||uy||1REF|) 1.9273253639242×102 9.425645768786×103 

In this section, we want to investigate the sensitivity of the previously shown errors to the moving grid local Mach number, which involves the simulation parameter ω. The moving grid maximum local Mach number MaMG,max is defined as follows:

MaMG,max=||vr,max||cs=rmax||ω||cs.
(64)

For all test cases studied in this work except the uniform flow past a rotating cylinder, a square moving mesh has been used, which means that the grid location leading to the maximum velocity is the corner of the moving grid: in other words, rmax=lMG2/2, where lMG is the side length of the square moving mesh. For this sensitivity study, the following values for ω have been chosen: ω1,ω2=2ω1,ω3=3ω1,ω4=4ω1,ω5=5ω1,ω6=6ω1. The corresponding MaMG,max values are displayed in Table VII.

TABLE VII.

Values of the moving grid maximum local Mach number MaMG,max for the four considered test cases.

MaMG,max,1MaMG,max,2MaMG,max,3MaMG,max,4MaMG,max,5MaMG,max,6
Uniform flow 0.029 0.034 0.043 0.057 0.086 0.17 
Poiseuille flow 0.029 0.034 0.043 0.057 0.086 0.17 
Acoustic pulse 0.0049 0.0058 0.0073 0.0097 0.015 0.029 
Convected vortex 0.049 0.059 0.074 0.098 0.15 0.29 
MaMG,max,1MaMG,max,2MaMG,max,3MaMG,max,4MaMG,max,5MaMG,max,6
Uniform flow 0.029 0.034 0.043 0.057 0.086 0.17 
Poiseuille flow 0.029 0.034 0.043 0.057 0.086 0.17 
Acoustic pulse 0.0049 0.0058 0.0073 0.0097 0.015 0.029 
Convected vortex 0.049 0.059 0.074 0.098 0.15 0.29 

In Fig. 18, we have plotted the evolution of the mean error ||EX||1 [as defined in Eq. (54)] against the local Mach number as defined in Eq. (64). It should be noted that in this figure, the error is expressed in % for the uniform and Poiseuille flow test cases, while absolute error values are given for the acoustic pulse and convected vortex test cases. Indeed, as explained in Sec. IV B 1, the reference datasets for these test cases show at some points very tiny values that in consequence would lead to very high error values with poor physical meaning if quantifying the error in %.

FIG. 18.

Evolution of the mean error ||EX||1 [as defined in Eq. (54)] against the local Mach number for density (left) and velocity (right) and for the four considered test cases (from top to bottom: uniform flow, Poiseuille flow, acoustic pulse, convected vortex).

FIG. 18.

Evolution of the mean error ||EX||1 [as defined in Eq. (54)] against the local Mach number for density (left) and velocity (right) and for the four considered test cases (from top to bottom: uniform flow, Poiseuille flow, acoustic pulse, convected vortex).

Close modal

First, it can be seen overall there is no significant influence of the rotating grid velocity on the calculated error on the macroscopic values. No specific pattern is observed in the evolution of the error, as depending on the test case it might decrease or increase depending on the moving grid velocity. The error always remains in the same order of magnitude, which means a fair robustness for the method.

The second observation is that, as shown in Sec. IV B 2, there is a significant gap between the error obtained with the original method and the one obtained using the direct coupling method with regard to the convected vortex test case. It can also be observed here that this difference is even more enhanced for high ω values. In other words, the higher the moving grid velocity, the more beneficial the use of the direct coupling method over the original one.

Finally, a mesh convergence study on the direct coupling method has been conducted for the convected vortex test case, based on the L2 and infinity norms, respectively, defined as

||EX||2=i(X(xi)Xref(xi))2iXref2(xi),
(65)
||EX||=maxi|X(xi)Xref(xi)|maxi|Xref(xi)|.
(66)

The L2 and infinity norms of both density and velocity for three different mesh resolutions (number of cells N1=9559,N2=4N1,N3=16N1) are shown in Fig. 19. A dashed line of slope 2 shows that the convergence is found to be second order, which was theoretically expected.

FIG. 19.

Evolution of the L2 norm ||EX||2 (left) and infinity norm ||EX|| (right)—as defined in Eqs. (65) and (66)—of the density and velocity fields computed by the direct coupling method against the number of cells Ncells, for the convected vortex test case. A dashed line of slope 2 shows that the convergence is found to be second order.

FIG. 19.

Evolution of the L2 norm ||EX||2 (left) and infinity norm ||EX|| (right)—as defined in Eqs. (65) and (66)—of the density and velocity fields computed by the direct coupling method against the number of cells Ncells, for the convected vortex test case. A dashed line of slope 2 shows that the convergence is found to be second order.

Close modal

In this paper, we present a new direct coupling scheme based on the overset technique to tackle moving boundary problems within the lattice Boltzmann framework. As opposed to previous studies on the overset technique, the approach here is based on the interpolation of distribution functions rather than moments at the borders between the moving and fixed grids. The present direct coupling approach also includes additional hypothesis ensuring mass and momentum conservation at the interface nodes between fixed and moving grids. The method has been validated with three aerodynamic test cases (an uniform flow, a Poiseuille flow, and an uniform flow past a rotating cylinder) and two aeroacoustic test cases (an acoustic pulse and a convected vortex). The results are compared to the ones obtained by a reference computation that does not contain any moving subdomain. It has been shown that the direct coupling method results compare very well with the original overset method results as well as the reference computation results. Additionally, it has been demonstrated that the direct coupling method allows an improvement of the accuracy of the lattice Boltzmann overset algorithm for aeroacoustics. In particular, the convected vortex test case reveals that a better ability to conserve the vortex structure over time, as well as a reduction of spurious acoustic distorsions at the fixed/moving interface, are provided by the direct coupling approach. Extensions to turbulent and thermal flows are the subject of further investigations.

This work was supported by ANR Industrial Chair ALBUMS (Grant No. ANR-18-CHIN-0003-01). The ProLB solver was used to perform the work.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
Aidun
,
C. K.
, and
Clausen
,
J. R.
, “
Lattice-Boltzmann method for complex flows
,”
Annu. Rev. Fluid Mech.
42
,
439
472
(
2010
).
2.
Astoul
,
T.
,
Wissocq
,
G.
,
Boussuge
,
J.-F.
,
Sengissen
,
A.
, and
Sagaut
,
P.
, “
Lattice Boltzmann method for computational aeroacoustics on non-uniform meshes: A direct grid coupling approach
,” arXiv:2004.14887 (
2020
).
3.
Bhatnagar
,
P. L.
,
Gross
,
E. P.
, and
Krook
,
M.
, “
A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems
,”
Phys. Rev.
94
,
511
(
1954
).
4.
Cai
,
S.
,
Degrigny
,
J.
,
Boussuge
,
J.
, and
Sagaut
,
P.
, “
Coupling turbulence wall models and immersed boundaries on Cartesian grids
,”
J. Comput. Phys.
429
,
109995
(
2021
).
5.
Coreixas
,
C. G.
, “
High-order extension of the recursive regularized lattice Boltzmann method
,” Ph.D. thesis (
Institute National Polytechnique de Toulouse
,
2018
).
6.
Far
,
E. K.
,
Geier
,
M.
, and
Krafczyk
,
M.
, “
Simulation of rotating objects in fluids with the cumulant lattice Boltzmann model on sliding meshes
,”
Comput. Math. Appl.
79
,
3
16
(
2020
).
7.
Far
,
E. K.
,
Geier
,
M.
,
Kutscher
,
K.
, and
Krafczyk
,
M.
, “
Simulation of micro aggregate breakage in turbulent flows by the cumulant lattice Boltzmann method
,”
Comput. Fluids
140
,
222
231
(
2016
).
8.
Favier
,
J.
,
Revell
,
A.
, and
Pinelli
,
A.
, “
A lattice Boltzmann-immersed boundary method to simulate the fluid interaction with moving and slender flexible objects
,”
J. Comput. Phys.
261
,
145
161
(
2014
).
9.
Feng
,
Y.
,
Boivin
,
P.
,
Jacob
,
J.
, and
Sagaut
,
P.
, “
Hybrid recursive regularized thermal lattice Boltzmann model for high subsonic compressible flows
,”
J. Comput. Phys.
394
,
82
99
(
2019
).
10.
Feng
,
Z.-G.
, and
Michaelides
,
E. E.
, “
The immersed boundary-lattice Boltzmann method for solving fluid–particles interaction problems
,”
J. Comput. Phys.
195
,
602
628
(
2004
).
11.
Filippova
,
O.
, and
Hänel
,
D.
, “
A novel lattice BGK approach for low Mach number combustion
,”
J. Comput. Phys.
158
,
139
160
(
2000
).
12.
Gariepy
,
M.
,
Trepanier
,
J.-Y.
, and
Malouin
,
B.
, “
Generalization of the far-field drag decomposition method to unsteady flows
,”
AIAA J.
51
,
1309
1319
(
2013
).
13.
Gendre
,
F.
,
Ricot
,
D.
,
Fritz
,
G.
, and
Sagaut
,
P.
, “
Grid refinement for aeroacoustics in the lattice Boltzmann method: A directional splitting approach
,”
Phys. Rev. E
96
,
023311
(
2017
).
14.
Guo
,
S.
,
Feng
,
Y.
, and
Sagaut
,
P.
, “
Improved standard thermal lattice Boltzmann model with hybrid recursive regularization for compressible laminar and turbulent flows
,”
Phys. Fluids
32
,
126108
(
2020
).
15.
Guo
,
Z.
,
Zheng
,
C.
, and
Shi
,
B.
, “
Discrete lattice effects on the forcing term in the lattice Boltzmann method
,”
Phys. Rev. E
65
,
046308
(
2002
).
16.
Jacob
,
J.
,
Malaspinas
,
O.
, and
Sagaut
,
P.
, “
A new hybrid recursive regularised Bhatnagar–Gross–Krook collision model for lattice Boltzmann method-based large eddy simulation
,”
J. Turbul.
19
,
1051
1076
(
2018
).
17.
Krüger
,
T.
,
Kusumaatmaja
,
H.
,
Kuzmin
,
A.
,
Shardt
,
O.
,
Silva
,
G.
, and
Viggen
,
E. M.
,
The Lattice Boltzmann Method
(
Springer International Publishing
,
2017
), Vol.
10
, pp.
4
15
.
19.
Lallemand
,
P.
, and
Luo
,
L.-S.
, “
Lattice Boltzmann equation with overset method for moving objects in two-dimensional flows
,”
J. Comput. Phys.
407
,
109223
(
2020
).
20.
Latt
,
J.
, and
Chopard
,
B.
, “
Lattice Boltzmann method with regularized pre-collision distribution functions
,”
Math. Comput. Simul.
72
,
165
168
(
2006
).
21.
Li
,
Y.
, “
An improved volumetric LBM boundary approach and its extension for sliding mesh simulation
,” Ph.D. Dissertation (Iowa State University, 2011).
22.
Ma
,
C.
,
Wu
,
J.
, and
Zhang
,
T.
, “
A high order spectral difference-based phase field lattice Boltzmann method for incompressible two-phase flows
,”
Phys. Fluids
32
,
122113
(
2020
).
23.
Malaspinas
,
O.
, “
Increasing stability and accuracy of the lattice Boltzmann scheme: Recursivity and regularization
,” arXiv:1505.06900 (
2015
).
24.
Meldi
,
M.
,
Vergnault
,
E.
, and
Sagaut
,
P.
, “
An arbitrary Lagrangian–Eulerian approach for the simulation of immersed moving solids with lattice Boltzmann method
,”
J. Comput. Phys.
235
,
182
198
(
2013
).
25.
Mittal
,
S.
, and
Kumar
,
B.
, “
Flow past a rotating cylinder
,”
J. Fluid Mech.
476
,
303
334
(
2003
).
26.
Peskin
,
C. S.
, “
Flow patterns around heart valves: A numerical method
,”
J. Comput. Phys.
10
,
252
271
(
1972
).
27.
Peskin
,
C. S.
, “
The immersed boundary method
,”
Acta Numer.
11
,
479
517
(
2002
).
28.
Saadat
,
M. H.
, and
Karlin
,
I. V.
, “
Arbitrary Lagrangian–Eulerian formulation of lattice Boltzmann model for compressible flows on unstructured moving meshes
,”
Phys. Fluids
32
,
046105
(
2020
).
29.
Sagaut
,
P.
, “
Toward advanced subgrid models for lattice-Boltzmann-based large-eddy simulation: Theoretical formulations
,”
Comput. Math. Appl.
59
,
2194
2199
(
2010
).
30.
Shan
,
X.
, and
Chen
,
H.
, “
Lattice Boltzmann model for simulating flows with multiple phases and components
,”
Phys. Rev. E
47
,
1815
(
1993
).
31.
Shan
,
X.
,
Yuan
,
X.-F.
, and
Chen
,
H.
, “
Kinetic theory representation of hydrodynamics: A way beyond the Navier–Stokes equation
,”
J. Fluid Mech.
550
,
413
441
(
2006
).
32.
Succi
,
S.
,
The Lattice Boltzmann Equation: For Fluid Dynamics and Beyond
(
Oxford University Press
,
2001
).
33.
Toubin
,
H.
, and
Bailly
,
D.
, “
Development and application of a new unsteady far-field drag decomposition method
,”
AIAA J.
53
,
3414
3429
(
2015
).
34.
Wilhelm
,
S.
,
Jacob
,
J.
, and
Sagaut
,
P.
, “
An explicit power-law-based wall model for lattice Boltzmann method-Reynolds averaged numerical simulations of their flow around airfoils
,”
Phys. Fluids
30
,
065111
(
2018
).
35.
Zhan
,
N.
,
Chen
,
R.
, and
You
,
Y.
, “
Meshfree method based on discrete gas-kinetic scheme to simulate incompressible/compressible flows
,”
Phys. Fluids
33
,
017112
(
2021
).
36.
Zhang
,
R.
,
Sun
,
C.
,
Li
,
Y.
,
Satti
,
R.
,
Shock
,
R.
,
Hoch
,
J.
, and
Chen
,
H.
, “
Lattice Boltzmann approach for local reference frames
,”
Commun. Comput. Phys.
9
,
1193
1205
(
2011
).
37.
Zhao
,
S.
,
Farag
,
G.
,
Boivin
,
P.
, and
Sagaut
,
P.
, “
Toward fully conservative hybrid lattice Boltzmann methods for compressible flows
,”
Phys. Fluids
32
,
126118
(
2020
).
38.
Zhu
,
L.
,
He
,
G.
,
Wang
,
S.
,
Miller
,
L.
,
Zhang
,
X.
,
You
,
Q.
, and
Fang
,
S.
, “
An immersed boundary method based on the lattice Boltzmann approach in three dimensions, with application
,”
Comput. Math. Appl.
61
,
3506
3518
(
2011
).