Abstract

Some numerical approaches to solve fluid structure interaction problems in blood flow are reviewed. Fluid structure interaction is the interaction between a deformable structure with either an internal or external flow. A discussion on why the compliant artery associated with fluid structure interaction should be taken into consideration in favor of the rigid wall model being included. However, only the Newtonian model of blood is assumed, while various structure models which include, amongst others, generalized string models and linearly viscoelastic Koiter shell model that give a more realistic representation of the vessel walls compared to the rigid structure are presented. Since there exists a strong added mass effect due to the comparable densities of blood and the vessel wall, the numerical approaches to overcome the added mass effect are discussed according to the partitioned and monolithic classifications, where the deficiencies of each approach are highlighted. Improved numerical methods which are more stable and offer less computational cost such as the semi-implicit, kinematic splitting, and the geometrical multiscale approach are summarized, and, finally, an appropriate structure and numerical scheme to tackle fluid structure interaction problems are proposed.

1. Introduction

Fluid structure interaction is defined as the interaction between deformable structures with an internal or surrounding fluid flow. Such deformation can either be stable or oscillatory. Problems involving fluid structure interaction are classified into the one-way problem which occurs when the movement of the structure controls the motion of fluid but the fluid’s motion does not influence the structure, or the two-way fluid structure interaction problem when the movement of the structure influences the motion of the fluid and vice versa [1].

Fluid structure interaction is more often considered in modelling biofluids because the interaction between the blood and vessel wall is of great clinical interest, for example, in studying cardiovascular diseases which are a major cause of death in developed countries [2].

The interaction between blood flow and vessel wall is often neglected because the coupled fluid and solid equations are complicated and difficult to solve [3]. Earlier numerical models used to predict blood flow are based on rigid geometries [4] in which only the arterial lumen needs to be reconstructed and discretized, yielding results that are reasonably accurate and can be obtained in a relatively short time [5]. However, there are still further considerations to be taken into account such as the elastic nature and stresses on the arterial wall that play crucial roles in arterial disease, as well as the material property alterations with the development of the atherosclerotic lesion [6].

Numerous studies had been carried out to compare the effect of the rigid and compliant wall on blood flow. Lee and Xu [7] indicated that the axial velocities at the center being of a rigid wall are higher compared to the ones in the compliant model. Mass conservation theory is utilized to explain such phenomena as the internal fluid pressure exerted on the vessel wall pushes the vessel wall outward consistently and slows the fluid flow due to the flow area expansion. Rigid wall simulation of blood flow through arteries also overpredicts the wall shears stress. These findings showed that incorporating fluid structure interaction has significant effects on blood flow characteristics [2, 7]. Siogkas et al. [5] considered a rigid wall assumption and fluid structure interaction simulations of blood flow in the arterial segments. They concluded that the computational time for the simulation of fluid structure interaction is longer than that of the rigid wall assumption. It was found that the required time for the simulation is 30–40 minutes in the case of rigid wall but takes almost 5 hours in the compliant vessel.

Fluid structure interaction describes the wave propagation in arteries driven by the pulsatile blood flow. From the theoretical point of view, such problems are complex and challenging due to the high nonlinearity of the problem. Not only the fluid equation exhibits nonlinearity, the structure displacement modifies the fluid domain which generates geometrical nonlinearities as well [8].

The generalized string model had been utilized as the structure of blood flow in compliant vessels and arteries [915]. Causin et al. [16] explained that the generalized string model is a structural model derived from the theory of linear elasticity for a cylindrical tube with small thickness. The reference configuration is a cylindrical surface of the base circle radius that is supposed to move in radially, neglecting the longitudinal and angular displacements. Nobile and Vergara [15] pointed out that the generalized string model neglects bending as well. According to Čanić et al. [1720], there are no analytical results which are able to prove the well posedness of fluid structure interaction problems without assuming the structure model that includes the higher order derivative terms, capturing the viscoelastic behavior, or the terms describing bending rigidity. In hemodynamics, there exists a strong added mass effect issue in which the fluid and structure have comparable densities. If the structure density is higher than the fluid density, such as in aeroelasticity, the added mass effect is negligible. Various structural models are discussed in Section 2.

Numerical approaches of fluid structure interaction which are discussed in Section 3 can be broadly classified into two: the partitioned approach and the monolithic approach. Partitioned approach can be further subdivided into the loosely and strongly coupled algorithms [2027]. In hemodynamics, the use of explicit partitioned algorithm turns out to be problematic where stability is concerned, particularly due to the added mass effect. In addition, the implicit partitioned algorithms are also affected by the added mass effect in terms of convergence. Special treatment of the interface conditions needs to be considered [8, 16, 2830]. To date, it seems that only the monolithic and implicit schemes are applicable in blood flow simulation involving fluid structure interaction. However, they are costly in terms of computational cost, computational time, and memory requirement [911, 3133]. In Section 4, the improved numerical methods which are stable but with low computational time are summarized.

