`ISRN Applied MathematicsVolume 2012 (2012), Article ID 391974, 14 pageshttp://dx.doi.org/10.5402/2012/391974`
Research Article

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

School of Mathematical Sciences, University of KwaZulu-Natal, Private Bag X01, Scottsville, Pietermaritzburg 3209, South Africa

Received 10 February 2012; Accepted 1 March 2012

Academic Editors: H. Homeier and C. Lu

Copyright © 2012 Josè C. Pedro and P. Sibanda. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

We present a staggered approach for the solution of the piston fluid-structure problem in a time-dependent 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.

#### 1. 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 Lefranç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 first-order 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.

#### 2. 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 and is attached to an external fixed point with a linear spring stiffness . The spring is characterized by three different lengths, namely, the unstretched length , the length at rest under pressure , and the length at any given time during the interaction process.

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

The displacement, velocity, and acceleration of the piston are given, respectively, by , , and 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 complemented by initial conditions where is the piston cross-section, is the fluid pressure applied to the piston, and is the atmospheric pressure.

##### 2.1. 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 -direction and is described by the one-dimensional nonlinear Euler equations where , , and 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 of any quantity (e.g., mass, momentum, and energy) is defined as the quantity flowing through a section per unit time. For a fixed section, it is regarded as the local fluid velocity: where defines the orientation vector of the section (along the -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.

Figure 2: Discretization of the one-dimensional flow.

