An Algorithm for the Strong-Coupling of the Fluid-Structure Interaction Using a Staggered Approach

We present a staggered approach for the solution of the piston fluid-structure problem in a timedependent domain. The one-dimensional fluid flow is modelled using the nonlinear Euler equations. We investigate the time marching fluid-structure interaction and integrate the fluid and structure equations alternately using separate solvers. The Euler equations are written in moving mesh coordinates using the arbitrary Lagrangian-Eulerian ALE approach and discretised in space using the finite element method while the structure is integrated in time using an implicit finite difference Newmark-Wilson scheme. The influence of the time lag is studied by comparing two different structural predictors.


Introduction
Multiphysics problems involving a coupling between two or more different interacting physical phenomena encountered in engineering include fluid-structure interactions in aerodynamics, vibroaeroacoustic problems, the modeling of solidification and melting processes, and soft tissue mechanics, see Soulaïmani et al. 1 .These problems are particularly challenging to solve since the coupling calls for the use of different solvers in different parts of the solution domain, and with different mesh requirements thereby increasing the complexity of the computational effort.The problem of fluid-structure interaction FSI where fluid flow induces forces and thermal fluxes on a solid structure is well studied in the literature.The motion of the two phases modifies not only the fluid domain but also the velocities and the temperature fields at the fluid-structure interface.There are several practical and challenging examples of FSI, including, for instance, the interaction between the sail of a boat or a plane and the surrounding aerodynamic flow, the interaction between a bridge and the wind, and the interaction between vessels and blood flow Guruswamy 2 and Lefranc ¸ois and Boufflet 3 .
The two approaches that are commonly used to solve FSI problems are the staggered approach see, e.g., Farhat et al. 4 and Farhat and Lesoinne 5 , where a separate physics solver is used for each sub-domain, and the monolithic approach Blom 6 , where the fluid and structure are integrated in time as a single system.The staggered approach is the most widely used of the two and was used, for example, by Farhat 7 to study the panel flutter.Test cases were studied by Farhat and Lin 8 , Blom and Leyland 9 , Piperno 10, 11 , Prananta 12, 13 and Rausch et al. 14 .Blom and Leyland 15 and Piperno 10 studied the flutter of airfoils using this method.
In this paper the staggered methodis applied to the piston problem where the fluid domain is separated from a solid domain by a moving surface.The piston problem has previously been the subject of several studies.It was studied, for example, by Blom 6 and by Piperno 10,16 in his investigation of subcycling techniques on coupling algorithms.
Other studies of the piston problem include Fazio and Leveque 17 who developed a one-dimensional moving mesh method for the hyperbolic system of conservation laws based on a high-resolution finite-volume wave propagation method, which was implemented using the CLAWPACK software package.They solved the modified system of hyperbolic conservation laws on a fixed uniform computational grid, with a grid mapping function computed simultaneously in such a way that in the physical space certain features are tracked by cell interfaces.The scheme was tested on, among other test problems, a moving piston whose motion was tracked by the moving mesh.Chertock and Kurganov 18 extended the interface tracking method and developed a simple Eulerian finite-volume method for compressible fluid in domains with moving boundaries.The scheme was tested on the same piston problem taken from Fazio and Leveque 17 .Borsche et al. 19 considered a general system of conservation laws coupled with a system of ordinary differential equations.One of the test problems they considered is the classical piston problem.They assumed the flow to be isentropic and described the system by using the Lagrangian formulation.For the fluid solver they used a local Lax-Friedrichs scheme while using an explicit Euler method for the structure.The coupling was done after each time step.
Bendiksen 20 used a coupling algorithm in the Runge-Kutta scheme where the interaction was also at intermediate levels.Prananta and Hounjet 21 used structural as well as aerodynamical predictors in the algorithm.Van Zuijlen and Bijl 22 investigated the computational efficiency of higher-order partitioned implicit explicit IMEX schemes by comparing them with the monolithic second-order BDF scheme.They integrated both the fluid and the structure by the explicit singly diagonal implicit Runge-Kutta ESDIRK scheme.
The staggered approach, however, has its limitations, and it was shown in 23 that using high-order implicit schemes is not sufficient because, as the fluid and the structure domains are advanced in time successively, there will always be lags between the fluid and the solid phase solutions.Recently, it was introduced and advocated the staggered schemes using structural predictor 24 .The prediction techniques can reduce considerably the order of the energy conservation errors 25 .However, choosing a predictor is not always a straightforward prospect since some predictors are designed to work in combination with specific time integration methods to achieve stable schemes 24, 25 .Blom 6 used a staggered scheme with a structure predictor for linear acoustic equations.He used the zero-and firstorder predictors for the structural velocity in order to compare their performance.To integrate the structure he used the Newmark method, and to integrate the fluid modelled by acoustic equations he used the finite volume method.The results showed that, for linear acoustic simulations, the first-order predictor performed better than the zero-order predictor.
In order to gain an understanding of the structure predictors we take, in this paper the zeroth-and first-order structure predictors are used for nonlinear Euler equations.The performance of the two structure predictors is compared.In order to combine the predictors with the time integration scheme for the structure, we use the Newmark method as a structural solver and the arbitrary Lagrangian Euler ALE formulation for the Euler equations combined with the finite element method on a dynamic mesh.

