Numerical Stability of Partitioned Approach in Fluid-Structure Interaction for a Deformable Thin-Walled Vessel

Abstract

Added-mass instability is known to be an important issue in the partitioned approach for fluid-structure interaction (FSI) solvers. Despite the implementation of the implicit approach, convergence of solution can be difficult to achieve. Relaxation may be applied to improve this implicitness of the partitioned algorithm, but this commonly leads to a significant increase in computational time. This is because the critical relaxation factor that allows stability of the coupling tends to be impractically small. In this study, a mathematical analysis for optimizing numerical performance based on different time integration schemes that pertain to both the fluid and solid accelerations is presented. The aim is to determine the most efficient configuration for the FSI architecture. Both theoretical and numerical results suggest that the choice of time integration schemes has a significant influence on the stability of FSI coupling. This concludes that, in addition to material and its geometric properties, the choice of time integration schemes is important in determining the stability of the numerical computation. A proper selection of the associated parameters can improve performance considerably by influencing the condition of coupling stability.

1. Introduction

Fluid-structure interaction (FSI) is used widely in biomechanical computer simulations. The modelling of pulsatile blood flow in elastic vessels requires a framework that can handle the blood-vessel interaction, and the implementation of FSI can solve the time dependent biofluid flow through its elastic structure. Useful information such as the severity of vessel damage by abnormal flow, degree of plaque growth or risk of its rupture in diseased arteries, and the aggravation of atherosclerosis can be generated for medical evaluation. In general, FSI is an architecture processing the interaction of a solid structure with a dynamic fluid that can be implemented by the monolithic and partitioned approaches. A configured FSI solver can use the former approach to solve a system of governing equations for the fluid and solid domains [1]. Although FSI gives strong coupling between the two domains, the limitations are as follows:(i)demanding computational power for solving a large system of equations; (ii)need for further development of preconditioning; (iii)lack of specialized capabilities that pertain to “legacy software” such as ABAQUS and ANSYS.

The partitioned approach is attractive due to its advantage of having software modularity that allows selection of an appropriate solver among the well-established solvers for each of the domains. Nevertheless, its efficiency is inferior to its monolithic counterpart due to the existence of an “added-mass instability,” which commonly occurs in problems involving large deformation and light weight structure. This causes divergence and failure before the final solution can be achieved. In order to handle this issue, small values for coupling relaxation factors of interface loads must be used in order to maintain stability. However, that will lead to a significant increase in computational time. As such, much work has been performed to search for techniques that can deal with this instability. For example, adaptive relaxation techniques that are based on using information of earlier iterations for approximating an appropriate relaxation value and providing stability have been employed to increase the speed of calculation [2]. It is well understood that FSI solution experiences this numerical instability when the following conditions [3] are observed:(i)stability of FSI solution tends to be critically severe when density ratio of fluid to solid is excessively high; (ii)increase in fluid viscosity leads to a decrease in stability of the FSI solver, and a corresponding increase in structural stiffness improves this stability; (iii)temporal discretisation schemes used for FSI calculation can influence the instability condition; (iv)decrease in time step size used for FSI calculation can give an earlier occurrence of its instability.

The observed behaviors of instability can be explained mathematically when conditions for stability of both explicit and implicit FSI solution of a flexible cylindrical vessel are demonstrated [4]. However, the impact of time integration schemes of the fluid and solid accelerations that is used in the FSI calculation is still not fully understood and therefore justifies the need for a thorough investigation in this paper. The aim of this work is to analyze performance of FSI using a partitioned approach based on different time integration schemes that pertain to the structural mechanics. We conduct the analysis on a simplified problem of pressure wave propagation along a flexible cylindrical vessel. Influence of the parameters relating to the FSI performance can be summarized as follows:(i)time integration schemes for solid and fluid accelerations; (ii)time step size; (iii)ratio of fluid to solid densities.

We explore the influence of time integration schemes on a partitioned approach for fluid-structure interaction problems by organising our work in a logical manner. In Section 2, we present a mathematical description of a simplified pressure wave propagation along a flexible cylindrical vessel, fundamental FSI conditions, and the discretization schemes used for deriving the results. Section 3 provides a mathematical analysis of the stability of implicit FSI based on the different time integration schemes that are discussed. In Section 4, numerical validations and results of the same problem (that are used in Section 3) are demonstrated in order to confirm the validity of our theoretical proofs. Finally, Section 5 summarizes the influence of the time integration schemes and the associated parameters on the stability of the FSI.

2. Mathematical Background of Fluid-Structure Interaction

2.1. Governing Equations

