3.1.1 Acquisition of clinical images

Angiography is an imaging technique successfully used to ‘identify’ the vessel lumen. It exploits the property that a liquid inside the vessel appears brighter than the vessel wall and the surrounding tissue. Angiographies are usually acquired as two-dimensional (2D) images, corresponding to different slices of the domain of interest, but three-dimensional (3D) acquisitions of volumes are also possible.

One of the most common techniques for obtaining an angiography is X-ray imaging, based on the projection of X-ray beams through the body onto suitable screens, and on the contrast produced in the 2D image by the different absorption properties of body structures. To highlight the vessel lumen, a radio-opaque dye is inserted into the bloodstream through the arterial system. To reconstruct tortuous geometries, a rotational angiography (RA) is performed, where X-ray sources and detectors are rapidly rotated around the patient, allowing one to acquire many projections within a few seconds. The excellent spatial resolution of projection angiography (about $0.2~\text{mm}$ , $0.4~\text{mm}$ for RA) makes this technique the gold standard for most vascular imaging applications. Another X-ray angiography technique, widely used for blood flow simulation, is based on computed tomography (CT) technology, where multiple X-ray sources and detectors are rotated rapidly around the patient, allowing one to acquire 3D images with excellent spatial resolution (less than $1~\text{mm}$ in computed tomography angiography, CTA). Unlike projection angiography, another advantage of CTA is the possibility of using intravenous rather than arterial injections. Recently, temporally resolved CTA (4D-CTA) has become feasible. This allows one to obtain several (15–20) 3D images during a heartbeat.

Difficulties may arise in the presence of metal artifacts due to metallic prostheses such as pacemakers, resulting in streaks on the images obscuring anatomical details; see Robertson et al. (Reference Robertson, Yuan, Wang and Vannier1997) and Faggiano, Lorenzi and Quarteroni (Reference Faggiano, Lorenzi and Quarteroni2014) for possible mathematical treatments.

Another widely used technique to obtain angiographies is magnetic resonance imaging (MRI), based on the different decay rates exhibited by body structures on exposure to radio frequency (RF) energy. This is called magnetic resonance angiography (MRA). The generated contrast in the images can be tuned by selecting different RF stimuli. This allows MRA to be suitably tuned to detect soft tissues. Another advantage of MRA is that angiography can be generated without using exogenous agents. However, usually an intravenous injection of a paramagnetic contrast agent is used to improve the blood signal and reduce the acquisition time (contrast-enhanced (CE)-MRA).

Finally, we mention ultrasound (US) imaging, based on the reflections of high-frequency sound waves (a few MHz) transmitted into the body. Ultrasound is the least expensive and invasive of the techniques discussed here, and allows real-time acquisition of 2D images. In contrast, its spatial resolution is the poorest. Recently it has even been possible to acquire 3D images (3D US) by reconstructing a 3D volume from 2D slices.

On the other hand, only a few techniques currently allow us to obtain images of vessel walls. Among these we cite black blood (BB)-MRA, by which the vessel wall and the surrounding tissue can also be viewed, and intravascular ultrasound (IVUS), which is however very invasive since the transducer is placed directly into the artery (typically a coronary artery) via a catheter.

No matter which technique is used, from a mathematical standpoint we can assume that at the end of the acquisition step we obtain a vector $\boldsymbol{I}^{\mathit{clin}}$ , whose $j$ th component, $I_{j}^{\mathit{clin}}$ , corresponds to the intensity of the image at the point $\boldsymbol{x}_{j}$ in grey-scale representing the contrast generated by the imaging technique. The collection of the points $\boldsymbol{x}_{j}$ , $j=1,\ldots ,N^{\mathit{clin}},$ forms the lattice ${\mathcal{L}}^{\mathit{clin}}$ , where $N^{\mathit{clin}}$ is the total number of acquisition points (pixels or voxels) where the image contrast has been evaluated. Here and below, a lattice is a simple collection of points determined by the point coordinates. It can be useful to associate a corresponding image intensity (scalar) function with the image intensity vector $\boldsymbol{I}^{\mathit{clin}}$ , which is typically obtained by interpolation, and will be denoted by $I^{\mathit{clin}}(\boldsymbol{x})$ .

3.1.2 Image enhancement

Medical images are often affected by noise and artifacts that may interfere with the quality of the final results of the preprocessing step. Thus, prior to the reconstruction of the 3D geometry, an imaging enhancement is usually performed.

One popular enhancement technique is resampling, consisting in suitably changing the resolution of the images in one or more directions. In practice, an interpolation of image intensity values $\boldsymbol{I}^{\mathit{clin}}$ onto a more refined lattice is performed. The most commonly used methods are constant interpolation, first-order composite Lagrangian interpolation, B-spline (Unser Reference Unser1999), and windowed sinc interpolation.

The noise in the medical images may be due to thermal effects in the signal processing electronics or to other undesired sources. Reduction of noise could be obtained by means of a smoothing filter, which does not require any prior information about the nature of the noise and has a regularizing effect on the image. This technique is the most commonly used both for CT and MRI images. A very popular filter is the Gaussian filter, consisting in performing a discrete convolution with a Gaussian kernel over the lattice ${\mathcal{L}}^{\mathit{clin}}$ of the image intensity $\boldsymbol{I}^{\mathit{clin}}$ . Unfortunately, together with the noise, smoothing could also filter significant high-frequency image contents. Moreover, since the image is separated from the background by sharp boundaries, characterized by high-frequency content, the smoothing filter could blur and move the boundaries. To prevent this, anisotropic diffusion filtering has been introduced (Perona and Malik Reference Perona and Malik1990): the heat equation is solved for a new image intensity function, with diffusion coefficient decreasing for increasing values of the gradient magnitude of intensity. By so doing, the filtering is not performed at the boundaries where the gradient is large.

Another technique, called multiscale vessel enhancement (Frangi et al. Reference Frangi, Niessen, Hoogeveen, van Walsum and Viergever1999), exploits the specific tubular shape of vascular geometries, and therefore assumes that the minimum modulus eigenvalue of the Hessian matrix of the image intensity function $I^{\mathit{clin}}$ is small, while the other two are large and of equal sign.

At the end of this substep we obtain a new image intensity vector $\boldsymbol{I}^{\mathit{en}}$ whose $j$ th component, $I_{j}^{\mathit{en}}$ , denotes the intensity of the enhanced image in grey-scale at the point $\boldsymbol{x}_{j}$ , $j=1,\ldots ,N^{\mathit{en}}$ , in the lattice ${\mathcal{L}}^{\mathit{en}}$ (and correspondingly an associated enhanced image intensity function $I^{\mathit{en}}(\boldsymbol{x})$ via interpolation). Here, $N^{\mathit{en}}$ is the total number of points in the enhanced image intensity vector. Usually, $N^{\mathit{en}}>N^{\mathit{clin}}$ .

3.1.3 Image segmentation

Image segmentation is the cornerstone of the preprocessing step. It consists in the construction of the shape of a vascular district from the image obtained after the enhancement substep. In particular, the segmentation allows one to detect those points of the lattice ${\mathcal{L}}^{\mathit{en}}$ which – presumably – belong to the boundary of the vessel lumen. The precise definition of the boundary of the lumen is a challenging task which generally requires considerable experience on the part of the user.

The first technique we describe is thresholding, consisting in selecting a threshold $k$ to identify the points $\boldsymbol{x}_{j}\in {\mathcal{L}}^{\mathit{en}}$ such that $I_{j}^{\mathit{en}}>k$ . This is motivated by the assumption that $k$ separates different anatomical structures, in our case the vessel lumen (characterized by intensity values larger than $k$ ) and the background, obtained by the collection of points for which $I_{j}^{\mathit{en}}\leq k$ . The value of $k$ is determined either manually or via a suitable algorithm. In the latter case, one commonly used strategy is full width at half maximum (FWHM), consisting in setting the threshold halfway between the peak intensity within the lumen and the intensity of the background. For the segmentation of special structures, such as calcifications or stents, higher-bound thresholds are used (Boldak, Rolland and Toumoulin Reference Boldak, Rolland and Toumoulin2003).

A more sophisticated class of segmentation methods than thresholding is given by front propagation methods, where the propagation of a suitable wavefront is tracked. The speed of the wave is small in regions where $\boldsymbol{I}^{\mathit{en}}$ changes rapidly and high for other regions, so the wavefront slows down when approaching the boundary. The most popular front propagation method is the fast marching method, which provides an efficient solution to the eikonal problem

where ${\mathcal{D}}^{\mathit{en}}\subset \mathbb{R}^{3}$ is a region that contains all $\boldsymbol{x}_{j}\in {\mathcal{L}}^{\mathit{en}}$ , and where suitable boundary conditions are prescribed on a selected boundary where the propagation starts (Zhu and Tian Reference Zhu and Tian2003). In the above equation, $V$ is the speed of the wavefront and $T(\boldsymbol{x})$ is the first arrival time at point $\boldsymbol{x}$ . In fact, $T$ consists of isocontours, denoting a collection of surfaces describing the shape of the waveform. The vessel boundary is then represented by the points $\boldsymbol{x}_{j}\in {\mathcal{L}}^{\mathit{en}}$ such that $T(\boldsymbol{x}_{j})=T^{b}$ (up to a given tolerance), where $T^{b}$ is a suitable value selected by the user.

Another class of segmentation methods is that of deformable models, where a suitable energy is minimized, allowing the deformation of the body (in our case the boundary of the vessel lumen) to reach a final state with smallest energy, accounting for external terms derived from the image and internal terms constraining the boundary to be regular. The most widely used class of deformable models is the level set method, where a deformable surface is represented implicitly as the zero level of a higher-dimensional embedding function (Sethian Reference Sethian1999). Deformable models, for example those based on cylindrically parametrized surface meshes, incorporate anatomical knowledge of the vessel shape (Frangi et al. Reference Frangi, Niessen, Hoogeveen, van Walsum and Viergever1999, Yim et al. Reference Yim, Cebral, Mullick, Marcos and Choyke2001).

For the segmentation of the vessel wall, Steinman et al. (Reference Steinman, Thomas, Ladak, Milner, Rutt and Spence2001), starting from BB-MRA images, segmented the vessel wall outer boundary using the same deformable model as used for the vessel lumen segmentation. Usually, BB-MRA or other images detecting the vessel wall are not routinely acquired in clinical practice. In this case, a reasonable approach to obtaining the vessel wall is to ‘extrude’ the reconstructed boundary lumen along the outward unit vector by using a suitable function specifying the vessel wall thickness in the different regions of the district of interest.

In those cases where the image intensity vectors $\boldsymbol{I}^{\mathit{clin}}$ and $\boldsymbol{I}^{\mathit{en}}$ refer to 2D slices, application of the above segmentation strategies leads to identification of several vessel boundaries (contours), one for each slice, which now need to be connected to obtain the 3D boundary surface. This operation is called surface reconstruction. A simple procedure is to connect successive contours with straight lines defining surface triangle edges. This strategy is not suitable in the presence of changes of shape such as in bifurcations. Better surface reconstruction is provided by least-squares fitting of polynomial surfaces to the contour set (Wang, Dutton and Taylor Reference Wang, Dutton and Taylor1999). This strategy is suited to managing bifurcations whose branches are fitted separately with a successive extension into the parent vessel. A variant of this approach has been proposed in Geiger (Reference Geiger1993), where contours are first filled with triangles which are then connected to the triangles of the adjacent contours by means of tetrahedra. The final lumen surface is then represented by the boundary of this tetrahedral mesh (formed by triangles). We also mention shape-based interpolation where, for each contour, a characteristic function with positive (resp. negative) values for points located inside (resp. outside) the contour is constructed. The final lumen boundary surface is then represented by the zero level set of the interpolation of all these characteristic functions (Raya and Udupa Reference Raya and Udupa1990). Finally, we briefly describe interpolation by means of radial basis functions (RBFs), which provide a flexible way of interpolating data in multi-dimensional spaces, even for unstructured data where interpolation nodes are scattered and/or do not form a regular grid, and for which it is often impossible to apply polynomial or spline interpolation (Carr, Fright and Beatson Reference Carr, Fright and Beatson1997, Fornefett, Rohr and Stiehl Reference Fornefett, Rohr and Stiehl2001). The coefficients in the linear combination with respect to the RBF basis are determined by solving a suitable linear system, which is invertible under very mild conditions (Peiró et al. Reference Peiró, Formaggia, Gazzola, Radaelli and Rigamonti2007).

A special mention must go to centreline reconstruction. The centreline is a one-dimensional curve centred inside the vessel lumen. Many segmentation tools use the centreline as the starting point, making the assumption that the shape of the section around each centreline location is known (O’Donnell, Jolly and Gupta Reference O’Donnell, Jolly and Gupta1998). Centreline reconstruction allows complete reconstruction of the computational domain when using one-dimensional modelling of blood flow: see Section 4.5.1.

In any case, at the end of the segmentation step we obtain the lattice ${\mathcal{L}}^{\mathit{surf}}$ which collects the points $\boldsymbol{x}_{j}$ , $j=1,\ldots ,N^{\mathit{surf}}$ , classified as belonging to the lumen vessel surface or to the outer wall, where $N^{\mathit{surf}}$ denotes the total number of points of the surface lattice.

3.1.4 Building the computational mesh

Once the final boundary lattice ${\mathcal{L}}^{\mathit{surf}}$ (made up of points on the lumen boundary) is made available, we are ready to build the volumetric mesh ${\mathcal{T}}^{\mathit{vol}}$ in the lumen. This mesh usually consists of unstructured tetrahedra, because of their flexibility in filling volumes of complex shape.

Unstructured volumetric meshes are constructed starting from an analytical expression, say ${\mathcal{S}}(\boldsymbol{x})$ , representing the surface associated with the boundary lattice ${\mathcal{L}}^{\mathit{surf}}$ . This expression can derive from an explicit representation, for instance a bivariate parametric function built as a collection of adjacent polygons. The latter are typically triangles, generated by Lagrangian shape functions, or patches, generated by high-degree polynomials such as NURBS (Sun et al. Reference Sun, Starly, Nam and Darling2005). Alternatively, the surface is represented implicitly as the isosurface of an embedding function. Note that some of the segmentation strategies described above, such as deformable models and those used for the surface reconstruction, directly provide an analytical expression ${\mathcal{S}}(\boldsymbol{x})$ of the lumen boundary surface.

For the construction of unstructured volumetric meshes ${\mathcal{T}}^{\mathit{vol}}$ , we mention two possible approaches. In the first, a boundary surface mesh ${\mathcal{T}}^{\mathit{surf}}$ is generated. To this end, we start from a lattice $\widetilde{{\mathcal{L}}}^{\mathit{surf}}$ (in principle different to ${\mathcal{L}}^{\mathit{surf}}$ ) composed of points of ${\mathcal{S}}$ . Then, the Voronoi diagram for $\widetilde{{\mathcal{L}}}^{\mathit{surf}}$ is constructed. This is a partition of ${\mathcal{S}}$ into non-overlapping regions, each one containing exactly one point (node) of $\widetilde{{\mathcal{L}}}^{\mathit{surf}}$ and composed of all the points of ${\mathcal{S}}$ that are closer to that node than to any other node. Starting from the Voronoi diagram, it is possible to generate a Delaunay mesh ${\mathcal{T}}^{\mathit{surf}}$ : see Thompson, Soni and Weatherill (Reference Thompson, Soni and Weatherill1999). We emphasize that the vertices of the mesh ${\mathcal{T}}^{\mathit{surf}}$ do not necessarily coincide with the points of the lattice $\widetilde{{\mathcal{L}}}^{\mathit{surf}}$ . Popular algorithms for generating a Delaunay mesh have been proposed by Watson (Reference Watson1981) and Weatherill and Hassan (Reference Weatherill and Hassan1994). Once a surface mesh ${\mathcal{T}}^{\mathit{surf}}$ is made available, the volumetric mesh ${\mathcal{T}}^{\mathit{vol}}$ is generated. The latter could be obtained by advancing front methods, where, starting from the triangles of the surface mesh, a front composed of internal nodes is generated. These new nodes allow us to identify tetrahedra, whose validity is verified by checking that they do not intersect the front (Bentley and Friedman Reference Bentley and Friedman1979).

The second approach relies on directly generating the volumetric mesh ${\mathcal{T}}^{\mathit{vol}}$ , for example by means of Delaunay 3D mesh generation, where a starting volumetric lattice $\widetilde{{\mathcal{L}}}^{\mathit{vol}}$ is obtained by locating the nodes in the volume ${\mathcal{V}}(\boldsymbol{x})$ contained in ${\mathcal{S}}(\boldsymbol{x})$ . One of the main problems related to this approach is that boundary meshing is often difficult, since the related surface triangulation could not be of Delaunay type. An alternative approach is given by octree mesh generation, where ${\mathcal{V}}(\boldsymbol{x})$ is embedded in a box and successive subdivisions are performed until the smallest cells permit accurate description of the boundary. Despite being faster, this strategy generates meshes with poor quality near the boundary.

When a volumetric mesh ${\mathcal{T}}^{\mathit{vol}}$ is obtained, a further step (mesh optimization) could be introduced prior to the generation of the final mesh, so as to improve its quality. This prevents mesh distortion, for example the presence of very small angles, which could reduce the convergence of algorithms for the solution of the PDE of interest and thus their accuracy. Mesh optimization is incorporated in the strategies described above; it leads to an optimal mesh, providing the best accuracy for a given number of nodes.

