A displacement-based finite element formulation for incompressible and nearly-incompressible cardiac mechanics.

The Lagrange Multiplier (LM) and penalty methods are commonly used to enforce incompressibility and compressibility in models of cardiac mechanics. In this paper we show how both formulations may be equivalently thought of as a weakly penalized system derived from the statically condensed Perturbed Lagrangian formulation, which may be directly discretized maintaining the simplicity of penalty formulations with the convergence characteristics of LM techniques. A modified Shamanskii-Newton-Raphson scheme is introduced to enhance the nonlinear convergence of the weakly penalized system and, exploiting its equivalence, modifications are developed for the penalty form. Focusing on accuracy, we proceed to study the convergence behavior of these approaches using different interpolation schemes for both a simple test problem and more complex models of cardiac mechanics. Our results illustrate the well-known influence of locking phenomena on the penalty approach (particularly for lower order schemes) and its effect on accuracy for whole-cycle mechanics. Additionally, we verify that direct discretization of the weakly penalized form produces similar convergence behavior to mixed formulations while avoiding the use of an additional variable. Combining a simple structure which allows the solution of computationally challenging problems with good convergence characteristics, the weakly penalized form provides an accurate and efficient alternative to incompressibility and compressibility in cardiac mechanics.

Introduction

The human heart is a remarkably complex organ, translating cellular ATP consumption into the systemic blood flow [1]. Over the last four decades, computational modeling of cardiac mechanics has evolved, incorporating biophysically-based hyperelastic strain energy laws [2-5], anisotropic tissue structure [6-8], patient-specific geometries [9] and cellular activation [10] to effectively simulate the myocardial behavior assuming basic Newtonian physics [11]. Based on tunable parameters [12,13], cardiac models provide a framework for studying and assessing heart function, offering spatiotemporally varying metrics– such as strain, stress, work and power– which are otherwise inaccessible clinically [14,15].

While cardiac modeling is capable of providing quantitative data of clinical relevance, a number of modeling questions remain actively pursued in the community. An issue commonly discussed in cardiac mechanics is the choice of modeling myocardial tissue as an incompressible [16,3,17,4,18,19,5,20,21] or nearly incompressible [14,22-24] material. While this choice is inherently based on tissue behavior which must be determined experimentally, both models continue to be used either to model incompressible/nearly incompressible behavior or, in some cases, for numerical convenience.

A range of relevant numerical schemes have been applied in heart models, one of the most popular being the penalty method [25,14,22,15,23]. An advantage of this approach is its simplified form, requiring only the solution of the tissue displacement. However, when applied in the finite element method (FEM) framework, displacement-based formulations near the incompressible limit exhibit locking leading to sub-optimal convergence rates and poor numerical approximations in classic elastic models [26-30]. Critically, the penalty method lacks monotonic convergence to the incompressible solution as the bulk modulus is increased, making it challenging to employ as an approximate model to an incompressible cardiac material model.

The development of numerical strategies circumventing these issues has been a field of significant research effort in the solid mechanics community. Among others, the B-Bar method introduced by Hughes [31] and its generalization to finite strains [32,33], the reduced or selective integration technique [34,30,35], the augmented Lagrangian method [36,37], have been successfully employed to enforce incompressibility while tackling the numerical difficulties and locking phenomena associated with the penalty formulation. An alternative approach used extensively in solid mechanics, also known to alleviate locking, is the class of multi-field variational principles, which gained popularity with the pioneering work of Herrmann on isotropic linear elasticity [38]. Herrmann’s principle was also extended to orthotropic materials by Taylor [39] and Key [40], to nonlinear formulations [41,42] and elasto-plastic applications [43].

The most common of these mixed formulations is the Lagrange Multiplier (LM) method, a two-field variational approach which has been used widely to enforce incompressibility of the myocardium by introducing a variable to respresent the hydrostatic pressure [16,17,4,19,44]. While the LM method is known to improve numerical convergence [45,46,29,47] and avoid locking phenomena, the use of an additional variable results in increased computational cost and enhanced complexity in the linear algebra involved, due to the indefinite nature of the resulting stiffness matrix [26,45].

The Perturbed Lagrangian (PL) formulation was introduced to address this issue, by augmenting the energy functional of the LM approach with a penalty/compressibility term [48-50]. The PL is a two-field variational approach suitable for the solution of nearly incompressible problems, where pressure and displacement are treated as independent variables. Sussman and Bathe introduced a generalized form of the PL approach, the formulation, which has been used extensively in the computational mechanics literature [51,52,48] and has also been applied in the myocardium [44]. Similarly, the well established three-field Hu-Washizu formulation by Simo et al. [33] extends the PL formulation by introducing pressure and dilatation as independent variables [37,46,51,53]. This approach has also been employed in cardiac mechanics [24] (though this procedure comes with the cost of computing an additional variable). The use of a separate interpolation for the independent variables, allows efficient and accurate approximations, alleviating the numerical difficulties associated with both the penalty and LM methods. The efficiency of these methods was also enhanced with the use of discontinuous interpolation for the pressure and dilatation fields (static condensation) [50,47,33,46] allowing the estimation of these fields on element level and leading to a generalized displacement-only formulation. Further, Bercovier [50] proved that, for Herrmann’s principle, the PL (and its statically condensed equivalent) converges monotonically to the incompressible problem as the bulk modulus is increased. Nevertheless, as suggested by Sussman and Bathe [47], static condensation may exhibit convergence difficulties during the Newton–Raphson procedure.

