A modified fractional step method for fluid–structure interaction problems

Abstract

We propose a Lagrangian fluid formulation particularly suitable for fluid–structure interaction (FSI) simulation involving free-surface flows and light-weight structures. The technique combines the features of fractional step and quasi-incompressible approaches. The fractional momentum equation is modified so as to include an approximation for the current-step pressure using the assumption of quasi-incompressibility. The volumetric term in the tangent matrix is approximated allowing for the element-wise pressure condensation in the prediction step. The modified fractional momentum equation can be readily coupled with a structural code in a partitioned or monolithic fashion. The use of the quasi-incompressible prediction ensures convergent fluid–structure solution even for challenging cases when the densities of the fluid and the structure are similar. Once the prediction was obtained, the pressure Poisson equation and momentum correction equation are solved leading to a truly incompressible solution in the fluid domain except for the boundary where essential pressure boundary condition is prescribed. The paper concludes with two benchmark cases, highlighting the advantages of the method and comparing it with similar approaches proposed formerly.

1

Introduction

Fluid–structure interaction (FSI) problems involving incompressible fluid flows and flexible structures are found in many civil and mechanical engineering applications. Active research has been carried out in the field of FSI over past two decades and multiple numerical models were developed (a review can be found, e.g. in [1]). For the problems involving light-weight structures interacting with free-surface flows quasi-incompressible Lagrangian fluid formulations coupled to the standard structural formulations proved to be very advantageous. The evolution of the free-surfaces and FSI interfaces could be tracked since it was naturally defined by the position of the moving Lagrangian mesh. On the other hand, quasi-incompressible formulations circumvent the added mass effect [2] typically encountered when standard truly incompressible fluid formulations were used [3]. This benefit was achieved due to the relaxation of the incompressibility constraint introduced by the assumption of slight compressibility. Quasi-incompressible fluid formulations have been widely used for FSI simulation both in the finite element method (FEM) [4–7] and the smooth particle hydrodynamics (SPH) contexts [8–11].

For designing monolithic FSI solvers (i.e. the ones that rely on the solution of the coupled problem using a single discrete system) quasi-incompressible formulations are particularly benefitial as they allow for pressure condensation in the fluid domain while maintaining the velocity/pressure coupling. This results in (a) better monolithic system conditioning due to elimination of the different variables scales (b) simplicity of coupling with the structure when both sub-domains are described using the same primary variable (displacement or velocity). In such case elements of the fluid and the structure simply share the same degrees of freedom at a contact node. Thus, a fluid–structure problem can be solved very similarly to a single-material one. Of course, in such approaches the use of fitting interface meshes is obligatory.

Under the assumption of quasi-incompressibility the pressure is related to the kinematic field (velocity or displacement) via a constitutive equation involving the compressibility constant, also called the bulk modulus. For high values of the bulk modulus quasi-incompressible formulations provide an acceptable approximation of the incompressible behavior. The bulk modulus must be sufficiently high to conserve mass in a satisfactory way and introduce the sound propagation speed at least 1 order of magnitude higher than the expected velocity of the bulk flow. For FSI problems the one can update the fluid pressure in the coupling step using the constitutive relation ensuring that the pressure accounts for the motion of the structure. As we shall see further, this pressure update does not involve linear system solution and is therefore computationally cheap.

