Cardiac Mechanics Modeling using svFSI

Modeling cardiac mechanics is a challenging task as the heart muscle is a complex fibrous structure that is predominantly incompressible [1], and undergoes large deformation during a cardiac cycle. svFSI is designed specifically to tackles all these challenges. It has a robust implementation of nonlinear FEM for large deformation problems and several advanced fiber-reinforced hyperelastic material models. In the following sections, we will demonstrate svFSI‘s capability in modeling cardiac mechanics using the benchmark example from Land et al. [2], i.e., the passive inflation of an idealized left ventricle. The complete example can be downloaded from here and you can find the accompanying YouTube tutorial here. We will also briefly explain the theory and implementation of nonlinear continuum mechanics in svFSI.

Mesh Generation

Mesh generation for cardiac mechanics modeling are quite standard and can be achieved following the general guide here.

On the other hand, the heart wall is a composite of layers (or sheets) of parallel myocytes, which are the predominant fiber types. These fiber and sheet directions enable defining a local orthonormal coordinate system inside the cardiac muscle. This local coordinate system is crucial for using a structurally-based constitutive relation for the cardiac muscle [9]. For 3D problems, the users are required to prescribe fiber and sheet directions in the computational domain for certain constitutive relations. Users are provided an option to define a constant fiber direction through the following input directives,

# Fiber direction
Fiber direction: (1.0, 0.0, 0.0)

# Sheet direction
Fiber direction: (0.0, 1.0, 0.0)

svFSI also supports specifying distributed fiber and sheet direction generated by rule-based methods [3] through the following input commands,

# Fiber direction
Fiber direction file path: fibersLong.vtu

# Sheet direction
Fiber direction file path: fibersSheet.vtu

In the latter approach, the vtu file should have the local coordinate system stored as a vector with the name “FIB_DIR” at the centroid of each element. An example is provided here.

Input File

svFSI requires a plain-text input file to specify simulation parameters. An overview of the syntax could be found here. The SimVascular GUI currently supports limited input configurations. To access more advanced functions of svFSI, users are recommended to create their own input file by modifying existing templates. Below is a template of the input file for modeling passive inflation of a left ventricle,

# File svFSI.inp
#----------------------------------------------------------------
# General simulation parameters

Continue previous simulation: f
Number of spatial dimensions: 3
Number of time steps: 100
Time step size: 0.01
Spectral radius of infinite time step: 0.50
Searched file name to trigger stop: STOP_SIM

Save results to VTK format: 1
Name prefix of saved VTK files: result
Increment in saving VTK files: 1
Start saving after time step: 1

Increment in saving restart files: 100
Convert BIN to VTK format: 0

Verbose: 1
Warning: 1
Debug: 0

#----------------------------------------------------------------
# Mesh data including volume (vtu) and surface (vtp) meshes,
# domains and fiber distributions

Add mesh: msh {
   Mesh file path:    <path to mesh-complete folder>/mesh-complete.mesh.vtu
   Add face: endo {
      Face file path: <path to mesh-complete folder>/mesh-surfaces/endo.vtp
   }
   Add face: epi {
      Face file path: <path to mesh-complete folder>/mesh-surfaces/epi.vtp
   }
   Add face: base {
      Face file path: <path to mesh-complete folder>/mesh-surfaces/base.vtp
   }
   Fiber direction: (1.0, 0.0, 0.0)
   Fiber direction: (0.0, 1.0, 0.0)
}

