Cardiac Mechanics Benchmark Problems

Overview

The objective of this case study is to test FEBio’s structural mechanics solver against benchmark problems in cardiac mechanics presented in the study by Land and co-workers ^Land-2015Land S, Gurev V, Arens S, et al. Verification of cardiac mechanics software: benchmark problems and solutions for testing active and passive material behaviour. Proc Math Phys Eng Sci. 2015;471(2184):20150641. doi:10.1098/rspa.2015.0641.

Land et al. (2015) describe three problems: (1) Deformation of a beam, (2) inflation of a ventricle, and (3) inflation and active contraction of a ventricle. In this case study, we reproduce Problems 2 and 3.

The models in this case study can be found on the FEBioStudio repository.

Constitutive Law

The constitutive model adopted in that benchmark study is that of Guccione et al. ^{Guccione-1995}Guccione JM, Costa KD, McCulloch AD. Finite element stress analysis of left ventricular mechanics in the beating dog heart. J Biomech. 1995;28(10):1167‐1177. doi:10.1016/0021-9290(94)00174-3, whose strain energy density is $W=\frac{c}{2}\left(e^{Q}-1\right)$, where $Q=b_{f}E_{11}^{2}+b_{t}\left(E_{22}^{2}+E_{33}^{2}+2E_{23}^{2}\right)+2b_{fs}\left(E_{12}^{2}+E_{13}^{2}\right)$, $b_f$, $b_t$ and $b_{fs}$ are unitless material coefficients, and $E_{ij}$ represent rectangular components of the Green-Lagrange strain tensor $\mathbf{E}$.

In FEBio, this constitutive relation can be reproduced using the Fung orthotropic constitutive model ^{Ateshian-2006}Ateshian GA, Costa KD. A frame-invariant formulation of Fung elasticity. J Biomech. 2009;42(6):781‐785. doi:10.1016/j.jbiomech.2009.01.015, with material properties that map into Guccione’s model according to Young’s moduli $E_1=c\, b_f$, $E_2=c\,b_t$, $E_3=c\,b_t$, shear moduli $G_{12}=\frac{c}{2}b_{fs}$, $G_{23}=\frac{c}{2}b_t$, $G_{13}=\frac{c}{2}b_{fs}$, and Poisson ratios $\nu_{12}=\nu_{23}=\nu_{31}=0$. Since the Fung orthotropic model uses an uncoupled formulation to simulate nearly-incompressible behavior, a bulk modulus $k$ also needs to be specified, which should be two to three orders of magnitude larger than the largest Young’s modulus.

Problem 2

For Problem 2, the material properties of the Fung orthotropic material were set to $E_{1}=E_{2}=E_{3}=10$, $G_{12}=G_{23}=G_{31}=5$, $\nu_{12}=\nu_{23}=\nu_{31}=0$, $c=10$ and $k=1000$. To better enforce the incompressibility constraint, Lagrangian augmentation was turned on, and set to converge with a tolerance of 0.01.

Problem 3

For Problem 3, the material properties of the FEBio model were set to $E_{1}=16$, $E_{2}=E_{3}=4$, $G_{12}=G_{31}=4$, $G_{23}=2$, $\nu_{12}=\nu_{23}=\nu_{31}=0$, $c=2$ and $k=1600$.

Active contraction was enforced in FEBio using the uncoupled prescribed uniaxial active contraction material model, lumped with the Fung orthotropic model into an uncoupled solid mixture. The contractile stress was ramped up linearly with time, from 0 to 60.

To better enforce the incompressibility constraint, Lagrangian augmentation was turned on within the uncoupled solid mixture material, and set to converge with a tolerance of 0.01.

FEBio Model

Geometry and Mesh

The details of the geometry can be found in the Land et al. (2015) paper. For implementation into FEBio, a custom C++ code was written to generate a meshed geometry with user-prescribed numbers of elements in the circumferential and axial directions, and through the wall thickness. For the FEBio simulations reported in this case study, 36 elements were used along the circumferential directions, 40 elements along the axial direction, and 10 elements through the wall thickness. These values produced 14040 8-node hexahedral elements and 360 6-node pentahedral elements (Figure 1).