Quasi-incompressible fluid formulation based on linear velocity-constant pressure finite elements was proposed in [4]. Element-wise constant pressure formulation facilitated pressure condensation at an elemental level, i.e. prior to assembly. This facilitated solving the entire fluid–solid problem using a unified approach with the velocity being the only primary variable. However, the drawback of the formulation was the volumetric locking phenomenon (well-known in constant-pressure elements) that manifested already at moderately high values of the bulk modulus. On the other hand, low values of bulk modulus led to poor approximations of the incompressible behavior. In [6,12] an alternative based on linear pressure interpolations was proposed. The formulation exhibited superior behavior in terms of volumetric locking. Nonetheless, the computational cost of the solver increased due to the impossibility of condensing pressure elementally when using linear pressure approximation. The global pressure condensation procedure had to be introduced. Moreover, when approaching incompressibility limit the pressure instability problems (inf-sup instability [13]) due to using equal order velocity-pressure interpolations manifested. In general, the ambiguity of the quasi-incompressible or penalty approaches can be expressed as follows: the compressibility constant must be large enough to approximate the incompressibility accurately, but at the same time it must be small enough not to lead to “stiff” governing systems. An improvement with respect to modeling the incompressible behavior can be found in [14], where an idea of combining the above-mentioned quasi-incompressible approaches with the fractional step strategy was proposed. The method consisted in using the momentum equation of a quasi-incompressible fluid as a prediction (fractional momentum equation). The subsequent solution of the pressure Poisson's equation and the momentum equation correction led to the truly incompressible solution. The method allowed for using relatively low values of bulk modulus in the quasi-incompressible momentum equation, since the truly incompressible solution was recovered at the correction step. The necessity of the computationally expensive global pressure condensation inherited from [6] due to the use of linear pressure interpolations defined the main drawback of the methodology.

In the present work we propose one further improvement of the methods’ family developed in [4,6,14]. Following the idea of combining the quasi-incompressible prediction with the fractional step method, we propose to use the approximation of the volumetric term in the tangent matrix that allows for computationally efficient elemental pressure condensation, defining a major advantage in comparison with [14]. We also introduce the fluid–structure interaction coupling strategy where the modified fractional momentum equation is solved together with the momentum equation of the structure monolithically, while the subsequent “incompressible correction” steps are carried out in the fluid domain exclusively.

The paper is organized as follows. We first introduce the modified fractional momentum equation using the quasi-incompressibility assumption. An approximate linearization of the volumetric term is introduced. Correction steps ensuring truly incompressible solution are specified next. Then the solution procedure for the FSI problems is outlined. The paper concludes with two challenging FSI benchmark examples.

2

Numerical model

2.1

Governing equations for the fluid

Let us consider a fluid domain Ω with the fixed boundary Γ_d. We shall consider viscous incompressible Newtonian fluids being the most common in the majority of the engineering applications. The governing system are therefore the Navier–Stokes equations equipped with the incompressibility condition. These can be written as:

where v is the velocity vector, p the pressure, t the time, g the body force, ρ the density, μ the dynamic viscosity and ϵ = (∇ v + ∇ ^Tv)/2 – the deviatoric strain rate.

At the fixed wall Γ_d, homogeneous Dirichlet boundary conditions are prescribed:

2.1.1

Finite element formulation

The equal order linear velocity/pressure interpolations over 3-noded triangles (2D) or 4-noded tetrahedra (3D) are used here for the space discretization of the governing equations Eqs. (1) and (2). We assume Backward Euler time integration scheme exclusively for the sake of simplicity. All the arguments presented in the paper are valid for any implicit time integration scheme. In the implementation carried out in this work the second order Newmark–Bossak scheme is used [6]. Being standard, the space and time discretization are not discussed here (see e.g. [15,16]). Pressure stabilization term is added due to the use of the equal order velocity-pressure formulation (Algebraic Sub-Grid Scales (ASGS) stabilization [17] is implemented here). Lagrangian description of the fluid is considered.

Given v¯n and p¯n at t_n, the time discrete problem consists in finding v¯n+1 and p¯n+1 at t_n+1 as the solution of

where M, L, G and S are the mass, the Laplacian, the gradient and the stabilization matrices, respectively. v¯ and p¯ are the velocity and pressure, respectively, and F¯ is the body force vector. Subindices indicate the time step. Note the absence of the convective term due to the use of the Lagrangian kinematic framework.

The matrices and vectors are assembled from the elemental contributions defined as

