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 𝐿𝑠0, the length at rest under pressure 𝐿se, and the length 𝐿𝑠(𝑡) at any given time during the interaction process.

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 𝑝𝑚̈𝑞(𝑡)+𝑘𝑞(𝑡)=𝐴𝑥=𝐿0+𝑞(𝑡)𝑝0,(2.1) complemented by initial conditions 𝑞(0)=𝑞0,̇𝑞(0)=0,(2.2) where 𝐴 is the piston cross-section, 𝑝(𝑥=𝐿0+𝑞(𝑡)) is the fluid pressure applied to the piston, and 𝑝0 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 𝜕𝜌+𝜕𝜕𝑡𝜌𝑢𝜌𝐸𝜕𝑥𝜌𝑢𝜌𝑢2[]+𝑝𝑢(𝜌𝐸+𝑝)=0,𝑥0,𝐿,𝑡0,(2.3) where 𝜌, 𝑢, 𝐸 and 𝑝 denote the density, velocity, total energy, and pressure, respectively. The equations are closed by the equation of state for a perfect gas: 1𝑝=(𝛾1)𝜌𝐸2𝜌𝑢2,(2.4) 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 𝜕𝜕𝜕𝑡𝑈+𝜌𝜕𝑥𝐹=0,where𝑈=𝜌𝑢𝜌𝐸,𝐹=𝜌𝑢𝜌𝑢2+𝑝𝑢(𝜌𝐸+𝑝).(2.5)

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: 𝐹(𝑄)=𝑆𝑄𝑢𝑛𝑑𝑆,(2.6) 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.

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 𝑥=0, [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 𝑡𝑛+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 𝑡𝑛+1.

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 𝐹(𝑄)=𝑆𝑄𝑢𝜔𝑥𝑛𝑑𝑆,(2.7) 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 𝜕𝜕𝑡𝐽𝑈𝑖𝜕+𝐽𝐹𝜕𝑥𝑖𝜔𝑥𝑈𝑖=0for𝑖=1,2,3,(2.8) 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 𝐽(𝑥,𝑡)=𝑑𝑥(𝜉,𝑡),𝑑𝜉𝑥(𝑡)𝑓(𝑥,𝑡)𝑑𝑥=𝜉𝑓(𝑥(𝜉,𝑡),𝑡)𝐽𝑑𝜉.(2.9) The integration is performed on the fixed space 𝜉[0,𝐿0]. The boundary conditions for this coupled structure problem are given by a zero-flow condition at 𝑥=0 and by ensuring kinematic compatibility between the fluid flow and piston velocity at 𝑥=𝐿(𝑡), that is, 𝑢(0,𝑡)=0,𝑢(𝐿(𝑡),𝑡)=̇𝑞(𝑡)for𝑡0.(2.10)

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 ̇𝑞; 𝑞𝑛+1=𝑞𝑛+Δ𝑡̇𝑞𝑛+Δ𝑡24̈𝑞𝑛+̈𝑞𝑛+1,̇𝑞𝑛+1=̇𝑞𝑛+Δ𝑡2̈𝑞𝑛+̈𝑞𝑛+1,(3.1) where 𝑛 and 𝑛+1 correspond to times 𝑡 and 𝑡+Δ𝑡, respectively, and Δ𝑡 is the time step between two successive solutions. From the first relation in (3.1) we have that ̈𝑞𝑛+1=4Δ𝑡24Δ𝑞Δ𝑡̇𝑞𝑛̈𝑞𝑛,(3.2) where the variation between two successive times, Δ𝑞=𝑞𝑛+1𝑞𝑛, is obtained by substituting (3.1) and (3.2) into (2.1). The variables are updated at time 𝑡𝑛+1: 1Δ𝑞=4𝑚+𝑘Δ𝑡2𝐴Δ𝑡2𝑝𝑛𝑝0+𝑚4Δ𝑡̇𝑞𝑛+Δ𝑡2̈𝑞𝑛+1Δ𝑡2𝑘𝑞𝑛.(3.3) This relation allows the new piston position, 𝑞𝑛+1, to be computed from 𝑞𝑛, ̇𝑞𝑛, and ̈𝑞𝑛. Thus, the structural position is updated according to 𝑞𝑛+1=𝑞𝑛+Δ𝑞,(3.4) and the velocity ̇𝑞𝑛+1 and acceleration ̈𝑞𝑛+1 are updated according to (3.1) and (3.2), respectively. For the first step 𝑞1 of (3.4) the initial conditions are 𝑞(0)=𝑞0, ̇𝑞(0)=̇𝑞0, and ̈𝑞(0)=̈𝑞0. The first two initial conditions are given in (2.1), and it is easy to show, also from (2.1), that ̈𝑞0=1𝑚𝑘𝑞0+𝐴𝑝(0)𝑝0,(3.5) where 𝑝(0) is the uniform pressure in the chamber resulting from an adiabatic variation entailed by the initial change in the piston position 𝑞0.

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), 𝐿00𝜕𝐽𝑈𝑖𝜕𝐹𝜕𝑡+𝐽𝑖𝜕𝑥𝜓(𝜉)𝑑𝜉=0,for𝑖=1,2,3,(3.6) where 𝜓(𝜉) is any test function of class 𝐂1. Equation (3.6) may be written as 𝐿00𝜓𝜕𝐽𝑈𝑖𝜕𝑡𝑑𝜉+𝐿00𝐽𝜕𝐹𝜕𝑥𝑖𝜓𝑑𝜉=0.(3.7) In the above equation, the temporary derivative is evaluated at constant 𝜉. Switching back to the time-varying elements, (3.7) can be written as 0𝐿(𝑡)𝜓𝜕𝑈𝑖𝜕𝑡𝑑𝑥+0𝐿(𝑡)𝜕𝐹𝜕𝑥𝑖𝜓𝑑𝑥=0.(3.8) Integrating by parts the last term yields 0𝐿(𝑡)𝜓𝜕𝑈𝑖𝜕𝑡𝑑𝑥0𝐿(𝑡)𝜕𝜓𝐹𝜕𝑥𝑖𝜓𝐹𝑑𝑥+𝑖0𝐿(𝑡)=0.(3.9) Finally, switching back again to constant 𝜉[0,𝐿0], we have 𝐿00𝜓𝜕𝐽𝑈𝑖𝜕𝑡𝑑𝜉𝐿00𝜕𝜓𝐹𝜕𝜉𝑖𝜓𝐹𝑑𝜉+𝑖𝐿00=0.(3.10)