2. Fluid Structure Interaction Problem Formulation

Fluid structure interaction problem can be divided into three parts: fluid problem, structure problem, and coupling condition. The model discussed here is based on the work of [9, 10]. Consider the flow of an incompressible, viscous Newtonian fluid in a two-dimensional symmetric channel with thin and deformable walls in Figure 1.

2.1. The Fluid Problem

Let and denote the horizontal and vertical coordinates, respectively. Assume that the fluid domain is supplanted by a symmetry boundary condition at the axis of symmetry. The fluid domain is denoted by , where with the lateral (top) boundary condition The fluid flow is governed by the Navier-Stokes equations: where is the fluid velocity, is the fluid pressure, is the fluid density, and is the fluid stress tensor. The fluid is assumed as Newtonian so that the fluid stress tensor is given by , where is the fluid viscosity and is the rate-of-strain tensor .

Blood is known as a suspension of red blood cells, white blood cells, and platelets in plasma. Although blood is not a Newtonian fluid, it is well accepted that, in medium-to-large arteries, the Newtonian assumption is acceptable. The non-Newtonian nature due to the particular rheology is relevant to the small arteries and capillaries where the diameter of the arteries and the size of the cell are comparable [18, 35]. For a critical review on blood flow, one can refer to [36] where the blood rheology, blood viscosity models, and conditions are listed. In this paper, only Newtonian fluid will be considered as in [911, 33].

2.2. The Structure Problem

Since the fluid structure interaction problem is complicated, the simplified model is used whenever possible. Previous studies indicated that the simplified mathematical model presenting the major physical characteristics is reasonable. A common set of simplifying models includes the use of two-dimensional models instead of the more realistic three-dimensional ones, cylindrical geometry of a section of an artery without branching, neglecting the viscoelastic term and the bending rigidity, and even a further reduction to one-dimensional models. The two-dimensional and three-dimensional models are rather complex while the one-dimensional models suffer from a serious drawback as they are not closed and oversimplifying the viscous fluid [1820].

Recent studies on the two-dimensional models with some simplification assumptions include that of Nobile and Vergara [15]. They assumed that the structure behaves as a membrane which implies that the structure is a thin elastic shell with no bending, whose thickness is neglected and which can be described by a two-dimensional manifold. A simple inertia-algebraic membrane model which considers small deformation is considered. The structure equation with initial conditions is as follows: where ,   are initial conditions. If in the particular case , it is known as the algebraic model.

Nobile [37] proposed a generalised string model derived from a cylindrical configuration. Let be the cylindrical reference surface of radius , while the longitudinal and angular displacement are neglected; thus the radial displacement is given by where account for shear deformation while introduced the viscoelastic behaviour. By neglecting the viscoelastic terms and the term of second derivatives in , the resulting equation (7) is so-called independent ring model Further simplification by neglecting the inertia term will result in the simple algebraic equation Generalised string model had been widely used as the structure of blood flow in compliant vessels and arteries [1216, 32]. It is derived from the theory of linear elasticity for a cylindrical tube with small thickness. The reference configuration is a cylindrical surface of the base circle radius that is supposed to move in radially, the longitudinal and angular displacements being neglected [16]. Causin et al. [16] suggested that the results will be more qualitative in the present example of a nonnegligible second-order term.

Guidoboni et al. [9, 10] proposed the generalized string model which includes the elastic and viscoelastic behavior. is assumed as a linearly viscoelastic thin shell, undergoing only transversal displacement :

Čanić et al. [1820] stated that there are no analytical results which are able to prove the well posedness of fluid structure interaction problems without assuming the structure model that includes the higher order derivative terms capturing the viscoelastic behavior or with the terms describing bending rigidity. They explained that the bending rigidity of the vessels walls which are being neglected might mean oversimplifying the physics. Thus, their motivation was to derive the Koiter shell equations in cylindrical coordinate. The linearly viscoelastic cylindrical Koiter shell model is given as There are two interesting features: bending rigidity plays a nonnegligible role in the approximation of the original problem, and the fluid viscous dissipation imparts long-term viscoelastic memory effects on the motion of the arterial walls.

2.3. The Coupling Condition

The coupling condition between both fluid and structure is The initial and boundary conditions for fluid velocity and the structure displacement are prescribed as