In this paper, we consider the statically condensed Perturbed Lagrangian formulation of Bercovier [50] and others [48,49], which may be conveniently thought of as a weakly penalized form with an optional choice of projection operator. In this generalized form, with an appropriate choice of the projection operator, we may choose to strengthen or weaken the constraint resulting in the PL, LM or penalty formulations. Using this generalization, we derive an estimate detailing the error convergence of these methods (in a linear setting) and introduce modifications to a Newton-Raphson scheme [54,55] to significantly improve nonlinear convergence properties for standard and weakly penalized formulations (particularly for high bulk modulus). The scheme is further augmented to take advantage of a Shamanskii-type Newton scheme [54,55] boosting computational performance by enabling re-use of the Jacobian matrix (and its inverses or preconditioners) estimated at previous time/loading steps. As this re-use is particularly sensitive to stiffness, we modify the scheme to effectively maintain nonlinear convergence behavior. Further, we examine the direct numerical discretization of the weakly penalized form which may be made efficient through the use of discontinuous projection operators. The weakly penalized form is then compared with the mixed variational formulations (LM, PL), as well as the penalty method, showing that the modified form maintains the convergence characteristics of the mixed variational forms and avoids locking behaviors observed in the penalty method. This comparison is performed on a model left ventricle, which to the best of our knowledge is the first application of this combination of the PL method and static condensation in cardiac mechanics. Further we verify the result proven for linear problems in [50], showing that the error between the weakly penalized formulation and the incompressible solution indeed decreases with a rate inversely proportional to the bulk modulus. As a result, the formulation enables modeling of the myocardium as nearly incompressible or incompressible (with an error proportional to , with k being the bulk modulus).

Below we expand on this approach to illustrate the general minimization problem (Section 2.1) and show how both penalty and LM formulations may be thought of equivalently as weakly penalized constraints in the continuous setting (Section 2.1.1). The basis for locking is then reviewed in Section 2.2, motivating the introduction of the weakly penalized approach. The different convergence behavior of the various schemes is also illustrated through their error estimates at the solution of a linear incompressible problem (Section 2.4). We then introduce modifications to the mechanical system to improve the nonlinear convergence behavior of the weakly penalized scheme for high values of the bulk modulus (Section 2.5.3). Moreover, we modify the SNR scheme for the weakly penalized and penalty formulations to enable better computational efficiency (Section 2.5.4). The numerical convergence of these different methods is then compared, showing optimal convergence and locking phenomena in the various schemes [27,56] for a two-dimensional problem and a cardiac model (Section 3). Finally, the LM, penalty and weakly penalized formulations are compared in terms of accuracy and convergence for whole-cycle cardiac mechanics. Our results suggest that the weakly penalized form provides an accurate and computationally efficient alternative to the LM, PL and penalty methods and can be applied successfully in the numerical implementation of incompressible and nearly incompressible cardiac models (Section 4).

Methods

In this section we show how LM and Penalty formulations can be viewed uniformly through a weakly penalized form (PL) (Section 2.1). Subsequently, we introduce the discretized forms, illustrating the deviation of the two schemes and resulting locking phenomena (Section 2.2). These motivate the use of an alternative discretization strategy, leading to a displacement-only formulation. The solution to this system is then demonstrated and optimized to accommodate the numerical stiffening due to the weakly penalized terms.

Continuous minimization problem

Problems of static (or quasi-static) solid mechanics involve finding the deformation, , of a body defined over a domain as shown in Fig. 1. Here, the body is under the action of a body force, , and some boundary traction .

Fig. 1

The undeformed and deformed body under consideration. Here represents the reference state of the body, its boundaries subject to Dirichlet and traction conditions, respectively.

The solution is a function commonly sought in an appropriate function space (usually or a space smooth enough to ensure existence and uniqueness of the solution of the minimization problem [61,48]), subject to the Dirichlet boundary condition , i.e.The displacement of the body may be obtained by the principle of virtual work, equivalent to the principle of stationary potential energy [45]. Following the principle of stationary potential energy we seek to find a minimizer of the total potential energy functional, , describing the total potential energy of the body under consideration (see Fig. 1). Assuming that the traction and body forces are not functions of the displacement , then under static equilibrium, the total potential energy for a hyperelastic body may be expressed as a sum of the internal and external potential energy as,where represents the strain energy function [45]. According to the principle of stationary potential energy, the body will deform in a way that minimizes its total potential energy . This problem can be expressed as,In the case of incompressibility, the deformation is required to preserve the determinant of the deformation gradient ,In this case, the solution is found [29], which satisfies,We note that as it is not, in general, straightforward to construct the space it is often preferable to seek the solution in the entire space.

Weakly penalized form and the penalty/LM/PL methods

For later comparisons, in this section we introduce a weakly penalized form of the mechanical problem and show its equivalence with both penalty and LM formulations. Here we introduce the projection operator which, for any function , denotes the orthogonal projection onto W, i.e.In this way, we may elect to represent g coarsely or finely by adjusting the selection of the space W (as we will discuss further in the following sections). We may then introduce the weakly penalized total potential energy functional,where an additional penalty term has been added, representing the growth in energy resulting from material compression as is typical for many penalty methods. However, the presence of the projection enables the selective weakening or strengthening of the constraint by allowing it to hold weakly through Eq. (6). Clearly, when for any (for example, when as is the case for the continuous mechanical system [48]), thenand Eq. (7) reduces to the classic total potential energy functional,where the k-dependent term denotes the volumetric penalty term often used in cardiac mechanics [25,23]. The Perturbed Lagrangian formulation may be derived by introducing an additional variable, with,Substituting the orthogonal projection with the added variable, and adding the Galerkin orthogonality condition (Eq. (6)) we arrive at the PL functional [50,47,48,57],

This general purpose formulation has been employed for the solution of nearly incompressible materials [47,50,48,49] and cardiac mechanics [44]. Note that as the previous formulation becomes the classic LM method [26,29].

In the continuous setting, there is a solution which satisfies,for all approaches and all values of k. However, this equivalence is often lost in the discrete setting as different strategies are applied to discretize the function space, , and the orthogonal projection, .

Finite element approximation

In the FEM framework used in the solution of Problems (12)–(14), the domain is subdivided into a mesh of non-overlapping elements [56]. The displacement is then interpolated with functions in , consisting of a set of piecewise polynomials () on the mesh , i.e.where denotes the order of interpolation used for the displacement and is the space of continuous vector functions. Letting be an vector function comprised of the basis functions , the resulting displacement solution is then expressed as the weighted sum,In all approaches, the minimization of the total potential energy occurs over .