N stands for the vector of standard linear FE shape functions, Ω_e is the element integration domain, τ is an algorithmic stabilization coefficient defined as τ=(((2||v¯||)/h)+(4ν/h2))−1, where h is the element size. Note also that the discrete operators given by Eqs. (6)–(11) correspond to the unknown current configuration X_n+1 according to updated Lagrangian approach [6,18]. Thus, the system is nonlinear and must be solved in an iterative manner and the discrete operators must be updated at every nonlinear iteration.

Fractional step split. Following the basic idea of the fractional step methods [19–21], the momentum equation is split into two parts by introducing the intermediate velocity v˜

where v˜ is the above-mentioned intermediate or “fractional” velocity and p¯n+1g is the guess of the end-of-step pressure. Eq. (12) is known as “fractional momentum” and Eq. (13) as “end-of-step momentum” equations.

The novelty of these equations in contrast to those of the standard second order fractional step approach consists in using p¯n+1g instead of p¯n in the fractional momentum equation. This idea was originally proposed in [14], where p¯n+1g was computed assuming slight compressibility. This is particularly advantageous for the fluid–structure interaction problems since it implies that the fractional momentum equation becomes equivalent to the momentum equation of a quasi-incompressible fluid, that can be successfully coupled to the momentum equation of the structure according to [6] or [4]. This type of coupling is completely free of the spurious added mass effect since the fluid pressure becomes coupled to the kinematic field (velocity or displacement) of the FSI problem via the pressure constitutive equation according to the quasi-incompressibility assumption. When the classical fractional step method is applied, the fluid pressure becomes completely segregated from the kinematic field in the fluid–structure coupling step [3].

The pressure Poisson's equation is obtained by applying the incompressibility condition Eq. (14) to the end-of-step momentum equation (Eq. (13)), leading to

Using the typical approximation DM⁻¹G ≈ L, we arrive at the final system to be solved:

Fractional momentum equation solution and pressure prediction. The momentum Eq. (16) is nonlinear due to the dependence of the discrete operators on the unknown current configuration X_n+1. Thus, it must be solved iteratively. For this reason, let us define the residual of the fractional momentum equation

According to the assumption of quasi-incompressibility the unknown current-step pressure can be computed by adding the term proportional to the divergence of velocity to the pressure of the previous step (see [6,14] for details):

where K is the bulk modulus of the fluid. The space-discrete form of the constitutive Eq. (20) using linear velocity-pressure finite elements reads

where M_p is the pressure mass matrix.

In order to avoid matrix inversion for obtaining the current step pressure, the pressure mass matrix M_p will be taken in the lumped format. Multiplying Eq. (21) by the inverse of the mass matrix and performing time integration one obtains:

where

Now the unknown pressure increment is expressed in terms of the fractional velocity. Since its computation involves multiplication of assembled (global) matrices M and D, it can be called the “global pressure condensation”.

Expressing the pressure increment as a function of velocity according to Eq. (23) allows us to solve the nonlinear equation

(with r˜m defined by Eq. (19)) exclusively for velocity:

with the subsequent velocity and pressure update: v˜i+1=v¯n+δv¯ and p¯n+1g=p¯n+KΔtMp−1Dv˜i+1, where “i” stands for the nonlinear iteration index at time t_n+1 and H is the tangent matrix defined as

As the non-constant pressure term is now included in the residual, it must be accounted for in the linearization. Using Eq. (22) permits to linearize the pressure gradient term with respect to velocity, giving

thus leading to the following expression for the dynamic tangent:

However, calculation of the volumetric term KΔtGMp−1D in the tangent matrix is a computationally tedious procedure and is feasible only if matrix-free methods are applied [6] for avoiding global matrices’ multiplication and storing the product. Instead, we propose an approximation of this term that corresponds to a element-wise constant pressure approximation. Note that this approximation will be used exclusively in the tangent matrix maintaining linear pressure elsewhere.

Discontinuous approximation for the pressure in the tangent matrix. In case of element-wise constant pressure approximation, the pressure increment (Eq. (20)) can be computed as:

where the operator B and volumetric constitutive matrix C_K are defined (in 2D) as