A simplified flexible cylindrical tube of radius 𝑟, length 𝐿, thickness ℎ_𝑠, density 𝜌_𝑠, Young’s modulus 𝐸, and Poisson’s coefficient 𝜈 is chosen for our mathematical analysis. It allows our mathematical and numerical examination to be performed and provides sufficient information for conducting a realistic simulation. We define a fluid domain Ω_𝐹 and a structural domain Ω_𝑆 that interact at the common boundary Γ. Deformation of cylindrical tube is allowed only in the radial and longitudinal directions. Inlet and outlet of the fluid domain are subjected to pressure boundary conditions. The schematic representation of our computational domain for a fluid through the tube is shown in Figure 1. Note that the dashed lines represent the axis of symmetry of the tube.

Figure 1 

Schematic representation of the computational domain for a flexible vessel transporting fluid.

2.1.1. Solid Governing Equation

We refer to the governing equation of deformation for a flexible cylindrical tube in [4] as follows: 𝜌_𝑠ℎ_𝑠𝜕²𝑑_Γ𝜕𝑡²+𝑎₀𝑑_Γ𝜕−𝑏²𝑑_Γ𝜕𝑥²=𝑃_Γ, (1) where 𝑎₀=𝐸ℎ_𝑠/𝑟²(1−𝜈²), 𝑏=𝐾_𝑇𝐵ℎ_𝑠. Note that 𝐾_𝑇 denotes the Timoshenko shear correction factor, 𝑃_Γ is structural load on the interface Γ due to the external forcing term from the fluid, 𝑑_Γ is displacement at the interface, 𝑥 is position in space, and 𝑡 is position in time.

2.1.2. Fluid Governing Equations

Considering a general variable property per unit mass that is denoted as 𝜙, the generic form of fluid governing equations in an Arbitrary Lagrangian-Eulerian (ALE) frame of reference is stated as 𝜕𝜌𝜙𝜕𝑡_Temporalterm+∇⋅𝜌𝐯−𝐯_𝑐𝜙_{Convectionterm}−∇⋅𝜌Υ_𝜙∇𝜙_{Diffusionterm}=𝑆_𝜙𝜙_Sourceterm. (2)

In the finite volume method, it is required that the governing equation is satisfied over the control volume 𝐕_𝑃 around a point 𝑃. Therefore, it can be rewritten in the integral form as _𝑡^𝑡+Δ𝑡𝜕𝜕𝑡_𝐕𝑃𝜌𝜙𝑑𝐕_{𝑃Temporalterm}+_𝐕𝑃∇⋅𝜌𝐯−𝐯_𝑐𝜙𝑑𝐕_{𝑃Convectionterm}−_𝐕𝑃∇⋅𝜌Υ_𝜙∇𝜙𝑑𝐕_{𝑃Diffusionterm}=𝑑𝑡_{𝑡𝑡+Δ𝑡𝐕𝑃}𝑆_𝜙𝜙𝑑𝐕_{𝑃Sourceterm}𝑑𝑡. (3)

2.1.3. The Effective FSI Governing Equation

In order to analyze stability of FSI calculation, an effective FSI governing equation can be expressed as 𝜌_𝑠ℎ_𝑠𝜕²𝑑_Γ𝜕𝑡²+𝑎₀𝑑_Γ𝜕−𝑏²𝑑_Γ𝜕𝑥²=𝑃_Γ. (4) The structural load that is influenced by the fluid load 𝑃_ext,Γ is given by 𝑃_Γ=𝑃_ext,Γ−𝜌_𝑓𝑀_𝑎𝜕²𝑑_Γ𝜕𝑡², (5) where 𝑀_𝑎 denotes the added-mass operator matrix.

By coupling the two domains, we can achieve the interaction of the fluid and structure based on 𝜌_𝑠ℎ_𝑠𝜕²𝑑_Γ𝜕𝑡²+𝜌_𝑓𝑀_𝑎𝜕²𝑑_Γ𝜕𝑡²+𝑎₀𝑑_Γ𝜕−𝑏²𝑑_Γ𝜕𝑥²=𝑃_ext,Γ. (6)

Since we are interested in studying the added-mass instability where the mass term dominates the stiffness term, some nonlinearity is neglected. This results in 𝜌_𝑠ℎ_𝑠𝜕²𝑑^𝑠_Γ𝜕𝑡²+𝜌_𝑓𝑀_𝑎𝜕²𝑑^𝑓_Γ𝜕𝑡²+𝑎₀𝑑^𝑠_Γ=𝑃_ext,Γ. (7)

2.2. Fundamental FSI Coupling Conditions

The partitioned approach can enable the physical integration of the fluid and solid domains that is demanded by the simulation of flow through an elastic vessel. The fluid-structure interface is enforced by iteration between the structural and fluid physics modules until convergence is reached. The arbitrary mesh motion at discrete time steps can be computed by the timeintegration of the set of partial-differential equations in Section 2.1 for both the solid and fluid structures that govern the mesh motion.

