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

Authors: Josè C. Pedro; P. Sibanda
Journal: ISRN Applied Mathematics
 (2012)

---

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


---

## Body

## 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 Figure1. 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 Ls0, the length at rest under pressure Lse, and the length Ls(t) 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, byq(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(2.1)mq¨(t)+kq(t)=A(p(x=L0+q(t))-p0),
complemented by initial conditions
(2.2)q(0)=q0,q˙(0)=0,
where A is the piston cross-section, p(x=L0+q(t)) is the fluid pressure applied to the piston, and p0 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 thex-direction and is described by the one-dimensional nonlinear Euler equations
(2.3)∂∂t{ρρuρE}+∂∂x{ρuρu2+pu(ρE+p)}=0,∀x∈[0,L],t≥0,
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:
(2.4)p=(γ-1)(ρE-12ρu2),
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
(2.5)∂∂tU+∂∂xF=0,whereU={ρρuρE},F={ρuρu2+pu(ρE+p)}.The fluxF(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:
(2.6)F(Q)=∫S‍Qu→⋅n→dS,
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.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 atx=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 Fi 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 tn and a circle symbol (∘) depicts the position at time tn+1.Figure 3
Moving physical spacex(t) representation.It is essential to be able to move nodes in order to avoid the traversing node effect, visible in Figure3 between nodes 4 and 5 wherein if nodes 4 and 5 are fixed, then they will be outside the problem domain at time tn+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 ωx, the flow rate given by (2.6) is corrected as
(2.7)F~(Q)=∫SQ(u→-ω→x)⋅n→dS,
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
(2.8)∂∂t(JUi)+J∂∂x(Fi-ωxUi)=0fori=1,2,3,
where F~=(Fi-ωxUi) 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
(2.9)J(x,t)=dx(ξ,t)dξ,∫x(t)f(x,t)dx=∫ξf(x(ξ,t),t)Jdξ.
The integration is performed on the fixed space ξ∈[0,L0]. 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,
(2.10)u(0,t)=0,u(L(t),t)=q˙(t)fort≥0.


## 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 thex-direction and is described by the one-dimensional nonlinear Euler equations
(2.3)∂∂t{ρρuρE}+∂∂x{ρuρu2+pu(ρE+p)}=0,∀x∈[0,L],t≥0,
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:
(2.4)p=(γ-1)(ρE-12ρu2),
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
(2.5)∂∂tU+∂∂xF=0,whereU={ρρuρE},F={ρuρu2+pu(ρE+p)}.The fluxF(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:
(2.6)F(Q)=∫S‍Qu→⋅n→dS,
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.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 atx=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 Fi 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 tn and a circle symbol (∘) depicts the position at time tn+1.Figure 3
Moving physical spacex(t) representation.It is essential to be able to move nodes in order to avoid the traversing node effect, visible in Figure3 between nodes 4 and 5 wherein if nodes 4 and 5 are fixed, then they will be outside the problem domain at time tn+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 ωx, the flow rate given by (2.6) is corrected as
(2.7)F~(Q)=∫SQ(u→-ω→x)⋅n→dS,
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
(2.8)∂∂t(JUi)+J∂∂x(Fi-ωxUi)=0fori=1,2,3,
where F~=(Fi-ωxUi) 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
(2.9)J(x,t)=dx(ξ,t)dξ,∫x(t)f(x,t)dx=∫ξf(x(ξ,t),t)Jdξ.
The integration is performed on the fixed space ξ∈[0,L0]. 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,
(2.10)u(0,t)=0,u(L(t),t)=q˙(t)fort≥0.


## 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 q and q˙;
(3.1)qn+1=qn+Δtq˙n+Δt24(q¨n+q¨n+1),q˙n+1=q˙n+Δt2(q¨n+q¨n+1),
where n and n+1 correspond to times t and t+Δt, respectively, and Δt is the time step between two successive solutions. From the first relation in (3.1) we have that
(3.2)q¨n+1=4Δt2Δq-4Δtq˙n-q¨n,
where the variation between two successive times, Δq=qn+1-qn, is obtained by substituting (3.1) and (3.2) into (2.1). The variables are updated at time tn+1:
(3.3)Δq=14m+kΔt2[AΔt2(pn-p0)+m(4Δtq˙n+Δt2q¨n+1)-Δt2kqn].
This relation allows the new piston position, qn+1, to be computed from qn, q˙n, and q¨n. Thus, the structural position is updated according to
(3.4)qn+1=qn+Δq,
and the velocity q˙n+1 and acceleration q¨n+1 are updated according to (3.1) and (3.2), respectively. For the first step q1 of (3.4) the initial conditions are q(0)=q0, q˙(0)=q˙0, and q¨(0)=q¨0. The first two initial conditions are given in (2.1), and it is easy to show, also from (2.1), that
(3.5)q¨0=1m(-kq0+A(p(0)-p0)),
where p(0) is the uniform pressure in the chamber resulting from an adiabatic variation entailed by the initial change in the piston position q0.


### 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),
(3.6)∫0L0[∂(JUi)∂t+J∂F~i∂x]ψ(ξ)dξ=0,fori=1,2,3,
where ψ(ξ) is any test function of class C1. Equation (3.6) may be written as
(3.7)∫0L0ψ∂(JUi)∂tdξ+∫0L0J∂∂xF~iψdξ=0.
In the above equation, the temporary derivative is evaluated at constant ξ. Switching back to the time-varying elements, (3.7) can be written as
(3.8)∫0L(t)ψ∂Ui∂tdx+∫0L(t)∂∂xF~iψdx=0.
Integrating by parts the last term yields
(3.9)∫0L(t)ψ∂Ui∂tdx-∫0L(t)∂ψ∂xF~idx+[ψF~i]0L(t)=0.
Finally, switching back again to constant ξ∈[0,L0], we have
(3.10)∫0L0ψ∂JUi∂tdξ-∫0L0∂ψ∂ξF~idξ+[ψF~i]0L0=0.The time integration betweentn and tn+1 gives
(3.11)∫0L0ψ(JUi)n+1dξ-∫0L0(JU)ndξ-Δt(∫0L0∂ψ∂ξF~in+1/2dξ+[ψF~in+1/2]0L0)=0.
The spatial discretization of the finite elements of the mesh followed by an assembling process gives
(3.12)[M]n+1{Ui}n+1-[M]n{Ui}(n)-Δ{Ri}n+1/2={0},fori=1,2,3,
where {Ui} 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 tn+1 and tn, respectively, and {Ri}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]
(3.13)F~i+1/2LW2=F~(Ui+1/2LW),
where the state Ui+1/2LW is computed from
(3.14)Ui+1/2LW=12(Ui+1n+Uin)+12ΔtΔx(F~(Uin)-F~(Ui+1n)),
is applied to solve the system (3.12) for each new time step tn+1. The temporal stability criteria are given by the Courant-Friedrichs-Lewy (CFL) condition as
(3.15)Δt=CFL×min⁡(Le|u+c0+wx|),withCFL<1,
where c0=γRT 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 V* velocities. The mesh nodes and mesh velocities for j=1,…,N are calculated from
(3.16)ωxjn=j-1N-1V*,xjn+1=xjn+Δtωxjn,
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])
(3.17)V*=12(q˙n+1+q˙n).