Since the fluid has a moving boundary, a node attached to a movable boundary, such as the piston at , 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 , [3]. However, since (2.3) is not valid for moving nodes, it is necessary to rewrite these equations such that the flow term takes into account the motion of the nodes. We then consider that any point of the fluid domain is movable with a given velocity , [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 and a circle symbol () depicts the position at time .

Figure 3: Moving physical space representation.

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 .

##### 2.2. 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 , 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 is the corrected flow with respect to the moving space coordinate, defines the ALE grid velocity, and denotes the Jacobian of the frame transformation , where the substitution rule between and is defined by The integration is performed on the fixed space . The boundary conditions for this coupled structure problem are given by a zero-flow condition at and by ensuring kinematic compatibility between the fluid flow and piston velocity at , that is,

#### 3. Numerical Schemes

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

##### 3.1. 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 and ; where and correspond to times and , respectively, and is the time step between two successive solutions. From the first relation in (3.1) we have that where the variation between two successive times, , is obtained by substituting (3.1) and (3.2) into (2.1). The variables are updated at time : This relation allows the new piston position, , to be computed from , , and . Thus, the structural position is updated according to and the velocity and acceleration are updated according to (3.1) and (3.2), respectively. For the first step of (3.4) the initial conditions are , , and . The first two initial conditions are given in (2.1), and it is easy to show, also from (2.1), that where is the uniform pressure in the chamber resulting from an adiabatic variation entailed by the initial change in the piston position .

##### 3.2. 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 . 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 , we have

The time integration between and gives The spatial discretization of the finite elements of the mesh followed by an assembling process gives where is the global vector of unknowns of the th equation in (2.8), and are the global mass matrices at times and , respectively, and is an 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 is computed from is applied to solve the system (3.12) for each new time step . The temporal stability criteria are given by the Courant-Friedrichs-Lewy (CFL) condition as where is the local speed of sound. The Courant-Fredrichs-Lewy (CFL) condition ensues convergence of the finite element method.

###### 3.2.1. 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 and the right-hand velocities. The mesh nodes and mesh velocities for are calculated from where is the piston velocity computed by the structure solver [6]. The boundary velocity in the fluid solver is calculated from (see Blom [6])

#### 4. 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 the state of the fluid, structure, and mesh are known. The next steps are taken to integrate the fluid-structure system from to , and we proceed as follows (see [3, 6]) (1)Predict the state of the structure at the end of the current time step . 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 [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.

##### 4.1. Prediction

To predict the state of the structure at the end of the current time different predictors are applied to the problem: prediction 1 that is given by the zero-order prediction and prediction 2 given by a first-order prediction The velocity 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.

#### 5. 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 is the spring rigidity, is the mass of the piston, is the length of the chamber at rest, is the initial displacement of the piston, is the initial temperature, is the universal gas constant and is the local speed of the sound.

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

Remark 5.1. The parameters in Table 1 are chosen according to the criteria given in Lefranç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  s. We have that for ,  Hz, which implies that  s. The characteristic times and 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  m and  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 . 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.

Figure 4: Amplitude of the piston at .
Figure 5: Amplitude of the piston by prediction 1 algorithm.
Figure 6: Amplitude of the piston by prediction 2 algorithm at and , respectively (a). Zoom to emphasize the damping (b).

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.

Figure 7: Density at , staggered algorithm, via prediction 1 (a) and prediction 2 (b).
Figure 8: Pressure at , staggered algorithm, via prediction 1 (a) and prediction 2 (b).

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

#### References

1. A. Soulaïmani, Z. Feng, and A. B. H. Ali, “Solution techniques for multi-physics problems with application to computational nonlinear aeroelasticity,” Nonlinear Analysis, Theory, Methods and Applications, vol. 63, no. 5–7, pp. e1585–e1595, 2005.
2. G. P. Guruswamy, “A review of numerical fluids/structures interface methods for computations using high-fidelity equations,” Computers and Structures, vol. 80, no. 1, pp. 31–41, 2002.
3. E. Lefrançois and J.-P. Boufflet, “An introduction to fluid-structure interaction: application to the piston problem,” SIAM Review, vol. 52, no. 4, pp. 747–767, 2010.
4. C. Farhat, K. G. van der Zee, and P. Geuzaine, “Provably second-order time-accurate loosely-coupled solution algorithms for transient nonlinear computational aeroelasticity,” Computer Methods in Applied Mechanics and Engineering, vol. 195, no. 17-18, pp. 1973–2001, 2006.
5. C. Farhat and M. Lesoinne, “Two efficient staggered algorithms for the serial and parallel solution of three-dimensional nonlinear transient aeroelastic problems,” Computer Methods in Applied Mechanics and Engineering, vol. 182, no. 3-4, pp. 499–515, 2000.
6. F. J. Blom, “A monolithical fluid-structure interaction algorithm applied to the piston problem,” Computer Methods in Applied Mechanics and Engineering, vol. 167, no. 3-4, pp. 369–391, 1998.
7. C. Farhat, “High performance simulation of coupled nonlinear transient aeroelastic problems,” AGARD Report R-807, Cosmase Course, Lausanne, Switzerland, 1995.
8. C. Farhat and T. Y. Lin, “A structure attached corotational fluid grid for transient aeroelastic compu- tations,” AIAA Journal, vol. 31, no. 3, pp. 597–599, 1993.
9. F. J. Blom and P. Leyland, “Analysis of fluid-structure interaction by means of dynamic unstructured meshes,” in Proceedings of the 4th International Mechanical Engineering Congress and Expositionon Fluid-Structure Interaction, Aeroelasticity, Flow-Induced Vibrations and Noise, Paudoussis et al., Ed., vol. 53-1, pp. 3–10, Dallas, Tex, USA, 1997.
10. S. Piperno, Simulation numéric de phénomènes d'interaction fluid-structure [Ph.D. thesis], Ecole Nationale Des Ponts et Chausses, 1995.
11. S. Piperno, “Two-dimensional Euler aeroelastic simulations with interface matching relaxation,” in Proceedings of the 2nd ECCOMAS Conference on Numerical Methods in Engineering, pp. 898–904, 1996.
12. B. B. Prananta, M. H. L. Hounjet, and R. J. Zwaan, “Thin layer Navier Stokes and its application for aeroelastic analysis of an airfoil in transonic flow,” in Proceedings of the International Forum on Aeroelasticity and Structural Dynamics, Manchester, UK, 1996.
13. B. B. Prananta and M. H. L. Hounjet, “Aeroelastic simulation with advanced CFD methods in 2D and 3D transonic flow,” in Proceedings of the Unsteady Aerodynamics Conference Royal Aeronautical Socuiety, London, UK, 1996.
14. R. D. Rausch, J. T. Batina, and H. T. Y. Yang, “Euler flutter analysis of airfoils using unstructured dynamic meshes,” Journal of Aircraft, vol. 27, no. 5, pp. 436–443, 1990.
15. F. J. Blom and P. Leyland, “Analysis of fluid-structure interaction on moving airfoils by means of an improved ALE-method,” AIAA Paper number 97-1770, 1997.
16. S. Piperno, “Staggered time-integration methods for a one-dimensional Euler aeroelastic problem,” Raport de Recherche CERMICS number 94-33, 1994.
17. R. Fazio and R. J. LeVeque, “Moving-mesh methods for one-dimensional hyperbolic problems using CLAWPACK,” Computers & Mathematics with Applications, vol. 45, no. 1–3, pp. 273–298, 2003.
18. A. Chertock and A. Kurganov, “A simple Eulerian finite-volume method for compressible fluids in domains with moving boundaries,” Communications in Mathematical Sciences, vol. 6, no. 3, pp. 531–556, 2008.
19. R. Borsche, R. M. Colombo, and M. Garavello, “On the coupling of systems of hyperbolic conservation laws with ordinary differential equations,” Nonlinearity, vol. 23, no. 11, pp. 2749–2770, 2010.
20. O. O. Bendiksen, “A New Approach to Computational Aeroelasticity,” AIAA Paper number 91-0939-CP, 1991.
21. B. B. Prananta and M. H. L. Hounjet, “Large time step aero-structural coupling procedures for aeroelastic simulation,” in Proceedings of the International Forum on Aeroelasticity and structural Dynamics, vol. 2, pp. 63–70, Rome, Italy, 1997.
22. A. H. Van Zuijlen and H. Bijl, “Implicit and explicit higher order time integration schemes for structural dynamics and fluid-structure interaction computations,” Computers and Structures, vol. 83, no. 2-3, pp. 93–105, 2005.
23. V. Jean-Mark, D. V. Pascal, H. Charles, and L. Benoit, “Strong coupling algorithm to solve fluid-structure interaction problems with a staggered approach,” Report, Open Engineering SA, 2009.
24. S. Piperno, “Explicit/implicit fluid/structure staggered procedures with a structural predictor and fluid subcycling for 2D inviscid aeroelastic simulations,” International Journal for Numerical Methods in Fluids, vol. 25, no. 10, pp. 1207–1226, 1997.
25. C. Michler, S. J. Hulshoff, E. H. van Brummelen, and R. de Borst, “A monolithic approach to fluid-structure interaction,” Computers and Fluids, vol. 33, no. 5-6, pp. 839–848, 2004.
26. C. Farhat, P. Geuzaine, and C. Grandmont, “The discrete geometric conservation law and the nonlinear stability of ALE schemes for the solution of flow problems on moving grids,” Journal of Computational Physics, vol. 174, no. 2, pp. 669–694, 2001.
27. R. Löhner, Applied Computational Fluid Dynamics Techniques: An Introduction Based on Finite Element Methods, John Wiley & Sons, 2nd edition, 2008.
28. L. Bardella and F. Genna, “Newmark's time integration method from the discretization of extended functionals,” Transactions of the ASME on Journal of Applied Mechanics, vol. 72, no. 4, pp. 527–537, 2005.
29. D. S. Burnett, Finite Element Analysis, Addison-Wesley Publishing Company, Reading, Mass, USA, 1987.
30. E. F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics, Springer, Berlin, Germany, 2nd edition, 1999.