A numerical framework for simulating fluid-structure interaction phenomena

In this paper, a numerical tool able to solve fluid-structure interaction problems is proposed. The lattice Boltzmann method is used to compute fluid dynamics, while the corotational finite element formulation together with the Time Discontinuous Galerkin method are adopted to predict structure dynamics. The Immersed Boundary method is used to account for the presence of an immersed solid in the lattice fluid background and to handle fluid-structure interface conditions, while a Volume-of-Fluid-based method is adopted to take trace of the evolution of the free surface. These ingredients are combined through a partitioned staggered explicit strategy, according to an efficient and accurate algorithm recently developed by the authors. The effectiveness of the proposed methodology is tested against two different cases. The former investigates the dam break phenomenon, involving the modeling of the free surface. The latter involves the vibration regime experienced by two highly deformable flapping flags obstructing a flow. A wide numerical campaign is carried out by computing the error in terms of interface energy artificially introduced at the fluid-solid interface. Moreover, the structure behavior is dissected by simulating scenarios characterized by different values of the Reynolds number. Present findings are compared to literature results, showing a very close agreement.


INTRODUCTION
omputational fluid dynamics (CFD) represents a set of scientific methods whose aim is to solve problems, which involve fluids by adopting computers, algorithms and numerical methods.CFD can be studied and analyzed at three different levels.The "top" one is the so-called macroscopic level, consisting of the solution of the Navier-Stokes equations, which govern flow behavior.The "bottom" one reduces a fluid flow problem to a microscopic one and it is solved by adopting laws arising from molecular dynamics.Between these approaches, in the last decades the lattice Boltzmann (LB) method [1,2] developed as a powerful alternative to standard approaches.It is characterized by a C mesoscopic point of view, where the flow behavior emerges by the space-time evolution of a function describing the density of a set of fictitious particles.These particles move upon an Eulerian grid which is kept fixed during the overall analysis.Among the possible advantages of the LB method over the macroscopic-based one, the most immediate is the greater computational efficiency, especially in interaction problems involving moving structures.It has been demonstrated that the solution given by the LB method is equivalent to the one of the Navier-Stokes equations for an incompressible flow with a second order of accuracy [3].In order to account for the presence of a solid body immersed in the lattice fluid background, the Immersed Boundary (IB) method is adopted [4,5].With respect to the well-consolidated bounce-back rule [6,7], such method is characterized by superior properties in terms of stability.According to a recent work carried out by the authors, it preserves a high level of accuracy and it leads to a quite general implementation of the numerical algorithm [8].A Volume-Of-Fluid-based approach (VOF) [9,10] is used to take trace of the liquid-gas interface movements due to the mass advection.In particular, the model uses the VOF to update the nodes fill levels (in terms of mass) and a pressure boundary to restore the missing distribution functions at the liquid-gas interface.When moving structures are considered, the non-linear structure dynamics is computed by adopting the Time Discontinuous Galerkin (TDG) scheme.This choice over standard Newmark algorithms is motivated by its superior properties in terms of stability, accuracy and convergence [11][12][13].In this paper, the LB, VOF, IB and TDG methods are combined within a staggered explicit coupling algorithm called FELBA (Finite Element Lattice Boltzmann Analysis), which has been validated by the authors for different applications, including mechanics [14,15], industry [16], flapping flight [17,18], and even shallow waters [19], among the others.In particular, the above ingredients are combined together in order to predict the flow physics arising in two different cases: a dam break and two solid structures obstructing the development of a flow in a channel.Structures are idealized by Euler-Bernoulli finite elements capable of undergoing large displacements, according to the corotational formulation [20,21].The accuracy of the coupling algorithm is evaluated in terms of the energy that is artificially introduced in the system at the fluid-structure interface.Here, scenarios characterized by different values of the Reynolds number are investigated.Findings in terms of the three components of the tip displacement are given for the above-mentioned configurations, together with considerations about the artificially introduced interface energy.As regards the simulation of the dam break phenomenon, findings from a numerical analysis are compared to literature results, showing the effectiveness of the proposed approach.

PROBLEM STATEMENT
he flow physics is computed by solving the LB equation, that is fi(x+ci t, t+t) = fi(x,t) + [fi(x,t)-fi eq (x,t)] with i=0,…,8, (1) where f i are the particle distribution functions, x is the position, t is the time,  is the relaxation frequency, t is the time step and ci are the lattice vectors in the adopted D2Q9 model [2].Notice that the relaxation frequency is strictly related to the fluid viscosity, i.e. =(1/-0.5)/3.The equilibrium distribution functions f i eq are computed as a second-order expansion in the local Mach number [2].The presence of an immersed body is accounted for by adding at the right-hand side of Eq. ( 1) the quantity 3c i w i g i (x,t), where g i (x,t) is a term computed via an implicit velocity-correction based IB method [22].Once the corrected Eq. ( 1) is solved, the macroscopic density  and fluid velocity v are computed as On the other hand, deformable solids are idealized as corotational beam finite elements.The resultant non-linear equation of the solid motion is integrated in time by adopting the TDG scheme according to the implementation proposed in [13].As discussed in [16], the TDG scheme has shown superior properties with respect to the adoption of standard Newmark algorithms for solving fluid-structure interaction problems.As mentioned above, in order to handle fluid-structure interface condition, the IB method is adopted.In this way, the solid body is represented by a set of Lagrangian points, independent from the Eulerian fluid mesh, thus allowing a quite general implementation of the algorithm of computation.If a free surface flow is considered, the difference in terms of density and viscosity between the two phases composing the system is usually very high.Thus, it is possible to affirm that the behavior of the overall system is governed by the phase with the higher density [9].The phase characterized by the lower density is considered only as a surface boundary condition, thus it is neglected in the rest of the domain.The adopted free surface model represents a combination of the LB and VOF methods, the latter being responsible of the interface tracking.This method requires a particular treatment T of the interface lattice sites.Since the part composed of the lower density is not affected by the solution of the LB equation, the missing values of the particle distribution functions should be computed by a pressure boundary, as sketched in Fig. 1.The interested reader can refer to [9,10] for further discussion about this approach.The present method, as a whole, is briefly described in the following steps: 1. initialize fluid and solid domain; 2. solve the fluid system, Eq. (1) (LB); 3. compute mass fluxes (VOF) and reconstruct a valid set of distribution functions at the interface (pressure boundary); 4. compute the fluid macroscopic variables and v, Eq. ( 2); 5. enforce the no-slip condition at the fluid-structure interface via the iterative IB algorithm (compute g i (x,t) and compute the forces on the structure); 6. correct the populations with the term g i (x,t) and compute, again, the fluid macroscopic variables and v, Eq. (2); 7. perform structure solution; 8. update nodes fill levels and re-initialize lattice state (i.e.emptied interface node become gas node, etc.); 9. advance in time going to step 2.

