key: cord-0045949-sjtj1knz
authors: Sun, Pengtao; Zhang, Chen-Song; Lan, Rihui; Li, Lin
title: Monolithic Arbitrary Lagrangian–Eulerian Finite Element Method for a Multi-domain Blood Flow–Aortic Wall Interaction Problem
date: 2020-05-25
journal: Computational Science - ICCS 2020
DOI: 10.1007/978-3-030-50436-6_5
sha: 91834bc7ab68ed5b4a2a401ab0ba65585461cfb0
doc_id: 45949
cord_uid: sjtj1knz

In this paper, an arbitrary Lagrangian–Eulerian (ALE) finite element method in the monolithic approach is developed for a multi-domain blood flow–aortic wall interaction problem with multiple moving interfaces. An advanced fully discrete ALE-mixed finite element approximation is defined to solve the present fluid–structure interaction (FSI) problem in the cardiovascular environment, in which two fields of structures are involved with two fields of fluid flow, inducing three moving interfaces in between for the interactions. Numerical experiments are carried out for a realistic cardiovascular problem with the implantation of vascular stent graft to demonstrate the strength of our developed ALE-mixed finite element method.

In this paper, we consider the problem of interactions between the incompressible free viscous fluid flow and the deformable elastic structure of multiple fields. This problem is of great importance in a wide range of applications, such as in the cardiovascular area of physiology, where a significant example of this type of problem can be described as the blood fluid-vessel interaction problem. During a cardiovascular cycle, any abnormal mechanical change on the elastic property of the aortic vessel could induce some severe cardiovascular diseases (CVDs), such as the vascular stenosis, the vascular rupture, the aneurysm, the aortic dissection, and etc. Hence, a comprehensive modeling study and an accurate numerical simulation on the interactions between the blood fluid flow and the elastic aortic wall become extremely important and urgent to help medical professionals understand how the CVDs occur, grow, and fatally affect the patients, as well as to help cardiovascular physicians to find out the best medical treatment.

In the past several decades, the vascular stent graft has become one of the most popular fixation devices to carry out surgery for endovascular aortic repair (EVAR) in order to rescue CVD patients. For instance, without the need for an open surgery, the implanted vascular stent graft may cure the vascular stenosis by expanding the constrictive aortic lumen back to normal and then supporting therein as the backbone of the aorta; can treat the aortic dissection (see the left of Fig. 1 ) by substituting for the dissected aortic wall that is induced by an intima tear; and can remedy the abdominal aortic aneurysm (AAA) (see the middle of Fig. 1 ) by replacing the bulgy part of the aortic wall that is caused by a lack of elasticity in the aorta. In all cases, because of the elastic expansion of the vascular stent and the residual elasticity energy remained in the diseased aortic wall, most of the blood fluid that is trapped in the void space between the stent graft and the pathologically changed aortic wall are supposed to be squeezed back to the blood vessel lumen through the stent graft. But, there might be still a small amount of blood fluid remaining in the corner of aorta and will be eventually solidified as the thrombi. Thus, all extra blood fluids that are trapped in the shrinking void spaces are supposed to be expelled. Then the diseased blood vessel would be gradually cured to a great extent, demonstrated as either bouncing back for the aortic aneurysm case, or no longer deteriorated and then eventually normalized in an improving cardiovascular environment for the aortic dissection case (see Fig. 1 , right). However, the EVAR surgery has been long lacking a quantitative study about how and under what circumstance the CVDs are given rise to and then grow, as well as how they are changed after the implantation of a vascular stent graft, i.e., an accurate and efficient mathematical modeling and numerical simulation of the hemodynamic FSI are still largely missing in the investigations before and after the EVAR surgery in the modern medical field of CVDs. Moreover, it has been increasingly recognized that a personalized vascular stent graft is crucially needed for CVD patients, which cannot be done either without a quantitative numerical study of the hemodynamic FSI on the feasibility of the vascular stent graft in the cardiovascular environment.

