Lattice Boltzmann method for fluid-structure interaction in compressible flow

We present a two-way coupled fluid-structure interaction scheme for rigid bodies using a two-population lattice Boltzmann formulation for compressible flows. Arbitrary Lagrangian-Eulerian formulation of the discrete Boltzmann equation on body-fitted meshes is used in a combination with polynomial blending functions. The blending function approach localizes mesh deformation and allows treating multiple moving bodies with a minimal computational overhead. We validate the model with several test cases of vortex induced vibrations of single and tandem cylinders and show that it can accurately describe dynamic behavior of these systems. Finally, in the fully compressible regime, we demonstrate that the proposed model accurately captures complex phenomena such as transonic flutter over an airfoil.


I. INTRODUCTION
Interaction between fluid and deformable or moving bodies is important area of research in computational mechanics. Behavior of solid structures under the influence of fluid flow, and vice versa, is of special interest in the fields of bio-mechanics, mechanical engineering and aerospace engineering with applications ranging from blood flow inside the heart and blood vessels [1], material erosion due to bubble cavitation [2,3], effect of shock or blast waves on solid structures [4,5], and vibration or deformation of aircraft wings at high speeds [6][7][8]. For aerospace vehicle design, the study of fluid structure interactions (FSI) at transonic and supersonic speeds is especially important in order to detect shock induced vibrations and deformations such as aeroelastic flutter. A robust and efficient model to accurately predict fluid loads on structures is therefore crucial.
The lattice Boltzmann method (LBM) has steadily been gaining prominence in the domain of computational fluid dynamics (CFD). Promising results have been shown in a wide array of flows ranging from turbulence [9], multiphase flows [10] and multi-component flows [11] to rarefied gas flows [12] and relativistic hydrodynamics [13], to name a few. LBM arrives at the macroscopic equations of fluid dynamics from a mesoscopic perspective, where the flow is described by discretized particle distribution functions (populations) f i (x,t), which are associated with a set of discrete velocities C = {c i , i = 0, ..., Q − 1}, forming the links of a space-filling lattice. The LBM algorithm eventually reduces to a simple and highly efficient stream-and-collide procedure, where the populations are advected along the discrete velocities and relaxed towards local equilibrium distribution functions, which are designed to recover the Navier-Stokes equations in the hydrodynamic limit.
However, despite a lot of progress in the development of various extensions of LBM, the classical LBM is limited to the incompressible flow regime. This is mainly due to the geometrical restrictions of standard lattices, such as D2Q9 or D3Q27 in D = 2 or D = 3 dimensions with Q = 9 and Q = 27 discrete velocities, respectively, which induce errors in the fluid stress * Corresponding author bdorschn@ethz.ch tensor and break Galilean invariance. While systematically increasing the number of velocities, leading to high-order lattices [14][15][16][17][18], was thought to be a path towards compressible flows, severe restrictions due to the added computational costs and tight bounds on the temperature range remain. Another approach is to retain standard lattices but introduce appropriate correction terms to counteract the anomalous terms in the stress tensor. In the literature, different proposals for the implementation of such correction terms from a variety of authors exist [19][20][21][22][23]. Recently, we have extended this approach within the two population setting, which has shown promising results for compressible flows on standard lattices [24,25].
Another intrinsic property of LBM is that it is based on uniform Cartesian grids. While this allows for a highly efficient and exact advection scheme, non-uniform body-fitted grids can be advantageous to resolve the thin boundary layers of complex geometries. Moreover, on non-uniform meshes, a variable Courant-Friedrichs-Levy (CFL) number can be achieved easily, while it is fixed to unity for the classical LBM. The class of so-called off-lattice Boltzmann methods were developed to overcome these limitations and extend LBM to non-uniform grids. A majority of these methods turn to purely Eulerian formulations such as finite-volumes [26][27][28], finitedifferences [29,30], or finite-elements [31,32] in order to solve the discrete Boltzmann equation. While this does allow for non-uniform grids and a flexible CFL number, a partial differential equation for each discrete velocity needs to be solved in order to advect the populations. This typically requires small times steps and repeated non-local evaluation of spatial gradients, which leads to a significant computational overhead and can result in prohibitively high costs [33]. Another avenue are semi-Lagrangian methods [33][34][35][36][37], which start from a characteristic equation to resolve these issues. Promising results have been shown for wall-bounded turbulent flows [38] and even compressible flows [37,39,40].
Here, we build on these results and continue our development of semi-Lagrangian LBM for fluid-structure interaction problems of compressible flows. In the literature, there exist a variety of FSI schemes also within the realm of LBM. However, the majority are based on Cartesian meshes (see, e.g., [41,42]) using Cartesian boundary conditions or a form of the immersed boundary method (see, e.g., [43][44][45]). For these schemes, the added costs of FSI coupling is low, since even for large displacements, no re-meshing or alike is required. However, as stated above, the uniformity of the mesh limits the attainable resolution at the boundary.
In order to keep the body-fitted nature of the grid while also allowing for mesh motion, we turn to the Arbitrary Lagrangian-Eulerian (ALE) scheme. The ALE approach allows for some nodes of the mesh to move with the boundary, as in a Lagrangian description, while nodes away from the boundary stay fixed, as in an Eulerian description. Using ALE, the need to re-mesh the entire domain after displacement of the body is reduced, and the model can accommodate much larger distortions of the domain than afforded by a purely Lagrangian description, with a higher resolution near the wall than possible with a purely Eulerian description [46].
Within the ALE description, care must be taken to accommodate multiple bodies in the domain, and bodies whose motion is not a rigid body motion. In such cases, it is necessary to limit the mesh deformation due to the moving body (or bodies) to its own vicinity in order to maintain mesh conformity throughout the domain. In the literature, this problem is solved by using either of two methods. One approach is the use of polynomial blending functions that smoothly distribute the motion of the mesh spatially, such that mesh nodes close to and attached to the body move with it, while the mesh nodes that are far away from the moving body remain stationary [46][47][48]. Another approach is to solve a differential equation, such as the diffusion or linear elasticity equation in order to generate the mesh deformation field [49,50]. While the latter is more robust, solving an additional differential equation adds a significant amount of computational costs compared to solving a polynomial equation.
In our previous works [39,40], it was demonstrated that by using an arbitrary Eulerian-Lagrangian (ALE) formulation in combination with a dual population LB formulation, highspeed compressible flows with moving geometries can be captured accurately and efficiently. However, in those works, only a single body was considered and the motion was prescribed. Here, we aim at extending this framework towards genuine two-way coupled FSI, including multiple bodies as well as un-prescribed motion. To that end, we start with a rigid body FSI using an ALE formulation, where solids do not deform but are coupled with the fluid through a spring-mass system. Multiple bodies with independent displacements are accounted for by using a blended mesh formulation. Finally, we demonstrate that the proposed scheme is able to accurately capture complex phenomena such as vortex induced vibrations (VIV) and aeroelastic flutter, which is of crucial importance to aerospace applications. To the best of our knowledge there does not exist a LBM model for fluid-structure interaction in the compressible flow regime.
The paper is organized as follows: In section Sec. II we review the ALE description of LB and its extension to the semi-Lagrangian framework. We also introduce concept of blended meshes, which is crucial for the treatment of multiple moving bodies with conforming meshes. Sec. III describes the fluidstructure coupling, the force computation, and the boundary conditions necessary for the LB model. The ordinary differential equations that define the rigid body motion of the bodies are defined in Sec. IV. Validation of the coupled system is presented in Sec. V, where we consider one degree of freedom (DOF) and two DOF systems for flows past oscillating cylinders as well as airfoil flutter. Simulations over a range of Reynolds number and Mach number regimes are presented. It is shown that the LB model together with a simple yet robust mesh deformation scheme ensures mesh conformity even under highly irregular motion. Finally, conclusions are summarized in Sec. VI.