The time integration between 𝑡𝑛 and 𝑡𝑛+1 gives 𝐿00𝜓𝐽𝑈𝑖𝑛+1𝑑𝜉𝐿00(𝐽𝑈)𝑛𝑑𝜉Δ𝑡𝐿00𝜕𝜓𝐹𝜕𝜉𝑖𝑛+1/2𝜓𝐹𝑑𝜉+𝑖𝑛+1/2𝐿00=0.(3.11) The spatial discretization of the finite elements of the mesh followed by an assembling process gives [𝑀]𝑛+1𝑈𝑖𝑛+1[𝑀]𝑛𝑈𝑖𝑅(𝑛)Δ𝑖𝑛+1/2={0},for𝑖=1,2,3,(3.12) where {𝑈𝑖} is the (𝑁×1) global vector of unknowns of the 𝑖th equation in (2.8), [𝑀]𝑛 and [𝑀]𝑛+1 are the global mass matrices at times 𝑡𝑛+1 and 𝑡𝑛, respectively, and {𝑅𝑖}𝑛+(1/2) is an (𝑁×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] 𝐹𝐿𝑊2𝑖+1/2=𝐹𝑈𝐿𝑊𝑖+1/2,(3.13) where the state 𝑈𝐿𝑊𝑖+1/2 is computed from 𝑈𝐿𝑊𝑖+1/2=12𝑈𝑛𝑖+1+𝑈𝑛𝑖+12Δ𝑡𝐹𝑈Δ𝑥𝑛𝑖𝐹𝑈𝑛𝑖+1,(3.14) is applied to solve the system (3.12) for each new time step 𝑡𝑛+1. The temporal stability criteria are given by the Courant-Friedrichs-Lewy (CFL) condition as 𝐿Δ𝑡=CFL×min𝑒||𝑢+𝑐0+𝑤𝑥||,withCFL<1,(3.15) where 𝑐0=𝛾𝐑𝐓 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 (0) and the right-hand 𝑉 velocities. The mesh nodes and mesh velocities for 𝑗=1,,𝑁 are calculated from 𝜔𝑥𝑗𝑛=𝑗1𝑉𝑁1,𝑥𝑗𝑛+1=𝑥𝑗𝑛+Δ𝑡𝜔𝑥𝑗𝑛,(3.16) where 𝑉 is the piston velocity computed by the structure solver [6]. The boundary velocity in the fluid solver is calculated from (see Blom [6]) 𝑉=12̇𝑞𝑛+1+̇𝑞𝑛.(3.17)

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 𝑡𝑛+1, and we proceed as follows (see [3, 6]) (1)Predict the state of the structure at the end of the current time step (𝑡=𝑡𝑛+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 𝑡=𝑡𝑛+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.

4.1. Prediction

To predict the state of the structure at the end of the current time 𝑡=𝑡𝑛+1 different predictors are applied to the problem: prediction 1 that is given by the zero-order prediction ̇𝑞𝑛+1={̇𝑞𝑛}(4.1) and prediction 2 given by a first-order prediction ̇𝑞𝑛+1={̇𝑞𝑛}+Δ𝑡{̈𝑞𝑛}.(4.2) The velocity ̇𝑞𝑛+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.

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, 𝐿0 is the length of the chamber at rest, 𝑞0 is the initial displacement of the piston, 𝐓0 is the initial temperature, 𝐑 is the universal gas constant and 𝑐0=𝛾𝐓0 is the local speed of the sound.

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, 𝒯𝑓char𝐿(𝑡)𝑐0,(5.1)(ii)for the structure it is the natural period of the piston, that is, 𝒯𝑠char=𝒯0.(5.2)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 𝒯𝑓char0.0036 s. We have that for 𝑚=0.8Kg, 𝑓0=563 Hz, which implies that 𝒯𝑠char=0.00177715 s. The characteristic times 𝒯𝑓char and 𝒯𝑠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 𝑞(0)=0.2 m and ̇𝑞(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.

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.