A mesh is deemed valid for blood flow simulations if it allows recovery of outputs of physical interest. In arteries, the mesh should be fine enough to capture wall shear stress (WSS) (Celik et al. Reference Celik, Ghia, Roache, Freitas, Coleman and Raad2008) and, to this end, the construction of a boundary layer mesh is essential, even at low Reynolds numbers (Bevan et al. Reference Bevan, Nithiarasu, van Loon, Sazonov, Luckraz and Garnham2010). Here $\mathit{WSS}$ expresses the magnitude of tangential viscous forces exerted by the fluid on the lumen boundary $\unicode[STIX]{x1D6F4}^{t}$ , defined by

where $\boldsymbol{v}$ is the fluid velocity, $\boldsymbol{n}$ is the outward unit vector, and $\boldsymbol{\unicode[STIX]{x1D70F}}^{(j)}$ , $j=1,2$ , represent the tangential unit vectors. Note that $\mathit{WSS}$ is a scalar function of $\boldsymbol{x}\in \unicode[STIX]{x1D6F4}^{t}$ and $t>0$ . In componentwise notation,

For the structure domain, hexahedral meshing is preferable so as to prevent the locking phenomenon, whereas tetrahedral meshes are used when conforming meshes at the boundary lumen interface are needed, in view of fluid–structure interaction problems (see Section 4.3). Usually, three or four layers of elements are enough to obtain an accurate result (Younis et al. Reference Younis, Kaazempur-Mofrad, Chan, Isasi, Hinton, Chau, Kim and Kamm2004).

For recent reviews on geometric reconstruction for blood flow simulation, see Antiga et al. (Reference Antiga, Peiró, Steinman, Formaggia, Quarteroni and Veneziani2009), Sazonov et al. (Reference Sazonov, Yeo, Bevan, Xie, van Loon and Nithiarasu2011) and Lesage et al. (Reference Lesage, Angelini, Bloch and Funka-Lea2009).

4.2.1 Modelling the structure as a 2D membrane

In some circumstances, because of the thinness of the vessel wall, a non-linear shell model has been proposed: see e.g. Leuprecht et al. (Reference Leuprecht, Perktold, Prosi, Berk, Trubel and Schima2002) and Zhao et al. (Reference Zhao, Xu, Hughes, Thom, Stanton, Ariff and Long2000). In this case, the structure problem is described by two-dimensional equations defined with respect to the middle surface, consisting of the computation of the deformation of this surface.

A simpler equation may be obtained if the structure is modelled as a 2D membrane whose position in space at any time exactly coincides with internal boundary $\unicode[STIX]{x1D6F4}^{t}$ , yielding the so-called generalized string model

(Quarteroni et al. Reference Quarteroni, Tuveri and Veneziani2000c ). Here $\unicode[STIX]{x1D6F4}$ denotes the reference membrane configuration, $d_{r}$ is the radial displacement and $H_{s}$ is the structure thickness. The tensor $\boldsymbol{P}$ accounts for shear deformations and, possibly, for pre-stress,

where $\unicode[STIX]{x1D70C}_{1}(\boldsymbol{x})$ and $\unicode[STIX]{x1D70C}_{2}(\boldsymbol{x})$ are the mean and Gaussian curvatures of $\unicode[STIX]{x1D6F4}$ , respectively (Nobile and Vergara Reference Nobile and Vergara2008), and $f_{s}$ is the forcing term, given by a measurement of the fluid pressure. Equation (4.17) is derived from Hooke’s law for linear elasticity under the assumptions of small thickness, plane stresses, and negligible elastic bending terms (Zienkiewicz and Taylor Reference Zienkiewicz and Taylor2005). To account for the effect of the surrounding tissue, the term $\unicode[STIX]{x1D712}$ in (4.17) needs to be augmented by the elastic coefficient of the tissue $\unicode[STIX]{x1D6FC}_{ST}$ (Formaggia, Quarteroni and Vergara Reference Formaggia, Quarteroni and Vergara2013).

A further simplification arises when $\unicode[STIX]{x1D6F4}$ denotes the lateral surface of a cylinder. By discarding any dependence on the circumferential coordinate, model (4.17) reduces to

where $k$ is the Timoshenko correction factor, $G$ is the shear modulus, $R_{0}$ is the initial cylinder radius and $z$ is the axial coordinate. Often, in the latter case, a viscoelastic term of the form $\unicode[STIX]{x1D6FE}_{v}(\unicode[STIX]{x2202}^{3}d_{r}/\unicode[STIX]{x2202}^{2}z\unicode[STIX]{x2202}t)$ is also added, with $\unicode[STIX]{x1D6FE}_{v}$ denoting a suitable viscoelastic parameter (Quarteroni et al. Reference Quarteroni, Tuveri and Veneziani2000c ).

4.4.1 Conjecturing velocity and pressure profiles

A common trick to effectively prescribe the flow rate condition (3.2) consists in prescribing a velocity profile

where $\boldsymbol{g}=\boldsymbol{g}(t,\boldsymbol{x})$ is chosen in such a way as to satisfy (3.2), that is,

The flow rate condition (3.2) is therefore replaced with the standard (vectorial) Dirichlet condition (4.24). A classical choice for $\boldsymbol{g}$ is a parabolic profile (e.g. for flow simulations in the carotid arteries: Campbell et al. Reference Campbell, Ries, Dhawan, Quyyumi, Taylor and Oshinski2012), a constant profile (often used for the ascending aorta: Moireau et al. Reference Moireau, Xiao, Astorino, Figueroa, Chapelle, Taylor and Gerbeau2012), or that obtained from the Womersley solution (He, Ku and Moore Reference He, Ku and Moore1993). Both the parabolic and Womersley profiles require a circular section to be prescribed, while non-circular sections require an appropriate morphing (He et al. Reference He, Ku and Moore1993).

In spite of its straightforward implementation, this choice has a major impact on the solution, particularly in the neighbourhood of the section $\unicode[STIX]{x1D6E4}^{t}$ and for elevated values of the Reynolds number (Veneziani and Vergara Reference Veneziani and Vergara2007). To reduce the sensitivity to the arbitrary choice of the profile, the computational domain can be artificially elongated by operating what is called a flow extension (Moore, Steinman and Ethier Reference Moore, Steinman and Ethier1997).

A similar approach could be applied to the mean normal Cauchy stress condition (3.3) as well. In the case at hand, we can postulate that the pressure on $\unicode[STIX]{x1D6E4}^{t}$ is constant and that the viscous stress can be neglected, that is, we can prescribe

Note that the above condition in particular satisfies the defective condition (3.3). Condition (4.26) is generally acceptable because the pressure changes in arteries mainly occur along the axial direction and the viscous stresses are negligible on orthogonal cross-sections.

Since $P\boldsymbol{n}$ plays the role of boundary normal Cauchy stress when implemented in the framework of finite element approximations, no further action is required beyond assembling the matrix for Neumann conditions. For this reason, the treatment has been called the ‘do-nothing’ approach (Heywood et al. Reference Heywood, Rannacher and Turek1996). As pointed out by Veneziani Reference Veneziani and Rannacher1998a , Reference Veneziani1998b , this procedure is in fact not completely ‘innocent’. The do-nothing approach corresponds to the following weak formulation (for simplicity we assume homogeneous Dirichlet conditions, $\boldsymbol{v}_{\mathit{up}}=\mathbf{0}$ ): for each $t>0$ , find $\boldsymbol{v}\in \widetilde{\boldsymbol{V}}^{t},\,\boldsymbol{v}=\boldsymbol{v}_{0}$ – chk for $t=0$ in $\unicode[STIX]{x1D6FA}_{f}$ , and $p\in L^{2}(\unicode[STIX]{x1D6FA}_{f}^{t})$ such that

for all

and $q\in L^{2}(\unicode[STIX]{x1D6FA}_{f}^{t})$ .

A do-nothing formulation for the flow rate conditions is possible too: see Heywood et al. (Reference Heywood, Rannacher and Turek1996) and Veneziani (Reference Veneziani1998b ).

Note that as an alternative to (3.3), other defective conditions involving the fluid pressure could be considered as well. This is the case for mean pressure conditions (Heywood et al. Reference Heywood, Rannacher and Turek1996), or conditions involving the total pressure (Formaggia et al. Reference Formaggia, Quarteroni and Vergara2013), defined by $p_{\mathit{tot}}=p+(\unicode[STIX]{x1D70C}_{f}/2)|\boldsymbol{v}|^{2}$ . For a comprehensive review of these conditions, we refer the interested reader to Quarteroni et al. (Reference Quarteroni, Veneziani and Vergara2016c ).

4.4.2 Augmented formulation

An alternative approach is to regard the flow rate boundary condition (3.2) as a constraint for the solution of the fluid problem and then enforce it via a Lagrange multiplier approach. As a scalar constraint, we need a scalar multiplier $\unicode[STIX]{x1D706}=\unicode[STIX]{x1D706}(t)$ at each time, resulting in the following weak formulation (we again consider the case of homogeneous Dirichlet conditions): for each $t>0$ , find $\boldsymbol{v}\in \widetilde{\boldsymbol{V}}^{t},\,\boldsymbol{v}=\boldsymbol{v}_{0}$ for $t=0$ in $\unicode[STIX]{x1D6FA}_{f},\,p\in L^{2}(\unicode[STIX]{x1D6FA}_{f}^{t})$ , and $\unicode[STIX]{x1D706}\in \mathbb{R}$ such that

for all $\boldsymbol{w}\in \widetilde{\boldsymbol{V}}^{t},\,q\in L^{2}(\unicode[STIX]{x1D6FA}_{f}^{t})$ , and $\unicode[STIX]{x1D713}\in \mathbb{R}$ , and where we have set

See Formaggia et al. (Reference Formaggia, Gerbeau, Nobile and Quarteroni2002) and Veneziani and Vergara (Reference Veneziani and Vergara2005), who also analyse the well-posedness of this problem.

Besides prescribing the flow rate condition (3.2), the above augmented formulation enforces at each time a constant-in-space normal Cauchy stress on $\unicode[STIX]{x1D6E4}^{t}$ aligned with its normal direction, which precisely coincides with the Lagrange multiplier $\unicode[STIX]{x1D706}$ , that is,

This method is particularly suitable when the artificial cross-section is orthogonal to the longitudinal axis, so that vector $\boldsymbol{n}$ is truly aligned along the axial direction.

Since the velocity spatial profile is not prescribed a priori, this technique has been used to improve the parabolic-based law implemented in the Doppler technology for the estimation of the flow rate starting from the peak velocity (Ponzini et al. Reference Ponzini, Vergara, Redaelli and Veneziani2006, Vergara et al. Reference Vergara, Ponzini, Veneziani, Redaelli, Neglia and Parodi2010, Ponzini et al. Reference Ponzini, Vergara, Rizzo, Veneziani, Roghi, Vanzulli, Parodi and Redaelli2010).

Extension of the augmented formulation to the case of compliant walls is addressed in Formaggia, Veneziani and Vergara (Reference Formaggia, Veneziani and Vergara2009b ) and to the quasi-Newtonian case in Ervin and Lee (Reference Ervin and Lee2007).

An augmented formulation has been proposed by Formaggia et al. (Reference Formaggia, Gerbeau, Nobile and Quarteroni2002) to prescribe condition (3.3) as well. However, as noticed by Formaggia et al. (Reference Formaggia, Gerbeau, Nobile and Quarteroni2002), in this case it yields at each time the condition

where $\unicode[STIX]{x1D706}=\unicode[STIX]{x1D706}(t)$ is again the constant-in-space Lagrange multiplier. This is a non-homogeneous Dirichlet boundary condition for the fluid velocity, which is generally incompatible with the no-slip condition $\boldsymbol{v}=\boldsymbol{\unicode[STIX]{x1D719}}$ prescribed at the physical boundary $\unicode[STIX]{x1D6F4}^{t}$ . For this reason, the approach is no longer used.

4.4.3 A control-based approach

A different strategy for the fulfilment of condition (3.2) is based on the minimization of the mismatch functional

constrained by the fact that $\boldsymbol{v}$ satisfies the incompressible Navier–Stokes equations (Formaggia, Veneziani and Vergara Reference Formaggia, Veneziani and Vergara2008). This PDE-constrained optimization – which can be regarded as the dual of the above augmented strategy – yields a system of optimality conditions to be fulfilled (also referred to as the Karush–Kuhn–Tucker (KKT) system): see Section 9.1.2 for further details. In particular, Formaggia et al. (Reference Formaggia, Veneziani and Vergara2008) used the normal component of the normal Cauchy stress on $\unicode[STIX]{x1D6E4}^{t}$ as the control variable for the minimization of the mismatch functional. This approach was considered by Formaggia et al. (Reference Formaggia, Veneziani and Vergara2009b ) for the compliant case, whereas Lee (Reference Lee2011), Galvin and Lee (Reference Galvin and Lee2013) and Galvin, Lee and Rebholz (Reference Galvin, Lee and Rebholz2012) address the non-Newtonian, quasi-Newtonian and viscoelastic cases.

The same approach can also be used to fulfil the defective condition (3.3) provided that a suitable functional to be minimized is introduced (Formaggia et al. Reference Formaggia, Veneziani and Vergara2008). This allows us to prescribe (3.3) on a section oblique with respect to the longitudinal axis too. In this case the control variable is the complete normal Cauchy stress vector, that is, the direction of the normal Cauchy stress is also a priori unknown.

Boundary data may also be lacking for the cross-section of the vessel wall. In this case we end up with defective BC issues for the vessel wall: see e.g. Quarteroni et al. (Reference Quarteroni, Veneziani and Vergara2016c ).

4.5.1 The 1D and 0D models

Numerical modelling of the entire cardiovascular system by means of 3D models is currently not affordable because of the complexity of the computational domain, which is composed of thousands of arteries and veins and billions of arterioles, capillaries and venules (Nichols and O’Rourke Reference Nichols and O’Rourke2005). In many applications, reduced-dimensional models are used instead, either as stand-alone models or coupled with the 3D ones.

The first one-dimensional (1D) model was introduced almost 250 years ago by Euler (Reference Euler1775). Subsequently this approach was brought into the engineering environment by Barnard et al. (Reference Barnard, Hunt, Timlake and Varley1966), Hughes (Reference Hughes1974) and Hughes and Lubliner (Reference Hughes and Lubliner1973). These models allow the description of blood flow in a compliant vessel where the only space coordinate is that of the vessel axis.

One-dimensional models may be derived from 3D models by making simplifying assumptions on the behaviour of the flow, the structure, and their interaction (Quarteroni and Formaggia Reference Quarteroni and Formaggia2004, Peiró and Veneziani Reference Peiró, Veneziani, Quarteroni, Formaggia and Veneziani2009). The starting fluid domain is a truncated cone: see Figure 4.2. Referring to cylindrical coordinates $(r,\unicode[STIX]{x1D711},z)$ , we make the following simplifying assumptions: (i) the axis of the cylinder is fixed; (ii) for any $z$ , the cross-section ${\mathcal{S}}(t,z)$ is a circle with radius $R(t,z)$ ; (iii) the solution of both fluid and structure problems does not depend on $\unicode[STIX]{x1D711}$ ; (iv) the pressure is constant over each section ${\mathcal{S}}(t,z)$ ; (v) the axial fluid velocity $v_{z}$ dominates the other velocity components; (vi) only radial displacements are allowed, so the structure deformation takes the form $\boldsymbol{d}=d\boldsymbol{e}_{r}$ , where $\boldsymbol{e}_{r}$ is the unit vector in the radial direction, and $d(t,z)=R(t,z)-R_{0}(z)$ , where $R_{0}(z)$ is the reference radius at the equilibrium; (vii) the fluid is assumed to obey Poiseuille’s law, so that the viscous effects are modelled by a linear term proportional to the flow rate; (viii) the vessel structure is modelled as a membrane with constant thickness.

Figure 4.2. Fluid domain for the derivation of the 1D model.

We introduce the following quantities:

As for the structure and its interaction with the fluid, we need to introduce a membrane law, which in fact prescribes a relation between the pressure and the lumen area (which is determined by $d_{r}$ ) of the form

where $\unicode[STIX]{x1D713}$ is a given function satisfying $\unicode[STIX]{x2202}\unicode[STIX]{x1D713}/\unicode[STIX]{x2202}A>0,\,\unicode[STIX]{x1D713}(A_{0})=0$ . Here $\boldsymbol{\unicode[STIX]{x1D6FD}}$ is a vector of parameters describing the mechanical properties of the membrane.

By integrating over the sections ${\mathcal{S}}$ the momentum fluid equation (4.19a ) in the $z$ -direction and the mass conservation law (4.19b ), we obtain the system

where $\boldsymbol{U}=[A\,\,Q]^{T}$ is the vector of the unknowns,

is the Coriolis coefficient, $K_{r}=-2\unicode[STIX]{x1D70B}\unicode[STIX]{x1D707}s^{\prime }(1)$ is the friction parameter,

while

represent the flux matrix and the dissipation vector term, respectively. A complete derivation of the model can be found in Pedley (Reference Pedley1980), Hughes (Reference Hughes1974) and Peiró and Veneziani (Reference Peiró, Veneziani, Quarteroni, Formaggia and Veneziani2009), for example. Classical choices of the velocity profile $s$ are flat ( $\unicode[STIX]{x1D6FC}=1$ ) and parabolic ( $\unicode[STIX]{x1D6FC}=4/3$ ).