3. Numerical Approaches for Fluid Structure Interaction Problems

In this paper, the numerical approaches in solving fluid structure interaction problems are classified into two: namely, the partitioned and the monolithic approach [2024, 26]. Other numerical approaches such as the conforming and nonconforming mesh associated with the immersed boundary method in solving fluid structure interactions have been reviewed in [27].

Partitioned approach treats the fluid and structure problems as two computational fields which can be solved using two distinct solvers. The interface conditions between the fluid and structure are solved through loosely or strongly coupled algorithms [2024, 26, 27]. Figure 2 shows that loosely coupled algorithms are known as explicit algorithms while strongly coupled algorithms are known as implicit algorithms.

Partitioned approach is based on the successive solution of three subproblems and allows one to reuse the existing codes. Monolithic approach treats the fluid and structure as a single system. In other words, flow and structure problem is solved with a single code. The interfacial conditions are implicit in the solution procedure. The solution procedures of the monolithic and partitioned approach are illustrated in Figure 3 where and denote the fluid and structure solution, respectively [21, 23, 27].

3.1. Partitioned Approach

According to Sieber [34], information in loosely coupling algorithms will be exchanged between the solvers only once per time step. This implies that the fluid and structure should be in equilibrium. Then the data can be exchanged only if both fluid and structure variables are constant within each time step. Before starting the iterations process, all materials, fluid properties parameters, fluid and structure variables, time step, and the convergence criteria should be initialized. However, convergence problems might increase due to the nature of explicit coupling algorithms. Thus the choice of time-step size was restricted and it was not suitable for large structural deformations problems. Figure 4 shows the comparison in terms of stability, generality, and programming efforts for both couplings.

Andersson and Ahl [26] summarized some issues about loosely and strongly coupled algorithms. For loosely coupled algorithms, instability issue increases with decreasing the mass density ratio. Besides, the decrease in time-step size further increases the instability, known as the artificial added mass effect. Errors in the predictions along with the added mass effect caused the incorrect coupling forces that led to the instability. For strongly coupled algorithms, it was more stable for low mass density ratio. On the other hand, due to more subiterations, the computational time increases when the ratio is reduced.

Deparis et al. [38] stated that standard loosely coupling algorithms solved the fluid, geometry, and the interface explicitly and the structure implicitly. The computational cost was cheap but unstable especially when the structure was light. Several suggestions had been made to overcome the stability issues. Nobile and Vergara [29] proposed Robin interface conditions to be enforced to solve fluid and structure subproblems. Burman and Fernández [28] proposed a stabilized explicit coupling for fluid structure interaction based on Nitsche’s scheme. Numerical simulation of fluid structure interaction problems involving a viscous compressible fluid and elastic structure was considered. The explicit coupling scheme without correction had given a stable approximation with poor accuracy. High order accuracy was achieved after a few correction iterations and the results were comparable to that with implicit scheme solution [28].

It has also been suggested that the geometry and interface coupling should be treated implicitly [18]. Vierendeels et al. [39] proposed a coupling method for strongly coupled fluid structure interaction problems with partitioned solvers. They solved the reduced order models for fluid and structure problem and a small number of coupling iterations. Commercial CFD software package Fluent 6.2 was used as the fluid solver and Abaqus 6.5 was used as the structural solver. This coupling method showed a satisfactory convergence behavior. Thus, it can be summarized that the explicit partitioned algorithms are not suitable for problems in hemodynamics. It was proved to be problematic with the stability issues as the added mass effect of the fluid on the structure [16, 25, 29]. Implicit partitioned algorithms were also affected by the added mass effect as they converge slowly. Special treatments of the interface conditions had to be considered.

3.2. Monolithic Approach

Deficiencies in the partitioned method had motivated the investigation on monolithic methods [21]. Hron and Turek [40] and Hron and Mádlík [41] stated that the monolithic approach which treated the problem as a single continuum with coupling automatically takes care of the internal interface. This gets rid of the problematic interface treatment when the fluid and structure are solved separately. The results computed using monolithic approaches were ten times more accurate, but the computational cost was three to four times higher than those of the partitioned methods as stated in Michler et al. [21].

Heil [22] explained that if the fluid is incompressible or the problem is steady, the solution of a large system of coupled nonlinear algebraic equations is needed. The solution of a nonlinear system by Newton’s method was utilized since it yielded a powerful and rapidly converging scheme. However, repeated assembly of the Jacobian matrix and the solutions associated with the linear systems for Newton corrections contributed to the increase in computational cost. Thus they developed an efficient preconditioning technique that allows the rapid iterative solution instead of applying the Newton method as in [20].