The primary point of departure between the penalty and PL formulations comes in the choice of orthogonal projection. In the case of the penalty method, the orthogonal projection remains on the continuous space , leaving the total potential energy unchanged. In contrast, the PL approach given in Eqs. (11) and (14) requires a numerical approximation of the pressure variable, by , it is hence natural to introduce a discrete function space , i.e.This is equivalent to introducing the projection operator which satisfies the Galerkin orthogonality condition (Eq. (6)) on . This change means that projects the incompressibility onto a discrete set of polynomials of degree , effectively relaxing strict satisfaction of the constraint. The weakening of the constraint through the projection operator is similar to reduced/selective integration techniques which also weaken the constraint on compressibility/incompressibility. However, while these techniques have been proven consistent for specific quadrature and element schemes [26], proof is required to ensure each scheme achieves an optimal rate of convergence. The spaces and for the weakly penalized approach are often selected to satisfy the inf–sup condition (where we note ), which ensures uniqueness in the multiplier for all k [58,26,59-61]. Additionally, for appropriately selected spaces – such as with – convergence rates in the energy norm are optimal (i.e. ).

It is well known that these two methods need not be equivalent in the discrete setting, as can also be confirmed by the presence of locking phenomena in penalty applications [26,27]. These facts can also be observed through the dependence of the methods on the penalty parameter k. As , the PL becomes the classic LM method and the approximation spaces reduce to the subsets,for the penalty and LM methods, respectively. These spaces are nested, i.e.illustrating the more restrictive subset of functions over which the minimization problem can be considered. As a consequence, the space is typically a small subset of and can be too restrictive. This occurs as a small violation of the incompressibility constraint can cause a significant increase in the strain energy even though the approximate solution may have a minor degree of error from the true solution.

In contrast, while the LM approach effectively weakens the satisfaction of the constraint, it also has proven optimal convergence rates when and are chosen to satisfy the inf–sup condition [26]. This circumvents over-constraining of the approximation space, but comes at the expense of computing an additional variable.

Discrete weakly penalized form

In the discrete setting, the projection operator introduced in Eq. (6) can be written as,where is a test function in denotes the projection of g on is the mass matrix,and is the weighted function over the test space ,Considering the introduced term in Eq. (7), we may writewhere here . Following from Eq. (19) and noting the requirement that the projection holds for is equivalent to requiring it hold for all (where is the dimension of the discrete space ), can be seen to satisfy the linear system,where is given in Eq. (20)Inverting in Eq. (23) and substituting into Eq. (22), the weakly penalized system (Eq. (7)) may be written in discrete form as,

This form, also used by Bercovier [50] and others [48,49], reduces the system into a single minimization problem on , eliminating the pressure variable. However, the presence of requires the solve of Eq. (23), incuring similar computational cost as computing the PL solution. Considering the discrete weak form and its solution (as we show later), this weakly penalized term requires matrix–matrix products which are (1) expensive, (2) nonlinearly dependent on the solution and (3) generally of a more dense sparsity than the standard penalty system.

These practical issues presented stem from the choice of global orthogonal projection, , which unfortunately complicates the computation. However, the choice of is, generally speaking, arbitrary and ideally should balance the need for accuracy with ease of computation. An alternative approach is to use a local orthogonal projection, , satisfyingwhere and are the mass matrix and the weighted constraint vectors on the element . That is, the local orthogonal projection, , satisfies Eq. (6) on the piecewise discontinuous space,

Using this locally continuous, but globally discontinuous interpolation space, the total potential energy for the body becomes,

This localized projection, also known as static condensation has been employed by Sussman and Bathe [47], Bercovier [50] and Simo et al. [33,37] to enhance the efficiency of the formulation, while avoiding over-constraining of the approximation space and locking phenomena. Indeed, by localizing the projection, computations remain on the element level, reducing the computational cost relative to Eq. (25). Moreover, localization of the penalty term preserves sparsity of the penalty system, significantly reducing sparsity to that resulting from Eq. (25). These practical improvements come at the cost of restricting the approximation space. Again, as we send , the approximation space of the weakly penalized formulation in Eq. (27) is restricted to the space,which we note is,Note that in a practical setting, as k can not be infinite, an augmented Lagrangian iterative scheme can be applied to iteratively increase k. In this way, the weakly penalized form can provide equivalent results to the incompressible LM method.

Though the weakly penalized form in Eq. (27) does not mandate the inf–sup condition, usage of inf–sup stable spaces (such as Nicolaides–Boland [60] or Crouzeix–Raviart [62] elements) with globally discontinuous pressure ensures optimal convergence. However, as we demonstrate, even for some spaces which are not inf–sup stable (for instance ), this weakening of the constraint is sufficient to restore convergence.1

Error estimates for the generalized weakly penalized form

Using the generalized weakly penalized form, Appendix B shows an error estimate in Lemma 1 in the case of a linear elastic model. Here, the projection and bulk modulus in the discrete model are left general, enabling extension to the methods (LM, PL, and Penalty) considered in the paper. From this analysis, we observe in Lemma 1 (when ) that the continuous model () and discrete model satisfy the estimate,(see Appendix B). Here, the estimate shows the approximation depends on three principle terms which relate the error incurred due to the projection as well as the error approximating divergence free and non-divergence free components of the solution. Importantly, the bound depends on the subspaces and which, in general do not exhibit straightforward convergence properties. However, as demonstrated in Corollary 1 of Appendix B, if projects to a discrete space which satisfies the inf–sup condition [50,59,26,61], the estimate can be extended so the infimum is taken over . This factor – which is exploited in the LM and PL formulations (as well as the weakly penalized form suggested above) – enables straightforward application of standard estimates derived from interpolation theory.

However, in the case of the penalty method, remains the continuous -projection, limiting the extension discussed. While some simplification can be made (see Corollary 2), the estimate requires use of which may be prohibitive, limiting the convergence behavior. While more straightforward estimates may be derived (see for example [36]), in these the scaling constant C depends on k.

These results are in agreement with the the numerical results presented later in this paper.

Discrete weak form of the weakly penalized formulation

This section deals with the weak form of the weakly penalized formulation and the modifications introduced in the residual evaluation to enable improved nonlinear convergence and better performance of the SNR scheme. The discrete weak form for the weakly penalized formulation can be obtained by requiring that the directional derivative of the total potential energy functional vanishes in all arbitrary directions at , i.e.with the homogeneous zero Dirichlet subspace of . Following this procedure, the discrete weak form can be written in operator notation as,