The term $\unicode[STIX]{x2202}A_{0}/\unicode[STIX]{x2202}z$ in $\boldsymbol{B}$ is typically non-positive, accounting for vessel ‘tapering’, that is, reduction of the area of the lumen when proceeding from proximal to distal arteries. The term $\unicode[STIX]{x2202}\boldsymbol{\unicode[STIX]{x1D6FD}}/\unicode[STIX]{x2202}z$ originates from possibly different mechanical properties along the vessel, to describe, for example, the presence of atherosclerotic plaques or vascular prostheses.

If $A>0$ , system (4.30) has two real distinct eigenvalues,

(see e.g. Quarteroni and Formaggia Reference Quarteroni and Formaggia2004), so it is strictly hyperbolic (see e.g. LeVeque Reference Leveque1992). Under physiological conditions, $c_{1}\gg \unicode[STIX]{x1D6FC}\bar{v}$ , yielding $\unicode[STIX]{x1D706}_{1}>0$ and $\unicode[STIX]{x1D706}_{2}<0$ , we thus have two waves travelling in opposite directions, associated with corresponding characteristic variables. An explicit expression of these variables as a function of the physical variables could, in general, be derived, that is,

A simple membrane law (4.29) can be obtained by the algebraic relation

(Formaggia, Lamponi and Quarteroni Reference Formaggia, Lamponi and Quarteroni2003, Formaggia et al. Reference Formaggia, Quarteroni and Vergara2013), where $\unicode[STIX]{x1D708}$ is the Poisson modulus of the membrane, $E$ is its Young’s modulus, and $H_{s}$ is its thickness, yielding

This simple law, stating that the membrane radial displacement $d_{r}$ is linearly proportional to the fluid pressure, has been considered successfully in many applications: see e.g. Steele et al. (Reference Steele, Wan, Ku, Hughes and Taylor2003), Matthys et al. (Reference Matthys, Alastruey, Peiró, Khir, Segers, Verdonck, Parker and Sherwin2007) and Grinberg et al. (Reference Grinberg, Cheever, Anor, Madsen and Karniadakis2010). Other laws have been proposed to account for additional features of arterial walls, such as viscoelasticity, wall-inertia and longitudinal pre-stress (Quarteroni et al. Reference Quarteroni, Tuveri and Veneziani2000c , Formaggia et al. Reference Formaggia, Lamponi and Quarteroni2003).

A further geometric reduction is represented by the so-called lumped parameter models, which are zero-dimensional (0D) models obtained by integrating the 1D problem over the axial direction. These are typically used to describe the peripheral part of the arterial and venous tree, such as the capillaries and the arterioles.

In this case only dependence on time is allowed, and nominal values of the unknowns are used as representative of the entire compartment. To this end, we introduce the average flow rate and pressure in the district at hand, defined by

respectively, where $z_{p}$ and $z_{d}$ are the proximal and longitudinal abscissas of the segment, respectively, $l$ is its length and $V$ its volume. The convective term is dropped since in the peripheral sites the velocity is small.

If we take the longitudinal average of the momentum equation given by the first of (4.30) and combine it with (4.34), we obtain the ordinary differential equation (ODE)

where $\widehat{P}_{d}$ and $\widehat{P}_{p}$ are the distal and proximal pressure, respectively. When taking the longitudinal average of the mass conservation law given by the second of (4.30) and using (4.34), we obtain

where $\widehat{Q}_{d}$ and $\widehat{Q}_{p}$ are the distal and proximal flow rate, respectively. See Peiró and Veneziani (Reference Peiró, Veneziani, Quarteroni, Formaggia and Veneziani2009).

The two ODEs (4.35, 4.36) can be regarded as the starting point for a 0D description of a compartment model of an arterial tract. In fact, the term $L(\text{d}\widehat{Q}/\text{d}t)$ , with $L=\unicode[STIX]{x1D70C}_{f}l/A_{0}$ , corresponds to the blood acceleration, $R\widehat{Q}$ , with $R=\unicode[STIX]{x1D70C}_{f}K_{R}l/A_{0}^{2}$ , stems from the blood resistance due to the viscosity, while $C(\text{d}\widehat{P}/\text{d}t)$ , with $C=\sqrt{A_{0}}l/\unicode[STIX]{x1D702}$ , is due to the compliance of the vessel wall. Usually an electrical analogy is used to easily interpret 0D models. In particular, the flow rate plays the role of the current whereas the pressure is the potential. Accordingly, the acceleration term is represented by an inductance, the viscosity term by a resistance, and the compliance term by a capacitance.

To close the system (4.35, 4.36) (featuring four unknowns) we also need to include the boundary conditions originally prescribed on the 1D model. For instance, we can assume a Dirichlet condition at the inlet and a Neumann condition at the outlet. Thus, we may localize the unknown pressure $\widehat{P}$ at the proximal section ( $\widehat{P}\approx \widehat{P}_{p}$ ), assuming that the distal pressure $\widehat{P}_{d}$ is given. Similarly we assume that the flow rate $\widehat{Q}$ is approximated by $\widehat{Q}_{d}$ and that the proximal flow rate $\widehat{Q}_{p}$ is given. Then, from (4.35, 4.36), we obtain

corresponding to the electrical circuit drawn in Figure 4.3. Other sequences corresponding to different boundary conditions and then to different state variables are also possible (see e.g. Quarteroni et al. Reference Quarteroni, Veneziani and Vergara2016c ). Even though these schemes are equivalent in terms of functionality, they play a different role when coupled with higher-dimensional models (see e.g. Quarteroni et al. Reference Quarteroni, Veneziani and Vergara2016c ).

Figure 4.3. Example of lumped parameter scheme for an arterial tract.

For a description of more complex vascular districts, we may combine several 0D elementary tracts, by gluing them via classical continuity arguments. However, the lumped parameter models have mainly been used to provide suitable boundary conditions at the distal artificial sections of 3D and 1D models. In this case, one simple compartment is enough to describe the entire arterial system downstream of the region of interest. Examples are provided by the windkessel model (Westerhof, Lankhaar and Westerhof Reference Westerhof, Lankhaar and Westerhof2009), featuring an average resistance and capacitance, the $3$ -element windkessel model (Westerhof et al. Reference Westerhof, Lankhaar and Westerhof2009), where a second resistance is added before the windkessel compartment, and the $4$ -element windkessel model (Stergiopulos et al. Reference Stergiopulos, Westerhof, Meister and Westerhof1996, Stergiopulos, Westerhof and Westerhof Reference Stergiopulos, Westerhof and Westerhof1999), where an inductance element is added to the 3-element windkessel model. A 0D model given simply by a resistance is used to provide absorbing boundary conditions at the outlets of the fluid domain in FSI simulations: see e.g. Nobile and Vergara (Reference Nobile and Vergara2008). Instead, more sophisticated approaches account for the propagative dynamics associated with the peripheral circulation, such as the structured tree model (Olufsen et al. Reference Olufsen, Peskin, Kim, Pedersen, Nadim and Larsen2000), which assumes an asymmetric self-similar structure for the peripheral network.

4.5.2 Geometric multiscale coupling

The geometric multiscale approach, first introduced in Quarteroni and Veneziani (Reference Quarteroni, Veneziani and Bristeau1997), consists in the coupling between the 3D, 1D and 0D models. The idea is to use higher-dimensional models in those regions where a very detailed description is required, and lower-dimensional models in the remaining part of the region of interest. This allows us to describe much of the circulatory system.

As discussed earlier, 0D models are typically used to provide boundary conditions for 3D and 1D models. For this reason, the coupling between 3D or 1D models with an extended 0D network has rarely been considered in applications. Instead, the 3D/1D coupling has received a great deal of attention. For this reason, we will detail only the latter case, while referring to Quarteroni, Ragni and Veneziani (Reference Quarteroni, Ragni and Veneziani2001), Quarteroni and Veneziani (Reference Quarteroni and Veneziani2003), Vignon-Clementel et al. (Reference Vignon-Clementel, Figueroa, Jansen and Taylor2006), Kim et al. (Reference Kim, Vignon-Clementel, Figueroa, Ladisa, Jansen, Feinstein and Taylor2009), Migliavacca et al. (Reference Migliavacca, Balossino, Pennati, Dubini, Hsia, de Leval and Bove2006) and Haggerty et al. (Reference Haggerty, Mirabella, Restrepo, de Zélicourt, Rossignac, Sotiropoulos, Spray, Kanter, Fogel, Yoganathan, Holzapfel and Kuhl2013) for the 3D/0D coupling, and Formaggia et al. (Reference Formaggia, Nobile, Quarteroni and Veneziani1999) and Fernández, Milisic and Quarteroni (Reference Fernández, Milisic and Quarteroni2005) for the 1D/0D coupling.

As shown in Figure 4.4, we consider a 3D/FSI problem (4.19) in a 3D cylindrical domain together with initial conditions and boundary conditions at the proximal boundaries and at the external structure.

Figure 4.4. Schematic representation of the reference 3D/1D coupled model.

At the distal boundaries, the 3D problem is coupled with the 1D model (4.30) written in the domain $z\in [0,L]$ , together with initial conditions and a boundary condition at the distal boundary. We let $\unicode[STIX]{x1D6E4}^{t}=\unicode[STIX]{x1D6E4}_{f}^{t}\cup \unicode[STIX]{x1D6E4}_{s}^{t}$ be the coupling interface from the 3D side, which corresponds to the point $z=0$ from the 1D side (see Figure 4.4).

A major mathematical issue is how to couple 3D and 1D models at the common interface. Several strategies can be pursued, yielding many (alternative) sets of interface conditions. For a rigorous derivation of the 3D/1D problem and a detailed discussion of the interface conditions, we refer the interested reader to the recent review article by Quarteroni et al. (Reference Quarteroni, Veneziani and Vergara2016c ).

Here we will follow the guiding principle described in Formaggia et al. (Reference Formaggia, Quarteroni and Vergara2013), where suitable interface conditions are derived from a global energy estimate. In particular, we introduce, together with the 3D energy (4.21), the 1D energy

(Formaggia et al. Reference Formaggia, Gerbeau, Nobile and Quarteroni2001), where $\unicode[STIX]{x1D712}(A)=\int _{A_{0}}^{A}\unicode[STIX]{x1D713}(\unicode[STIX]{x1D70F})\text{d}\unicode[STIX]{x1D70F}$ and $\unicode[STIX]{x1D713}$ is the vessel law (see (4.29)). Note that the stand-alone 1D problem satisfies bounds for this energy functional as proved in Formaggia et al. (Reference Formaggia, Gerbeau, Nobile and Quarteroni2001).

Let $P_{\mathit{tot}}=\unicode[STIX]{x1D713}(A)+(\unicode[STIX]{x1D70C}_{f}/2)\bar{v}^{2}$ be the total pressure for the 1D model, and let $p_{\mathit{tot}}=p+(\unicode[STIX]{x1D70C}_{f}/2)|\boldsymbol{v}|^{2}$ be that of the 3D model. For the interface coupling conditions holding at $\unicode[STIX]{x1D6E4}^{t}$ , let us assume that the following inequality holds:

Then, for all $t>0$ , the coupled 3D/1D problem (4.19, 4.30) with homogeneous boundary conditions satisfies the energy decay property

For the proof see Formaggia et al. (Reference Formaggia, Quarteroni and Vergara2013).

The above result provides an indication of how to find suitable interface conditions for the 3D/1D coupled problem. In particular, for inequality (4.38) to be fulfilled it is sufficient that the interface conditions

hold for the fluid, together with

for the structure (Formaggia, Moura and Nobile Reference Formaggia, Moura and Nobile2007, Formaggia et al. Reference Formaggia, Quarteroni and Vergara2013). Similarly, inequality (4.38) holds if relation (4.40) is replaced by

The interface conditions (4.39) prescribe the continuity of the flow rate (kinematic condition) and a dynamic condition involving the total pressures. Note that (4.40) and (4.41) are in fact independent of the 1D model, resulting in boundary conditions for the 3D structure problems only. This allows a discontinuity to manifest in the displacement between the 3D and 1D models.

Usually, dynamic interface conditions involving the pressure (instead of the total pressure) are considered in place of (4.39b ), such as

However, the above condition does not satisfy the compatibility condition (4.38). More precisely, in this case we have

Even though the right-hand side is not necessarily (always) negative, numerical evidence indicates that condition (4.42) leads to stable results for haemodynamic applications: see Malossi et al. (Reference Malossi, Blanco, Crosetto, Deparis and Quarteroni2013). This interface condition is indeed the most commonly used among dynamic ones.

4.6.1 Numerical methods for the fluid problem

We start by reviewing some numerical methods for the fluid problem (4.1) together with its initial and boundary conditions.

As for the time discretization, we usually consider implicit methods with a semi-implicit treatment of the convective term and (for a moving domain) of the fluid domain. The problem is solved at discrete time $t^{n+1}$ in the domain $\unicode[STIX]{x1D6FA}_{f}^{\ast }$ and with convective term $\unicode[STIX]{x1D70C}_{f}(\boldsymbol{v}^{\ast }\cdot \unicode[STIX]{x1D6FB})\boldsymbol{v}^{n+1}$ , where $\unicode[STIX]{x1D6FA}_{f}^{\ast }$ and $\boldsymbol{v}^{\ast }$ are suitable extrapolations of $\unicode[STIX]{x1D6FA}_{f}^{n+1}$ and $\boldsymbol{v}^{n+1}$ of the same order as the time discretization. This choice introduces a CFL-like (Courant–Friedrichs–Lewy) restriction on the time step to preserve absolute stability (Quarteroni and Valli Reference Quarteroni and Valli1994). However, this condition is very mild in haemodynamic applications, since, for accuracy purposes, the pulsatility of the blood signal and the quick dynamics around systole can only be accommodated by choosing a small $\unicode[STIX]{x1D6E5}t$ . Usually, a second-order approximation is considered a good choice in haemodynamics; in this respect, the second-order backward difference formula (BDF2) and Crank–Nicolson are the most widely used methods (Quarteroni, Sacco and Saleri Reference Quarteroni, Sacco and Saleri2000a ).

The first class of methods we present is based on a decomposition of the semi-discrete problem at the spatial continuous level (differential splitting or projection methods). The basic idea underlying these methods is to split the computation of velocity and pressure, with a final step aiming to recover the incompressibility constraint. In what follows we detail the Chorin–Teman method (Chorin Reference Chorin1968, Temam Reference Temam1969), originally proposed for homogeneous Dirichlet conditions and fixed domain, which is the progenitor of these methods. We only detail the case of the backward Euler discretization.

Chorin–Temam method

For $n\geq 0$ , at time $t^{n+1}$ :

1. (1) Solve in $\unicode[STIX]{x1D6FA}_{f}$ the advection-reaction-diffusion problem with homogeneous Dirichlet condition for the intermediate unknown velocity $\widetilde{\boldsymbol{v}}^{n+1}$ :

   $$\begin{eqnarray}\unicode[STIX]{x1D70C}_{f}\frac{\widetilde{\boldsymbol{v}}^{n+1}-\boldsymbol{v}^{n}}{\unicode[STIX]{x1D6E5}t}-\unicode[STIX]{x1D707}(\unicode[STIX]{x1D6FB}\widetilde{\boldsymbol{v}}^{n+1}+(\unicode[STIX]{x1D6FB}\widetilde{\boldsymbol{v}}^{n+1})^{T})+\unicode[STIX]{x1D70C}_{f}(\boldsymbol{v}^{n}\cdot \unicode[STIX]{x1D6FB})\widetilde{\boldsymbol{v}}^{n+1}=\mathbf{0}.\end{eqnarray}$$

2. (2) Solve in $\unicode[STIX]{x1D6FA}_{f}$ the pressure problem with homogeneous Neumann conditions:

   $$\begin{eqnarray}\triangle p^{n+1}=\frac{\unicode[STIX]{x1D70C}_{f}}{\unicode[STIX]{x1D6E5}t}\unicode[STIX]{x1D6FB}\cdot \widetilde{\boldsymbol{v}}^{n+1}.\end{eqnarray}$$

3. (3) Correct the velocity:

   $$\begin{eqnarray}\boldsymbol{v}^{n+1}=\widetilde{\boldsymbol{v}}^{n+1}-\frac{\unicode[STIX]{x1D6E5}t}{\unicode[STIX]{x1D70C}_{f}}\unicode[STIX]{x1D6FB}p^{n+1}.\end{eqnarray}$$

This splitting method is based on the Ladyzhenskaya theorem (Girault and Raviart Reference Girault and Raviart1986), which states that a vector function belonging to $[L^{2}(\unicode[STIX]{x1D6FA}_{f})]^{3}$ can always be decomposed as the sum of a solenoidal part and a gradient term. In fact, the correction step corresponds to projecting the intermediate velocity onto

Thus, it is possible to show that $\boldsymbol{v}^{n+1}$ and $p^{n+1}$ are in fact solutions of the original semi-discrete problem. The Chorin–Temam method is very effective since it overcomes the saddle-point nature of the problem and solves two standard uncoupled elliptic problems. However, it suffers from inaccuracies at the boundary. In particular, the tangential velocity cannot be controlled (see the definition of $\boldsymbol{H}$ ) and spurious pressure values appear as a consequence of the artificial Neumann condition for the pressure problem (Rannacher Reference Rannacher1992). This has an effect on the accuracy of the semi-discrete solution. In particular, the following error estimate holds true (Rannacher Reference Rannacher1992):

The use of higher-order time approximations leads to the same accuracy.