II. LATTICE BOLTZMANN MODEL
A. Arbitrary Lagrangian-Eulerian Description of LB In this section we briefly review the ALE formulation for LBM as proposed in [40]. To that end, we consider the Boltzmann equation, where the populations f i (x,t) are defined on the deformed physical domain, (x,t) and Ω i is the collision operator. In the ALE description, all computations take place on the undeformed domain (see Fig. 1). Hence, we wish to obtain a mapping G x (X,t) between the undeformed domain, (X,t 0 ) and the deformed domain, (x,t). We also introduce a unit reference frame, ξ which is fixed, and the transformation G ξ (X, t) between the undeformed domain and the unit reference domain. For now, we concern ourselves with the deformed and the undeformed domains, we will return to the unit reference domain in Sec. II B . In the simplest case of prescribed motion of a rigid body, the mapping from the deformed domain to the undeformed domain is a rigid transformation, consisting of a combination of translation and rotation. For coupled problems, where such a mapping is not known a-priori, it can be constructed at each time step from the motion of the body. Further, in case of multiple moving bodies, the mapping will not be a rigid transformation. This scenario is addressed in Sec. II D.
Assuming the existence of such a continuous mapping between the deformed and the undeformed initial domain x = G x (X,t), that maps a location X in the undeformed domain to a location x in the deformed domain, we have for the temporal derivative of this mapping, the mapping velocity V G , For the time derivative of the populations f i in the deformed domain we can write, The space derivative can similarly be transformed to the undeformed domain, Here g is the Jacobian of the ALE mapping and its inverse is given by where the mapping metrics ∂ X x, ∂ Y x, ∂ X y and ∂ Y y are the par- Finally, from Eq. (1), (3), and (4), we have the ALE description of the Boltzmann equation in the undeformed domain, where we can see that the discrete velocity c i has been replaced byĉ i , which is simply the transformation of the discrete velocity set to the initial undeformed domain, The resulting equation is in the same form as the discrete Boltzmann equation in the deformed domain, with only the discrete velocities different. These transformed discrete velocities cannot be assumed to be integers anymore and an offlattice propagation scheme is now necessary.