The conditions required when solving FSI pertain to the kinematic and dynamic nature as suggested by numerous works (see [5] and reference therein). The kinematic condition ensures the compatibility of displacement across FSI interface and can be written as 𝑑^𝑓_Γ=𝑑^𝑠_Γ, (8) where 𝑑 is displacement, 𝑑^𝑠_Γ and 𝑑^𝑓_Γ are the displacements of the solid and fluid interface, respectively.

Assuming that if no-slip condition is used on the fluid side of the FSI interface, this condition leads to a relationship between fluid velocity 𝐯_Γ and rate of change of displacement, which can be written as 𝑛⋅𝑢_Γ=𝑛⋅𝑑̇^𝑠_Γ. (9)

The dynamic condition ensures the compatibility of traction across FSI interface and gives rise to 𝑛⋅𝜏^𝑠_Γ=𝑛⋅𝜏^𝑓_Γ, (10) where 𝜏_Γ represents stress on the interface.

These two conditions are normally utilized in FSI codes that adopt the partitioned approach. By implementation of the kinematic condition, fluid nodes on the FSI interface are updated according to their corresponding solid nodes. By doing the same for the dynamic condition, the equilibrium of stress on FSI interface is ensured and the fluid pressure is integrated into a fluid force, which is used in applying to the solid nodes along the interface.

2.3. Time Integration Schemes

2.3.1. Discretization of Structural Acceleration

The nonlinear version of the generalized-𝛼 time integration scheme that was introduced by [6] can be used for discretization of solid governing equation [7] and is given by 𝑑̇^𝑛+1=𝑑̇^𝑛+(1−𝛿)𝑑̈^𝑛+𝛿𝑑̈^𝑛+1𝑑Δ𝑡,_𝑛+1=𝑑^𝑛+𝑑̇^𝑛1Δ𝑡+2−𝛼𝑑̈^𝑛+𝛼𝑑̈^𝑛+1Δ𝑡², (11) where Δ𝑡 denotes the discrete time step interval.

After some mathematical manipulations, (11) can be recasted into the following equations as: 𝑑̇^𝑛+1=𝑑̇^𝑛+Δ𝑡(1−𝛿)𝑑̈^𝑛+Δ𝑡𝛿𝑑̈^𝑛+1,𝑑̈^𝑛+1=1𝛼Δ𝑡²𝑑^𝑛+1+𝑑^𝑛−1̇𝑑𝛼Δ𝑡^𝑛−1̈𝑑2𝛼−1^𝑛. (12) The value of structural displacement, velocity, and acceleration is interpolated between time level as 𝑑̈^𝑛+1,𝛼=1−𝛼_𝑚𝑑̈^𝑛+1+𝛼_𝑚𝑑̈^𝑛,𝑑̇^𝑛+1,𝛼=1−𝛼_𝑓𝑑̇^𝑛+1+𝛼_𝑓𝑑̇^𝑛,𝑑^𝑛+1,𝛼=1−𝛼_𝑓𝑑^𝑛+1+𝛼_𝑓𝑑^𝑛. (13) In order to maintain unconditional stability and second order accuracy of time integration, the following criteria must be satisfied [6]: 1𝛿≥21,𝛼≥21𝛿,𝛿=2−𝛼_𝑚+𝛼_𝑓,𝛼_𝑚≤𝛼_𝑓≤12. (14) For numerical damping, [8] suggests that these parameters can be written in terms of amplitude decay factor 𝛾 as 1𝛼=4(1+𝛾)²1,𝛿=2𝛼+𝛾,_𝑓=0,𝛼_𝑚=−𝛾. (15)

When minimum numerical damping is applied, the values of 𝛼, 𝛿, 𝛼_𝑚, and 𝛼_𝑓 are set to be 1/4, 1/2, 0, and 0, respectively. These values represent time integration that has zero numerical damping, which can be achieved by setting 𝛾=0. The equations for 𝑑̈^𝑛+1 and 𝑑̇^𝑛+1 can be stated as 𝑑̇^𝑛+1=𝑑̇^𝑛+12Δ𝑡𝑑̈^𝑛+12Δ𝑡𝑑̈^𝑛+1,𝑑̈^𝑛+1=4Δ𝑡²𝑑^𝑛+1+𝑑^𝑛−Δ𝑡𝑑̇^𝑛−Δ𝑡²4𝑑̈^𝑛. (16) After mathematical manipulation, the structural acceleration can be written in terms of deformation at different time levels as 𝑑̈^𝑛+1,𝛼≈1Δ𝑡²4𝑑^𝑛+1−16𝑑^𝑛+32𝑑^𝑛−1−20𝑑^𝑛−2−12̇𝑑Δ𝑡^𝑛−2−𝑑̈^𝑛−2. (17) This time integration scheme cannot be put in terms of displacement only, and it has the fully recursive characteristics, which means that the calculation of time step 𝑛+1 utilizes information of all previous time steps down to the intial step.