The Piston Problem
We consider a gas contained in a one-dimensional chamber, closed on its right side by a moving piston and on its left by a fixed wall.The configuration is depicted in Figure 1.The piston has a mass m and is attached to an external fixed point with a linear spring stiffness k.The spring is characterized by three different lengths, namely, the unstretched length L s0 , the length at rest under pressure L se , and the length L s t at any given time during the interaction process.
The displacement, velocity, and acceleration of the piston are given, respectively, by q t , q t , and q t with regard to its position at rest.
For the movement of the piston we use the undamped differential equation of motion with one degree of freedom, namely, the classical mass-spring system where A is the piston cross-section, p x L 0 q t is the fluid pressure applied to the piston, and p 0 is the atmospheric pressure.

Modeling the Fluid Phase
The fluid phase is modelled using the inviscid compressible Euler equations.The fluid motion is assumed to vary only in the x-direction and is described by the one-dimensional nonlinear Euler equations where ρ, u, E and p denote the density, velocity, total energy, and pressure, respectively.The equations are closed by the equation of state for a perfect gas: where γ is the ratio of the specific heat.These three equations correspond to the conservation laws of mass, momentum, and total energy, respectively, and are usually combined in a vectorial form such as The flux F Q of any quantity Q e.g., mass, momentum, and energy is defined as the quantity flowing through a section S per unit time.For a fixed section, it is regarded as the local fluid velocity: where n defines the orientation vector of the section along the x-axis in this case .
The boundary delimiting the fluid domain moves in time, and it is necessary to solve the flow problem on a dynamic mesh 9, 15, 26 .A popular formulation for solving flow problems on dynamic meshes is the arbitrary Lagrangian Euler ALE 26 method.In this approach, the numerical scheme for solving the flow equations on moving grids typically incurs the computation of some geometric quantities involving the grid positions and velocities.The evaluation of these quantities as well as the time integration of fluxes on moving grids is accomplished using the discrete geometric conservation law DGCL 26 .We use the finite element method to discretize the one-dimensional flow in space.This involves computing the solution at discrete nodes within the fluid domain, such that two successive nodes form a finite element.The calculation mesh is shown in Figure 2.
Since the fluid has a moving boundary, a node attached to a movable boundary, such as the piston at x L t , must follow it 3, 6, 22 .In order to prevent nodes impinging or traversing, interior nodes must be moved, except for the node attached to the fixed boundary located at x 0, 3 .However, since 2.3 is not valid for moving nodes, it is necessary to rewrite these equations such that the flow term F i takes into account the motion of the nodes.We then consider that any point of the fluid domain is movable with a given velocity ω x , 6, 22 .This concept of moving coordinates is illustrated in Figure 3, where we consider a mesh composed of five nodes located at regular intervals along the domain and indexed from 1 to 5. Here, the bullet symbol • depicts the position at time t n and a circle symbol • depicts the position at time t n 1 .
It is essential to be able to move nodes in order to avoid the traversing node effect, visible in Figure 3 between nodes 4 and 5 wherein if nodes 4 and 5 are fixed, then they will be outside the problem domain at time t n 1 .

Arbitrary Lagrangian-Euler (ALE) Formulation
A combination of the Eulerian and Lagrangian fluid flow formulations, the arbitrary Lagrangian-Euler ALE formulation, is often preferred for general cases of fluid flow, 3, 6, 22, 27 .Correctly calculating the flow passing through a moving section is vital to ensure the conservation of fluid properties.For a section moving with velocity ω x , the flow rate given by 2.6 is corrected as to take into account the section motion.Applying the same corrective strategy, the conservative form of 2.5 may be written in arbitrary Lagrangian-Eulerian ALE form as where F F i − ω x U i is the corrected flow with respect to the moving space coordinate, ω x x, t defines the ALE grid velocity, and J x, t denotes the Jacobian of the frame transformation x t → ξ, where the substitution rule between x t and ξ is defined by The integration is performed on the fixed space ξ ∈ 0, L 0 .The boundary conditions for this coupled structure problem are given by a zero-flow condition at x 0 and by ensuring kinematic compatibility between the fluid flow and piston velocity at x L t , that is, u 0, t 0, u L t , t q t for t ≥ 0. 2.10