B. Semi-Lagrangian Framework
In this paper, we turn towards semi-Lagrangian propagation as our choice of off-lattice scheme, which was recently extended to compressible flows using a two-population LBM on standard lattice and body-fitted meshes [39]. This leads to a simpler and more efficient treatment of complex boundaries while keeping the D2Q9 lattice structure, even at high Reynolds and Mach numbers. An example of a body fitted mesh over a NACA64A010 airfoil, with refinement zones to adequately resolve shock waves, is shown in Fig. 3.
The off-lattice nature of the model implies that the departure points of the characteristic lines do not necessarily coincide with a lattice node anymore and therefore, the reconstruction of the populations at these points require interpolation. Fig. 2 describes such a scenario that results from a non-uniform grid.
In order to obtain the values of the populations f i at the offlattice departure points, we need to interpolate based on the FIG. 2. Schematic of a second-order finite-element mesh, the semi-Lagrangian advection along the discrete velocityĉ and the mapping to reference cell [39]. known values of the populations. Here, this is done using a second-order finite-element reconstruction, where the summation is carried out over all collocation points ξ s , s = 1, . . . , 9 of a reference cell with ξ = (ξ , η), (−1 ≤ ξ , η ≤ 1). N s are standard second-order Lagrangian finiteelement shape functions and ξ d p is the departure point on the reference cell. Here, we use the total Lagrangian formulation which means that all computations take place on the finiteelement reference cell, and hence, as a first step we need to compute the mapping, G ξ (X dp ,t) of the departure point in the undeformed frame, X dp = X −ĉ i δt, to the departure point in unit reference frame, ξ d p . This involves solving the following system for ξ d p : where N s are standard first-order Lagrangian finite-element shape functions. The finite-element interpolation here is computed over a first-order element, instead of second-order element in order to save computational time. Reducing the interpolation order for the departure point computation has negligible effect on the accuracy of the model (see also [39]). The computation of spatial gradients, is done using the finite-element approximation of the first derivative. For a variable C we have, where C s are the values of the variable at collocation points and J −1 is the inverse of the Jacobian of the transformation matrix computed as, and is the determinant of the Jacobian. The metrics of transformation ∂ ξ x, ∂ η x, ∂ ξ y, ∂ η y are computed with To summarize, the semi-Lagrangian propagation on moving grids can therefore be split into two steps, (i) solution of Eq. (10) to determine the departure point in the reference frame, and (ii) finite-element interpolation of the population at the departure point using Eq. (9).