An improvement of the above method is given by the rotational incremental variant of the Chorin–Temam scheme (Timmermans, Minev and van de Vosse Reference Timmermans, Minev and van de Vosse1996). Below we detail the case of BDF2 and second-order extrapolation of the convective term, since the first-order approximation does not lead to any improvement.

Rotational–incremental Chorin–Temam method

For $n\geq 0$ , at time $t^{n+1}$ :

1. (1) Solve in $\unicode[STIX]{x1D6FA}_{f}$ the advection-reaction-diffusion problem with homogeneous Dirichlet condition in the intermediate unknown velocity $\widetilde{\boldsymbol{v}}^{n+1}$ :

   $$\begin{eqnarray}\displaystyle & & \displaystyle \displaystyle \unicode[STIX]{x1D70C}_{f}\frac{3\widetilde{\boldsymbol{v}}^{n+1}-4\boldsymbol{v}^{n}+\boldsymbol{v}^{n-1}}{2\unicode[STIX]{x1D6E5}t}-\unicode[STIX]{x1D707}(\unicode[STIX]{x1D6FB}\widetilde{\boldsymbol{v}}^{n+1}+(\unicode[STIX]{x1D6FB}\widetilde{\boldsymbol{v}}^{n+1})^{T})\nonumber\\ \displaystyle & & \displaystyle \qquad +\,\unicode[STIX]{x1D70C}_{f}((2\boldsymbol{v}^{n}-\boldsymbol{v}^{n-1})\cdot \unicode[STIX]{x1D6FB})\widetilde{\boldsymbol{v}}^{n+1}+\unicode[STIX]{x1D6FB}p^{n}=\mathbf{0}.\nonumber\end{eqnarray}$$

2. (2) Solve the pressure problem:

   $$\begin{eqnarray}\displaystyle \begin{array}{@{}ll@{}}\displaystyle \triangle p^{n+1}=\triangle p^{n}+\biggl(\frac{3\unicode[STIX]{x1D70C}_{f}}{2\unicode[STIX]{x1D6E5}t}-\unicode[STIX]{x1D707}\biggr)\unicode[STIX]{x1D6FB}\cdot \widetilde{\boldsymbol{v}}^{n+1}, & \displaystyle \boldsymbol{x}\in \unicode[STIX]{x1D6FA}_{f},\\ \displaystyle \frac{\unicode[STIX]{x2202}p^{n+1}}{\unicode[STIX]{x2202}\boldsymbol{n}}=\unicode[STIX]{x1D707}(\unicode[STIX]{x1D6FB}\times \unicode[STIX]{x1D6FB}\times \widetilde{\boldsymbol{v}}^{n+1})\cdot \boldsymbol{n}, & \displaystyle \boldsymbol{x}\in \unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}_{f}.\end{array} & & \displaystyle \nonumber\end{eqnarray}$$

3. (3) Correct the velocity:

   $$\begin{eqnarray}\boldsymbol{v}^{n+1}=\widetilde{\boldsymbol{v}}^{n+1}-\frac{2}{3}\frac{\unicode[STIX]{x1D6E5}t}{\unicode[STIX]{x1D70C}_{f}}\unicode[STIX]{x1D6FB}(p^{n+1}-p^{n}+\unicode[STIX]{x1D707}\unicode[STIX]{x1D6FB}\cdot \widetilde{\boldsymbol{v}}^{n+1}).\end{eqnarray}$$

Unlike the classical Chorin–Temam scheme, in the above method the boundary conditions for the pressure problem are consistent, and no numerical boundary layer for the pressure is observed. This is confirmed by the improved error estimate

(Guermond and Shen Reference Guermond and Shen2006).

The above two methods belong to the general class of pressure-correction methods; see also Codina and Badia (Reference Codina and Badia2006), Guermond, Minev and Shen (Reference Guermond, Minev and Shen2006) and Guermond and Quartapelle (Reference Guermond and Quartapelle1998), for example. A different class is obtained by switching the role of velocity and pressure in the splitting, that is, the viscous term is now ignored or treated explicitly in the first step and the velocity is then corrected accordingly (velocity-correction schemes: Orszag, Israeli and Deville Reference Orszag, Israeli and Deville1986, Karniadakis, Israeli and Orszag Reference Karniadakis, Israeli and Orszag1991). These schemes are characterized by the same non-optimal error estimates as the pressure-correction schemes due to artificial Neumann conditions for the pressure problem. Again, an improvement could be obtained by considering a rotational-incremental variant (Guermond and Shen Reference Guermond and Shen2003).

In haemodynamics it is often the case that Neumann boundary conditions are prescribed at some artificial section. The extension of the differential splitting methods to this case is addressed in Guermond, Minev and Shen (Reference Guermond, Minev and Shen2005): on the Neumann boundary we have an artificial Dirichlet condition for the pressure, which again diminishes the optimal rate of convergence with respect to $\unicode[STIX]{x1D6E5}t$ .

Given the next methods we are going to review, it is convenient to introduce the algebraic problem arising from the application of a Galerkin-like method to the semi-discrete-in-time problem. First of all we notice that the solvability of the discretized-in-space problem is guaranteed by a suitable compatible choice of the approximation spaces for the velocity and the pressure in order to satisfy the discrete inf-sup stability condition (Quarteroni and Valli Reference Quarteroni and Valli1994). As is well known, an example of finite elements for a tetrahedral mesh is provided by piecewise polynomials of order 2 for the velocity approximation and of order 1 for the pressure approximation. This choice guarantees the existence and uniqueness of the solution to the linearized fully discrete problem and is often used to provide a stable solution in haemodynamics. In this case, we have the optimal error estimate

provided that $\boldsymbol{v}^{n}$ and $p^{n}$ are sufficiently regular. (For other stable choices see Quarteroni and Valli Reference Quarteroni and Valli1994 and Boffi, Brezzi and Fortin Reference Boffi, Brezzi and Fortin2013.) Alternatively, suitable stabilization terms could be added to the problem, circumventing the inf-sup condition and allowing the use of polynomials of equal order. In this case, additional terms are added to the mass conservation equation and, if needed, to the momentum conservation equation. Usually, these techniques also stabilize convected-dominated problems arising when the Reynolds number is high, for example in the aorta or in stenotic carotid arteries. One technique is streamline upwind/Petrov–Galerkin (SUPG) (Muller et al. Reference Muller, Sahni, Lia, Jansen, Shephard and Taylor2005). A generalization of SUPG is the variational multiscale (VMS) method (Hughes Reference Hughes1995, Hughes, Mazzei and Jansen Reference Hughes, Mazzei and Jansen2000), which is based on decomposition of the unknown into two terms, one accounting for the large scales and the other for the small scales. The same decomposition is used for the test functions, so that a system of two coupled problems is obtained. The VMS method is also useful since it allows modelling of transition to turbulence effects, which may occur in some pathological conditions such as stenoses (Akkerman et al. Reference Akkerman, Bazilevs, Calo, Hughes and Hulshoff2008). For an application to haemodynamics, see Forti and Dede’ (Reference Forti and Dede’2015).

Here we introduce the algebraic problem related to the fully discretized linearized problem. For the sake of exposition, we limit ourselves to the cases without stabilization terms. For the more general case, we refer the interested reader to Elman, Silvester and Wathen (Reference Elman, Silvester and Wathen2005) and Benzi, Golub and Liesen (Reference Benzi, Golub and Liesen2005), for example. At each time step we have

where $\boldsymbol{V}$ and $\boldsymbol{P}$ are the vectors collecting the velocity and pressure unknowns, $A=\unicode[STIX]{x1D70C}_{f}(\unicode[STIX]{x1D6FC}/\unicode[STIX]{x1D6E5}t)M_{f}+\unicode[STIX]{x1D70C}_{f}N(\boldsymbol{V}^{\ast })+\unicode[STIX]{x1D707}K$ ( $M_{f}$ is the mass matrix, $N$ is the matrix related to the linearized convective term, $K$ is the stiffness matrix), $\boldsymbol{F}_{f}$ accounts for non-homogeneous Dirichlet and Neumann conditions and the terms coming from time discretization, $\unicode[STIX]{x1D6FC}$ depends on the time discretization scheme, and where we have omitted the current temporal index  $^{n+1}$ .

The above non-symmetric linear system can be solved by means of a Krylov method, for example the GMRES method. Suitable preconditioners are required. A classical choice is given by block preconditioners, which again split the solution of the velocity and of the pressure:

If $P_{A}=A$ and $P_{\unicode[STIX]{x1D6F4}}=\unicode[STIX]{x1D6F4}=BA^{-1}B^{T}$ (the Schur complement), then the solution is achieved in three GMRES iterations (Elman et al. Reference Elman, Silvester and Wathen2005). In fact, this choice is equivalent to formally solving the momentum equation for the velocity and substituting its expression in the mass equation. However, in practice this preconditioner is not efficient, since the linear system involving the Schur complement is too onerous, $\unicode[STIX]{x1D6F4}$ being a full matrix whose explicit construction requires knowledge of $A^{-1}$ . Efficient preconditioners can be obtained by approximating $\unicode[STIX]{x1D6F4}$ (and, if needed, $A$ ). For low Reynolds numbers (less than 10, say), an effective choice is given by $P_{\unicode[STIX]{x1D6F4}}=\unicode[STIX]{x1D707}^{-1}M_{P}$ , where $M_{P}$ is the pressure mass matrix (or even its diagonal: Elman and Silvester Reference Elman and Silvester1996). Therefore this is a good choice in haemodynamics for small vessels. For increasing Reynolds numbers, the convergence properties of this preconditioner deteriorates since it does not account for the convective term. A better choice for medium and large vessels is given by $P_{\unicode[STIX]{x1D6F4}}=A_{P}F_{P}^{-1}M_{P}$ , where $A_{P}$ is the pressure stiffness matrix, and $F_{P}=\unicode[STIX]{x1D707}A_{P}+\unicode[STIX]{x1D70C}_{f}N_{P}(\boldsymbol{V}^{\ast })$ , where $N_{P}$ is the matrix related to the convective term defined on the pressure space (pressure convection-diffusion preconditioner: Elman et al. Reference Elman, Silvester and Wathen2005, Benzi et al. Reference Benzi, Golub and Liesen2005). For the solution of the velocity problem, suitable preconditioners for the advection-reaction-diffusion problem could be introduced. Alternatively, fast solutions such as V-cycle multigrid can be considered (Turek Reference Turek1999).

Another class of preconditioners is obtained by an inexact block LU factorization of the fluid matrix. The starting point is the exact factorization

Again, different preconditioners are obtained by suitable approximations $\widehat{A}_{1}$ and $\widehat{A}_{2}$ of $A$ and $\widehat{\unicode[STIX]{x1D6F4}}$ of $\unicode[STIX]{x1D6F4}$ , leading to

A convenient choice is $\widehat{A}_{1}=\widehat{A}_{2}=D_{A}$ and $\widehat{\unicode[STIX]{x1D6F4}}=BD_{A}^{-1}B^{T}$ , where $D_{A}$ is the diagonal of $A$ (SIMPLE preconditioner: Patankar and Spalding Reference Patankar and Spalding1972, Li and Vuik Reference Li and Vuik2004). Note that in this case $\widehat{\unicode[STIX]{x1D6F4}}$ is sparse and could be explicitly assembled. This is an effective choice when the fluid matrix is diagonally dominant, that is, when small values of $\unicode[STIX]{x1D6E5}t$ are used. Another choice is the Yosida preconditioner, where

(Veneziani Reference Veneziani2003). Again, the efficiency deteriorates for increasing $\unicode[STIX]{x1D6E5}t$ . The Yosida preconditioner was originally introduced as a solver in Veneziani (Reference Veneziani1998b ); see also Quarteroni, Saleri and Veneziani Reference Quarteroni, Saleri and Veneziani1999, Reference Quarteroni, Saleri and Veneziani2000b . This led to a splitting of the velocity and of the pressure computation, which could be seen as the algebraic counterpart of the Chorin–Temam method (algebraic pressure-correction methods). In particular, we have the following steps:

Again, an incremental version of the algebraic pressure-correction methods could also be considered (Quarteroni et al. Reference Quarteroni, Saleri and Veneziani2000b ). An extension to spectral methods is provided in Gervasio, Saleri and Veneziani (Reference Gervasio, Saleri and Veneziani2006).

For a recent comparison of the performance of different preconditioners (all described above) used for haemodynamic applications, see Deparis, Grandperrin and Quarteroni (Reference Deparis, Grandperrin and Quarteroni2014).

In healthy conditions, blood flow is mainly laminar. Transitional flow may develop in some pathological instances, or with the assistance of devices. In these circumstances, suitable mesh refinement is often employed, possibly accompanied by the use of turbulence models. We mention the case of stenotic carotid arteries, for which Stroud, Berger and Saloner (Reference Stroud, Berger and Saloner2002) and Grinberg, Yakhot and Karniadakis (Reference Grinberg, Yakhot and Karniadakis2009) use Reynolds-averaged Navier–Stokes (RANS) models, Lee et al. (Reference Lee, Lee, Fischer, Bassiouny and Loth2008), Fischer et al. (Reference Fischer, Loth, Lee, Lee, Smith and Bassiouny2007) and Cheung et al. (Reference Cheung, Wong, Yeoh, Yang, Tu, Beare and Phan2010) use direct numerical simulation (DNS), and Rayz, Berger and Saloner (Reference Rayz, Berger and Saloner2007) and Lancellotti et al. (Reference Lancellotti, Vergara, Valdettaro, Bose and Quarteroni2015) use large eddy simulation (LES). Bazilevs et al. (Reference Bazilevs, Gohean, Hughes, Moser and Zhang2009) use the VMS formulation to describe transitional effects in the ascending aorta under the influence of the left ventricular assist device (LVAD).

In Figure 4.5 we show some examples of numerical results obtained in four real geometries reconstructed from radiological images (see the caption for details). These results highlight the complex pattern of blood flow induced by the geometry and by the heart pulsatility. To highlight the transitional effects in stenotic carotid arteries, Figure 4.5(c) plots the $Q$ criterion, defined by

where $\boldsymbol{S}=\unicode[STIX]{x1D6FB}\boldsymbol{u}+(\unicode[STIX]{x1D6FB}\boldsymbol{u})^{T}$ and $\boldsymbol{\unicode[STIX]{x1D6FA}}=\unicode[STIX]{x1D6FB}\boldsymbol{u}-(\unicode[STIX]{x1D6FB}\boldsymbol{u})^{T}$ (Lee et al. Reference Lee, Lee, Fischer, Bassiouny and Loth2008). Positive values of $Q$ indicate locations where rotations dominate strain and shear.

Figure 4.5. (a) Velocity vectors in the aneurysm of an abdominal aorta (CT images from the Vascular Surgery and Radiology Divisions at Fondazione IRCSS Cà Granda, Ospedale Maggiore Policlinico, Milan, Italy). (b) Velocity streamlines in a stenotic carotid artery (MRI images from the Vascular Surgery and Radiology Divisions at Ospedale Maggiore Policlinico, Milan). (c) Coherent vortical structures by $Q$ criterion in a stenotic carotid artery (we show only the regions with $Q>50\,000$ shaded by the velocity magnitude; CT images from the Vascular Surgery and Radiology Divisions at Ospedale Maggiore Policlinico, Milan). (d) Wall shear stress in an ascending aorta (MRI images from the Cardio Surgery and Radiology Divisions at Ospedale Borgo Trento, Verona, Italy). Numerical results were obtained using the finite element library LifeV, P2/P1 finite elements, the backward Euler scheme for the time discretization with a semi-implicit treatment of the non-linear term, and the Yosida preconditioner. For the stenotic carotid arteries an LES model has been used.

4.6.2 Numerical methods for the vessel wall problem

In this section we review some of the most commonly used numerical approaches for the solution of problem (4.7), endowed with its initial and boundary conditions.

For the time discretization, a popular family of schemes is that of Newmark, which is characterized by two parameters $\unicode[STIX]{x1D703}$ and $\unicode[STIX]{x1D701}$ (Newmark Reference Newmark1959). The special combination $\unicode[STIX]{x1D703}=1/2$ and $\unicode[STIX]{x1D701}=1/4$ yields the following semi-discrete form of (4.16) (for simplicity we set $P_{\mathit{ext}}=0$ and $\boldsymbol{h}_{s}=\mathbf{0}$ ):

where $\widehat{\boldsymbol{w}}^{n+1}$ and $\widehat{\boldsymbol{a}}^{n+1}$ denote approximations of vessel wall velocity and acceleration, respectively. This method is unconditionally absolutely stable and second-order accurate with respect to $\unicode[STIX]{x1D6E5}t$ . An extension of Newmark schemes is provided by the generalized-alpha method (Chung and Hulbert Reference Chung and Hulbert1993); see e.g. Isaksen et al. (Reference Isaksen, Bazilevs, Kvamsdal, Zhang, Kaspersen, Waterloo, Romner and Ingebrigtsen2008) for an application to haemodynamics.

Space discretization is typically based on finite elements. Whatever the implicit temporal scheme chosen, a system of non-linear algebraic equations is obtained after space and time discretization, that is,

where $\unicode[STIX]{x1D6FD}$ depends on the time discretization (e.g.  $\unicode[STIX]{x1D6FD}=4$ for the Newmark method (4.43)), $\boldsymbol{D}$ is the vector collecting the vessel wall displacement unknowns, $M_{s}$ is the mass matrix, $M_{s}^{\mathit{ext}}$ is the boundary mass matrix related to $\unicode[STIX]{x1D6F4}_{\mathit{ext}}$ , $\boldsymbol{\unicode[STIX]{x1D6E4}}$ is the non-linear operator defined by