#----------------------------------------------------------------
# Equations solved
# Here we use mixed formulation, ustruct, with stabilization
# displacement-based formulation can be invoked through struct
Add equation: ustruct {
   # Define min and max number of iterations, and convergence
   # tolerance for the nonlinear solver (Newton method)
   Coupled: 1
   Min iterations: 4
   Max iterations: 10
   Tolerance: 1e-6
   Use Taylor-Hood type basis: f

   # Define constitutive model and its parameters
   Constitutive model: Guccione {
      C:   1.0e4
      bf:  1.0
      bt:  1.0
      bfs: 1.0
   }

   # Define a small density value for quasi-steady (static)
   # simulation
   Density: 1e-3

   # Define Poisson ratio
   Poisson ratio: 0.5

#==================================================================================
#  This block is for struct; remove if ustruct is used
   # Define elasticity modulus for the volumetric part of the
   # strain energy function
   Elasticity modulus: 1.0
   # Penalty method to enforce incompressibility
   Dilational penalty model: ST91
   # Use a penalty parameter, if different from bulk modulus
   Penalty parameter: 1.0E6
#==================================================================================

#==================================================================================
#  This block is for ustruct; remove if struct is used
   # Define elasticity modulus to calculate tauM/tauC
   Elasticity modulus: 2.0e5
   # Stabilization paramters
   Momentum stabilization coefficient: 1e-3
   Continuity stabilization coefficient: 1e-3
#==================================================================================

   # Define variables for output
   Output: Spatial {
      Displacement: t
      Velocity: t
      Jacobian: t
   }

   # Linear solver parameter
   # Jacobi preconditioner (FSILS) is the default.
   # For ustruct, trilinos-ilu preconditioner is recommended.
   LS type: GMRES
   {
      Preconditioner:      FSILS
      Tolerance:           1D-6
      Max iterations:      1000
      Krylov space dimension: 50
   }

   # Apply zero displacement BC at the base
   Add BC: base {
      Type: Dir
      Value: 0.0
      Impose on state variable integral: t
   }

   # Add pressure load at the endo surface
   # Here we define a ramp function to linearly
   # increase pressure load from zero to the
   # desired value to increase solver stability.
   # Follower pressure load is set as the direction
   # of the applied load follows deformation.
   Add BC: endo {
      Type: Neu
      Time dependence: Unsteady
      Temporal values file path: <path to load.dat file>
      Ramp function: t
      Follower pressure load: t
   }
}

The applied load on the endocardial surface is provided in a file, load.dat, defining the ramp function. In this file, the first line specifies that there are two data points and the value will change linearly. The second and third line specify the time and the value at that time. Note that for a ramp function, the code expects data at two time points only. Should the simulation go beyond the last time stamp (t=1.0 in the current example), a constant value equal to the last extrema (1.0e4) will be applied.

2    1
0.0  0.0
1.0  1.0e4
# content of load.dat

Theory and Implementation

This section briefly explains the theory and implementation of the nonlinear solid dynamics solver in svFSI. Two types of formulations are provided in svFSI for modeling solid mechanics: STRUCT and uSTRUCT. STRUCT uses a pure displacement-based formulation of the balance laws (mass and momentum conservation principles), while uSTRUCT uses a mixed (velocity-pressure) formulation of the governing equations. The latter approach uses either a stabilized equal-order discretization for velocity and pressure function spaces, or the inf-sup-conforming Taylor-Hood-type finite element discretization.

Kinematics

During a cardiac cycle, the heart undergoes large deformations which can no longer be described by linear elasticity. Since the domain is continuously deforming, we introduce the concepts of a reference configuration (denoted $\Omega_{0}$) and a current configuration (denoted $\Omega_{t}$). The reference configuration is fixed for all times, and often refers to the initial geometry when $t=0$. We will use the vector $\mathbf{X}$ to denote the physical coordinates of the geometry in the reference configuration, and define the following relation:

$$\mathbf{x}(\mathbf{X}, t)=\mathbf{X}+\mathbf{u}(\mathbf{X}, t)$$

where $\mathbf{x}$ describes the physical coordinates of the geometry in the current configuration and $\mathbf{u}$ is the time-varying displacement vector field on $\Omega_{0}$ which acts as a mapping between the reference and current configurations such that $\mathbf{u}: \Omega_{0} \Rightarrow \Omega_{t} .$ We also define the following relationships and tensors, which are standard for describing nonlinear structural mechanics:

$$ \ddot{\mathbf{u}}=\frac{d^{2} \mathbf{u}}{d t^{2}}, \quad \mathbf{F}=\frac{\partial \mathbf{x}}{\partial \mathbf{X}}, \quad J=\operatorname{det}(\mathbf{F}) $$