In this paper, we study the effect of using an elastic material model to describe the deformable aortic wall and stent graft, and their interactions with the blood fluid flow, where, due to the loss of elasticity of the aortic wall, which is however still within the range of linear elasticity, the blood vessel lumen is inflated to form an abdominal aortic aneurysm, and the vascular stent graft is implanted into the aneurysm to separate the blood fluid into two parts. Thus four domains, i.e., two fluid domains and two structure domains, are separated by three interfaces. To model the blood fluid, we consider the Navier-Stokes equations, under the assumptions of incompressibility and Newtonian rheology, which is conventionally defined in Eulerian description. The elastic structure equation, which is conventionally defined in Lagrangian description, can be generally defined by various elastic materials such as the classical Hookean-type material model. In addition, in such a set of governing equations of FSI problem, the time and space dependence of the primary unknowns and of the moving interfaces between the fluid and the elastic structure play a significant role for the dynamic interactions in between, where, the well defined interface conditions among four domains are crucial as well.

Regarding the numerical methodology to be studied in this paper, we develop an advanced numerical method based upon the arbitrary Lagrangian-Eulerian (ALE) mapping within a monolithic approach to solve the present FSI problem. In the first place, we prefer the monolithic approach [10] , in view of its unconditional stability and the immunity of any systematic error in the implementation of interface conditions for any kind of FSI problem. Moreover, a high-performance preconditioning linear algebraic solver can also be developed and parallelized for the monolithic system without doing an alternating iteration by subdomains [1] . In contrast, the partitioned approach [2] , which decouples the FSI system and iteratively solves the fluid and the structure equations via an iteration-by-subdomain, is conditionally stable and conditionally convergent under a particular range of the physical parameters of FSI model, e.g., when fluid and structural densities are of the same order. This is called the added-mass effect [4] specifically induced by the partitioned approach. Unfortunately, the hemodynamic FSI problems are within the particular range of the added-mass effect, making the partitioned approach very difficult to use. Hence again, the monolithic approach is the primarily reliable method to be studied in this paper.

On the top of the monolithic approach, we develop an ALE method for the present FSI problem. As a type of body-fitted mesh method, ALE technique [10] has become the most accurate and also the most popular approach for solving FSI problems and other general interface problems [6, 8] , where the mesh on the interface is accommodated to be shared by both the fluid and the structure, and thus to automatically satisfy the interface conditions across the interface.

The structure of the paper is the following: in Sect. 2 we define a type of dynamic FSI model. The ALE mapping and the weak form in ALE description are defined in Sect. 3. We further develop the monolithic ALE finite element method for the present FSI model in semi-and fully discrete schemes in Sect. 4 . Numerical experiments for a realistic cardiovascular problem with the implantation of vascular stent graft are carried out in Sect. 5.

, respectively represent two fluid domains: Ω f1 t is embraced by the implanted stent graft while Ω f2 t occupies the interior of the aortic aneurysm, both of which are filled with an incompressible and viscous blood fluid, and, two elastic structure domains: Ω s1 t is the aortic wall while Ω s2 t is formed by the implanted vascular stent graft. As illustrated in Fig. 2 , four subdomains are separated by three interfaces:

also change with t ∈ I and are termed as the current (Eulerian) domains with respect to x k , in contrast to their initial (reference/Lagrangian) domains, where, Ω f 1 t : the blood vessel lumen, Ω f 2 t : the blood fluid squeezed inside the aneurysm, Ω s 1 t : the aortic wall, Ω s 2 t : the vascular graft, and three interfaces Γ 12

In this paper, we are interested in studying a pressure-driven flow through the deformable channel with a two-way coupling between the incompressible fluid and the elastic material, where the fluid motion is defined in Eulerian description and the structure motion is defined in Lagrangian description as follows

Elastic structure motion:

where, v f is the velocity of the free fluid defined in Ω f t ,û s is the displacement of elastic structure that leads to ∂ûs

, ρ f and ρ s are the constant densities of the incompressible free fluid and the elastic structure, respectively. D· Dt denotes the classical concept of material derivative as