where R is the residual function and the operators , and C are defined as,where and are the discrete deformation gradient and second Piola stress tensors, respectively [45]. Here represents the local basis coefficients for on the element and the element matrix denotes the linearized constraint derived from the PL functional, i.e.with . Note that and are identical to those in Eq. (20) and (24) with .

The weak forms for the penalty, LM and PL formulations can be derived similarly as outlined in Appendix A.

Nonlinear solution for the weakly penalized problem

In order to solve the mechanical system introduced in Eq. (30) (as well as the others discussed in the Appendix A), we look to use the global Shamanskii–Newton–Raphson (SNR) method [54]. This method has been shown to be effective for problems in fluid–structure interaction [55], enabling faster computation by re-using the Jacobian matrix over multiple time/load steps. Following the procedure outlined in [55], on the SNR iteration we update each subsequent guess of the solution, , using the iterative formula,Here denotes the basis weights at the iteration, the update vector, and and are the residual and Jacobian, respectively, defined in Eq. (33).The key distinction between the Newton–Raphson method, global Newton–Raphson method and SNR method introduced in Eq. (32) are the parameters and . In classical Newton–Raphson, these parameters play no role in the update process (in this case ). The global Newton–Raphson scheme, however, uses the parameter to scale the descent direction to minimize and ensure that , i.e. that the residual decays (in this case ).

In constrast, the SNR method shown in Algorithm 1 combines the minimization procedure with a re-use scheme, where denotes the use of a Jacobian matrix determined at a previous load step/iteration. The selection of is based on the convergence characteristics governed by and ITER.2 Here governs the degree of residual decrease required to continue using the currently stored Jacobian. ITER may also signal re-computation and enables a cap on the degree of re-use. Finally, the parameters is used to cope with stiff problems which, due to Dirichlet conditions, may see initial increases in the residual computation.

Jacobian and residual construction for the weakly penalized problem

As both penalty and LM formulations have been outlined elsewhere, here we focus on the Jacobian and residual evaluations for the weakly penalized approach. Note that and may be written as,where is the FEM assembly operator, denotes specific elements, the subscript denotes vector, matrices or operators constructed on the element and and are element-level residual contributions stemming from the weakly penalized term, , and all other terms, . Introducing a short-hand notation, i.e. , we can express the element-level Jacobian contributions and as,The element-level matrix denotes those terms resulting from the elasticity stress/boundary contributions3 and is evaluated using central finite differencing (typically , where h is the mesh size). denotes those terms which result from the weakly penalized form. Here we assume that , defined in Eq. (31), is independent of when we linearize the C operator introduced in Section 2.5. This linearization does not seem to impact convergence of the Newton scheme and preserves symmetric positive semi-definiteness of the weakly penalized matrix term. As we see, is comprised of the local element mass matrix and the linearized constraint equation introduced in the previous section. We note that and its inverse are linear and thus may be computed once for the entire simulation. On the other hand, the linearized constraint must be re-computed due to its nonlinear dependence on the solution. However, computing this matrix is quick as it does not require differencing.

Residual modifications for the nonlinear solve of the weakly penalized system

As mentioned previously, the weakly penalized system can be thought of as a generalized formulation which can result in the PL, penalty or statically condensed PL formulations depending on the choice of the projection operator. The equivalence of these methods under the weakly penalized regime, allows us to combine and take advantage of the good characteristics of each method. For instance, the weakly penalized formulation combines the simplified structure of the penalty method with the convergence characteristics of the PL formulation. However, due to the stiffness of the linear system at high values of the bulk modulus, the penalized formulations (classic penalty/weakly penalized) exhibit deteriorated nonlinear convergence. This stands in stark contrast to the PL method which (for inf–sup stable schemes) exhibits fast convergence even for high bulk modulus. However, we observe that, when the choice of provides equivalence with the discrete PL method, poor nonlinear convergence is observed though, in principle, the convergence should be similar. Examining the update formulae for both weakly penalized and PL approaches (see Appendix C), we observe that deteriorated convergence stems from: (1) initial residual amplification, and (2) the amplification of the residual.

The first factor, mentioned in Section 2.5.1, results from non-monotonicity in the residual. This manifests particularly early during the nonlinear solve, where the initial residual becomes amplified after the first iteration. This is particularly evident with Dirichlet conditions on a stiff material, where the norm of the boundary displacement (for example) may be much smaller than the updated residual due to stiffness in the material. This issue which is not observed in the PL solution, is circumvented in the Newton–Raphson procedure and the SNR procedure outlined in Algorithm 1 by enabling an initial amplification of the residual and relaxing strict requirements on monotonicity.

Less trivial issue to address is the amplification of the residual resulting from large k, which can lead to poor convergence or stalling in the iterative solve. This increased residual, due to strong dependence on the weakly penalized term, often does not imply divergence but rather results from a k-dependent scaling of the projection problem and its nonlinear dependence on as shown in Appendix C. Based on Eq. (C.7) (and choosing ), we may re-write the weakly penalized residual contribution in terms of the linearized guess and a remainder,Here the first term approximates the current linearized guess of the projection, while the second examines how well the current guess satisfies the projection problem. In other words, the first term denotes the hydrostatic contribution to momentum, while the second represents the error remaining in the projection problem amplified by the bulk modulus k. For large k, the later term can become disproportionately scaled, making nonlinear convergence more challenging. Moreover, this issue is avoided in the PL formulation where the hydrostatic constraint is not scaled by k as can be seen in the residual terms of Eq. (C.3).

To circumvent this, we modify the Newton–Raphson scheme to measure convergence of instead of . Clearly, as ; however, measuring convergence of avoids issues due to high bulk modulus. Further, extending this form to the Penalty formulation, we must select the projection . As, in this case, the rank of is no longer finite dimensional, we may instead write the modified Penalty form in its equivalent integral form, i.e.

Similarly to the modifications introduced in the weakly penalized approach, in the modified penalty formulation we measure the convergence of instead of . This avoids the amplification of the error in the second term of the residual, allowing better nonlinear convergence of the scheme.

Residual modifications for the SNR solve of the weakly penalized system