The pressure gradient can be approximated as

Thus, the linearization of the pressure gradient with respect to velocity can be expressed as:

and the resulting tangent matrix is

Now the linearization of the gradient of pressure increment involves elemental matrices B and C_K exclusively. Thus, using element-wise constant pressure approximation permitted to condense pressure at the elemental level minimizing the associated the computational cost.

Pressure Poisson's equation and the correction step. The next step to be carried out is the correction of the pressure, i.e. obtaining the end-of-step incompressible pressure using Eq. (17). Solution of Eq. (17) requires to impose the pressure boundary conditions due to the presence of the Laplacian L. According to the methodology presented in [14], p¯n+1=p¯n+1g can be used as an essential boundary condition for the pressure necessary for solving the Poisson's equation. The quality of this approximation depends exclusively on the value of the bulk modulus K used in the prediction step. Having the pressure fixed to the predicted value p¯n+1g at the free surface (or at the interface with the structure in the FSI problems), pressure Poisson's equation is solved elsewhere in the domain to give the end-of-step pressure p¯n+1.

This step can be thus viewed as a correction of the predicted pressure p¯n+1g to the correct end-of-step one everywhere except for the free surface (or FSI interface), where the “slightly compressible” pressure is maintained. Consequently, the projection step is carried out according to Eq. (18) and returns the end-of-step divergence-free velocity everywhere in the domain except for the pressure boundary, where the divergence-free velocity is approximated. The implementation procedure of the modified fractional step scheme is summarized in Table 1.

Table 1.

Implementation procedure of the modified fractional step scheme.

1.   Solve the modified fractional momentum Eq. (24) using the iterative procedure according to Eq. (25) and tangent matrix defined by Eq. (34)

•  Update pressure according to Eq. (20) and add it to the fractional momentum residual

•  Compute the new nodal position Xn+1 and move the nodes and update discrete operators

•  Repeat until convergence in terms of velocity is achieved

2.  Solve the pressure Poisson's Eq. (17), using the predicted pressure as the boundary condition

3.  Solve the end-of-step momentum Eq. (18)

3

Coupling with a structure

The present approach can be incorporated into the FSI strategies proposed in our works on quasi-incompressible Lagrangian fluids [6,7]. According to these approaches, a unique discretization is applied to the entire domain containing both the fluid and the structure. A single monolithic FSI system of equations is solved. Thus, the interaction becomes an intrinsic feature of the method and does not involve iterative boundary conditions’ exchange between the fluid and the structure sub-domains. The fractional momentum equation for the fluid and the momentum equation of the structure can be assembled into a single system of equations. This step completely coincides with the procedure proposed in [6,7]. The coupling can be applied to any structure provided that its only nodal degree of freedom is velocity (see e.g. [4] for the velocity-based formulation for solids). Alternatively, the equations for the fluid can be rewritten in terms of displacements instead of velocity in order to facilitate the coupling with displacement-based structural formulations. This ensures the compatibility of the degrees of freedom in the nodes shared by the fluid and the solid elements.

The coupled solution involves a monolithic step where fractional momentum of the fluid and the momentum equation of the structure are solved in a single system of equations and a subsequent correction step carried out in the fluid domain exclusively. The correction step consists in solving the pressure Poisson's equation and end-of-step momentum equation. It ensures that the true incompressibility in the fluid domain is fulfilled.

The solution strategy of the coupled system is summarized in Table 2.

Table 2.

Implementation procedure of the monolithic solution for the FSI problems using the proposed fluid formulation.

1.  Start the nonlinear loop

•  Assemble the monolithic FSI system (consisting of the modified fractional momentum for the fluid and the structural momentum equation) in a standard FE manner

•  Solve the monolithic FSI system for the primary kinematic variable (velocity or displacement)

•  At every nonlinear iteration update the guess for the fluid pressure using the quasi-incompressible fluids’ constitutive relation

•  Repeat until convergence in velocity is achieved

