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.

## I. INTRODUCTION

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).

## II. THE LATTICE BOLTZMANN METHOD USING AN OVERSET TECHNIQUE AND A HYBRID RECURSIVE REGULARIZED COLLISION MODEL

### A. The hybrid recursive regularized collision model

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

where $\Omega i(x,t)$ is the collision operator and $\Delta 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

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

and

where *w _{i}* and

*c*are, respectively, the lattice weights and the speed of sound of the lattice, the $Hi(n)$ are the discrete Hermite polynomials of order

_{s}*n*, and $a0(n)$ and $a1(n)$ are the Hermite coefficients of the equilibrium and off-equilibrium distribution, respectively,

where *q* is the number of discrete lattice velocities.

In addition, $\Pi (1)$ is the deviatoric stress tensor, defined as

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,

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

with

and the off-equilibrium Hermite coefficients write

### B. The lattice Boltzmann overset method

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).

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

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

where $\omega $ 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)

Thus, $uc$ is equal to the velocity ** u** shifted by $F/(2\rho )$, 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

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.

### C. The border condition problem in the lattice Boltzmann overset framework

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 *f _{i}* 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?

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.

### D. The original overset-HRR method

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**, $\Pi (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*

where $\varphi $ 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:

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\xd7\omega $, where $\omega $ 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

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

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:

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

where *S _{i}* is the source term standing for the Coriolis and centrifugal forces related to the non-inertial rotating frame, namely,

Also, the $\u0302$ symbol means that $f\u0302if$ 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\xd7\omega $ 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)(28)$\u2200k\u2208[1..8],\u2003uPm(xkm)=um(xkm)+vr(xkm),$
- the fixed frame being the inertial reference frame, there is no source term
*S*and_{i}(29)$fif=fi,regf=g(0)(\rho f,uf)+g(1)(\rho f,uf,\Pi f,(1))=fi\u0302m.$

## III. A NEW LATTICE BOLTZMANN DIRECT COUPLING OVERSET METHOD FOR THE MOVING SUBDOMAIN BOUNDARY PROBLEM

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.

### A. Rotation and interpolation of distribution functions in the context of the overset technique

#### 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

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

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

where $\varphi $ is the bilinear interpolation operator and the *I* index means that it is an interpolated value.

The macroscopic variables can then be reconstructed as

It should be noted that the deviatoric stress tensor is reconstructed by interpolation and rotation, exactly as in the original method, that is, $\Pi m,(1)(xm)=\Pi 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

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 *S _{i}* 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

Also, the $\u0302$ symbol means that $f\u0302if$ 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\xd7\omega $ 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)(38)$fi;Pm,(0)=g(0)(xkm;\rho km,ukm+vr),$(39)$fi;Pm,(1)=g(1)(xkm;\rho km,ukm),$
- the fixed frame being the inertial reference frame, there is no source term
*S*and_{i}(40)$fif=fi;P,R,If,(0)+fi;P,R,If,(1)=f\u0302im.$

### B. Description of the direct coupling overset method

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

where $fi\u2009fine,(1)$ and $fi\u2009coarse,(1)$ are the distribution functions at the interface related, respectively, to the fine and coarse grids, and $R=0.5\u2009\tau fine/\tau 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:

where $\Pi 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)

where $\varphi $ is the bilinear interpolation operator and the *I* index means that it is an interpolated value. In Eq. (43), $\Pi f,(1)\u2009(xkf)$ is defined as

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

where $\Phi 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)