where $\widehat{\boldsymbol{e}}_{i}$ is the $i$ th basis function, and $\boldsymbol{G}_{s}$ is the vector related to the right-hand side of the time-discretized equation. Notice that we have omitted the temporal index, which is understood. The above system is linearized by means of the Newton method, obtaining at each time step a sequence of linear systems of the form

where $k\geq 1$ is the Newton iteration index, to be solved until convergence occurs. Here $\unicode[STIX]{x1D6FF}\boldsymbol{D}_{(k)}=(\boldsymbol{D}_{(k)}-\boldsymbol{D}_{(k-1)})$ , $T$ is the matrix related to the linearization of the first Piola–Kirchhoff tensor, that is,

where $D_{F}$ is the Gâteaux derivative with respect to $\boldsymbol{F}$ .

For the solution of the above linear system, domain decomposition (DD) methods are often used as efficient preconditioners for iterative Krylov methods. Since the matrix $T$ is symmetric, the conjugate gradient method is usually considered for iterations. Among DD preconditioners, finite element tearing and interconnect (FETI) methods (Farhat and Roux Reference Farhat and Roux1991) are very often used in structural mechanics. In particular, all floating FETI methods have been considered for vessel wall problems, for example in Augustin, Holzapfel and Steinbach (Reference Augustin, Holzapfel and Steinbach2014). As in classical FETI methods, Lagrange multipliers are introduced to glue the solution together at the subdomain interfaces. In addition, Lagrange multipliers are also used to prescribe Dirichlet boundary conditions. This simplifies the implementation of the FETI method since all the subdomains are treated in the same way. A variant successfully used for arterial vessel walls is the so-called dual–primal FETI method: see e.g. Balzani et al. (Reference Balzani, Brands, Klawonn, Rheinbach and Schroder2010). Finally, we mention yet another class of DD methods considered for this problem, a two-level overlapping Schwarz method with an energy minimization coarse space: see Dohrmann and Widlund (Reference Dohrmann and Widlund2009).

In Figure 4.6 we illustrate some examples of numerical results obtained in real geometries reconstructed from radiological images (see the caption for details). These results highlight the anisotropic internal stresses characterizing vascular vessel walls.

Figure 4.6. (a) Von Mises internal stresses in a carotid artery (MRI images from the Vascular Surgery and Radiology Divisions at Ospedale Maggiore Policlinico, Milan, Italy). (b) Von Mises stresses in an abdominal aortic aneurysm (mesh from www.vascularmodel.com/sandbox/doku.php?id=start). Numerical results were obtained using LifeV (carotid artery) and the finite element library redbKIT v2.1 (github.com/redbKIT/redbKIT/releases) (AAA), P2 finite elements, a Newmark unconditionally stable scheme for the time discretization and an exponential vessel wall law.

4.6.3 Numerical methods for the fluid–structure interaction problem

The numerical solution of the coupled FSI problem (4.19) requires us to manage three sources of non-linearities, namely,

1. (i) the fluid domain is unknown (geometric non-linearity),

2. (ii) the fluid subproblem is non-linear (fluid constitutive non-linearity),

3. (iii) the vessel displacement subproblem is non-linear (structure constitutive non-linearity),

together with two different kinds of coupling,

1. (iv) the displacement of the fluid domain at the FS interface needs to match the displacement of the vessel wall (geometric adherence: see (4.19f )),

2. (iiv) the fluid and vessel displacement subproblems are coupled by means of the kinematic and dynamic conditions (4.19c , 4.19d ) (physical coupling).

Explicit scheme for geometric coupling and adherence

For $n\geq 1$ , at time step $t^{n}$ :

1. (1) Solve the FSI problem:

   (4.45a ) $$\begin{eqnarray}\displaystyle \displaystyle \frac{\unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D6E5}t}\boldsymbol{v}^{n}+\unicode[STIX]{x1D70C}_{f}((\boldsymbol{v}^{\ast }-\boldsymbol{v}_{f}^{\ast })\cdot \unicode[STIX]{x1D6FB})\boldsymbol{v}^{n}-\unicode[STIX]{x1D6FB}\cdot \boldsymbol{T}_{f}(\boldsymbol{v}^{n},p^{n})=\boldsymbol{g}_{f}^{n} & & \displaystyle \displaystyle \text{in }\unicode[STIX]{x1D6FA}_{f}^{\ast },\end{eqnarray}$$

   (4.45b ) $$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6FB}\cdot \boldsymbol{v}^{n}=0 & & \displaystyle \displaystyle \quad \text{in }\unicode[STIX]{x1D6FA}_{f}^{\ast },\end{eqnarray}$$

   (4.45c ) $$\begin{eqnarray}\displaystyle \displaystyle \boldsymbol{v}^{n}=\frac{\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D6E5}t}\boldsymbol{d}^{n}+\boldsymbol{g}_{fs}^{n} & & \displaystyle \displaystyle \text{on }\unicode[STIX]{x1D6F4}^{\ast },\end{eqnarray}$$

   (4.45d ) $$\begin{eqnarray}\displaystyle \displaystyle \boldsymbol{T}_{s}(\boldsymbol{d}^{n})\boldsymbol{n}^{\ast }=\boldsymbol{T}_{f}(\boldsymbol{v}^{n},p^{n})\boldsymbol{n}^{\ast } & & \displaystyle \displaystyle \text{on }\unicode[STIX]{x1D6F4}^{\ast },\end{eqnarray}$$

   (4.45e ) $$\begin{eqnarray}\displaystyle \displaystyle \frac{\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D6FD}}{\unicode[STIX]{x1D6E5}t^{2}}\widehat{\boldsymbol{d}}^{n}-\unicode[STIX]{x1D6FB}\cdot \widehat{\boldsymbol{T}}_{s}(\widehat{\boldsymbol{d}}^{n})=\widehat{\boldsymbol{g}}_{s}^{n} & & \displaystyle \displaystyle \text{in }\unicode[STIX]{x1D6FA}_{s},\end{eqnarray}$$
   where $\unicode[STIX]{x1D6FC}$ still depends on the time discretization scheme, and
   $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{g}_{f}^{n}=\boldsymbol{g}_{f}^{n}(\boldsymbol{v}^{n-1},\boldsymbol{v}^{n-2},\ldots ), & \displaystyle \nonumber\\ \displaystyle & \displaystyle \boldsymbol{g}_{s}^{n}=\boldsymbol{g}_{s}^{n}(\boldsymbol{d}^{n-1},\boldsymbol{d}^{n-2},\ldots ), & \displaystyle \nonumber\\ \displaystyle & \displaystyle \boldsymbol{g}_{fs}^{n} & \displaystyle =\boldsymbol{g}_{fs}^{n}(\boldsymbol{d}^{n-1},\boldsymbol{d}^{n-2},\ldots )\nonumber\end{eqnarray}$$
   account for the terms at previous time steps coming from time discretizations of order $p$ of the corresponding equations (4.44a ), (4.44e ) and (4.44c ), respectively.

2. (2) Then, solve the fluid geometry problem,

   (4.46a ) $$\begin{eqnarray}\displaystyle \displaystyle -\triangle \widehat{\boldsymbol{d}}_{f}^{n}=\mathbf{0} & & \displaystyle \displaystyle \text{in }\unicode[STIX]{x1D6FA}_{f},\end{eqnarray}$$

   (4.46b ) $$\begin{eqnarray}\displaystyle \displaystyle \widehat{\boldsymbol{d}}_{f}^{n}=\widehat{\boldsymbol{d}}^{n} & & \displaystyle \displaystyle \text{on }\unicode[STIX]{x1D6F4},\end{eqnarray}$$
   and build $\unicode[STIX]{x1D6FA}_{f}^{\ast }$ accordingly.

In the above substeps, the FSI problem (4.45) is still coupled by means of the physical coupling given by the interface conditions (4.45c , 4.45d ): see issue (v) above. For the solution of this problem, both partitioned and monolithic procedures have been successfully considered so far in haemodynamics. In partitioned schemes, the fluid and vessel wall subproblems are solved separately, one or more times per time step. Each of the two problems is equipped with a suitable boundary condition at the FS interface $\unicode[STIX]{x1D6F4}^{\ast }$ derived by splitting the physical interface conditions (4.45c , 4.45d ).

Explicit Dirichlet–Neumann scheme

For $n\geq 1$ , at time step $t^{n}$ :

1. (1) Solve the fluid Oseen problem with a Dirichlet condition at the FS interface:

   (4.47a ) $$\begin{eqnarray}\displaystyle \displaystyle \frac{\unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D6E5}t}\boldsymbol{v}^{n}+\unicode[STIX]{x1D70C}_{f}((\boldsymbol{v}^{\ast }-\boldsymbol{v}_{f}^{\ast })\cdot \unicode[STIX]{x1D6FB})\boldsymbol{v}^{n}-\unicode[STIX]{x1D6FB}\cdot \boldsymbol{T}_{f}(\boldsymbol{v}^{n},p^{n})=\boldsymbol{g}_{f}^{n} & & \displaystyle \displaystyle \text{in }\unicode[STIX]{x1D6FA}_{f}^{\ast },\end{eqnarray}$$

   (4.47b ) $$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6FB}\cdot \boldsymbol{v}^{n}=0 & & \displaystyle \displaystyle \text{in }\unicode[STIX]{x1D6FA}_{f}^{\ast },\end{eqnarray}$$

   (4.47c ) $$\begin{eqnarray}\displaystyle \displaystyle \boldsymbol{v}^{n}=\frac{\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D6E5}t}\boldsymbol{d}^{n-1}+\boldsymbol{g}_{fs}^{n-1} & & \displaystyle \displaystyle \text{on }\unicode[STIX]{x1D6F4}^{\ast }.\end{eqnarray}$$

2. (2) Then, solve the non-linear vessel wall problem with a Neumann condition at the FS interface:

   (4.48a ) $$\begin{eqnarray}\displaystyle \displaystyle \frac{\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D6FD}}{\unicode[STIX]{x1D6E5}t^{2}}\widehat{\boldsymbol{d}}^{n}-\unicode[STIX]{x1D6FB}\cdot \widehat{\boldsymbol{T}}_{s}(\widehat{\boldsymbol{d}}^{n})=\widehat{\boldsymbol{g}}_{s}^{n} & & \displaystyle \displaystyle \text{in }\unicode[STIX]{x1D6FA}_{s},\end{eqnarray}$$

   (4.48b ) $$\begin{eqnarray}\displaystyle \displaystyle \widehat{\boldsymbol{T}}_{s}(\widehat{\boldsymbol{d}}^{n})\widehat{\boldsymbol{n}}=\widehat{\boldsymbol{T}}_{f}(\widehat{\boldsymbol{v}}^{n},\widehat{p}^{n})\widehat{\boldsymbol{n}} & & \displaystyle \displaystyle \text{on }\unicode[STIX]{x1D6F4}.\end{eqnarray}$$

Note that the time-discretized kinematic condition (4.47c ) differs from (4.45c ) since now we are considering an explicit Dirichlet condition for the fluid subproblem, so the right-hand side is computed at the previous time step.

A ‘parallel’ version of this scheme is obtained by replacing $\widehat{\boldsymbol{T}}_{f}(\widehat{\boldsymbol{v}}^{n},\widehat{p}^{n})\widehat{\boldsymbol{n}}$ with $\widehat{\boldsymbol{T}}_{f}(\widehat{\boldsymbol{v}}^{n-1},\widehat{p}^{n-1})\widehat{\boldsymbol{n}}$ in (4.48b ). Notice that in the monolithic FSI problem (4.45), the dynamic continuity condition (4.45d ) is written in the current configuration $\unicode[STIX]{x1D6F4}^{\ast }$ , whereas for the structure subproblem alone (4.48) is written in the reference configuration $\unicode[STIX]{x1D6F4}$ . Accordingly, in what follows the structure interface quantities will be written in the current configuration in monolithic FSI problems and in the reference configuration when the structure problem is uncoupled in view of a partitioned scheme.

Unfortunately, the explicit Dirichlet–Neumann scheme can be unconditionally absolutely unstable. In particular, Causin, Gerbeau and Nobile (Reference Causin, Gerbeau and Nobile2005) prove that this happens if the fluid and structure densities are comparable, which is precisely the case for haemodynamics. This is due to the so-called high added mass effect. See also Forster, Wall and Ramm (Reference Forster, Wall and Ramm2007) and Nobile and Vergara (Reference Nobile and Vergara2012).

Stable loosely coupled algorithms have been introduced recently. To this end, in the FSI problem (4.45) the interface conditions (4.45c , 4.45d ) are replaced by two linear independent combinations

where $\unicode[STIX]{x1D70E}_{f}\neq \unicode[STIX]{x1D70E}_{s}$ are, in general, two functions of space and time. This naturally leads to the following.

Explicit Robin–Robin scheme

For $n\geq 1$ , at time step $t^{n}$ :

1. (1) Solve the Oseen problem with a Robin condition at the FS interface:

   (4.50a ) $$\begin{eqnarray}\displaystyle \displaystyle \frac{\unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D6E5}t}\boldsymbol{v}^{n}+\unicode[STIX]{x1D70C}_{f}((\boldsymbol{v}^{\ast }-\boldsymbol{v}_{f}^{\ast })\cdot \unicode[STIX]{x1D6FB})\boldsymbol{v}^{n} & & \displaystyle \nonumber\\ \displaystyle \displaystyle \qquad -\,\unicode[STIX]{x1D6FB}\cdot \boldsymbol{T}_{f}(\boldsymbol{v}^{n},p^{n})=\boldsymbol{g}_{f}^{n} & & \displaystyle \displaystyle \text{in }\unicode[STIX]{x1D6FA}_{f}^{\ast },\end{eqnarray}$$

   (4.50b ) $$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6FB}\cdot \boldsymbol{v}^{n}=0 & & \displaystyle \displaystyle \text{in }\unicode[STIX]{x1D6FA}_{f}^{\ast },\end{eqnarray}$$

   (4.50c ) $$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D70E}_{f}\boldsymbol{v}^{n}+\boldsymbol{T}_{f}(\boldsymbol{v}^{n},p^{n})\boldsymbol{n}^{\ast } & & \displaystyle \nonumber\\ \displaystyle \qquad \unicode[STIX]{x1D70E}_{f}\biggl(\frac{\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D6E5}t}\boldsymbol{d}^{n-1}+\boldsymbol{g}_{fs}^{n-1}\biggr)+\boldsymbol{T}_{s}(\boldsymbol{d}^{n-1})\boldsymbol{n}^{\ast } & & \displaystyle \displaystyle \text{on }\unicode[STIX]{x1D6F4}^{\ast }.\end{eqnarray}$$

2. (2) Then, solve the non-linear vessel wall problem with a Robin condition at the FS interface:

   (4.51a ) $$\begin{eqnarray}\displaystyle \displaystyle \frac{\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D6FD}}{\unicode[STIX]{x1D6E5}t^{2}}\widehat{\boldsymbol{d}}^{n}-\unicode[STIX]{x1D6FB}\cdot \widehat{\boldsymbol{T}}_{s}(\widehat{\boldsymbol{d}}^{n})=\widehat{\boldsymbol{g}}_{s}^{n} & & \displaystyle \displaystyle \text{in }\unicode[STIX]{x1D6FA}_{s},\end{eqnarray}$$

   (4.51b ) $$\begin{eqnarray}\displaystyle \displaystyle \frac{\unicode[STIX]{x1D70E}_{s}\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D6E5}t}\widehat{\boldsymbol{d}}^{n}+\widehat{\boldsymbol{T}}_{s}(\widehat{\boldsymbol{d}}^{n})\widehat{\boldsymbol{n}}=\unicode[STIX]{x1D70E}_{s}\widehat{\boldsymbol{v}}^{n}+\widehat{\boldsymbol{T}}_{f}(\widehat{\boldsymbol{v}}^{n},\widehat{p}^{n})\widehat{\boldsymbol{n}}-\unicode[STIX]{x1D70E}_{s}\widehat{\boldsymbol{g}}_{fs}^{n} & & \displaystyle \displaystyle \text{on }\unicode[STIX]{x1D6F4}.\end{eqnarray}$$

Burman and Fernández (Reference Burman and Fernández2009) propose a discontinuous Galerkin (DG)-like mortaring of the interface coupling conditions (4.45c , 4.45d ). The corresponding block Gauss–Seidel explicit algorithm could be reinterpreted as an explicit Dirichlet–Robin scheme, where in (4.50c , 4.51b ) we have $\unicode[STIX]{x1D70E}_{f}=+\infty$ , $\unicode[STIX]{x1D70E}_{s}=-\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D707}/h$ , where $\unicode[STIX]{x1D6FE}$ is the DG penalty parameter. Burman and Fernández show that the DG interface penalty and the viscous dissipation are not able to control the pressure fluctuations at the FS interface appearing in the discrete energy estimate. For this reason, they propose adding to the fluid problem a consistent stabilization term penalizing the pressure fluctuations, which is proved to be absolutely stable under a CFL-like condition. Banks, Henshaw and Schwendeman (Reference Banks, Henshaw and Schwendeman2014) have introduced a stable explicit Robin–Robin scheme, setting the parameters in the Robin interface conditions (4.50c , 4.51b ) after analysing the incoming and outgoing characteristic variables of the vessel wall problem. In particular, the above parameters are defined in terms of the outgoing characteristic variable for the fluid subproblem and in terms of the incoming characteristic variable for the vessel wall subproblem. For linear elasticity, the following values are obtained: $\unicode[STIX]{x1D70E}_{f,\text{norm}}=\sqrt{\unicode[STIX]{x1D70C}_{s}(\unicode[STIX]{x1D707}_{1}+2\unicode[STIX]{x1D707}_{2})}$ and $\unicode[STIX]{x1D70E}_{f,\text{tang}}=\sqrt{\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D707}_{2}}$ in the normal and tangential direction, respectively ( $\unicode[STIX]{x1D707}_{1}$ and $\unicode[STIX]{x1D707}_{2}$ are the Lamé constants), and $\unicode[STIX]{x1D70E}_{s}=-\unicode[STIX]{x1D70E}_{f}$ . This choice allows the travelling information contained in the characteristic variables to provide a tighter coupling of the fluid and structure problems than that enforced by (4.45c , 4.45d ).