$$ \mathbf{E}=\frac{1}{2}(\mathbf{C}-\mathbf{I}), ~~ \mathbf{C}=\mathbf{F}^{T} \mathbf{F} $$

where $\ddot{\mathbf{u}}$ refers to the structural acceleration and the time derivative operator, $d^{2} / d t^{2}$, is applied on the reference configuration. $\mathbf{F}$ denotes the deformation gradient tensor and $J$ is the the determinant of $\mathbf{F}$ tensor that denotes the Jacobian of the transformation. $\mathbf{C}$ is the Cauchy-Green deformation tensor while $\mathbf{E}$ is the Green-Lagrange strain tensor.

STRUCT: displacement-based nonlinear elastodynamics

The STRUCT equation solves equations for nonlinear structural dynamics using finite element formulation. We start with the function spaces and weak form. We require that the trial and weighting functions satisfy their respective properties on the current domain. The strong form of the momentum balance is

$$ \rho\ddot{\mathbf{u}} = \nabla_{x}\cdot{\sigma} + \rho\mathbf{f_b}, $$ $$ \sigma\cdot \mathbf{n} = \mathbf{h},~\mathrm{on}~\left(\Gamma_{t}\right)_{h}. $$

where $\sigma$ is the stress term in the current configuration.

The resulting weak form of the structural dynamics equations is

$$ \int_{\Omega_{t}} \mathbf{w} \cdot \rho \ddot{\mathbf{u}} d \Omega+\int_{\Omega_{0}} \nabla_{X} \mathbf{w}:(\mathbf{F} \mathbf{S}) d \Omega-\int_{\Omega_{t}} \mathbf{w} \cdot \rho \mathbf{f_b} d \Omega-\int_{\left(\Gamma_{t}\right)_{h}} \mathbf{w} \cdot \mathbf{h} d \Gamma_{h}=0 $$

The acceleration term (i.e. $\int_{\Omega_{t}} \mathbf{w} \cdot \rho \ddot{\mathbf{u}} d \Omega$), body forcing term (i.e. $\int_{\Omega_{t}} \mathbf{w} \cdot \mathbf{f_b} d \Omega$), and the natural boundary condition term (i.e. $\int_{\left(\Gamma_{t}\right)_{h}} \mathbf{w} \cdot \mathbf{h} d \Gamma_{h}$) are all evaluated in the current configuration. These terms are commonly referred to as external work done on the structure by body forces and surface tractions. The remaining stress term in the equation is often referred to as internal work done on the structure by internal stresses, which we will treat differently here. We rewrite this in the reference configuration in terms of the deformation gradient and the second Piola-Kirchhoff stress tensor, which is commonly denoted as $\mathbf{S}$.

svFSI uses a splitting approach where the strain energy and the resulting second Piola-Kirchhoff stress tensor, $\mathbf{S}$, are decomposed into deviatoric (or isochoric, $\mathbf{S}_{iso}$) and dilational (or volumetric, $\mathbf{S}_{vol}$) components. The specific form of $\mathbf{S}$ will depend on the chosen constitutive model (isochoric and volumetric). More information on these Material models can be found in the literature [1]. It is noted that the symbol $E$ denotes the elastic modulus of the structure and is not to be confused with $\mathbf{E}$, which denotes the Green-Lagrange strain tensor, and $\nu$ represents Poisson’s ratio. Some key material parameters can then be defined as,

$$ \lambda=\frac{E v}{(1+v)(1-2 v)}, \quad \mu=\frac{E}{2(1+v)}, \quad \kappa=\lambda + \frac{2}{3} \mu $$

where, $\lambda$ and $\mu$ are the Lame’s first and second parameters, respectively, and $\kappa$ is the bulk modulus. The second Piola-Kirchhoff stress tensors for a few standard constitutive models are given as,