When maximum numerical damping is applied, the values of 𝛼, 𝛿, 𝛼_𝑚, and 𝛼_𝑓 are set to be 1, 3/2, −1, and 0, respectively. This can be achieved alternatively by setting 𝛾=1. Therefore, equations for 𝑑̈^𝑛+1 and 𝑑̇^𝑛+1 can be written as 𝑑̇^𝑛+1=𝑑̇^𝑛+12Δ𝑡𝑑̈^𝑛+32Δ𝑡𝑑̈^𝑛+1,𝑑̈^𝑛+1=1Δ𝑡²𝑑^𝑛+1+𝑑^𝑛−Δ𝑡𝑑̇^𝑛−Δ𝑡²2𝑑̈^𝑛. (18) Subsequently, this leads to 𝑑̈^𝑛+1,𝛼=1Δ𝑡²2𝑑^𝑛+1−5𝑑^𝑛+4𝑑^𝑛−1−𝑑^𝑛−2. (19)

2.3.2. Discretization of Fluid Acceleration

Two backward Euler schemes are used for the discretization of fluid acceleration. The first order backward Euler scheme can be written as 𝑢^𝑛+1−𝑢^𝑛=𝑢̇Δ𝑡^𝑛+1, (20) while the second order backward Euler scheme is 𝑢^𝑛+1−𝑢^𝑛=1Δ𝑡3𝑢^𝑛−𝑢^𝑛−1+2Δ𝑡3𝑢̇^𝑛+1. (21)

In this paper, only the zeroth order structural predictor is used in order to estimate the fluid deformation corresponding to structural displacement of previous time step 𝑑^𝑛. This structural predictor can be written as 𝑑^𝑛+1_Γ,𝑃=𝑑^𝑛_Γ. (22)

3. Mathematical Analysis

A stability analysis of FSI calculation based on implicit coupling is presented. By adopting this implicit coupling, several iterations are needed for each time step, and relaxation is applied to maintain the stability of calculation.

3.1. Numerical Procedure of Implicit Coupling

An interface code processes the data transfer between the fluid and solid domains and mesh association across their respective processors. The following steps are used for achieving implicit coupling. Initial guess 𝑑₀^𝑛+1 is given and coupling iteration 𝑚={0,1,2,3,…}.(1)Estimate interface deformation according to previous time or coupling step by using (22), and update the interior fluid mesh. (2)Execute fluid solver, 𝐹(𝑑), such that 𝑃^𝑛+1_𝑚+1=𝐹𝑑^𝑛+1_𝑚+1. (23)(3)Execute solid solver, 𝑆(𝑃), such that 𝑑^𝑛+1_𝑚+1=𝑆𝑃^𝑛+1_𝑚+1. (24)(4)Apply relaxation such that 𝑑^𝑛+1_𝑚+1=(1−𝜔)𝑑^𝑛+1_𝑚+1+𝜔𝑑_𝑚^𝑛+1. (25) Note that the critical relaxation factor is denoted by 𝜔.(5)Check convergence. The solution is converged if the following conditions are applied: 𝑟^𝑛+1,𝑠<𝜀₀,𝑟^𝑛+1,𝑠=𝑑^𝑛+1_𝑚+1−𝑑_𝑚^𝑛+1,𝑟^𝑛+1,𝑓<𝜀₀,𝑟^𝑛+1,𝑓=𝑃^𝑛+1_𝑚+1−𝑃_𝑚^𝑛+1. (26)

3.2. Stability Condition When Minimum Numerical Damping Is Applied

In this section, we discuss the analysis of explicit FSI coupling when minimum numerical damping is implemented. The discretized forms of the structural acceleration and fluid acceleration are given by 𝑑̈_Γ^𝑛+1,𝑠≈1Δ𝑡²4𝑑_Γ^𝑛+1−16𝑑^𝑛_Γ+32𝑑_Γ^𝑛−1−20𝑑_Γ^𝑛−2−12̇𝑑Δ𝑡^Γ𝑛−2−𝑑̈_Γ^𝑛−2,𝑑̈_Γ^𝑛+1,𝑓=1Δ𝑡²𝑑_Γ^𝑛+1−2𝑑^𝑛_Γ+𝑑_Γ^𝑛−1. (27)

Contrary to explicit coupling, the latest information used for calculation of FSI interface acceleration of fluid domain comes from the same time level as the one that is used for the structural acceleration. This is due to the application of relaxation factor. Then, discretization of (7) is achieved by substituting (27) to give 𝜌_𝑠ℎ_𝑠Δ𝑡²4𝑑_Γ^𝑛+1−16𝑑^𝑛_Γ+32𝑑_Γ^𝑛−1−20𝑑_Γ^𝑛−2−12Δ𝑡³𝑑̇_Γ^𝑛−2−1Δ𝑡²𝑑̈_Γ^𝑛−2+𝜌_𝑓𝑀_𝑎Δ𝑡²𝑑_Γ^𝑛+1−2𝑑^𝑛_Γ+𝑑_Γ^𝑛−1+𝑎₀𝑑_Γ^𝑛+1=𝑃^𝑛+1_ext,Γ. (28)