Numerical Schemes
In this section we describe the numerical schemes for the structure and for the fluid flow.

Structure Solver
We apply an implicit finite difference Newmark-Wilson scheme for time integration of 2.1 .
Further information about this method can be found in 3, 6, 28 .This scheme is based on the time series expansion of q and q; q n 1 q n Δt qn Δt where the variation between two successive times, Δq q n 1 − q n , is obtained by substituting 3.1 and 3.2 into 2.1 .The variables are updated at time t n 1 : This relation allows the new piston position, q n 1 , to be computed from q n , qn , and qn .Thus, the structural position is updated according to q n 1 q n Δq, 3.4 and the velocity qn 1 and acceleration qn 1 are updated according to 3.1 and 3.2 , respectively.For the first step q 1 of 3.4 the initial conditions are q 0 q 0 , q 0 q0 , and q 0 q0 .The first two initial conditions are given in 2.1 , and it is easy to show, also from 2.1 , that q0 1 m −kq 0 A p 0 − p 0 , 3.5 where p 0 is the uniform pressure in the chamber resulting from an adiabatic variation entailed by the initial change in the piston position q 0 .

Fluid Solver
We apply the finite element method 29 for the spatial discretization and a Lax-Wendroff scheme for the temporal resolution.The finite element method requires a variational form of 2.8 , where ψ ξ is any test function of class C 1 .Equation 3.6 may be written as In the above equation, the temporary derivative is evaluated at constant ξ.Switching back to the time-varying elements, 3.7 can be written as Integrating by parts the last term yields Finally, switching back again to constant ξ ∈ 0, L 0 , we have The time integration between t n and t n 1 gives The spatial discretization of the finite elements of the mesh followed by an assembling process gives where {U i } is the N × 1 global vector of unknowns of the ith equation in 2.8 , M n and M n 1 are the global mass matrices at times t n 1 and t n , respectively, and {R i } n 1/2 is an N × 1 global residual vector calculated at the halfway time step.Equation 3.12 is a system of equations to be solved at each time step.For numerical approximation of the flux, the explicit two-stage Lax-Wendroff procedure 30 where the state U LW i 1/2 is computed from is applied to solve the system 3.12 for each new time step t n 1 .The temporal stability criteria are given by the Courant-Friedrichs-Lewy CFL condition as where c 0 γRT is the local speed of sound.The Courant-Fredrichs-Lewy CFL condition ensues convergence of the finite element method.

Fluid Mesh Deformation Technique
In order to ensure kinematic compatibility between the fluid domain and the piston position and to prevent traversing by fluid nodes near the piston, fluid mesh deformation at each time step is necessary.Therefore, the mesh velocity is calculated by a linear interpolation of the left 0 and the right-hand V * velocities.The mesh nodes and mesh velocities for j 1, . . ., N are calculated from where V * is the piston velocity computed by the structure solver 6 .The boundary velocity in the fluid solver is calculated from see Blom 6 V * 1 2 qn 1 qn .3.17

Fluid-Structure Coupling: The Staggered Algorithm
The scheme chosen to couple the fluid and the structure is the staggered algorithm as described in Section 1.The staggered algorithm is applied with a structure predictor.At time t t n the state of the fluid, structure, and mesh are known.The next steps are taken to integrate the fluid-structure system from t n to t n 1 , and we proceed as follows see 3, 6 .
1 Predict the state of the structure at the end of the current time step t t n 1 .Here two different predictions are applied to the problem. 2 Integrate the fluid to the next time level using the predicted state of the structure.
The fluid is spatially discretized by the finite element method and integrated in time using the explicit two-stage Lax-Wendroff scheme.3 Update the structure to the next time level using the fluid pressures on the boundary.The structure is updated by the scheme described in Section 3.1.
The most important part of this algorithm is the prediction of the structural state at the end of the current time step t t n 1 23 .The velocity in a time step must correspond to the distance covered by the same time step 6 .The structural velocity is linear in time for the constant average acceleration method.The distance covered is calculated using the trapezoidal integration of the velocity.

Prediction
To predict the state of the structure at the end of the current time t t n 1 different predictors are applied to the problem: prediction 1 that is given by the zero-order prediction qn 1 qn 4.1 and prediction 2 given by a first-order prediction qn 1 qn Δt qn .

4.2
The velocity qn 1 in 3.17 is calculated using the structural prediction 4.1 or 4.2 .In the literature the performance of these predictors has been investigated for linear acoustic equations.In this paper we investigate their performance when nonlinear Euler equations are considered.

