Keywords

1 Introduction

Fluid-structure interaction (FSI) problems remain among the most challenging problems in computational mechanics and computational fluid dynamics. The difficulties in the simulation of these problems stem from the fact that they rely on the coupling of models of different nature: Lagrangian in the solid and Eulerian in the fluid. The so-called Arbitrary Lagrangian Eulerian (ALE) method [7, 9] copes with this difficulty by adapting the fluid solver along the deformations of the solid medium in the direction normal to the interface. These methods allow us to accurately account for continuity of stresses and velocities at the interface, where the body-fitted conforming mesh is used and the surface mesh is accommodated to be shared between the fluid and solid, and thus to automatically satisfy the interface kinematic condition. In this paper, we use the ALE method to reformulate and solve fluid equations on a moving fluid mesh which moves along with structure motions through the interface.

In contrast to a large amount of FSI literatures which are dedicated to either a rigid structure [8, 11], or a non-rotational structure [4, 7], or a stationary fluid domain (thus the interface) [5, 13], few works are contributed to FSI problems with a rotational and deformable structure. This is mainly due to the fact that the corresponding mathematical and mechanical model is lacking. Some relevant models arise from the field of graphics and animation applications for the co-rotational linear elasticity [14], but still far away from our need. To simulate such a FSI problem with a more realistic operating condition, we first derive a linearized constitutive model of structure combining the elastic deformation with the rotation, then present its weak formulation and develop a specific ALE mapping technique to deal with the co-rotational fluid and structure. Furthermore, we describe an efficient monolithic iterative algorithm and ALE-type mixed finite element method to solve the hydrodynamic/hemodynamic FSI problem.

The application of our developed monolithic ALE finite element method in this paper is to study the numerical performance of an artificial heart pump running in the hemodynamic FSI system. It has been always significant to study efficient and accurate numerical methods for simulating hemodynamic FSI problems in order to help patients being recovered from a cardiovascular illness, especially from a heart failure. The statistical data tell us about 720,000 people in the U.S. suffer heart attacks each year, in which about 2,000–2,300 heart transplants are performed annually in the United States [1]. The vast majority of these patients are unsuitable to take the heart transplantation, or, are waiting for the proper human heart to transplant. Under such circumstance, the artificial heart pump is thus the only choice to sustain their lives. Since 1980s, the left ventricle auxiliary device (LVAD), also usually called the artificial heart pump, has become an effective treatment by breaking the vicious cycle of heart failure. Currently, although a massive amount of heart failure patients need the artificial heart to assistant their treatment, an individually designed artificial heart for a specially treated patient is still lacking since such individual design heavily relies on an accurate modeling and solution of the artificial heart-bloodstream-vascular interactions, which may tell the doctors the accurate shape and size of the artificial heart pump to be manufactured, and the exact location to implant it into the cardiovascular system in order to help on saving the patients’ lives from their heart failure illnesses [2].

The goal of this paper is to develop advanced modeling and novel numerical techniques for hemodynamic FSI problems in order to effectively perform stable, precise and state of the art simulation for the cardiovascular interactions between the artificial heart pump and blood flow. In our case, the artificial heart pump is equipped with a rotating impeller that may bear a small deformation under the bloodstream impact but large rotations, simultaneously, interacting with the incompressible blood flow through the interface transmissions.

2 Model Description of FSI Problems

Let us consider a deformable elastic structure in \(\varOmega _{\mathrm s}\) that is immersed in an incompressible fluid domain \(\varOmega _{\mathrm f}\). \(\varOmega _{\mathrm f}\cup \varOmega _{\mathrm s}=\varOmega \in R^d,\ \varOmega _{\mathrm f}\cap \varOmega _{\mathrm s}=\emptyset \). The interface of fluid and structure is \(\varGamma =\partial \varOmega _{\mathrm s}\). As shown in Fig. 1, the Eulerian coordinates in the fluid domain are described with the time invariant position vector \({\varvec{x}}\); the solid positions in the initial configuration \(\hat{\varOmega }_{\mathrm s}\) and the current configuration \(\varOmega _{\mathrm s}\) are represented by \({\varvec{X}}\) and \({\varvec{x}}({\varvec{X}},t)\), respectively. We use \({\varvec{v}}_{\mathrm f}\) to denote the velocity of fluid, \(p_{\mathrm f}\) the pressure of fluid, \({\varvec{v}}_{\mathrm s}\) the velocity of solid structure, and \({\varvec{u}}_{\mathrm s}\) the displacement of structure. Thus, \({\varvec{u}}_{\mathrm s}({\varvec{x}}({\varvec{X}},t),t)=\hat{\varvec{u}}_{\mathrm s}({\varvec{X}},t) = {\varvec{x}}-{\varvec{X}}\), and \({\varvec{v}}_{\mathrm s}=\frac{d{\varvec{x}}}{dt}=\frac{\partial {\varvec{u}}_{\mathrm s}}{\partial t}\). The Jacobian matrix \({\varvec{F}}=\frac{\partial {\varvec{x}}}{\partial {\varvec{X}}}\) describes the structure deformation gradient.

Fig. 1.
figure 1

Schematic domain of FSI problem.

Full size image