While the Shamanskii–Newton–Raphson scheme can significantly enhance performance of the PL scheme, acceleration in the SNR approach for penalized methods is minimal or even worse than standard Global Newton–Raphson. This deterioration in performance is due, predominantly, to the stiffness of the system for high k and the inevitable inaccuracies introduced in the descent direction by the re-used Jacobian. In contrast, this deterioration in performance is not observed in both PL and LM formulations, which can take significant advantage of matrix re-use (as we will show). Once more, by examining the equivalence between weakly penalized and PL formulations (see Appendix C), we observe this deterioration may be circumvented using the modified form,Note that, as , the difference in the constraint residual and as a result,which represents the standard residual resulting from Eq. (30). That is, as we converge, the modified residual converges to that given by evaluating C. Further, we note that if , the first term in the definition drops away, leaving us with the expected residual.

The modifications we introduced to the residual of the weakly penalized formulation can also be applied in the penalty method to improve its nonlinear convergence behavior by once more selecting the projection operator and using the equivalence of the weakly penalized/penalty forms. Expressing these modifications in integral form, we may write the modified Penalty form with residuals,With the introduced modifications in the residual monitored for convergence, the modified penalty method is able to significantly exploit matrix re-use, and substantially improve its computational efficiency.

Results

In this section we study the convergence behavior of the LM, PL, penalty and weakly penalized formulations outlined in Section 2, in the solution of incompressible mechanics problems, for varying values of k. Furthermore, the PL, penalty and weakly penalized approaches are compared on nearly-incompressible solid mechanics problems, assuming various values for the bulk modulus.

Mechanical tests

Elongation of a two-dimensional square domain

The convergence behavior of the LM, penalty and weakly penalized methods was compared on the simple case of the stretch of a square domain (Fig. 2(a)). The body was assumed to be made of a Neo-Hookean material, described by the deviatoric strain energy/Second-Piola stress tensor [45],where the material parameter is analogous to the shear modulus of linear elasticity, is the right Cauchy–Green deformation tensor, is the unity tensor, and are the first and third invariants of , and d is the dimension of the domain (in our case ).

Fig. 2

Discretization of the two-dimensional square domain: (a) A Neohookean material in a square domain () under no slip (bottom edge), no traction (side walls), and vertical displacement of (top edge). The shear modulus of was used. (b) Number of elements and degrees of freedom (DOFS) in each discretization for the () LM/PL and () weakly penalized/penalty methods.

The domain was discretized using six different meshes of inf–sup stable Taylor–Hood quadrilateral elements [26]. For the weakly penalized formulation, a quadratic interpolation was used for the displacement field and a discontinuous linear interpolation was used for the pressure. The actual solution was approximated using the LM and PL solution (for the incompressible and nearly-incompressible comparison respectively) on a finer mesh (mesh7), with a cubic interpolation for the displacement and a quadratic interpolation for the pressure. Fig. 2(b) presents the number of elements in all seven meshes, and the corresponding degrees of freedom when the LM, PL, penalty and weakly penalized methods were used.

Cardiac mechanics in the left ventricle

The three methods were tested on a model of the passive inflation of a left ventricle under diastolic loading conditions. The left ventricle (LV) was modeled as a thick-walled truncated ellipsoid (Fig. 3(a)). A standard generic heterogeneous fiber field was used to represent the structure of the tissue [7], where the fiber angle varied linearly between and , from endocardium to epicardium [63].

Fig. 3

Discretization of the cardiac model: (a) The idealized LV was modeled as a thick-walled ellipsoid truncated at of the total height. Typical cardiac dimensions were used (semi-major axis=, semi-minor axis=, wall thickness= at the apex, at the base). The red and blue curves denote the epicardium () and endocardium () fiber directions [63,7], respectively. Zero traction condition was applied on the epicardial surface, and the base was held fixed. (b) Number of elements and degrees of freedom (DOFS) in each discretization and error of the three methods when used in the cardiac cycle test ( for the penalty (PEN) and weakly penalized (WP) methods). (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this article.)

Several hyperelastic constitutive laws have been proposed to model the myocardial tissue. In this work, the myocardium was modeled using the transversely isotropic exponential law introduced by Guccione et al. [3], a model used frequently in the literature.4 This constitutive law is defined with respect to the fiber coordinate system where a coordinate system aligned with local tissue structure is defined everywhere in the material by orthogonal fiber , sheet , and sheet normal unit vectors [5,11]. Letting,denote the components of strain aligned with the local tissue structure directions, then the fiber oriented Green strain is,with . The strain energy and second Piola stress tensor may then be written as [3],where stores the material parameters governing the stress response to strain in fiber/sheet/normal directions, i.e.The parameters used were [25]. The endocardial surface of the ventricular model was passively loaded to (), to cover normal and pathological LV functions at end diastole. The LV was inflated using 150 equal load steps, by setting the boundary traction equal to the product of the pressure and the deformed surface normal.

In order to simulate a cardiac cycle, the cardiac model was modified to include myocardial contraction through an active tension generation model [10]. Active tension generation was incorporated into the cardiac model by the addition of the active stress in the fiber direction of the stress tensor. The cardiac model was also coupled to a Windkessel model representing the systemic circulation using the parameters given by Korakianitis et al. [64]. The coupling was enabled through the use of a Lagrange Multiplier which enforces the same rate of change of LV volume in the two models [65].

The LV was discretized using four different meshes of hexahedral elements (Fig. 3(b)). On the first three discretizations, a quadratic interpolation was used for the displacement. The pressure field was interpolated using linear continuous (for the LM/PL methods) and discontinuous (for the weakly penalized formulation) Lagrange polynomials. In the cardiac tests the results on the first three meshes were compared with the LM and PL solution on mesh4 (for the incompressible and compressible comparison respectively), where a quadratic-linear interpolation scheme was employed.

Numerical solution

The solid mechanics tests presented in Section 3.1 were used to test the convergence behavior of the methods discussed. The convergence rate of each method was acquired by observing the change in the error between a high resolution benchmark solution and the approximate solution, with mesh refinement. As compressible methods may be selected as an approximation to incompressible behavior, we tested the convergence characteristics of the penalty and weakly penalized approaches to the incompressible LM solution (i.e. in Eq. (11)). We also examined the ability of these approaches to model compressible behavior, comparing the results with a fine grid compressible solution(s) (PL solution(s)). The error tolerance for these tests was set to .