C. Kinetic Equations
Here, we briefly describe the two-population model with correction terms as proposed in [24], where populations f are used to conserve mass and momentum, and populations g are responsible for energy conservation. We limit ourselves to two dimensions and the standard D2Q9 lattice with the discrete velocities c i = (c ix , c iy ), c iα ∈ {−1, 0, 1}. The kinetic equations read: where δt is the time-step. The local equilibria f i eq and g i eq , and the quasi-equilibrium g i * satisfy the local conservation laws for the macroscopic quantities: density ρ, velocity u, and total energy E, The temperature T is defined as, where C v and C p are the specific heat at constant volume and pressure, respectively, related by C v = C p − R, and R is the adiabatic gas constant, which here is set to 1. The adiabatic exponent is denoted by γ = C p /C v . The Mach number is defined as ratio of velocity to the speed of sound Ma = u/c s , where c s = √ γRT is the speed of sound. The equilibrium populations must also satisfy the Maxwell-Boltzmann (MB) relations in order to recover the full Navier-Stokes-Fourier equations in the hydrodynamic limit.
The equilibrium populations, f i eq , can be written in product form as, where, where α ∈ {x, y}. Equilibrium populations g i eq and quasiequilibrium g i * are constructed to ensure that the higher order moments of energy are reproduced [51] and read, with the weights and q eq is the equilibrium heat flux vector, where the enthalpy is defined as H = (C v + 1)/T + u 2 /2. R eq is the 4th order Maxwell-Boltzmann moment defined as, The term φ i in Eq. (15) is the correction term for the momentum equation, which is defined as (31) Another correction term, Ψ i , for the energy equation, appears in the definition of the equilibrium and quasiequilibrium populations, g i * and g i eq (Eq. (25) and (26)), which reads, Here, is the symbol for the Hadamard product or the component-wise product.
These correction terms are designed to cancel out the errors in the fluid pressure tensor that are originating from the restriction of the D2Q9 lattice. While we only present the resulting expressions here, a detailed derivation can be found in [24].
From the Chapman-Enskog analysis of the kinetic equations together with the equilibrium populations, we are able to recover the full Navier-Stokes-Fourier equations: with the heat flux q = −κ∇T and the viscous stress tensor, where the pressure is given by p = ρT . The viscosity µ and the thermal conductivity κ are related to the relaxation parameters ω and ω 1 by

D. Concept of Blended Mesh
In the case of a single body undergoing a prescribed motion, the mapping x = G x (X,t) is known a-priori. In such cases, the ALE transformation is obtained trivially by a combined translation and rotation transformation of the entire domain. The entire domain assumes a rigid motion that mirrors the motion of the moving body. For coupled motion, this mapping is not available analytically and is instead constructed at each time step from the values of the body displacement. The time derivatives in this case are then calculated using first-order finite-differences.
Here, we focus on the polynomial blended function approach for maintaining mesh conformity, for reasons explained below. In this approach, a moving body is given a region of influence, inside of which the motion of the mesh is influenced by the motion of the body. This region is divided into two parts: the rigid region and the blended region (see Fig. 4). Inside the rigid region (dark gray in Fig. 4), the mesh motion coincides with the motion of the body. This ensures mesh conformity and quality close to the walls. Inside the blended region, the influence of the moving body on mesh deformation is gradually reduced, and vanishes at the boundary of the blended region (light gray in Fig. 4). Beyond the blended region, the mesh is stationary and no special treatment of the nodes is necessary. The choice of how big to make the blended and rigid regions are dependent on the size of the body in question and the amplitude of the motion.
Inside the blended region, the distribution of mesh motion is done using blending functions. In this work, we use a 5th order polynomial function: where s R is the normalized position of the mesh node inside the blended region.
With Eq. (38), the new position of a node in the physical domain, x, can be computed from the rigid mapping, G x (X,t), and position of the node in the undeformed domain, X, as: where r blend is the blending factor, and G x (X) is the discrete map from physical to moving domain, computed at each time-step from the output of the structural solver.
In this first approach to moving bodies, we use radial distances from the center of the body for defining the region of influence. While this approach works well, it may not be optimal for slender bodies. In that case, an approach based on the distance field from each point on the moving body might be more efficient.
The blended function approach can be extended easily to multiple bodies by assigning a region of influence to each individual moving body. Given the approach as described above, it can be seen that it is not suitable for multiple moving bodies that are placed in very close proximity to each other, since the regions of influence cannot overlap.
For the simulations in this paper, we only consider motion of one or more rigid bodies. In case of multiple bodies, there is sufficient distance between them to ensure that each can be assigned an exclusive region of influence. As such we do not stand to benefit from the expensive solution of differential equations, at each time-step, to obtain the ALE map. However, for cases involving relative motion of closely placed bodies, or any arbitrary non-rigid motion of bodies, solving an additional differential equation might be a robust and attractive alternative.