For simplification, the analysis can be done by considering only one eigenvector, 𝐯_𝑖, of solution on the interface. This is because the added-mass operator 𝑀_𝐴 is a real positive matrix. Therefore, 𝑑_Γ=Σ𝑑_𝑖𝐯_𝑖. Notice that the added-mass operator 𝑀_𝐴 can be represented by the 𝑖th eigenvalues of 𝑀_𝐴, 𝜇_𝑖. This results in 𝜌_𝑠ℎ_𝑠Δ𝑡²4𝑑_𝑖^𝑛+1−16𝑑^𝑛_𝑖+32𝑑_𝑖^𝑛−1−20𝑑_𝑖^𝑛−2−12Δ𝑡³𝑑̇_𝑖^𝑛−2−1Δ𝑡²𝑑̈_𝑖^𝑛−2⋯+𝜌_𝑓𝜇_𝑖Δ𝑡²𝑑_𝑖^𝑛+1−2𝑑^𝑛_𝑖+𝑑_𝑖^𝑛−1+𝑎₀𝑑_𝑖^𝑛+1=Δ𝑡²𝜌_𝑠ℎ_𝑠𝑃^𝑛+1_ext,𝑖. (29) By substituting (25) into (29), we obtain 𝜌_𝑠ℎ_𝑠Δ𝑡²4𝑑_𝑖^𝑛+1𝜔−16𝑑^𝑛_𝑖+32𝑑_𝑖^𝑛−1−20𝑑_𝑖^𝑛−2−12Δ𝑡³𝑑̇_𝑖^𝑛−2−1Δ𝑡²𝑑̈_𝑖^𝑛−2+𝑎₀𝑑^𝑛+1_𝑚+1,𝑖𝜔⋯=1−𝜔𝜔2𝜌_𝑠ℎ_𝑠Δ𝑡²+𝑎₀𝑑^𝑛+1_𝑚,𝑖−𝜌_𝑓𝜇_𝑖Δ𝑡²𝑑^𝑛+1_𝑚,𝑖−2𝑑^𝑛_𝑖+𝑑_𝑖^𝑛−1+Δ𝑡²𝜌_𝑠ℎ_𝑠𝑃^𝑛+1_ext,𝑖. (30) By means of Von Neumann stability analysis [9], 1𝜔4𝜌_𝑠ℎ_𝑠Δ𝑡²+𝑎₀𝑑^𝑛+1_𝑚+1,𝑖=1−𝜔𝜔4𝜌_𝑠ℎ_𝑠Δ𝑡²+𝑎₀−𝜌_𝑓𝜇_𝑖Δ𝑡²𝑑^𝑛+1_𝑚,𝑖⋯+𝑔𝑑^𝑛_𝑖,𝑑_𝑖^𝑛−1,𝑑_𝑖^𝑛−2,𝑑̇_𝑖^𝑛−2,𝑑̈_𝑖^𝑛−2,𝑃^𝑛+1_ext,𝑖,𝑑^𝑛+1_𝑚+1,𝑖𝑑^𝑛+1_𝑚,𝑖≈((1−𝜔)/𝜔)4𝜌_𝑠ℎ_𝑠/Δ𝑡²+𝑎₀−𝜌_𝑓𝜇_𝑖/Δ𝑡²(1/𝜔)4𝜌_𝑠ℎ_𝑠/Δ𝑡²+𝑎₀𝜔=1−4𝜌_𝑠ℎ_𝑠+𝑎₀Δ𝑡²+𝜌_𝑓𝜇_𝑖4𝜌_𝑠ℎ_𝑠+𝑎₀Δ𝑡²_{Growthfactor,𝐺}. (31)

The absolute value of the growth factor has to be less than unity for the solution to be stable. That is, |𝐺|=𝜔1−4𝜌_𝑠ℎ_𝑠+𝑎₀Δ𝑡²+𝜌_𝑓𝜇_max4𝜌_𝑠ℎ_𝑠+𝑎₀Δ𝑡²<1,0<𝜔<8𝜌_𝑠ℎ_𝑠+2𝑎₀Δ𝑡²4𝜌_𝑠ℎ_𝑠+𝑎₀Δ𝑡²+𝜌_𝑓𝜇_max. (32)(33)

From (33), we see that not only does the material and geometrical properties (such as fluid and solid densities, Young’s modulus of the structure, and the maximum eigenvalue of added mass matrix) influence the allowable relaxation factor, but they also affect the time step size of the calculation. It can be further deduced that if the time step size is significantly small and approaches zero, (33) becomes 80<𝜔<4+𝜌_𝑓𝜇_max/𝜌_𝑠ℎ_𝑠. (34)