The stress tensor of each phase is defined as

where, p f denotes the fluid pressure, μ f , μ s and λ s are constant physical parameters representing the fluid viscosity, the shear modulus, and the Lamé constant of elastic structures, respectively. Here different elastic structures may have different elastic parameters, i.e., μ s Ωs k = μ s k , λ s Ωs k = λ s k ,ρ s Ωs k =ρ s k (k = 1, 2), inducing different structure displacementsû s Ωs k =û s k (k = 1, 2).

The equation of elastic structure can be also reformulated in terms of the structure velocity, v s , as follows by substituting ∂ûs ∂t =v s andû s = u 0

whereσ s (v s ) = 2μ s ε t 0v s dτ + λ s∇ · t 0v s dτ I .

The core of ALE method is the introduction of ALE mapping, which is essentially a type of time-dependent bijective affine mapping, X t ∈ (W 1,∞ (Ω f )) d , defined from the initial (Lagrangian) domainΩ f ⊂ R d to the current (Eulerian) domain Ω f t ⊂ R d , as follows [8, 9] X t :

and

The key reason why this mapping works well for FSI problems and other standard interface problems is because the following Proposition 1 holds for any , t) , t) and its ALE time derivative,

Proposition 1 [9] . v ∈ H 1 (Ω t ) and dv

, and vice versa.

Only the current (Eulerian) domain that involves a moving boundary/interface needs the ALE mapping to produce a moving mesh therein which can accommodate the motion of the boundary/interface without breaking the mesh connectivity along the time. In practice, one way to define the affinetype ALE mapping, X t :Ω f → Ω f t , is the harmonic extension technique. For instance, to compute the ALE mapping for the fluid mesh, X t , in the present FSI problem, we solve the following Laplace equation

Once the ALE mapping X t is computed, then the fluid mesh and its moving velocity, w f , can be updated, accordingly, as:

In view of the definition of ALE time derivative dv f dt x , the interface condition (1(k)) and the time differentiation of the harmonic ALE Eq. (5), the weak form of FSI problem (1)-(3) can be defined as follows.

Introduce the following trilinear function

Applying the Green's theorem, the incompressibility (1(b)), the interface conditions (1(j)) and (5(b)), we obtain

Then (7) can be reformulated as

We first introduce finite element spaces to discretize the functional spaces defined in (6) . Conventionally, the finite element approximation to the linear elasticity equation is defined in a Lagrange-type piecewise polynomial (finite element) space,V s h ⊂V s b , to accommodate its numerical solutionv s,h . In order to satisfy the interface condition (1(j)), we need not only the meshes in Ω f t and Ω s t are conforming through Γ t , but also the finite element space of fluid (Navier-Stokes) equations, 

For any t ∈ [0, T ], we consider a discretization of the mapping X t by means of piecewise linear Lagrangian finite elements, denoted by X h,t :

So, differentiating (5) with time, we can similarly define a discrete ALE mapping forŵ f,h as follows

Choosing P 2 P 1 mixed element, P 2 element and P 1 element to construct the finite element spaces

h , respectively, we define the following semi-discrete ALE finite element approximation to FSI problem (1)-(3) based on the weak form (9). Find (v f,h ,v s,h , p f,h ,ŵ 

f,h ) ∈ (H 2 ∩L ∞ )(0, T ; V f t,h ×V s h ×Q f t,h × W f h ) such that ρ f ∂v f,h ∂t ĥ x f , ψ f,h Ω f t + β(v f,h , v f,h , ψ f,h ) Ω f t − β(ŵ f,h • X −1 h,t , v f,h , ψ f,h ) Ω f t + 1 2 v f,h ∇ · (ŵ f,h • X −1 h,t ), ψ f,h Ω f t + 2μ f D(v f,h ), D(ψ f,h ) Ω f t − (p f,h , ∇ · ψ f,h ) Ω f t + (∇ · v f,h , q f,h ) Ω f t + (∇ŵ f,h ,∇ξ f,h )Ω f +ρ s ( ∂v s,h ∂t ,ψ s,h )Ω s + 2μ s ε t 0v s,h dτ , ε(ψ s,h ) Ω s + λ s ∇ · t 0v s,h dτ ,∇ ·ψ s,h Ω s = ρ f (f f , ψ f,h ) Ω f t +ρ s (f s ,ψ s,h )Ω s − 2μ s ε(û 0 s ), ε(ψ s,h ) Ωs − λ s ∇ ·û 0 s ,∇ ·ψ s,h Ωs , ∀(ψ f,h ,ψ s,h , q f,h ,ξ f,h ) ∈ V f t,h,0 ×V s h,0 × Q f t,h ×Ŵ f h,0 .(11)