III. FLUID STRUCTURE COUPLING
In this work, we use a weak coupling between solid and fluid with a staggered time-stepping approach, where fluid and solid time steps are carried out sequentially and information is exchanged in between. In particular, after a fluid time-step the LB solver computes the forces acting on the body and passes it to the structural solver, which uses the forces to compute the displacement, before passing the displacement back to the fluid solver. While this coupling is first-order accurate in time and does not enforce fluid-solid interface conditions exactly, it is computationally efficient. It is further known that such loose coupling induces artificial energy at the interface due to its staggered nature [52]. This can lead to severe instabilities for small solid-fluid density ratios, a phenomenon known as the added-mass effect. Alternatively, strong coupling approaches using sub-iterations to converge to the solid-fluid interface condition can be used to mitigate these issues. However, the added-mass effect is proportional to the time step size for compressible flows and has negligible influence for the cases considered here and no stability issues have been observed [53]. Moreover, the LB time step size is well within the stability limits that are imposed by structural solver that is used here.
The force F (x j ,t) = [F x , F y ] T acting on the jth element on the wall is given by where ds is the area of the wall element, n is the outwardpointing wall normal vector, and the fluid stress tensor is given by The non-equilibrium stress tensor can be obtained by where the non-equilibrium populations f Finally, the total force acting on the body can then be obtained via a summation of the element-wise forces over the wall, To impose wall boundary conditions in the fluid solver we follow the approach of [39,54] and replace the missing populations at the wall nodes with Grad's approximation (see where we use the target values ρ tgt , u tgt , and T tgt that are left to be specified in order to impose the boundary conditions. The equilibrium populations f eq i and g eq i are obtained from Eq. (21) and (25). The non-equilibrium population are defined using non-equilibrium values of the higher order moments. f with the non-equilibrium pressure tensor, P (1) , and the population g (1) i is given as, with and Here, S is the strain rate tensor defined as, Spatial gradients are obtained using a first order finite-element approximation of the derivative as given in Eq. (11).
To enforce no slip conditions at the wall ( point A in Fig. 5), the target values of the velocity is set to zero, u tgt = 0. The target values for density and temperature at the wall are obtained by applying a zero-gradient condition normal to the wall, For any scalar, this can be achieved by taking a first-order finite-difference approximation at a location B, which lies a unit normal distance away from the wall (see Fig. 5). In the present model with semi-Lagrangian propagation, this distance can be at maximum equal to the time-step size δt.
The value of the scalar at this point, B, is obtained by interpolation using second-order finite-elements. Using this value, together with the zero-gradient boundary conditions, we have For nodes that do not lie on the wall but lie close to the wall, and as a result also have missing populations, the target values are assigned as the value of the variable at the previous timestep.

Inlet and Outlet
At the inlet, we use equilibrium populations to impose the inlet target values directly.
At the outlet, we also impose equilibrium conditions and the target values are obtained in the same manner as for the wall boundary condition using a zero gradient condition for density, temperature as well as for velocity. Until now we have described the fluid part of the fluid structure interaction. We now describe the structural model, which takes fluid forces on the body as an input and provides the resultant displacement of the body as an output. We present two generalized models, for 1 degree of freedom and 2 degrees of freedom motion, respectively.

A. 1DOF Spring-Mass System
A single degree of freedom system, with displacement h and external driving force F y , as shown in Fig. 6, can be modeled using the Newton's second law of motion as: where m, and K h are the mass of the system and the spring constant, respectively. Here, we omit the damping term and restrict ourselves to undamped motion. The two degree-of-freedom system with coupled pitching and plunging displacements, as shown in Fig. 7, can be described using the second law of motion as before, which yields: Here, m is the total mass of the system, h and α are the plunge and pitch displacements, I α is the moment of inertia about the elastic axis (EA), and S α represents the static moment about the EA. K h and K α are the spring constants for the linear (plunge) and torsion (pitch) spring, respectively. The transverse force (lift) and the moment about the EA are represented by F y and M EA , respectively. Once again, we only consider undamped systems in this work. The airfoil half chord b is used to normalize distances, and x a is the normalized distance of the mass center from the elastic, axis which is positive if the center of mass behind the EA (towards the trailing edge) and negative if it is ahead.
The resulting system of ordinary differential equations are solved using a 4th-order Runge-Kutta method.