## 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 q and q˙;
(3.1)qn+1=qn+Δtq˙n+Δt24(q¨n+q¨n+1),q˙n+1=q˙n+Δt2(q¨n+q¨n+1),
where n and n+1 correspond to times t and t+Δt, respectively, and Δt is the time step between two successive solutions. From the first relation in (3.1) we have that
(3.2)q¨n+1=4Δt2Δq-4Δtq˙n-q¨n,
where the variation between two successive times, Δq=qn+1-qn, is obtained by substituting (3.1) and (3.2) into (2.1). The variables are updated at time tn+1:
(3.3)Δq=14m+kΔt2[AΔt2(pn-p0)+m(4Δtq˙n+Δt2q¨n+1)-Δt2kqn].
This relation allows the new piston position, qn+1, to be computed from qn, q˙n, and q¨n. Thus, the structural position is updated according to
(3.4)qn+1=qn+Δq,
and the velocity q˙n+1 and acceleration q¨n+1 are updated according to (3.1) and (3.2), respectively. For the first step q1 of (3.4) the initial conditions are q(0)=q0, q˙(0)=q˙0, and q¨(0)=q¨0. The first two initial conditions are given in (2.1), and it is easy to show, also from (2.1), that
(3.5)q¨0=1m(-kq0+A(p(0)-p0)),
where p(0) is the uniform pressure in the chamber resulting from an adiabatic variation entailed by the initial change in the piston position q0.