Now we develop the fully discrete monolithic ALE finite element approximation for the FSI model (1) . Introduce a uniform partition 0 = t 0 < t 1 < · · · < t N = T with the time-step size Δt = T /N. Set t n = nΔt, ϕ n = ϕ(x n , t n ) for n = 0, 1, · · · , N. Define the following backward Euler time differences based on the discrete ALE mapping X h,t :

where, X m,n = X h,tn • (X h,tm ) −1 . By the Taylor's expansion in a subtle way, we can obtain the first-order convergence with respect to Δt for (12(a)) [5] . Based on the semi-discrete scheme (11), now we can define a fully discrete ALE finite element method (FEM) for (1) as follows. diameter and continues to grow, or begins to cause symptoms, the patient may need surgery to repair the artery before the aneurysm bursts. An endovascular aortic repair (EVAR) surgery involves replacing the weakened section of the vessel with an artificial tube, called a vascular stent graft. A computational mesh for an aneurysm CVD patient implanted with a stent graft is shown in Fig. 4 , where the vascular stent graft, that is represented by an artificial blood vessel, expands inside the artery and eventually sticks to the inlet end of the aortic wall. Hence, a large blood fluid cavity is formed between the aortic aneurysm and the stent graft. Note that the blood fluid is thus separated by the stent graft into two parts: one inside the stent graft lumen driven by the incoming blood flow, while the other part is squeezed into the aneurismal cavity, waiting for being expelled due to the residual elasticity energy that remains in the diseased aorta.

To numerically model this multi-domain cardiovascular FSI problem in an accurate and effective fashion, an important technical issue needs to be clarified as follows. As shown on the right of Fig. 4 , we extend the outlet end of the aortic wall further in order to connect two blood fluid fields together through the gap between the aortic wall and the stent graft. Thus, two blood fluid fields which respectively flow through the stent graft lumen and the aneurysm cavity, are eventually contained in one single connected fluid body and share the same inlet and outlet. This is crucial for the success of our FSI simulation, since the blood fluid that is contained inside the aneurysm cavity no longer owns an inlet thus no more incoming flow after the stent graft expands and squeezes onto the Time=0.001s

Time=0.05s inlet end of the aortic wall, therefore the incompressible fluid equations that is adopted to model the blood fluid flow will be violated in this fluid part if it is not connected with the main blood flow stream through the stent graft lumen.

Another important technical issue is that, to make ALE method successfully work for this multi-domain FSI simulation all the time, we need to avoid two structure parts, the aortic wall and the stent graft, to come into contact with each other. Otherwise, there is no space left for the blood fluid to flow through the contacting places of these two structure parts, which will induce a failure of ALE method since the ALE fluid mesh in these contacting area will lose many fluid mesh cells due to the physical contact of two structure parts, recalling that the preservation of mesh connectivity all the time is a unique feature and a warrant for success of ALE method. To deal with such a possible situation during the multi-domain FSI simulation, we set up a threshold for the distance in the normal direction between the aortic wall and the stent graft, namely, we do not let them physically get in touch with each other but just treat they "contact" at a pseudo zero distance in the normal direction once the threshold is reached. In addition, they also do not make any relative shift in the tangential direction since they are merged together on the inlet end of the blood vessel, as shown in Fig. 4 . By that way, we do not have to handle the complicated contacting problem occurring between two structures at this time. Simultaneously, the ALE fluid mesh is still feasible within the threshold gap, thus our developed ALE mixed finite element method still works well for the multi-domain cardiovascular FSI problem. Such idea to avoid the contact between multiple structures is also similarly adopted in [3] where a thin layer of fluid is introduced around structures, such that there is no real "contact" between both phases. In the realistic numerical simulation, the forces are still transferred via the remaining small layer of fluid. Our approach is stable under refinement of the temporal and spatial discretization. We will carry out a more comprehensive study in the future for this approach by comparing with experiments and numerical benchmarks, in depth.

