key: cord-0045960-akvs0nlk
authors: Wang, Yongxing; Jimack, Peter K.; Walkley, Mark A.
title: A Three-Dimensional, One-Field, Fictitious Domain Method for Fluid-Structure Interactions
date: 2020-05-25
journal: Computational Science - ICCS 2020
DOI: 10.1007/978-3-030-50436-6_3
sha: d061c15836e80dc044d79f8ed8577e1e9a99aeab
doc_id: 45960
cord_uid: akvs0nlk

In this article we consider the three-dimensional numerical simulation of Fluid-Structure Interaction (FSI) problems involving large solid deformations. The one-field Fictitious Domain Method (FDM) is introduced in the framework of an operator splitting scheme. Three-dimensional numerical examples are presented in order to validate the proposed approach: demonstrating energy stability and mesh convergence; and extending two dimensional benchmarks from the FSI literature. New three dimensional benchmarks are also proposed.

Numerical simulation of Fluid-Structure Interaction (FSI) problems is a computational challenge due to its strong nonlinearity, especially in the case of large solid deformations. This challenge is exacerbated in three dimensions due to the need for efficient numerical algorithms to handle the large number of degrees of freedom that are inevitably required. In this paper we generalize our recent onefield Fictitious Domain Method (FDM) [22, 23] from two to three dimensions, enhance the efficiency and robustness of the proposed time-stepping scheme, and demonstrate the resulting algorithm's capabilities on a number of challenging test problems. We also provide potential benchmark problems to allow results to be compared against those from other schemes in the future.

Lagrangian and Arbitrary Lagrangian-Eulerian (ALE) methods are widely adopted when considering a relatively small solid deformation [6, 11, 15] . Discrete remeshing can be used for large deformations [10, 18] , however this can be very costly in the case of three dimensions and can present challenges for mass conservation. The cut finite element method (cutFEM) [4, 5] may also be applied to solve FSI problems [14, 19] , although it is not trivial to deal with the discontinuous integral across the elements cut by the moving fluid-solid interface, especially in three dimensional cases. The Fictitious Domain Method (FDM) [1, 3, 8, 12, 13] uses two meshes to represent the fluid and solid separately, which can easily handle large deformation of the solid. However the FDM approach solves a very large equation system: both the velocity in the whole domain (fluid and solid) and the displacement in the solid domain, coupled via a distributed Lagrange multiplier (DLM) which is also an unknown variable. The Immersed Finite Element Method (IFEM) [2, 17, 24, 27] also uses two meshes but only solves for velocity in the whole domain, while the solid information is assembled on the right-hand side of the fluid equation as a prescribed force term. This IFEM approach achieves FSI behaviour through this forcing term, and is therefore relatively efficient in three dimensional simulations compared to the DLM approach. Performance of the IFEM method depends strongly on the fluid and solid properties and usually works well when the solid behaves similarly to the fluid (such as a relatively soft solid) [20] . It has been successfully used, for example, in the area of biomechanics [16, 27] .

The one-field FDM approach [22, 23] similarly only solves for one velocity field in the whole domain. However, this proposed method assembles the solid equations and implicitly includes them in the equation system. The one-field FDM approach has the same generality and robustness as the FDM/DLM: both of them solve the fluid equations and solid equations as one system. However the former needs to solve only for one velocity field while the latter solves for fluid velocity, solid displacement and Lagrange multiplier. The proposed onefield FDM may also be regarded as a special linearisation of the implicit IFEM, which however is more robust compared with explicit IFEM and more efficient compared with the implicit IFEM [24] , allowing a wide range of solid parameters to be considered and naturally dealing with the case of different densities between fluid and solid [22, 23] . In short, the one-field FDM combines the FDM/DLM advantage of robustness and the classical/explicit IFEM advantage of efficiency. The scheme has been validated through comparison with idealised two-dimensional test cases and against experimental data and simulation results drawn from the literature [21] .

In this article, the one-field FDM is extended, implemented and validated in three dimensions for the first time through the use of a newly applied operator splitting scheme. The paper is organized as follows. The control equations and a general finite element weak formulation are introduced in Sect. 2.1 and 2.2 respectively, followed by time discretization in Sect. 2.3. The operator splitting scheme is introduced in Sect. 2.4, followed by the linearisation (implementation detail) in Sect. 2.5 and the numerical algorithm for the final linear equation system in Sect. 2.6. Several three-dimensional numerical tests are given in Sect. 3, and conclusions are presented in Sect. 4.

In this section, we review the one-field fictitious domain method [22] and develop it further based upon a three-step operator splitting scheme and the case novel block-matrix preconditioners. The system is described in a manner that is independent of the spatial dimensions, thus ensuring its capability in three dimensions, which is the primary purpose of this paper. 