$$ \mathbf{S}^{StVK} = 2 \mu \mathbf{E} + \lambda \operatorname{tr}(\mathbf{E}) \mathbf{I}, ~~~ \textrm{St. Venant-Kirchhoff} $$ $$ \mathbf{S}^{mStVK} = \kappa \operatorname{ln}(J) \mathbf{C}^{-1} + \mu(\mathbf{C} - \mathbf{I}), ~~~ \textrm{Modified St. Venant-Kirchhoff} $$ $$ \mathbf{S_{iso}}^{nHK} = \mu J^{-2/3} \left(\mathbf{I} - \frac{1}{3} \operatorname{tr}(\mathbf{C}) \mathbf{C}^{-1} \right), ~~ \textrm{Neo-Hookean} $$

where $\mathbf{I}$ is the identity matrix. For the Neo-Hookean and other hyperelastic constitutive models, the $\mathbf{S}$ tensor is computed as $\mathbf{S} = \mathbf{S_{iso}} + \mathbf{S_{vol}}$, where $\mathbf{S_{vol}} = p J \mathbf{C}^{-1}$, and $p$ is the hydrostatic pressure computed based on the chosen dilational strain-energy function. See the section on Material models and the corresponding references for the available dilational penalty models in svFSI .

uSTRUCT: mixed formulation nonlinear elastodynamics

svFSI allows solving nonlinear elastodynamics using a mixed formulation where the structure’s velocity and pressure are the unknown degrees of freedom [4]. Two variants are available within this feature: (a) first, a novel variational multiscale stabilized (VMS) formulation that allows equal-order interpolation of velocity and pressure bases using a unified framework [4]; (b) second, using the classical inf-sup stable Taylor-Hood type elements where the velocity basis is derived from a function space that is one order higher relative to the pressure basis. In the displacement-based formulation, a hyperelastic material model assumes the existence of a Helmholtz free energy potential. However, uSTRUCT postulates hyperelasticity using Gibbs free energy potential [4] and takes the following additive decoupled form as, $$ G(\overline{\mathbf{C}},p,T) = G_{iso}(\overline{\mathbf{C}},T) + G_{vol}(p,T) $$

where $G_{vol}(p,T)=\int \rho^{-1}dp$, $\overline{\mathbf{C}}=J^{-2/3}\mathbf{C}$, $p$ is the pressure and $T$ is the temperature. It is worth mentioning that the free energy above is the specific free energy, i.e. the energy per unit mass. The free energy per unit volume is $G^R=\rho_0G$, where $\rho_0$ is the density of the reference configuration. The Helmholtz free energy per unit volume can be obtained by a Legendre transformation of $-G^R$ as,

$$ H^R(\overline{\mathbf{C}},J,T)=\sup_p\left(-pJ+G^R(\overline{\mathbf{C}},p,T) \right). $$

and due to the additive splitting of the Gibbs free energy, we have,

$$ H^R(\overline{\mathbf{C}},J,T)=G_{iso}^R(\overline{\mathbf{C}},T) + \sup_p\left(-pJ+G_{vol}^R(p,T) \right). $$

It is noted that Gibbs free energy naturally introduces pressure into the stress term. The governing equations for the motion of a continuum body in the current configuration are,

$$ \frac{\mathrm{d} \mathbf{u}}{\mathrm{d} t} - \mathbf{v} = \mathbf{0} $$ $$ \beta(p) \frac{\mathrm{d} p}{\mathrm{d} t} + \nabla_x \cdot \mathbf{v} = 0 $$ $$ \rho(p) \frac{\mathrm{d} \mathbf{v}}{\mathrm{d} t} - \nabla_x \cdot \mathbf{\sigma_{dev}} + \nabla_x p - \rho(p) \mathbf{f_b} = \mathbf{0}. $$

The above system of equations represent the kinematic relation between displacement and velocity, balance of mass and linear momentum. $\sigma_{dev}$ is the deviatoric Cauchy stress, while $\rho$ and $\beta$ are the density and isothermal compressibility coefficient, respectively, defined as functions of pressure. The expressions for these quantities are given as follows,