Semi-implicit pressure-vessel wall coupled scheme

For $n\geq 1$ , at time step  $t^{n}$ :

1. (1) Solve the ALE-advection-diffusion problem with a Dirichlet condition at the FS interface:

   (4.52a ) $$\begin{eqnarray}\displaystyle \displaystyle \frac{\unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D6E5}t}\widetilde{\boldsymbol{v}}^{n}+\unicode[STIX]{x1D70C}_{f}((\boldsymbol{v}^{\ast }-\boldsymbol{v}_{f}^{\ast })\cdot \unicode[STIX]{x1D6FB})\widetilde{\boldsymbol{v}}^{n} & & \displaystyle \nonumber\\ \displaystyle \qquad -\,\unicode[STIX]{x1D6FB}\cdot \unicode[STIX]{x1D707}(\unicode[STIX]{x1D6FB}\widetilde{\boldsymbol{v}}^{n}+(\unicode[STIX]{x1D6FB}\widetilde{\boldsymbol{v}}^{n})^{T})=\boldsymbol{g}_{f}^{n} & & \displaystyle \displaystyle \text{in }\unicode[STIX]{x1D6FA}_{f}^{\ast },\end{eqnarray}$$

   (4.52b ) $$\begin{eqnarray}\displaystyle \displaystyle \widetilde{\boldsymbol{v}}^{n}=\frac{\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D6E5}t}\boldsymbol{d}^{n-1}+\boldsymbol{g}_{fs}^{n-1} & & \displaystyle \displaystyle \text{on }\unicode[STIX]{x1D6F4}^{\ast }.\end{eqnarray}$$

2. (2) Then, solve the coupled pressure-vessel wall problem. To this end, introduce the following iterations on index $k\geq 1$ .

   1. (a) Solve the pressure problem with a Neumann condition at the FS interface:

      $$\begin{eqnarray}\displaystyle \begin{array}{@{}ll@{}}\displaystyle \triangle p_{(k)}^{n}=\frac{\unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D6E5}t}\unicode[STIX]{x1D6FB}\cdot \widetilde{\boldsymbol{v}}_{(k)}^{n} & \displaystyle \quad \text{in }\unicode[STIX]{x1D6FA}_{f}^{\ast },\\ \displaystyle \frac{\unicode[STIX]{x2202}p_{(k)}^{n}}{\unicode[STIX]{x2202}\boldsymbol{n}^{\ast }}=\frac{\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D6E5}t}\boldsymbol{d}_{(k-1)}^{n}+\boldsymbol{g}_{fs}^{n} & \displaystyle \text{on }\unicode[STIX]{x1D6F4}^{\ast }.\end{array} & & \displaystyle \nonumber\end{eqnarray}$$

   2. (b) Then, solve the non-linear vessel wall problem with a Neumann condition at the FS interface:

      $$\begin{eqnarray}\displaystyle \begin{array}{@{}ll@{}}\displaystyle \frac{\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D6FD}}{\unicode[STIX]{x1D6E5}t^{2}}\widehat{\boldsymbol{d}}_{(k)}^{n}-\unicode[STIX]{x1D6FB}\cdot \widehat{\boldsymbol{T}}_{s}(\widehat{\boldsymbol{d}}_{(k)}^{n})=\widehat{\boldsymbol{g}}_{s}^{n} & \displaystyle \text{in }\unicode[STIX]{x1D6FA}_{s},\\ \displaystyle \widehat{\boldsymbol{T}}_{s}(\widehat{\boldsymbol{d}}_{(k)}^{n})\widehat{\boldsymbol{n}}=\widehat{\boldsymbol{T}}_{f}(\widehat{\widetilde{\boldsymbol{v}}}^{n},\widehat{p}_{(k)}^{n})\widehat{\boldsymbol{n}} & \displaystyle \text{on }\unicode[STIX]{x1D6F4}.\end{array} & & \displaystyle \nonumber\end{eqnarray}$$

Astorino, Chouly and Fernández (Reference Astorino, Chouly and Fernández2009) apply the DG mortaring approach to this projection scheme, leading to a Robin–Robin-like scheme. An algebraic version of the projection scheme proposed by Fernández et al. (Reference Fernández, Gerbeau and Grandmont2007) is introduced in Badia et al. (Reference Badia, Quaini and Quarteroni2008c ).

Implicit Robin–Robin scheme

For $n\geq 1,\,k\geq 1$ , at time step $t^{n}$ /iteration  $k$ :

1. (1) Solve the Oseen problem with a Robin condition at the FS interface:

   (4.55a ) $$\begin{eqnarray}\displaystyle \displaystyle \frac{\unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D6E5}t}\boldsymbol{v}_{(k)}^{n}+\unicode[STIX]{x1D70C}_{f}((\boldsymbol{v}^{\ast }-\boldsymbol{v}_{f}^{\ast })\cdot \unicode[STIX]{x1D6FB})\boldsymbol{v}_{(k)}^{n} & & \displaystyle \nonumber\\ \displaystyle \displaystyle \qquad -\,\unicode[STIX]{x1D6FB}\cdot \boldsymbol{T}_{f}(\boldsymbol{v}_{(k)}^{n},p_{(k)}^{n})=\boldsymbol{g}_{f}^{n} & & \displaystyle \displaystyle \text{in }\unicode[STIX]{x1D6FA}_{f}^{\ast },\end{eqnarray}$$

   (4.55b ) $$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6FB}\cdot \boldsymbol{v}_{(k)}^{n}=0 & & \displaystyle \displaystyle \text{in }\unicode[STIX]{x1D6FA}_{f}^{\ast },\end{eqnarray}$$

   (4.55c ) $$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D70E}_{f}\boldsymbol{v}_{(k)}^{n}+\boldsymbol{T}_{f}(\boldsymbol{v}_{(k)}^{n},p_{(k)}^{n})\boldsymbol{n}^{\ast } & & \displaystyle \nonumber\\ \displaystyle \displaystyle \qquad =\unicode[STIX]{x1D70E}_{f}\biggl(\frac{\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D6E5}t}\boldsymbol{d}_{(k-1)}^{n}+\boldsymbol{g}_{fs}^{n}\biggr)+\boldsymbol{T}_{s}(\boldsymbol{d}_{(k-1)}^{n})\boldsymbol{n}^{\ast } & & \displaystyle \displaystyle \text{on }\unicode[STIX]{x1D6F4}^{\ast }.\end{eqnarray}$$

2. (2) Then, solve the non-linear vessel wall problem with a Robin condition at the FS interface:

   (4.56a ) $$\begin{eqnarray}\displaystyle \displaystyle \frac{\unicode[STIX]{x1D70C}_{s}\unicode[STIX]{x1D6FD}}{\unicode[STIX]{x1D6E5}t^{2}}\widehat{\boldsymbol{d}}_{(k)}^{n}-\unicode[STIX]{x1D6FB}\cdot \widehat{\boldsymbol{T}}_{s}(\widehat{\boldsymbol{d}}_{(k)}^{n})=\widehat{\boldsymbol{g}}_{s}^{n} & & \displaystyle \displaystyle \text{in }\unicode[STIX]{x1D6FA}_{s},\end{eqnarray}$$

   (4.56b ) $$\begin{eqnarray}\displaystyle \displaystyle \frac{\unicode[STIX]{x1D70E}_{s}\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D6E5}t}\widehat{\boldsymbol{d}}_{(k)}^{n}+\widehat{\boldsymbol{T}}_{s}(\widehat{\boldsymbol{d}}_{(k)}^{n})\widehat{\boldsymbol{n}} & & \displaystyle \nonumber\\ \displaystyle \displaystyle \qquad =\unicode[STIX]{x1D70E}_{s}\widehat{\boldsymbol{v}}_{(k)}^{n}+\widehat{\boldsymbol{T}}_{f}(\widehat{\boldsymbol{v}}_{(k)}^{n},\widehat{p}_{(k)}^{n})\widehat{\boldsymbol{n}}-\unicode[STIX]{x1D70E}_{s}\widehat{\boldsymbol{g}}_{fs}^{n} & & \displaystyle \displaystyle \text{on }\unicode[STIX]{x1D6F4}.\end{eqnarray}$$

As proved in Causin et al. (Reference Causin, Gerbeau and Nobile2005), a small relaxation parameter is needed to achieve convergence in the implicit Dirichlet–Neumann scheme (corresponding to setting $\unicode[STIX]{x1D70E}_{f}=+\infty$ , $\unicode[STIX]{x1D70E}_{s}=0$ in (4.55c , 4.56b )). In practice, an Aitken relaxation procedure is often used to dynamically estimate an efficient relaxation parameter (Deparis Reference Deparis2004, Kuttler and Wall Reference Kuttler and Wall2008). A better situation is obtained by properly selecting the parameters in the Robin interface conditions (4.55c , 4.56b ). In particular, the choice

(where, as usual, $E$ and $\unicode[STIX]{x1D708}$ are Young’s modulus and the Poisson ratio for the vessel material at small deformations, and $H_{s}$ and $R$ are the representative thickness and radius of the vessel) yields fast convergence without any relaxation (Robin–Neumann scheme: Badia, Nobile and Vergara Reference Badia, Nobile and Vergara2008a , Nobile, Pozzoli and Vergara Reference Nobile, Pozzoli and Vergara2013). Gerardo-Giorda, Nobile and Vergara (Reference Gerardo-Giorda, Nobile and Vergara2010) and Gigante and Vergara (Reference Gigante and Vergara2015) have characterized the optimal value of $\unicode[STIX]{x1D70E}_{s}$ , leading to further improvement in the convergence history. Yu, Baek and Karniadakis (Reference Yu, Baek and Karniadakis2013) have derived a Dirichlet–Robin scheme by means of a generalized fictitious method, where the coefficients of the fluid pressure and vessel wall acceleration are changed to account for the added mass effect. This again allows one to obtain good convergence properties for haemodynamic parameters without any relaxation. Another class of implicit methods with good convergence properties for high added mass effect is based on adding a suitable interface artificial compressibility (IAC) consistent term to the fluid problem, proportional to the jump of pressure between two successive iterations (Degroote et al. Reference Degroote, Swillens, Bruggeman, Segers and Vierendeels2010). Degroote (Reference Degroote2011) showed that for a finite volume approximation, the IAC method based on Dirichlet–Neumann iterations is equivalent to a Robin–Neumann scheme for a suitable choice of the parameter $\unicode[STIX]{x1D70E}_{f}$ .

Explicit Dirichlet–Neumann scheme for the FSI problem with membrane structure

Given the quantities at previous time steps, at time step $t^{n}$ :

1. (1) Solve the Oseen problem with a Dirichlet condition at the FS interface:

   $$\begin{eqnarray}\displaystyle \begin{array}{@{}ll@{}}\displaystyle \frac{\unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D6E5}t}\boldsymbol{v}^{n}+\unicode[STIX]{x1D70C}_{f}((\boldsymbol{v}^{\ast }-\boldsymbol{v}_{f}^{\ast })\cdot \unicode[STIX]{x1D6FB})\boldsymbol{v}^{n}-\unicode[STIX]{x1D6FB}\cdot \boldsymbol{T}_{f}(\boldsymbol{v}^{n},p^{n})=\boldsymbol{g}_{f}^{n} & \displaystyle \text{in }\unicode[STIX]{x1D6FA}_{f}^{\ast },\\ \displaystyle \unicode[STIX]{x1D6FB}\cdot \boldsymbol{v}^{n}=0 & \displaystyle \text{in }\unicode[STIX]{x1D6FA}_{f}^{\ast },\\ \displaystyle \boldsymbol{v}^{n}\cdot \boldsymbol{n}^{\ast }=\frac{\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D6E5}t}d_{r}^{n-1}+g_{fs}^{n-1} & \displaystyle \text{on }\unicode[STIX]{x1D6F4}^{\ast },\\ \displaystyle \boldsymbol{v}^{n}-(\boldsymbol{v}^{n}\cdot \boldsymbol{n}^{\ast })\boldsymbol{n}^{\ast }=\mathbf{0} & \displaystyle \text{on }\unicode[STIX]{x1D6F4}^{\ast }.\end{array} & & \displaystyle \nonumber\end{eqnarray}$$

2. (2) Then, solve the membrane problem:

   $$\begin{eqnarray}\frac{\unicode[STIX]{x1D70C}_{s}H_{s}\unicode[STIX]{x1D6FD}}{\unicode[STIX]{x1D6E5}t^{2}}\widehat{d}_{r}^{n}-\unicode[STIX]{x1D6FB}\cdot (\boldsymbol{P}\unicode[STIX]{x1D6FB}\widehat{d}_{r}^{n})+\unicode[STIX]{x1D712}H_{s}\widehat{d}_{r}^{n}=-\widehat{\boldsymbol{T}}_{f}(\widehat{\boldsymbol{v}}^{n},\widehat{p}^{n})\widehat{\boldsymbol{n}}\cdot \widehat{\boldsymbol{n}}+\widehat{g}_{s}^{n}\quad \text{on }\unicode[STIX]{x1D6F4}.\end{eqnarray}$$

As in the case of scheme (4.47, 4.48), unfortunately this scheme is unconditionally absolutely unstable in the haemodynamic regime (Causin et al. Reference Causin, Gerbeau and Nobile2005).

Different algorithms are obtained by linearly combining the interface conditions (4.57c ) and (4.57e ) and by substituting the new condition in (4.57c ). In this case, we have to solve a coupled problem consisting of the following.

Robin–Neumann coupling for the FSI problem with membrane structure

Given the quantities at previous time steps, at time step $t^{n}$ :

1. (1) Solve the Oseen problem with a Robin condition at the FS interface:

   (4.59a ) $$\begin{eqnarray}\displaystyle \displaystyle \frac{\unicode[STIX]{x1D70C}_{f}\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D6E5}t}\boldsymbol{v}^{n}+\unicode[STIX]{x1D70C}_{f}((\boldsymbol{v}^{\ast }-\boldsymbol{v}_{f}^{\ast })\cdot \unicode[STIX]{x1D6FB})\boldsymbol{v}^{n}-\unicode[STIX]{x1D6FB}\cdot \boldsymbol{T}_{f}(\boldsymbol{v}^{n},p^{n})=\boldsymbol{g}_{f}^{n} & & \displaystyle \displaystyle \text{in }\unicode[STIX]{x1D6FA}_{f}^{\ast },\end{eqnarray}$$

   (4.59b ) $$\begin{eqnarray}\displaystyle \displaystyle \unicode[STIX]{x1D6FB}\cdot \boldsymbol{v}^{n}=0 & & \displaystyle \displaystyle \text{in }\unicode[STIX]{x1D6FA}_{f}^{\ast },\end{eqnarray}$$

   (4.59c ) $$\begin{eqnarray}\displaystyle \displaystyle (\unicode[STIX]{x1D70E}_{f}\boldsymbol{v}^{n}+\boldsymbol{T}_{f}(\boldsymbol{v}^{n},p^{n})\boldsymbol{n}^{\ast })\cdot \boldsymbol{n}^{\ast }\unicode[STIX]{x1D70E}_{f}\biggl(\frac{\unicode[STIX]{x1D6FC}}{\unicode[STIX]{x1D6E5}t}d_{r}^{n}+g_{fs}^{n}\biggr) & & \displaystyle \nonumber\\ \displaystyle \qquad -\,\biggl(\frac{\unicode[STIX]{x1D70C}_{s}H_{s}\unicode[STIX]{x1D6FD}}{\unicode[STIX]{x1D6E5}t^{2}}d_{r}^{n}-\unicode[STIX]{x1D6FB}\cdot (\boldsymbol{P}\unicode[STIX]{x1D6FB}d_{r}^{n})+\unicode[STIX]{x1D712}H_{s}d_{r}^{n}\biggr) & \displaystyle \text{on }\unicode[STIX]{x1D6F4}^{\ast }, & \displaystyle\end{eqnarray}$$

   (4.59d ) $$\begin{eqnarray}\displaystyle \displaystyle \boldsymbol{v}^{n}-(\boldsymbol{v}^{n}\cdot \boldsymbol{n}^{\ast })\boldsymbol{n}^{\ast }=\mathbf{0} & & \displaystyle \displaystyle \text{on }\unicode[STIX]{x1D6F4}^{\ast }.\!\end{eqnarray}$$

