ARTICLE TYPE

# Left heart hemodynamics simulations with fluid-structure interaction and reduced valve modeling

Oscar Ruz<fnref n="1" /> | Jérôme Diaz<fnref n="2" /> | Marina Vidrascu<fnref n="1" /> | Philippe Moireau<fnref n="2" /> | Dominique Chapelle<fnref n="2" /> | Miguel A. Fernández<fnref n="1" />

**Correspondence**

Corresponding author M.A. Fernández,

Email: miguel.fernandez@inria.fr

**Abstract**

The combination of reduced models of cardiac valve dynamics with a one-way kinematic uncoupling of blood flow and electromechanics is a widespread approach for reducing the complexity of cardiac hemodynamics simulations. This comes however with a number of shortcomings: artificial pressure oscillations, missing isovolumetric phases and valve laws without precise continuous formulation. This paper is aimed at overcoming these three difficulties while still mitigating computational cost. A novel reduced model of valve dynamics is proposed in which unidirectional flow is enforced in a mathematically sound fashion. Artificial pressure oscillations are overcome by considering a fluid-structure interaction model, which couples bi-ventricular electromechanics and blood flow in the left cavities. The interface coupling is solved in a partitioned fashion via an unconditionally stable loosely coupled scheme. A priori energy estimates are derived for both the continuous coupled problem and its numerical approximation. The benefits and limitations of the proposed approaches are illustrated in a comprehensive numerical study.

**KEYWORDS**

Cardiac hemodynamics, fluid-structure interaction, reduced valve modeling, loosely coupled scheme.

1 | INTRODUCTION

The numerical simulation of cardiac hemodynamics (that is, of blood flows in the heart chambers) opens many fascinating perspectives for enhanced device design, treatment planning and clinical decision (see, e.g., <fnref n="1" />,<fnref n="2" /> for recent reviews), notably because computer simulations provide direct access to mechanical quantities (such as local pressures and shear stresses) which can hardly be measured in a non-invasive fashion. The genuine multi-physics nature of the heart function makes these simulations a permanent challenge in computational modeling and scientific computing. Indeed, a complete mathematical model of cardiac hemodynamics must involve also the interaction with different additional complex phenomena, such as cardiac electromechanics and valve dynamics.

The simulation of all the fluid-structure interaction phenomena involved in the heart is a complex problem, with only a few studies reported in the literature so far (see ³,⁴,⁵,⁶,⁷). Heart valves are definitely a bottleneck of the problem, particularly due to their large deflections, fast dynamics and contact phenomena at high pressure-drops. Fortunately, not all cardiac hemodynamics simulations require a complete multi-physics model of the heart, and alternative reduced modeling options can be considered depending on the relevant physics driving the phenomena of interest (see, e.g., <fnref n="1" />,<fnref n="2" />).

With the purpose of reducing the computational complexity of the problem, several reduced models of cardiac hemodynamics have been proposed in the literature. One of the most widespread approaches consists in combining a reduced modeling of the valve dynamics with a one-way kinematic uncoupling of blood flow and cardiac electromechanics. In short, motion data, coming either from electromechanical simulations (see, e.g., ⁸,⁹,¹⁰,¹¹,¹²,¹³) or from imaging measurements (see, e.g., ¹⁴,¹⁵,¹⁶,¹⁷,¹⁸,¹⁹,²⁰,²¹,²²), is enforced on the boundaries of the fluid cavities, while the opening and closing of the valves is described by switching boundary conditions (see, e.g., ²³,²⁴) or via resistive immersed methods (see ²⁵,²⁶,¹⁰,¹²).

<footnote n="1">¹Sorbonne Université, CNRS (LJLL), Inria, Paris, France</footnote>

<footnote n="2">²Inria & LMS, École Polytechnique, CNRS, Institut Polytechnique de Paris, France</footnote>
These modeling simplifications come however at a price. First, unphysical pressure oscillations are observed during the opening/closing phases of the valves. This often requires the prescription of artificial refractory periods. Second, the isovolumetric phases are neglected in the cardiac cycle, since the fluid pressure within the ventricle is not well defined in a fully enclosed domain with Dirichlet data. In¹⁰ (see also¹²), this issue has been overcome by leveraging the natural leakage induced by the resistive immersed valve model (for moderate values of the resistance), but it requires an a priori knowledge of the ventricle pressure. And third, the opening/closing valve laws generally lack mathematical soundness. They are indeed defined in a time semi-discrete framework (as a discrete algorithm rather than as a mathematical model), without a clear continuous counterpart. The contributions of this paper are aimed at overcoming these three issues.

We first introduce a novel reduced valve model in which unidirectional flow is enforced in a mathematically sound fashion through an inequality constraint on the fluid velocity across the immersed valve. This naturally avoids the need of artificially switching between different interface conditions or between different values of the resistance. Related unilateral constraints have been proposed in another context with the purpose of deriving energy estimates for the Navier-Stokes equations with outflow boundary conditions (see²⁷,²⁸,²⁹). In the context of kinematical uncoupling, numerical evidence indicates that this approach reduces the magnitude of the spurious pressure oscillations but it does not remove them.

We postulate that these artificial oscillations are due to the lack of compatibility between the kinematical data enforced on the cardiac cavities and the status of the reduced valve model. We investigate this hypothesis by proposing a cardiac hemodynamics model which combines the above mentioned reduced valve dynamics models with a fluid-structure interaction model in the left ventricle. Alternative fluid-structure interaction models for the left-ventricle have been proposed in³⁰,³¹,³²,³³,³⁴. A salient feature of our proposed approach is that it couples a bi-ventricular model of active cardiac mechanics with the Navier-Stokes system in the left heart. Furthermore, in order to mitigate computational complexity, notably with respect to the kinematically uncoupling approach, the interface coupling is numerically solved via a Robin-Robin loosely coupled scheme reported in³⁵ (see also³⁶,³⁷), which is expected to deliver nearly optimal accuracy. An a priori energy estimate, guaranteeing unconditional energy stability, is derived for the resulting fully discrete method. A comprehensive numerical study, involving a simplified and a realistic left-heart geometry, illustrates the accuracy and robustness of the proposed approach. Numerical evidence indicates, in particular, that fluid-structure interaction removes the above mentioned spurious pressures oscillations, irrespectively of the considered reduced valve model.

The rest of the paper is organized as follows. Section 2 presents the different mathematical models considered in this paper: from active cardiac mechanics, blood flow and reduced valve modeling to kinematical uncoupling and fluid-structure interaction. The numerical approximation of the proposed coupled fluid-structure interaction system is addressed in Section 3. We introduce a loosely coupled scheme and derive an energy stability estimate. Section 4 is devoted to the description and discussion of the numerical examples. Finally, a summary of the main results and some directions of future work are drawn in Section 5.

## 2 | PROBLEM SETTING AND MATHEMATICAL MODELS

In this section, we introduce a distributed fluid-structure interaction model of the left ventricle. The myocardium is described as a passive hyperelastic material with an active law along the muscle fibers (Section 2.2), whereas fluid flow within the left cavities is modeled by the incompressible Navier-Stokes equations with a reduced modeling of the valve dynamics (Section 2.3).

### 2.1 | Geometrical setting

We denote by $\Omega \subset \mathbb{R}^3$ the reference computational domain, partitioned as $\Omega_f \cup \Omega_s$, where $\Omega_f$ stands for the left heart fluid cavities and $\Omega_s$ represents the ventricles myocardium (see Figure 1). The fluid domain can be split into three regions $\Omega_f^{la} \cup \Omega_f^{ao} \cup \Omega_f^{lv}$, which respectively correspond to the fluid cavities in the left atria, ascending aorta and left ventricle. For the solid domain we consider the partition $\Omega_s^{lv} \cup \Omega_s^{rv}$, where $\Omega_s^{lv}$ and $\Omega_s^{rv}$ respectively denote the left and right ventricular myocardiums. As mentioned above, the present study focuses on the fluid-structure interaction between the left ventricular flow and the myocardium mechanics. This interaction takes place at the left ventricle endocardium interface, denoted by $\Gamma_1^{\text{def}} = \partial\Omega_s^{lv} \cap \partial\Omega_f^{lv}$.
<figure><img src="image_1.png" /><figcaption>FIGURE 1 Geometrical description of the reference computational domain $\Omega$.</figcaption></figure>

## 2.2 | Cardiac mechanics model

The mechanical response of the myocardium is described by the fundamental law of dynamics, formulated in a total Lagrangian approach, in terms of the deformation map $\phi: \Omega_s \times \mathbb{R}^+ \to \mathbb{R}^3$. This mapping describes the position of each point in the reference configuration at any given time instant. The solid displacement field $y: \Omega_s \times \mathbb{R}^+ \to \mathbb{R}^3$ is defined as $y(x, t) \stackrel{\text{def}}{=} \phi(x, t) - x$ for all $x \in \Omega_s$. We introduce the deformation gradient of the solid $F_s \stackrel{\text{def}}{=} \nabla\phi = I + \nabla y$, where $I$ denotes the identity matrix. The Jacobian of the solid system is given by the relation $J_s \stackrel{\text{def}}{=} \det F_s$. The boundary of $\Omega_s$ (see Figure 1) is assumed to be partitioned as $\Gamma_r \cup \Gamma_a \cup \Gamma_b \cup \Gamma_v \cup \Gamma_1$, where $\Gamma_r$ denotes the endocardium of the right ventricle, $\Gamma_a$ the apex of the heart, $\Gamma_b$ the base of the heart, $\Gamma_v$ the annulus of the mitral and tricuspid valves, and $\Gamma_1$ the endocardium of the left ventricle. For simplicity, the solid is assumed to be clamped on $\Gamma_v$, and a viscoelastic boundary condition is enforced on $\Gamma_a \cup \Gamma_b$, which emulates the support provided by the surrounding tissue (e.g., as the pericardium).

The mechanical response of the cardiac tissue is described in terms of the second Piola-Kirchhoff stress tensor $\Sigma$, which is assumed to be composed of a passive contribution $\Sigma_p$, and an active contribution $\Sigma_a$. We combine both contributions using an extension of the classical Hill-Maxwell rheology scheme, which takes into account the nonlinear response, as outlined in$^{38,39,40}$. Therefore, the overall second Piola-Kirchhoff stress tensor of the cardiac tissue can be expressed as $\Sigma = \Sigma_p + \Sigma_a$. The passive behavior, primarily influenced by the collagen and elastin matrix surrounding the heart fibers, is fully three-dimensional (3D) and is composed of hyperelastic and viscous contributions:

$$ \Sigma_p = \frac{\partial W_e(E)}{\partial E} + \frac{\partial W_v(\dot{E})}{\partial \dot{E}}, \qquad (1) $$

where $E$ stands for the Green-Lagrange strain tensor, defined by the relations

$$ E \stackrel{\text{def}}{=} \frac{1}{2}(C - I), \quad C \stackrel{\text{def}}{=} F_s^T F_s. $$
The cardiac tissue is modeled as nearly incompressible and transversely isotropic, with a privileged direction along the muscle fibers, represented by their unit tangential vector **a**. The hyperelastic potential is inspired from the Holzapfel–Ogden constitutive law (see$^{41}$) and defined as

$$W_e = \mu_1(J_1 - 3) + \mu_2(J_2 - 3) + C_0 e^{C_1(J_1-3)^2} + C_2 e^{C_3(J_4-1)^2} + \kappa(J_s - 1) - \kappa \ln(J_s),$$

where $J_1 \stackrel{\text{def}}{=} \text{trC}$, $J_2 \stackrel{\text{def}}{=} \frac{1}{2}((\text{trC})^2 - \text{trC}^2)$, $J_4 \stackrel{\text{def}}{=} a^T C a$ are the reduced invariants and $\mu_1, \mu_2, C_0, C_1, C_2, C_3, \kappa$ denote model parameters adjusted on the specific regions of the heart tissue. On the other hand, we consider the following simple relation for the viscous pseudo-potential

$$W_v = \frac{\nu}{2} I : (\dot{E})^2, \qquad (2)$$

where $\nu$ denotes the passive viscosity coefficient.

The active behavior, driven by the chemical reactions within myocytes that induce the contraction and relaxation of sarcomeres (resulting in the subsequent contraction and relaxation of the myocardium along the fibers of the heart), is modeled following the approach outlined in$^{38,39}$. The dynamics of the active stiffness $k_c$ and of the contractile stress $\tau_c$ produced by the sarcomeres are hence given by

$$\begin{cases} \dot{k}_c = -(|u_e| + \alpha_d |\dot{e}_c|) k_c + n_0 k_0 |u_e|, \\ \dot{\tau}_c = -(|u_e| + \alpha_d |\dot{e}_c|) \tau_c + \dot{e}_c k_c + n_0 \sigma_0 |u_e|, \end{cases} \qquad (3)$$

where $e_c$ represents the extension along the fibers, which is further elaborated below, and $u_e$ denotes the chemical activation (or electrical input) and $|x|_+ \stackrel{\text{def}}{=} \max\{x, 0\}$. The term $\alpha_d |\dot{e}_c|$ accounts for the bridges destruction upon rapid length changes. The parameter $0 \le n_0 \le 1$ is a reduction factor that accounts for the Frank-Starling effect, which describes the relation of the maximum active tension with the current strain of the sarcomeres. The constant $k_0$ is the maximum stiffness, while $\sigma_0$ accounts for the solid contractibility.

<figure><img src="image_2.png" /><figcaption>FIGURE 2 Rheological model from$^{38,40}$.</figcaption></figure>

Following$^{38,40}$, the passive and active elements of the myocardium tissue are assembled using the rheological scheme of Figure 2. The total stress acting along the fiber results from a balance of contractile and elastic stresses, respectively denoted $\sigma_c$ and $\sigma_s$. The contractile stress is given by the constitutive relation

$$\sigma_c = \tau_c + \mu_c \dot{e}_c, \qquad (4)$$
where $\mu_c > 0$ is a dissipation coefficient. On the other hand, the elastic stress produced in the sarcomeres (due to the elastic response of the Z disks) can be modeled as

$$ \sigma_s = E_s e_s, \qquad (5) $$

where $E_s$ represents the Young modulus of the elastic component of the sarcomeres and $e_s$ is the elastic extension of these components. The overall stress along the fibers, denoted by $\sigma_a$, can be related with the contractile and elastic stresses according to the following relations (see$^{38}$):

$$ \sigma_a = \frac{\sigma_c}{1 + 2e_s} = \frac{\sigma_s}{1 + 2e_c}. \qquad (6) $$

On the other hand, the extension along the fibers, denoted by $e_a$, is related with the contractile and elastic extensions as

$$ 1 + 2e_a = (1 + 2e_s)(1 + 2e_c). \qquad (7) $$

Additionally, $e_a$ is given in terms of the Green-Lagrange strain tensor as

$$ e_a \stackrel{\text{def}}{=} a^T E a. \qquad (8) $$

The final relation for the active behaviour is described as

$$ \Sigma_a = \sigma_a a \otimes a. $$

Note that from, (7), we get

$$ e_s = \frac{e_a - e_c}{1 + 2e_c}, \qquad (9) $$

so that this expression can be used to eliminate $e_s$ in the above relations. In particular, owing to (5) and (6), we have

$$ \sigma_a = E_s \frac{e_a - e_c}{(1 + 2e_c)^2}. \qquad (10) $$

Similarly, by combining (4), (6) and (9), we get the following dynamics for $e_c$

$$ \mu_c \dot{e}_c + \tau_c = E_s \frac{(e_a - e_c)(1 + 2e_a)}{(1 + 2e_c)^3}. $$

We can then use (8) to express $e_a$ in terms of $E$ in the two previous relations, which together with (3) closes the active modeling.

A lumped parameter (0D) model is used to describe the blood flow in the right cardiac cavities (see, e.g.,$^{40}$). The coupling with ventricular tissue is operated by the relations