V. RESULTS
In the following section we present several test cases that were validated against what exist in the literature in order to show the suitability of our model to capture the dynamics and non-linearity arising from the coupling between the fluid and structural solvers. We build up from the well studied case of 1DOF vortex induced vibrations (VIV) of a single cylinder in incompressible flow, to the more complex dynamics of multiple bodies. Finally, to demonstrate a fully compressible case with two degree of freedom, we consider the case of flutter of an airfoil at transonic speeds.
A. 1DOF System: Oscillating Cylinder The motion of an elastically mounted cylinder immersed in cross-flow is widely studied and is used as a model problem for the analysis of vortex induced vibrations. The schematic of the system is shown in Fig. 6.
Here, we compare our results to those obtained by Qiu et al [55] where they used a fourth-order spectral difference method. A second-order finite-element mesh with 23,116 el- Non-dimensionalising Eq. (58) in the same way as Zhao et al [56], we have,ḧ where the lift coefficient is given by C L = 2F y /(ρU 2 ∞ D) and m * is the non-dimensional mass or the mass ratio, V R is the reduced velocity, and f N is the natural frequency of the system. h FIG. 11. 1DOF cylinder in cross-flow. Top: non-dimensional frequency of vortex shedding (red) and non-dimensional natural frequency of structural system (black). Middle: normalized amplitude of oscillation: present (red) and literature [55] (black). Bottom: amplitude of lift coefficient against: present (red) and literature (black). [55] The non-dimensional parameters are defined as: Fixing the mass ratio m * = 10, and the reduced velocity V R = 0.06Re, allows us to set the natural frequency f N of the structural system according to Eq. (63). The Reynolds number is defined as Re = ρU ∞ D/µ. The free-stream Mach number is set to Ma ∞ = U ∞ / √ γT ∞ = 0.1, T ∞ = 1/3 and γ = 1. The Strouhal number is given by St = f D/U ∞ , where f is the vortex shedding frequency. From the literature [57], it is well known that the maximum amplitudes of oscillation are characterized by the synchronisation of the vortex shedding frequency and the natural frequency of the oscillating system. This phenomenon, known as lock-in, presents itself physically as a sudden jump in the oscillation amplitude and the lift co-efficient.
Simulations were carried out at Reynolds numbers between Re = 70 and Re = 110, in the vicinity of the expected lock-in regime. Fig 11 shows the evolution of the oscillations over time, in the the lock-in regime. Snapshots of vorticity for Re = 90 (Fig. 9) and Re = 110 (Fig. 10) show the vortex shedding from the oscillating cylinders. The natural frequency of the structural system, from Eq. (63), is inversely proportional to the Reynolds number, given a fixed free-stream velocity. A clear picture of the onset of the lock-in regime can be seen in fig 11 (top). It can be observed that at approximately Re = 90, the vortex shedding frequency approaches the natural frequency of the oscillating system. A sharp jump in the amplitude of oscillation, a result of lock-in, is observed between Re = 80 and Re = 90. The corresponding jump in the maximum lift coefficient results in values that are approximately twice as large compared to the case of a stationary cylinder. The values of maximum amplitude and lift coefficient that were measured in our simulations are plotted in Fig. 11 (middle and bottom). The agreement with literature is very good, both qualitatively and quantitatively, as can be seen in Fig. 11. The simulations were run for a sufficiently long time to avoid accounting for the initial transient (see Fig. 8) . The addition of another body to the system leads to the emergence of interesting dynamic behavior, even without an additional degree of freedom. Fig. 12 represents this spring mass system, with two spherical bodies that are rigidly connected by a rod, which causes both bodies to oscillate as one. was defined in the previous section. The ratio of the body mass to the mass of the displaced fluid was set to m * = 2.