2.  Prescribing the pressure to the predicted value at the FSI interface, solve the pressure Poisson's equation for the end-of-step pressure in the fluid domain.

3.  Perform the correction of the momentum equation in the fluid domain

4.  Go to next time step

4

Examples

In this section two fluid–structure interaction examples are simulated using the formulation proposed. These benchmarks are characterized by the density ratio of the fluid and the structure close to 1, which defines a challenging setting for conventional fluid–structure interaction coupling schemes.

4.1

Deformation of a rubber seal

This example studies the deformation of an elastic plate subjected to water pressure. It was proposed and studied in detail in [8,9]. A water container of width A = 0.1 m with water level L = 0.14 m is closed by a rubber cover of height H = 0.079 m and width s = 0.005 m, which is fixed at the top (see Fig. 1). The cover is released and exposed to the water column, which induces deformation. The rubber cover is modeled with the following properties: density ρ = 1100 kg/m³, Young's modulus E = 6 MPa and Poisson's ratio ν = 0.4

Fig. 1.

Elastic seal subjected to water pressure.

Rubber is a nonlinear material, however in the present study linear elastic approximation is used. The value of the Young's modulus is approximated here as 0.5(E₀ + E₁₀₀), where E₀ is the Young's modulus of the virgin material and E₁₀₀ is the value that corresponds to the deformation of 100% (these values are taken from [9]).

Water is modeled using actual properties of water. The bulk modulus value used in the prediction step was K = 10 KPa. The unstructured uniform mesh with the elemental size of 0.001 m was.

Qualitative comparison among the results of the present simulation and the experimental data published in [8] is shown in Fig. 2. The domain configurations at t = 0.04, 0.12, 0.28 and 0.4 s are displayed. One can see a good agreement between the numerical and experimental data both in terms of the free surface of the fluid and the seal deformation.

Fig. 2.

Elastic seal subjected to water pressure.

Quantitative comparison is shown in Fig. 3. Fig. 3a and b shows the evolution of horizontal and vertical displacements of the middle of the seal tip, respectively. One can see that the maximum horizontal and vertical displacements (around 0.042 and 0.017 m, respectively) are very close to the experimental values. Certain discrepancy is observed once the maximum deformation is obtained, which suggests that the hyper-elastic effects of the real rubber are considerable. However, in the elastic region the numerical and the experimental results match well. Overall, one can see a good agreement among the results obtained with the present simulation and the ones presented in [8]. In particular, our simulation could predict the slight “reopening” of the seal, i.e. the increment of displacements observed around 0.35 s. The difference in the maximum value is less than 10/vertical directions.

Fig. 3.

Displacement of the rubber seal: numerical results vs. experimental data.

4.2

Vertical elastic beam in a shallow oil sloshing

This example was studied in [5], where experimental data and simulation results are provided. It models rotational motion of a rectangular container filled with liquid and a vertical elastic beam clamped at the bottom. The geometry of the model is shown in Fig. 4a. The tank has a length L = 0.609 m and a height H = 0.3445 m. The real set-up has a width of 0.039 m. The container moves around a fixed point located in the mid-point of the bottom wall (x = 0.3045 m, y = 0 m). The motion with an amplitude of 4 degrees and a period of 1.21 s is prescribed to the container walls. The beam is made of polyurethane resin with the following properties: the density is 1100 kg/m³, the Young modulus measured in a traction test amounts approximately 6 MPa. The beam has a thickness b = 0.004 m and width of 0.0332 m which is enough to simulate a 2D flow without touching the lateral walls.

Fig. 4.

Oil sloshing in the container with an elastic beam.

The tank was filled with sunflower oil, with the density of 917 kg/m³ and the kinematic viscosity of 5e−5 m²/s. The original free surface level of the liquid coincided with the beam height (h = 0.1148 m). It is important to note that in the experiment, when the motor is started there is a transition from the rest state to the harmonic motion due to inertia. The numerical simulation accounted for this shift by introducing a delay of 0.25 s in the onset of the tank motion. Uniform unstructured mesh with size of 0.003 m was used.