This means that if the time step size is sufficiently small, the Young modulus is no longer a factor that determines the criteria for stability of an implicit algorithm. Therefore, the critical value of relaxation factor converges and does not vary with a further decrease in time step size. Moreover, it can also be concluded that if 𝜌_𝑓𝜇_max=𝜌_𝑠ℎ_𝑠, the relaxation factor has to be strictly less than 8/5 to allow convergence of solution.

3.2.1. Stability Condition When Maximum Numerical Damping Is Applied

The discretized forms of structural acceleration and fluid acceleration are 𝑑̈_Γ^𝑛+1,𝑠=1Δ𝑡²2𝑑_Γ^𝑛+1−5𝑑^𝑛_Γ+4𝑑_Γ^𝑛−1−𝑑_Γ^𝑛−2,𝑑̈_Γ^𝑛+1,𝑓=1Δ𝑡²𝑑_Γ^𝑛+1−2𝑑^𝑛_Γ+𝑑_Γ^𝑛−1. (35) Unlike explicit coupling, the latest information used for calculation of FSI interface acceleration of fluid domain comes from the same time level as that used for structural acceleration. This is due to the application of the relaxation factor. The discretized form of (7) is achieved by substituting (35) and (45) to give 𝜌_𝑠ℎ_𝑠Δ𝑡²2𝑑_Γ^𝑛+1−5𝑑^𝑛_Γ+4𝑑_Γ^𝑛−1−𝑑_Γ^𝑛−2+𝜌_𝑓𝑀_𝑎Δ𝑡²𝑑_Γ^𝑛+1−2𝑑^𝑛_Γ+𝑑_Γ^𝑛−1+𝑎₀𝑑_Γ^𝑛+1=𝑃^𝑛+1_ext,Γ,𝜌_𝑠ℎ_𝑠Δ𝑡²2𝑑_𝑖^𝑛+1−5𝑑^𝑛_𝑖+4𝑑_𝑖^𝑛−1−𝑑_𝑖^𝑛−2+𝜌_𝑓𝜇_𝑖Δ𝑡²𝑑_𝑖^𝑛+1−2𝑑^𝑛_𝑖+𝑑_𝑖^𝑛−1+𝑎₀𝑑_𝑖^𝑛+1=Δ𝑡²𝜌_𝑠ℎ_𝑠𝑃^𝑛+1_ext,𝑖. (36)(37) By substituting (25), (37) can be written as 𝜌_𝑠ℎ_𝑠Δ𝑡²2𝑑^𝑛+1_𝑚+1,𝑖𝜔−5𝑑^𝑛_𝑖+4𝑑_𝑖^𝑛−1−𝑑_𝑖^𝑛−2+𝑎₀𝑑^𝑛+1_𝑚+1,𝑖𝜔=1−𝜔𝜔2𝜌_𝑠ℎ_𝑠Δ𝑡²+𝑎₀𝑑^𝑛+1_𝑚,𝑖−𝜌_𝑓𝜇_𝑖Δ𝑡²𝑑^𝑛+1_𝑚,𝑖−2𝑑^𝑛_𝑖+𝑑_𝑖^𝑛−1+Δ𝑡²𝜌_𝑠ℎ_𝑠𝑃^𝑛+1_ext,𝑖,1𝜔2𝜌_𝑠ℎ_𝑠Δ𝑡²+𝑎₀𝑑^𝑛+1_𝑚+1,𝑖=1−𝜔𝜔2𝜌_𝑠ℎ_𝑠Δ𝑡²+𝑎₀−𝜌_𝑓𝜇_𝑖Δ𝑡²𝑑^𝑛+1_𝑚,𝑖+𝑔𝑑^𝑛_𝑖,𝑑_𝑖^𝑛−1,𝑑_𝑖^𝑛−2,𝑃^𝑛+1_ext,𝑖. (38) By means of Von Neumann stability analysis, 𝑑^𝑛+1_𝑚+1,𝑖𝑑^𝑛+1_𝑚,𝑖≈((1−𝜔)/𝜔)2𝜌_𝑠ℎ_𝑠/Δ𝑡²+𝑎₀−𝜌_𝑓𝜇_𝑖/Δ𝑡²(1/𝜔)2𝜌_𝑠ℎ_𝑠/Δ𝑡²+𝑎₀𝜔=1−2𝜌_𝑠ℎ_𝑠+𝑎₀Δ𝑡²+𝜌_𝑓𝜇_𝑖2𝜌_𝑠ℎ_𝑠+𝑎₀Δ𝑡²_{Growthfactor,𝐺}. (39)