## 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),
(3.6)∫0L0[∂(JUi)∂t+J∂F~i∂x]ψ(ξ)dξ=0,fori=1,2,3,
where ψ(ξ) is any test function of class C1. Equation (3.6) may be written as
(3.7)∫0L0ψ∂(JUi)∂tdξ+∫0L0J∂∂xF~iψdξ=0.
In the above equation, the temporary derivative is evaluated at constant ξ. Switching back to the time-varying elements, (3.7) can be written as
(3.8)∫0L(t)ψ∂Ui∂tdx+∫0L(t)∂∂xF~iψdx=0.
Integrating by parts the last term yields
(3.9)∫0L(t)ψ∂Ui∂tdx-∫0L(t)∂ψ∂xF~idx+[ψF~i]0L(t)=0.
Finally, switching back again to constant ξ∈[0,L0], we have
(3.10)∫0L0ψ∂JUi∂tdξ-∫0L0∂ψ∂ξF~idξ+[ψF~i]0L0=0.The time integration betweentn and tn+1 gives
(3.11)∫0L0ψ(JUi)n+1dξ-∫0L0(JU)ndξ-Δt(∫0L0∂ψ∂ξF~in+1/2dξ+[ψF~in+1/2]0L0)=0.
The spatial discretization of the finite elements of the mesh followed by an assembling process gives
(3.12)[M]n+1{Ui}n+1-[M]n{Ui}(n)-Δ{Ri}n+1/2={0},fori=1,2,3,
where {Ui} 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 tn+1 and tn, respectively, and {Ri}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]
(3.13)F~i+1/2LW2=F~(Ui+1/2LW),
where the state Ui+1/2LW is computed from
(3.14)Ui+1/2LW=12(Ui+1n+Uin)+12ΔtΔx(F~(Uin)-F~(Ui+1n)),
is applied to solve the system (3.12) for each new time step tn+1. The temporal stability criteria are given by the Courant-Friedrichs-Lewy (CFL) condition as
(3.15)Δt=CFL×min⁡(Le|u+c0+wx|),withCFL<1,
where c0=γRT 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 V* velocities. The mesh nodes and mesh velocities for j=1,…,N are calculated from
(3.16)ωxjn=j-1N-1V*,xjn+1=xjn+Δtωxjn,
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])
(3.17)V*=12(q˙n+1+q˙n).


## 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 V* velocities. The mesh nodes and mesh velocities for j=1,…,N are calculated from
(3.16)ωxjn=j-1N-1V*,xjn+1=xjn+Δtωxjn,
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])
(3.17)V*=12(q˙n+1+q˙n).


## 4. Fluid-Structure Coupling: The Staggered Algorithm

The scheme chosen to couple the fluid and the structure is the staggered algorithm as described in Section1. The staggered algorithm is applied with a structure predictor. At time t=tn the state of the fluid, structure, and mesh are known. The next steps are taken to integrate the fluid-structure system from tn to tn+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=tn+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 Section3.1.The most important part of this algorithm is the prediction of the structural state at the end of the current time step t=tn+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 timet=tn+1 different predictors are applied to the problem: prediction 1 that is given by the zero-order prediction
(4.1){q˙n+1}={q˙n}
and prediction 2 given by a first-order prediction
(4.2){q˙n+1}={q˙n}+Δt{q¨n}.
The velocity q˙n+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.


## 4.1. Prediction

To predict the state of the structure at the end of the current timet=tn+1 different predictors are applied to the problem: prediction 1 that is given by the zero-order prediction
(4.1){q˙n+1}={q˙n}
and prediction 2 given by a first-order prediction
(4.2){q˙n+1}={q˙n}+Δt{q¨n}.
The velocity q˙n+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 Table1, where k is the spring rigidity, m is the mass of the piston, L0 is the length of the chamber at rest, q0 is the initial displacement of the piston, T0 is the initial temperature, R is the universal gas constant and c0=γℛT0 is the local speed of the sound.Table 1
Numerical values of the parameters for the piston problem.
k
m
L
0
q
0
L
s
0
A
γ
p
0
T
0
ℛ
c
0
107 N/m
0.8 Kg
1 m
0.20 m
1.2 m
1
1.4
105 Pa
300 K
287
334.7 m/sRemark 5.1.
The parameters in Table1 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 acharacteristic 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,
(5.1)Tcharf≈L(t)c0,(ii)
for the structure it is the natural period of the piston, that is,
(5.2)Tchars=T0.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𝒯charf≈0.0036 s. We have that for m=0.8 Kg, f0=563 Hz, which implies that 𝒯chars=0.00177715 s. The characteristic times 𝒯charf and 𝒯chars 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 Section3.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. Figure4 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.Amplitude of the piston atCFL=0.80.
(a)
(b)Amplitude of the piston by prediction 1 algorithm.
(a)
(b)Amplitude of the piston by prediction 2 algorithm atCFL=0.85 and CFL=0.99999, respectively (a). Zoom to emphasize the damping (b).
(a)
(b)
(c)
(d)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. Figures7 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.Densityρ at Tmax=3×𝒯0, staggered algorithm, via prediction 1 (a) and prediction 2 (b).
(a)
(b)Pressurep at Tmax=3×𝒯0, staggered algorithm, via prediction 1 (a) and prediction 2 (b).
(a)
(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.



---

*Source: 391974-2012-06-20.xml*