Fluid Motion in Eulerian Description

$$\begin{aligned} \rho _{\mathrm f}\left( \frac{\partial \varvec{v}_{\mathrm f}}{\partial t}+\varvec{v}_{\mathrm f}\cdot \nabla \varvec{v}_{\mathrm f}\right)= & {} \nabla \cdot \varvec{\sigma }_{\mathrm f}+\rho _{\mathrm f}\varvec{f}_{\mathrm f}, \quad \text{ in } \varOmega _{\mathrm f}\end{aligned}$$
(1)
$$\begin{aligned} \nabla \cdot \varvec{v}_{\mathrm f}= & {} 0, \quad \text{ in } \varOmega _{\mathrm f}\end{aligned}$$
(2)
$$\begin{aligned} \varvec{\sigma }_{\mathrm f}= & {} {-p_{\mathrm f}\varvec{I}+\mu _{\mathrm f}\left( \nabla \varvec{v}_{\mathrm f}+\left( \nabla \varvec{v}_{\mathrm f}\right) ^T\right) },\end{aligned}$$
(3)
$$\begin{aligned} \varvec{v}_{\mathrm f}= & {} \varvec{v}_{\mathrm B}, \quad \text{ on } \partial \varOmega \end{aligned}$$
(4)
$$\begin{aligned} \varvec{v}_{\mathrm f}(\varvec{x},0)= & {} \varvec{v}_{{\mathrm f0}}, \quad \text{ in } \varOmega _{\mathrm f} \end{aligned}$$
(5)

Solid Motion in Lagrangian Description

$$\begin{aligned} \rho _{\mathrm s}\frac{\partial ^2\hat{\varvec{u}}_{\mathrm s}}{\partial t^2}= & {} \nabla \cdot (\varvec{J}\varvec{\sigma }_{\mathrm s}\varvec{F}^{-T})+\rho _{\mathrm s}\hat{\varvec{f}}_{\mathrm s},\quad \text{ in } \hat{\varOmega }_{\mathrm s}\end{aligned}$$
(6)
$$\begin{aligned} \hat{\varvec{u}}_{\mathrm s}(\varvec{X},0)= & {} \hat{\varvec{u}}_{\mathrm s0},\ \text{ and } \ \frac{\partial \hat{\varvec{u}}_{\mathrm s}}{\partial t}(\varvec{X},0)\ =\ \hat{\varvec{v}}_{\mathrm s0},\quad \text{ in } \hat{\varOmega }_{\mathrm s} \end{aligned}$$
(7)

where the stress tensor \(\varvec{\sigma }_{\mathrm s}\) in the reference configuration \(\varvec{X}\), denoted by \(\hat{\varvec{\sigma }}_{\mathrm s}\), is defined by the first Piola-Kirchhoff stress tensor, \(\varvec{P}\), as \(\hat{\varvec{\sigma }}_s=\varvec{P}=\varvec{J}\varvec{\sigma }_{\mathrm s}\varvec{F}^{-T}=\varvec{F}\varvec{S}\), where \(\varvec{S}=\varvec{J}\varvec{F}^{-1}\varvec{\sigma }_{\mathrm s}\varvec{F}^{-T}\) denotes the second Piola-Kirchhoff stress tensor. For a compressible St. Venant-Kirchhoff (STVK) material [10], \( \varvec{\sigma }_{\mathrm s}(x)=\frac{1}{\varvec{J}}\varvec{F}(2\mu _{\mathrm s} \varvec{E}+\lambda _{\mathrm s} (tr\varvec{E})\varvec{I})\varvec{F}^T\), thus \(\varvec{S}=2\mu _s\varvec{E}+\lambda _{\mathrm s} (tr\varvec{E})\varvec{I}\), where \(\varvec{E}\) is the Green-Lagrangian finite strain tensor, defined as \( \varvec{E}=(\varvec{F}^T\varvec{F}-\varvec{I})/2=[\nabla _{\varvec{X}} \varvec{u}_{\mathrm s}+(\nabla _{\varvec{X}} \varvec{u}_{\mathrm s})^T+(\nabla _{\varvec{X}} \varvec{u}_{\mathrm s})^T\nabla _{\varvec{X}} \varvec{u}_{\mathrm s}]/2 \). \(\lambda _{\mathrm s},\,\mu _{\mathrm s}\) are the Lam\(\acute{\text{ e }}\)’s constants.

Interface Conditions

$$\begin{aligned} \varvec{v}_{\mathrm f}= & {} \varvec{v}_{\mathrm s}, \quad \ \ \ \ \ \text {on } \varGamma \end{aligned}$$
(8)
$$\begin{aligned} \varvec{\sigma }_{\mathrm f}\cdot \varvec{n}= & {} \varvec{\sigma }_{\mathrm s}\cdot \varvec{n}, \quad \text {on } \varGamma \end{aligned}$$
(9)

which are called the kinematic and dynamic lateral interface conditions, respectively, describing the continuities of velocity and of normal stress.

ALE Mapping Technique. Define \(\mathcal {A}(\cdot ,t)=\mathcal {A}_t:\hat{\varOmega }(\varvec{X})\mapsto \varOmega (\varvec{x})\) is an arbitrary diffeomorphism, satisfying the following conventional ALE (harmonic) mapping:

$$\begin{aligned} \left\{ \begin{array}{rll} -\varDelta \mathcal {A}&{}=0, &{}\text{ in } \hat{\varOmega }_{\mathrm {f}},\\ \mathcal {A}&{}= 0,&{} \text{ on } \partial \hat{\varOmega }_{\mathrm {f}}\backslash \hat{\varGamma },\\ \mathcal {A}&{}= \hat{\varvec{u}}_{\mathrm s},&{} \text{ on } \hat{\varGamma }, \end{array} \right. \end{aligned}$$
(10)

where \(\mathcal A\) is the displacement of fluid grid which equals to the structure displacement on the interface, resulting in a moving fluid grid according to the structure motion and further fulfilling the kinematic interface condition.

Fluid Motion in ALE Description. The momentum equation of fluid motion (1) in ALE description is given as

$$\begin{aligned} \rho _{\mathrm f}\partial ^{\mathcal A}_t \varvec{v}_{\mathrm f}+\rho _{\mathrm f}((\varvec{v}_{\mathrm f}-\varvec{w})\cdot \nabla )\varvec{v}_{\mathrm f}=\nabla \cdot \varvec{\sigma }_{\mathrm f}+\rho _{\mathrm f}\varvec{f}_{\mathrm f}, \end{aligned}$$
(11)

where \(\varvec{w}=\frac{\partial \mathcal A}{\partial t}\circ \mathcal A_t^{-1}\) denotes the velocity of fluid grid and

$$\partial ^{\mathcal A}_t \varvec{v}_{\mathrm f}|_{(x,t)}:=\left. \frac{\partial \left[ \varvec{v}_{\mathrm f}\left( \mathcal A(\varvec{X},t),t\right) \right] }{\partial t}\right| _{(\varvec{X},t)=\left( \mathcal {A}_t^{-1}(x),t\right) }=\frac{\partial \varvec{v}_{\mathrm f}}{\partial t}+(\varvec{w}\cdot \nabla ) \varvec{v}_{\mathrm f}$$

is the ALE time derivative. For the simplicity of notation, in what follows, we suppose \(\varvec{f}_{\mathrm f}=\varvec{f}_{\mathrm s}=0\) by assuming no external force is acted on fluid and structure.

3 A Monolithic Weak Formulation of FSI

Introduce \(\hat{\varvec{V}}_{\mathrm s}=\{\hat{\varvec{v}}_{\mathrm s}\in (H^1(\hat{\varOmega }_{\mathrm s}))^d|\hat{\varvec{v}}_{\mathrm s}(\varvec{X})=\varvec{v}_{\mathrm f}(\varvec{x}(\varvec{X},t))\circ \mathcal A\ \text{ on }\ \hat{\varGamma }\}\), \(\varvec{V}_{\mathrm f}=\{\varvec{v}_{\mathrm f}\in (H^1(\varOmega _{\mathrm f}))^d|\varvec{v}_{\mathrm f}=\varvec{v}_{\mathrm B}\ \text{ on }\ \partial \varOmega \}\), \(\varvec{W}_{\mathrm f}=L^2(\varOmega _{\mathrm f})\), \(\hat{\varvec{Q}}_{\mathrm f}=\{\mathcal A\in (H^1(\hat{\varOmega }_{\mathrm f}))^d|\mathcal A=0\ \text{ on }\ \partial \hat{\varOmega }_{\mathrm f}\cap \partial \hat{\varOmega }, \text{ and } \mathcal A=\hat{\varvec{u}}_{\mathrm s} \ \text{ on }\ \hat{\varGamma }\}\), we define a monolithic ALE weak formulation of FSI model as follows: find \((\hat{\varvec{v}}_{\mathrm s},\varvec{v}_{\mathrm f},p,\mathcal A)\in (\hat{\varvec{V}}_{\mathrm s} \oplus \varvec{V}_{\mathrm f}\oplus \varvec{W}_{\mathrm f}\oplus \hat{\varvec{Q}}_{\mathrm f})\times L^2([0,T])\) such that

$$\begin{aligned} \left\{ \begin{aligned} \left( \rho _{\mathrm s}\frac{\partial \hat{\varvec{v}}_{\mathrm s}}{\partial t},\phi \right) +\left( \hat{\varvec{\sigma }}_{\mathrm s}\left( \hat{\varvec{u}}_{\mathrm s}^0+\int _0^t\hat{\varvec{v}}\left( \tau \right) d\tau \right) ,\varepsilon \left( \phi \right) \right)&\\ +\left( \rho _{\mathrm f}\partial ^{\mathcal A}_t{\varvec{v}}_{\mathrm f},\psi \right) +\left( \rho _{\mathrm f}\left( {\varvec{v}}_{\mathrm f}-\varvec{w}\right) \cdot \nabla {\varvec{v}}_{\mathrm f},\psi \right) +\left( {\varvec{\sigma }}_{\mathrm f},\varepsilon \left( \psi \right) \right)&=0,\\ \left( \nabla \cdot {\varvec{v}}_{\mathrm f},q\right)&=0,\\ \left( \nabla \mathcal {A},\nabla \xi \right)&=0, \\ \forall \phi \in \hat{\varvec{V}}_{\mathrm s},\ \psi \in \varvec{V}_{\mathrm f},\ q\in \varvec{W}_{\mathrm f},\ \xi \in \hat{\varvec{Q}}_{\mathrm f},&\\ \end{aligned} \right. \end{aligned}$$
(12)

where \(\varepsilon (\varvec{u})=\frac{1}{2}\left( \nabla \varvec{u}+\left( \nabla \varvec{u}\right) ^T\right) \). The boundary integrals that arise from the integration by parts for fluid and structure equations are cancelled due to the continuity of normal stress condition (9). To efficiently implement the kinematic condition (8), we employ the master-slave relation technique, namely, two sets of grid nodes are defined on the interface, one belongs to the fluid grid and the other one belongs to the structure grid, both share the same position and the same degrees of freedom of velocity on the interface. The usage of structure velocity \(\hat{\varvec{v}}_{\mathrm s}\) instead of the displacement \(\hat{\varvec{u}}_{\mathrm s}\) as the principle unknown in (12) has an advantage to precisely apply the master-slave relation, where, \(\hat{\varvec{v}}_{\mathrm s}=\frac{\partial \hat{\varvec{u}}_{\mathrm s}}{\partial t}\) or \(\hat{\varvec{u}}_{\mathrm s}=\hat{\varvec{u}}_{\mathrm s}^0+\int _0^t\hat{\varvec{v}}_{\mathrm s}(\tau )d\tau \), then (8) becomes \(\hat{\varvec{v}}_{\mathrm s}(\varvec{X})=\varvec{v}_{\mathrm f}(\varvec{x}(\varvec{X},t))\circ \mathcal A\ \text{ on }\ \hat{\varGamma }\).

4 FSI Problem Involving a Rotating Elastic Structure

Based on the constitutive law of STVK material, we rewrite the structure equation in Lagrangian description (6) as follows

$$\begin{aligned} \rho _{\mathrm s} \frac{\partial ^2 \hat{\varvec{u}}_{\mathrm s}}{\partial t^2}=\nabla \cdot \left( \lambda _{\mathrm s}(tr{\varvec{E}}){\varvec{F}}+2\mu _{\mathrm s}{\varvec{F}}{\varvec{E}}\right) . \end{aligned}$$
(13)

Suppose \({\varvec{X}}_0\) is the structure centroid in statics, but the center of mass in dynamics, we define the structure displacement with regard to \(\varvec{X}_0\) as [6]

$$\begin{aligned} {\varvec{x}}-{\varvec{X}}_0=\varvec{R}({\varvec{X}}+\hat{\varvec{u}}_{\mathrm d}-{\varvec{X}}_0), \end{aligned}$$
(14)

where \(\hat{\varvec{u}}_{\mathrm d}\) is the deformation displacement in the local coordinate system whose origin is \(\varvec{X}_0\), \(\varvec{R}\) is the rotator from the base configuration \(\mathcal {C}^0\) to the corotated configuration \(\mathcal {C}^R\), given by \(\varvec{R}=\varvec{T}_R\varvec{T}_0^T\), here \(\varvec{T}_0\) and \(\varvec{T}_R\) are the global-to-local displacement transformations [6]. \(\varvec{T}_R,\ \varvec{T}_0\) and \(\varvec{R}\) are all orthogonal matrices. We specialize one case of which the rotation of axis is Z-axis, resulting in

$$\begin{aligned} \varvec{R}(\theta )=\left( \begin{array}{r@{\quad }r@{\quad }r} \cos (\theta -\theta _0)&{}-\sin (\theta -\theta _0)&{}0\\ \sin (\theta -\theta _0)&{}\cos (\theta -\theta _0)&{}0\\ 0&{}0&{}1\\ \end{array} \right) , \end{aligned}$$
(15)

where \(\theta _0\) is a fixed initial angle from the globally-aligned configuration \(\mathcal {C}^G\) in \(\mathcal {C}^0\) and \(\theta =\theta (t)\) is a time-dependent angle from \(\mathcal {C}^G\) in \(\mathcal {C}^R\). Without loss of generality, we let \(\theta _0=0\). Hence, the total structure displacement, \(\hat{\varvec{u}}_s\), can be defined as

$$\begin{aligned} \hat{\varvec{u}}_{\mathrm s}={\varvec{x}}-{\varvec{X}}=(\varvec{R}-{\varvec{I}})({\varvec{X}}-{\varvec{X}}_0)+\varvec{R}\hat{\varvec{u}}_{\mathrm d}=\hat{\varvec{u}}_\theta +\varvec{R}\hat{\varvec{u}}_{\mathrm d}, \end{aligned}$$
(16)

where \(\hat{\varvec{u}}_\theta =(\varvec{R}-{\varvec{I}})({\varvec{X}}-{\varvec{X}}_0)\) is the rotational displacement of the structure.

From (14), we obtain \( {\varvec{F}}= \varvec{R}({\varvec{I}}+\varPhi ), \text{ where } \varPhi =\nabla \hat{\varvec{u}}_{\mathrm d} \). Then, \( {\varvec{E}}=({\varvec{F}}^T{\varvec{F}}-{\varvec{I}})/2=(\varPhi +\varPhi ^T+\varPhi ^T\varPhi )/2. \) Thus, the first Piola-Kirchhoff stress \(\varvec{P}={\varvec{F}}(\lambda _{\mathrm s}(tr{\varvec{E}})+2\mu _{\mathrm s}{\varvec{E}})\), leading to a function of \(\varPhi \) as follows

$$ P(\varPhi )=\varvec{R}({\varvec{I}}+\varPhi )\left[ \frac{1}{2}\lambda _{\mathrm s}\left( tr\varPhi +tr\varPhi ^T+tr\left( \varPhi ^T\varPhi \right) \right) +\mu _{\mathrm s}\left( \varPhi +\varPhi ^T+\varPhi ^T\varPhi \right) \right] . $$

The following linear approximation is then derived by Taylor expansion at \(\varPhi =0\),

$$ P\approx \frac{1}{2}\lambda _{\mathrm s}\varvec{R}\left( tr\varPhi +tr\varPhi ^T\right) +\mu _{\mathrm s}\varvec{R}\left( \varPhi +\varPhi ^T\right) =\varvec{R}\left( \lambda _{\mathrm s} tr\left( \varepsilon \left( \hat{\varvec{u}}_{\mathrm d}\right) \right) {\varvec{I}}+2\mu _{\mathrm s}\varepsilon \left( \hat{\varvec{u}}_{\mathrm d}\right) \right) , $$

where \(\varepsilon (\hat{\varvec{u}}_{\mathrm d})=\frac{1}{2}(\nabla \hat{\varvec{u}}_{\mathrm d}+(\nabla \hat{\varvec{u}}_{\mathrm d})^T)\) is the linear approximation of \({\varvec{E}}(\hat{\varvec{u}}_{\mathrm d})\). We may equivalently rewrite \(\hat{\varvec{\sigma }}_{\mathrm s}=\varvec{R}{\varvec{D}}\varepsilon (\hat{\varvec{u}}_{\mathrm d})\), where \({\varvec{D}}\) denotes the constitutive matrix in terms of the Young’s modulus E and Poisson’s ratio \(\nu \) [3]. Thus, (13) is approximated by the following co-rotational linear elasticity model

$$\begin{aligned} \rho _{\mathrm s} \frac{\partial ^2 \hat{\varvec{u}}_{\mathrm s}}{\partial t^2}= & {} \nabla \cdot \left( \varvec{R}{\varvec{D}}\varepsilon \left( \hat{\varvec{u}}_{\mathrm d}\right) \right) . \end{aligned}$$
(17)

The weak formulation of (17) is defined as: find \(\hat{\varvec{u}}_{\mathrm s}\in L^2([0,T];\hat{\varvec{V}}_{\mathrm s})\) such that

$$\begin{aligned} \begin{aligned} \left( \rho _{\mathrm s} \frac{\partial ^2 \hat{\varvec{u}}_{\mathrm s}}{\partial t^2},\phi \right) + \left( {\varvec{D}}\varepsilon \left( \hat{\varvec{u}}_{\mathrm d}\right) ,\varepsilon \left( \varvec{R}^T\phi \right) \right) =&\int _{\partial \varOmega _{\mathrm s}}{\varvec{\sigma }}_{\mathrm s}{\varvec{n}}_{\mathrm s}\phi \ ds,\quad \forall \phi \in \hat{\varvec{V}}_{\mathrm s}. \end{aligned} \end{aligned}$$
(18)

Due to (16), we have \(\varvec{R}\hat{\varvec{u}}_{\mathrm d}=\hat{\varvec{u}}_{\mathrm s}-\hat{\varvec{u}}_\theta \), thus \(\varepsilon (\hat{\varvec{u}}_{\mathrm d})=\varepsilon (\varvec{R}^T\hat{\varvec{u}}_{\mathrm s})-\varepsilon (\varvec{R}^T\hat{\varvec{u}}_\theta )\). Note that \(\hat{\varvec{u}}_\theta =(\varvec{R}-{\varvec{I}})({\varvec{X}}-{\varvec{X}}_0)\), then \({\varvec{D}}\varepsilon (\varvec{R}^T\hat{\varvec{u}}_\theta )={\varvec{D}}\varepsilon ((I-\varvec{R}^T)({\varvec{X}}-{\varvec{X}}_0))\). Therefore, (18) is rewritten as

$$\begin{aligned} \left( \rho _{\mathrm s} \frac{\partial ^2 \hat{\varvec{u}}_{\mathrm s}}{\partial t^2},\phi \right) + \left( {\varvec{D}}\varepsilon \left( \varvec{R}^T\hat{\varvec{u}}_{\mathrm s}\right) ,\varepsilon \left( \varvec{R}^T\phi \right) \right) =&\left( {\varvec{D}}\varepsilon ((I-\varvec{R}^T)({\varvec{X}}-{\varvec{X}}_0)),\varepsilon (\varvec{R}^T\phi )\right) \nonumber \\&+\int _{\partial \varOmega _{\mathrm s}}{\varvec{\sigma }}_{\mathrm s}{\varvec{n}}_{\mathrm s}\phi \ ds,\ \forall \phi \in \hat{\varvec{V}}_{\mathrm s}, \end{aligned}$$
(19)

where, the stiffness term on the L.H.S. is symmetric positive definite, the first term on the R.H.S. is the body force contribution from the rotation, the boundary integral term shall be canceled later due to the continuity condition of normal stress (9). The remaining unknown quantity in (19) is the rotational matrix \(\varvec{R}\) if the structure rotation is passive. In the following we introduce two propositions to the computation of \(\varvec{R}\) by means of the structure velocity \(\hat{\varvec{v}}_{\mathrm s}\), both proofs are relatively easy and thus are omitted here.

Proposition 1

The angular velocity \(\omega (t)=(\omega _1,\omega _2,\omega _3)^T\) satisfies

$$\begin{aligned} I\omega =\int _{\varOmega _{\mathrm s}}\rho _{\mathrm s}{\varvec{r}}\times \hat{\varvec{v}}_{\mathrm s}\ d{\varvec{X}}, \end{aligned}$$
(20)

where, \(I=\int _{\varOmega _{\mathrm s}}\rho _{\mathrm s}{\varvec{r}}^2\ d{\varvec{X}}\) is the inertia, \(\varvec{r}\) is the position vector to the rotation of axis which has the following relation with \(\omega \)

$$\begin{aligned} \frac{d{\varvec{r}}}{dt}=\omega \times {\varvec{r}}=\tilde{\omega }{\varvec{r}}=\left( \begin{array}{ccc} 0&{}-\omega _3&{}\omega _2\\ \omega _3&{}0&{}-\omega _1\\ -\omega _2&{}\omega _1&{}0\\ \end{array} \right) \left( \begin{array}{c} r_1\\ r_2\\ r_3\\ \end{array} \right) ,\ \text{ and } \ \hat{\varvec{v}}_{\mathrm s}=\frac{\partial \hat{\varvec{u}}_{\mathrm s}}{\partial t}=\frac{d{\varvec{r}}}{dt}. \end{aligned}$$
(21)

Proposition 2

\(\varvec{R}\) satisfies the following O.D.E.

$$\begin{aligned} \frac{d\varvec{R}}{dt}=\tilde{\omega }\varvec{R}, \quad \text{ and }\quad \varvec{R}(0)=\varvec{R}_0. \end{aligned}$$
(22)

Now, by iteratively solving (19), (20), and (22) together, we are able to sequentially obtain the structure velocity \(\hat{\varvec{v}}_{\mathrm s}\), the angular velocity \(\omega \), then the rotational matrix \(\varvec{R}\), further, the structure displacement \(\hat{\varvec{u}}_{\mathrm s}\).

5 Algorithm Description

We first split \(\varOmega _{\mathrm f}=\varOmega _{\mathrm {sf}}\cup \varOmega _{\mathrm {rf}}\), where \(\varOmega _{\mathrm {rf}}\) is an artificial buffer zone containing the rotating elastic structure inside, which could be a disk in 2D or a cylinder in 3D, as shown in Fig. 2. And, the size of \(\varOmega _{\mathrm {rf}}\) which is characterized by the radius of its cross section is usually taken as the middle between the outer boundary of the flow channel and the rotating structure in order to guarantee the mesh quality inside and outside of the artificial buffer zone \(\varOmega _{\mathrm {rf}}\). Both \(\varOmega _{\mathrm {rf}}\) and \(\varOmega _{\mathrm s}\) suppose to rotate about the same rotation of axis with the same angular velocity under the circumstance of a small strain arising from the structure. Suppose all the necessary solution data from the last time step are known: \({\varvec{v}}_{\mathrm f}^{n-1}\), \(\hat{\varvec{v}}_{\mathrm s}^{n-1},\ \hat{\varvec{u}}_{\mathrm s}^{n-1}\), \(\varvec{R}^{n-1}\), \(w^{n-1}\), \(\mathcal {A}^{n-1}\), and the mesh on the last time step, \(\mathbb {T}_h^{n-1}=\mathbb {T}_{{\mathrm s \mathrm f},h}\cup \mathbb {T}_{{\mathrm {rf}},h}^{n-1}\cup \hat{\mathbb {T}}_{{\mathrm s},h}\), where \(\hat{\mathbb {T}}_{{\mathrm s},h}\) is the Lagrangian structure mesh in \(\hat{\varOmega }_{\mathrm s}\) which is always fixed, \(\mathbb {T}_{{\mathrm {sf}},h}\) is the mesh in the stationary fluid domain \(\varOmega _{\mathrm {sf}}\) which is also fixed, and \(\mathbb {T}_{{\mathrm {rf}},h}\) is the mesh in the rotational fluid domain \(\varOmega _{\mathrm {rf}}\) which needs to be computed all the time. Thus, we decompose the displacement of the rotational fluid mesh \(\mathbb {T}_{{\mathrm {rf}},h}\) to two parts: the rotational part \(\varvec{u}_\theta \) and the deformation part \(\mathcal A\) which is attained from a specific ALE mapping defined in (26), where \(\mathcal A\) is only subject to a structure deformation displacement \(\hat{\varvec{u}}_{\mathrm s}-\varvec{u}_\theta \) on \(\hat{\varGamma }\), and the local adjustment of the mesh nodes on \(\partial \varOmega _{\mathrm {rs}}=\partial \varOmega _{\mathrm {rf}}\cap \partial \varOmega _{\mathrm {sf}}\) for the sake of a conforming fluid mesh across \(\partial \varOmega _{\mathrm {rs}}\). Such specific ALE mapping always guarantees a shape-regular fluid mesh in \(\varOmega _{\mathrm {rf}}\), and still conforms with the stationary fluid mesh in \(\varOmega _{\mathrm {sf}}\) through \(\partial \varOmega _{\mathrm {rs}}\).

Fig. 2.
figure 2

An illustration of the rotational part (a buffer zone) of the fluid (green) separates the structure (blue) from the stationary part of the fluid (red) in 2D & 3D. (Color figure online)

Full size image

In the following, we define an implicit relaxed fixed-point iterative scheme for the FSI simulation at the current n-th time step, i.e., we first iteratively solve the implicit nonlinear momentum equations of both fluid and structure on a known fluid mesh \({\mathbb {T}}_{{\mathrm f},h}\) obtained from the previous step and the fixed structure mesh \(\hat{\mathbb {T}}_{{\mathrm s},h}\) until the convergence, then compute a new fluid mesh \({\mathbb {T}}_{{\mathrm f},h}\) based on the newly obtained structure velocity, and start another inner iteration to solve the momentum equations on the new fluid mesh. Continue this fixed-point iteration until the fluid mesh is converged, then march to the next time step. The detailed algorithms are described in Algorithms 1–3.

figure a
figure b

In Algorithm 2, the ALE time derivative, \(\partial ^{\mathcal A}_t \varvec{v}_{\mathrm f}\), is discretized as

$$\begin{aligned} \partial ^{\mathcal A}_t \varvec{v}_f\approx \frac{\varvec{v}_f({\varvec{x}},t_n)-\varvec{v}_f(\mathcal A_{n-1}\circ \mathcal {A}_{n}^{-1}({\varvec{x}}),t_{n-1})}{\varDelta t}, \end{aligned}$$
(24)

where \(\mathcal A_{n-1}\circ \mathcal A_{n}^{-1}({\varvec{x}})\) is on the corresponding grid node at \(t=t_{n-1}\) as long as \({\varvec{x}}\) is on one grid node at \(t=t_n\). Thus the interpolation between different time levels is avoided due to the mesh connectivity that is guaranteed by the ALE mapping.

figure c

6 Application to the Artificial Heart Pump

To test our model and numerical method, in this section we study the hemodynamic interaction of blood flow with an artificial heart pump which is embedded into the blood vessel and immersed in the blood flow. The artificial heart pump to be studied in this paper consists of three parts: the rotating rotor in the middle and the unmoving head guide and tail guide on two terminals. It locates close to the inlet of the vascular lumen and its rotation of axis is fixed, as illustrated in Fig. 3. In principle, the pump rotor plays a role of impeller to increase the blood pressure when the blood flows through it, the head/tail guides are used to stabilize the incoming and outgoing blood flow, altogether helping the failing human heart to propel a stable blood flow through the entire human body.

Fig. 3.
figure 3

Computational meshes – Left: the artificial heart pump (3 parts from the left end to the right end: head guide, rotor and tail guide); Right: the blood flow mesh in a vascular lumen, where the rotational part that immerses the pump rotor is separated from the stationary part by two discs between the pump rotor and the tail/head guide of the pump. The meshes on two discs are shown in Fig. 4.

Full size image

To apply our ALE finite element method to the above specific artificial heart pump, we need to separate the entire computational domain shown on the right of Fig. 3 to three parts: the rotational part containing the pump rotor and the surrounding fluid area, and two stationary parts including the unmoving head/tail guide of the pump and surrounding fluid regions. Two discs with particular meshes shown in Fig. 4 are made between the head/tail guide and the rotor to fulfill this purpose. The mesh on each disc is made along a series of concentric circles, by which the local adjustment of the mesh motion on the interface of the rotational fluid region (\(\varOmega _{\mathrm {rf}}\)) and the stationary fluid region (\(\varOmega _{\mathrm {sf}}\)), \(\partial \varOmega _{\mathrm {rs}}\), can then be easily calculated.

Fig. 4.
figure 4

Interface meshes on \(\partial \varOmega _{\mathrm {rs}}\) between the stationary fluid and the rotational fluid regions.

Full size image

To start the artificial heart pump, the rotor is always given an initial angular velocity, \(\omega \), in the unit of revolution per minute (rpm), starting from \(\omega =1000\) rpm then running up to \(\omega =8000\) rpm to work as an impeller of the blood flow by increasing the blood pressure to normal level. In turn, the artificial heart pump itself bears a relatively tiny deformation due to its nearly rigid structure material (the Young’s modulus is up to \(1.95\times 10^{11}\) Pa), thus the developed co-rotational linear elasticity model (17) works well for the artificial heart pump. We will simulate such a hemodynamic FSI process by carrying out Algorithms 1–3 with a stable Stokes-pair (e.g. \(P^2P^1\) or MINI mixed finite element) or a stabilized scheme (e.g. \(P^1P^1\) mixed element with pressure-stabilization).

It is worth mentioning that our numerical simulations are carried out on Tianhe-1A which locates in the National Supercomputing Center, Tianjin, China, and each node of which has two Intel Xeon X5670 CPUs of six cores and 24 GB memory. The mesh we used in the computation has 473,403 vertices, 2,614,905 cells, and 10 boundary layers are attached on the rotor blades and the wall near them. With \(P^2P^1\) mixed finite element, the total number of DOFs is 11,296,137. The time step is chosen to be \(5\times 10^{-5}\) s to reach a rotational steady state at t = 0.13 s. The linear system is solved by the additive Schwartz method, and on each processor, the local linear system is solved with the direct solver MUMPS. As for the timing issue, 1024 cores is used in our computation, and each time step cost about 47 s, the total time of the simulation cost about 33 h.

Considering the Reynolds number of blood flow near the high-speed rotating rotor surface is higher than elsewhere, behaving like a turbulence flow, we then also employ the Reynolds-Averaged Simulation (RAS) turbulence model [12] to replace the laminar fluid model in (1) if \(\omega > 4000\) rpm. We conduct a comparison computation and show a large difference between the RAS model and the laminar model, as elucidated in Fig. 5, the convergence history of the pressure drop from the inlet to the outlet obtained by the RAS model is much stabler while the result of laminar model is oscillatory. Under a prescribed angular velocity \(\omega =6000\) rpm of the rotor and the incoming flow velocity 3 L/min (about 0.2 m/s) at the inlet, we obtain numerical results as shown in Figs. 6 and 7 for the artificial heart pump simulation with both rotating rotor and unmoving guides are interacting with the surrounding blood flow. We can see that the pressure is greatly increased after the blood flow passes through the pump, inducing the blood flow being propelled further forward along with the time marching as shown in Fig. 7 with a reasonable velocity vector field near and behind the high-speed rotating pump rotor.

Fig. 5.
figure 5

Comparison of convergence history of pressure drop between RAS turbulence model and the laminar model with higher angular velocity \(\omega \).

Full size image
Fig. 6.
figure 6

Cross-sectional results of velocity magnitude and pressure in flow direction developing at 0.02 s, 0.09 s and 0.13 s from the left to the right.

Full size image
Fig. 7.
figure 7

The velocity vector field.

Full size image
Fig. 8.
figure 8

Comparisons of velocity magnitude (top) and pressure (bottom) between our method (left) and the commercial package (right).

Full size image
Fig. 9.
figure 9

Comparison of convergence history of pressure drop vs iteration steps between our method and the commercial package (Star-CCM+).

Full size image

To validate our numerical results, we attempt to rebuild the above FSI simulation in a commercial CFD package (Star-CCM+), and compare the obtained numerical results with ours. However, the commercial CFD package cannot deal with the fluid-rotating structure interaction problem with the co-existing rotational fluid/structure regions and stationary fluid/structure regions, which is incomparable with our developed numerical method. In order to compare with the commercial CFD package, we therefore simplify the original setup of the artificial heart pump to let it consist of the rotor only, and remove all unmoving parts of the pump. Then, we let the blood flow inside the entire vascular lumen rotates together with the pump rotor. Figure 8 illustrates that numerical results obtained from both the commercial package and our method are comparable in regard to the profile and the value range of the velocity magnitude and the pressure. Though, a detailed observation over Fig. 8 shows a slight difference between the results of our method and the commercial package, e.g., our velocity is smaller near the head guide and tail guide but larger near the rotor, our pressure is a bit larger near the head guide. That is because some numerical techniques used in our method are still different from the commercial CFD package in many ways. For instance, we use the streamline-diffusion finite element method versus the upwind finite volume method in the CFD package, and a different way to generate the boundary layer meshes, etc.

In addition, we also compare the convergence history of the pressure-drop field with the commercial CFD package, the comparison is shown in Fig. 9 which illustrates that our result matches well with that of the commercial package, moreover, our method even shows a smoother and stabler convergence process. On the other hand, considering that our ALE finite element method is also able to deal with a realistic artificial heart pump that contains both rotating rotor and unmoving head/tail guide and their surrounding blood flow which, however, the commercial CFD package cannot do, we can say that our well developed method is more capable and more powerful in realistic applications to many kinds of fluid-rotating structure interaction problems.

7 Conclusions

We build a Eulerian-Lagrangian model for fluid-structure interaction problem with the rotating elastic body based on the arbitrary Lagrangian-Eulerian approach. The kinematic interface condition is easily dealt with by adopting the velocity as the principle unknown of structure equation and the master-slave relation technique. The dynamic interface condition vanishes in the finite element approximation of a monolithic scheme. Our well developed iterative algorithm demonstrates a satisfactory numerical result for a hemodynamic FSI problem involving a rotating artificial heart pump, where, the RAS turbulence model is employed in our FSI simulation to tackle the turbulent flow behavior near the high-speed rotating rotor. Comparisons with the commercial CFD package validate our numerical results and further our numerical methods, also illustrates that the developed numerical techniques are efficient and flexible to explore the interaction between fluid and a rotational elastic structure.