In the following context, Ω f t ⊂ R d and Ω s t ⊂ R d (d = 3 in this article) denote the moving fluid and solid domain respectively, with the moving interface Γ t = Ω f t ∩ Ω s t as shown schematically in Fig. 1 .

with Γ D and Γ N being Dirichlet and Neumann boundaries respectively. We denote by X the reference coordinates of the solid, by x = x(·, t) the current coordinates of the solid, and by x 0 the initial coordinates of the solid. Notice that we choose X to be the stress-free configuration, which may be different to the initial configuration x 0 .

Let ρ, μ, u, σ and g denote the density, viscosity, velocity, stress tensor and acceleration due to gravity respectively. We assume both an incompressible fluid in Ω f t and incompressible solid in Ω s t . The conservation of momentum and conservation of mass therefore take the same form as follows. Momentum equation:

and continuity equation:

An incompressible Newtonian constitutive equation in Ω f t can be expressed as:

with τ f = μ f Du f being the deviatoric part of stress σ f , and Du = ∇u + ∇ T u. An incompressible neo-Hookean solid with viscosity μ s is assumed in Ω s t [3] , and the constitutive equation may be expressed as:

with τ s = c 1 FF T − I + μ s Du s being the deviatoric part of stress σ s , and F = ∂x ∂X = ∂x ∂x0 ∂x0 ∂X = ∇ 0 x∇ X x 0 being the deformation tensor of the solid, and c 1 is a solid material parameter. Finally the system is completed with continuity of velocity u f = u s and normal stress σ f n s = σ s n s on interface Γ t , and standard Dirichlet/Neumann boundary (on Γ D /Γ N ) and initial conditions.

In the following context, let L 2 (ω) be the square integrable functions in domain ω, and H 1 (ω) = u : u, ∂u ∂xi ∈ L 2 (ω) for i = 1, · · · , d . We also denote by H 1 0 (ω) the subspace of H 1 (ω) whose functions have zero values on the Dirichlet boundary of ω, and denote by L 2 0 (ω) the subspace of L 2 (ω) whose functions have zero mean value.

Given v ∈ H 1 0 (Ω) d , we perform the following symbolic operations:

Integrating the stress terms by parts, the above operations, using constitutive equations (3) and (4), give:

whereh denotes the prescribed normal stress on Γ N . Note that the integrals on the interface Γ t are cancelled out due to the continuity of normal stress: σ f n s = σ s n s , because they are internal forces for the whole FSI system. Combining with the following symbolic operations for q ∈ L 2 (Ω),

leads to the weak form of the FSI system as follows.

, the following equation holds:

where ρ δ = ρ s − ρ f and μ δ = μ s − μ f , and the integral over Ω s t , d(·) dt is the time derivative with respect to a frame moving with the solid velocity u s = u| Ω s t .

Using the backward Euler method to discretise in time, Eq. (6) may be approximated as follows.

Given u n , p n and Ω s n , find u n+1 ∈ H 1 (Ω) d , p n+1 ∈ L 2 (Ω) and Ω s n+1 , such that for ∀v ∈ H 1 0 (Ω) d , ∀q ∈ L 2 (Ω), the following equation holds:

and Ω s n+1 is updated from Ω s n by the following formula:

The formulation of (7) is implicit. However we shall solve it semi-implicitly via the following operator spitting scheme which is based upon [9] .

(1) Convection step:

(2) Diffusion step:

(3) Pressure step:

In the above, ∇ n (·) represents the divergence in the current coordinates at t = t n and D n = ∇ n + ∇ T n . Note that the variables u n+1/3 and u n+2/3 are just intermediate values, not specifically the velocity at time t = t n + Δt 3 or t = t n + 2Δt 3 . The notation F n+1/3 or F n+2/3 is interpreted as follows:

with t = n + 1/3 or n + 2/3. Using this splitting scheme, standard approaches can be taken to solve the pure convection equation (9) (see [9] ), and iterative methods with an efficient preconditioner can be applied to solve the "degenerate" Stokes Equations (11) (see Sect. 2.6 of [7, 21] ). The main challenge is in how to approximate the term F n+2/3 F T n+2/3 − I in Eq. (10), which is nonlinearly related to the solid displacement and hence to the solid velocity. In the following subsection we focus on expressing and linearising F n+2/3 F T n+2/3 − I in terms of velocity u n .

The specific choice of linearisation is the core of this proposed one-field FDM approach, and is what makes it distinctive from all other schemes. Let F n+2/3 F T n+2/3 − I be denoted by F t F T t − I = s t with t = n + 2/3, then s t can be computed as follows:

Using the chain rule, this last equation can also be expressed as:

or