RESULTS AND DISCUSSION
n this section, two cases are investigated.First, the dam break phenomenon is simulated.Then, the dynamics of two slender deformable beams invested by a flow is examined.

Dam break
Making reference to Fig. 2, a water column of height H=0.3m and width B=0.6m is placed in the right part of a tank of length L=1.61m, height D=6m.At t=0 the wall placed at x=0 is removed.The system described has been discretized using a relatively coarse grid, composed of 400 lattice nodes in x-direction and 150 in z-direction.During the test, the height of the water column is recorded at 4 checkpoints, H1, H2, H3 and H4.In the following, findings from a LB simulation are compared to the experimental data given in [23][24][25].In Fig. 3, the time evolution of the simulated and experimental free-surface profiles are compared.Moreover, in Tab. 1, the simulated and experimental time instants at which the water column crosses the checkpoints H2, H3 and H4 and impacts the boundary are reported.As it is possible to observe, a very close agreement is experienced.Finally, in Fig. 4, the present findings in terms of water height are compared with literature results [23][24][25] at the above mentioned checkpoints.As it can be noted, also in this case a good agreement is obtained.

Two slender beams obstructing a flow
Making reference to Fig. 5, two identical beams (blue lines) of length L=25 are placed in a channel of width W=200 and height H=60.As it is possible to observe, as Re grows, the amplitude of the displacement tends to reduce.Such behavior has to be addressed to the fact that the viscous-dependent shear component of the stress induced by the flow upon the beam tends to vanish for progressively higher values of Re.On the other hand, a more turbulent flow induces more complex oscillation for the beam.In order to assess the effectiveness of the coupling algorithm to properly enforce the interface conditions, the authors define an interface energy rate at the fluid-structure interface that is computed from both the fluid, J f , and the structure, J s , sides as the product between forces and velocities computed from the fluid and solid sides, respectively.As shown in Fig. 7, the two rates lie on the diagonal of the graph, thus showing values relatively close each other.Only a part of the graph moves apart from the diagonal, which corresponds to the beginning of the simulation when the velocity profile suddenly impacts the beams.Finally, some snapshots of the velocity field are given in Fig. 8.

CONCLUSIONS
n this paper, a numerical framework able to simulate fluid-structure interaction problems has been presented.Specifically, the fluid domain is solved by the LB method, whereas structure dynamics is predicted via the non-linear TDG method.The IB method is adopted to handle fluid-structure interface conditions.The moving free surface is treated by a VOF-based method.All the ingredients described are coupled in a global algorithm that has been applied in the two different test cases presented.In the former, a dam break phenomenon is simulated by adopting a free surface model.Present findings are compared to experimental data, showing a very close agreement.In the latter, a constant uniform rightward velocity profile invests two slender beams and scenarios characterized by different values of the Reynolds number are dissected.The effectiveness and the accuracy of the proposed approach are assessed by computing the interface energy rates from the fluid and solid sides, which show very close values.Based on numerical results, the I authors conclude that the proposed approach appears to be feasible and highly promising for simulating a wide range of fluid-structure interaction phenomena.

Figure 1 :
Figure 1: Handling of the missing distribution functions.Red dotted arrows represent the distribution that need to be reconstructed.

Figure 3 :
Figure 3: Dam break: numerical and experimental [23] time evolution of the free-surface profile.

Figure 5 :
Figure 5: Two slender beams obstructing a flow: sketch of the problem set-up.Moreover, the distance D is set to 100 and it is chosen as the characteristic length of the problem.At the west-most section, a constant uniform rightward velocity profile V=0.01 is prescribed, whereas at the east-most section outflow boundary conditions are enforced.The fluid possesses density =1, viscosity  and is initially at rest.The zero-velocity condition imposed at the bottom and top walls complete the definition of the problem.Notice that all the quantities are given in LB units.Beams possess elastic modulus E=0.01, cross sectional area A=16 and inertia moment J=21.33,

Figure 6 :
Figure 6: Two slender beams obstructing a flow: time history of the horizontal component of the tip displacement at different Reynolds number, Re=10 (red), 25 (green), 50 (blue).

Table 1 :
Dam break: time [ms] required to reach the checkpoints and impact the boundary.