B. 1DOF System: Tandem Linked Oscillating Cylinders
Despite the fact that the motion is limited by the linkage, the effect of the downstream body, in the wake of the upstream body, reveals interesting dynamics. Zdravkovich [58] studied the flow and dynamic response of multiple cylinders in VIV and classified them into different regimes based on the dis-tance between bodies. Based on this classification, proximity interference and wake interference regimes are defined. Proximity interference occurs when the distance between the bodies is sufficiently small, while wake interference occurs when one body is submerged in the wake of another. The overlap of these regimes, proximity and wake interference, is also possible and leads to complex dynamic behavior. Here, we present test cases of VIV of two tandem cylinders at Re = 150 exhibiting all these regimes, i.e., L/D = 2, 4, 6. We also demonstrate the affect of the reduced frequency and compare our results to the work of Zhao et al [59]. In the following tests, the free-stream Mach number is set to Ma ∞ = 0.1 with γ = 1 and T = 1/3.
For L/D = 2, Zhao et. al reported beating behavior in the range of reduced velocity 4.5 ≤ V R ≤ 6. Beats are caused by interference between waves which differ by a small amount in frequency, leading to a periodic oscillation where amplitude rises slowly to its peak value and then drops rapidly before repeating again. This behaviour is well captured by our model and Fig. 13 shows the vorticity contours and the time evolution of the amplitude for V R = 5. According the classification by Zdravkovich, L/D = 2 falls under the proximity and wake interference regime with intermittent vortex shedding from the upstream cylinder. We can observe that the vortex shedding occurs for both the upstream and downstream cylinders. It is interesting to note that this phenomenon does not occur for stationary tandem cylinders at L/D = 2, i.e., vortex shedding from the upstream cylinder only occurs due to the motion of the cylinders. In the stationary case, the flow is largely stagnant in the gap between the cylinders, resulting in much lower values of the lift coefficient on the upstream cylinder. This stagnation flow in the gap also presents itself momentarily in the case of the elastic cylinders, and is responsible for the beating. In Fig. 13(6), we can see such a case when the gap is mostly stagnant, leading to a drop in lift, and causing a sharp drop in oscillation amplitude. The frequency domain is shown in Fig. 14(a) where the secondary modes causing the beats can be visualized.
At V R = 6, the natural frequency of the system is approximately equal to the vortex shedding frequency of a single stationary cylinder. The resulting synchronization leads to a regular periodic motion of the cylinders, as shown in Fig. 15 for L/D = 2, 4, and 6. Even with dual cylinders the lock-in regime centers approximately around the vortex shedding frequency of a stationary cylinder. Due to a hard lock-in, the power spectral density as shown in Fig. 14(b) exhibits a single peak, which is close to the natural frequency of the system, in this regime.
At higher values of reduced velocity, V R ≥ 9 for L/D = 6, Zhao et al [59] reported transverse oscillations with two or more dominant frequencies. Fig. 16 shows the time series of transverse displacement for L/D = 6 and V R = 9 which is able to reproduce the data reported in the literature. The frequency domain Fig. 14(c) shows two significant modes, which are larger than the natural frequency and a minor mode at the vortex shedding frequency.
The combined results of all cases, over a range of reduced velocities and distances between cylinders, are reported in Table I. It shows good quantitative agreement with what is reported in the literature while the previously shown plots are able to reproduce the dynamics reported in literature qualitatively.  [59] for two linked cylinders oscillating in cross-flow. By mounting each cylinder on separate springs, we introduce another degree of freedom, transverse to the direction of flow of the fluid. This spring-mass system with two independently oscillating bodies is shown in Fig. 17. We use a The time-step is set to δt = 1/3. For the implementation of the blended mesh, the radius of the rigid zone around each body, r rigid = 50 lattice units, and the radius of the blended zone, r blend = 400 lattice units.
Once again, there is no change to the equation of motion Eq. (60), however, we now have to solve one differential equation for each body.
The cylinders are mounted in-line with a L/D = 5, Re = 1000, and a free-stream Mach number of Ma ∞ = 0.1 with γ = 1 and T ∞ = 1/3. We simulate the response of a wide range of reduced velocities, from 4 to 10, and compare our results to the work of Jester and Kallinderis [60]. The mass ratio is set to m * = 4.
Different from the previous case of linked cylinders, the forces experienced by each body act solely on it and we therefore see distinct responses for the upstream and the downstream cylinder. According to the classification by Zdravkovich, this is also a case of wake-interference. We expect to see the upstream cylinder behave like an isolated cylinder, with a maximum amplitude peak near the range of reduced velocities that synchronize with the vortex shedding frequency. Away from this lock-in region, the amplitude of oscillation quickly drops.
The downstream cylinder, in the wake of the upstream cylinder, experiences a much larger amplitudes, over a much larger range of reduced velocities. Thus the dynamic response of the downstream cylinder is very different from the response of an isolated cylinder, while the upstream cylinder behaves much in the same way. The results agree qualitatively with previous experimental findings of [61] as well quantitatively with numerical results obtained in [60].
A visualization of the wake interference phenomenon is shown in Fig. 18, where we plot a snapshot of vorticity contours along with the displacement of both cylinders with time, for reduced velocity V R = 6. It can be observed that there is a phase shift of π/2 between the upstream and the downstream cylinder. We can further see that the gap between the two cylinders is large enough for vortex shedding to occur and the vortices, generated from the upstream cylinder, impinge on the downstream cylinder. For instance, in subfigure on the downstream cylinder, which causes the downstream cylinder to have a higher amplitude oscillations, even at much higher values of reduced velocity. In Fig. 19 we visualise the grid, deformed using blended function approach, around the cylinders at two time instances which correspond to the maximum displacement of the downstream cylinder. Animations of the mesh deformation and the vorticity can be found in the supplementary material. Fig. 20 shows the maximum amplitude of both cylinders plotted against the reduced velocity. The solid line represents the results of [60], which are in very close agreement with our results.