For Problem 3, active contraction was prescribed along fibers whose angle $\alpha$ ranged from -90° at the epicardial surface to +90° at the endocardial surface. In FEBio, the C++ custom code also generated local element material axes such that the first axis in this orthonormal basis was aligned with the local fiber direction (Figure 2).

Boundary Conditions

As described in the Land et al. (2015) paper, the base plane was fixed in all directions. A pressure $p$ was applied to the endocardial surface using FEBio’s pressure load. For Problem 2 $p$ was set to 10, whereas Problem 3 used $p=15$. To achieve best numerical convergence, FEBio’s settings for this pressure used the nonlinear, non-symmetric version of this surface load implementation.

The FEBio model files can be found in Problem 2 and Problem 3.

Results

Problem 2

The FEBio model for Problem 2 had 46365 degrees of freedom. The analysis ran in 91 s on a desktop computer with 18 cores. The mesh in its reference and final configurations is shown in Figure 3. Near-incompressibility was enforced effectively, as assessed by the value of $J=\det\mathbf{F}$ (the determinant of the deformation gradient). For this analysis, it was found that $0.999\le J\le1.002$.

Land et al. (2015) reported the apex location at the end of the ventricle inflation analysis, for the various software packages tested in their benchmark study in their Figure 6, which we reproduce here with the inclusion of FEBio results.

The deformation of a line located in the middle of the ventricular wall was reported in their Figure 7, also reproduced here with the inclusion of FEBio results.

Problem 3

The FEBio analysis for Problem 3 also had 46365 degrees of freedom; it ran in 90 s on a desktop computer with 18 cores. Near-incompressibility was enforced effectively nearly everywhere ($0.997\lesssim J\lesssim1.003$) except for a few elements at the apex, where $J$ dropped to $0.986$ in the worst case. For this problem, the combination of inflation and anisotropic active contraction produced a twisting motion which could be observed using a cutaway view as shown in Figure 4.

The position of the apex was presented in Figure 9 of Land et al. (2015), which is reproduced here with the inclusion of FEBio results.

The deformation of a line located in the middle of the ventricular wall was reported in their Figure 10 and Figure 11, also reproduced here with the inclusion of FEBio results.

Discussion

The objective of this case study was to compare the performance of FEBio against two cardiac mechanics benchmark problems presented in the study by Land et al. (2015).

Results for the ventricle inflation analysis in Problem 2 show that FEBio predicts the apex location within the cluster of results that are most consistent with each other (Figure 6). Similarly, the deformation of a line at the center of the ventricular wall agrees with the great majority of comparison cases presented in Land et al. (2015) (Figure 7).

Results for the ventricle inflation, superposed with active anisotropic contraction (Problem 3), show similarly good agreement for the apex location (Figure 9) and the deformation of a line in the ventricle wall (Figure 10 and Figure 11). In all these cases, the FEBio results are in good agreement with the most clustered results of other open-source software packages examined in the study by Land et al. (2015).

In summary, this case study demonstrates that FEBio’s implementation of finite deformation hyperelasticity and active contraction are well-suited for the examination of problems in cardiac mechanics.

Footnotes

• Land-2015
  Land S, Gurev V, Arens S, et al. Verification of cardiac mechanics software: benchmark problems and solutions for testing active and passive material behaviour. Proc Math Phys Eng Sci. 2015;471(2184):20150641. doi:10.1098/rspa.2015.0641
• Guccione-1995
  Guccione JM, Costa KD, McCulloch AD. Finite element stress analysis of left ventricular mechanics in the beating dog heart. J Biomech. 1995;28(10):1167‐1177. doi:10.1016/0021-9290(94)00174-3
• Ateshian-2006
  Ateshian GA, Costa KD. A frame-invariant formulation of Fung elasticity. J Biomech. 2009;42(6):781‐785. doi:10.1016/j.jbiomech.2009.01.015