Heil et al. [20] studied the fluid structure interaction in collapsible channel with monolithic and partitioned approaches. Both approaches were competitive in the test case involving steady problems. In unsteady problems, strongly coupled partitioned solvers suffered from severe convergence problems and an under-relaxation parameter needs to be applied in stabilizing the solution procedure. Monolithic solvers become more essential in unsteady problems but required an efficient precondition for the large problems, particularly in three-dimensional problems [20].

Razzaq et al. [42, 43] presented numerical simulation of fluid structure interaction in hemodynamics with monolithic approach. They restricted the research on two-dimensional models which allow the systematic tests of the proposed methods. The corresponding monolithic treatment of the fluid structure interaction problems suggested that a stable second-order time stepping scheme as well as the same finite elements for fluid and structure should be utilized. Hron and Turek [40] and Hron and Mádlík [41] applied different types of discretization in space and time. They solved the simplified two-dimensional examples with finite element and Crank-Nicolson for the space and time discretization, respectively. The resulting nonlinear algebraic system was solved by an approximate Newton’s method. The results obtained had high accuracy and robustness.

4. Improved Numerical Methods

Although stability and accuracy of partitioned approach can be improved through prediction techniques, their error remains larger than monolithic solutions [20, 21]. To date, monolithic and implicit schemes seem to be applicable in fluid structure interactions in blood flow. However, subiterations that are performed at each time step increase the computational time and computational costs [9, 10, 3033]. Several approaches based on the coupling algorithms had been proposed in recent research works, such as the semi-implicit, kinematic splitting algorithm, and the geometrical multiscale approach.

4.1. Semi-Implicit Approach

Fernández et al. [30, 31] proposed a semi-implicit scheme to solve the numerical simulation of fluid structure interaction problems involving strong added mass effect, particularly in hemodynamics. The idea of the semi-implicit scheme was to treat the added mass effect implicitly while other contributions such as geometrical nonlinearities, viscous, and convective effects are solved explicitly. Such explicit-implicit splitting can be naturally performed using a Chorin-Temam projection scheme in the fluid. The authors claimed that this scheme was numerically stable, given in theoretical and numerical evidence for a wide range of physical and discrete parameters.

However, Astorino et al. [44] stated that the scheme proposed in [30, 31] had computing limitations such that (i)it was assumed that the fluid problem is to be solved with a projection scheme; (ii)the energy was not perfectly balanced. Astorino et al. [44] then modified the scheme in [30, 31] by treating the explicit part of the coupling with Nitsche-based mortaring. The authors claimed that their scheme was independent of the added mass effect.

Badia et al. [8] proposed a similar previous semi-implicit approach which was based on the inexact block-LU factorization of the linear system. The linear system was obtained after the space-time discretization and linearization of the fluid structure interaction problems. The idea presented was to decouple the fluid velocity computation of the strongly coupled fluid structure system. Only pressure and structure unknowns were involved, with the advantage of reducing the computational costs and maintaining stability. Since the pressure was still coupled to the structure, the stability of the scheme was independent of the added mass effect.

4.2. Kinematic Splitting Algorithm

Guidoboni et al. [9, 10, 33] proposed a new version of the loosely coupled-type algorithm. The algorithm which is known as the kinematically coupled scheme is aligned with the crucial role of the kinematic condition for the proposed algorithm. This scheme applied the hypothesis that the arterial wall was modelled as a thin shell so that such scheme does not suffer from instabilities related to the high nonlinear interfacial coupling between the flow and structure. The idea of the kinematically coupled scheme was presented as follows.(1)Use operator splitting for time-discretization.(2)No iterations between the fluid and structure subproblems were required.(3)Impose the kinematic condition in strong form in order to maintain the tight link between fluid and structure in each sub-problem.(4)The fluid stress at the interface did not have to be computed explicitly.

Kinematically coupled scheme splits the structure into two parts: the hydrodynamic load exerted by the fluid on the structure and the purely elastic part without the hydrodynamic load. The hydrodynamic part, consisting of the fluid stress acting on the interface and the viscoelastic terms, is treated together with the fluid. By adding the hydrodynamic part of the structure equation to the fluid equation and by utilizing the kinematic interface condition, they deal with the inertia of both fluid and structure at the same time, thereby getting around the difficulty associated with the added mass effect. The elastic part was treated separately and this enabled the use of a wide range of structural models [9, 10].