Fig. 4b displays the evolution of the horizontal displacement d_x of the beam's upper left corner. The results obtained with the present method are compared with the experimental data and the numerical simulation [5,22]. One can see a good agreement with the experimental data and an almost exact match with the numerical results.

Let us examine next the benefits of the methodology proposed in the present paper in comparison with the former Lagrangian FSI formulations based on quasi-incompressibility assumption, such as [4,5] or [6]. The quality of the quasi-incompressible approximations strongly depends upon the value of the bulk modulus. Fig. 5a shows the evolution of the horizontal and vertical displacement of the elastic beam for different values of the bulk modulus K obtained using the quasi-incompressible formulation. One can see that for the value of K = 1000 KPa a perfectly smooth solution is obtained. When K = 100 Pa is used spurious oscillations appear and the solution deteriorates. For K = 10 KPa the spurious oscillations become considerably large and the fluid behavior cannot be considered incompressible anymore. These oscillations represent the compressibility effects that manifest when the bulk modulus is not large enough and the corresponding pressure waves travel with low-enough velocity to create a series of compressions–extensions in the medium.

Fig. 5.

Vertical and horizontal displacements of the elastic beam's tip in sloshing problem.

The present methodology, on the contrary, exhibits perfectly smooth solutions even for small values of the bulk modulus. Fig. 5b shows the solutions obtained for K = 100 KPa and K = 10 KPa. One can see that the graphs practically coincide. No spurious oscillations are observed. This proves that the correction step carried out after the quasi-incompressible prediction considerably improves the solution quality in the fluid domain. Note that even in the limiting case of K = 0 one would simply recover the classical fractional step solution. However, non-zero values of bulk modulus must be used in order to ensure convergent fluid–structure interaction solution. Fig. 6 shows average computational time per time step versus bulk modulus value (the simulations were executed on Intel i7 (4 cores, 2.80 GHz) PC). One can see that for the problem at hand using K = 1000 KPa leads to computational cost practically five times higher than that of K = 10 KPa. This comparison does not pretend to represent any exact estimate of the computational benefit of the proposed method. It merely indicates the gain in the computational efficiency. It is worth noting that in 3D the improvement in the computational efficiency due to the possibility of using low bulk modulus values is even higher.

Fig. 6.

Average computational cost per time step in sloshing simulation for different values of the bulk modulus.

5

Summary and conclusions

This paper presented a fluid formulation that combined the features of fractional step and quasi-incompressible approaches. It consisted of the prediction for the velocity and the pressure under the assumption of quasi-incompressibility and their subsequent correction by means of pressure Poisson's equation and end-of-step momentum equation. The formulation maintained the attractive features of the formerly proposed quasi-incompressible formulations for the fluid–structure interaction coupling. However, it led to truly incompressible solutions even when using low bulk modulus values in the prediction step which defines an important benefit of the proposed approach. Moreover, the proposed approximation of the volumetric term in the tangent matrix permitted elemental pressure condensation at the prediction step, defining an important improvement of our former formulation, where computationally expensive global pressure condensation was mandatory [6].

Using the presented formulation we achieved:

• •

  Truly incompressible solution in the fluid domain

• •

  Possibility of using low values of bulk modulus in the prediction step

• •

  Straight-forward and stable coupling to the structure

• •

  Possibility of using conventional pressure stabilization techniques (in the Poisson's equation)

• •

  Computationally efficient solution strategy

• •

  No added mass effect

The permissible time step in the proposed method is restricted by the non-negativity requirement of the elements’ Jacobian (i.e. no element can be inverted within a time step). This restriction can be possibly alleviated by using innovative streamline time integration techniques proposed in [23]. This defines one clear research line for the future.

Acknowledgement

The author of this work expresses his gratitude to the Spanish Ministerio de Economia y Competitividad for its generous support via Europa Investigacion (EUIN2015-62730) and FPDI-2013-18471 grants.