In this work, for the first time we successfully apply our monolithic ALE FEM [7, 10] and FSI codes to this multi-domain (artery, stent graft, blood fluid in the vessel lumen and the aneurysmal cavity) FSI problem with three moving interfaces in a much subtle fashion. Given a time-dependent incoming velocity profile on the inlet to simulate the pulsatile blood flow, some numerical results are gained for the velocity field of this multi-domain FSI problem in its early stage in forms of magnitude contour, vector and streamline, as shown in Fig. 5 .

We can clearly see that the larger velocity magnitudes occur in the stent graft lumen rather than in the aneurysm cavity, which is because the incoming blood fluid continues to flow through the inlet to the blood vessel lumen while the contained blood fluid inside the aneurysm cavity is less active due to the lack of residual elasticity inside the diseased aortic wall, but after a realistically long term run, the aortic wall of the aneurysm, driven by the residual elasticity, shall be able to be shrunk back to some extent, and then squeeze the remaining blood out of the aneurysm cavity through the small gap, as illustrated in Fig. 5 after 0.05 s in the computation.

Scalable parallel methods for monolithic coupling in fluidstructure interaction with application to blood flow modeling

Added-mass effect in the design of partitioned algorithms for fluid-structure problems

Eulerian techniques for fluid-structure interactions: Part II -applications

Fluid-structure interaction problems with strong added-mass effect

A monolithic arbitrary Lagrangian-Eulerian finite element analysis for a Stokes/parabolic moving interface problem

Mixed finite element analysis for an elliptic/mixedelliptic coupling interface problem with jump coefficients

Numerical simulation of an immersed rotating structure in fluid for hemodynamic applications

Convergence of a finite element/ALE method for the Stokes equations in a domain depending on time

A stability analysis for the arbitrary Lagrangian Eulerian formulation with finite elements

Modeling and simulation for fluidrotating structure interaction

such that for n = 0, 1, · · · , N − 1,

Assume that the mesh at the time step t n is denoted by T n h = T n f,h ∪T s,h , whereT s,h is the Lagrangian structural mesh inΩ s which is always fixed, T f,h is the ALE fluid mesh in Ω f t which needs to be updated all the time through the discrete ALE mapping, i.e., T n f,h =T f,h + X h,tn , where, X h,tn is only subject to the structure displacementû n s,h = Δt 2 (v n s,h +v n−1 s,h ) +û n−1 s,h onΓ for the sake of a conforming mesh across Γ t , or, X h,tn = ΔtDue to a small deformation displacement of the structure, such specific ALE mapping always guarantees a shape-regular fluid mesh in Ω f t . Suppose all the necessary solution variables at the time step t n are given as:In the following, we define an implicit iterative scheme for the FSI simulation at the current (n + 1)-th time step.

1. Let j = 1, χ n+1,0 = χ n , and T n+1,0Solve the following linearized mixed finite element equation for χ n+1,j with χ n+1,j−1 given at the j − 1 iteration:For a given tolerance ε, determine whether or not the following stopping criteria hold for the relatively iterative errors: 

A hypothetical CVD patient with a growing aneurysm is used as an illustration example. The patient's CT imaging data (See Fig. 3 ) are collected to construct a geometrical model (See the right of Fig. 2) . The cardiovascular parameters of this patient are estimated based on the patient's routine test results such as blood work, blood pressure, diabetes. Usually, if an aneurysm becomes about 2 in. in