$$ \rho(p) = \left( \frac{\mathrm{d} G_{vol}}{\mathrm{d} p} \right)^{-1}, $$ $$ \beta(p) = \frac{1}{\rho} \frac{\mathrm{d} \rho}{\mathrm{d} p} = -\frac{\partial^2 G_{vol}}{\partial p^2} / \frac{\partial G_{vol}}{\partial p}, $$ $$ \mathbf{\sigma_{dev}} = J^{-1} \bar{\mathbf{F}} \left( \mathbb{P}:\bar{\mathbf{S}} \right) \bar{\mathbf{F}}^T, $$ $$ \bar{\mathbf{S}} = 2 \frac{\partial G_{iso}^R}{\partial \bar{\mathbf{C}}} = 2 \frac{\partial (\rho_0 G_{iso})}{\partial \bar{\mathbf{C}}}, $$

where $\mathbb{P} = \mathbb{I} - \frac{1}{3}\mathbf{C} \otimes \mathbf{C}^{-1}$ is the projection tensor, $\bar{\mathbf{F}} = J^{-1/3} \mathbf{F}$ and $\bar{\mathbf{C}} = J^{-2/3} \mathbf{C}$.

This mixed finite element problem is stabilized using variational multiscale method to allow using equal-order interpolating functions for velocity and pressure unknowns, employ linear elements and handle material incompressibility without suffering from locking issues. Defining an appropriate mixed function space, the stabilized weak form can then be written in the current configuration as,

$$ \mathbf{B}_k := \int_{\Omega_x} \mathbf{w}_\mathbf{u} \cdot \left( \frac{\mathrm{d} \mathbf{u}}{\mathrm{d} t} - \mathbf{v} \right) \mathrm{d} \Omega_x = \mathbf{0} $$

$$ \mathbf{B}_p := \int_{\Omega_x} w_p \left( \beta(p) \frac{\mathrm{d} p}{\mathrm{d} t} + \nabla_x \cdot \mathbf{v} \right) \mathrm{d} \Omega_\mathbf{x} \ + \sum_e \int_{\Omega_x^e} \tau_M^e \nabla_x w_p \cdot \left( \rho(p)\frac{\mathrm{d} \mathbf{v}}{\mathrm{d} t} - \nabla_x \cdot \mathbf{\sigma_{dev}} + \nabla_x p - \rho(p)\mathbf{f_b} \right) \mathrm{d} \Omega_x^e = 0 $$

$$ \mathbf{B}_m := \int_{\Omega_x} \left( \mathbf{w}_\mathbf{v} \cdot \rho(p) \frac{\mathrm{d} \mathbf{v}}{\mathrm{d} t} + \nabla_x \mathbf{w}_\mathbf{v} : \mathbf{\sigma_{dev}} - \nabla_x \cdot \mathbf{w}_\mathbf{v} p - \mathbf{w}_\mathbf{v} \cdot \rho(p)\mathbf{f_b} \right) \mathrm{d} \Omega_x \ -\int_{\Gamma_x^h} \mathbf{w}_\mathbf{v} \cdot \mathbf{h} \mathrm{d} \Gamma_x + \sum_e \int_{\Omega_x^e} \tau_C \left(\nabla_x \cdot \mathbf{w}_\mathbf{v} \right) \left( \beta(p) \frac{\mathrm{d} p}{\mathrm{d} t} + \nabla_x \cdot \mathbf{v} \right) \mathrm{d} \Omega_x^e = 0. $$

The stabilization parameters are chosen as,

$$ \mathbf{\tau}_M = \tau_M\mathbf{I}_{nd}, \quad \tau_M = c_m \frac{\Delta x^e}{c\rho}, \quad \tau_C = c_c c\Delta x^e \rho $$

where, $\Delta x^e$ is the diameter of the circumscribing sphere of the tetrahedral element, $c_m$ and $c_c$ are two non-dimensional parameters, and $c$ is the maximum wave speed in the solid body. For compressible materials, $c$ is the bulk wave speed. Assuming isotropic small-strain linear elastic material, the bulk wave speed can be approximated as, $c=\sqrt{ \left( \lambda+2\mu \right) / \rho_0}$, where $\lambda$ and $\mu$ are the Lame’s parameters. For incompressible materials, $c = \sqrt{\frac{\mu}{\rho_0}}$ is the shear wave speed. Further details about the formulation, finite element discretization and time integration could be found in Liu et al. [4].

Material models

Below is the list of available material constitutive models in svFSI :