The problems under consideration were implemented in CHeart –a multi-physics software tool based on [66-68,55] and expanded by the CHeart team at KCL. All problems were solved on a Dell OPTIPLEX 990, quad-core (Intel® Core i7–2600 CPU @ 3.40 GHz), on an 2.1 GHz AMD Opteron Interlagos 32 processor and on an SGI with 640 2.67 GHz processors (Intel® Xeon® CPU E7–8837).

Numerical results for the convergence rates

The LM, penalty and weakly penalized formulations were initially compared to the incompressible formulation of the elongation problem (Section 3.1.1). Fig. 4 compares the error of the penalty and weakly penalized methods in the solution of the incompressible elongation problem measured over the entire domain as well as a horizontal patch excluding the corners (where singularities in the solution occur). Finally, the importance of interpolation order is highlighted in Fig. 5, where linear interpolation was used for both penalty and weakly penalized formulations (where the local orthogonal projection was selected as the set of piecewise-discontinuous constants).

Fig. 4

Comparison of the convergence behavior of the 3 methods on the two-dimensional elongation problem (the slope of these curves is denoted by ): The error between the (a) penalty and (b) weakly penalized approaches () and the incompressible fine grid solution () for six different values of the bulk modulus k. Convergence of the LM method is shown in black for comparison, whereas the red line represents the highest value of k. (c) Illustration of the error for penalty/weakly penalized () approaches measured over a subset of the domain excluding the region around the four corners of the square.

Fig. 5

Comparison of the convergence behavior of the penalty and weakly penalized formulations when a lower order interpolation scheme is used: The errors between the (red) penalty and (black) weakly penalized forms () and the fine grid solution to the incompressible () elongation problem using linear interpolations are illustrated. (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this article.)

Similar results can be observed in the passive inflation problem detailed in Section 3.1.2. Here convergence of the -norm displacement error in the different methods is shown in Fig. 6 for both approximations to the fine grid incompressible (for ) and compressible passive inflation problems (with and ). The LM, penalty and weakly penalized methods were also compared on the cardiac cycle model showing consistent results to those illustrated in the passive inflation test. Representative results of this comparison are illustrated in Fig. 3(b), while Fig. 7(a) illustrates convergence of the weakly penalized pressure–volume loops with mesh refinement.

Fig. 6

Comparison of the norm error of the various formulations compared to fine grid (a, b) incompressible and (c, d) compressible solutions of the passive inflation problem: (Top) Comparison as an approximation to the incompressible problem (LM solution) for (a) penalty and (b) weakly penalized methods for different k values. (Bottom) Comparison as an approximation to the compressible problem for PL, penalty and weakly penalized (WP) forms with (c) and (d) . denotes the slope of the LM and PL curves.

Fig. 7

Comparison over the cardiac cycle: (a) Display of the pressure–volume loops of the weakly penalized solution of the cardiac cycle on the first three meshes. These pressure–volume loops are converging to the pressure–volume loop of the LM solution on mesh4. (b) The norm (top) and semi-norm (bottom) comparison of the displacement error between the penalty and weakly penalized formulations () in different phases of the cardiac cycle on mesh2. Letters A, B, C, D map the time in cycle (b) to cardiac phase on the pressure–volume loop (a).

Examining the different effects of these methods on cardiac mechanics, we solved each model over a single cardiac cycle, comparing the differences in their behavior, using the LM method as the point of reference. The cardiac cycle was solved on an intermediate mesh (mesh2), using a quadratic-linear interpolation scheme for the displacement and pressure variables. Fig. 7 illustrates the pressure–volume loop derived from the coupled Windkessel-ventricle model as well as the differences between the LM model and both penalty and weakly penalized methods (with ) throughout the cycle. Fig. 9 shows the fiber strain, of the LM method, at the point where the most significant variations are observed ( after end diastole) as well as the absolute difference in fiber strain between the LM solution on every mesh and the respective penalty and weakly penalized approximations.

Fig. 9

Illustration of the absolute difference (error) in fiber strain () at , between the Penalty and LM method (a, d, g) and the weakly penalized and LM method (b, e, h) on the first three meshes (colors representing values of the error between 0 (blue) and (red)). Additionally, the fiber strain for the LM method on each mesh is displayed (c, f, i), with colors ranging from blue () to red (). Figures created in CMGUI [69]. (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this article.)

Numerical results for the efficiency of the different formulations

The PL, penalty and weakly penalized formulations were compared in terms of their nonlinear convergence behavior. Representative values for the number of iterations of the Newton–Rapson scheme, the number of Jacobian matrix computations (and their respective times) along with the linear solution time and total time are presented in Table 1. The improvement in the efficiency of the different methods when the SNR scheme is applied is presented as well. Finally, the effect of the modifications we introduced in the SNR scheme for the penalty method can be deduced as the table compares the application of the SNR scheme to the penalty method with and without the introduced modifications. Note that although not presented here, the nonlinear behavior of the weakly penalized system when the SNR is applied without the introduced modifications is similar to that of the penalty method without the introduced modifications (PEN). Similar observations can be made using Fig. 8 which compares the number of Jacobian and residual computations when the classic Newton–Raphson and the SNR scheme are used for the different methods.

Table 1

Comparison of average number of Newton–Raphson iterations and Jacobian computations per load step as well as their respective average times between Newton–Raphson and Shamanskii–Newton–Raphson schemes. The efficiency of the SNR scheme with (PEN-MOD) and without (PEN) the introduced modifications on the penalty method is presented as well. The total solve time per load step and the total time per load step are illustrated as well. This comparison was performed on the passive inflation test (Section 3.1.2) on mesh2 ( for the penalty, weakly penalized and Perturbed Lagrangian (PL) methods), the simulations were run on a single processor and a direct solver was used.

	J compute time [s]a	J computationsa	R compute time [s]a	R Computations a	Solve time [s]a	Total time [s]a
Newton–Raphson
PEN	181.55	4	1.93	4	47.51	231.13
WP	242.05	3.81	2.42	3.81	41.46	286.06
PL	246.77	3.88	2.49	3.88	45.36	294.75
Shamanskii–Newton–Raphson
PEN	25.80	0.46	11.67	10.82	7.2	44.81
PEN-MOD	10.14	0.113	14.07	9.51	1.76	26.08
WP-MOD	3.95	0.047	13.19	9.79	1.18	18.47
PL	4.66	0.053	13.43	9.88	1.37	19.61