2. (2) Solve the membrane problem:

   (4.60) $$\begin{eqnarray}\frac{\unicode[STIX]{x1D70C}_{s}H_{s}\unicode[STIX]{x1D6FD}}{\unicode[STIX]{x1D6E5}t^{2}}\widehat{d}_{r}^{n}-\unicode[STIX]{x1D6FB}\cdot (\boldsymbol{P}\unicode[STIX]{x1D6FB}\widehat{d}_{r}^{n})+\unicode[STIX]{x1D712}H_{s}\widehat{d}_{r}^{n}=-\widehat{\boldsymbol{T}}_{f}(\widehat{\boldsymbol{v}}^{n},\widehat{p}^{n})\widehat{\boldsymbol{n}}\cdot \widehat{\boldsymbol{n}}+\widehat{g}_{s}^{n}\quad \text{on }\unicode[STIX]{x1D6F4}.\end{eqnarray}$$

The above problem could be solved either monolithically or by means of a block Gauss–Seidel method that in fact introduces sub-iterations splitting the solution of (4.59) and (4.60). When $\boldsymbol{P}=\mathbf{0}$ , the special choice

introduced in (4.59c ), yields

At this stage, this is a Robin condition for the fluid problem without any explicit dependence on $d_{r}^{n}$ . Thus, the monolithic problem given by (4.59, 4.60, 4.61) is equivalent to the stand-alone fluid problem (4.59a , 4.59b , 4.59d , 4.62): see Nobile and Vergara (Reference Nobile and Vergara2008). The solution of this fluid problem can then be used to feed the right-hand side of (4.60) and to get the structure displacement $\widehat{d}_{r}^{n}$ . In this way, the fluid and structure problems are in fact decoupled, even if the coupling conditions are treated implicitly. This provides a smart and efficient way to solve the monolithic problem (4.57) exactly, at the expense of a single fluid problem solve (note that the membrane problem (4.60) is solved very cheaply).

Starting from this result, Guidoboni et al. (Reference Guidoboni, Glowinski, Cavallini and Canic2009) have proposed a stable Robin–Neumann scheme based on an operator splitting, for a general membrane law ( $\boldsymbol{P}\neq \mathbf{0}$ ). The inertial vessel wall term is treated implicitly as in the previous case, leading to a Robin boundary condition for the fluid with $\unicode[STIX]{x1D70E}_{f}=\unicode[STIX]{x1D70C}_{s}H_{s}\unicode[STIX]{x1D6FD}/\unicode[STIX]{x1D6FC}\unicode[STIX]{x1D6E5}t$ , whereas the elastic and algebraic contributions are treated explicitly. Fernández (Reference Fernández2013) proposes an incremental version of this scheme, where the elastic and algebraic parts of the membrane law are included in the Robin condition for the fluid problem by means of a suitable extrapolation from previous time steps. Finally, we mention Colciago et al. (Reference Colciago, Deparis and Quarteroni2014), who treat the whole membrane law implicitly, leading to a generalized Robin condition, which however requires an ad hoc implementation.

4.6.4 Numerical methods for defective boundary problems

For the numerical solution of the augmented formulation (4.27), one could rely either on a monolithic strategy, where the full augmented matrix is constructed and solved, or on splitting techniques. For the latter, Formaggia et al. (Reference Formaggia, Gerbeau, Nobile and Quarteroni2002) and Veneziani and Vergara (Reference Veneziani and Vergara2005) propose writing the Schur complement equation with respect to the Lagrange multiplier of the linearized and discretized (in time and space) augmented formulation. This is a linear system whose dimension is equal to the number of flow rate conditions, say $m\geq 1$ . By using the GMRES method to solve this system iteratively, the exact solution is reached after exactly $m$ iterations (in exact arithmetic). At each iteration, the solution of a standard fluid problem with Neumann conditions is needed (exact splitting technique). The solution of a further standard fluid problem is required to compute the initial residual in the GMRES algorithm. This approach is quite expensive, even for the case $m=1$ , which requires the solution of two fluid problems per time step. However, it preserves modularity; indeed, it can be implemented using available standard fluid solvers in a black box mode. This is an interesting property when applications to cases of real interest are addressed: see Viscardi et al. (Reference Viscardi, Vergara, Antiga, Merelli, Veneziani, Puppini, Faggian, Mazzucco and Luciani2010), Vergara et al. (Reference Vergara, Viscardi, Antiga and Luciani2012), Piccinelli et al. (Reference Piccinelli, Vergara, Antiga, Forzenigo, Biondetti and Domanin2013) and Guerciotti et al. Reference Guerciotti, Vergara, Azzimonti, Forzenigo, Buora, Biondetti and Domanin2015, Reference Guerciotti, Vergara, Ippolito, Quarteroni, Antona and Scrofani2017.

To reduce the computational time required by the exact splitting approach, Veneziani and Vergara (Reference Veneziani and Vergara2007) propose a different (inexact) splitting procedure, requiring the solution of $m$ steady problems out of the temporal loop and of one unsteady null flow rate problem at each time step. This strategy introduces an error near the section which is smaller than the one based on conjecturing the velocity profile in the original (non-null) flow rate problem: see Section 4.4.1.

For the numerical solution of the control-based approach described in Section 4.4.3, Formaggia et al. Reference Formaggia, Veneziani and Vergara2008, Reference Formaggia, Veneziani and Vergara2009b have considered classical iterative methods for the solution of the resulting KKT system.

Recently, a numerical approach based on the Nitsche method has been considered to prescribe a flow rate condition. In particular, the original idea of prescribing Dirichlet conditions with a consistent penalization approach (Nitsche Reference Nitsche1970/71) has been extended by Zunino (Reference Zunino2009) to the case of flow rate boundary conditions. This strategy does not introduce further variables other than those of the original problem, but it does require careful tuning of a penalization parameter. In addition, it deals with non-standard bilinear forms that need ad hoc implementation. However, it should be very effective if flow rate conditions are implemented in a DG code. Vergara (Reference Vergara2011) has considered a similar approach to fulfilling the mean pressure condition (3.3) and the FSI case; see also Porpora et al. (Reference Porpora, Zunino, Vergara and Piccinelli2012).

For a more comprehensive overview of numerical strategies for defective boundary problems, we refer the reader to Formaggia and Vergara (Reference Formaggia and Vergara2012).

4.6.5 Numerical methods for the geometric reduced models and multiscale approach

For the numerical solution of the 1D reduced model (4.30), in principle any convenient approximation method for non-linear hyperbolic equations can be used. The peculiar feature of this model, however, is the lack of discontinuous solutions. A common approach relies on the finite element version of the Lax–Wendroff scheme, thanks to its excellent dispersion properties (Formaggia et al. Reference Formaggia, Gerbeau, Nobile and Quarteroni2001). As this scheme is explicit, a CFL-like condition is required to ensure absolute stability. In the presence of a viscoelastic term, the 1D model is usually discretized by means of a splitting procedure where the solution is split into two components, one satisfying the pure elastic problem and the other satisfying the viscoelastic correction (Formaggia et al. Reference Formaggia, Lamponi and Quarteroni2003, Malossi, Blanco and Deparis Reference Malossi, Blanco and Deparis2012). High-order methods are suitable for capturing the (physical) reflections at bifurcations induced by the vessel tapering: see e.g. Sherwin et al. (Reference Sherwin, Franke, Peiró and Parker2003b ) and Sherwin et al. (Reference Sherwin, Formaggia, Peiró and Franke2003a ) for a discontinuous Galerkin discretization and Muller and Toro (Reference Muller and Toro2014) for a finite volume scheme.

Regarding 0D models, they are in general described by systems of differential and algebraic equations, possibly non-linear due to the presence of diodes to represent the valves (Formaggia et al. Reference Formaggia, Quarteroni and Veneziani2009a ). Usually, for haemodynamic applications, these systems can be reduced to classical Cauchy problems and solved by classical Runge–Kutta methods.

Figure 4.8. Pressure wave propagation in an ascending aorta (3D model) and in the 1D model of the systemic circulation. Numerical results were obtained using LifeV; the Newton method was used for the interface equation. Courtesy of C. Malossi.

As for the solution of the 3D/1D coupled problems described in Section 4.5.2, we can in principle identify three different strategies, namely partitioned schemes, monolithic approaches, and methods based on the solution of an interface equation. In partitioned schemes, the 3D and 1D problems are solved separately in an iterative framework. The coupling interface conditions can be enforced in many different ways. For example, we can prescribe the flow rate condition (4.39a ) to the 3D problem and the pressure condition (4.42) to the 1D problem. Different algorithms are obtained by switching the role of the interface conditions in the iterative algorithm or by considering other interface conditions (e.g. (4.39b )). This is also the case when one of the two interface conditions is replaced by a condition expressing the continuity of the characteristic variable $W_{1}$ entering the 1D domain (Formaggia et al. Reference Formaggia, Gerbeau, Nobile and Quarteroni2001, Papadakis Reference Papadakis2009), that is, according to (4.33),

In any case, each of these approaches yields a 3D problem with a defective boundary condition, which could be tackled by one of the strategies described in Section 4.4. Formaggia et al. (Reference Formaggia, Moura and Nobile2007) and Papadakis (Reference Papadakis2009) have successfully considered explicit algorithms based on the solution of the 3D and 1D problems only once per time step. These algorithms enforce a limitation on $\unicode[STIX]{x1D6E5}t$ , which, however, is milder with respect to that imposed by the numerical scheme adopted for the 1D model. As an alternative, Blanco, Feijóo and Urquiza (Reference Blanco, Feijóo and Urquiza2007) and Blanco et al. (Reference Blanco, Pivello, Urquiza and Feijóo2009) have introduced iterative methods applied directly to the monolithic linearized system. A different approach to solving the 3D/1D coupled problem relies on writing an interface equation involving only the 3D/1D interface unknowns. We can interpret this equation as the geometric heterogeneous counterpart of the Schur complement equation. For its numerical solution, Leiva, Blanco and Buscaglia (Reference Leiva, Blanco and Buscaglia2011), Malossi et al. (Reference Malossi, Blanco, Crosetto, Deparis and Quarteroni2013) and Blanco, Deparis and Malossi (Reference Blanco, Deparis and Malossi2013) have used either the Broyden or the Newton method, in combination with GMRES. Methods relying on the numerical solution of the interface equation are simple to implement in the case of multiple interfaces, such as those arising in complex arterial networks.

In Figure 4.8 we give a numerical result obtained by the coupling between the 3D model of an ascending aorta and a 1D model of the systemic circulation. This result highlights the suitability of the 1D model in providing absorbing conditions for the 3D model and in propagating the pressure wave along the whole network (Malossi Reference Malossi2012).

7.1.1 The bidomain model

As observed in Section 5.3, the electrical activation of the heart is the result of two processes: at the microscopic scales, the generation of ionic currents through the cellular membrane producing a local action potential, and at the macroscopic scales, the travelling of the action potential from cell to cell allowed by the presence of the gap junctions. The former is a discrete process, in the sense that there is a delay between the depolarization of a cardiomyocyte and its neighbours, whereas the latter can be assimilated to a smooth process (Durrer et al. Reference Durrer, van Dam, Freud, Janse, Meijler and Arzbaecher1970).

Figure 7.1. Electrical circuit for the sequence of two cardiac cells. Each consists of a capacitor and a series of resistances, one for each ionic current (here only sodium and potassium ionic channels are depicted). In the intracellular region, two adjacent cells are connected by a resistance representing a gap junction. However, the latter is not explicitly modelled at the macroscopic scales: instead its effect is hidden in the conductivity tensor (see the text).

At the macroscopic level, the propagation of the potentials is described by means of partial differential equations, suitably coupled with ordinary differential equations modelling the ionic currents in the cells. In particular, the single membrane cell can be modelled as a capacitor separating charges that accumulate at its intracellular and extracellular surfaces. Moreover, as observed in Section 5, ionic currents cross the membrane through channels which open and close during excitation. A suitable model can therefore be expressed via the simple electric circuit depicted in Figure 7.1, for which

where $I_{m}$ is the membrane current per unit volume, $C_{m}$ is the membrane capacitance, $\unicode[STIX]{x1D712}_{m}$ is the surface area-to-volume ratio (see Section 6.3 for a quantification of the latter two), $V_{m}(t,\boldsymbol{x})$ is the trans-membrane potential, and $I_{\mathit{ion}}(t,\boldsymbol{x})$ are the ionic currents. Due to the conservation of current and charge, this current should equal the divergence of both the intracellular and extracellular current fluxes $\boldsymbol{j}_{i}$ and $\boldsymbol{j}_{e}$ :

Ohm’s law in the intracellular and extracellular regions gives

where $\boldsymbol{\unicode[STIX]{x1D6F4}}_{i},\boldsymbol{\unicode[STIX]{x1D6F4}}_{e}$ are the conductivity tensors and $\unicode[STIX]{x1D719}_{i}(t,\boldsymbol{x}),\unicode[STIX]{x1D719}_{e}(t,\boldsymbol{x})$ are the intracellular and extracellular potentials, so that

Note that, due to the anisotropy of the cardiac tissue induced by the presence of fibres and sheets, each conductivity tensor is in general expressed in terms of three scalar quantities representing the conductivities along the fibre direction $\boldsymbol{a}_{f}(\boldsymbol{x})$ , the direction $\boldsymbol{a}_{s}(\boldsymbol{x})$ orthogonal to $\boldsymbol{a}_{f}$ and tangential to sheets, and the direction $\boldsymbol{a}_{n}(\boldsymbol{x})$ orthogonal to sheets, that is,

Putting together all the above equations and using a homogenization procedure (see e.g. Colli Franzone et al. Reference Colli Franzone, Pavarino and Scacchi2014 for a rigorous derivation), we obtain for each $t>0$ the following system of two partial differential equations called the parabolic–parabolic (PP) formulation of the bidomain equations:

where $I_{i}^{\mathit{ext}}(t,\boldsymbol{x}),\,I_{e}^{\mathit{ext}}(t,\boldsymbol{x})$ are applied currents per unit volume.

Thanks to (7.1), $\unicode[STIX]{x1D6FB}\cdot (\boldsymbol{j}_{i}+\boldsymbol{j}_{e})=0$ ; thus, using (7.2) and (7.3), we obtain the following parabolic–elliptic (PE) formulation of the bidomain equations:

Due to the homogenization procedure, the effect of the gap junctions, which at the cellular level contributes to determining the current flux $\boldsymbol{j}_{i}$ , is hidden in the conductivity tensor $\boldsymbol{\unicode[STIX]{x1D6F4}}_{i}$ . We also notice that both bidomain problems (7.5) and (7.6) hold in the entire computational domain $\unicode[STIX]{x1D6FA}_{\mathit{mus}}$ given by the union of the myocardium with the endocardium and epicardium: see Figures 6.1 and 7.2. Indeed, again because of the homogenization procedure, no geometric distinction is made between the intracellular and extracellular regions, even if their different functionality is maintained in the bidomain models.

Figure 7.2. (a) Left ventricular myocardial domain obtained by the cut at the base (corresponding to $\unicode[STIX]{x1D6F4}_{b}$ ), and (b) corresponding fluid cavity domain.

7.1.2 Cardiac cell models

In order to close equations (7.5) and (7.6), we need to provide a model for the ionic current $I_{\mathit{ion}}$ . In what follows we briefly describe three families of models, featuring different levels of complexity and accuracy.

The first family, the so-called reduced ionic models, only provide a description of the action potential and disregard sub-cellular processes:

where $f,\boldsymbol{g}_{w}$ are suitable functions, while $\boldsymbol{w}:[0,T]\times \unicode[STIX]{x1D6FA}_{\mathit{mus}}\rightarrow \mathbb{R}^{M}$ collects the so-called gating variables which represent the percentage of open channels per unit area of the membrane. The most celebrated reduced model for ventricular cells is the FitzHugh–Nagumo model (FitzHugh Reference Fitzhugh1961), where $f(V_{m},w)=-kV_{m}(V_{m}-a)(V_{m}-1)-w$ and $g_{w}(V_{m},w)=\unicode[STIX]{x1D716}(V_{m}-\unicode[STIX]{x1D6FE}w)$ , for suitable constant parameters $k,\,a,\,\unicode[STIX]{x1D6FE}$ . In this case the gating variable $w$ in fact plays the role of a recovery function, which allows us to model the refractoriness of cells. Other, more sophisticated ventricular cell models of this family include the Rogers–McCulloch model (Rogers and McCulloch Reference Rogers and Mcculloch1994), the Aliev–Panfilov model (Aliev and Panfilov Reference Aliev and Panfilov1996), the Fenton–Karma model (Fenton and Karma Reference Fenton and Karma1998) and the Bueno-Orovio model (Bueno-Orovio, Cherry and Fenton Reference Bueno-Orovio, Cherry and Fenton2008). Like FitzHugh–Nagumo, the first two models are characterized by the dynamics of one gating variable and by a cubic non-linear expression of the ionic current. In contrast, the Fenton–Karma model and its Bueno-Orovio variant, specifically applied to human ventricular cells, include two and three gating variables, respectively, and a more complex non-linearity in the ionic current expression. These simple models are very appealing, especially because their parameters have a direct physical interpretation, such as the action potential duration, allowing easy setting of the model properties. For example, the Aliev–Panfilov model has been used successfully in the initial simulations of ventricular fibrillation in a real geometry (Panfilov Reference Panfilov1999). However, they are not able to describe any process occurring at the level of the ionic channels or the cell, so they are recommended when one is only interested in electrical activity of the heart.

The second family of ventricular cell models we consider is that of the so-called first-generation models. Unlike reduced models, which express the ionic current by means of the sole function $f$ , they allow explicit description of the kinetics of different ionic currents by using several gating variables. They are given by