Volumetric constitutive models for the struct/ustruct equations
Volumetric Model Input Keyword
Quadratic model “quad”, “Quad”, “quadratic”, “Quadratic”
Simo-Taylor91 model[5] “ST91”, “Simo-Taylor91”
Miehe94 model[6] “M94”, “Miehe94”
Isochoric constitutive models for the struct/ustruct equations
Isochoric Model Input Keyword
Saint Venant-Kirchhoff$^\dagger$ “stVK”, “stVenantKirchhoff”
modified Saint Venant-Kirchhoff$^\dagger$ “m-stVK”, “modified-stVK”, “modified-stVenantKirchhoff”
Neo-Hookean model “nHK”, “nHK91”, “neoHookean”, “neoHookeanSimo91”
Mooney-Rivlin model “MR”, “Mooney-Rivlin”
Holzapfel-Gasser-Ogden model [7] “HGO”
Guccione model [8] “Guccione”, “Gucci”
Holzapfel-Ogden model [9] “HO”, “Holzapfel”

$^\dagger$ These models are not available for ustruct.

References

[1] Holzapfel, G. A. (2002). Nonlinear solid mechanics: a continuum approach for engineering science. Wiley.

[2] Land, S., Gurev, V., Arens, S., Augustin, C. M., Baron, L., Blake, R., Bradley, C., Castro, S., Crozier, A., Favino, M., Fastl, T. E., Fritz, T., Gao, H., Gizzi, A., Griffith, B. E., Hurtado, D. E., Krause, R., Luo, X., Nash, M. P., … Niederer, S. A. (2015). Verification of cardiac mechanics software: benchmark problems and solutions for testing active and passive material behaviour. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Science, 471 (2184), 20150641. https://doi.org/10.1098/rspa.2015.0641

[3] Bayer, J. D., Blake, R. C., Plank, G., & Trayanova, N. A. (2012). A Novel Rule-Based Algorithm for Assigning Myocardial Fiber Orientation to Computational Heart Models. Annals of Biomedical Engineering, 40(10), 2243–2254. https://doi.org/10.1007/s10439-012-0593-5

[4] Liu, J., & Marsden, A. L. (2018). A unified continuum and variational multiscale formulation for fluids, solids, and fluid–structure interaction. Computer Methods in Applied Mechanics and Engineering, 337, 549–597. https://doi.org/10.1016/J.CMA.2018.03.045

[5] Simo, J. C., & Taylor, R. L. (1991). Quasi-incompressible finite elasticity in principal stretches. Continuum basis and numerical algorithms. Computer Methods in Applied Mechanics and Engineering, 85(3), 273–310. https://doi.org/10.1016/0045-7825(91)90100-K

[6] Miehe, C. (1994). Aspects of the formulation and finite element implementation of large strain isotropic elasticity. International Journal for Numerical Methods in Engineering, 37(12), 1981–2004. https://doi.org/10.1002/nme.1620371202

[7] Gasser, T. C., Ogden, R. W., & Holzapfel, G. A. (2006). Hyperelastic modelling of arterial layers with distributed collagen fibre orientations. Journal of The Royal Society Interface, 3(6), 15–35. https://doi.org/10.1098/rsif.2005.0073

[8] Guccione, J. M., McCulloch, A. D., & Waldman, L. K. (1991). Passive material properties of intact ventricular myocardium determined from a cylindrical model. Journal of Biomechanical Engineering, 113(February), 42–55. https://doi.org/10.1115/1.2894084

[9] Holzapfel, G. A, & Ogden, R. W. (2009). Constitutive modelling of passive myocardium: a structurally based framework for material characterization. Philosophical Transactions of the Royal Society Series A, 367(1902), 3445–3475. https://doi.org/10.1098/rsta.2009.0091

[10] Nolan, D. R., Gower, A. L., Destrade, M., Ogden, R. W., & McGarry, J. P. (2014). A robust anisotropic hyperelastic formulation for the modelling of soft tissue. Journal of the Mechanical Behavior of Biomedical Materials, 39, 48–60. https://doi.org/10.1016/J.JMBBM.2014.06.016