Then s t can be expressed based on the previous coordinate x n as follows:

Using x t = x n + Δtu t (see (12) ), this can finally be expressed as:

There are two nonlinear terms in this equation, which may be linearised as

and ∇ n u t s n ∇ T n u t = ∇ n u t s n ∇ T n u n + ∇ n u n s n ∇ T n u t − ∇ n u n s n ∇ T n u n .

Substituting s n+2/3 = F n+2/3 F T n+2/3 − I, using expressions (17), (18) and (19) , into the diffusion step (10) , and neglecting terms of O Δt 2 , after some algebra we produce the one-field FDM formulation:

In this section, we shall discuss the numerical algorithms in order to solve the final linear equations from the diffusion step (20) and pressure step (11) . For convection step, we use the Taylor-Galerkin method in this paper [9] . Let us write Eq. (20) in an operator matrix form as follows:

where A n = M/Δt + K + P T n (M s n /Δt + K s n ) P n , and b n = f + P T n f s n + Mu n+1/3 /Δt + P T n M s n P n u n /Δt. The above matrix operators are defined as:

We use the finite element interpolation to approximate P n after discretisation in space.

A preconditioned Conjugate Gradient method can efficiently solve Eq. (21) . We use the incomplete Cholesky decomposition of matrix M/Δt + K as a preconditioner in order to solve Eq. (21) . Very good convergence performance can be observed from our numerical tests (although the precise performance of the linear algebraic solver is not the topic of this article).

Similarly, the "degenerate Stokes" problem (pressure step (11) ) can also be expressed in an operator matrix form:

where, ∀v in H 1 (Ω) d and q ∈ L 2 (Ω), Ω (Bv) q = − Ω q∇ · v. We use the MinRes algorithm [7] to solve the system with the following preconditioner:

We justify this since we can derive a Schur complement in the form of S = B T M −1 B. The operators that are discretised in this form imply that S will be spectrally equivalent to a discrete Laplacian. Hence we expect this preconditioner will be effective for this system, similarly to analysis for Stokes equation [7] .

In the following numerical tests, the convection and diffusion steps are discretised with quadratic finite elements (tri-quadratic hexahedra and quadratic tetrahedra), and the pressure step is discretised with the Taylor-Hood element. For stability it is sufficient that μ δ ≥ 0 [21, 22] , however for simplicity, and to be consistent with [2, 17, 27] for example, we assume μ δ = 0 in these tests.

In this section, we consider a 3D oscillating ball which is an extension of the 2D disc in [23, 28] . We use this example to test stability of the proposed approach by investigating the evolution of total energy:

The four different energy contributions/terms in the above equation The initial velocities of x and y components are the same as that used in [22, 28] , which are prescribed by the stream function Φ = Φ 0 sin(ax)sin(by), with Φ 0 = 5.0×10 −2 and a = b = 2π. The z component of velocity is initially set to be 0. In this test, ρ f = 1, μ f = 0.01, ρ s = 1.5 and c 1 = 1. In order to visualise the mesh and deformation of the solid, a snapshot of fluid velocity and pressure are presented in Fig. 2 , and the corresponding deformed solid is displayed in Fig. 3(a) . It can be seen from Fig. 3(b) that the total energy is nonincreasing, which is an indication of stability. In addition, the total energy converges to the initial system energy as we reduce the size of the time step, which shows the desired energy conservation property of the proposed scheme.

In this test we consider a cylindrical pillar oscillating in a cuboid channel as shown in Fig. 4 , which is a 3D extension of the 2D leaflet in [13, 22, 25] . We use this example to test the mesh convergence of the proposed scheme. The size of the cuboid is: length L = 3, height H = 1 and width W = 0.5. The cylinder is located at the center of the cuboid's base, with radius of r = 0.05 and height h = 0.8. We use a symmetry boundary condition on the top, front and back Table 1 . Material properties for the oscillating cylinder and oscillating tri-leaflets.

Fluid Solid ρ f = 100 kg/m 3 ρ s = 100 kg/m 3 μ f = 10 N · s/m 2 c1 = 10 7 N/m 2 surfaces of the cuboid. All the velocity components are fixed to be zero at the bottom of the cuboid, and the inlet and outlet flow are defined by:

We use the same material properties as used in [13, 22, 25] for the 2D leaflet (see Table 1 ), which is a natural extension of the corresponding 2D problem with similar boundary conditions. We use a tri-quadratic hexahedras fluid mesh of size 10×20×60 (width×height×length) for a coarse mesh, 16×32×96 for a medium mesh and 20 × 40 × 120 for a fine mesh. We use a linear tetrahedral solid mesh of 10304 elements with 2675 vertices for a coarse mesh, 19040 elements with 4786 vertices for a medium mesh and 38080 elements with 8883 vertices for a fine mesh. A stable small time step Δt = 1.0 × 10 −4 is adopted for all the cases. In order to visualise the results of this simulation, snapshots of the velocity norm and stream lines in the background domain and the solid deformations are presented in Fig. 5 and 6 respectively. The displacement of initial point (1.55, 0.8, 0.5) (the top of the cylinder) for three different meshes is plotted in Fig. 7 as a function of time, from which mesh convergence with regard to the displacement is observed (the medium and fine mesh results are almost indistinguishable in these plots). 

In this section we consider a 3D circular tube with flexible, opening tri-leaflets. A similar case has been studied in [26] . The computational domain is shown in Fig. 8 with L = 2 and R = 0.5 in this test. Note that there is a small gap (with the angle α = 0.4 • as shown in Fig. 8(b) ) between the three parts of the tri-leaflets in order to avoid contact, which is not currently included in our which is an extension of formula (26) . The material properties for the fluid and solid are presented in Table 1 . Two different meshes are used to test this problem: a coarse mesh of 12750 tri-quadratic hexahedra with 106151 nodes and a fine mesh of 90000 tri-quadratic hexahedra with 735861 nodes in order to discretise the cube; a coarse mesh of 7390 linear tetrahedra with 2657 vertices and a fine mesh of 27460 linear tetrahedra with 8917 vertices in order to discretise the trileaflets. The density of the background mesh can be observed in Fig. 9 , which also presents a snapshot of the velocity norm and stream lines. The solid mesh can be observed in Fig. 11 , which also shows the deformation of the tri-leaflets with horizontal velocity (x component) at different stages in order to visualise the pattern of the oscillation. The displacement at one of the tri-leaflet tips is plotted as a function of time in Fig. 10 , from which it can be seen that the coarse mesh leads to small oscillation, however it is not present in the fine mesh simulation. The maximal fluid velocity is u x = 18.2 at the centre of the channel when the tri-leaflets are completely open, and the maximal solid velocity at the leaflet tip is u x = 18.2 when the tri-leaflets are completely close. 

In this article the one-field Fictitious Domain Method (FDM) [22] is extended in three ways: through an efficient operator splitting scheme, the implementation of block-matrix preconditioning and into three space dimensions. One numerical example is presented in order to validate the energy conservation, a second is used to test mesh convergence, and the last numerical example is extended from a two-dimensional benchmark for comparison and, we believe, can act as a 3D benchmark for future comparison. It can be seen from these tests that the onefield FDM approach may be adopted to simulate a variety of FSI problems with large solid deformation in three space dimensions. We know from our 2D tests that, for soft solids, execution times are comparable with an IFEM implementation: and considerably faster as the solid becomes more stiff (capable of using a larger time step). Consequently we propose the one-field FDM as a general approach that combines that robustness of FDM/DLM and the computational efficiency of IFEM.

A fictitious domain/mortar element method for fluid-structure interaction

The finite element immersed boundary method with distributed Lagrange multiplier

A fictitious domain approach with Lagrange multiplier for fluid-structure interactions

Fictitious domain methods using cut elements: III. A stabilized Nitsche method for Stokes problem

Cut finite element methods for partial differential equations on embedded manifolds of arbitrary codimensions

Performance of a new partitioned procedure versus a monolithic procedure in fluid-structure interaction

Finite Elements and Fast Iterative Solvers

A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: application to particulate flow

Finite Element Methods for Incompressible Viscous Flow

An energy stable monolithic Eulerian fluid-structure finite element method

Solvers for large-displacement fluid-structure interaction problems: segregated versus monolithic approaches

A mortar approach for fluid-structure interaction problems: immersed strategies for deformable and rigid bodies

A fictitious domain/distributed Lagrange multiplier based fluid-structure interaction scheme with hierarchical B-Spline grids

A stabilised immersed framework on hierarchical b-spline grids for fluid-flexible structure interaction with solid-solid contact

Fixed-point fluid-structure interaction solvers with dynamic relaxation

Immersed finite element method and its applications to biological systems

The immersed boundary method

Numerical study of a monolithic fluid-structure formulation

Monolithic cut finite element based approaches for fluid-structure interaction

Interpolation functions in the immersed boundary and finite element methods

A one-field fictitious domain method for fluid-structure interactions

A one-field monolithic fictitious domain method for fluid-structure interactions

Energy analysis for the one-field fictitious domain method for fluid-structure interactions

A theoretical and numerical investigation of a family of immersed finite element methods

A DLM/FD method for fluid/flexible-body interactions

A three-dimensional fictitious domain method for the simulation of fluid-structure interactions

Immersed finite element method

A fixed-mesh method for incompressible flowstructure systems with finite solid deformations