Using the relation $fi(1)=fi\u2212fi(0)(X)$, where $X=(\rho ,ux,uy,uz)$ is the vector of relevant macroscopic variables, it yields

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\u0302if$ of Subsection III A. This reads

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\u0302if$ are known. The only unknown is ** X**. Thus, solving Eq. (47) allows to obtain all macroscopic variables $X=(\rho ,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

### C. Resolution of the non-linear system

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

where $G(X)=\u2211i\Phi i\u2009(f\u0302if\u2212fim,(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 $\delta X$ a small variation around the root $X0$. Then, the Newton–Raphson method provides a good approximation for the root as

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

This algorithm can then be applied iteratively to obtain

The convergence criteria are set at the user's convenience; here, we consider the method as converged when $||\delta X||<10\u221210$.

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.

## IV. NUMERICAL VALIDATION OF THE NEW DIRECT COUPLING OVERSET-HRR ALGORITHM

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.

### A. Aerodynamic test cases

#### 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 $\rho 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 ($\rho outlet=\rho ref$). The Reynolds number is *Re *=* *90. The computational domain is a box of dimensions $Lx=100\xd7Ly=100$. It contains a fixed grid of the same dimensions and a moving grid of dimensions $Lxm=14\xd7Lym=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 $\omega 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.

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

and

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 $||\rho ||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.

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\rho ||1=2.86\xb710\u22127$ % for the density and $||Eux||1=6.54\xb710\u22126$ % 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 ($|\u2009||\rho ||1Original\u2212||\rho ||1DC\u2009|tfinal=|\u2009||ux||1Original\u2212||ux||1DC\u2009|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

being a constant in $\mathbb{R}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)].*X*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 $10\u221218$ 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\xd710\u22126$ % for density and $1.490\xd710\u22124$ % for velocity.

x = 15
. | Original method . | DC method . |
---|---|---|

$Mean\u2009(||E\rho ||1)$ (%) | $1.668\u2009800\u2009001\u2009297\u2009654\xd710\u22126$ | $1.668\u2009800\u2009001\u2009297\u2009654\xd710\u22126$ |

$Max\u2009(||E\rho ||1)$ (%) | $5.960\u2009000\u2009004\u2009634\u2009480\xd710\u22126$ | $5.960\u2009000\u2009004\u2009634\u2009480\xd710\u22126$ |

$Mean\u2009(||Eux||1)$ (%) | $1.362\u2009759\u2009979\u2009612\u2009352\xd710\u22125$ | $1.362\u2009759\u2009979\u2009612\u2009352\xd710\u22125$ |

$Max\u2009(||Eux||1)$ (%) | $2.235\u2009999\u2009966\u2009892\u2009892\xd710\u22125$ | $2.235\u2009999\u2009966\u2009892\u2009892\xd710\u22125$ |

x = 15
. | Original method . | DC method . |
---|---|---|

$Mean\u2009(||E\rho ||1)$ (%) | $1.668\u2009800\u2009001\u2009297\u2009654\xd710\u22126$ | $1.668\u2009800\u2009001\u2009297\u2009654\xd710\u22126$ |

$Max\u2009(||E\rho ||1)$ (%) | $5.960\u2009000\u2009004\u2009634\u2009480\xd710\u22126$ | $5.960\u2009000\u2009004\u2009634\u2009480\xd710\u22126$ |

$Mean\u2009(||Eux||1)$ (%) | $1.362\u2009759\u2009979\u2009612\u2009352\xd710\u22125$ | $1.362\u2009759\u2009979\u2009612\u2009352\xd710\u22125$ |

$Max\u2009(||Eux||1)$ (%) | $2.235\u2009999\u2009966\u2009892\u2009892\xd710\u22125$ | $2.235\u2009999\u2009966\u2009892\u2009892\xd710\u22125$ |

x = 85
. | Original method . | DC method . |
---|---|---|

$Mean\u2009(||E\rho ||1)$ (%) | $1.311\u2009200\u2009001\u2009019\u2009586\xd710\u22126$ | $1.311\u2009200\u2009001\u2009019\u2009586\xd710\u22126$ |

$Max\u2009(||E\rho ||1)$ (%) | $5.960\u2009000\u2009004\u2009634\u2009480\xd710\u22126$ | $5.960\u2009000\u2009004\u2009634\u2009480\xd710\u22126$ |

$Mean\u2009(||Eux||1)$ (%) | $4.357\u2009639\u2009934\u2009587\u2009714\xd710\u22125$ | $4.357\u2009639\u2009934\u2009587\u2009714\xd710\u22125$ |

$Max\u2009(||Eux||1)$ (%) | $1.489\u2009999\u2009977\u2009559\u2009619\xd710\u22124$ | $1.489\u2009999\u2009977\u2009559\u2009619\xd710\u22124$ |

x = 85
. | Original method . | DC method . |
---|---|---|

$Mean\u2009(||E\rho ||1)$ (%) | $1.311\u2009200\u2009001\u2009019\u2009586\xd710\u22126$ | $1.311\u2009200\u2009001\u2009019\u2009586\xd710\u22126$ |

$Max\u2009(||E\rho ||1)$ (%) | $5.960\u2009000\u2009004\u2009634\u2009480\xd710\u22126$ | $5.960\u2009000\u2009004\u2009634\u2009480\xd710\u22126$ |

$Mean\u2009(||Eux||1)$ (%) | $4.357\u2009639\u2009934\u2009587\u2009714\xd710\u22125$ | $4.357\u2009639\u2009934\u2009587\u2009714\xd710\u22125$ |

$Max\u2009(||Eux||1)$ (%) | $1.489\u2009999\u2009977\u2009559\u2009619\xd710\u22124$ | $1.489\u2009999\u2009977\u2009559\u2009619\xd710\u22124$ |

#### 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\xd7Ly=100$ and $Lxm=14\xd7Lym=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 $\omega 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):

where $ux,max=0.1$.

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

The temporal evolutions of the spatial L1 norm of the density field $||\rho ||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, $|\u2009||\rho ||1Original\u2212||\rho ||1DC\u2009|tfinal=3.216\xd710\u221210$ and $|\u2009||ux||1Original\u2212||ux||1DC\u2009|tfinal=4.894\xd710\u221211$. In the end, the following error values are found when compared to the reference computation: $||E\rho ||1=7.89\xd710\u22125$ % for the density and $||Eux||1=7.10\xd710\u22125$ % 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.

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\xd710\u22124$ % for density and $2.430\xd710\u22122$ % for velocity.

x = 15
. | Original method . | DC method . |
---|---|---|

$Mean\u2009(||E\rho ||1)$ (%) | $1.176\u2009583\u2009906\u2009798\u2009957\xd710\u22124$ | $1.174\u2009208\u2009598\u2009866\u2009126\xd710\u22124$ |

$Max\u2009(||E\rho ||1)$ (%) | $1.546\u2009827\u2009095\u2009575\u2009228\xd710\u22124$ | $1.428\u2009284\u2009107\u2009822\u2009683\xd710\u22124$ |

$Mean\u2009(||Eux||1)$ (%) | $7.881\u2009696\u2009528\u2009406\u2009885\xd710\u22124$ | $7.880\u2009878\u2009689\u2009514\u2009660\xd710\u22124$ |

$Max\u2009(||Eux||1)$ (%) | $3.638\u2009914\u2009170\u2009708\xd710\u22123$ | $3.638\u2009914\u2009170\u2009708\xd710\u22123$ |

x = 15
. | Original method . | DC method . |
---|---|---|

$Mean\u2009(||E\rho ||1)$ (%) | $1.176\u2009583\u2009906\u2009798\u2009957\xd710\u22124$ | $1.174\u2009208\u2009598\u2009866\u2009126\xd710\u22124$ |

$Max\u2009(||E\rho ||1)$ (%) | $1.546\u2009827\u2009095\u2009575\u2009228\xd710\u22124$ | $1.428\u2009284\u2009107\u2009822\u2009683\xd710\u22124$ |

$Mean\u2009(||Eux||1)$ (%) | $7.881\u2009696\u2009528\u2009406\u2009885\xd710\u22124$ | $7.880\u2009878\u2009689\u2009514\u2009660\xd710\u22124$ |

$Max\u2009(||Eux||1)$ (%) | $3.638\u2009914\u2009170\u2009708\xd710\u22123$ | $3.638\u2009914\u2009170\u2009708\xd710\u22123$ |

x = 85
. | Original method . | DC method . |
---|---|---|

$Mean\u2009(||E\rho ||1)$ (%) | $2.955\u2009730\u2009755\u2009523\u2009385\xd710\u22125$ | $2.955\u2009730\u2009755\u2009523\u2009385\xd710\u22125$ |

$Max\u2009(||E\rho ||1)$ (%) | $7.157\u2009658\u2009657\u2009838\u2009396\xd710\u22125$ | $7.157\u2009658\u2009657\u2009838\u2009396\xd710\u22125$ |

$Mean\u2009(||Eux||1)$ (%) | $1.364\u2009851\u2009243\u2009599\u20090\xd710\u22122$ | $1.364\u2009823\u2009061\u2009568\u20091\xd710\u22122$ |

$Max\u2009(||Eux||1)$ (%) | $2.430\u2009288\u2009413\u2009781\u20090\xd710\u22122$ | $2.430\u2009288\u2009413\u2009781\u20090\xd710\u22122$ |

x = 85
. | Original method . | DC method . |
---|---|---|

$Mean\u2009(||E\rho ||1)$ (%) | $2.955\u2009730\u2009755\u2009523\u2009385\xd710\u22125$ | $2.955\u2009730\u2009755\u2009523\u2009385\xd710\u22125$ |

$Max\u2009(||E\rho ||1)$ (%) | $7.157\u2009658\u2009657\u2009838\u2009396\xd710\u22125$ | $7.157\u2009658\u2009657\u2009838\u2009396\xd710\u22125$ |

$Mean\u2009(||Eux||1)$ (%) | $1.364\u2009851\u2009243\u2009599\u20090\xd710\u22122$ | $1.364\u2009823\u2009061\u2009568\u20091\xd710\u22122$ |

$Max\u2009(||Eux||1)$ (%) | $2.430\u2009288\u2009413\u2009781\u20090\xd710\u22122$ | $2.430\u2009288\u2009413\u2009781\u20090\xd710\u22122$ |

#### 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\xd7Ly=798$. It includes a circular moving grid of radius 30. The moving grid rotates at the angular velocity $\omega 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\u221e=0.017,0.0)$, while at the outlet $\rho ref=1$. The viscosity is imposed as $\nu =1.70\xd710\u22123$, 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\omega z/U\u221e$), 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.

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 *C _{p}* 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).

. | Original method . | DC method . | Mittal and Kumar (2003) . |
---|---|---|---|

C _{L} | 1.168 043 262 926 71 | 1.168 034 110 111 57 | 1.068 |

C _{D} | 2.476 433 333 358 73 | 2.476 384 452 553 85 | 2.4 |

. | Original method . | DC method . | Mittal and Kumar (2003) . |
---|---|---|---|

C _{L} | 1.168 043 262 926 71 | 1.168 034 110 111 57 | 1.068 |

C _{D} | 2.476 433 333 358 73 | 2.476 384 452 553 85 | 2.4 |

Figure 12 shows the density and velocity vertical profiles upstream ($x/D=\u221210$) 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\xd710\u22123$ % 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.

$x/D=\u221210$ . | Original method . | DC method . |
---|---|---|

$Mean\u2009(||E\rho ||1)$ (%) | $2.980\u2009363\u2009525\u2009742\u2009923\u20096\xd710\u22124$ | $2.980\u2009363\u2009525\u2009742\u2009923\u20096\xd710\u22124$ |

$Max\u2009(||E\rho ||1)$ (%) | $8.333\u2009402\u2009778\u2009411\u2009081\xd710\u22124$ | $8.333\u2009402\u2009778\u2009411\u2009081\xd710\u22124$ |

$Mean\u2009(||E||u||||1)$ (%) | $3.799\u2009049\u2009107\u2009179\u2009938\xd710\u22122$ | $3.799\u2009106\u2009136\u2009119\u2009404\u20096\xd710\u22122$ |

$Max\u2009(||E||u||||1)$ (%) | $6.264\u2009432\u2009073\u2009432\u2009659\xd710\u22122$ | 6.264 431 322 722 924 $\xd710\u22122$ |

$x/D=\u221210$ . | Original method . | DC method . |
---|---|---|

$Mean\u2009(||E\rho ||1)$ (%) | $2.980\u2009363\u2009525\u2009742\u2009923\u20096\xd710\u22124$ | $2.980\u2009363\u2009525\u2009742\u2009923\u20096\xd710\u22124$ |

$Max\u2009(||E\rho ||1)$ (%) | $8.333\u2009402\u2009778\u2009411\u2009081\xd710\u22124$ | $8.333\u2009402\u2009778\u2009411\u2009081\xd710\u22124$ |

$Mean\u2009(||E||u||||1)$ (%) | $3.799\u2009049\u2009107\u2009179\u2009938\xd710\u22122$ | $3.799\u2009106\u2009136\u2009119\u2009404\u20096\xd710\u22122$ |

$Max\u2009(||E||u||||1)$ (%) | $6.264\u2009432\u2009073\u2009432\u2009659\xd710\u22122$ | 6.264 431 322 722 924 $\xd710\u22122$ |

$x/D=10$ . | Original method . | DC method . |
---|---|---|

$Mean\u2009(||E\rho ||1)$ (%) | $5.436\u2009273\u2009590\u2009404\u2009883\xd710\u22124$ | $5.436\u2009273\u2009590\u2009404\u2009883\xd710\u22124$ |

$Max\u2009(||E\rho ||1)$ (%) | $4.166\u2009805\u2009560\u2009179\u2009329\xd710\u22123$ | $4.166\u2009805\u2009560\u2009179\u2009329\xd710\u22123$ |

$Mean\u2009(||E||u||||1)$ (%) | $5.677\u2009881\u2009987\u2009115\u200941\xd710\u22121$ | $5.679\u2009093\u2009365\u2009767\u2009347\xd710\u22121$ |

$Max\u2009(||E||u||||1)$ (%) | 6.812 550 348 059 371 | 6.813 554 108 265 357 |

$x/D=10$ . | Original method . | DC method . |
---|---|---|

$Mean\u2009(||E\rho ||1)$ (%) | $5.436\u2009273\u2009590\u2009404\u2009883\xd710\u22124$ | $5.436\u2009273\u2009590\u2009404\u2009883\xd710\u22124$ |

$Max\u2009(||E\rho ||1)$ (%) | $4.166\u2009805\u2009560\u2009179\u2009329\xd710\u22123$ | $4.166\u2009805\u2009560\u2009179\u2009329\xd710\u22123$ |

$Mean\u2009(||E||u||||1)$ (%) | $5.677\u2009881\u2009987\u2009115\u200941\xd710\u22121$ | $5.679\u2009093\u2009365\u2009767\u2009347\xd710\u22121$ |

$Max\u2009(||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: $|\u2009||\rho ||1Original\u2212||\rho ||1DC\u2009|tfinal=1.141\xd710\u22129$ and $|\u2009||u||1Original\u2212||u||1DC\u2009|=4.852\xd710\u22126$.

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.

### B. Aeroacoustic test cases

#### 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\xd7Ly=100$. It contains a fixed grid of the same dimensions and a moving grid of dimensions $Lxm=39\xd7Lym=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 $\omega z=0.000\u200959$, and the kinematic viscosity is $\nu =1\xd710\u22124$. The acoustic pulse is initialized at the center of the box as follows:

where *A *=* *0.001, $xc=yc=50$, and *R _{c}* = 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.

Figure 14 shows the horizontal profiles located at the center of the domain (*y *=* *50) of density, both components of velocity *u _{x}* and

*u*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.

_{y}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\xb710\u22126$ for density and $6.367\xb710\u22127$ 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.

. | Original method . | DC method . |
---|---|---|

$Mean\u2009|\u2009\Vert \rho \Vert 1\u2212\Vert \rho \Vert 1REF|$ | $8.126\u2009333\u2009333\u2009202\u2009674\xd710\u22127$ | $8.126\u2009333\u2009333\u2009202\u2009674\xd710\u22127$ |

$Max\u2009| \Vert \rho \Vert 1\u2212\Vert \rho \Vert 1REF|$ | $1.192\u2009500\u2009000\u2009080\u2009088\xd710\u22126$ | $1.192\u2009500\u2009000\u2009080\u2009088\xd710\u22126$ |

$Mean\u2009(|\u2009\Vert u\Vert 1\u2212||u||1REF\u2009|)$ | $3.783\u2009709\u2009888\u2009991\u2009750\xd710\u22127$ | $3.783\u2009664\u2009374\u2009708\u2009644\xd710\u22127$ |

$Max\u2009(|\u2009\Vert u\Vert 1\u2212||u||1REF\u2009|)$ | $6.367\u2009026\u2009105\u2009555\u2009673\xd710\u22127$ | $6.366\u2009852\u2009250\u2009907\u2009429\xd710\u22127$ |

. | Original method . | DC method . |
---|---|---|

$Mean\u2009|\u2009\Vert \rho \Vert 1\u2212\Vert \rho \Vert 1REF|$ | $8.126\u2009333\u2009333\u2009202\u2009674\xd710\u22127$ | $8.126\u2009333\u2009333\u2009202\u2009674\xd710\u22127$ |

$Max\u2009| \Vert \rho \Vert 1\u2212\Vert \rho \Vert 1REF|$ | $1.192\u2009500\u2009000\u2009080\u2009088\xd710\u22126$ | $1.192\u2009500\u2009000\u2009080\u2009088\xd710\u22126$ |

$Mean\u2009(|\u2009\Vert u\Vert 1\u2212||u||1REF\u2009|)$ | $3.783\u2009709\u2009888\u2009991\u2009750\xd710\u22127$ | $3.783\u2009664\u2009374\u2009708\u2009644\xd710\u22127$ |

$Max\u2009(|\u2009\Vert u\Vert 1\u2212||u||1REF\u2009|)$ | $6.367\u2009026\u2009105\u2009555\u2009673\xd710\u22127$ | $6.366\u2009852\u2009250\u2009907\u2009429\xd710\u22127$ |

#### 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\xd7Ly=200$. It contains a fixed grid of the same dimensions and a moving grid of dimensions $Lxm=79\xd7Lym=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 $\omega z=0.006$, and the kinematic viscosity is $\nu =1\xb710\u221210$. The convected vortex is initialized at the center of the box as follows ($Re=4\xb71010$, *Ma *=* *0.17):

where $\rho 0=1,\u2009Ux=0.1,\u2009\u03f5=0.15,\u2009xc=yc=100$, and *R _{c}* = 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.

Figure 16 shows the horizontal profiles located at the center of the domain (*y *=* *50) of density, both components of velocity *u _{x}* and

*u*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

_{y}*t*=

*0, $t=10\u2009tadv$, and $t=20\u2009tadv$, where*

*t*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.

_{adv}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\xd7,\u20096.18\xd7,\u20092.04\xd7$ 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\xd7$ and $2.11\xd7$ lower using the direct coupling approach over the original method, for the left one (*x *=* *25) and the right one (*x *=* *76), respectively.

. | Original method . | DC method . |
---|---|---|

$Mean\u2009|\u2009\Vert \rho \Vert 1\u2212\Vert \rho \Vert 1REF|$ | $9.333\u2009193\u2009049\u2009999\u2009910\xd710\u22124$ | $4.074\u2009471\u2009959\u2009999\u2009959\xd710\u22124$ |

$Max\u2009| \Vert \rho \Vert 1\u2212\Vert \rho \Vert 1REF|$ | $4.724\u2009383\u2009400\u2009000\xd710\u22123$ | $2.678\u2009275\u2009100\u2009000\xd710\u22123$ |

$Mean\u2009(|\u2009\Vert ux\Vert 1\u2212||ux||1REF\u2009|)$ | $1.530\u2009553\u2009475\u2009055\xd710\u22123$ | $3.368\u2009238\u2009869\u2009740\u2009786\xd710\u22124$ |

$Max\u2009(|\u2009\Vert ux\Vert 1\u2212||ux||1REF\u2009|)$ | $4.819\u2009702\u2009046\u2009274\xd710\u22123$ | $7.792\u2009670\u2009588\u2009818\u2009718\xd710\u22124$ |

$Mean\u2009(|\u2009\Vert uy\Vert 1\u2212||uy||1REF\u2009|)$ | $4.765\u2009100\u2009393\u2009857\xd710\u22123$ | $2.212\u2009539\u2009263\u2009029\xd710\u22123$ |

$Max\u2009(|\u2009\Vert uy\Vert 1\u2212||uy||1REF\u2009|)$ | $1.927\u2009325\u2009363\u2009924\u20092\xd710\u22122$ | $9.425\u2009645\u2009768\u2009786\xd710\u22123$ |

. | Original method . | DC method . |
---|---|---|

$Mean\u2009|\u2009\Vert \rho \Vert 1\u2212\Vert \rho \Vert 1REF|$ | $9.333\u2009193\u2009049\u2009999\u2009910\xd710\u22124$ | $4.074\u2009471\u2009959\u2009999\u2009959\xd710\u22124$ |

$Max\u2009| \Vert \rho \Vert 1\u2212\Vert \rho \Vert 1REF|$ | $4.724\u2009383\u2009400\u2009000\xd710\u22123$ | $2.678\u2009275\u2009100\u2009000\xd710\u22123$ |

$Mean\u2009(|\u2009\Vert ux\Vert 1\u2212||ux||1REF\u2009|)$ | $1.530\u2009553\u2009475\u2009055\xd710\u22123$ | $3.368\u2009238\u2009869\u2009740\u2009786\xd710\u22124$ |

$Max\u2009(|\u2009\Vert ux\Vert 1\u2212||ux||1REF\u2009|)$ | $4.819\u2009702\u2009046\u2009274\xd710\u22123$ | $7.792\u2009670\u2009588\u2009818\u2009718\xd710\u22124$ |

$Mean\u2009(|\u2009\Vert uy\Vert 1\u2212||uy||1REF\u2009|)$ | $4.765\u2009100\u2009393\u2009857\xd710\u22123$ | $2.212\u2009539\u2009263\u2009029\xd710\u22123$ |

$Max\u2009(|\u2009\Vert uy\Vert 1\u2212||uy||1REF\u2009|)$ | $1.927\u2009325\u2009363\u2009924\u20092\xd710\u22122$ | $9.425\u2009645\u2009768\u2009786\xd710\u22123$ |

### C. Sensitivity study to the moving grid local Mach number and convergence rate

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 $\omega $. The moving grid maximum local Mach number $MaMG,max$ is defined as follows:

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 *l _{MG}* is the side length of the square moving mesh. For this sensitivity study, the following values for $\omega $ have been chosen: $\omega 1,\u2004\omega 2=2\omega 1,\omega 3=3\omega 1,\u2004\omega 4=4\omega 1,\u2004\omega 5=5\omega 1,\u2004\omega 6=6\omega 1$. The corresponding $MaMG,max$ values are displayed in Table VII.

. | $MaMG,max,1$ . | $MaMG,max,2$ . | $MaMG,max,3$ . | $MaMG,max,4$ . | $MaMG,max,5$ . | $MaMG,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,1$ . | $MaMG,max,2$ . | $MaMG,max,3$ . | $MaMG,max,4$ . | $MaMG,max,5$ . | $MaMG,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 %.

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 $\omega $ 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

The L2 and infinity norms of both density and velocity for three different mesh resolutions (number of cells $N1=9559,\u2009N2=4\u2009N1,\u2009N3=16\u2009N1$) 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.

## V. CONCLUSION

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.

## ACKNOWLEDGMENTS

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.

## DATA AVAILABILITY

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