$$ \left\{ \begin{array}{l} F_s \Sigma n_s = -P_r n_s \quad \text{on } \Gamma_r, \\ \displaystyle \int_{\Gamma_r} J_s u_s \cdot F_s^{-T} n_s = Q_r (P_r, P_{ar,r}, P_{at,r}), \end{array} \right. $$

The first enforces the dynamic balance on the fluid-solid interface $\Gamma_r$, in terms of the unknown homogeneous fluid pressure $P_r$, while the second corresponds to the balance of mass in the right ventricular cavity, with $Q_r$ denoting the net flux across the valves given by the function

$$ Q_r(P_r, P_{ar,r}, P_{at,r}) \stackrel{\text{def}}{=} \begin{cases} K_{at,r}(P_r - P_{at,r}) & \text{if } P_r \le P_{at,r}, \\ 0 & \text{if } P_{at,r} \le P_r \le P_{ar,r}, \\ K_{ar,r}(P_r - P_{ar,r}) & \text{if } P_r \ge P_{ar,r}. \end{cases} \qquad (11) $$

Here, $P_{ar,r}$, $P_{at,r}$ respectively denote the right atrium and pulmonary artery pressures, respectively, while $K_{at,r}$, $K_{ar,r}$ stand for the conductances of the valves in the open configurations. Furthermore, a two-stage Windkessel system is used to model the impact of the external circulation on the dynamics of the unknown arterial pressure $P_{ar,r}$ (see, e.g.,$^{42}$), viz.,

$$ \left\{ \begin{array}{l} C_{p,r} \dot{P}_{ar,r} + \frac{P_{ar,r} - P_{d,r}}{R_{p,r}} = |Q_r|_+, \\ C_{d,r} \dot{P}_{d,r} + \frac{P_{d,r} - P_{ar,r}}{R_{p,r}} = \frac{P_{vs,r} - P_{d,r}}{R_{d,r}}, \end{array} \right. \qquad (12) $$
where $P_{d,r}$ stands for the unknown distal pressure, $P_{vs,r}$ for a given remote venous pressure and $R_{p,r}$, $C_{p,r}$ and $R_{d,r}$, $C_{d,r}$ respectively denote the proximal and distal parameters.

The resulting active cardiac mechanics system reads therefore as follows: find the displacement $\mathbf{y} : \Omega_s \times \mathbb{R}^+ \to \mathbb{R}^3$, the velocity $\mathbf{u}_s : \Omega_s \times \mathbb{R}^+ \to \mathbb{R}^3$, with $\mathbf{u}_s = \partial_t \mathbf{y}$, the extension along the fibers $e_c : \Omega_s \times \mathbb{R}^+ \to \mathbb{R}$, the active stiffness $k_c : \Omega_s \times \mathbb{R}^+ \to \mathbb{R}^+$, the active tension $\tau_c : \Omega_s \times \mathbb{R}^+ \to \mathbb{R}$, the right ventricular pressure $P_r : \mathbb{R}^+ \to \mathbb{R}$, the right proximal pressure $P_{ar,r} : \mathbb{R}^+ \to \mathbb{R}$ and the right distal pressure $P_{d,r} : \mathbb{R}^+ \to \mathbb{R}$ such that

$$
\left\{
\begin{array}{ll}
\rho_s \partial_t \mathbf{u}_s - \mathbf{div}(\mathbf{F}_s \boldsymbol{\Sigma}(\mathbf{E}, e_c)) = \mathbf{0} & \text{in } \Omega_s, \\
\boldsymbol{\Sigma} = \frac{\partial W_e}{\partial \mathbf{E}} + \frac{\partial W_v}{\partial \dot{\mathbf{E}}} + \frac{E_s(\mathbf{a}^\mathrm{T} \mathbf{E} \mathbf{a} - e_c)}{(1+2e_c)^2} \mathbf{a} \otimes \mathbf{a} & \text{in } \Omega_s, \\
\mu_c \dot{e}_c + \tau_c = E_s \frac{(\mathbf{a}^\mathrm{T} \mathbf{E} \mathbf{a} - e_c)(1+2\mathbf{a}^\mathrm{T} \mathbf{E} \mathbf{a})}{(1+2e_c)^3} & \text{in } \Omega_s, \\
\dot{k}_c = -(|u_e| + \alpha_d |\dot{e}_c|) k_c + n_0 k_0 |u_e|_{l_4} & \text{in } \Omega_s, \\
\dot{\tau}_c = -(|u_e| + \alpha_d |\dot{e}_c|) \tau_c + \dot{e}_c k_c + n_0 \sigma_0 |u_e|_{l_4} & \text{in } \Omega_s, \\
\mathbf{y} = \mathbf{0} & \text{on } \Gamma_v, \\
\mathbf{F}_s \boldsymbol{\Sigma} \mathbf{n}_s + \mathbf{ay} + \mathbf{bu}_s = \mathbf{0} & \text{on } \Gamma_a \cup \Gamma_b, \\
\mathbf{F}_s \boldsymbol{\Sigma} \mathbf{n}_s = -P_r \mathbf{n}_s & \text{on } \Gamma_r, \\
\int_{\Gamma_r} J_s \mathbf{u}_s \cdot \mathbf{F}_s^{-\mathrm{T}} \mathbf{n}_s = Q_r(P_r, P_{ar,r}, P_{at,r}), \\
C_{p,r} \dot{P}_{ar,r} + \frac{P_{ar,r} - P_{d,r}}{R_{p,r}} = |Q_r|_{l_4}, \\
C_{d,r} \dot{P}_{d,r} + \frac{P_{d,r} - P_{ar,r}}{R_{p,r}} = \frac{P_{vs,r} - P_{d,r}}{R_{d,r}},
\end{array}
\right.
\quad (13)
$$

where $\rho_s$ stands for the solid density and $a, b > 0$ are given coefficients of the viscoelastic support boundary condition on $\Gamma_a \cup \Gamma_b$. The type of boundary conditions imposed on the left endocardium wall $\Gamma_1$ depends on the type of model coupling approach considered with the fluid in the left ventricle. This will be discussed in Section 2.4 below.

## 2.3 | Fluid model with reduced valve dynamics

The dynamics of the fluid flow within the left cavities are described by the incompressible Navier-Stokes equations in arbitrary Lagrangian-Eulerian (ALE) formulation. The current configuration of the fluid domain, denoted by $\Omega_f(t)$, is hence parametrized by the so-called ALE map $\mathcal{A} : \Omega_f \times \mathbb{R}^+ \to \mathbb{R}^3$, with $\mathcal{A}_t \stackrel{\text{def}}{=} \mathcal{A}(\cdot, t)$ one-to-one and so that $\Omega_f(t) = \mathcal{A}(\Omega_f, t)$, for all $t > 0$. The ALE map $\mathcal{A}$ can be given in terms of the displacement of the fluid domain $\mathbf{y}_f : \Omega_f \times \mathbb{R}^+ \to \mathbb{R}^3$, as $\mathcal{A}(\mathbf{x}, t) \stackrel{\text{def}}{=} \mathbf{x} + \mathbf{y}_f(\mathbf{x}, t)$. Note that both $\mathcal{A}$ and $\mathbf{y}_f$ are arbitrarily defined in the interior of $\Omega_f$. With the ALE map, we can define the computational domain velocity, also called ALE velocity, as $\mathbf{w} = \partial_t \mathcal{A}$, the gradient of deformation $\mathbf{F} \stackrel{\text{def}}{=} \nabla \mathcal{A}$ and the Jacobian $J \stackrel{\text{def}}{=} \det \mathbf{F}$.

The boundary of the fluid domain $\Omega_f$ is assumed to be partitioned as $\Gamma_{la} \cup \Gamma_{aw} \cup \Gamma_{pv} \cup \Gamma_{ao} \cup \Gamma_1$ (see Figure 1), where $\Gamma_{la} \subset \partial\Omega_f^{la}$ denotes the left atrium wall, $\Gamma_{aw} \subset \partial\Omega_f^{ao}$ the wall of the ascending aorta, $\Gamma_{pv} \subset \partial\Omega_f^{la} \setminus \Gamma_{la}$ the inlet of the pulmonary veins left atrium and $\Gamma_{ao} \subset \partial\Omega_f^{la} \setminus \Gamma_{aw}$ the outlet of the ascending aorta. Furthermore, the closed configuration of the mitral and aortic valves are represented through two immersed surfaces, denoted by $\Gamma_{mv}$ and $\Gamma_{av}$, respectively, oriented with the normal vector $\mathbf{n}_f$ along the flow direction.

**Remark 1.** For a given field $f$ defined in the current configuration $\Omega_f(t)$, we shall use the notation $\hat{f}$ to defined its corresponding ALE description in $\Omega_f$, by composing with the ALE map, namely, $\hat{f}(\mathbf{x}, t) \stackrel{\text{def}}{=} f(\mathcal{A}(\mathbf{x}, t), t) \quad \forall \mathbf{x} \in \Omega_f$. Conversely, a given field defined in the reference configuration $\Omega_f$ can also be transported into the current configuration by composition with $\mathcal{A}_t^{-1}$. For the sake of simplicity, this composition is omitted. For instance, $\mathbf{w}$ in (14)$_1$ below has to be read as $\mathbf{w} \circ \mathcal{A}_t^{-1}$.
The fluid problem in the left cavities with a reduced model of the valve dynamics consists in finding the fluid velocity $\hat{u} : \Omega_f \times \mathbb{R}^+ \to \mathbb{R}^3$, and the fluid pressure $\hat{p} : \Omega_f \times \mathbb{R}^+ \to \mathbb{R}$, such that

$$ 
\left\{ 
\begin{array}{ll} 
\rho_f (\partial_t \mathbf{u}|_{\mathcal{A}} + (\mathbf{u}-\mathbf{w}) \cdot \nabla \mathbf{u}) - \mathbf{div} \sigma(\mathbf{u}, p) + \sum_{i \in \{\text{mv,av}\}} f_{\Gamma_i}(\mathbf{u}-\mathbf{w}, p) \delta_{\Gamma_i} = \mathbf{0} & \text{in } \Omega_f(t), \\ 
\mathbf{div} \mathbf{u} = 0 & \text{in } \Omega_f(t), \\ 
\mathbf{u} = \mathbf{0} & \text{on } \Gamma_{\text{la}} \cup \Gamma_{\text{aw}}, \\ 
\sigma(\mathbf{u}, p) \mathbf{n}_f + \frac{\rho_f}{2} |\mathbf{u} \cdot \mathbf{n}_f|_- \mathbf{u} = -P_{\text{pv,ao}} \mathbf{n}_f & \text{on } \Gamma_{\text{pv}} \cup \Gamma_{\text{ao}}, 
\end{array} 
\right. 
\quad (14) 
$$

with $\rho_f$ standing for the fluid density, $\delta_{\Gamma_i}$ for the Dirac measure on the immersed surface $\Gamma_i$ and $\sigma(\mathbf{u}, p)$ for the fluid Cauchy stress tensor, given by:

$$ \sigma(\mathbf{u}, p) \stackrel{\text{def}}{=} -p\mathbf{I} + 2\mu_f \varepsilon(\mathbf{u}), \quad \varepsilon(\mathbf{u}) \stackrel{\text{def}}{=} \frac{1}{2} (\nabla \mathbf{u} + \nabla \mathbf{u}^T), $$

where $\mu_f > 0$ denotes the dynamic viscosity of the fluid. Furthermore we have used the general notation $f \stackrel{\text{def}}{=} \hat{f} \circ \mathcal{A}_t^{-1}$, for any field $\hat{f}$ defined in the reference configuration $\Omega_f$.

As boundary conditions, we impose a non-slip condition on $\Gamma_{\text{la}} \cup \Gamma_{\text{aw}}$, while we consider normal tractions on $\Gamma_{\text{pv}} \cup \Gamma_{\text{ao}}$. The enforced pressure $P_{\text{pv,ao}}(t)$ is analytically given in $\Gamma_{\text{pv}}$, while in $\Gamma_{\text{ao}}$ it is determined through a one-stage Windkessel model. Furthermore, in order to circumvent potential backflow numerical instabilities on Neumann boundaries (see, e.g.,$^{43}$), we include the following dissipative term, $\partial_t \frac{\rho_f}{2} |\mathbf{u} \cdot \mathbf{n}_f|_- \mathbf{u}$, with $|x| \stackrel{\text{def}}{=} \max\{0, -x\}$. Finally, the term $f_{\Gamma_i}(\mathbf{u}-\mathbf{w}, p)$ represents the fluid stress jump $\llbracket \sigma(\mathbf{u}, p) \mathbf{n}_f \rrbracket$ induced by the opening and closing dynamics of the valve $\Gamma_i$, with the notation

$$ \llbracket f \rrbracket \stackrel{\text{def}}{=} f^+ - f^-, \quad \llbracket f \mathbf{n} \rrbracket \stackrel{\text{def}}{=} \llbracket f \rrbracket \mathbf{n}, \quad f^\pm(\mathbf{x}) \stackrel{\text{def}}{=} \lim_{\epsilon \to 0^\pm} f(\mathbf{x} + \epsilon \mathbf{n}_f) \quad \forall \mathbf{x} \in \Gamma_i. $$

The term $f_{\Gamma_i}(\mathbf{u}-\mathbf{w}, p)$ can be modeled in different ways. In this work, we consider the resistive immersed surface (RIS) approach and we propose as well a novel reduced valve modeling approach, with superior mathematical foundations.

### 2.3.1 | Resistive immersed surface valve model

In the RIS model, introduced in$^{25}$ (see also$^{26,10}$), the expression of $f_{\Gamma_i}(\mathbf{u}-\mathbf{w}, p)$ is given by a distributed Ohm's law, namely,

$$ f_{\Gamma_i}(\mathbf{u}-\mathbf{w}, p) = R_{\Gamma_i}(t)(\mathbf{u}-\mathbf{w}) \quad \text{on } \Gamma_i, \quad (15) $$

where $R_{\Gamma_i}(t)$ plays the role of a time-dependent resistance which regulates the fluid flow across the valve $\Gamma_i$. Typically, $R_{\Gamma_i}(t)$ vanishes when the valve is open and takes a large value $R_{\text{closed}}$ whenever the valve is closed. The transition between these two states of the valve depends on both the current status of the valve (closed or open) and on the nature of the flow across $\Gamma_i$. Roughly, the basic idea consists in closing the valve when the relative flow rate becomes negative (or smaller than a given tolerance) and opening the valve when the pressure jump across the valve reaches a positive value. The dependence of $R_{\Gamma_i}$ on $\mathbf{u}-\mathbf{w}$ and $p$ is hence rather intricate. In fact, instead of providing an explicit continuous expression of $R_i$ in terms of $\mathbf{u}-\mathbf{w}$ and $p$, the RIS model is often formulated in a time semi-discrete framework by setting the valve resistance at time level $n$, $R_{\Gamma_i}^n$, as a function of $R_{\Gamma_i}^{n-1}$, $\mathbf{u}^{n-1}-\mathbf{w}^{n-1}$ and $p^{n-1}$ via the procedure reported in Algorithm 1 (see$^{25}$). Furthermore, refractory times and continuous transitions of the resistance are often necessary to mitigate the impact of pressure oscillations observed during the opening and closing of the valves.

### 2.3.2 | Velocity-constrained valve model

We introduce here a novel reduced model for the valves, which, unlike the RIS model, has a sound mathematical formulation. The main idea of this modeling approach consists in enforcing unidirectional flow via an inequality constraint on the normal relative velocity at the interface $\Gamma_i$. This yields the optimality conditions

$$ (\mathbf{u}-\mathbf{w}) \cdot \mathbf{n}_f \ge 0, \quad \llbracket \sigma(\mathbf{u}, p) \mathbf{n}_f \rrbracket \cdot \mathbf{n}_f \ge 0, \quad (\llbracket \sigma(\mathbf{u}, p) \mathbf{n}_f \rrbracket \cdot \mathbf{n}_f)((\mathbf{u}-\mathbf{w}) \cdot \mathbf{n}_f) = 0 \quad \text{on } \Gamma_i, \quad (16) $$
**Algorithm 1** Explicit resistance update in the RIS valve model.

```python
if R_{\Gamma_i}^{n-1} = R_{\text{closed}} (\text{closed valve}) \text{ then}
    if \int_{\Gamma_i} \|p^{n-1}\| > 0 \text{ then}
        R_{\Gamma_i}^n = 0
    else
        R_{\Gamma_i}^n = R_{\text{closed}}^n
end if
else if R_{\Gamma_i}^{n-1} = 0 (\text{open valve}) \text{ then}
    if \int_{\Gamma_i} (u^{n-1} - w^n) \cdot n_f \le 0 \text{ then}
        R_{\Gamma_i}^n = R_{\text{closed}}^n
    else
        R_{\Gamma_i}^n = 0
end if
end if
```

with $i \in \{\text{ao, mv}\}$. Note that the stress jump $\|\sigma(u, p)n_f\| \cdot n_f$ across the interface $\Gamma_i$ corresponds to the (scalar) Lagrange multiplier associated with the inequality constraint $(u-w) \cdot n_f \ge 0$. Thus, the second relation simply enforces that it has to be positive (when the valve is closed) or zero (when the valve is open). Finally, the third corresponds to the so-called complementary condition.

<figure><img src="image_3.png" /><figcaption>**FIGURE 3** Graphical representation of the penalized velocity constraint (19) as a regularization of the ideal graph (dotted).</figcaption></figure>

The optimality conditions (16) can equivalently be rewritten as the following non-linear identity (see, e.g., ⁴⁴ Chapter 6):

$$ \|\sigma(u, p)n_f\| \cdot n_f = \frac{1}{\epsilon} |(u-w) \cdot n_f - \epsilon \|\sigma(u, p)n_f\| \cdot n_f|_{-} \quad \text{on } \Gamma_i \quad (17) $$

for all $\epsilon > 0$. In this work, we consider a penalized version of (17), which amounts to taking $\epsilon$ sufficiently small, so that the stress contribution $\epsilon \|\sigma(u, p)n_f\| \cdot n_f$ in the right-hand side of (17) can be neglected. Furthermore, we assume no tangential friction, viz.,

$$ \|\sigma(u, p)n_f\| \cdot \tau_f = 0 \quad \text{on } \Gamma_i, \quad (18) $$

where $\tau_f$ denotes the tangential basis vectors on $\Gamma_i$. We hence consider the following expression for $f_{\Gamma_i}$:

$$ f_{\Gamma_i}(u - w, p) = -\frac{1}{\epsilon} |(u - w) \cdot n_f|_{-} n_f \quad \text{on } \Gamma_i, \quad (19) $$

which is graphically represented in Figure 3.

There are several fundamental differences between the reduced valve models given by (15) and (18)-(19). First, the open and closing dynamics in (19) are entirely driven by the local sign of $(u-w) \cdot n_f$ on $\Gamma_i$, without the need of the explicit switch between resistances of Algorithm 1. Second, in constrast to (15), the relation (19) can be straightforwardly combined with an implicit time-stepping scheme. This introduces a genuine non-linear term in the discrete problem, which can be handled in the fluid solver with a few Newton iterations (as in the numerical results of Section 4). Finally, it is worth noting that (owing to (18)) the
fluid can slip tangentially along the valve whenever it is closed, which is not the case of (15). This could be avoided by adding a suitable penalty term on the tangential velocity in (18). This variant is not explored in this work.

### 2.3.3 | Open valve configuration

In order to obtain a physiological flow pattern within the left ventricle (see, e.g., 45,10), the open configuration of the mitral valve, denoted here by $\Gamma_{mv}^{open}$ (see Figure 4), must be integrated within the blood flow model (14). To this purpose, a term similar to (15), supported on $\Gamma_{mv}^{open}$ and with associated resistance $R_{\Gamma_{mv}^{open}}$, is added to the momentum conservation equation (14)1.

<figure><img src="image_4.png" /><figcaption>FIGURE 4 Geometry of the mitral valve: closed configuration (in green) and open configuration (in purple).</figcaption></figure>

In the case of the RIS valve model of Section 2.3.1, the value of $R_{\Gamma_{mv}^{open}}$ is simply given in terms of the $R_{\Gamma_{mv}}$ (see, e.g., 25,10), that is,

$$R_{\Gamma_{mv}^{open}} = \begin{cases} 0 & \text{if } R_{\Gamma_{mv}} = R_{\text{closed}}, \\ R_{\text{closed}} & \text{if } R_{\Gamma_{mv}} = 0. \end{cases}$$

In the velocity-constrained model (Section 2.3.2), since there is no switch between resistances, the open configuration is activated in terms of the relative flow rate across the closed configuration, $Q_{mv} \stackrel{\text{def}}{=} \int_{\Gamma_{mv}} (u - w) \cdot n_f$, as follows:

$$R_{\Gamma_{mv}^{open}} = \begin{cases} 0 & \text{if } Q_{mv} \le Q_{\text{min}}, \\ R_{\text{closed}} & \text{if } Q_{mv} > Q_{\text{min}}, \end{cases}$$

where $Q_{\text{min}}$ stands for an user-defined minimal flow tolerance which characterizes the closed valve configuration.

## 2.4 | Model coupling approaches

In this work, cardiac hemodynamics are simulated by considering two different types of modeling coupling approaches on the left ventricle interface $\Gamma_1$: kinematical uncoupling (or one-way) approach and fluid-structure interaction (FSI).

### 2.4.1 | Kinematically uncoupled model

In this modeling approach, the fluid system (14) is solved with a prescribed motion of the fluid cavity $\Omega^f(t)$ through the relations

$$\begin{cases} \mathbf{y}_f = \mathcal{L}(\mathbf{y}_{\text{data}}), & \mathbf{w} = \partial_t \mathbf{y}_f, & \mathbf{A} = \mathbf{I}_{\Omega_f} + \mathbf{y}_f, & \Omega_f(t) = \mathcal{A}(\Omega_f, t), \\ \mathbf{u} = \mathbf{w} & \text{on } \Gamma_1. \end{cases} \qquad (20)$$

Here, $\mathbf{y}_{\text{data}}: \Gamma_1 \times \mathbb{R}^+ \to \mathbb{R}^3$ denotes a given displacement of the left ventricle endocardium $\Gamma_1$ and the symbol $\mathcal{L}$ stands for a (possibly nonlinear) lifting operator from $\Gamma_1$ to $\Omega_f$. In practice, $\mathbf{y}_{\text{data}}$ can be obtained from imaging measurements (see, e.g., 14,15,16,17,18,19,20,21,22) or from electromechanical simulations (see, e.g., 8,9,10,11,12,13). In the latter case, this amounts to solve
<figure><img src="image_5.png" /><figcaption>FIGURE 5 One-way kinematical uncoupling approach in the left ventricle.</figcaption></figure>

(13) with a lumped parameter description of the fluid also in the left heart, namely, by adding the relations

$$ \left\{ \begin{array}{l} F_s \Sigma n_s = -P_1 n_s \quad \text{on } \Gamma_1, \\ \int_{\Gamma_1} J_s u_s \cdot F_s^{-T} n_s = Q_1(P_1, P_{\text{ar},1}, P_{\text{at},1}), \end{array} \right. \qquad (21) $$

where $P_1$, $P_{\text{at},1}$ and $P_{\text{ar},1}$ respectively stand for the pressure in the left ventricle, left atria and aorta, and $Q_1$ denotes a net flux valve function similar to (11). Finally, a two-stage Windkessel model, similar to (12), is considered to describe the dynamics of the unknown pressure $P_{\text{ar},1}$. This approach is schematized in Figure 5.

## 2.4.2 | Coupled fluid-structure interaction model

In this modeling approach, the dynamics of the solid (13) and of the fluid (14) are coupled through the following standard interface coupling conditions on the endocardial interface $\Gamma_1$:

$$ \left\{ \begin{array}{l} y_f = \mathcal{L}(y|_{\Gamma_1}), \quad w = \partial_t y_f, \quad \mathcal{A} = I_{\Omega_f} + y_f, \quad \Omega_f(t) = \mathcal{A}(\Omega_f, t), \\ \mathbf{u} = \mathbf{u}_s \circ \mathcal{A}_f^{-1} \quad \text{on } \Gamma_1(t), \\ F_s \Sigma n_s = -J \hat{\sigma}(\mathbf{u}, p) F^{-T} \hat{\mathbf{n}}_f \quad \text{on } \Gamma_1. \end{array} \right. \qquad (22) $$

The first relation enforces the so-called geometrical compatibility condition, the second corresponds to the kinematic continuity while the third enforces the dynamic equilibrium at the fluid-solid interface $\Gamma_1$. This approach is schematized in Figure 6.

We conclude this section by providing an a priori energy estimate for the coupled system given by (13), (14) and (22).
<figure><img src="image_6.png" /><figcaption>FIGURE 6 FSI coupling approach in the left ventricle.</figcaption></figure>

2.4.2.1 | A priori energy estimate.

We define the total energy of the system as

$$ \mathcal{E}(t) \stackrel{\text{def}}{=} \frac{\rho_f}{2} \|\mathbf{u}\|_{0,\Omega_f(t)}^2 + \frac{\rho_s}{2} \|\mathbf{u}_s\|_{0,\Omega_s}^2 + \int_{\Omega_s} W_e(y) + \frac{1}{2} \int_{\Omega_s} E_s |e_s|^2 + \int_{\Omega_s} U_c + \frac{a}{2} \|y\|_{\Gamma_a \cup \Gamma_b}^2 + \frac{C_{p,r}}{2} |P_{ar,r}|^2 + \frac{C_{d,r}}{2} |P_{d,r}|^2, $$

where the first two terms represent the instantaneous kinetic energy of the fluid and the solid, respectively. The third term corresponds to the elastic energy of the solid and the fourth term to the (passive) elastic energy of the heart fibers. The fifth term accounts for the microscopic elastic energy of the actin-myosin bridges, with $U_c$ denoting the density of elastic energy stored in the sarcomeres and whose dynamics are given by the following relation (see$^{38}$):

$$ \dot{U}_c = -(|u_e| + \alpha_d |\dot{e}_c|) U_c + \tau_c \dot{e}_c + n_0 U_0 |u_e|_+, \quad (23) $$

where $U_0 > 0$ stands for the initial elastic energy in the sarcomeres.

Remark 2. As outlined in$^{38}$, from the relation (3) it follows that $k_c \in [0, k_0]$, provided that $k_c(0) \in [0, k_0]$. Furthermore, it follows from (23) that $U_c \ge \tau_c^2/(2k_c) \ge 0$.

We also introduce the total energy dissipation rate of the system as

$$ \mathcal{D}(t) \stackrel{\text{def}}{=} 2\mu_f \|\boldsymbol{\varepsilon}(\boldsymbol{u})\|_{0,\Omega_f(t)}^2 \, \mathrm{d}s + \sum_{i \in \{\mathrm{mv}, \mathrm{av}\}} \int_{\Gamma_i} f_{\Gamma_i} \cdot \boldsymbol{u} + \int_{\Omega_s} (|u_e| + \alpha_d |\dot{e}_c|) U_c + \int_{\Omega_s} \mu_c |\dot{e}_c|^2 + \int_{\Omega_s} \frac{\partial W_v}{\partial \dot{\boldsymbol{E}}} : \dot{\boldsymbol{E}} \\ + b \|\boldsymbol{u}_s\|_{\Gamma_a \cup \Gamma_b}^2 + \left( K_{\mathrm{ar},r} |P_r - P_{\mathrm{ar},r}|_+^2 + K_{\mathrm{at},r} |P_r - P_{\mathrm{at},r}|_-^2 + \frac{|P_{\mathrm{d},r} - P_{\mathrm{ar},r}|^2}{R_{\mathrm{p},r}} + \frac{|P_{\mathrm{d},r}|^2}{R_{\mathrm{d},r}} \right). $$
Here, the first and second terms correspond to the dissipation induced, respectively, by the viscosity of the fluid and by the reduced modeling of the valves. Note, in particular, that the second term is non negative whenever the valves are fixed (i.e., $w|_{\Gamma_i} = 0$). The third term represents the energy dissipated during the binding and unbinding mechanisms and the fourth the passive dissipation within the sarcomeres. Finally, the last two terms correspond respectively to the passive viscosity of the cardiac tissue and the dissipation induced by the reduced valve and circulation models.

**Theorem 1.** Let $(y, u_s, e_c, k_c, \tau_c, y_f, u, p, P_r, P_{ar,r}, P_{d,r})$ be solution of the coupled system (13), (14) and (22). The following energy identity holds

$$ \mathcal{E}(t) + \int_0^t \mathcal{D}(t) = \mathcal{E}(0) + \int_0^t \int_{\Omega_s} n_0 U_0 |u_e|_+ - \sum_{i \in \{\text{pv,ao}\}} \int_0^t \int_{\Gamma_i} P_i u \cdot n_f + \int_0^t \left( \frac{P_{\text{vs,r}} P_{\text{d,r}}}{R_{\text{d,r}}} - K_{\text{at,r}} P_{\text{at,r}} |P_r - P_{\text{at,r}}|_- \right), \quad (24) $$

where the second term in the right-hand side represents the work produced by the actin-myosin engines and the third and fourth terms the work delivered by the external pressures.

*Proof.* The proof follows by combining energy arguments from$^{46}$ and$^{38}$. We first proceed by taking the scalar product of (13)$_1$ with $u_s$ and of (14)$_1$ with $u$. By integrating the resulting expressions over $\Omega_s$ and $\Omega_f(t)$, respectively, and by adding them and using the interface conditions (22) and the boundary conditions in (13) and (14), we get the following fundamental identity:

$$ \frac{\int_{\Omega_f(t)} \rho_f \partial_t \boldsymbol{u} |_{\boldsymbol{\mathcal{A}}} \cdot \boldsymbol{u} + \rho_f \int_{\Omega_f(t)} (\boldsymbol{u} - \boldsymbol{w}) \cdot \nabla \boldsymbol{u} \cdot \boldsymbol{u} + 2\mu \int_{\Omega_f(t)} |\boldsymbol{\epsilon}(\boldsymbol{u})|^2 + \sum_{i \in \{\text{mv,av}\}} \int_{\Gamma_i} f_{\Gamma_i} \cdot \boldsymbol{u} + \sum_{i \in \{\text{pv,ao}\}} \frac{\rho_f}{2} \int_{\Gamma_i} |\boldsymbol{u} \cdot \boldsymbol{n}_f|_- |\boldsymbol{u}|^2}{T_1} \\ + \rho_s \underbrace{\int_{\Omega_s} \partial_t \boldsymbol{u}_s \cdot \boldsymbol{u}_s + \int_{\Omega_s} \boldsymbol{\Sigma} : \mathbf{D}_y \boldsymbol{E} \cdot \boldsymbol{u}_s + \int_{\Gamma_a \cup \Gamma_b} (ay + bu_s) \cdot \boldsymbol{u}_s + P_r \underbrace{\int_{\Gamma_r} J_s \boldsymbol{F}_s^{-\mathrm{T}} \boldsymbol{n}_s \cdot \boldsymbol{u}_s}_{T_3} = - \sum_{i \in \{\text{pv,ao}\}} \int_{\Gamma_i} P_i \boldsymbol{u} \cdot \boldsymbol{n}_f}_{T_2} = - \sum_{i \in \{\text{pv,ao}\}} \int_{\Gamma_i} P_i \boldsymbol{u} \cdot \boldsymbol{n}_f, \quad (25) $$

where the linear operator $\mathbf{D}_y \boldsymbol{E}$ is defined as

$$ \mathbf{D}_y \boldsymbol{E} \cdot \boldsymbol{v}_s \stackrel{\text{def}}{=} \frac{1}{2} (\boldsymbol{F}_s^{\mathrm{T}} \nabla \boldsymbol{v}_s + (\nabla \boldsymbol{v}_s)^{\mathrm{T}} \boldsymbol{F}_s). $$

We now proceed by estimating each term $T_i$, $i = 1, 2, 3$, separately. Using a change of variables, the identity $\partial_t J = J \widehat{\text{div}} \boldsymbol{w}$, integration by parts in the convective term and (22), it follows that (see, e.g.,$^{46}$):

$$ T_1 = \frac{\rho^f}{2} \frac{d}{dt} \int_{\Omega_f(t)} |\boldsymbol{u}|^2 + 2\mu \int_{\Omega_f(t)} |\boldsymbol{\epsilon}(\boldsymbol{u})|^2 + \sum_{i \in \{\text{mv,av}\}} \int_{\Gamma_i} f_{\Gamma_i} \cdot \boldsymbol{u}. $$

As regards the solid contributions, we first note that owing to (13)$_2$ and to (8) we have

$$ \int_{\Omega_s} \boldsymbol{\Sigma} : \dot{\boldsymbol{E}} = \int_{\Omega_s} \left( \frac{\partial W_e}{\partial \boldsymbol{E}} + \frac{\partial W_v}{\partial \dot{\boldsymbol{E}}} \right) : \dot{\boldsymbol{E}} + \int_{\Omega_s} \sigma_a \boldsymbol{a}^{\mathrm{T}} \dot{\boldsymbol{E}} \boldsymbol{a}. \quad (26) $$

On the other hand, by differentiating in time (8) and by using (7) and then (4) and (5), we get

$$ \int_{\Omega_s} \sigma_a \boldsymbol{a}^{\mathrm{T}} \dot{\boldsymbol{E}} \boldsymbol{a} = \int_{\Omega_s} \sigma_a \dot{\boldsymbol{e}}_a = \int_{\Omega_s} \sigma_s \dot{\boldsymbol{e}}_s + \int_{\Omega_s} \sigma_c \dot{\boldsymbol{e}}_c = \int_{\Omega_s} E_s e_s \dot{\boldsymbol{e}}_s + \int_{\Omega_s} \tau_c \dot{\boldsymbol{e}}_c + \int_{\Omega_s} \mu_c \dot{\boldsymbol{e}}_c^2. $$

By inserting this relation into (26) and using (23), we have

$$ T_2 = \frac{d}{dt} \left( \frac{\rho^s}{2} \int_{\Omega_s} |\boldsymbol{u}_s|^2 + \int_{\Omega_s} W_e + \frac{1}{2} \int_{\Omega_s} E_s |e_s|^2 + \int_{\Omega_s} U_c \right) + \int_{\Omega_s} \frac{\partial W_v}{\partial \dot{\boldsymbol{E}}} : \dot{\boldsymbol{E}} \\ + \frac{a}{2} \|\boldsymbol{y}\|_{\Gamma_a \cup \Gamma_b}^2 + b \|\boldsymbol{u}_s\|_{\Gamma_a \cup \Gamma_b}^2 + \int_{\Omega_s} \mu_c \dot{\boldsymbol{e}}_c^2 + \int_{\Omega_s} (|u_e| + \alpha_d |\dot{\boldsymbol{e}}_c|) U_c - \int_{\Omega_s} n_0 U_0 |u_e|_+ . $$

Finally, for the third term, using (13)$_9$ we get

$$ T_3 = P_r Q_r = P_r |Q_r|_+ - P_r |Q_r|_- = (P_r - P_{\text{ar,r}}) |Q_r|_+ + P_{\text{ar,r}} |Q_r|_+ - P_r |Q_r|_- \quad (27) $$
Owing to the expression of the valve function (11), we have

$$ (P_r - P_{ar,r}) |Q_r|_+ = K_{ar,r} |P_r - P_{ar,r}|_+^2. $$

On the other hand, using a standard energy argument in (12), we obtain

$$ P_{ar,r} |Q_r|_+ = \frac{d}{dt} \left( \frac{C_{p,r}}{2} |P_{ar,r}|^2 + \frac{C_{d,r}}{2} |P_{d,r}|^2 \right) + \frac{|P_{ar,r} - P_{d,r}|^2}{R_{p,r}} + \frac{|P_{d,r}|^2}{R_{d,r}} - \frac{P_{vs,r} P_{d,r}}{R_{d,r}}. $$

For the last term in (27), using (11) again, we have

$$ -P_r |Q_r|_+ = -(P_r - P_{at,r}) |Q_r|_+ + P_{at,r} |Q_r|_+ = K_{at,r} |P_r - P_{at,r}|_+^2 + K_{at,r} P_{at,r} |P_r - P_{at,r}|_+ $$

Finally, the energy equality (24) follows by gathering all the above relations into (25) and by integrating the resulting expression over the time interval (0, t). This completes the proof. □

*Remark 3.* It should be noted that by using the Young's, Cauchy-Schwarz, trace and Korn's inequalities in the right-hand side of (24), we can absorb the velocity and Windkessel pressures contributions into $D(t)$, which yields the a priori energy bound in terms of the problem data

$$ \mathcal{E}(t) + \frac{1}{2} \int_0^t \mathcal{D}(t) \le \mathcal{E}(0) + \int_0^t \int_{\Omega_s} n_0 U_0 |u_e|_+ + \frac{C_T C_K}{\mu} \sum_{i \in \{\text{pv,ao}\}} \int_0^t \int_{\Gamma_i} |P_i|^2 + \frac{1}{2} \int_0^t (|P_{vs,r}|^2 + K_{at,r} |P_{at,r}|^2), $$

where $C_T, C_K > 0$ denote the constants of the trace and Korn's inequalities, respectively. In particular, since $w|_{\Gamma_i} = 0$ for $i \in \{\text{mv, av}\}$, we have $f_{\Gamma_i} \cdot u \ge 0$, so that $D(t) \ge 0$.

# 3 | NUMERICAL APPROXIMATION

This section is devoted to the numerical approximation of the non-linear coupled system (13), (14) and (22).

## 3.1 | Time semi-discretization: Robin-based loosely coupled scheme

The time interval $(0, T)$ is discretized into $N$ sub-intervals $(t_n, t_{n+1})$, where $t_n \stackrel{\text{def}}{=} \tau n, T = t_N$, and $\tau$ denotes the time-step length. In what follows, we shall make use of the following standard notation:

$$ \partial_{\tau} f^{n+1} \stackrel{\text{def}}{=} \frac{1}{\tau} (f^{n+1} - f^n), \quad f^{n+\frac{1}{2}} \stackrel{\text{def}}{=} \frac{f^n + f^{n+1}}{2}, $$

for the first-order backward finite difference and the midpoint value, respectively.

The fluid subproblem (14) is discretized in time using a first-order backward Euler scheme with a semi-implicit treatment of the convective term. For the solid subproblem (13), a modified mid-point scheme is considered for the momentum conservation equation (13)$_1$ with the so-called Gonzalez correction of the stresses, given by (see$^{47}$):

$$ \Sigma_g^{n+\frac{1}{2}} \stackrel{\text{def}}{=} \frac{\partial_{\tau} W_e(E^{n+1}) - \frac{\partial W_e}{\partial E}(E^{n+\frac{1}{2}}) : \dot{E}^{n+\frac{1}{2}}}{|\dot{E}^{n+\frac{1}{2}}|^2} \dot{E}^{n+\frac{1}{2}}, \qquad (28) $$

with the definitions

$$ E^{n+\frac{1}{2}} \stackrel{\text{def}}{=} E(y^{n+\frac{1}{2}}), \quad \dot{E}^{n+\frac{1}{2}} \stackrel{\text{def}}{=} \partial_{\tau} E^{n+1}, \qquad (29) $$

Owing to (1), the discrete passive stress tensor is given by

$$ \Sigma_p^{n+\frac{1}{2}} \stackrel{\text{def}}{=} \frac{\partial W_e}{\partial E}(E^{n+\frac{1}{2}}) + \frac{\partial W_v}{\partial \dot{E}}(\dot{E}^{n+\frac{1}{2}}), $$
where, using the chain rule and (2), the viscous-pseudo potential contribution yields

$$ \frac{\partial W_v}{\partial \dot{\mathbf{E}}} (\dot{\mathbf{E}}^{n+\frac{1}{2}}) \stackrel{\text{def}}{=} \nu D_y \mathbf{E}^{n+\frac{1}{2}} \cdot \mathbf{u}_s^{n+\frac{1}{2}}. \qquad (30) $$

Following³⁸, the dynamics of the active component are discretized in time using a first-order backward Euler scheme in terms of the discrete variables $e_c^{n+1}$, $\varrho^{n+1} \stackrel{\text{def}}{=} \tau_c^{n+1}/\sqrt{k_c^{n+1}}$ and $\gamma^{n+1} \stackrel{\text{def}}{=} \sqrt{k_c^{n+1}}$, which guarantees the positiveness of $k_c^{n+1}$ and the bound on the discrete elastic energy in the sarcomeres (see Remark 4 below).

For the time-discretization of the interface conditions (22), we build on the loosely coupled approach proposed in³⁵ (see also⁴⁸,³⁷). This basically consists in combining an explicit treatment for the geometrical compatibility (22)¹ with the following Robin-Robin discretization of the kinematic/dynamic coupling conditions (22)²,³:

$$ \left\{ \begin{array}{l} F_s^{n+\frac{1}{2}} \overline{\Sigma}^{n+\frac{1}{2}} n_s + \alpha \mathbf{u}_s^{n+\frac{1}{2}} = \alpha \hat{\mathbf{u}}^n - J^n \hat{\sigma}(\mathbf{u}^n, p^n) (\mathbf{F}^n)^{-\mathrm{T}} \hat{\mathbf{n}}_f \quad \text{on } \Gamma_1, \\ \\ J^{n+1} \hat{\sigma}(\mathbf{u}^{n+1}, p^{n+1}) (\mathbf{F}^{n+1})^{-\mathrm{T}} \hat{\mathbf{n}}_f + \alpha \hat{\mathbf{u}}^{n+1} = \alpha \mathbf{u}_s^{n+\frac{1}{2}} + J^n \hat{\sigma}(\mathbf{u}^n, p^n) (\mathbf{F}^n)^{-\mathrm{T}} \hat{\mathbf{n}}_f \quad \text{on } \Gamma_1, \end{array} \right. \qquad (31) $$

where $\alpha > 0$ is the so-called Robin parameter (user-defined). A salient feature of (31), notably with respect to alternative loosely coupled schemes reported in the literature (see, e.g.,⁴⁹,⁵⁰,⁵¹,⁵²,⁵³,³³), is that the scheme can be proved to deliver unconditional stability and nearly optimal accuracy (see³⁵).

The proposed time semi-discrete approximation of the coupled problem (13), (14) and (22) reads therefore as follows: For $n \ge 0$:

1. Solve active solid: Find $\mathbf{y}^{n+1} : \Omega^s \to \mathbb{R}^3$, $\mathbf{u}_s^{n+1} : \Omega^s \to \mathbb{R}^3$, $e_c^{n+1} : \Omega^s \to \mathbb{R}$, $\varrho^{n+1} : \Omega^s \to \mathbb{R}$, $\gamma^{n+1} : \Omega^s \to \mathbb{R}^+$ and $(P_r^{n+1}, P_{ar,r}^{n+1}, P_{d,r}^{n+1}) \in \mathbb{R}^3$ such that $\mathbf{u}_s^{n+\frac{1}{2}} = \partial_\tau \mathbf{y}^{n+1}$ and

$$ \left\{ \begin{array}{l} \rho_s \partial_\tau \mathbf{u}_s^{n+1} - \mathbf{div} \left( \mathbf{F}_s^{n+\frac{1}{2}} \overline{\mathbf{\Sigma}^{n+\frac{1}{2}}} \right) = \mathbf{0} \quad \text{in } \Omega_s, \\ \\ \overline{\mathbf{\Sigma}}^{n+\frac{1}{2}} = \mathbf{\Sigma}_p^{n+\frac{1}{2}} + \mathbf{\Sigma}_g^{n+\frac{1}{2}} + \frac{E_s (a^\mathrm{T} \mathbf{E}^{n+\frac{1}{2}} \mathbf{a} - e_c^{n+\frac{1}{2}})}{(1 + 2e_c^{n+\frac{1}{2}})^2} \mathbf{a} \otimes \mathbf{a} \quad \text{in } \Omega_s, \\ \\ \gamma^{n+1} \varrho^{n+1} + \mu_c \partial_\tau e_c^{n+1} = E_s \frac{(a^\mathrm{T} \mathbf{E}^{n+\frac{1}{2}} \mathbf{a} - e_c^{n+\frac{1}{2}})(1 + 2a^\mathrm{T} \mathbf{E}^{n+\frac{1}{2}} \mathbf{a})}{(1 + 2e_c^{n+\frac{1}{2}})^3} \quad \text{in } \Omega_s, \\ \\ \partial_\tau (\gamma^{n+1})^2 = -(\|u_c^{n+1}\| + \alpha_d |\partial_\tau e_c^{n+1}|)(\gamma^{n+1})^2 + n_0 k_0 \|u_c^{n+1}\|_+ \quad \text{in } \Omega_s, \\ \\ \partial_\tau \varrho^{n+1} = \gamma^{n+1} \partial_\tau e_c^{n+1} + \frac{n_0 \sigma_0 \|u_c^{n+1}\|_+}{\gamma^{n+1}} \left(1 - \frac{k_0}{2\sigma_0 \gamma^{n+1}} \varrho^{n+1}\right) - \frac{1}{2} \left(\|u_c^{n+1}\|_+ + |\partial_\tau e_c^{n+1}|\right) \varrho^{n+1} \quad \text{in } \Omega_s, \\ \\ \mathbf{F}_s^{n+\frac{1}{2}} \overline{\mathbf{\Sigma}}^{n+\frac{1}{2}} \mathbf{n}_s = -(ay^{n+\frac{1}{2}} + bu_{s,h}^{n+\frac{1}{2}}) \quad \text{on } \Gamma_a \cup \Gamma_b, \qquad (32) \\ \\ \mathbf{y}^{n+1} = \mathbf{0} \quad \text{on } \Gamma_v, \\ \\ \mathbf{F}_s^{n+\frac{1}{2}} \overline{\mathbf{\Sigma}}^{n+\frac{1}{2}} \mathbf{n}_s = -P_r^{n+\frac{1}{2}} \mathbf{n}_s \quad \text{on } \Gamma_r, \\ \\ \int_{\Gamma_r} J_s^{n+\frac{1}{2}} \mathbf{u}_s^{n+\frac{1}{2}} \cdot (\mathbf{F}_s^{n+\frac{1}{2}})^{-\mathrm{T}} \mathbf{n}_s = Q_r (P_r^{n+\frac{1}{2}}, P_{ar,r}^{n+\frac{1}{2}}, P_{at,r}^{n+\frac{1}{2}}), \\ \\ \mathbf{F}_s^{n+\frac{1}{2}} \overline{\mathbf{\Sigma}}^{n+\frac{1}{2}} \mathbf{n}_s + \alpha \mathbf{u}_s^{n+\frac{1}{2}} = \alpha \hat{\mathbf{u}}^n - J^n \hat{\sigma}(\mathbf{u}^n, p^n) (\mathbf{F}^n)^{-\mathrm{T}} \hat{\mathbf{n}}_f \quad \text{on } \Gamma_1, \\ \\ C_{p,r} \partial_\tau P_{ar,r}^{n+1} + \frac{P_{ar,r}^{n+\frac{1}{2}} - P_{d,r}^{n+\frac{1}{2}}}{R_{p,r}} = |Q_r^{n+\frac{1}{2}}|_+, \\ \\ C_{d,r} \partial_\tau P_{d,r}^{n+1} + \frac{P_{d,r}^{n+\frac{1}{2}} - P_{ar,r}^{n+\frac{1}{2}}}{R_{p,r}} = \frac{P_{vs,r}(t_{n+\frac{1}{2}}) - P_{d,r}^{n+\frac{1}{2}}}{R_{d,r}}. \end{array} \right. $$

2. Fluid domain update:

$$ \mathbf{y}_f^{n+1} = \mathcal{L}(\mathbf{y}^{n+1}|_{\Gamma_1}), \quad \mathbf{w}^{n+1} = \partial_\tau \mathbf{y}_f^{n+1}, \quad \mathcal{A}^{n+1} = I_{\Omega_f} + \mathbf{y}_f^{n+1} \quad \text{in } \Omega_f^{n+1} = \mathcal{A}^{n+1}(\Omega_f). \qquad (33) $$
3. Solve fluid: Find $\hat{u}^{n+1} : \Omega_f \to \mathbb{R}^3$ and $\hat{p}^{n+1} : \Omega_f \to \mathbb{R}$ such that

$$
\left\{
\begin{array}{ll}
\rho_f \partial_\tau \boldsymbol{u}^{n+1}|_{\boldsymbol{\mathcal{A}}} + \rho_f(\boldsymbol{u}^n - \boldsymbol{w}^{n+1}) \cdot \boldsymbol{\nabla} \boldsymbol{u}^{n+1} - \mathbf{div} \boldsymbol{\sigma}(\boldsymbol{u}^{n+1}, p^{n+1}) + \sum_{i \in \{\text{mv,av}\}} f_{\Gamma_i}^{n+1} \delta_{\Gamma_i} = \mathbf{0} & \text{in } \Omega_f^{n+1}, \\
\mathbf{div} \boldsymbol{u}^{n+1} = \mathbf{0} & \text{in } \Omega_f^{n+1}, \\
\boldsymbol{u}^{n+1} = \mathbf{0} & \text{in } \Gamma_{\text{la}} \cup \Gamma_{\text{aw}}, \\
\boldsymbol{\sigma}(\boldsymbol{u}^{n+1}, p^{n+1}) \boldsymbol{n}_f + \frac{\rho_f}{2} |\boldsymbol{u}^n \cdot \boldsymbol{n}_f|_{-} \boldsymbol{u} = -P_{\text{pv,ao}}^{n+1} \boldsymbol{n}_f & \text{on } \Gamma_{\text{pv}} \cup \Gamma_{\text{ao}}, \\
J^{n+1} \boldsymbol{\hat{\sigma}}(\boldsymbol{u}^{n+1}, p^{n+1})(\boldsymbol{F}^{n+1})^{-T} \boldsymbol{\hat{n}}_f + \alpha \hat{\boldsymbol{u}}^{n+1} = \alpha \boldsymbol{u}_s^{n+1/2} + J^n \boldsymbol{\hat{\sigma}}(\boldsymbol{u}^n, p^n)(\boldsymbol{F}^n)^{-T} \boldsymbol{\hat{n}}_f & \text{on } \Gamma_1.
\end{array}
\right.
$$

Remark 4. As shown in $^{38}$, the scheme (32) preserves the sign properties of $k_c^{n+1} \in (0, n_0k_0)$ as well as the bound $U_c^{n+1} \ge (\tau_c^{n+1})^2/(2k_c^{n+1}) \ge 0$.

## 3.2 | Finite element approximation: fully discrete scheme

For the spatial approximation of the time semi-discrete scheme introduced in the previous section, we adopt a finite element methodology. To this purpose, we let $\mathcal{T}_{s,h}$ and $\mathcal{T}_{f,h}$ be triangulations of the reference domains $\Omega_s$ and $\Omega_f$, respectively, where the subscript $h > 0$ indicates the level of spatial refinement. We assume that $\mathcal{T}_{s,h}$ and $\mathcal{T}_{f,h}$ match on the fluid-solid interface $\Gamma_1$, and that $\mathcal{T}_{f,h}$ is also fitted to the immersed surfaces $\Gamma_{mv}$, $\Gamma_{av}$ and $\Gamma_{mv}^{open}$. We define the following finite element spaces for the fluid velocity, fluid pressure and fluid domain displacement, respectively,

$$V_{f,h} \stackrel{\text{def}}{=} \{ \hat{\mathbf{v}}_{f,h} \in [C^0(\Omega_f)]^3 : \hat{\mathbf{v}}_{f,h} |_{K} \in [\mathbb{P}_1(K)]^3 \quad \forall K \in \mathcal{T}_{f,h} \},$$

$$Q_{f,h} \stackrel{\text{def}}{=} \{ \hat{q}_h \in C^0(\Omega_f \setminus (\Gamma_{mv} \cup \Gamma_{av} \cup \Gamma_{mv}^{\text{open}})) : \hat{q}_h |_{K} \in \mathbb{P}_1(K) \quad \forall K \in \mathcal{T}_{f,h} \},$$

$$W_{f,h} \stackrel{\text{def}}{=} \{ \mathbf{y}_{f,h} \in V_{f,h} : \mathbf{y}_{f,h} = \mathbf{0} \quad \text{on} \quad \Gamma_{mv} \cup \Gamma_{av} \cup \Gamma_{mv}^{\text{open}} \}.$$

Note that the discrete pressure space $Q_{f,h}$ is made of piece-wise affine functions which are globally continuous in $\Omega_f$ except across the immersed valves $\Gamma_{mv}$, $\Gamma_{av}$ and $\Gamma_{mv}^{open}$. This is fundamental in order to guarantee interfacial mass conservation across the valves. For the approximation of the solid displacement and velocity, we consider the finite elements space

$$V_{s,h} \stackrel{\text{def}}{=} \{ \mathbf{v}_{s,h} \in [C^0(\Omega_s)]^3 : \mathbf{v}_{s,h} |_{K} \in [\mathbb{P}_1(K)]^3 \quad \forall K \in \mathcal{T}_{s,h}, \quad \mathbf{v}_{s,h} |_{\Gamma_v} = \mathbf{0} \}. $$

We also introduce the interface matching trace space

$$\mathbf{\Lambda}_h \stackrel{\text{def}}{=} \{ \mathbf{v}_{s,h} |_{\Gamma_1} : \mathbf{v}_{s,h} \in V_{s,h} \} = \{ \hat{\mathbf{v}}_{f,h} |_{\Gamma_1} : \hat{\mathbf{v}}_{f,h} \in V_{f,h} \},$$

which will be used in the formulation of the spatial discrete counterpart of (31). Finally, we denote by $\mathcal{L}_h : \mathbf{\Lambda}_h \to W_{f,h}$ a given discrete counterpart of the lifting operator $\mathcal{L}$ involved in (33).

The dynamics of the active mechanics variables $e_c^{n+1}$, $\varrho^{n+1}$ and $\gamma^{n+1}$ in (32) are discretized in space (element-wise) using a collocation scheme at each quadrature point $q_l, l = 1, \dots, N_{q,h}$, where $N_{q,h} \stackrel{\text{def}}{=} |\mathcal{T}_{s,h}|n_q$ and $n_q$ denotes the number of quadrature points in each element of the triangulation $\mathcal{T}_{s,h}$. This yields a total number of $N_{q,h}$ degrees of freedom for each discrete unknown, viz., given by the arrays

$$e_{c,h}^{n+1} \stackrel{\text{def}}{=} \{ e_{c,l}^{n+1} \}_{l=1}^{N_{q,h}}, \quad \varrho_h^{n+1} \stackrel{\text{def}}{=} \{ \varrho_l^{n+1} \}_{l=1}^{N_{q,h}}, \quad \gamma_h^{n+1} = \{ \gamma_l^{n+1} \}_{l=1}^{N_{q,h}}. \quad (34)$$

The proposed fully discrete approximation of the non-linear coupled problem (13), (14) and (22) is detailed in Algorithm 2. In the solid step (35), the symbols $a_l$, $E_l^{n+1/2}$, $\Sigma_{p,l}^{n+1/2}$, $\Sigma_{g,l}^{n+1/2}$, $\bar{\Sigma}_l^{n+1/2}$ indicate the evaluation of the corresponding field at the quadrature point $q_l, l = 1, \dots, N_{q,h}$. Note that the stress along the fibers at $q_l$ is evaluated in terms of the discrete unknowns (34), which is then integrated via quadrature formula in the discrete principle of virtual work of the solid. It is also worth noting that, at each Newton iteration of (35), the increments associated to the unknowns (34) can be eliminated from the solid tangent system by a static condensation procedure (see, e.g., $^{40}$).
**Algorithm 2** Fully discrete Robin-based loosely coupled scheme.

For $n \ge 0$:

1. Solve active solid: Find $(\mathbf{y}_h^{n+1}, \mathbf{u}_{s,h}^{n+1}, \mathbf{e}_{c,h}^{n+1}, \mathbf{\varrho}_h^{n+1}, \gamma_h^{n+1}, P_r^{n+1}, P_{ar,r}^{n+1}, P_{d,r}^{n+1}) \in \mathbf{V}_{s,h} \times \mathbf{V}_{s,h} \times \mathbb{R}^{N_{q,h}} \times \mathbb{R}^{N_{q,h}} \times (\mathbb{R}^+)^{N_{q,h}} \times \mathbb{R} \times \mathbb{R} \times \mathbb{R}$ with $\mathbf{u}_{s,h}^{n+\frac{1}{2}} = \partial_{\tau} \mathbf{y}_h^{n+1}$ and such that

$$
\left\{
\begin{aligned}
& \rho_s \int_{\Omega_s} \partial_{\tau} \mathbf{u}_{s,h}^{n+1} \cdot \mathbf{v}_{s,h} + \int_{\Omega_s} \bar{\mathbf{\Sigma}}^{n+\frac{1}{2}} : D_y \mathbf{E}^{n+\frac{1}{2}} \cdot \mathbf{v}_{s,h} + \alpha \int_{\Gamma_1} \mathbf{u}_{s,h}^{n+\frac{1}{2}} \cdot \mathbf{v}_{s,h} + \int_{\Gamma_a \cup \Gamma_b} (a_y^{n+\frac{1}{2}} + b \mathbf{u}_{s,h}^{n+\frac{1}{2}}) \cdot \mathbf{v}_{s,h} = \\
& \qquad + \alpha \int_{\Gamma_1} \hat{\mathbf{u}}_h^n \cdot \mathbf{v}_{s,h} - \int_{\Gamma_1} \lambda_h^n \cdot \mathbf{v}_{s,h} - \int_{\Gamma_1} J_s^{n+\frac{1}{2}} P_r^{n+\frac{1}{2}} (\mathbf{F}_s^{n+\frac{1}{2}})^{-T} \mathbf{v}_{s,h} \cdot \mathbf{n}_s, \\
& \bar{\mathbf{\Sigma}}_l^{n+\frac{1}{2}} = \mathbf{\Sigma}_{p,l}^{n+\frac{1}{2}} + \mathbf{\Sigma}_{g,l}^{n+\frac{1}{2}} + \frac{E_s (\mathbf{a}_l^T \mathbf{E}_l^{n+\frac{1}{2}} \mathbf{a}_l - \mathbf{e}_{c,l}^{n+\frac{1}{2}})}{(1 + 2\mathbf{e}_{c,l}^{n+\frac{1}{2}})^2} \mathbf{a}_l \otimes \mathbf{a}_l, \\
& \gamma_l^{n+1} \varrho_l^{n+1} + \mu_c \partial_{\tau} \mathbf{e}_{c,l}^{n+1} = E_s \frac{(\mathbf{a}_l^T \mathbf{E}_l^{n+\frac{1}{2}} \mathbf{a}_l - \mathbf{e}_{c,l}^{n+\frac{1}{2}}) (1 + 2\mathbf{a}_l^T \mathbf{E}_l^{n+\frac{1}{2}} \mathbf{a}_l)}{(1 + 2\mathbf{e}_{c,l}^{n+\frac{1}{2}})^3}, \\
& \partial_{\tau} (\gamma_l^{n+1})^2 = - (|\mathbf{u}_{c,l}^{n+1}| + \alpha_d |\partial_{\tau} \mathbf{e}_{c,l}^{n+1}|) (\gamma_l^{n+1})^2 + n_0 k_0 |\mathbf{u}_{c,l}^{n+1}|_+, \\
& \partial_{\tau} \varrho_l^{n+1} = \gamma_l^{n+1} \partial_{\tau} \mathbf{e}_{c,l}^{n+1} + \frac{n_0 \sigma_0 |\mathbf{u}_{c,l}^{n+1}|_+}{\gamma_l^{n+1}} \left(1 - \frac{k_0}{2\sigma_0 \gamma_l^{n+1}} \varrho_l^{n+1}\right) - \frac{1}{2} (|\mathbf{u}_{c,l}^{n+\frac{1}{2}}| + |\partial_{\tau} \mathbf{e}_{c,l}^{n+1}|) \varrho_l^{n+1}, \\
& \int_{\Gamma_1} J_s^{n+\frac{1}{2}} \mathbf{u}_{s,h}^{n+\frac{1}{2}} \cdot (\mathbf{F}_s^{n+\frac{1}{2}})^{-T} \mathbf{n}_s = Q_r (P_r^{n+\frac{1}{2}}, P_{ar,r}^{n+\frac{1}{2}}, P_{at,r}^{n+\frac{1}{2}}), \\
& C_{p,r} \partial_{\tau} P_{ar,r}^{n+1} + \frac{P_{ar,r}^{n+\frac{1}{2}} - P_{d,r}^{n+\frac{1}{2}}}{R_{p,r}} = |Q_r^{n+\frac{1}{2}}|_+, \\
& C_{d,r} \partial_{\tau} P_{d,r}^{n+1} + \frac{P_{d,r}^{n+\frac{1}{2}} - P_{ar,r}^{n+\frac{1}{2}}}{R_{p,r}} = \frac{P_{vs,r}(t_{n+\frac{1}{2}}) - P_{d,r}^{n+\frac{1}{2}}}{R_{d,r}},
\end{aligned}
\right.
\quad (35)
$$

for all $\mathbf{v}_{s,h} \in \mathbf{V}_{s,h}$ and $l = 1, \dots, N_{q,h}$.

2. Update fluid domain:

$$ \mathbf{y}_{f,h}^{n+1} = \mathcal{L}_h(\mathbf{y}_h^n|_{\Gamma_1}), \quad \mathbf{w}_h^{n+1} = \partial_{\tau} \mathbf{y}_{f,h}^{n+1}, \quad \mathcal{A}_h^{n+1} = \mathbf{I}_{\Omega_f} + \mathbf{y}_{f,h}^{n+1}, \quad \Omega_f^{n+1} = \mathcal{A}_h^{n+1}(\Omega_f). \quad (36) $$

3. Solve fluid: Find $(\hat{\mathbf{u}}_h^{n+1}, \hat{P}_h^{n+1}) \in \mathbf{V}_{f,h} \in \mathcal{Q}_{f,h}$ such that

$$
\rho_f \partial_{\tau} \int_{\Omega_f^{n+1}} \mathbf{u}_h^{n+1} \cdot \mathbf{v}_{f,h} - \rho_f \int_{\Omega_f^{n+1}} (\mathrm{div} \, \mathbf{w}_h^{n+1}) \mathbf{u}_h^{n+1} \cdot \mathbf{v}_{f,h} + a_{\Omega_f^{n+1}} (\mathbf{u}_h^n, \mathbf{w}_h^{n+1}; (\mathbf{u}_h^{n+1}, \mathbf{p}_h^{n+1}), (\mathbf{v}_{f,h}, q_h)) \\
+ s_{GCL} (\mathbf{u}_h^n, \mathbf{w}_h^{n+1}; \mathbf{u}_h^{n+1}, \mathbf{v}_{f,h}) + \sum_{i \in \{mv, av\}} \int_{\Gamma_i} f_{\Gamma_i} (\mathbf{u}_h^{n+1} - \mathbf{w}_h^{n+1}, \mathbf{p}^{n+1}) \cdot \mathbf{v}_{f,h} + \frac{\rho_f}{2} \int_{\Gamma_{pv,ao}} |\mathbf{u}_h^n \cdot \mathbf{n}_f|_- \mathbf{u}_h^{n+1} \cdot \mathbf{v}_{f,h} \\
+ \alpha \int_{\Gamma_1} \hat{\mathbf{u}}_h^{n+1} \cdot \hat{\mathbf{v}}_{f,h} = \alpha \int_{\Gamma_1} \mathbf{u}_{s,h}^{n+\frac{1}{2}} \cdot \hat{\mathbf{v}}_{f,h} + \int_{\Gamma_1} \lambda_h^n \cdot \hat{\mathbf{v}}_{f,h} - \int_{\Gamma_{pv} \cup \Gamma_{ao}} P_{pv,ao}^{n+1} \mathbf{v}_{f,h} \cdot \mathbf{n}_f \quad (37)
$$

for all $(\hat{\mathbf{v}}_{f,h}, \hat{q}_h) \in \mathbf{V}_{f,h} \times \mathcal{Q}_{f,h}$.

4. Interfacial fluid stress update: Set $\lambda_h^{n+1} \in \mathbf{\Lambda}_h$ with

$$ \lambda_h^{n+1} = \lambda_h^n + \alpha (\mathbf{u}_{s,h}^{n+\frac{1}{2}} - \hat{\mathbf{u}}_h^{n+1}) \quad \text{on } \Gamma_1. \quad (38) $$

As in$^{35}$, the fluid trilinear form the fluid sub-problem (37) is given by the relation

$$
\begin{align}
a_{\Omega_f^{n+1}} (\mathbf{u}_h^n, \mathbf{w}_h^{n+1}; (\mathbf{u}_h^{n+1}, \mathbf{p}_h^{n+1}), (\mathbf{v}_{f,h}, q_h)) &\stackrel{\text{def}}{=} \rho_f \int_{\Omega_f^{n+1}} (\mathbf{u}_h^n - \mathbf{w}_h^{n+1}) \cdot \nabla \mathbf{u}_h^{n+1} \cdot \mathbf{v}_{f,h} + \frac{\rho_f}{2} \int_{\Omega_f^{n+1}} (\mathrm{div} \, \mathbf{u}_h^n) \mathbf{u}_h^{n+1} \cdot \mathbf{v}_{f,h} \nonumber \\
&\quad - \frac{\rho_f}{2} \int_{\Gamma_1^{n+1}} (\mathbf{u}_h^n - \mathbf{u}_s^{n+\frac{1}{2}}) \cdot \mathbf{n}_f \mathbf{u}_h^{n+1} \cdot \mathbf{v}_{f,h} \nonumber \\
&\quad + 2\mu_f \int_{\Omega_f^{n+1}} \epsilon(\mathbf{u}_h^{n+1}) : \epsilon(\mathbf{v}_{f,h}) - \int_{\Omega_f^{n+1}} \mathbf{p}^{n+1} \mathrm{div} \, \mathbf{v}_{f,h} + \int_{\Omega_f^{n+1}} q_h \mathrm{div} \, \mathbf{u}_h^n \nonumber \\
&\quad + s_{f,h} (\mathbf{u}_h^n, \mathbf{w}_h^{n+1}; (\mathbf{u}_h^{n+1}, \mathbf{p}^{n+1}), (\mathbf{v}_h^f, q_h)). \tag{39}
\end{align}
$$
The second term in the right-hand side corresponds to the so-called Temam's trick (see, e.g.,⁵⁴), to cope with the fact that the discrete fluid velocities are not divergence-free. The third term is a strongly consistent interfacial stabilization, introduced in³⁵, which serves to balance the contributions of the convective term on the interface, by noting that $\mathbf{u}^n|_{\Gamma_1} \neq \mathbf{w}^{n+1}|_{\Gamma_1} = \mathbf{u}_s^{n+\frac{1}{2}}|_{\Gamma_1}$. Finally, the last term in (39) corresponds to the Streamline Upwind Petrov-Galerkin (SUPG)/Pressure-Stabilized Petrov-Galerkin (PSPG) stabilization (see, e.g.,⁵⁵), which guarantees robustness for high local Reynolds numbers and with respect to the lack of inf-sup stability in the fluid velocity/pressure discrete spaces $V_{f,h}/Q_{f,h}$.

On the other hand, to cope with the stability issues related to the lack of discrete geometric conservation law (see, e.g.,⁵⁶), the fluid sub-problem (37) includes the following weakly consistent stabilization term (from⁵⁷ Chapter 5):

$$s_{\text{GCL}}^{n+1}(\mathbf{u}_h^n, \mathbf{w}_h^{n+1}; \mathbf{u}_h^{n+1}, \mathbf{v}_{f,h}) \stackrel{\text{def}}{=} -\frac{\rho_f}{2\tau} \left( \int_{\Omega_f^{n+1}} \mathbf{u}_h^{n+1} \cdot \mathbf{v}_{f,h} - \int_{\Omega_f^n} \mathbf{u}_h^{n+1} \cdot \mathbf{v}_{f,h} \right) + \frac{\rho_f}{2} \int_{\Omega_f^{n+1}} (\text{div } \mathbf{w}^{n+1}) \mathbf{u}_h^{n+1} \cdot \mathbf{v}_{f,h}. \quad (40)$$

Finally, it is also worth noting that the relation (38) in Algorithm 2 is nothing but the spatial discrete counterpart of (31)$_2$, in which the intermediate variable $\lambda_h^{n+1}$ approximates the interfacial fluid stress on the reference configuration $J^{n+1}\bar{\sigma}(\mathbf{u}^{n+1}, p^{n+1})(\mathbf{F}^{n+1})^{-T}\hat{\mathbf{n}}_f$ as a variational residual (see³⁵,³⁷).

## 3.3 | Discrete energy bound: unconditional stability

The purpose of this section is to derive an energy bound for the numerical approximation of the coupled problem (13), (14) and (22) provided by Algorithm 2. For $n \ge 1$, the discrete energy is defined as

$$\mathcal{E}_h^n \stackrel{\text{def}}{=} \frac{\rho_f}{2} \|\mathbf{u}_h^n\|_{0,\Omega_f^n}^2 + \frac{\rho_s}{2} \|\mathbf{u}_{s,h}^n\|_{0,\Omega_s}^2 + \int_{\Omega_s} W_c(\mathbf{E}^n) + \frac{\tau}{2} \left( \alpha \|\hat{\mathbf{u}}_h^n\|_{0,\Gamma_1}^2 + \frac{1}{\alpha} \|\lambda_h^n\|_{0,\Gamma_1}^2 \right) + \int_{\Omega_s} \frac{E_s}{2} |e_s^n|^2 + \int_{\Omega_s} U_c^n \\ + \frac{a}{2} \|\mathbf{y}_h^n\|_{\Gamma_a \cup \Gamma_b}^2 + \frac{C_{p,r}}{2} |P_{ar,r}^n|^2 + \frac{C_{d,r}}{2} |P_{d,r}^n|^2,$$

and the discrete dissipation rate as

$$\mathcal{D}_h^n \stackrel{\text{def}}{=} 2\mu_f \|\epsilon(\mathbf{u}_h^n)\|_{0,\Omega_f^n}^2 + \sum_{i \in \{\text{mv,av}\}} \int_{\Gamma_i} f_{\Gamma_i}^n \cdot \mathbf{u}_h^n + \alpha \|\mathbf{u}_{s,h}^{n-\frac{1}{2}} - \hat{\mathbf{u}}_h^{n-1}\|_{0,\Gamma_1}^2 + \int_{\Omega_s} \frac{\partial W_v}{\partial \dot{\mathbf{E}}} (\dot{\mathbf{E}}^{n-\frac{1}{2}}) : \dot{\mathbf{E}}^{n-\frac{1}{2}} \\ + \int_{\Omega_s} (|u_e^n| + \alpha_d |\partial_\tau e_c^n|) U_c^n + \int_{\Omega_s} \mu_c |\partial_\tau e_c^n|^2 + b \|\mathbf{u}_s^n\|_{\Gamma_a \cup \Gamma_b}^2 + K_{\text{ar,r}} |P_r^{n-\frac{1}{2}} - P_{\text{ar,r}}^{n-\frac{1}{2}}|^2 \\ + K_{\text{at,r}} |P_r^{n-\frac{1}{2}} - P_{\text{at,r}}(t_{n-\frac{1}{2}})|^2_+ + \frac{|P_{d,r}^{n-\frac{1}{2}} - P_{\text{ar,r}}^{n-\frac{1}{2}}|^2}{R_{p,r}} + \frac{|P_{d,r}^{n-\frac{1}{2}}|^2}{R_{d,r}}.$$

The following result establishes that, up to some numerical dissipation, Algorithm 2 satisfies an energy estimate similar to the one provided by Theorem 1 for the continuous coupled problem (13), (14) and (22).

**Theorem 2.** Let $\{(\mathbf{y}_h^{n+1}, \mathbf{u}_{s,h}^{n+1}, \mathbf{e}_c^{n+1}, \boldsymbol{\theta}^{n+1}, \gamma^{n+1}, \mathbf{u}_h^{n+1}, p_h^{n+1}, \lambda_h^{n+1})\}_{n=0}^{N-1} \subset V_{s,h} \times V_{s,h} \times \mathbb{R}^{N_{q,h}} \times \mathbb{R}^{N_{q,h}} \times \mathbb{R}^{N_{q,h}} \times V_{f,h} \times Q_{f,h} \times \Lambda_h$ be given by Algorithm 2. The following energy estimate holds:

$$\mathcal{E}_h^n + \sum_{m=0}^{n-1} \tau \mathcal{D}_h^{m+1} \le \mathcal{E}_h^0 + \sum_{m=0}^{n-1} \tau \int_{\Omega} U_0 |u_e^{m+1}|_+ - \sum_{i \in \{\text{pv,ao}\}} \sum_{m=0}^{n-1} \tau \int_{\Gamma_i} P_i(t_{m+1}) \mathbf{u}_h^{m+1} \cdot \mathbf{n}_f \\ + \sum_{m=0}^{n-1} \tau \frac{P_{\text{vs,r}}(t_{m+\frac{1}{2}}) P_{\text{d,r}}^{m+\frac{1}{2}}}{R_{\text{d,r}}} - \sum_{m=0}^{n-1} \tau K_{\text{at,r}} P_{\text{at,r}}(t_{m+\frac{1}{2}}) |P_r^{m+\frac{1}{2}} - P_{\text{at,r}}(t_{m+\frac{1}{2}})|_-. \quad (41)$$
Proof. The proof builds on the energy arguments reported in ³⁵ and ³⁸. We first proceed by taking $v_{s,h} = u_{s,h}^{n+\frac{1}{2}}$ in (35). Owing to the relation $D_y E^{n+\frac{1}{2}} \cdot u_{s,h}^{n+\frac{1}{2}} = \dot{E}^{n+\frac{1}{2}} \stackrel{\text{def}}{=} \partial_\tau E^{n+1}$ and by using (28), the stiffness term yields

$$ 
\begin{aligned}
\int_{\Omega_s} \bar{\Sigma}^{n+\frac{1}{2}} : D_y E^{n+\frac{1}{2}} \cdot u_{s,h}^{n+\frac{1}{2}} &= \int_{\Omega_s} \frac{\partial W_e}{\partial E} (E^{n+\frac{1}{2}}) : \dot{E}^{n+\frac{1}{2}} + \int_{\Omega_s} \Sigma_g^{n+\frac{1}{2}} : \dot{E}^{n+\frac{1}{2}} + \int_{\Omega_s} \frac{\partial W_v}{\partial \dot{E}} (\dot{E}^{n+\frac{1}{2}}) : \dot{E}^{n+\frac{1}{2}} \\
&\quad + \int_{\Omega_s} \sigma_a^{n+\frac{1}{2}} \mathbf{a} \otimes \mathbf{a} : \partial_\tau E^{n+1}, \qquad (42) \\
&= \int_{\Omega_s} \partial_\tau W_e (E^{n+1}) + \int_{\Omega_s} \frac{\partial W_v}{\partial \dot{E}} (\dot{E}^{n+\frac{1}{2}}) : \dot{E}^{n+\frac{1}{2}} + \int_{\Omega_s} \sigma_a^{n+\frac{1}{2}} \partial_\tau e_a^{n+1},
\end{aligned}
$$

where, according to (10) and (8), we have used the notations

$$ \sigma_{a,l}^{n+\frac{1}{2}} \stackrel{\text{def}}{=} E_s \frac{e_{a,l}^{n+\frac{1}{2}} - e_{c,l}^{n+\frac{1}{2}}}{(1 + 2e_{c,l}^{n+\frac{1}{2}})^2}, \quad e_{a,l}^{n+\frac{1}{2}} \stackrel{\text{def}}{=} a_l^T \mathbf{E}_i^{n+\frac{1}{2}} \mathbf{a}_l $$

for $l = 1, \dots, N_{q,h}$. It hence only remains to estimate the last term in (42) coming from the active contribution. It should be noted that a little abuse of notation is made in this term, where the integral over $\Omega_s$ has to be read as a quadrature formula. To this purpose, we also introduce the notations

$$ e_{s,l}^{n+1} \stackrel{\text{def}}{=} \frac{e_{a,l}^{n+1} - e_{c,l}^{n+1}}{1 + 2e_{c,l}^{n+1}}, \quad \sigma_{s,l}^{n+\frac{1}{2}} \stackrel{\text{def}}{=} E_s e_{s,l}^{n+\frac{1}{2}}, \quad \sigma_{s,l}^{n+\frac{1}{2}} \stackrel{\text{def}}{=} \gamma_l^{n+1} \varrho_l^{n+1} + \mu_c \dot{e}_{c,l}^{n+\frac{1}{2}}, $$

which yields the stress identities

$$ \sigma_{a,l}^{n+\frac{1}{2}} = \frac{\sigma_{c,l}^{n+\frac{1}{2}}}{1 + 2e_{s,l}^{n+\frac{1}{2}}} = \frac{\sigma_{s,l}^{n+\frac{1}{2}}}{1 + 2e_{c,l}^{n+\frac{1}{2}}}. $$

Additionally, we also introduce the collocated discrete elastic energy stored in the sarcomeres

$$ \partial_{\tau} U_{c,l}^{n+1} = \gamma_l^{n+1} \varrho_l^{n+1} \partial_{\tau} e_{c,l}^{n+1} - (|u_{c,l}^{n+1}| + \alpha_d |\partial_{\tau} e_{c,l}^{n+1}|) U_{c,l}^{n+1} + n_0 U_0 |u_{c,l}^{n+1}|_+ . \quad (43) $$

By combining the above relations, we have

$$ 
\begin{aligned}
\sigma_{a,l}^{n+\frac{1}{2}} &= \sigma_{s,l}^{n+\frac{1}{2}} \partial_{\tau} e_{s,l}^{n+1} + \sigma_{c,l}^{n+\frac{1}{2}} \partial_{\tau} e_{c,l}^{n+1} \\
&= \frac{E_s}{2} \partial_{\tau} |e_{s,l}^{n+1}|^2 + \gamma_l^{n+1} \varrho_l^{n+1} \partial_{\tau} e_{c,l}^{n+1} + \mu_c |\partial_{\tau} e_{c,l}^{n+1}|^2.
\end{aligned}
$$

Therefore, by using (43), we finally get

$$ \int_{\Omega_s} \sigma_a^{n+\frac{1}{2}} \partial_\tau e_a^{n+1} = \int_{\Omega_s} \frac{E_s}{2} \partial_\tau |e_s^{n+1}|^2 + \int_{\Omega_s} \mu_c |\partial_\tau e_s^{n+1}|^2 \\ + \int_{\Omega_s} \left( \partial_\tau U_c^{n+1} + (|u_e^{n+1}| + \alpha_d |\partial_\tau e_c^{n+1}|) U_c^{n+1} - n_0 U_0 |u_e^{n+1}|_+ \right). $$

By gathering all the above contributions, the energy balance of the solid system yields

$$ 
\begin{aligned}
& \frac{\rho_s}{2} \partial_\tau \|\mathbf{u}_{s,h}^{n+1}\|_{0,\Omega_s}^2 + \partial_\tau \int_{\Omega_s} W_e(\mathbf{E}^{n+1}) + \partial_\tau \int_{\Omega} \left(\frac{E_s}{2}|e_s^{n+1}|^2 + U_c^{n+1}\right) + \int_{\Omega_s} \frac{\partial W_v}{\partial \dot{E}}(\dot{E}^{n+\frac{1}{2}}) : \dot{E}^{n+\frac{1}{2}} + \int_{\Omega_s} \mu_c |\partial_\tau e_s^{n+1}|^2 \\
& \qquad + \int_{\Omega} (|u_e^{n+1}| + \alpha_d |\partial_\tau e_c^{n+1}|) U_c^{n+1} + \alpha \int_{\Gamma_1} (\mathbf{u}_{s,h}^{n+1} - \hat{\mathbf{u}}_h^n) \cdot \dot{\mathbf{y}}_h^{n+1} - \int_{\Gamma_1} \lambda_h^n \cdot \dot{\mathbf{y}}_h^{n+1} \\
&= \int_{\Omega} n_0 U_0 |u_e^{n+1}|_+ - \int_{\Gamma_\tau} \int_{s}^{n+\frac{1}{2}} P_r^{n+\frac{1}{2}}(\mathbf{F}_s^{n+\frac{1}{2}})^{-T} \mathbf{u}_{s,h}^{n+1} \cdot \mathbf{n}_s.
\end{aligned}
\quad (44) $$

Finally, for the last term on the right endocardial boundary, we have

$$ \int_{\Gamma_r} J_s^{n+\frac{1}{2}} P_r^{n+\frac{1}{2}} (\mathbf{F}_s^{n+\frac{1}{2}})^{-T} \mathbf{u}_{s,h}^{n+1} \cdot \mathbf{n}_s = P_r^{n+\frac{1}{2}} Q_r^{n+\frac{1}{2}}, $$
so that, by proceeding in a similar fashion as in the proof of Theorem 1, we have

$$ P_r^{n+\frac{1}{2}} Q_r^{n+\frac{1}{2}} = K_{ar,r} |P_r^{n+\frac{1}{2}} - P_{ar,r}^{n+\frac{1}{2}}|_+^2 + \frac{C_{p,r}}{2} \partial_\tau |P_{ar,r}^{n+1}|^2 + \frac{C_{d,r}}{2} \partial_\tau |P_{d,r}^{n+1}|^2 + \frac{|P_{ar,r}^{n+\frac{1}{2}} - P_{d,r}^{n+\frac{1}{2}}|^2}{R_{p,r}} + \frac{|P_{d,r}^{n+\frac{1}{2}}|^2}{R_{d,r}} \\ + K_{at,r} |P_r^{n+\frac{1}{2}} - P_{at,r}^{n+\frac{1}{2}}|_-^2 + K_{at,r} P_{at,r}^{n+\frac{1}{2}} |P_r^{n+\frac{1}{2}} - P_{at,r}^{n+\frac{1}{2}}|_- - \frac{P_{vs,r}^{n+\frac{1}{2}} P_{d,r}^{n+\frac{1}{2}}}{R_{d,r}}. $$

The energy balance of the fluid subproblem is obtained by taking $(\hat{w}_{f,h}, \hat{q}_h) = (\hat{u}_h^{n+1}, \hat{p}_h^{n+1})$ in (37). By integrating by parts in the convective term, using the relation $w_h^{n+1}|_{\Gamma_1} = u_{s,h}^{n+\frac{1}{2}}|_{\Gamma_1}$ and the fact that $s_{f,h} \ge 0$, we have

$$ a_{\Omega^{f,n}} (\mathbf{u}_h^{n-1}, \mathbf{w}_h^n; (\mathbf{u}_h^n, \mathbf{p}_h^n), (\mathbf{u}_h^n, \mathbf{p}_h^n)) \geq 2\mu \|\epsilon(\mathbf{u}_h^n)\|_{0,\Omega^{f,n}}^2 + \frac{\rho^f}{2} \int_{\Omega^{f,n}} (\operatorname{div} \mathbf{w}_h^n) |\mathbf{u}_h^n|^2. $$

By combining this bound with (40) and (38), we finally get the following energy estimate for the fluid:

$$ \frac{\rho_f}{2} \partial_\tau \|\mathbf{u}_h^{n+1}\|_{0,\Omega_f^{n+1}}^2 + 2\mu_f \|\epsilon(\mathbf{u}_h^{n+1})\|_{0,\Omega_f^{n+1}}^2 + \sum_{i \in \{\text{mv,av}\}} \int_{\Gamma_i} f_{\Gamma_i}(\mathbf{u}_h^{n+1} - \mathbf{w}_h^{n+1}, p^{n+1}) \cdot \mathbf{u}_h^{n+1} \\ \leq \int_{\Gamma_1} \lambda_h^{n+1} \hat{\mathbf{u}}_h^{n+1} - \int_{\Gamma_{pv} \cup \Gamma_{ao}} P_{pv,ao}^{n+1} \mathbf{u}_h^{n+1} \cdot \mathbf{n}_f. \quad (45) $$

We now proceeding by adding the relations (44) and (45), and after some mathematical manipulations, it only remains to control the following interfacial terms on the right-hand side:

$$ \underbrace{\alpha \langle \hat{\mathbf{u}}_h^n - \mathbf{u}_{s,h}^{n+\frac{1}{2}}, \mathbf{u}_{s,h}^{n+\frac{1}{2}} \rangle - \langle \lambda_h^n, \mathbf{u}_{s,h}^{n+\frac{1}{2}} \rangle + \langle \lambda_h^{n+1}, \mathbf{u}_h^{n+1} \rangle}_{T_1}. $$

These terms can be handled as in$^{35}$, which yields

$$ T_1 = \frac{\alpha}{2} \left( \|\hat{\mathbf{u}}_h^n\|_{0,\Gamma_1}^2 - \|\hat{\mathbf{u}}_h^{n+1}\|_{0,\Gamma_1}^2 \right) + \frac{1}{2\alpha} \left( \|\lambda_h^n\|_{0,\Gamma_1}^2 - \|\lambda_h^{n+1}\|_{0,\Gamma_1}^2 \right) - \alpha \|\mathbf{u}_{s,h}^{n+\frac{1}{2}} - \hat{\mathbf{u}}_h^n\|_{0,\Gamma_1}^2. $$

The energy bound (41) then follows by multiplying by $\tau$ the resulting expression and then summing over $m = 0, \dots, n-1$, which completes the proof. $\square$

We conclude this section by noting that, using an argument as in Remark 3, the following energy bound can be inferred from (41):

$$ \mathcal{E}_h^n + \sum_{m=0}^{n-1} \frac{\tau}{2} \mathcal{D}_h^{m+1} \leq \mathcal{E}_h^0 + \sum_{m=0}^{n-1} \tau \int_{\Omega} U_0 |\mathbf{u}_e^{m+1}|_+ + \frac{C_T C_K}{\mu} \sum_{i \in \{\text{pv,ao}\}} \sum_{m=0}^{n-1} \tau \|P_i(t_{m+1})\|_{0,\Gamma_i}^2 + \sum_{m=0}^{n-1} \frac{\tau}{2} |P_{vs,r}(t_{m+\frac{1}{2}})|^2 \\ + \sum_{m=0}^{n-1} \frac{\tau K_{at,r}}{2} |P_{at,r}(t_{m+\frac{1}{2}})|^2, $$

which guarantees the unconditional stability of Algorithm 2 in the energy norm.

## 4 | NUMERICAL RESULTS

In this section, we present the numerical results of the left heart hemodynamics obtained using both the kinematical uncoupling approach of Section 4.1 and the fluid-structure interaction model of Section 4.2.2. The computational geometries used in the numerical simulations were obtained from the Zygote‡ heart model. For the solid we isolate the right and left ventricles and the endocardial surfaces are smoothed, with the 3-matic§ software, by removing trabeculae and papillary muscles (see Figure 7(a)).

<footnote>‡ https://www.zygote.com</footnote>

<footnote>§ https://www materialise.com</footnote>
For the fluid, the endocardial surface of the left ventricle is first extracted and a volume mesh of the left ventricular cavity is constructed, which ensures the interface matching nature of the fluid and solid meshes (see Figure 7(b)). The endocardial surfaces of the left atrium and of the ascending aorta are then extracted from the Zygote heart model, and volume meshes of these cavities are then generated, by incorporating the closed configuration of the immersed mitral $Γ_{mv}$ and aortic $Γ_{av}$ valves, and the open configuration of the mitral valve $Γ_{mv}^{open}$ in the fluid domain. These immersed surfaces are duplicated in the fluid mesh (cracked mesh), which facilitates the introduction of strong discontinuities in the pressure approximation across the valves. The solid and fluid volume meshes are made of tetrahedra, and obtained using the 3D Gmsh meshing software (see$^{58}$). More specifically, the solid mesh is composed of 45,584 vertices and 106,747 tetrahedra, while the fluid mesh is composed of 28,090 vertices and 148,333 tetrahedra. Fiber directions in the solid bi-ventricular mesh are synthetically generated using the automatic procedure

<figure><img src="image_7.png" /><figcaption>**FIGURE 8** Illustration of the fiber distribution in the right and left ventricles.</figcaption></figure>

reported in$^{59,60}$, where distance maps are computed in order to prescribe fiber elevation angles varying linearly from -65 degrees (epicardium) to +65 degrees (endocardium) across the thickness. The resulting fiber distribution is illustrated in Figure 8.

The numerical computations have been performed using the fluid solver of the FELiScE¶ library and the cardiac mechanics solver of the MoReFEM# library. Owing to the genuine partitioned nature of Algorithm 2, the FSI coupling is performed by linking the two softwares in a segregated fashion, via an external coupling interface based on the ZeroMQ|| message passing library.

<footnote>¶ https://gitlab.inria.fr/felisce</footnote>

<footnote># https://gitlab.inria.fr/MoReFEM</footnote>

<footnote>|| https://zeromq.org</footnote>
## 4.1 Kinematical uncoupling

In this section, we use the kinematical uncoupling modeling approach to delve into the influence on the left heart hemodynamics of the reduced valve models introduced in Sections 2.3.1 and 2.3.2. First, the 3D-0D cardiac mechanics problem (13) and (21) is solved with an analytical electrical input $u_e$, which mimics the FitzHugh-Nagumo dynamics (see$^{61,62}$), and initiates at the apex of the heart. Afterwards, the fluid problem (14) and (20) is solved separately, using the left endocardial motion provided by the 3D-0D cardiac mechanics solution, for both the RIS and the velocity-constrained reduced valve models.

<table><caption>TABLE 1 Passive physical parameters of the cardiac mechanics model.</caption><thead><tr><th>Parameter</th><th>Left Ventricle</th><th>Right Ventricle</th></tr></thead><tbody><tr><td>μ<sub>1</sub> (Pa)</td><td>0.0</td><td>0.0</td></tr><tr><td>μ<sub>2</sub> (Pa)</td><td>0.0</td><td>0.0</td></tr><tr><td>C<sub>0</sub> (Pa)</td><td>1.9 · 10<sup>3</sup></td><td>1.9 · 10<sup>3</sup></td></tr><tr><td>C<sub>1</sub></td><td>1.1 · 10<sup>-1</sup></td><td>1.1 · 10<sup>-1</sup></td></tr><tr><td>C<sub>2</sub> (Pa)</td><td>1.9 · 10<sup>3</sup></td><td>1.9 · 10<sup>3</sup></td></tr><tr><td>C<sub>3</sub></td><td>1.1 · 10<sup>-1</sup></td><td>1.1 · 10<sup>-1</sup></td></tr><tr><td>ν (N s m<sup>-1</sup>)</td><td>70</td><td>70</td></tr><tr><td>κ (Pa)</td><td>2 · 10<sup>5</sup></td><td>2 · 10<sup>5</sup></td></tr></tbody></table>

<table><caption>TABLE 2 Active physical parameters of the cardiac mechanics model.</caption><thead><tr><th>Parameter</th><th>Left Ventricle</th><th>Right Ventricle</th></tr></thead><tbody><tr><td>α<sub>d</sub></td><td>1.5</td><td>1.5</td></tr><tr><td>σ<sub>0</sub> (Pa)</td><td>4.6 · 10<sup>5</sup></td><td>3.72 · 10<sup>5</sup></td></tr><tr><td>E<sub>s</sub> (Pa)</td><td>3 · 10<sup>7</sup></td><td>3 · 10<sup>7</sup></td></tr><tr><td>k<sub>0</sub> (Pa)</td><td>1 · 10<sup>5</sup></td><td>1 · 10<sup>5</sup></td></tr><tr><td>η (N s m<sup>-1</sup>)</td><td>70</td><td>70</td></tr></tbody></table>

<table><caption>TABLE 3 Parameters of the Windkessel model.</caption><thead><tr><th>Parameter</th><th>Left Ventricle</th><th>Right Ventricle</th></tr></thead><tbody><tr><td>C<sub>p</sub> (m<sup>3</sup> Pa<sup>-1</sup>)</td><td>2 · 10<sup>-9</sup></td><td>2 · 10<sup>-10</sup></td></tr><tr><td>R<sub>p</sub> (Pa s m<sup>-3</sup>)</td><td>2 · 10<sup>7</sup></td><td>1 · 10<sup>7</sup></td></tr><tr><td>C<sub>d</sub> (m<sup>3</sup> Pa<sup>-1</sup>)</td><td>5.05 · 10<sup>-9</sup></td><td>1 · 10<sup>-8</sup></td></tr><tr><td>R<sub>d</sub> (Pa s m<sup>-3</sup>)</td><td>2 · 10<sup>8</sup></td><td>3 · 10<sup>7</sup></td></tr><tr><td>P<sub>vs</sub> (Pa)</td><td>1 · 10<sup>3</sup></td><td>1 · 10<sup>3</sup></td></tr></tbody></table>

<table><caption>TABLE 4 Physical parameters of the fluid model.</caption><thead><tr><th>Parameter</th><th>Left Ventricle</th></tr></thead><tbody><tr><td>ρ<sub>f</sub> (kg m<sup>-3</sup>)</td><td>1006</td></tr><tr><td>μ<sub>f</sub> (Pa s)</td><td>0.004</td></tr><tr><td>R<sub>closed</sub> (kg m<sup>-2</sup> s<sup>-1</sup>)</td><td>10<sup>7</sup></td></tr><tr><td>ε (m<sup>2</sup> s kg<sup>-1</sup>)</td><td>10<sup>-7</sup></td></tr></tbody></table>

The values of the considered physical parameters for the cardiac and fluid models are reported in Tables 1–4. The ventricles are assumed to be initially in a pre-loaded configuration, corresponding to the beginning of the diastolic phase. This initial configuration is obtained by solving a static version of (13) (i.e., time derivatives are set to zero) with a preload pressure of $P_r = 25$ Pa in the right ventricle and of $P_l = 500$ Pa in the left ventricle. The initial configuration of the fluid domain is then retrieved by lifting the initial solid displacement via $(20)_1$. The fluid is assumed to be initially at rest.
The electromechanical and fluid simulations cover a complete heart beat with a duration of $t = 0.85$ s. We consider a fixed time step of $\tau = 5 \cdot 10^{-4}$ seconds for the solid domain, and $\tau = 10^{-3}$ seconds for the fluid domain. The smaller time step in the solid facilitates the convergence of the Newton iterations during the ventricular contraction phase.

<figure><img src="image_8.png" /><figcaption>FIGURE 9 Deformed left heart cavity configurations and snapshots of the endocardial displacement magnitude at three different time instants.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIGURE 10 Long axis snapshots of the pressure and velocity fields at three different time instants.</figcaption></figure>

Figure 9 provides some snapshots of the endocardial displacement magnitude. During the filling phase (Figure 9(a)), the fluid domain primarily displaces laterally. In contrast, during the ejection phase (Figure 9(b)), a significant displacement is observed
<figure><img src="image_10.png" /><figcaption>FIGURE 11 Simulated Wiggers diagrams.</figcaption></figure>

near the apex of the left ventricle, representing its contraction. The isovolumic relaxation phase (Figure 9(c)) shows a slight decrease in the displacements due to the isovolumetric constraint.

In Figure 10 we have reported the fluid pressure and velocity fields obtained with the velocity constrained model on a long axis cut plane of the left heart at three different time instants. Similar results are obtained with the RIS model, which are omitted here for the sake of conciseness. At the filling phase (Figure 10(a)), the emergence of a single vortex that rotates clockwise is clearly visible, which is in agreement with available physiological observations (see, e.g.,$^{63,64}$). At this stage, the pressure in the left atrium is slightly larger than in the left ventricle, which induces blood flow from the left atrium to the left ventricle. On the other hand, the pressure in the aorta is higher than in the left ventricle, indicating that the aortic valve is closed, and there is no flux in that region. During the ejection phase (Figure 10(b)), there is a flux from the left ventricle to the aorta, with a slightly higher pressure in the former region. At this stage, the pressure in the left atrium is lower than in the left ventricle, indicating that the mitral valve is closed, and there is no flux in this region. Finally, during the iso-volumetric relaxation phase, the pressure inside the left ventricle is higher than in the left atrium but lower than in the aorta, indicating that both the mitral and aortic valves are closed.

Figure 11 presents the numerical Wiggers diagrams obtained with the RIS model and the velocity constrained reduced valve models. For comparison purposes, the 0D left ventricular pressure provided by the cardiac mechanics model (13) is also shown. One can clearly see that both models produce strong spurious oscillations in the ventricular pressure, which are not observed in cardiac mechanics simulation. The peak pressure occurs at $t = 225$ and $t = 226$ ms with values of $1.7 \cdot 10^7$ Pa and 44101.3 Pa for the RIS and the velocity constrained model, respectively. Lower oscillations are thus observed in the latter. Moreover, the RIS model exhibits a more pronounced delay in the onset of iso-volumetric relaxation. However, the oscillations persist irrespective of the valve model considered. These results suggest that the presence of spurious oscillations may be more related to the type of model coupling approach between the fluid and solid than to the reduced model considered for the valves.

## 4.2 | Fluid-structure interaction

The purpose of this section is twofold. We first consider a representative simplified two-dimensional setting in which numerical evidence on the accuracy of Algorithm 2 is provided (Section 4.2.1). The impact of the free Robin coefficient $\alpha$ and of the number of correction iterations is addressed in this example. Finally, in the context of the coupled FSI problem (13), (14) and (22), we delve into the capabilities of Algorithm 2 to provide physiological simulations of the left heart hemodynamics (Section 4.2.2), notably when compared with kinematical uncoupling results of the previous section.
## 4.2.1 | Academic test case

We consider a simplified two-dimensional geometrical setting which retains the main fundamental ingredients of the cardiac cycle phases in the left ventricle. The fluid domain $\Omega_f$ is assumed to be partitioned in three disjoint sub-domains $\Omega_{f,1} \stackrel{\text{def}}{=} [0, 0.025] \times [0, 0.03]$, $\Omega_{f,2} \stackrel{\text{def}}{=} [0.025, 0.115] \times [0, 0.03]$ and $\Omega_{f,3} \stackrel{\text{def}}{=} [0.115, 0.14] \times [0, 0.03]$ (dimensions in m), which respectively correspond to the atrium, the left ventricle and the aorta. The mitral and aortic valves are described by the immersed interfaces $\Gamma_{12} \stackrel{\text{def}}{=} \partial\Omega_{f,1} \cap \partial\Omega_{f,2}$ and $\Gamma_{23} \stackrel{\text{def}}{=} \partial\Omega_{f,2} \cap \partial\Omega_{f,3}$, respectively. The boundaries of the three fluid domains are partitioned as follows: $\partial\Omega_{f,1} = \Gamma_{f,1}^{\text{in}} \cup \Gamma_{f,1}^{\text{t}} \cup \Gamma_{f,1}^{\text{b}} \cup \Gamma_{12}$, $\partial\Omega_{f,2} = \Gamma_{1} \cup \Gamma_{f,2}^{\text{b}} \cup \Gamma_{12} \cup \Gamma_{23}$ and $\partial\Omega_{f,3} = \Gamma_{f,3}^{\text{out}} \cup \Gamma_{f,3}^{\text{t}} \cup \Gamma_{f,3}^{\text{b}} \cup \Gamma_{12}$, see Figure 12.

<figure><img src="image_11.png" /><figcaption>FIGURE 12 Geometrical configuration.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIGURE 13 Snapshots of the fluid pressure and velocity fields obtained with Algorithm 2 for $\tau = 10^{-4}$, $\alpha = 10^4$ at four different time instants.</figcaption></figure>

Flow in the fluid cavities is described by the Navier-Stokes system (14) with the RIS model (15) for the valves. The dynamics of the myocardial wall $\Gamma_1$ are described by a non-linear Timoshenko beam model (see, e.g.,$^{65}$), with the following time-periodic source term which replicates the effect of the active contraction/relaxation:

$$f_a(t) \stackrel{\text{def}}{=} \begin{cases} 5075 \tanh(100(t - mt_p - 0.1)) + 5075 & \text{if } t \in (mt_p, 0.18 + mt_p), \\ -5075 \tanh(100(t - mt_p - 0.2)) + 5075 & \text{if } t \in (0.18 + mt_p, (m+1)t_p), \end{cases}$$

where $t_p \stackrel{\text{def}}{=} 0.25$ is the length of the period and $m \in \mathbb{N}$ indicates the current period number. As regards the boundary conditions, the solid domain is clamped at its extremities $\partial\Gamma_1$. A non-slip condition is enforced on the fluid upper wall $\Gamma_{f,1}^t \cup \Gamma_{f,3}^t$, a symmetry
condition on the fluid bottom wall $\Gamma_{f,1}^b \cup \Gamma_{f,2}^b \cup \Gamma_{f,3}^b$ and the following conditions are enforced on the inlet and outlet fluid boundaries:

$$ \left\{ \begin{array}{ll} \sigma(u,p)n_f + \frac{\rho_f}{2} |u \cdot n_f|_- u = -100n_f & \text{on } \Gamma_{f,1}^{\text{in}}, \\ \sigma(u,p)n_f + \frac{\rho_f}{2} |u \cdot n_f|_- u = -(10^4 + R \int_{\Gamma_{f,3}^{\text{out}}} u \cdot n_f)n_f & \text{on } \Gamma_{f,3}^{\text{out}}. \end{array} \right. $$

Finally, coupling conditions similar to (22) are imposed on the fluid-solid interface $\Gamma_1$. The fluid and the solid are assumed to be initially at rest. $kg m^{-3}$, as Young's modulus $E = 10^5$ Pa and as Poisson ratio $\nu = 0.5$.

The fluid density is set to $\rho_f = 1006$ $kg m^{-3}$, the dynamic viscosity to $\mu = 0.004$ Pa s and the outlet resistance to $R = 1.5 \cdot 10^4$ Pa s $m^{-3}$. For the solid we take as density $\rho_s = 1200$ The resulting coupled problem is approximated numerically

<figure><img src="image_13.png" /><figcaption>FIGURE 14 Temporal evolution of the area $|\Omega_{f,2}^n|$ obtained with the strongly coupled scheme and Algorithm 2 (3 correction iterations and different $\alpha$).</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>FIGURE 15 Time history of pressure obtained with the strongly coupled scheme and Algorithm 2 for $\alpha = 10^4$ and different time-step lengths.</figcaption></figure>

using Algorithm 2 as in Section 3, with the exception that the Timoshenko beam is discretized in space using linear MITC (Mixed Interpolation of Tensorial Components) beam elements (see, e.g.,$^{65}$), available in the FELiScE library. The fluid mesh is made of 2202 triangles and the solid mesh of 40 segments.

For illustration purposes, we have reported in Figure 13 some snapshots of the fluid pressure and velocity fields obtained with Algorithm 2 for $\tau = 10^{-4}$ s, $\alpha = 10^4$ kg m$^{-2}$ s$^{-2}$ and by performing 3 correction iterations. We recall that each correction
<figure><img src="image_15.png" /><figcaption>FIGURE 16 Time history of the area $|\Omega_{f,2}^n|$ obtained with the strongly coupled scheme and Algorithm 2 for $\alpha = 10^4$ and different time-step lengths.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>FIGURE 17 Pressure-area loops obtained with the strongly coupled scheme and Algorithm 2 for $\alpha = 10^4$ and different time-step lengths.</figcaption></figure>

iteration amounts to perform once more Steps 1 and 3 of Algorithm 2 with updated Robin conditions, which enhances the accuracy of the approximation. In fact, iterating in this way until convergence yields a strongly coupled scheme. During the filling phase (Figure 13(a)), the valve $\Gamma_{12}$ is open while $\Gamma_{23}$ remains closed, allowing fluid to fill the ventricular domain $\Omega_{f,2}$. Subsequently, the pressure in $\Omega_{f,2}$ exceeds that in $\Omega_{f,1}$ but remains lower than in $\Omega_{f,3}$, so that both valves $\Gamma_{12}$ and $\Gamma_{23}$ are closed (see Figure 13(b))). During this phase the pressure increases smoothly with no area change in $\Omega_{f,2}^n$. In contrast, during the ejection phase (Figure 13(c)), the pressure in $\Omega_{f,2}$ exceeds that in $\Omega_{f,3}$, which triggers the opening of the valve $\Gamma_{23}$. Finally, Figure 13(d) corresponds to the situation in which both valves are closed, as the pressure in $\Omega_{f,2}$ is higher than in $\Omega_{f,1}$ but lower than in $\Omega_{f,3}$. Here, the pressure decreases smoothly with the area of $\Omega_{f,2}^n$ remaining constant. In summary, this 2D FSI example reproduces the four main phases of the cardiac cycle in the left ventricle: passive filling, isovolumetric contraction, ejection, and isovolumetric relaxation.

In order to assess the accuracy of Algorithm 2 (for different $\tau$, $\alpha$ and correction iterations), we compare the results obtained with those of the strongly coupled scheme. The time histories of the ventricular area $|\Omega_{f,2}^n|$ obtained with Algorithm 2 with $\alpha \in \{6 \cdot 10^3, 8 \cdot 10^3, 10^4, 1.2 \cdot 10^4\}$, $\tau = 10^{-4}$ and 3 correction iterations are reported in Figure 14. We can clearly see that the accuracy of Algorithm 2 is sensitive to $\alpha$. In particular, small and large values of $\alpha$ can compromise accuracy. These results are also in agreement with the error estimates provided in$^{35}$, which postulate an optimal value of $\alpha$ scaling as $\sqrt{\rho_s E}$.

Figure 15 compares the evolution of the average pressure within the ventricular cavity $\Omega_{f,2}^n$ for different time-step lengths $\tau \in \{10^{-3}, 2 \cdot 10^{-4}, 10^{-4}\}$ and correction iterations in Algorithm 2. The value of the Robin parameter is fixed and set to $\alpha = 10^4$,
based on the results of Figure 14. Once more, for comparison purposes we also show the results obtained with the strongly coupled scheme, for lowest time-step length $\tau = 10^{-4}$. As the number of correction iterations increases or the time-step length $\tau$ decreases, the pressure approximation provided by Algorithm 2 converges to the reference provided by the strongly coupled scheme. For large $\tau$, further correction iterations are necessary for accuracy. However, for moderate and small values of $\tau$, the approximation is close to the reference solution irrespectively of the number of correction iterations. Regardless of the improvement in accuracy when increasing the number of correction iterations or decreasing the time-step lengths, there are no significant differences in the evolution of the pressure with respect to these user-defined parameters.

The temporal evolution of the area of the fluid domain $\Omega_{f,2}^n$ is reported in Figure 16. We can clearly observe that the area is more sensitive than the pressure to the time-step length and the number of correction iterations. Particularly, the evolution of the area is underestimated for large values of $\tau$ or a small number of correction iterations. A salient feature of the results of Figure 16, is that it provides numerical evidence on that Algorithm 2 can accurately capture phases of constant areas (the isovolumetric phases in the cardiac hemodynamics setting) by reducing $\tau$ or by increasing the number of correction iterations. Similar conclusions can be inferred from Figure 17, which provides the pressure-area loops obtained with the strongly coupled scheme and Algorithm 2 for $\alpha = 10^4$ and different time-step lengths. The maximum pressure remains similar for all cases, regardless of the time-step length or the number of correction iterations. Nevertheless, the impact of these user-defined parameters on the area is evident. In general, the smaller the time-step lengths, the fewer correction iterations are necessary to capture the iso-area phases.

### 4.2.2 | Left heart hemodynamics

The purpose of this section is to illustrate numerically the capabilities of Algorithm 2 to provide meaningful left heart hemodynamics simulations, via the approximation of the coupled FSI problem (13), (14) and (22). We also compare the impact on these simulations of the two reduced valve dynamics models presented in Section 2.3. The geometrical setting, initial conditions and computational meshes are the same as in Section 4.1. The considered physical parameters are also those reported in Tables 1–4. We simulate three heartbeats, each with a duration of $t = 0.85$ s. Based on the results of Algorithm 2 in the simplified cardiac example of Section 4.2.1, we set $\tau = 2.5 \cdot 10^{-4}$ s, $\alpha = 10^4$ kg m$^{-2}$ s$^{-2}$ and two correction iterations are performed at each time-step, to accurately capture the isovolumetric phases. In the following, we only present the results of the last heartbeat, that is, for the time interval $[1.7, 2.55]$.

<figure><img src="image_17.png" /><figcaption>FIGURE 18 Pressure and volume time-histories obtained with Algorithm 2, with $\tau = 2.5 \cdot 10^{-4}$, $\alpha = 10^4$ and two correction iterations.</figcaption></figure>

Figure 18 presents the pressure and volume time-histories of the left ventricle and the associated pressure-volume loops obtained with Algorithm 2 and the two reduced valve models. Both reduced models yield similar results, which is also consistent with the results presented in Section 4.1. However, unlike the kinematic uncoupling approach, here no spurious pressure oscillations are seen. This provides numerical evidence on that, as postulated in the introduction, the artificial pressure oscillation observed with the kinematic uncoupling approach are due the lack of compatibility between the kinematical constraint and the
<figure><img src="image_18.png" /><figcaption>**FIGURE 19** Snapshots of the solid displacement magnitude (deformed configuration) at different time instants, obtained with Algorithm 2, with $\tau = 2.5 \cdot 10^{-4}$, $\alpha = 10^4$ and two correction iterations.</figcaption></figure>

status of the reduced valve model. It is also worth noting that the results of Figure 18 show a good reproduction of isovolumic phases.

For illustration purposes, snapshots of the displacement, velocity and pressure fields of the solid and the fluid are presented in Figures 19, 20 and 21, respectively. The cycle starts with the ventricles entering a phase of passive filling. During this phase, the mitral and tricuspid valves are open, allowing blood flow from the atria to the ventricles (see Figure 20(a)). At this stage, the fluid forms a jet immediately after the mitral valve, leading to the generation of a vortex ring that rotates clockwise, located almost at the intermediate section of the left ventricle, as shown in Figure 21(a). This is consistent with biomedical literature (see, e.g.,$^{64,63}$). This flow induces an increase in volume in the left ventricle, and also in the right ventricle due to the lumped parameter fluid coupled with it (see Figure 19(a)). The rate of volume rise remains nearly constant from the mid-to-late diastolic phase. The filling phase continues until the ventricles begin to contract, which marks the end of the diastolic phase. During this phase, the flow from the left atrium to the left ventricle remains relatively constant at a rate of $57.56 \text{ mL s}^{-1}$. The end-diastolic volume (EDV) is $EDV_l = 105.85 \text{ mL}$ for the right ventricle and $EDV_f = 125.72 \text{ mL}$ for the left ventricle, consistently with physiological values reported in literature (see, e.g.,$^{66,67,68}$).

At the end of the diastolic phase, the cardiac cycle transitions into the systolic phase, marked by the ventricular contraction. This contraction leads to an increase in the intraventricular pressure of the right and left ventricles. The increase in the intraventricular pressure induces a backflow from the ventricles to the atrium, leading to the closure of the tricuspid and mitral valves in the right and left ventricle, respectively. This marks the beginning of the isovolumic contraction phase, during which
<figure><img src="image_19.png" /><figcaption>**FIGURE 20** Snapshots of the fluid velocity magnitude at different time instants (deformed configuration) obtained with Algorithm 2, with $\tau = 2.5 \cdot 10^{-4}$, $\alpha = 10^4$ and two correction iterations.</figcaption></figure>

both cavities are closed while the ventricles contract. During this stage, the intraventricular pressure increases but remains lower than the aortic pressure, as seen in Figure 21(b). Figure 20(b) illustrates the velocity field, which becomes more homogeneous due to the contraction of the myocardium.

The isovolumic contraction phase persists until the pulmonary and aortic valves open, as the intraventricular pressures increases to a value higher than the pressures in the pulmonary artery and aorta. This results in the opening of the pulmonary and aortic valves, marking the beginning of the ejection phase. This phase is characterized by the ejection of the blood from the right ventricle into the pulmonary artery and from the left ventricle into the aorta. This phenomenon is evidenced by the substantial deformation in both the right and left ventricles, as it can be appreciated in Figures 19(c) and 19(d). During the first part of this phase, the intraventricular pressures continue to increase reaching maximum values of $p_{l,max} = 22849$ Pa and $p_{r,max} = 4494.78$ Pa at $t = 2.004$ s and $t = 1.99825$ s, respectively, consistently with the physiological ranges reported in the literature (see, e.g.,$^{69,70,71}$). The fluid volume in the left ventricle decreases as the ventricle contracts, reducing gradually the intraventricular pressure and the blood flux from the left ventricle into the aorta. Figures 20(c), 20(c) illustrate the unsteady nature of the flow, with a maximum flow rate of 319.13 mL s$^{-1}$ at $t = 2.00875$ s followed by a constant decrease towards the end of the ejection phase, in correspondence with the small deformations of the ventricles.

The ejection fraction (EF) represents the percentage of blood leaving the heart during contraction, given by $EF \stackrel{\text{def}}{=} SV/EDV$ where SV represents the stroke volume, computed as $SV \stackrel{\text{def}}{=} EDV - ESV$, with EDV representing end-diastolic volume, and
<figure><img src="image_20.png" /><figcaption>**FIGURE 21** Snapshots of the fluid pressure and velocity fields at different time instants (deformed configuration) obtained with Algorithm 2, with $\tau = 2.5 \cdot 10^{-4}$, $\alpha = 10^4$ and two correction iterations.</figcaption></figure>

ESV representing end-systolic volume. The ejection fraction of the left ventricle is $EF = 0.457$, within the physiological range reported in the literature (see, e.g.,$^{66,72}$).

During the last part of the ejection phase, the ventricles start to repolarize. This electrical change triggers the relaxation of the heart muscles, gradually easing the contraction of the ventricles. As the ventricular contraction diminishes, the pressure inside the ventricles drops below the pressure in the pulmonary artery and aorta (see Figure 21(e)). This drop in pressure causes the pulmonary and aortic valves to close, while the tricuspid and mitral valves remain closed. This phenomenon signifies the start of the isovolumetric relaxation phase, during which both ventricular chambers remain closed as the ventricles continue to relax and return to their original configuration, as seen in Figure 19(e). Similar to the isovolumetric contraction stage, the velocity field becomes more homogeneous during this phase (see Figure 20(e)). During the isovolumetric relaxation phase, the intraventricular pressure continues to decrease until its value is lower than the pressure in the atrium. Subsequently, the tricuspid valve and mitral valve open, initiating the filling stage and starting a new heart cycle.

5 | CONCLUSION

In this paper, we have introduced a novel reduced model of valve dynamics which circumvents the limitations of traditional resistive immersed surface (RIS) approaches, where valve laws are given in a pure algorithmic fashion. The basic idea consists
in enforcing unidirectional flow through a penalized inequality constraint on the normal velocity of the fluid. The price to pay (notably with respect to RIS) is a semi-smooth non-linearity in the fluid solver, which can however be efficiently handled with a few Newton iterations. The numerical comparisons between RIS and the velocity constrained models indicate that, though different in essence, they deliver rather similar results. In particular, artificial pressure oscillations are observed in both when combined with a one-way kinematic uncoupling approach.

An alternative left heart hemodynamics modeling approach which combines the above mentioned reduced valve dynamics models with a fluid-structure interaction (FSI) coupling in the left ventricle has also been proposed. In order to mitigate the computational complexity of the coupled system, the interface coupling is solved with a Robin-Robin loosely coupled scheme (Algorithm 2). A salient feature of the proposed method is that an a priori energy analysis (Theorem 2) guarantees the unconditional stability of the fully discrete approximations. The left heart hemodynamics simulations showed that the above mentioned spurious pressure oscillations are removed by the FSI coupling, irrespectively of the reduced valve model considered. The results also showed that the proposed loosely coupled scheme is able to reproduce more physiological isovolumetric phases with a temporal grid of around 4000 time-steps per cardiac cycle (i.e., $τ = 2 \cdot 10^{-4}$s), a couple of correction iterations and a judicious choice of the Robin coefficient $α$ (previously tuned in a simplified cardiac setting). Finally, simulated bio-markers such as ejection fractions, end-diastolic volume and min/max ventricular pressures where within the physiological ranges.

Extensions of this work can explore several directions. Ongoing work focuses on the electromechanical simulation of the atria via a 3D-shell modeling framework, which enables the use of general constitutive laws and facilitates the coupling with the surrounding media (e.g., the left ventricle). A more comprehensive left heart hemodynamics model could hence be envisioned, by incorporating FSI coupling with the atrium. Finally, a natural extension of the present work is also to tackle the open problem of second-order accuracy in the loosely coupled scheme (see<sup>73</sup>, for preliminary investigations in this direction).

## ACKNOWLEDGMENTS

This work was partially funded by the French National Research Agency (ANR), through the SIMR project (ANR-19-CE45-0020). The first author was financially supported by the Vicerrectoria de Investigación, Desarrollo e Innovación (VRIDEI) of the Universidad de Santiago de Chile, and the National Agency for Research and Development (ANID) Doctorado Becas Chile/2020-72210101 of the Government of Chile.

## CONFLICT OF INTEREST STATEMENT

The authors declare no potential conflict of interests.

## DATA AVAILABILITY STATEMENT

Data sharing not applicable to this article as no datasets were generated or analyzed during the current study.

## References

1. Mittal M, Seo J, Vedula V, et al. Computational modeling of cardiac hemodynamics: Current status and future outlook. *Journal of Computational Physics*. 2016;305:1065-1082.

2. Peirlinck M, Costabal Fi, Yao J, et al. Precision medicine in human heart modeling. *Biomechanics and Modeling in Mechanobiology*. 2021;20(3):803-831.

3. Okada JI, Washio T, Sugiura S, Hisada T. Clinical and pharmacological application of multiscale multiphysics heart simulator, UT-Heart. *Korean J. Physiol. Pharmacol.* 2019;23(5):295-303. doi: 10.4196/kjpp.2019.23.5.295

4. Gao H, Feng L, Qi N, Berry C, Griffith B, Luo X. A coupled mitral valve—left ventricle model with fluid-structure interaction. *Medical Engineering & Physics*. 2017;47:128-136.

5. Feng L, Gao H, Griffith B, Niederer S, Luo X. Analysis of a coupled fluid-structure interaction model of the left atrium and mitral valve. *International Journal for Numerical Methods in Biomedical Engineering*. 2019;35(11):e3254.

6. Viola F, Meschini V, Verzicco R. Fluid-Structure-Electrophysiology interaction in the left-heart: A multi-way coupled computational model. *European Journal of Mechanics - B/Fluids*. 2020;79:212-232.

7. Viola F, Del Corso G, De Paulis R, Verzicco R. GPU accelerated digital twins of the human heart open new routes for cardiovascular research. *Scientific Reports*. 2023;13(1):8230.

8. Augustin C, Crozier A, Neic A, et al. Patient-specific modeling of left ventricular electromechanics as a driver for haemodynamic analysis. *EP Europace*. 2016;18(4):121-129.
9. Karabelas E, Gsell M, Augustin C, et al. Towards a Computational Framework for Modeling the Impact of Aortic Coarctations Upon Left Ventricular Load. *Front Physiol.* 2018;9:538.

10. This A, Boilevin-Kayl L, Fernández M, Gerbeau JF. Augmented resistive immersed surfaces valve model for the simulation of cardiac hemodynamics with isovolumetric phases. *International journal for numerical methods in biomedical engineering*. 2020;36(3):e3223.

11. Zingaro A, Fumagalli I, Dede' L, et al. A geometric multiscale model for the numerical simulation of blood flow in the human left heart. *Discrete and Continuous Dynamical Systems - S*. 2022;15(8):2391-2427.

12. Zingaro A, Bucelli M, Fumagalli I, Dede' L, Quarteroni A. Modeling isovolumetric phases in cardiac flows by an Augmented Resistive Immersed Implicit Surface method. *International Journal for Numerical Methods in Biomedical Engineering*. 2023;39(12):e3767.

13. Zingaro A, Bucelli M, R. P. Regazzoni F, Dede' L, Quarteroni A. An electromechanics-driven fluid dynamics model for the simulation of the whole human heart. *Journal of Computational Physics*. 2024;504:112885.

14. Mihalef V, Ionasec R, Sharma P, et al. Patient-specific modelling of whole heart anatomy, dynamics and haemodynamics from four-dimensional cardiac CT images. *Interface Focus*. 2011;1(3):286-296.

15. Mangual J, Kraigher-Krainer E, A. D, et al. Comparative numerical study on left ventricular fluid dynamics after dilated cardiomyopathy. *Journal of Biomechanics*. 2013;46(10):1611-1617.

16. Seo V, Abraham T, Lardo A, Dawoud F, Luo H, Mittal R. Effect of the mitral valve on diastolic flow patterns. *Physics of Fluids*. 2014;26(12):121901.

17. Chnafa C, Mendez S, Nicoud F. Image-based large-eddy simulation in a realistic left heart. *Computers & Fluids*. 2014;94:173-187.

18. This A, Morales H, Bonnefous O, Fernández M, Gerbeau JF. A pipeline for image based intracardiac CFD modeling and application to the evaluation of the PISA method. *Computer Methods in Applied Mechanics and Engineering*. 2020;358:112627.

19. Karabelas E, Longobardi S, Fuchsberger J, et al. Global Sensitivity Analysis of Four Chamber Heart Hemodynamics Using Surrogate Models. *IEEE Transactions on Biomedical Engineering*. 2022;69(10):3216-3223.

20. Bennati L, Giambruno V, Renzi F, et al. Turbulent blood dynamics in the left heart in the presence of mitral regurgitation: a computational study based on multi-series cine-MRI. *Biomechanics and Modeling in Mechanobiology*. 2023;22(6):1829-1846.

21. Bennati L, Vergara C, Giambruno V, et al. An Image-Based Computational Fluid Dynamics Study of Mitral Regurgitation in Presence of Prolapse. *Cardiovascular Engineering and Technology*. 2023;14(3):457-475.

22. Zingaro A, Ahmad Z, Kholmovski E, et al. A comprehensive stroke risk assessment by combining atrial computational fluid dynamics simulations and functional patient data. *Scientific Reports*. 2024;14(1):9515.

23. Bucelli M, Dede' L, Quarteroni A, Vergara C. Partitioned and monolithic algorithms for the numerical solution of cardiac fluid-structure interaction. *Commun. Comput. Phys.* 2022;32(5):1217-1256.

24. Bakir A, Al Abed A, Stevens M, Lovell N, Dokos S. A Multiphysics Biventricular Cardiac Model: Simulations With a Left-Ventricular Assist Device. *Front Physiol.* 2018;9:1259.

25. Astorino M, Hamers J, Shadden S, Gerbeau JF. A robust and efficient valve model based on resistive immersed surfaces. *International Journal for Numerical Methods in Biomedical Engineering*. 2012;28(9):937-959.

26. Fedele M, Faggiano E, Dedè L, Quarteroni A. A patient-specific aortic valve model based on moving resistive immersed implicit surfaces. *Biomechanics and Modeling in Mechanobiology*. 2017;16(5):1779-1803.

27. Zhou G, Saito N. The Navier-Stokes equations under a unilateral boundary condition of Signorini's type. *J. Math. Fluid Mech.* 2016;18(3):481-510.

28. Saito N, Sugitani Y, Zhou G. Unilateral problem for the Stokes equations: the well-posedness and finite element approximation. *Appl. Numer. Math.* 2016;105:124-147.

29. Saito N, Sugitani Y, Zhou G. Energy inequalities and outflow boundary conditions for the Navier-Stokes equations. In: , , Model. Simul. Sci. Eng. Technol. Birkhäuser/Springer, 2016:307-317.

30. Cheng Y, Oertel H, Schenkel T. Fluid-Structure Coupled CFD Simulation of the Left Ventricular Flow During Filling Phase. *Annals of Biomedical Engineering*. 2005;33(5):567-576.

31. Nordsletten D, McCormick M, Kilner P, Hunter P, Kay D, Smith N. Fluid-solid coupling for the investigation of diastolic and systolic human left ventricular function. *International Journal for Numerical Methods in Biomedical Engineering*. 2011;27(7):1017-1039.

32. Santiago A, Aguado-Sierra J, Zavala-Aké M, et al. Fully coupled fluid-electro-mechanical model of the human heart for supercomputers. *International Journal for Numerical Methods in Biomedical Engineering*. 2018;34(12):e3140.

33. Bucelli M, Gabriel M, Quarteroni A, Gigante G, Vergara C. A stable loosely-coupled scheme for cardiac electro-fluid-structure interaction. *J. Comput. Phys.* 2023;490:Paper No. 112326, 24.
34. Bucelli M, Zingaro A, Africa P, Fumagalli I, Dede' L, Quarteroni A. A mathematical model that integrates cardiac electrophysiology, mechanics, and fluid dynamics: Application to the human left heart. *International Journal for Numerical Methods in Biomedical Engineering*. 2023;39(3):e3678.

35. Burman E, Durst R, Fernández M, Guzmán J, Ruz O. Robin-Robin loose coupling for incompressible fluid-structure interaction: non-linear setting and nearly-optimal error analysis. https://inria.hal.science/hal-04258861; 2023.

36. Burman E, Durst R, Guzmán J. Stability and error analysis of a splitting method using Robin-Robin coupling applied to a fluid-structure interaction problem. *Numer. Methods Partial Differential Equations*. 2022;38(5):1396-1406.

37. Burman E, Durst R, Fernández M, Guzmán J. Fully discrete loosely coupled Robin-Robin scheme for incompressible fluid-structure interaction: stability and error analysis. *Numer. Math.* 2022;151(4):807-840.

38. Chapelle D, Le Tallec P, Moireau P, Sorine M. Energy-preserving muscle tissue model: formulation and compatible discretizations. *International Journal for Multiscale Computational Engineering*. 2012;10(2).

39. Kimmig F, Chapelle D, Moireau P. Thermodynamic properties of muscle contraction models and associated discrete-time principles. *Advanced Modeling and Simulation in Engineering Sciences*. 2019;6(1):1-36.

40. Sainte-Marie J, Chapelle D, Cimrman R, Sorine M. Modeling and estimation of the cardiac electromechanical activity. *Computers & structures*. 2006;84(28):1743-1759.

41. Holzapfel G, Ogden R. Constitutive modelling of passive myocardium: a structurally based framework for material characterization. *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences*. 2009;367(1902):3445-3475.

42. Hart J. Normal resting pulse rate ranges. *Journal of Nursing Education and Practice*. 2015;5(8):95-98.

43. Esmaily Moghadam M, Bazilevs Y, Hsia TY, Vignon-Clementel I, Marsden A, Congenital Hearts Alliance (MOCHA) oM. A comparison of outlet boundary treatments for prevention of backflow divergence with relevance to blood flow simulations. *Computational Mechanics*. 2011;48(3):277-291.

44. Chouly F, Hild P, Renard Y. *Finite element approximation of contact and friction in elasticity*. 48 of *Advances in Mechanics and Mathematics*. Birkhäuser/Springer, 2023.

45. Pedrizzetti G, Domenichini F. Left Ventricular Fluid Mechanics: The Long Way from Theoretical Models to Clinical Applications. *Annals of Biomedical Engineering*. 2015;43(1):26-40.

46. Fernández MA, Gerbeau JF. Algorithms for fluid-structure interaction problems. In: ... 1 of MS&A. *Model. Simul. Appl*. Springer, 2009:307-346.

47. Gonzalez O. Exact energy and momentum conserving algorithms for general models in nonlinear elasticity. *Computer Methods in Applied Mechanics and Engineering*. 2000;190(13):1763-1783.

48. Seboldt A, Bukač M. A non-iterative domain decomposition method for the interaction between a fluid and a thick structure. *Numer. Methods Partial Differential Equations*. 2021;37(4):2803-2832.

49. Banks JW, Henshaw WD, Schwendeman DW. An analysis of a new stable partitioned algorithm for FSI problems. Part I: Incompressible flow and elastic solids. *J. Comput. Phys.*. 2014;269:108-137. doi: 10.1016/j.jcp.2014.03.006

50. Fernández M, Mullaert J, Vidrascu M. Generalized Robin-Neumann explicit coupling schemes for incompressible fluid-structure interaction: stability analysis and numerics. *Internat. J. Numer. Methods Engrg.*. 2015;101(3):199-229.

51. Serino DA, Banks JW, Henshaw WD, Schwendeman DW. A Stable Added-Mass Partitioned (AMP) Algorithm for Elastic Solids and Incompressible Flow: Model Problem Analysis. *SIAM Journal on Scientific Computing*. 2019;41(4):A2464-A2484.

52. Gigante G, Vergara C. On the stability of a loosely-coupled scheme based on a Robin interface condition for fluid-structure interaction. *Comput. Math. Appl.*. 2021;96:109-119. doi: 10.1016/j.camwa.2021.05.012

53. Gigante G, Vergara C. On the Choice of Interface Parameters in Robin-Robin Loosely Coupled Schemes for Fluid-Structure Interaction. *Fluids*. 2021;6(6).

54. Temam R. *Navier-Stokes equations*. 2 of *Studies in Mathematics and its Applications*. North-Holland Publishing Co., Amsterdam-New York. revised ed., 1979.

55. Tezduyar T. Stabilized Finite Element Formulations for Incompressible Flow Computations. In: ... 28 of *Advances in Applied Mechanics*. Elsevier, 1991:1-44.

56. Farhat C, Geuzaine P, Grandmont C. The discrete geometric conservation law and the nonlinear stability of ALE schemes for the solution of flow problems on moving grids. *Journal of Computational Physics*. 2001;174(2):669-694.

57. Smaldone S. *Numerical analysis and simulations of coupled problems for the cardiovascular system*. PhD thesis. Université Pierre et Marie Curie, Paris; 2014.

58. Geuzaine C, Remacle JF. Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities. *International Journal for Numerical Methods in Engineering*. 2009;79(11):1309-1331.
59. Chapelle D, Fernández M, Gerbeau JF, Moireau P, Sainte-Marie J, Zemzemi N. Numerical simulation of the electromechanical activity of the heart. In: Ayache N, Delingette H, Sermesant M., eds. *Functional Imaging and Modeling of the Heart*, Springer Berlin. 2009:357-365.

60. Moireau P. *Assimilation de données par filtrage pour les systèmes hyperboliques du second ordre - Applications à la mécanique cardiaque*. PhD thesis. Ecole Polytechnique X, Palaiseau; 2008.

61. FitzHugh R. Impulses and physiological states in theoretical models of nerve membrane. *Biophysical journal*. 1961;1(6):445-466.

62. Nagumo J, Arimoto S, Yoshizawa S. An active pulse transmission line simulating nerve axon. *Proceedings of the IRE*. 1962;50(10):2061-2070.

63. Di Labbio G, Kadem L. Jet collisions and vortex reversal in the human left ventricle. *Journal of biomechanics*. 2018;78:155-160.

64. Domenichini F, Querzoli G, Cenedese A, Pedrizzetti G. Combined experimental and numerical analysis of the flow structure into the left ventricle. *Journal of biomechanics*. 2007;40(9):1988-1994.

65. Bathe K. *Finite Element Procedures*. Prentice Hall, 1996.

66. Clay S, Alfakih K, Radjenovic A, Jones T, Ridgway J. Normal range of human left ventricular volumes and mass using steady state free precession MRI in the radial long axis orientation. *Magnetic Resonance Materials in Physics, Biology and Medicine*. 2006;19:41-45.

67. Gilbert J, Glantz SA. Determinants of left ventricular filling and of the diastolic pressure-volume relation. *Circulation research*. 1989;64(5):827-852.

68. Lorenz C, Walker E, Morgan V, Klein S, Graham T. Normal human right and left ventricular mass, systolic function, and gender differences by cine magnetic resonance imaging. *Journal of Cardiovascular Magnetic Resonance*. 1999;1(1):7-21.

69. Andersen O, Smiseth O, Dokainish H, et al. Estimating Left Ventricular Filling Pressure by Echocardiography. *Journal of the American College of Cardiology*. 2017;69(15):1937-1948.

70. Burns A, La Gerche A, D'hooge J, MacIsaac A, Prior D. Left ventricular strain and strain rate: characterization of the effect of load in human subjects. *European Journal of Echocardiography*. 2010;11(3):283-289.

71. Pinsky MR. The right ventricle: interaction with the pulmonary circulation. *Critical Care*. 2016;20:1-9.

72. Wood PW, Choy JB, Nanda NC, Becher H. Left ventricular ejection fraction and volumes: it depends on the imaging method. *Echocardiography*. 2014;31(1):87-100.

73. Burman E, Durst R, Fernández M, Guzmán J, Liu S. A second-order correction method for loosely coupled discretizations applied to parabolic-parabolic interface problems. *IMA J. Numer. Anal.* 2024. DOI:10.1093/imanum/drae075.