Absolute value of the growth factor 𝐺 has to be less than unity if the solution is to be stable. That is, |𝐺|=𝜔1−2𝜌_𝑠ℎ_𝑠+𝑎₀Δ𝑡²+𝜌_𝑓𝜇_max2𝜌_𝑠ℎ_𝑠+𝑎₀Δ𝑡²<1,0<𝜔<4𝜌_𝑠ℎ_𝑠+2𝑎₀Δ𝑡²2𝜌_𝑠ℎ_𝑠+𝑎₀Δ𝑡²+𝜌_𝑓𝜇_max. (40)(41)

From (41), the material and geometrical properties such as the densities of both the fluid and solid, Young’s modulus, and eigenvalue of added-mass matrix have influence on both the allowable relaxation factor and also the time step size of the calculation. This is similar to the case when the minimum numerical damping is applied. As the time step size approaches zero, (41) becomes 40<𝜔<2+𝜌_𝑓𝜇_max/𝜌_𝑠ℎ_𝑠. (42)

Here, it is proven that if time step size is very small, the Young modulus is no longer a factor determining the criteria for stability of implicit algorithm. This also means that the critical value of relaxation factor converges and does not vary with a further decrease in time step size. Moreover, it can also be concluded that if 𝜌_𝑓𝜇_max=𝜌_𝑠ℎ_𝑠, the relaxation factor has to be strictly less than 4/3 to allow convergence of solution.

4. Numerical Results

4.1. Validation of Computational Model

Pressure pulse velocity has been widely used as an indicator of blood vessel health. In the past, only physical experiment and analytical solutions were available for obtaining the pulse velocity. Although these analytical solutions can only be used for simple geometries of blood vessel, they can serve as solid verification for numerical simulation, which are becoming popular nowadays. Therefore, a simulation of fluid-structure interaction for calculating pressure pulse velocity in a compliant vessel is conducted to verify our proposed technique in order to compare with analytical solutions such as the Moens-Korteweg equation. The results presented in this section are obtained by using the same method as that used in the next section where numerical results are presented to confirm our theoretical results.

The geometrical model of a three-dimensional tube is shown in Figure 1. The cylindrical domain has a radius of 𝑅_𝑠=0.0005 m and a total length of 𝐿=0.06 m. The working fluid is modelled as incompressible fluid with fluid viscosity and density of 𝜇=0.01 Pa·s and 𝜌_𝑓 = 1,000 kg·m³, respectively. On the other hand, the compliant vessel is modelled as isotropic material with Poisson’s ratio of 𝜈=0.3 and density of 𝜌_𝑠 = 1,000 kg·m³, respectively. Young’s modulus of the solid structure is varied between 𝐸=50 and 200 kPa.

Figure 2 shows the computational mesh of both solid and fluid domains that consist of 1,800 and 24,674 elements, respectively. In the fluid mesh, the laminar boundary layer has 5 layers of thin hexagonal elements, while 4 layers are used for the thickness of solid domain. The mesh in the fluid domain is Geometric Conservation Law (GCL) compliant and the architecture generates arbitrary mesh motion at discrete time steps. Now, our main aim is to validate the theoretical results.

Figure 2 

Computational meshes of the fluid and solid domains for a flexible vessel.

Then, analytical sets of equations developed for pressure wave velocity are discussed next. The Moens-Korteweg equation [10, 11] is 𝑐=𝐸ℎ_𝑠2𝑅_𝑠𝜌_𝑠, (43) where 𝑐 represents the pressure wave velocity.

The Wylie and Streeter equation [12] is 𝑐=𝜑𝑐^∗_𝑓, (44) where 𝜑 is the stress factor and 𝑐^∗_𝑓 is the corrected pressure wave velocity as given by 𝜑=1−𝜈²𝐸1−𝐾ℎ_𝑠2𝑟1−𝐾𝜌_𝑓𝐸𝜌_𝑠⁻¹,𝑐^∗_𝑓=𝐸𝜌_𝑓2𝑅_𝑠+2ℎ_𝑠²ℎ_𝑠2𝑟+ℎ_𝑠−2(1−𝜈)⁻¹. (45)

Note that 𝐾 is the bulk modulus of elasticity of the tube walls. These two equations will be utilised for the validation of numerical method used in this work. The relationship between pressure wave velocity and Young’s modulus of flexible vessel illustrates adequate accuracy of the method (Figure 3).

Figure 3 

Relationship between pressure wave velocity and Young modulus of a flexible vessel.

The pressure pulse velocity is calculated by measuring the location of maximum deformation at two different times based on the sixth and eleventh millisecond. Distance between the two locations is divided by a time difference of five milliseconds in order to obtain the velocity of the pressure pulse. This location is used for the measurement as it represents the centre of the pressure wave. Figure 4 shows the contour plot of pressure wave propagation along the flexible cylindrical vessel at different Young’s moduli and simulation time levels.

Figure 4 

Pressure wave propagation along a flexible cylindrical vessel at various Young moduli and time levels.