Times/iterations given as the average per load step.

Fig. 8

Comparison of the number of Jacobian and residual computations for (a) the penalty, (b) the weakly penalized and (c) the PL formulation. The methods are compared over the passive inflation simulation (Section 3.1.2) on mesh 2 ( for the penalty, weakly penalized and Perturbed Lagrangian (PL) methods). The solid black line presents the number of Jacobian and residual computations for the classic Newton–Raphson scheme and the dotted lines present the number of Jacobian (black line) and residual (red line) computations when the SNR scheme is applied. (For interpretation of the references to colour in this figure caption, the reader is referred to the web version of this article.)

Discussion

Comparison of the methods for modeling incompressibility

Examining the ability of all approaches to model imcompressibility, as we noted in the introduction the most straightforward is the incompressible LM method which enforces weak incompressibility. However, we could consider the compressible penalty and weakly penalized approaches as approximations to the incompressible system. With this in mind, from Eq. (11), the error for any method should converge to zero with a rate proportional to . For the penalty method, we see from Fig. 4(a) and Fig. 6(a) that the error relative to the incompressible solution was generally higher and increased with increasing k, as a result of the well-known locking phenomena associated with displacement-only formulations. In the cardiac model, with the error was actually uniformly worse than all other values of the parameter at almost all levels of refinement, making the selection of an appropriate k to model incompressible behavior non-trivial.

As discussed in Section 2.2, we initially hypothesized that issues affiliated with the penalty method could be circumvented by projecting the constraint using an orthogonal projection operator, , resulting in the displacement-based formulation suggested by Bercovier [50] and others [48,49]. Indeed, from Fig. 4(b) and Fig. 6(b), we see that as k increased, the error in the approximation decreased proportionally to and became indistinguishable from the convergence of the LM method itself. The existence of a k-dependent error bound for the weakly penalized approach enables regulation of the error by an appropriate choice of k. Moreover, due to its dependence on the discretization, the error showed plateauing behavior for values of k which incured an error larger than the error associated with the discretization.

The locking behavior of the penalty method is observed to worsen with lower order as shown in Fig. 5, where linear elements were used. In this case, as the bulk modulus increased, the rate of convergence observed in the penalty method deteriorated to nearly zero. In contrast, the weakly penalized approach exhibited consistent linear convergence for .

We notice that for both elongation/cardiac problems, the rates of convergence from all methods were not optimal as we would expect based on the error estimates [61]. As the sub-optimal convergence rates appear in the application of all methods, we can assume that this is not a method-dependent issue. We believe that it is due to singularities in the two problems which limit convergence. In the elongation problem, singularities occur at all corners of the domain. Measuring convergence in a horizontal patch excluding corners as shown in Fig. 4(c), we see that for the rate of convergence in the weakly penalized method is restored to the expected order (for the semi-norm). In contrast, due to locking, no improvement to the rate of convergence is observed in the penalty approach. In the cardiac model, sub-optimal convergence is due to fixing the base plane of the model and the singularity in the fiber field near the apex. Even though the specific boundary condition and fiber field incur singularities, they were chosen because of their frequent use in cardiac models.

Comparison of the methods for modeling compressibility

Similar conclusions can be deduced by the application of the PL, penalty and weakly penalized formulations in the solution of compressible problems. In compressible problems, the three formulations should provide consistent results for low and moderate values of k. This is observed in Fig. 6(c). While we observed flattening of the convergence behavior to the incompressible solution for in both penalty and weakly penalized methods (and, though not shown, for PL), we see uniform and consistent convergence to the compressible solution.

Increasing k, however, caused deterioration in the convergence of the penalty method as shown in Fig. 6(d). Here, the error increased by almost an order of magnitude, while convergence remains consistent between the PL and weakly penalized methods. We note that identical behavior was observed in the 2D elongation problem and, while measurement of the error excluding corners, restored optimal convergence rates in the PL and weakly penalized methods, no improvement in the rate was observed in the penalty method.

Comparison over the cardiac cycle

In Fig. 7(b), we compare all methods over a cardiac cycle plotting the difference between both penalty and weakly penalized approaches (with ) and the incompressible LM method on the same discretization. This comparison was performed on an intermediate mesh, mesh2, consisting of 448 elements. Here we see that the LM and penalty methods differ in the semi-norm (which is indicative of errors we could expect in strain) by up to , while the peak difference between weakly penalized and LM approaches remains below . These differences in strain occured primarily during the systolic phase, with decreased error through the rest of the cardiac cycle. Similar conclusions can be deduced from the bulk behavior of these models as seen in Fig. 9, presenting a larger difference in fiber strain between the penalty and LM methods than the weakly penalized and LM methods. Interestingly, in this case the error of both the penalty and weakly penalized formulations decreased with mesh refinement.

The influence of these effects is heavily dependent on the k chosen for the model. Considering convergence (i.e. mesh3 with fine grid mesh4) of the compressible model over the cardiac cycle (even though not shown here), the maximum error for was ∼1% for the weakly penalized and LM methods and ∼10% for the penalty method. However, for the error for weakly penalized and LM methods remained around ∼1% while the error observed in the penalty method dropped to ∼3%. As the bulk modulus represents the tissue’s resistance to compression, its value is tied to the other cardiac constitutive parameters. Thus the influence of locking in the penalty method depends on the level of compressibility which is acceptable in the model. In general it seems that as , where C is the bulk scaling on the strain energy in Section 3.1.2, locking becomes increasingly more predominant.

Comparison of linearized systems and solution

In Section 2.5.2, we outlined the solution procedure for the weakly penalized formulation illustrating that the linearized system involves only the body displacement, . Considering the Jacobian for the LM (), and penalty methods (), shown in Eq. (56), the LM formulation has an indefinite saddle point structure while the penalty method adds to the principle block,Similar to , the Jacobian of the weakly penalized formulation shown in Eq. (34) also augments the block with a matrix, which, by construction, is symmetric positive semi-definite. Further, the Jacobian of the PL method augments the zero block matrix with a k-dependent term, avoiding the indefinite nature of the LM Jacobian (the Jacobians of the different formulations are outlined in Appendix A).