D. 2DOF System: NACA 64A010 Flutter Analysis
Flutter refers to an aeroelastic instability that occurs due to a coupling between aerodynamic, elastic, and inertial forces. Flutter vibrations are typically seen in slender bodies, especially for airfoils at transonic speeds due to the transient loads and shock interactions. If not controlled or actively damped, this can cause significant structural damage to the system. Below a certain threshold, the oscillations are damped out and pose no danger to the structure, but beyond, a divergent oscillating response known as flutter is observed. A neutrally stable response can also be observed, which it is known as a limit cycle oscillation (LCO).
Here, we demonstrate that our model can capture these dynamic phenomena by simulating the transonic flow over a NACA64A010 airfoil with 2 DOFs and at free-stream Mach number Ma ∞ = 0.85 with γ = 1.4 and at T ∞ = 0.2. As sketched in Fig. 7, the airfoil is mounted on a torsional and a linear spring, which allows a 2DOF motion of the airfoil.
In order to obtain the non-dimensionalized form of the 2DOF system (Eq. (59)), we introduce the nondimensionalised plunge,h = h/b, and non-dimensionalised time, τ = tω α , where ω α is the pitch frequency. The angular displacement, α is measured in radians. The nondimensionalised system, can then be written in the form: where and Here, we introduce the mass ratio m * , the reduced velocity V R and the uncoupled structural plunge frequency ω h respectively. We define the center of mass as the point where the mass of the body can be assumed to be concentrated. The radius of gyration is the distance from the mass center, where the moment of inertia of the concentrated mass equals the moment of inertia of the actual mass distribution of the body. We further denote x α as the non-dimensional distance of the center of mass from the elastic axis and r α is the radius of gyration, which is also measured from the elastic axis. The lift coefficient is defined as C L = 2F y /(ρU 2 ∞ 2b) and pitching moment coefficient is given by C M = 2M EA /(ρU 2 ∞ (2b) 2 ). These non-dimensional quantities and system parameters are defined as: These parameters are set as follows: x α = 1.8 , r α = 1.87, ω h /ω α = 1 , m * = 60, and V R = 0.5025. The elastic axis of the airfoil is located half a chord length in front of the leading edge of the airfoil. For our simulations, we use a secondorder finite-element mesh with 44,880 elements and 180,192 degrees of freedom, resulting in a near-wall resolution of 0.25 lattice units and 200 lattice units for the airfoil chord. The time-step was set to δt = 0.25/3.
Initially, a forced sinusoidal pitching perturbation with 1 • amplitude was applied for 2 full periods. Subsequently, the airfoil was allowed to oscillate freely under aerodynamic loads. The simulations are performed at Re = 1 × 10 4 , 1.25 × 10 4 , and 1.5 × 10 4 , where the Reynolds number based on the chord is defined as Re = ρU ∞ 2b/µ. We keep the Reynolds number moderate since the solver is two-dimensional and would produce nonphysical results for larger Reynolds numbers, where the flow would indeed be three-dimensional. Nevertheless, even at lower Reynolds numbers, the model can capture the complex dynamic coupling which leads to flutter. Fig. 21 shows the normalized transverse displacement and the angular displacement of the airfoil. From left to right, we can see that the airfoil exhibits damped, limit-cycle, and divergent oscillations, leading to flutter . The oscillations are damped at Re = 1 × 10 4 , exhibit LCO at Re = 1.25 × 10 4 , and flutter starts at Re = 1.5 × 10 4 . This dependence of flutter oscillations on the Reynolds number has been observed previously, both experimentally and numerically [62] [63].
It is important to note that the oscillations are of low amplitude, which is a result of keeping the Reynolds numbers low in order to maintain two-dimensional flow and warrant the use of a two-dimensional solver. This was also reported in the literature, see, e.g., [62], where the authors observed an increase of the oscillation amplitude of two orders of magnitude between Re = 1.6 × 10 4 and Re = 2.39 × 10 4 . The contours of Mach number around the airfoil are shown in Fig. 22, where we can observe the expansion of the subsonic flow over the airfoil leading to a locally supersonic region.
The two snapshots show the movement of the location of the resulting shock, as it oscillates with time.

VI. CONCLUSIONS
In this paper, we have proposed a two-way coupled fluidstructure interaction scheme for compressible flows using the semi-Lagrangian lattice Boltzmann method on unstructured meshes. Using an ALE formulation together with blending functions, the dynamics of multiple moving bodies can be described accurately. Thorough validation for vortex induced vibrations of single and tandem cylinders have shown that complex non-linear dynamics can be captured robustly and accurately, in excellent agreement with the literature. Finally, we presented a simulation of transonic flow over an airfoil that can freely pitch and plunge under aerodynamic loads. The non-linear coupling between the aerodynamic and structural forces leads to flutter and LCO regimes, which are captured by our model and agrees with what is reported in the literature. While we focused on two-dimensional flows in this paper, the extension to three dimensions is conceptually straightforward and will be focus of future works.