Next, we present a set of numerical results with different structural time integration scheme by varying the amplitude decay factor 𝛾. These results can be used to confirm our theoretical experiments. By our default configuration, which is used to as a reference to compare with other cases, we set vessel radius 𝑅_𝑠=0.005 m, vessel length 𝐿=0.06 m, vessel thickness ℎ_𝑠=0.001 m, solid density 𝜌_𝑠 = 30,000 kg·m³, Young’s modulus 𝐸 = 750,000 Pa, Poisson’s ratio 𝜈=0.5, fluid density 𝜌_𝑓 = 1,000 kg·m³, and fluid dynamic viscosity 𝜇_𝑓=0.01 Pa·s. At the initial condition, the fluid is assumed to be at rest and a pressure pulse with peak of 2,000 Pa is imposed at the inlet. The total duration of our observation 𝑇_max for instability is 0.02 s. For majority of the test cases, the instability occurs at the beginning of calculation before reaching 𝑇_max. In our program, the fluid acceleration is discretized by the second order accurate backward Euler scheme.

4.2. Influence of Simulation Parameters on 𝛾 versus 𝜔 Curve

We observe a dependence of the critical relaxation factor 𝜔 on the amplitude decay factor 𝛾. In particular, the required relaxation factors for maintaining stability become greater as the amplitude decay factor approaches zero and the structural time integration scheme becomes the average acceleration method. The numerical tests agree well with our theoretical results, and they confirm that the critical relaxation factor increases with the decrease in the structural amplitude decay factor. Moreover, it is observed that the required values of relaxation factor when using amplitude decay factor of zero can be approximately 30 percent higher than that when using amplitude decay factor of one as shown in Figure 5. By varying parameters and configurations of the FSI architecture such as (a) fluid time integration scheme, (b) time step size, and (c) fluid-solid density ratio, we can obtain information on their effect on the graph of critical relaxation factor versus amplitude decay factor. For all the graphs, it can be demonstrated that the value of critical relaxation factor is inversely proportional to the value of amplitude decay factor.

(a)

(b)

(c)

Figure 5 

Influence of amplitude decay factor on the critical relaxation factor based on (a) fluid time integration scheme; (b) time step size; and (c) fluid-solid density ratio.

First, we prove numerically the influence of the fluid time integration scheme on the FSI stability condition. By changing the discretization scheme used for fluid acceleration from the second to the first order accurate backward Euler scheme, the values of critical relaxation factor based on each value of 𝛾 increase considerably (Figure 5(a)). With 𝛾=1, the solution will be stable if the relaxation factor is set at lower than a conservative value of 0.7, while this threshold increases up to 0.95 when 𝛾 is set to zero.

Then, the influence of time step size on the stability of FSI calculation is tested. Here, we consider the same domain of calculation and physical parameter as before, and only time step size used for calculation is reduced from 0.001 to 0.0001 s. It is found that time step size also has considerable impact on the performance of FSI simulation. From Figure 5(b), it is shown that the values of critical relaxation factor corresponding with various values of amplitude decay factor decrease considerably when time step size is reduced. Moreover, when the time step size is further reduced to 0.00001 s, the change in critical relaxation factor remains almost unchanged. This observation agrees well with our theoretical results that when the time step size approaches zero, the value of critical relaxation factor converges as presented in Figure 5(b).

The influence of the fluid and solid structure densities on the performance of FSI calculation is considered next. The relationship between fluid-solid density ratio and the critical relaxation factor is based on different values of amplitude decay factor. Figure 5(c) illustrates the impact that the density ratio has on variation of critical relaxation factor, which is observed to be high only when the density ratio is relatively small. Its gradient tends to vanish as we increase the density ratio value. For fluid-solid density ratio of 0.033, the variation of the critical relaxation factor in the range of 𝛾=0 to 1 is considerably high, while it is negligible at a density ratio of 1.

5. Conclusion

We summarize that the stability condition of FSI solution can be influenced significantly by the material densities and the relaxation factor to be implemented. Its computational cost can be greatly reduced by implementation of an appropriate structural time integration scheme. In particular, the critical relaxation factor is higher when using the structural time integration scheme that does not contain numerical damping. Furthermore, the choice of time integration schemes for discretizing fluid acceleration has a strong impact on the stability condition. Our results suggest that more accurate schemes such as the second order accurate backward Euler can lead to a more stringent condition for the stability of FSI. Another important parameter is the time step size. Then, smaller time step size results in a more restrictive condition, and as the time step size approaches zero, the value of the critical relaxation factor converges to a specific value. Another factor worth mentioning is the density ratio of fluid to solid structure. It is found that this value has a considerable impact on FSI performance, and the stability condition is invariant to the choice of structural time integration schemes in the case of high fluid-solid density ratios. However, for the low density ratios, this value will mainly depend on the schemes used.