Guidoboni et al. [9, 10] considered the incompressible, viscous Newtonian fluid in a two-dimensional channel with thin, deformable walls in the generalized string model. Time discretization via Lie’s operator splitting was applied through the scheme. Since the operator splitting was developed only for the first-order formulation, the kinematic boundary condition was applied into the structure equation to transform the second-order formulation to the first-order formulation. The overall structure of the scheme was to solve the four subproblems with different numerical schemes. Existing fluid and structure solvers can be used as “black boxes.” Numerical results of kinematically coupled scheme showed excellent agreement with those obtained using an implicit scheme [9, 10].

Bukač et al. [11] extended the work of [9, 10] by replacing the generalized string model with linearly viscoelastic cylindrical Koiter shell model. The authors tried to capture the radial and longitudinal displacement of the linearly viscoelastic Koiter shell for the underlying fluid structure interaction problem. In addition, they aimed to increase the accuracy of kinematically coupled scheme with the modified Lie’s scheme. The modified scheme was named as kinematically coupled -scheme. The results were comparable with the monolithic scheme in [8]. Such a scheme was modular and easy to implement and had low computational cost.

The idea of kinematic splitting algorithm inspired Lukáčová-Medvid’ová et al. [45] to propose a similar technique to solve the fluid structure interaction problems of non-Newtonian fluids. Lukáčová-Medvid’ová et al. [45] claimed that their approach was more general than [9, 10, 33] because they allowed the use of second-order splitting method and non-Newtonian rheology. They applied implicit backward Euler discretization to the fluid and second-order Newmark scheme for the structure. The results were conditionally stable.

4.3. Geometrical Multiscale Approach

Formaggia et al. [13, 32] mentioned that although the coupling algorithm for fluid structure interaction should be implicit, it is difficult to simulate large regions. The simulation of three-dimensional fluid structure interaction suffered a pressure wave that had been generated and reflected at the flow section. Thus, geometrical multiscale approach was proposed by coupling the detailed three-dimensional fluid structure interaction model with a one-dimensional reduced model as shown in Figure 5. They applied an implicit coupling on the three-dimensional fluid structure interaction problem and the Lax-Wendroff scheme on the one-dimensional model. The explicit numerical algorithm was proposed for the geometrical multiscale coupling. Formaggia et al. [13, 32] attempted to eliminate the spurious reflection at the flow section through geometrical multiscale approach by implying one-dimensional reduced model as the absorbing boundary condition. The results showed that the pressure wave is quite well absorbed by the one-dimensional model.

Janela et al. [46] stated that as the flow is driven by a pressure pulse generated by a constant pressure, the vessel inflates initially near the inflow boundary. The motion propagates along the vessel until it reaches the outflow section and is reflected back. Such issue can be solved through geometrical multiscale approach as proposed by [32]. Janela et al. [47] proposed several absorbing boundary conditions in order to cope with the spurious reflection. The numerical approximation of three-dimensional and one-dimensional coupling was performed through a staggered algorithm, iterating the three-dimensional fluid structure interaction and one-dimensional model. The coupling can be performed implicitly, comprising subiterations at each time step or explicitly, with no subiterations at each time step. The proposed linear absorbing boundary conditions had been proven to be effective in absorbing the pressure wave.

Both Formaggia et al. [32] and Janela et al. [46, 47] mentioned that the homogenous boundary condition will lead to energy decay property. Standard homogenous boundary condition introduces the spurious reflections of the pressure wave which will cause the structure to continue oscillating. Proper boundary conditions should be chosen in order to cope with the reflection issues caused by the three-dimensional fluid structure interaction.

5. Conclusion

Fluid structure interaction needs to be included in models of blood flow as blood interacts mechanically with the vessel wall. It is suggested that the linearly viscoelastic Koiter shell model should be adopted to model the structure of the vessel wall since it takes into account the elastic and viscoelastic behavior with bending rigidity.

The main issue in the fluid structure interaction model of blood flow model is on how to get rid of the added mass effect so that the numerical solution will be stable and the computational cost is low. The monolithic scheme has been the most commonly used approach, but it is expensive in terms of computational cost and memory requirement. To get around this problem, various ways to improve on the partitioned approach have been sought.

Classical partitioned approach considers a problem separately as fluid, structure, and interface. Problem arises when the interface is solved separately. In the kinematically coupled scheme, which is a loosely coupled partitioned-type algorithm, an operator splitting is applied instead of the problem being split into the fluid and structure subproblems. Such splitting algorithm offers the flexibility of applying any suitable numerical methods in solving each subproblem. As the computational cost is measured according to the number of iterations, the computational cost of the kinematically coupled scheme is lower, with the results obtained being as accurate as those obtained from the implicit schemes.

Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgment

Financial supports provided by Vot 03H34 and Vot 01G31, Research University Grant Scheme, Universiti Teknologi Malaysia, are gratefully acknowledged.