The structure of these systems has a significant impact on their solution. While the actual system sizes (shown in Fig. 2(b) and 3(b)) are not substantially different, the indefinite structure of makes it more challenging to solve, requiring direct methods, “sophisticated” preconditioners or splitting schemes [70]. In contrast, the penalty, PL and weakly penalized strategies are more straightforward in structure, making them more ammenable to classic preconditioning strategies. However, as the bulk modulus k increases, care must be taken to deal with the conditioning of the linear system.

In addition to having contrasting linear structure, the methods also exhibited differing convergence behavior in the Newton–Raphson scheme outlined in Section 2.5.1.5 In general, the non-linear convergence of the weakly penalized formulation averaged ∼3.81 iterations per load step when the classic Newton–Raphson scheme was employed (Table 1). The modifications we introduced in the weakly penalized form (Section 2.5.2), enhanced the numerical ability of the scheme, which exhibited marginally better non-linear convergence behavior than that of the PL method.

Furthermore, the PL and weakly penalized forms were able to exploit the Jacobian re-use strategy (Shamanskii–Newton–Raphson scheme), leading to approximately decrease in the computational time of the Jacobian matrix (build and solution) and a reduction in the total time per loading step (for the weakly penalized form). By modifying the weakly penalized scheme to avoid the high sensitivity to the bulk modulus associated with displacement formulations, the weakly penalized formulation allows efficient re-use of the Jacobian matrix, whereas the performance of the penalty formulation (PEN) is not significantly improved when the SNR scheme is applied. When these modifications were also extended to the penalty method (PEN-MOD), they resulted in significant improvements in both the computational time of the Jacobian matrix ( decrease) and the total time per loading step ( decrease) compared to the classic Newton–Raphson scheme.

Similar conclusions can be deduced from Fig. 8 which compares the Jacobian and residual computations over the iteration number, with and without the SNR scheme. Clearly, the SNR scheme significantly reduces the number of Jacobian computations for all methods. It is important to note that these observations are consistent in all formulations, indicating that the modifications we introduced in the SNR scheme for both the weakly penalized and penalty formulations were able to significantly improve the performance of the methods.

Weakly penalized formulation

In Section 2.1.1 we illustrated how the energy functional for a hyperelastic solid can be written consistently for both penalty and LM methods by choosing both an appropriate space of solutions, , and orthogonal projection, . In the finite element context, we showed that in the LM method both and were necessarily discretized as both the displacement and pressure variables need to be computed, while in the penalty formulation the orthogonal projection is not necessarily discretized as the only unknown variable is the displacement (Section 2.2). As we have shown, this discretization of the orthogonal projection can restrict the approximation space for high values of k.

To circumvent this issue while retaining the single field approach, in Section 2.3 we apply a displacement-only formulation introduced by Bercovier [50] and others [48,49] which uses a localized discrete orthogonal projection operator, . Similar to augmented Lagrangian and reduced integration techniques [26], the aim of the discrete projection is to weaken the compressible/incompressible constraint, thereby enhancing the approximation space in the limit as k gets large. Furthermore, by appropriate restructuring of the weakly penalized system, we avoid the poor nonlinear convergence for high bulk modulus associated with displacement-only formulations. As shown in Figs. 4 and 6, the weakly penalized formulation restores convergence behavior while maintaining the simplicity of a single field approach. Finally, viewing the various methods under the same generalized framework allows us to extend the SNR modifications of the weakly penalized form to the penalty approach, significantly improving the computational performance of the scheme.

A convenient feature of the weakly penalized approach is that it enables more straightforward analysis by tapping into known finite element spaces. Though for uniqueness inf–sup stability is not necessary, this condition ensures optimal convergence in the null space of for linear problems, for appropriately chosen spaces. In the examples presented here, the projection was chosen to be one polynomial order less and piece-wise discontinuous. Though this pairing is not inf–sup stable, for the quadrilateral and hexahedral elements considered, this restored convergence. Another convenient choice are Nicolaides–Boland [60] elements, which give consistent results to those presented here.

Conclusion

In this paper we compared the use of different methods for approximating incompressible and compressible tissue mechanics in the heart. Noting that the choice of model is governed by both model validity and numerical considerations, we assessed the use of Lagrange (LM and PL) and penalty methods as models of both incompressible and compressible behavior. Motivated by the classic locking phenomena observed for linear mechanics [26,27], we apply an enhancement of the Bercovier [50] formulation which enables the single field approach while providing similar convergence behavior to the LM method. To the best of our knowledge this is the first application of this approach on heart models.

Observing the convergence behavior of these methods on a simple solid mechanics test and on cardiac models, we highlight the fact that the LM and penalty methods, although often used equivalently in cardiac mechanics, may present significant variations in results. This is due to the well-known deterioration of the convergence behavior of the penalty method for large values of the bulk modulus. Indeed in both the 2D elongation problem and the cardiac models, the penalty method generally has a larger error than the other two methods for all values of the bulk modulus. In contrast, the single field weakly penalized approach provides both improved rates of convergence and avoids issues associated with locking phenomena over these test problems. Further modifications introduced in this work enhance the computational performance of the numerical scheme, by allowing efficient application of the SNR re-use strategy and significantly reducing the computational time. The weakly penalized formulation can therefore provide an accurate and computationally inexpensive method that can be used to deal with incompressibility and near incompressibility in problems of cardiac mechanics and solid mechanics in general.

Algorithm 1. Shamanskii–Newton–Raphson (SNR)
Given initial guess U^0,β, and S_0. Compute R(U^0), and set n=0.
whileS_n‖R(U^n)‖>TOL
if (β=n) {ComputeJ(U^n)and[J(U^n)]^-1 (or preconditioner)}
SolveJ(U^β)δU^′=-R(U^n).
Findmin{‖R(U^n+α^nδU^′)‖,α^n∈[0,1]}.
SetδU^n=α^nU^′,U^n+1=U^n+δU^n.
if (‖R(U^n+1)‖>γS_n‖R(U^n)‖ orn> ITER) {Setβ=n+1}
Setn=n+1, andS_n=1
end while