Results and Discussion
In Blom 6 the fluid domain was discretised in space by the finite volume method and integrated in time using an implicit upwind scheme.It is found that in the case of the acoustic equations, the predictor 2 algorithm given by 4.2 performed better than the prediction 1 algorithm given by 4.1 .We compared the performance of these two predictors by discretizing the fluid domain using the finite element method and by integrating in time using the explicit two-stage Lax-Wendroff procedure.
The numerical values for the parameters used are given in Table 1, where k is the spring rigidity, m is the mass of the piston, L 0 is the length of the chamber at rest, q 0 is the initial displacement of the piston, T 0 is the initial temperature, R is the universal gas constant and c 0 γRT 0 is the local speed of the sound.
Remark 5.1.The parameters in Table 1 are chosen according to the criteria given in Lefranc ¸ois and Boufflet 3 in order to ensure strong coupling effects.We introduce the notion of a characteristic time which is defined as follows 3 : i for the fluid it is the time required for a pressure wave to cross the chamber from one side to the other, that is, ii for the structure it is the natural period of the piston, that is, If the two characteristic times are similar, the coupling is strong.In the case where one of the characteristic times significantly exceeds the other, the dynamics in question fluid or structure can be considered as quasisteady and the coupling is weak 3 .
Taking into account the criteria given above and considering that T f char ≈ 0.0036 s.We have that for m 0.8 Kg, f 0 563 Hz, which implies that T s char 0.00177715 s.The characteristic times T f char and T s char are of the same order; therefore the coupling is strong 3 .In order to test the performance of the two coupling approaches with respect to the amplitude of the piston, the fluid flow is solved by means of the fluid solver described in Section 3.2.The fluid is discretised with 100 finite elements and 101 nodes.The initial conditions for the piston are taken as q 0 0.2 m and q 0 20 m/s.We compared the amplitude of the piston for the two coupling algorithms and different CFL numbers.Figure 4 shows the amplitude of the piston for the two coupling algorithms at CFL 0.80.Apparently, there is not too much difference between the solutions, but   one can see that the second curve, obtained by applying the predictor 2 algorithm, is slightly more damped than the first one.The difference between the two algorithms becomes more noticeable as the CFL number increases.Figures 5 and 6 show the amplitude of the piston computed by predictions 1 and 2 algorithms, respectively, as the CFL number increases.It is found that as the CFL number increases the curve becomes more damped.The damping is more pronounced in the curves shown in Figure 6, which are obtained by the prediction 2 algorithm.From the structure solver given in Section 3.1 we find that there is no damping, so the damping of the signal comes from the fluid solver.As the CFL number increases the deviation of the solution from the equilibrium position also increases, and this deviation is larger for the prediction 1 algorithm since it is less accurate.
We could find no discussion in the literature regarding the performance of the above coupling approaches with respect to the fluid.We, therefore, found it appropriate to test the performance of the two prediction schemes in this study.In order to do so, the fluid was discretized with 100 finite elements and then with 200 finite elements, respectively.Figures 7  and 8 show the results for the density and pressure, respectively, obtained via prediction 1, left-side, and via prediction 2, right side.It is possible to see that the results obtained by applying prediction 2 with 100 finite elements compared favorably with the results obtained by prediction 1 with 200 finite elements, showing the efficiency of prediction 2.

Conclusions
An analysis of the time marching fluid-structure interaction algorithm has been presented.The relatively simple piston problem was chosen in order to gain an understanding of the coupling algorithms.The one-dimensional fluid is modeled using the Euler equations, which were presented in moving mesh coordinates using the arbitrary Lagrangian Eulerian ALE approach, discretised in space by the finite element method and integrated in time using an explicit method.The structure was integrated in time by an implicit finite difference Newmark-Wilson scheme.The fluid and the structure were integrated in time using separate solvers.The coupling between the fluid and structure solvers was realized by applying the staggered approach.Since the staggered approach suffers from a time lag, the influence of the time lag was studied by comparing two different predictions for the structure.The computations show that the differences between the two coupling algorithms become noticeable as the CFL number increases.The prediction 2 algorithm gave a higher accuracy.

Figure 1 :
Figure 1: A gas enclosed in a chamber with a moving piston.

Figure 2 :
Figure 2: Discretization of the one-dimensional flow.

Figure 3 :
Figure 3: Moving physical space x t representation.

Figure 4 :
Figure 4: Amplitude of the piston at CFL 0.80.

Figure 8 :
Figure 8: Pressure p at T max 3 × T 0 , staggered algorithm, via prediction 1 a and prediction 2 b .

Table 1 :
Numerical values of the parameters for the piston problem.