where $N$ is the total number of ionic currents, $M$ is the total number of gating variables, $V_{k}$ is the Nernst potential of the $k$ th ion (a constant value corresponding to the thermodynamic equilibrium of the ion at hand), $I_{k}$ is the current related to the $k$ th ion, $p_{j_{k}}$ accounts for the influence of the $j$ th gating variable on the $k$ th ionic current (possibly vanishing) and $G_{k}$ is the maximal conductance of the $k$ th ion. The $M$ components of $\boldsymbol{g}$ usually have the expression

where $w_{j}^{\infty }$ is the equilibrium state and $\unicode[STIX]{x1D70F}_{j}$ the characteristic time constant. The most famous model of this family, the Hodgkin–Huxley (HH) model (Hodgkin and Huxley Reference Hodgkin and Huxley1952), depends on three ionic currents, namely the sodium, potassium and leakage currents, and three gating variables:

Although introduced to describe the action potentials in nerves, the HH model inspired all the following models introduced specifically for the ventricle. Among these, we cite the Beeler–Reuter model (Beeler and Reuter Reference Beeler and Reuter1977), the Luo–Rudy I model (Luo and Rudy Reference Luo and Rudy1991) and the ten Tusscher–Panfilov model (ten Tusscher and Panfilov Reference ten Tusscher and Panfilov2006). These models have been widely used to study specific features of the ventricular electrical activation, such as re-entry and fibrillation (Xie et al. Reference Xie, Qu, Yang, Baher, Weiss and Garfinkel2004).

Finally, we mention the family of second-generation ventricular cell models, such as the Luo–Rudy dynamic model (Luo and Rudy Reference Luo and Rudy1994a , Reference Luo and Rudy1994b ), which, unlike first-generation models, provide a detailed description of ion concentration variables $\boldsymbol{c}$ and of many processes allowing for the study of channelopathies and drug action, for example. However, due to their increased complexity, the computational time required is huge for a complete heart model, and tuning of parameters is often very demanding. We refer to Clayton and Panfilov (Reference Clayton and Panfilov2008) for a discussion of second-generation models and Rudy and Silva (Reference Rudy and Silva2006) for a general review of cardiac cell models.

Although most research has focused on ventricular cell models like those mentioned above, specific models have also been introduced for the atrial cells (e.g. Hilgemann and Noble Reference Hilgemann and Noble1987) and sinoatrial node cells (e.g. Yanagihara, Noma and Irisawa Reference Yanagihara, Noma and Irisawa1980).

All the cardiac cell models belonging to the three families described above, used in combination with the bidomain problem (7.5) or (7.6), lead to a system of two PDEs coupled with two systems of ODEs, that is, the equations for the gating variables and ion concentrations at each point $\boldsymbol{x}$ . The general form of this coupled problem is as follows (we only describe the PP case in detail). At each time $t>0$ , find the potentials $V_{M},\,\unicode[STIX]{x1D719}_{i}$ and $\unicode[STIX]{x1D719}_{e}$ , the gating variable $\boldsymbol{w}$ and the ion concentrations $\boldsymbol{c}$ such that

where, together with the notation introduced above for the reduced and first-generation models, $\boldsymbol{c}:[0,T]\times \unicode[STIX]{x1D6FA}_{\mathit{mus}}\rightarrow \mathbb{R}^{S}$ collects the $S$ ionic concentration variables and $\boldsymbol{g}_{c}$ is a suitable function: see e.g. Colli Franzone et al. (Reference Colli Franzone, Pavarino and Scacchi2014). We observe, in general, the dependence of the Nernst potential $V_{k}$ on the variable $\boldsymbol{c}$ . Well-posedness results of the previous coupled problem are provided in Colli Franzone and Savaré (Reference Colli Franzone, Savaré, Lorenzi and Ruf2002), where the existence and uniqueness of the solution of the PP formulation coupled with the FitzHugh–Nagumo model is proved, and in Bourgault, Coudière and Pierre (Reference Bourgault, Coudière and Pierre2006), where a Faedo–Galerkin technique is applied to the PE formulation coupled with a general first-generation cell model.

The ODE systems modelling the gating variables and the ionic concentration variables are in general stiff, since the Jacobians $\unicode[STIX]{x2202}\boldsymbol{g}_{w}/\unicode[STIX]{x2202}\boldsymbol{w}$ and $\unicode[STIX]{x2202}\boldsymbol{g}_{c}/\unicode[STIX]{x2202}\boldsymbol{c}$ exhibit a wide range of eigenvalues.

For each $t>0$ , the weak formulation of the bidomain model (7.9) together with homogeneous Neumann conditions and initial conditions (see Section 7.1.4) reads as follows. Given $I_{i}^{\mathit{ext}}(t),\,I_{e}^{\mathit{ext}}(t)\in L^{2}(\unicode[STIX]{x1D6FA}_{\mathit{mus}})$ , find $V_{m}(t),\unicode[STIX]{x1D719}_{e}(t),\unicode[STIX]{x1D719}_{i}(t)\in H^{1}(\unicode[STIX]{x1D6FA}_{\mathit{mus}})$ , $\boldsymbol{w}\in [L^{2}(\unicode[STIX]{x1D6FA}_{\mathit{mus}})]^{M}$ and $\boldsymbol{c}\in [L^{2}(\unicode[STIX]{x1D6FA}_{\mathit{mus}})]^{S}$ such that

for all $z\in H^{1}(\unicode[STIX]{x1D6FA}_{\mathit{mus}})$ , $\boldsymbol{y}\in [L^{2}(\unicode[STIX]{x1D6FA}_{\mathit{mus}})]^{M}$ and $\boldsymbol{\unicode[STIX]{x1D701}}\in [L^{2}(\unicode[STIX]{x1D6FA}_{\mathit{mus}})]^{S}$ , together with (7.9c , 7.9d ).

7.1.3 Reduced continuous models: the monodomain and eikonal equations

To reduce the complexity of the bidomain models (7.5) and (7.6), an assumption of proportionality between the intracellular and extracellular conductivities is introduced, that is, $\boldsymbol{\unicode[STIX]{x1D6F4}}_{e}=\unicode[STIX]{x1D709}\boldsymbol{\unicode[STIX]{x1D6F4}}_{i}$ for a suitable constant $\unicode[STIX]{x1D709}$ . Substituting this relation in (7.6b ), eliminating $\boldsymbol{\unicode[STIX]{x1D6F4}}_{e}$ and substituting the corresponding relation for $\boldsymbol{\unicode[STIX]{x1D6F4}}_{i}$ in (7.6a ), we obtain the following monodomain equation. For each $t>0$ , find the trans-membrane potential $V_{m}$ such that

where

is the effective conductivity and

Again, a ventricular cell model is needed to provide the ionic current $I_{\mathit{ion}}$ . The same models discussed above for the coupling with the bidomain problem are used in combination with the monodomain problem too. Once the trans-membrane potential $V_{m}$ has been computed, the extracellular potential $\unicode[STIX]{x1D719}_{e}$ could be computed as a postprocessing by solving the elliptic problem (7.6b ).

For each $t>0$ , the weak formulation of the monodomain problem (7.11) together with homogeneous Neumann conditions and initial conditions (see Section 7.1.4) reads as follows. Given $I^{\mathit{ext}}(t)\in L^{2}(\unicode[STIX]{x1D6FA}_{\mathit{mus}})$ , find $V_{m}(t)\in H^{1}(\unicode[STIX]{x1D6FA}_{\mathit{mus}})$ , $\boldsymbol{w}\in [L^{2}(\unicode[STIX]{x1D6FA}_{\mathit{mus}})]^{M}$ and $\boldsymbol{c}\in [L^{2}(\unicode[STIX]{x1D6FA}_{\mathit{mus}})]^{S}$ such that

for all $z\in H^{1}(\unicode[STIX]{x1D6FA}_{\mathit{mus}})$ , $\boldsymbol{y}\in [L^{2}(\unicode[STIX]{x1D6FA}_{\mathit{mus}})]^{M}$ and $\boldsymbol{\unicode[STIX]{x1D701}}\in [L^{2}(\unicode[STIX]{x1D6FA}_{\mathit{mus}})]^{S}$ , together with (7.9c , 7.9d ).

Although the hypothesis underlying the monodomain model, that is, the proportionality between the internal and external conductivities, is not physiological, as shown by experiment, in some cases this model provides a very accurate solution compared to the bidomain model. In particular, this is true when there is no injection of current in the extracellular region (Colli Franzone, Pavarino and Taccardi Reference Colli Franzone, Pavarino and Taccardi2005, Potse et al. Reference Potse, Dubé, Richer, Vinet and Gulrajani2006). In contrast, when an external current is injected such as in defibrillation, the monodomain solution is no longer accurate, and the bidomain model is mandatory since the unequal anisotropy is fundamental to successful description of these scenarios (Trayanova Reference Trayanova2006).

A further simplification is provided by the eikonal equation. Starting from the bidomain model coupled with a simplified representation of the ionic current, which does not consider any gating variable and allows for the description only of the depolarization phase, Colli Franzone, Guerri and Rovida (Reference Colli Franzone, Guerri and Rovida1990) derive the following eikonal-diffusion equation:

where $\unicode[STIX]{x1D713}(\boldsymbol{x})$ is the unknown activation time (see Section 6.2), $c_{o}$ denotes the velocity of the depolarization wave along the fibre direction for a planar wavefront, and $\boldsymbol{M}=\boldsymbol{\unicode[STIX]{x1D6F4}}/(\unicode[STIX]{x1D712}C_{m})$ . A different derivation has been provided in Keener (Reference Keener1991), leading to the following eikonal-curvature equation:

These are both steady equations providing information on the activation of each cell. The contours of $\unicode[STIX]{x1D713}(\boldsymbol{x})$ give the position of the wavefront at time $t=\unicode[STIX]{x1D713}$ . The eikonal-diffusive model (7.14) is an elliptic equation, where the propagation speed is influenced by the tissue surrounding the wavefront. Once the activation time $\unicode[STIX]{x1D713}$ has been computed, it is possible to obtain an approximate value of the extracellular potential $\unicode[STIX]{x1D719}_{e}$ by solving at each time step a suitable elliptic problem: see Colli Franzone and Guerri (Reference Colli Franzone and Guerri1993).

In contrast, the eikonal-curvature model (7.15) is of parabolic type since the ‘diffusive’ term lacks the second derivative in the direction of propagation. This term is also proportional to an anisotropic generalization of the mean curvature (Tomlinson, Hunter and Pullan Reference Tomlinson, Hunter and Pullan2002). This implies that the propagation is faster when the wavefront is concave. This is in accordance with the diffusion of charge, which allows for faster depolarization in regions close to already depolarized tissues.

The eikonal equations are unsuitable for recovering the action potential and the ionic currents. However, they provide accurate results about the activation of cells even in complex scenarios such as front–front collision: see e.g. Colli Franzone and Guerri (Reference Colli Franzone and Guerri1993) for the eikonal-diffusive model.

The eikonal models are, however, very appealing from the computational point of view. First of all, they consist of a single steady PDE. Although non-linear, they do not require coupling with ODE systems. More importantly, the activation time, unlike the trans-membrane potential, does not feature any internal or boundary layer, so no special restriction on the mesh is needed in this case (see Section 6.1).

7.1.4 Boundary conditions and Purkinje network models

We discuss here the initial and boundary conditions of the problems introduced above. The bidomain and monodomain equations and the ODE systems for the gating variables and ionic concentrations need to be equipped with suitable initial conditions, that is,

for given functions $V_{m,0}(\boldsymbol{x}),\boldsymbol{w}_{0}(\boldsymbol{x}),\,\boldsymbol{c}_{0}(\boldsymbol{x})$ .

For the boundary conditions for the bidomain, monodomain and eikonal-diffusion problems, a homogeneous Neumann condition is commonly prescribed at the external surface $\unicode[STIX]{x1D6F4}_{\mathit{epi}}$ of the epicardium and, for a ventricular domain only, at the base $\unicode[STIX]{x1D6F4}_{b}$ (see Figures 6.1 and 7.2(a)) to prescribe null outgoing current fluxes. In particular, the following conditions have to be prescribed on $\unicode[STIX]{x1D6F4}_{\mathit{epi}}\cup \unicode[STIX]{x1D6F4}_{b}$ :

Moreover, for the bidomain problems (7.5) and (7.6), they force the following compatibility conditions on the applied external currents:

On the internal surface $\unicode[STIX]{x1D6F4}_{\mathit{endo}}$ of the endocardium, again Neumann conditions are prescribed. In this case, however, they could be non-homogeneous at specific stimulation points (e.g. the atrioventricular node and the points of the bundle of His). For the eikonal problem, Dirichlet data on the activation time could be prescribed at some specific locations if they are available, thanks to the NavX system, for example; see Section 6.2 (Sermesant et al. Reference Sermesant, Chabiniok, Chinchapatnam, Mansi, Billet, Moireau, Peyrat, Wong, Relan, Rhode, Ginks, Lambiase, Delingette, Sorine, Rinaldi, Chapelle, Razavi and Ayache2012). When redundant (e.g. when the electrical problem in the myocardium is coupled with the Purkinje network: see below) these data have been used to solve inverse problems, for example to estimate the conduction velocity in the myocardium (Sermesant et al. Reference Sermesant, Chabiniok, Chinchapatnam, Mansi, Billet, Moireau, Peyrat, Wong, Relan, Rhode, Ginks, Lambiase, Delingette, Sorine, Rinaldi, Chapelle, Razavi and Ayache2012) or to obtain personalized Purkinje networks (Vergara et al. Reference Vergara, Palamara, Catanzariti, Nobile, Faggiano, Pangrazzi, Centonze, Maines, Quarteroni and Vergara2014, Palamara et al. Reference Palamara, Vergara, Catanzariti, Faggiano, Centonze, Pangrazzi, Maines and Quarteroni2014, Reference Palamara, Vergara, Faggiano and Nobile2015).

If the mathematical model accounts for the presence of the Purkinje network, interface conditions on $\unicode[STIX]{x1D6F4}_{\mathit{endo}}$ describing the continuity of the current and of the potential at the Purkinje muscle junctions (PMJs) are implicitly provided for the bidomain and monodomain problems by the solution of the coupled muscle region/Purkinje network problem (Vigmond and Clements Reference Vigmond and Clements2007, Vergara et al. Reference Vergara, Lange, Palamara, Lassila, Frangi and Quarteroni2016). For the sake of exposition, we will not detail the bidomain and monodomain models for the Purkinje network, instead referring interested readers to Vigmond and Clements (Reference Vigmond and Clements2007), Bordas et al. (Reference Bordas, Gillow, Gavaghan, Rodríguez and Kay2012) and Vergara et al. (Reference Vergara, Lange, Palamara, Lassila, Frangi and Quarteroni2016). We will only note that, unlike the muscular case, in the network the gap junctions connecting two consecutive Purkinje cells are often explicitly modelled by means of resistances. Specific Purkinje cell models have also been developed, with the same structure as those developed for the muscular cells: see e.g. DiFrancesco and Noble (Reference Difrancesco and Noble1985). However, we will describe the mechanisms of coupling; in particular, we will refer to the coupled problem obtained by considering the monodomain problem both in the muscle region and in the Purkinje network (Vergara et al. Reference Vergara, Lange, Palamara, Lassila, Frangi and Quarteroni2016). We consider $N$ Purkinje muscle junctions located at $\boldsymbol{x}=\boldsymbol{s}_{j}$ and we assume that each of them could be modelled by means of a resistance $R_{\mathit{PMJ}}$ . Then, the monodomain/monodomain coupled problem reads as follows. For each $t>0$ , find $V_{p},\,V_{m},\,\boldsymbol{w}_{p},\,\boldsymbol{w}$ and $\unicode[STIX]{x1D6FE}_{j},\,j=1,\ldots ,N$ , such that

where $P_{m}(V_{m},\boldsymbol{w},F)=0$ denotes the monodomain problem in the myocardium with source term $F$ , $P_{p}(V_{p},\boldsymbol{w}_{p},\boldsymbol{\unicode[STIX]{x1D702}})=0$ is the monodomain problem in the Purkinje network with Neumann conditions with data $\unicode[STIX]{x1D702}_{j}$ at the Purkinje muscle junctions, $V_{p}$ and $\boldsymbol{w}_{p}$ are the trans-membrane potential and the gating variables in the Purkinje network, $\unicode[STIX]{x1D6FE}_{j}$ are the PMJ currents which are determined by Ohm’s laws (7.16c ), ${\mathcal{I}}_{Y}$ is the characteristic function related to the region $Y\subset \unicode[STIX]{x1D6FA}_{\mathit{mus}}$ , ${\mathcal{B}}_{r}(\boldsymbol{s}_{j})$ is the ball of radius $r$ centred at the point $\boldsymbol{s}_{j}$ , and $A_{r}$ is the volume of this ball. We observe that the two monodomain problems are coupled by means of the PMJ currents $\unicode[STIX]{x1D6FE}_{j}$ : for the 3D problem the latter act as source terms distributed in balls of radius $r$ , whereas for the network they act as Neumann conditions (Bordas et al. Reference Bordas, Gillow, Gavaghan, Rodríguez and Kay2012). A similar approach could be considered for the bidomain problems as well.

The coupling between eikonal muscular and Purkinje network problems has been addressed by Vergara et al. (Reference Vergara, Palamara, Catanzariti, Nobile, Faggiano, Pangrazzi, Centonze, Maines, Quarteroni and Vergara2014) for normal propagation and Palamara et al. (Reference Palamara, Vergara, Catanzariti, Faggiano, Centonze, Pangrazzi, Maines and Quarteroni2014) for pathological propagations.