A Block-Coupled Finite Volume Method for Incompressible Hyperelastic Solids

1. Introduction

A considerable range of engineering applications involve materials capable of sustaining large-strain elastic deformations while exhibiting incompressible or nearly incompressible mechanical behavior. Rubber is a canonical example that is capable of undergoing substantial stretching or bending without an appreciable change in volume. Many soft biological tissues, including muscle, skin, and arterial walls, are commonly idealized as incompressible because of their high water content, which effectively suppresses volumetric change. Hydrogels and silicone gels exhibit similar characteristics, and many engineering elastomers (such as polyurethane and silicone rubber) are also modeled as incompressible in mechanical analyses.

It is well established that numerical stress analysis of solids approaching the incompressible limit presents significant challenges, particularly when a displacement-based formulation is employed. In particular, fully integrated, displacement-based lower-order finite elements often exhibit volumetric locking, a phenomenon frequently accompanied by pressure oscillations. Both volumetric locking and pressure oscillations have been extensively studied, leading to the development of various strategies aimed at mitigating these issues. The most widely used approach to prevent volumetric locking within the Finite Element (FE) framework is the mixed displacement-pressure formulation [1]. This method introduces pressure as an independent variable alongside displacement, allowing the problem to be formulated in a way that enforces the incompressibility constraint without locking.

In Computational Solid Mechanics (CSM), the FE method remains the most widely used approach. However, the Finite Volume (FV) method - long established in Computational Fluid Dynamics - has been steadily expanding its reach into CSM. Its growing popularity is due to its conceptual simplicity and strong conservation properties. Over the years, the FV method has been applied to a broad spectrum of CSM problems, emerging in several distinct variants based on different control volume arrangements, most notably cell-centred and vertex-centred schemes. A comprehensive review is provided in [2].

Only a limited number of FV methods have been proposed for incompressible bodies [3]. The present study builds on the work of Professor I. Demirdžić and collaborators [4,5,6], who employed a mixed displacement–pressure formulation to prevent volumetric locking. The mathematical model consists of the momentum equation, the mass conservation equation (incompressibility condition) and the constitutive equation. Spatial discretisation is performed using a cell-centred FV method applied to meshes with control volumes of completely arbitrary topology. The resulting set of coupled linear algebraic equations is solved using a segregated procedure incorporating a SIMPLE-based algorithm [7], which ensures the required pressure–displacement coupling through the use of Rhie–Chow interpolation [8] to evaluate displacement at the control volume faces. In the present study, the following improvements to the above models are proposed:

• Block-coupled solution procedure – Instead of using a segregated solution procedure, in which displacement increment components and the pressure increment are solved separately, a block-coupled system of linear algebraic equations (representing the discretised momentum and pressure increment equations) is solved for displacement and pressure increments simultaneously. This leads to a more robust and efficient solution process, as there is no need to apply user-defined under-relaxation to the displacement and pressure increments. The proposed block-coupled solution procedure is based on the procedure proposed by Professor M. Darwish and collaborators [9,10] for the solution of laminar incompressible flow problems on collocated unstructured FV grids. Later, a similar algorithm was implemented [11,12] in the OpenFOAM framework [13], which serves as the starting point for the implementation of an algorithm for solving incompressible elastic/hyperelastic solid deformations in this work.
• Temporal discretisation - Emphasis is placed on the temporal accuracy of the model, as the ultimate objective is its application to fluid–structure interaction simulations in vascular flows. Four commonly used temporal discretisation schemes are implemented and tested, each incorporating a temporally consistent Rhie–Chow interpolation [14].
• Improved treatment of traction boundaries - Along the portions of the discretized spatial domain boundary where traction is specified, the displacement increment is calculated using the cell-face normal component obtained from the solution of the continuity (pressure increment) equation, while the remaining components are reconstructed using the selected constitutive equation. This approach provides substantially higher accuracy at the traction boundaries compared to the second-order extrapolation used in [4].
• Extended material applicability – The proposed FV solver is applicable to both linear elastic bodies and nonlinear hyperelastic materials described by arbitrary constitutive equations, with emphasis on constitutive relations used for modeling arterial walls and heart tissue (for exampe, Holzapfel-Gasser-Ogden (HGO) model [15] and Guccione model [16]).

The proposed FV solver for incompressible elastic/hyperelastic bodies is implemented within the solids4Foam toolbox [17,18], which extends the OpenFOAM framework to allow the straightforward implementation of advanced solid mechanics and fluid–structure interaction (FSI) procedures. In this way, the proposed FV solver provides the foundation for developing a self-contained FSI solver for vascular flows in the next step.

In the next section, the governing equations and constitutive relations are presented, followed by a brief description of the FV discretisation procedure and the solution algorithm. Finally, the capabilities of the method are demonstrated through a series of test cases.

2. Mathematical Model

The deformation of an arbitrary body of volume $\mathsf{\Omega }$, bounded by a surface $\mathsf{\Gamma }$, is governed by the law of conservation of linear momentum, which can be expressed in the following integral form:

$${\displaystyle \frac{\partial }{\partial t}}\underset{\mathsf{\Omega }}{\int }\rho {\displaystyle \frac{\partial \mathit{u}}{\partial t}}d\mathsf{\Omega }=\underset{\mathsf{\Gamma }}{\int }\mathit{n}·\sigma d\mathsf{\Gamma },$$

(1)

where t is time, $\rho$ is density, $\mathit{u}$ is the displacement vector, $\mathit{n}$ is the outward-pointing normal vector of the surface $\mathsf{\Gamma }$, $\sigma$ is the Cauchy stress tensor and the body force is omitted due to clarity. Equation (1) is integrated over the body in the deformed configuration. In this study, the alternative formulation of the linear momentum conservation law is used where integration is performed over the body in the initial undeformed configuration:

$${\displaystyle \frac{\partial }{\partial t}}\underset{{\mathsf{\Omega }}^0}{\int }{\rho }^0{\displaystyle \frac{\partial \mathit{u}}{\partial t}}d{\mathsf{\Omega }}^0=\underset{{\mathsf{\Gamma }}^0}{\int }\left(J{\mathit{F}}^{-\mathrm{T}}·{\mathit{n}}^0\right)·\sigma d{\mathsf{\Gamma }}^0,$$

(2)

where the superscript 0 indicates that a quantity refers to the initial configuration ($t=0\mathrm{s}$); for example, ${\mathsf{\Omega }}^0$ denotes the initial undeformed volume bounded by the corresponding initial undeformed surface ${\mathsf{\Gamma }}^0$. In the surface integral on the right-hand side of Equation (2), Nanson’s formula is used to transform the elementary area vector from the initial to the current configuration, where $\mathit{F}=\mathit{I}+{\left({\nabla }_0\mathit{u}\right)}^{\mathrm{T}}$ is the deformation gradient tensor, $J=det\mathit{F}$ is the Jacobian of the deformation gradient, $\mathit{I}$ is the unit (identity) tensor, and ${\nabla }_0$ denotes the gradient operator with respect to the initial configuration1. To close the model, it is necessary to specify the constitutive relation between the Cauchy stress tensor $\sigma$ and the displacement vector $\mathit{u}$, which is presented in the following section.

2.1. Constitutive Equations

The ultimate goal of this study is to develop a finite-volume solver for the simulation of large-strain dynamic deformation in elastic bodies composed of incompressible material, with a particular focus on applications in vascular system modelling. The solver currently handles incompressible elasticity/hyperelasticity but can be readily extended to compressible materials, following the same numerical approach as in [19]. Moreover, it consistently recovers the limit of linear elasticity, corresponding to small strains and rotations.

Constitutive relation for an isotropic linear elastic solid valid only under the small strain and small rotation assumption is Hook’s law:

$$\sigma =2\mu ϵ+\lambda \left(\nabla ·\mathit{u}\right)\mathit{I},$$

(3)

where $\mu$ and $\lambda$ are Láme coefficients and $\mathit{psy}\mathit{e}$ is the linear strain tensor:

$$ϵ={\displaystyle \frac{1}{2}}\left[\nabla \mathit{u}+{\left(\nabla \mathit{u}\right)}^{\mathrm{T}}\right].$$

(4)

The Láme coefficients $\mu$ and $\lambda$ as well as the bulk elasticity modulus K can be defined through Young’s elasticity modulus E and Poisson’s ratio $\nu$ by the following relations:

$$\mu ={\displaystyle \frac{E}{2(1+\nu )}},\lambda ={\displaystyle \frac{E\nu }{(1+\nu )(1-2\nu )}},K={\displaystyle \frac{E}{3(1-2\nu )}}.$$

(5)

It can be observed that Equation (3) is not well defined in the incompressible limit sice for $\nu \to {\displaystyle \frac{1}{2}}$, $\lambda \to \infty$ and $\nabla ·\mathit{u}\to 0$. To avoid this issue, pressure can be introduced as an additional independent variable, defined by the following equation:

$$p=-K\nabla ·\mathit{u}.$$

(6)

Taking into account Equation (6), Hook’s law (3) can be rewritten in a mixed (pressure-displacement) formulation as follows:

$$\begin{array}{c}\sigma =2\mu ϵ-{\displaystyle \frac{\lambda }{K}}p\mathit{I},\end{array}$$

(7)

$$\begin{array}{c}\nabla ·\mathit{u}=-{\displaystyle \frac{p}{K}},\end{array}$$

(8)

In the limit of incompressibility $\nu \to {\displaystyle \frac{1}{2}}$ while $\lambda /K\to 1$ and $1/K\to 0$, hence Equations (7) and (8) reduce to

$$\begin{array}{c}\sigma =2\mu ϵ-p\mathit{I},\end{array}$$

(9)

$$\begin{array}{c}\nabla ·\mathit{u}=0,\end{array}$$

(10)

which represents well defined constitutive relation for incompressible linear elastic material. Equation (10) represents the continuity equation for an incompressible material and is also referred to as the incompressibility condition.

Nonlinear elastic behaviour can be represented by various hyperelastic models. Here, the Neo-Hookean law is used as a reference, and the solver is structured to easily incorporate alternative formulations.

Deformation of a compressible Neo-Hookean solid (Ogden variant) can be described with the following constitutive relation in the mixed (pressure-displacement) formulation:

$$\sigma ={\displaystyle \frac{\mu }{J}}(\mathit{B}-\mathit{I})-{\displaystyle \frac{\lambda }{K}}p\mathit{I},$$

(11)

where $\mathit{B}=\mathit{F}·{\mathit{F}}^{\mathrm{T}}$ is the left Cauchy-Green tensor. A mixed formulation is again used in order to facilitate modelling of incompressible materials. In the limit of incompressibility Equation (11) reduces to:

$$\sigma =\mu (\mathit{B}-\mathit{I})-p\mathit{I}.$$

(12)

The constitutive relations described above are expressed in terms of the total displacement $\mathit{u}$, which is consistent with the total displacement formulation of the linear momentum conservation law (2). In this study, incremental displacement formulation of the linear momentum conservation law is applied since this approach can have positive effect on the efficiency and stability of an numerical procedure. Hance, it is necessary to express stress increment in terms of displacement increment $\delta \mathit{u}=\mathit{u}-{\mathit{u}}^n$, where superscript n refers to the known quantities from the previous time step2. Accordingly, one can express Cauchy stress increment tensor in terms of displacement increment and pressure increments for the incompressible Hook’s law as follows:

$$\delta \sigma =\sigma -{\sigma }^n=\mu \left[\nabla \delta \mathit{u}+{(\nabla \delta \mathit{u})}^{\mathrm{T}}\right]-\delta p\mathit{I},$$

(13)

where $\delta p$ is the pressure increment. Cauchy stress increment for incompressible Neo-Hookean solid reads as follows:

$$\delta \sigma =\mu \left[{\mathit{F}}^n·\left({\nabla }_0\delta \mathit{u}\right)+{\left({\nabla }_0\delta \mathit{u}\right)}^{\mathrm{T}}·{\left({\mathit{F}}^n\right)}^{\mathrm{T}}\right]-\delta p\mathit{I},$$

(14)

where the nonlinear quadratic term ${\left({\nabla }_0\delta \mathit{u}\right)}^{\mathrm{T}}·{\nabla }_0\delta \mathit{u}$ is omitted. One can note that expression (14) reduces to expression (13) if assumption of linear elasticity is satisfied ($\mathit{F}\approx \mathit{I}$).

2.2. Resulting Set of Equations for Incompressible Solid

Resulting set of governing equation which are solved numerically consists of linear momentum conservation law and mass conservation law (incompressibility conditions). As is mentioned earlier, in order to enhance numerical efficiency, the incremental formulation of the conservation laws is used.

Incremental formulation of the conservation law for total linear momentum reads as follows:

$${\displaystyle \frac{\partial }{\partial t}}\underset{{\mathsf{\Omega }}^0}{\int }{\rho }^0{\displaystyle \frac{\partial (\delta \mathit{u}+{\mathit{u}}^n)}{\partial t}}d{\mathsf{\Omega }}^0=\underset{{\mathsf{\Gamma }}^0}{\int }\left[{\left({\mathit{F}}^n\right)}^{-\mathrm{T}}·{\mathit{n}}^0\right]·\left[\delta \sigma +{\sigma }^n\right]d{\mathsf{\Gamma }}^0.$$

(15)

By expressing the Cauchy stress increment $\delta \sigma$ in the above equation in terms of the unknown displacement increment $\delta \mathit{u}$ and pressure increment $\delta p$ via the incremental constitutive relation (14), the balance of linear momentum takes the following final form:

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial }{\partial t}}\underset{{\mathsf{\Omega }}^0}{\int }{\rho }^0{\displaystyle \frac{\partial (\delta \mathit{u}+{\mathit{u}}^n)}{\partial t}}d{\mathsf{\Omega }}^0=\underset{{\mathsf{\Gamma }}^0}{\int }\delta \mathit{T}d{\mathsf{\Gamma }}^0+\underset{{\mathsf{\Gamma }}^0}{\int }{\mathit{T}}^{n}d{\mathsf{\Gamma }}^0,\end{array}$$

(16)

where the traction from the previous time instance ${\mathit{T}}^n$ and the traction increment $\delta \mathit{T}$ are defined as follows:

$${\mathit{T}}^n=\left[{\left({\mathit{F}}^n\right)}^{-\mathrm{T}}·{\mathit{n}}^0\right]·{\sigma }^n,$$

(17)

$$\begin{array}{cc}\hfill \delta \mathit{T}=\left[{\left({\mathit{F}}^{\mathrm{n}}\right)}^{-\mathrm{T}}·{\mathit{n}}^0\right]·\delta \sigma ={\mu }_{\mathrm{eff}}& \left\{({\mathit{n}}^0·{\nabla }_0\delta \mathit{u})+\left[{\mathit{n}}^0·{\left({\nabla }_0\delta \mathit{u}\right)}^{\mathrm{T}}\right]\right\}-\hfill \\ \hfill & -\left[{\left({\mathit{F}}^n\right)}^{-\mathrm{T}}·{\mathit{n}}^0\right]\delta p.\hfill \end{array}$$

(18)

Equation (16), with tractions defined by Equations (17) and (18), is not exact even for the Neo-Hookean model because, in its derivation, the current deformation gradient is approximated by that from the previous time step ($\mathit{F}\approx {\mathit{F}}^n$) and the quadratic term in the incremental constitutive relation (14) is neglected. In addition, the effective or instantaneous shear modulus ${\mu }_{\mathrm{eff}}$ is used in Equation (18) instead of the shear modulus $\mu$ of the Neo-Hookean model3. The procedure used to calculate aproximate value of instantaneous shear modulus is described in Appendix A. To compensate for these approximations, the following traction-correction term is added to the right-hand side of Equation (16):

$${\mathit{T}}^{\prime }=\left[{\left({\mathit{F}}^{*}\right)}^{-\mathrm{T}}·{\mathit{n}}^0\right]·\sigma ({\mathit{F}}^{*},p^{*})-(\delta {\mathit{T}}^{*}+{\mathit{T}}^n),$$

(19)

where superscript * represents quantities obtained in the previous iteration of the iterative solution procedure4. In such a way, the Cauchy stress tensor $\sigma ({\mathit{F}}^{*},p^{*})$ is evaluated using selected hyperelastic constitutive model based on the values of the unknown displacement (displacement gradient) and pressure fields from the previous iteration. Taking into account traction correction term defined by Equation (19), final form of linear momentum equation reads as follows:

$$\begin{array}{c}\hfill {\displaystyle \frac{\partial }{\partial t}}\underset{{\mathsf{\Omega }}^0}{\int }{\rho }^0{\displaystyle \frac{\partial (\delta \mathit{u}+{\mathit{u}}^n)}{\partial t}}d{\mathsf{\Omega }}^0=\underset{{\mathsf{\Gamma }}^0}{\int }\delta \mathit{T}d{\mathsf{\Gamma }}^0+\underset{{\mathsf{\Gamma }}^0}{\int }{\mathit{T}}^{n}d{\mathsf{\Gamma }}^0+\underset{{\mathsf{\Gamma }}^0}{\int }{\mathit{T}}^{\prime }d{\mathsf{\Gamma }}^0.\end{array}$$

(20)

In the case of linear-elastic material, the traction correction term vanishes and Equations (17) and (18) simplifies by taking into account that ${\mathit{F}}^n=\mathit{I}$.

In order to close the model which consist of two unknown fields (displacement increment and pressure increment) it is necessary to take into account Equation (10) transformed into incremental integral form:

$$\underset{\mathsf{\Gamma }}{\int }\mathit{n}·\delta \mathit{u}d\mathsf{\Gamma }=0,$$

(21)

where integration is performed over the volume of the body in the current (deformed) configuration. Equation (21) does not represent governing equation for pressure increment but additional constrain for displacement increment. The so called pressure increment equation will be derived later by combining discretised counterparts of Equations (20) and (21).

2.3. Initial and Boundary Conditions

To complete the mathematical model, initial and boundary conditions must be specified. For the numerical solution of transient problems, it is necessary to prescribe the displacement and velocity at the initial time and at all computational points of the spatial solution domain. Boundary conditions must be specified at all times along all spatial domain boundaries. The most common boundary conditions in solid mechanics are prescribed displacement (Dirichlet) and prescribed traction (Neumann). These conditions may be constant or time-dependent.

3. Numerical Model

The mathematical model is discretised in space using a second-order accurate collocated (cell-centered) unstructured FV method, while numerical integration in time is performed using various implicit schemes. The discretisation procedure is divided into two parts: discretisation of the solution domain and discretisation of the governing equations.

3.1. Solution Domain Discretization

The solution domain consist of temporal solution domain (ie. total solution time) and spatial solution domain. The total solution time is divided into a finite number of time steps $\Delta t$, and the equations are solved in a time-marching manner. In the unstructured FV discretisation, the spatial solution domain is partitioned into a finite number of convex polyhedral control volumes (CVs) or cells, each bounded by convex polygonal faces. The cells do not overlap and completely fill the spatial domain. Figure 1 illustrates a simple polyhedral control volume ${\mathsf{\Omega }}_P$ around the computational point P located at its centroid, the face f, the face area ${\mathsf{\Gamma }}_f$, the unit normal vector ${\mathit{n}}_f$ to the face, and the centroid N of the neighbouring CV sharing the face f. In general case, the interpolation point $f^{\prime }$ does not coincide with the face centre point f and distance between them is defined by the skewness correction vector ${\mathit{k}}_f$. The geometry of a CV is fully determined by the positions of its vertices.

3.2. Governing Equations Discretization

The second-order FV discretisation of an integral conservation equation transforms the surface integrals into sums of face integrals and approximates them and the volume integrals to the second-order accuracy by using the midpoint rule. Temporal discretisation is carried out by numerical integration of the governing equation in time from the old time instance $t^n$ to the new time instance $t^{n+1}=t^n+\Delta t$.

3.2.1. Momentum Equation Discretisation

The spatially discretized momentum equation (20) for the control volume ${\mathsf{\Omega }}_P$ can be written as:

$${\rho }_P^0{\displaystyle \frac{\partial }{\partial t}}\left[{\displaystyle \frac{\partial (\delta {\mathit{u}}_P+{\mathit{u}}_P^n)}{\partial t}}\right]{\mathsf{\Omega }}_P^0=\underset{\mathrm{implicit}}{\underbrace{\sum _f\delta {\mathit{T}}_f{\mathsf{\Gamma }}_f^0}}+\underset{\mathrm{explicit}}{\underbrace{\sum _f{\mathit{T}}_f^n{\mathsf{\Gamma }}_f^0+\sum _f{\mathit{T}}_f^{\prime }{\mathsf{\Gamma }}_f^0,}}$$

(22)

where the subscripts P and f denote the cell-centre and face-centre values, respectively, with the summation taken over all faces bounding cell P. In Equation (22), the first term on the right-hand side is treated implicitly after discretisation, whereas the last two terms are explicit and evaluated using displacement and pressure fields (or their increments) from the previous time step or iteration.

The traction increment $\delta {\mathit{T}}_f$ in the implicit term on the right-hand side of Equation (22) is obtained by discretizing Equation (18) at the centers of internal faces as follows:

$$\begin{array}{cc}\hfill \delta {\mathit{T}}_f={\left({\mu }_{\mathrm{eff}}\right)}_f& {({\mathit{n}}^0·{\nabla }_0\delta \mathit{u})}_f+{\left({\mu }_{\mathrm{eff}}\right)}_f{\left[{\mathit{n}}^0·{\nabla }_0(\delta \mathit{u}·{\mathit{n}}^0{\mathit{n}}^0)\right]}_f+\hfill \\ \hfill & +{\left({\mu }_{\mathrm{eff}}\right)}_f{\left[{\nabla }_{0t}(\delta \mathit{u}·{\mathit{n}}^0)\right]}_f^{*}-\left[{\left({\mathit{F}}_f^n\right)}^{-\mathrm{T}}·{\mathit{n}}_f^0\right]\delta p_f,\hfill \end{array}$$

(23)

where the second and third terms on the right-hand side are obtained by decomposition of the gradient transpose term as follows:

$${\mathit{n}}^0·{\left({\nabla }_0\delta \mathit{u}\right)}^{\mathrm{T}}={\nabla }_0(\delta \mathit{u}·{\mathit{n}}^0)={\mathit{n}}^0·{\nabla }_0(\delta \mathit{u}·{\mathit{n}}^0{\mathit{n}}^0)+{\nabla }_{0t}(\delta \mathit{u}·{\mathit{n}}^0),$$

(24)

where ${\nabla }_{0t}=(\mathit{I}-{\mathit{n}}^0{\mathit{n}}^0)·{\nabla }_0$ is the tangential gradient operator.

The face normal derivative of displacement increment in Equation (23) is discretized using the central scheme with non-orthogonal and skewness correction as follows (see Figure 1):

$${({\mathit{n}}^0·{\nabla }_0\delta \mathit{u})}_f={\displaystyle \frac{\delta {\mathit{u}}_N-\delta {\mathit{u}}_P}{d_{f,n}^0}}+{\displaystyle \frac{{\mathit{k}}_N^0·{\left({\nabla }_0\delta \mathit{u}\right)}_N^{*}-{\mathit{k}}_P^0·{\left({\nabla }_0\delta \mathit{u}\right)}_P^{*}}{d_{f,n}^0}},$$

(25)

where $d_{f,n}^0={\mathit{n}}_f^0·{\mathit{d}}_f^0$, ${\mathit{k}}_P^0=(\mathit{I}-{\mathit{n}}_f^0{\mathit{n}}_f^0)·({\mathbf{R}}_f^0-{\mathbf{R}}_P^0)$ and ${\mathit{k}}_N^0=(\mathit{I}-{\mathit{n}}_f^0{\mathit{n}}_f^0)·({\mathbf{R}}_N^0-{\mathbf{R}}_f^0)$. The first term on the right-hand side of Equation (25) is treated implicitly, while the correction term is treated explicitly. The discretisation scheme in Equation (25) is also applied to the face-normal derivative of the normal displacement increment, ${[{\mathit{n}}^0·{\nabla }_0(\delta \mathit{u}·{\mathit{n}}^0{\mathit{n}}^0)]}_f$, where $\delta \mathit{u}$ is replaced by its normal component, $\delta \mathit{u}·\left({\mathit{n}}^0{\mathit{n}}^0\right)$.

The face-centre value of the pressure increment $\delta p_f$ in Equation (23) is calculated using linear interpolation of the neighbouring cell-centre values with the skewness correction:

$$\delta p_f={{}^0\overline{\overline{\left(\delta p_P\right)}}}_f={{}^0\overline{\left(\delta p_P\right)}}_f+{\mathit{k}}_f^0·{{}^0\overline{\left({\nabla }_0\delta p\right)}}_f^{*},$$

(26)

where ${\mathit{k}}_f^0$ is the skewness correction vector between the cell-face interpolation point $f^{\prime }$ and the cell-face centre f, as illustrated in Figure 1. The linear interpolation procedure is then defined as follows:

$${{}^0\overline{\left(\delta p_P\right)}}_f=g_x^0\delta p_P+(1-g_x^0)\delta p_N,$$

(27)

here $g_x^0={\left(\overline{f^{\prime }N}\right)}^0/{\left(\overline{PN}\right)}^0$ is the interpolation factor (see Figure 1), and obtained cell-face value corresponds to the value at the interpolation point $f^{\prime }$, which coincides with the cell-face centre f for a mesh without skewness5. The skewness correction term in Equation (26) is treated explicitly, meaning it is evaluated using the pressure increment field obtained from the previous iteration.

Fully discretised momentum equation is obtained by numerical integration of Equation (22) over the time interval $[t^n,t^{n+1}]$, where approximation of the time integrals is carried out using implicit two-level Euler scheme, Crank-Nicolson (trapezoidal) scheme and implicit three-level backward scheme. Backward and trapezoidal schemes are also combined to form so called composite scheme proposed by Bathe [20]. The fully discretized momentum equation obtained by applying the above listed temporal discretisation schemes reads as follows:

• Euler scheme

  $${\rho }_P^0{\displaystyle \frac{{\mathit{v}}_P^{n+1}-{\mathit{v}}_P^n}{\Delta t}}{\mathsf{\Omega }}_P^0=\sum _f\delta {\mathit{T}}_f^{n+1}{\mathsf{\Gamma }}_f^0+\sum _f{\mathit{T}}_f^n{\mathsf{\Gamma }}_f^0+\sum _f{\left({\mathit{T}}_f^{\prime }\right)}^{*}{\mathsf{\Gamma }}_f^0,$$

  (28)

  $${\mathit{v}}_P^{n+1}={\displaystyle \frac{{\mathit{u}}_P^{n+1}-{\mathit{u}}_P^n}{\Delta t}}={\displaystyle \frac{\delta {\mathit{u}}_P^{n+1}}{\Delta t}},$$

  (29)
• Crank-Nicolson scheme

  $$\begin{array}{cc}\hfill {\rho }_P^0{\displaystyle \frac{{\mathit{v}}_P^{n+1}-{\mathit{v}}_P^n}{\Delta t}}{\mathsf{\Omega }}_P^0=& {\displaystyle \frac{1}{2}}\sum _f\delta {\mathit{T}}_f^{n+1}{\mathsf{\Gamma }}_f^0+{\displaystyle \frac{1}{2}}\sum _f{\mathit{T}}_f^n{\mathsf{\Gamma }}_f^0+\hfill \\ \hfill & {\displaystyle \frac{1}{2}}\sum _f{\left({\mathit{T}}_f^{\prime }\right)}^{*}{\mathsf{\Gamma }}_f^0+{\displaystyle \frac{1}{2}}\sum _f{\mathit{T}}_f^n{\mathsf{\Gamma }}_f^0,\hfill \end{array}$$

  (30)

  $${\mathit{v}}_P^{n+1}=2{\displaystyle \frac{{\mathit{u}}_P^{n+1}-{\mathit{u}}_P^n}{\Delta t}}-{\mathit{v}}_P^n=2{\displaystyle \frac{\delta {\mathit{u}}_P^{n+1}}{\Delta t}}-{\mathit{v}}_P^n,$$

  (31)
• backward scheme

  $$\begin{array}{cc}\hfill {\rho }_{\mathrm{P}}^0& {\displaystyle \frac{3{\mathit{v}}_P^{n+1}-4{\mathit{v}}_P^n+{\mathit{v}}_P^{n-1}}{2\Delta t}}{\mathsf{\Omega }}_{\mathrm{P}}=\hfill \\ \hfill & \sum _f\delta {\mathit{T}}_f^{n+1}{\mathsf{\Gamma }}_f^0+\sum _f{\mathit{T}}_f^n{\mathsf{\Gamma }}_f^0+\sum _f{\left({\mathit{T}}_f^{\prime }\right)}^{*}{\mathsf{\Gamma }}_f^0,\hfill \end{array}$$

  (32)

  $${\mathit{v}}_P^{n+1}={\displaystyle \frac{{\mathit{u}}_P^{n+1}-4{\mathit{u}}_P^n+{\mathit{u}}_P^{n-1}}{2\Delta t}}={\displaystyle \frac{\delta {\mathit{u}}_P^{n+1}-3{\mathit{u}}_P^n+{\mathit{u}}_P^{n-1}}{2\Delta t}}.$$

  (33)
• composite scheme [21]: a simplified version is implemented here, where the Crank-Nicolson and backward scheme are applied alternately.

In the above expressions, the velocity is introduced as an auxiliary variable, defined as the temporal derivative of the displacement, $\mathit{v}={\displaystyle \frac{\partial \mathit{u}}{\partial t}}$. The velocity at the current time step is subsequently expressed in terms of the displacement and its increment by employing the same temporal discretisation scheme.

The fully discretised momentum equation can now be written in the form of a linear algebraic equation, which for cell P takes the form:

$$\begin{array}{cc}\hfill {\mathbf{A}}_P^{\delta \mathit{u}}·\delta {\mathit{u}}_P^{n+1}+& \sum _N{\mathbf{A}}_N^{\delta \mathit{u}}·\delta {\mathit{u}}_N^{n+1}+\hfill \\ \hfill & {\mathbf{A}}_P^{\delta \mathit{u},\delta p}·\delta p_P^{n+1}+\sum _N{\mathbf{A}}_N^{\delta \mathit{u},\delta p}·\delta p_N^{n+1}={\mathbf{R}}_P^{\delta \mathit{u}},\hfill \end{array}$$

(34)

where ${\mathbf{A}}_P^{\delta \mathit{u}}$ and ${\mathbf{A}}_N^{\delta \mathit{u}}$ denote the central (diagonal) and neighbour (off-diagonal) tensorial coefficients associated with the displacement increment field, while ${\mathbf{A}}_P^{\delta \mathit{u},\delta p}$ and ${\mathbf{A}}_N^{\delta \mathit{u},\delta p}$ denote the corresponding vectorial coefficients associated with the pressure increment field. The source vector ${\mathbf{R}}_P^{\delta \mathit{u}}$ on the right-hand side depends on the solution due to the explicit treatment of certain terms (e.g., traction, non-orthogonal, and skewness corrections). The central and neighbour coefficients in Equation (34) as well as the right-hand side vector are defined in Appendix B for the first-order implicit Euler scheme.

3.2.2. Discretised Pressure Equation

Discretised momentum equation represented by Equation (34) must be complemented by an additional system of algebraic equation for unknown cell-centre pressure increment field. In this study, it is applied approach usually used to numerically solve incompressible fluid flow by applying collocated (cell-centered) finite volume method, where discretised pressure equation is derived by combining discretised momentum and continuity equations using Rhie-Chow momentum interpolation procedure [8].

For cell P, the discrete form of the incompressibility condition (21) reads:

$$\sum _f{\mathit{n}}_f^{1+1/2}·\delta {\mathit{u}}_f^{n+1}{\mathsf{\Gamma }}_f^{1+1/2}=0,$$

(35)

where the integration is performed on the cell faces mapped to the configuration at $t=t^{n+1/2}$, which represents the intermediate state between the previous and current configurations and will hereafter be denoted by an m subscript or superscript, that is, $t^m=t^{n+1/2}$. This approach yields the highest accuracy, as shown in [4]. Taking into account above introduced convention, simplified representation of the discrete form of the incompressibility condition now reads:

$$\sum _f{\mathit{n}}_f^m·\delta {\mathit{u}}_f^{n+1}{\mathsf{\Gamma }}_f^m=0.$$

(36)

The cell-face volume increment $\delta {\mathsf{\Omega }}_f^{n+1}={\mathit{n}}_f^m·\delta {\mathit{u}}_f^{n+1}{\mathsf{\Gamma }}_f^m$ in Equation (35) is calculated using the Rhie-Chow interpolation procedure [8] as follows:

$$\begin{array}{cc}\hfill \delta {\mathsf{\Omega }}_f^{n+1}=& {\mathit{n}}_f^m·{{}^m\overline{\overline{\left(\delta \mathit{u}\right)}}}_f^{n+1}{\mathsf{\Gamma }}_f^m+\delta {\mathsf{\Omega }}_f^{\mathrm{CRC}}-\hfill \\ \hfill & -{\left[{{}^m\overline{\left({\displaystyle \frac{A_P^{\delta \mathit{u}}}{{\mathsf{\Omega }}_P^0}}\right)}}_f\right]}^{-1}\left[{\widehat{\mathit{d}}}_f^m·{\left({\nabla }_m\delta p\right)}_f^{n+1}-{\widehat{\mathit{d}}}_f^m·{{}^m\overline{\left({\nabla }_m\delta p\right)}}_f^{*}\right]{\mathsf{\Gamma }}_f^m,\hfill \end{array}$$

(37)

where the interpolated cell-face displacement increment ${{}^m\overline{\overline{\left(\delta \mathit{u}\right)}}}_f^{n+1}$ is corrected by a term proportional to the difference between the pressure increment derivative interpolated to the face and the derivative computed directly at the face. This correction acts as a low-pass filter to prevent pressure-displacement decoupling and suppress oscillations in the pressure increment field. The diagonal coefficient $A_P^{\delta \mathit{u}}$ is scalar quantity and consists only contributions from the temporal and normal derivative term in the discretised momentum equation [first two therms at the right-hand side of Equation (A1)]6. Following the block-coupled solution procedure proposed for incompressible fluid flow in [9], only the interpolated gradient of pressure increment term is treated explicitly after discretisation, while the two reminding terms are implicit. The cell-face derivative of pressure increment ${\widehat{\mathit{d}}}_f^m·{\left({\nabla }_m\delta p\right)}_f^{n+1}$ in Equation (37) is discretised using central differencing formula:

$${\widehat{\mathit{d}}}_f^m·{\left({\nabla }_m\delta p\right)}_f^{n+1}={\displaystyle \frac{\delta p_N^{n+1}-\delta p_P^{n+1}}{|{\mathit{d}}_f^m|}},$$

(38)

where ${\widehat{\mathit{d}}}_f^m={\mathit{d}}_f^m/\left|{\mathit{d}}_f^m\right|$ is the unit vector along the line $\overline{PN}$ (see Figure 1). Here, it should be noted that the pressure derivative is taken along the line $\overline{PN}$ connecting two neighbouring cell centers instead along cell-face normal direction. According to [14], this simplification does not violate the efficiency of the low-pass filter.

Rhie-Chow momentum-weighted interpolation is widely used to prevent pressure-velocity decoupling in simulations of incompressible flow on meshes with a collocated variable arrangement. The original Rhie-Chow interpolation [8] was proposed for steady-state fluid flow simulations, and its direct application to transient problems can result in temporally inconsistent behaviour. Various modifications to the original formulation have been proposed for transient problems (see, for example, [22] and the comprehensive review in [14]). Here, we propose a temporally consistent Rhie-Chow interpolation that can be applied to transient simulations of incompressible solid deformation described by the governing equations in the pressure-displacement formulation. The modification is implemented through an explicit correction term $\delta {\mathsf{\Omega }}_f^{\mathrm{CRC}}$ in Equation (37), which is provided for first-order accurate implicit Euler temporal discretisation schemes in Appendix C.

Using Equation (37), the discretised incopressiblity condition (35) can be transformed into a discretised pressure increment equation, which can be expressed as a linear algebraic equation for the pressure increment of any cell P:

$$\begin{array}{cc}\hfill A_P^{\delta p}·\delta p_P^{n+1}+& \sum _N{A}_N^{\delta p}·\delta p_N^{n+1}+\hfill \\ \hfill & {\mathbf{A}}_P^{\delta p,\delta \mathit{u}}·\delta {\mathit{u}}_P^{n+1}+\sum _N{\mathbf{A}}_N^{\delta p,\delta \mathit{u}}·\delta {\mathit{u}}_N^{n+1}=r_P^{\delta p},\hfill \end{array}$$

(39)

where $A_P^{\delta p}$ and $A_N^{\delta p}$ are diagonal and off-diagonal scalar coefficients related to the unknown pressure increment field, ${\mathbf{A}}_P^{\delta p,\delta \mathit{u}}$ and ${\mathbf{A}}_N^{\delta p,\delta \mathit{u}}$ are diagonal and off-diagonal vectorial coefficients related to the unknown displacement increment field and $r_P^{\delta p}$ is the right-hand side of linear equation. The central and neighbour coefficients in Equation (39) as well as the right-hand side term are defined in Appendix B.

As will be shown later, Equations (34) and (39) are used to assemble block linear system which is solved for cell-centre displacement increment $\delta {\mathit{u}}_P$ and cell-centre pressure increment $\delta p_P$. After cell-centre unknowns are obtained, Rhie-Chow interpolation formula (37) is used to calculate cell-face volume increment (volume sweep by the cell-faces), which is after that used to calculate "conservative" cell-face displacement increment as follows:

$$\delta {\mathit{u}}_f^{n+1}=[\mathit{I}-{\left(\mathit{n}\mathit{n}\right)}_f^m]·{{}^m\overline{\overline{\left(\delta {\mathit{u}}_P\right)}}}_f^{n+1}+{\displaystyle \frac{\delta {\mathsf{\Omega }}_f^{n+1}}{{\mathsf{\Gamma }}_f^m}}{\mathit{n}}_f^m.$$

(40)

In this way, the calculated displacement increment is guaranteed to satisfy the discrete incompressibility condition (35).

3.2.3. Calculation of Gradients

The very important component of the solution procedure is the calculation of gradients of unknown fields. Pressure increment gradient is required at cell-centres, while displacement increment gradient has to be calculated at the cell-centres as well as at the cell-face centres. The cell-centre gradient of all fields is calculated using the least-squares linear fit [23]. This method produces a second-order accurate gradient irrespective of local mesh quality.

The face-centre gradient of displacement increment is used to evaluate explicit part at the right hand side of Equation (22). Full cell-centre gradient is composed from normal and tangential component as follows:

$${\left({\nabla }_0\delta \mathit{u}\right)}_f={({\mathit{n}}^0{\mathit{n}}^0·{\nabla }_0\delta \mathit{u})}_f+{\left({\nabla }_{0t}\delta \mathit{u}\right)}_f,$$

(41)

where the normal derivative of displacement increment ${({\mathit{n}}^0·{\nabla }_0\delta \mathit{u})}_f$ is calculated using Equation (25), while tangential face-centre gradient is calculated by taking tangential component of interpolated cell-centre gradient:

$${\left({\nabla }_{0t}\delta \mathit{u}\right)}_f=(\mathit{I}-{\mathit{n}}_f^0{\mathit{n}}_f^0)·{{}^0\overline{\left({\nabla }_0\delta \mathit{u}\right)}}_f,$$

(42)

where ${\mathit{n}}_f^0$ is the cell-face unit normal vector. The tangential gradient at the boundary cell face is evaluated by second-order accurate linear extrapolation of the neighbouring cell-centre gradients.

3.2.4. Mesh Vertices Displacement

Before the pressure-increment equation is discretized, the mesh is moved to an intermediate configuration by updating the positions of the vertices from their reference locations to the intermediate ones. The corresponding vertex displacements are reconstructed from the cell-center displacements of the surrounding cells using a weighted least-squares method with a linear fitting function (see details in [24]).

3.2.5. Initial and Boundary Condition

To start the transient simulation, it is necessary to specify the displacement and velocity at the initial time for all cell centres in the computational mesh. Boundary conditions must be applied to the cell faces that lie on the domain boundary.

Up to this point, the discretisation procedure has been described only for internal cell faces. To discretise the governing equations at boundary cell faces, it is necessary to take into account the specified boundary conditions. This is explained below for two typical boundary conditions: specified displacement increment (Dirichlet boundary condition) and specified traction (Neumann boundary condition).

If a specified displacement increment boundary condition is used, the normal derivative of the displacement increment in Equation (23) is approximated using the scheme defined by Equation (25), where $\delta {\mathit{u}}_N$ is replaced by the specified boundary cell-face displacement increment $\delta {\mathit{u}}_b$. At the same boundary, zero normal-derivative boundary condition is applied for the pressure increment field, so that Equation (37) yields a boundary-face volume increment corresponding to the specified displacement increment. The zero normal-derivative boundary condition is also used to update the pressure increment at the corresponding boundary cell-face at the end of each outer iteration of the iterative solution procedure.

For boundary faces where specified traction increment boundary condition is applied, the traction increment term $\delta {\mathit{T}}_f$ in Equation (22) is added to the source term (right hand side), and traction correction term $\delta {\mathit{T}}_f^{\prime }$ is set to zero. At the same boundary faces, the specified pressure increment boundary conditions is applied, where pressure increment value is calculated from the constitutive equation using gradient of displacement increment at the boundary from the previous iteration. The cell-face displacement increment term in Equation (37) is defined implicitly at the boundary cell-face b as follows:

$${{}^m\overline{\overline{\left(\delta \mathit{u}\right)}}}_b^{n+1}=\delta {\mathit{u}}_P^{n+1}+{\left({\displaystyle \frac{\partial \delta \mathit{u}}{\partial n}}\right)}_b^{*}{\delta }_{bn}^m,$$

(43)

where $\delta {\mathit{u}}_P^{n+1}$ is the displacement increment at the centre of the cell next to the considered boundary face b, ${\delta }_{bn}^m$ is the normal distance between boundary face centre b and adjacent cell-centre P, the normal derivative of the displacement increment is defined as follows:

$${\left({\displaystyle \frac{\partial \delta \mathit{u}}{\partial n}}\right)}_b^{*}={\displaystyle \frac{\delta {\mathit{u}}_b^{*}-\delta {\mathit{u}}_P^{*}}{{\delta }_{bn}^m}}.$$

(44)

The superscript * denotes quantities obtained in the previous iteration, and the boundary face displacement increment $\delta {\mathit{u}}_b^{*}$ is defined by Equation (40). In this equation, the tangential component of the displacement increment (first term on the right-hand side) is computed to satisfy the constitutive relation at traction-specified boundary faces.

3.2.6. Final Block-Coupled System of Linear Algebraic Equations

By combining the discretised momentum equation (34) with the discretised pressure-increment (incompressibility condition) equation (39), one obtains the following coupled system of algebraic equations for each cell P in the computational mesh:

$$\begin{array}{c}\hfill {\mathbf{A}}_P·{\mathbf{\Phi }}_P+\sum _N{\mathbf{A}}_N·{\mathbf{\Phi }}_N={\mathbf{R}}_P,\end{array}$$

(45)

where ${\mathbf{A}}_P$ and ${\mathbf{A}}_N$ are the $4\times 4$ dense matrices representing the diagonal and off-diagonal coefficients of the coupled linear system, obtained by combining the corresponding coefficients of the discretised momentum equation (34) and the discretised pressure equation (39), as is ilustrated for diagonal coefficient ${\mathbf{A}}_P$ in Figure 2. For each cell P, it is defined the $4\times 1$ vector $\mathbf{\Phi }$, which consists of the cell-centre unknown fields: the three components of the cell-centre displacement increment vector and the cell-centre pressure increment. For example, for cell P this vector of unknowns is given by ${\mathbf{\Phi }}_P={\left[{\left(\delta {\mathit{u}}_P\right)}_x{\left(\delta {\mathit{u}}_P\right)}_y{\left(\delta {\mathit{u}}_P\right)}_z\delta p_P\right]}^{\mathrm{T}}$. The $4\times 1$ vector ${\mathbf{R}}_P$ on the right-hand side of Equation (45) is composed of the corresponding right-hand side vectors of the linear systems (34) and (39), as follows: ${\mathbf{R}}_P={\left[{\left({\mathbf{R}}_P^{\delta \mathit{u}}\right)}^x{\left({\mathbf{R}}_P^{\delta \mathit{u}}\right)}^y{\left({\mathbf{R}}_P^{\delta \mathit{u}}\right)}^z{r}_P^{\delta p}\right]}^{\mathrm{T}}$. After composing the system of linear algebraic equations (45) for all cells in the computational mesh, the corresponding global linear system is obtained:

$$\mathbf{A}·\mathbf{\Phi }=\mathbf{R},$$

(46)

where $\mathbf{A}$ is the square block matrix, $\mathbf{\Phi }$ is the block vector of unknown fields and $\mathbf{R}$ is the block vector of the right-hand side, and their dimension is equal to the number of cells in the computational mesh. The block linear system (46) is solved using the Block-Jacobi preconditioner and the BiSGSTAB solver.

3.3. Solution Procedure

The overall solution procedure is summarized in Algorithm 1 and can be categorized as the fixed-point/Picard iteration method. The main part of the solution procedure is the outer iteration loop that starts with the assembly of the global linear system (46), where the elements of the vector $\mathbf{R}$ in the current iteration (k) are calculated using the values of unknown fields obtained in the previous outer iteration ($k-1$). The outer iteration loop repeats until the norm of the initial residual of the system (46) is greater than the predefined tolerance: $|{\mathbf{A}}^{k-1}·{\mathbf{\Phi }}^{k-1}-{\mathbf{R}}^{k-1}|>ϵ$.

Algorithm 1 Block coupled solution procedure.

while$t^{n+1}<t_{\mathrm{end}}$do     Initialize fields: $\delta p_{\mathrm{P}}^{n+1}$ and $\delta {\mathit{u}}_{\mathrm{P}}^{n+1}$     $k=0$     while $|{\mathbf{A}}^{k-1}·{\mathbf{\Phi }}^{k-1}-{\mathbf{R}}^{k-1}|>ϵ$ do         $k=k+1$         Calculate cell displacement increment gradient ${(\nabla \delta \mathit{u})}_P^{n+1}$         Calculate face displacement increment gradient ${(\nabla \delta \mathit{u})}_f^{n+1}$         Update total displacement gradients and corresponding deformation gradients         Update ${\sigma }_P=\sigma ({\mathit{F}}_P^{k-1},p_P^{k-1})$ and ${\sigma }_f=\sigma ({\mathit{F}}_f^{k-1},p_f^{k-1})$ using …         … selected constitutive equation         Assemble linear system (34)         Move mesh to "intermediate" configuration         Assemble linear system (39)         Assemble global linear system ${\mathbf{A}}^{k-1}·{\mathbf{\Phi }}^k={\mathbf{R}}^{k-1}$ and solve it for $\delta {\mathit{u}}_P^{n+1}$ and $\delta p_P^{n+1}$         Calculate cell-face sweep volume using Rhie-Chow interpolation (37)         Calculate conservative cell-face displacement increment $\delta {\mathit{u}}_f^{n+1}$ using Equation (39)         Move mesh to initial configuration         Update total displacement and total pressure:         ${\mathit{u}}_P^{n+1}={\mathit{u}}_P^n+\delta {\mathit{u}}_P^{n+1}$         ${\mathit{u}}_f^{n+1}={\mathit{u}}_f^n+\delta {\mathit{u}}_f^{n+1}$         $p_P^{n+1}=p_P^n+\delta p_P^{n+1}$     end while     Calculate cell-centre velocity ${\mathit{v}}_P$ using selected temporal scheme     $t^{n+1}=t^n+\Delta t$ end while

4. Results

In this section, seven test cases are presented to demonstrate the solver’s performance across different deformation regimes. The infinite plate with a circular hole is first used to assess accuracy in the linear elastic regime. Next, uniaxial stretch and simple shear tests highlight the nonlinear capabilities of the solver under elementary deformation modes. Inflation of a pressurized cylinder illustrates both steady-state and transient behaviour, including temporal accuracy analysis for Euler, backward, trapezoidal, and composite discretisation schemes. The heart tissue beam test evaluates the solver’s accuracy and stability in transient problems involving large rotations and displacements. Solver performance with the HGO constitutive model is demonstrated through inflation of a rat carotid artery. Finally, ventricle inflation serves as a realistic three-dimensional benchmark. In all cases, results obtained with the FV method are compared with analytical or FE solutions obtained using FEBio software suit [25,26].

4.1. Infinite Plate with a Circular Hole Subjected to Uniform Tension

A plate with a circular hole in its centre is loaded with uniform tension in the axial direction (see Figure 3). The problem is considered as the plane strain. If the radius of the hole is small compared to the dimensions of the plate, then the analytical solution obtained for an infinitely large plate describes the stress and displacement distribution around the hole [27]:

$$\begin{array}{cc}\hfill {\sigma }_{xx}(r,\theta )& =t_x\left[1-{\displaystyle \frac{a^2}{r^2}}\left({\displaystyle \frac{3}{2}}cos2\theta +cos4\theta \right)+{\displaystyle \frac{3}{2}}{\displaystyle \frac{a^4}{r^4}}cos4\theta \right],\hfill \\ \hfill {\sigma }_{yy}(r,\theta )& =t_x\left[-{\displaystyle \frac{a^2}{r^2}}\left({\displaystyle \frac{1}{2}}cos2\theta -cos4\theta \right)-{\displaystyle \frac{3}{2}}{\displaystyle \frac{a^4}{r^4}}cos4\theta \right],\hfill \\ \hfill {\sigma }_{xy}(r,\theta )& =t_x\left[-{\displaystyle \frac{a^2}{r^2}}\left({\displaystyle \frac{1}{2}}sin2\theta +sin4\theta \right)+{\displaystyle \frac{3}{2}}{\displaystyle \frac{a^4}{r^4}}sin4\theta \right],\hfill \end{array}$$

$$\begin{array}{cc}\hfill u_x(r,\theta )& ={\displaystyle \frac{a{t}_x}{8\mu }}\left\{{\displaystyle \frac{r}{a}}(k+1)cos\theta +2{\displaystyle \frac{a}{r}}\left[(k+1)cos\theta +cos3\theta \right]-2{\displaystyle \frac{a^3}{r^3}}cos3\theta \right\},\hfill \\ \hfill u_y(r,\theta )& ={\displaystyle \frac{a{t}_x}{8\mu }}\left\{{\displaystyle \frac{r}{a}}(k-3)sin\theta +2{\displaystyle \frac{a}{r}}\left[(1-k)sin\theta +sin3\theta \right]-2{\displaystyle \frac{a^3}{r^3}}sin3\theta \right\},\hfill \end{array}$$

where r and $\theta$ are polar coordinates, a is hole radius, $\mu$ is shear modulus, and $k=3-4\nu$ is the Kolosov constant for plane strain deformation.

The plate material is assumed to be linear elastic, isotropic, and incompressible (Poisson’s ratio $\nu =0.5$), with Young’s modulus $E={10}^9\mathrm{Pa}$. Owing to the problem’s symmetry, only a quarter of the plate is modelled, as shown in Figure 4 (shaded area). Following earlier studies of the same problem, e.g., [19], the exact solution corresponding to $t_x=10\mathrm{kPa}$ is prescribed on all traction-specified boundaries (top and right dashed edges) to eliminate finite-domain effects.

Spatial accuracy is demonstrated by comparing results obtained using the proposed FV solver with the corresponding analytical solution. This comparison was first conducted for three meshes of similar resolution, each consisting of a different cell shape: quadrilateral, triangular, and polygonal7, as shown in Figure 4. Figure 5 presents a comparison of the x-component of the displacement vector and the $xx$-component of the stress tensor with analytical results along the $AB$ boundary of the plate. Similarly, Figure 6 shows a comparison of the y-component of the displacement vector and the $yy$-component of the stress tensor along the $DE$ boundary. Very good agreement is observed between the analytical solution and the FV solver results for all three cell shapes.

The order of spatial accuracy of the proposed FV solver is evaluated by performing simulations on four quadrilateral meshes with successively increasing numbers of cells (1000, 4000, 16,000, 64,000), with the coarsest mesh shown in Figure 4(a). The cell-centre errors are calculated for displacement magnitude and $xx$-component of the Cauchy stress in respect to the available analytical solution and magnitude of the errors are computed using standard norms ($L^1$, $L^2$, and $L^{\infty }$). Figure 7 presents the error norms as a function of the average cell size. Second-order accuracy is achieved for the displacement across all three norms, whereas for the $xx$-component of the Cauchy stress tensor, first-order accuracy is observed with respect to the $L^{\infty }$ norm, while second-order accuracy is attained for the remaining two norms.

4.2. Uniaxial Extension and Simple Shear Tests

To test the solvers’ nonlinear capabilities for fundamental types of deformation, an uniaxial extension test and a simple shear test are performed on a cube with dimensions $1\mathrm{m}\times 1\mathrm{m}\times 1\mathrm{m}$. The cube is modeled as an incompressible material ($\nu =0.5$). Both tests are solved using the Neo-Hookean and Guccione material laws. For the Neo-Hookean material, the shear modulus of the cube is $\mu ={10}^6\mathrm{Pa}$. For the Guccione law, the material properties are $C=2\times {10}^6\mathrm{Pa}$, with $b_f=b_t=b_{fs}=1$. A detailed description of the incompressible Guccione constitutive relation is provided in Appendix D.2.

4.2.1. Uniaxial Extension Test

The simulation of uniaxial extension along the x-direction is performed on one quarter of the unit cube by applying symmetry boundary conditions on the $xz$- and $xy$-planes (see Figure 8). The $yz$-plane is also modeled as a symmetry plane to allow free contraction of the cube in the lateral directions. On the remaining boundaries, specified traction boundary conditions are applied, with the traction calculated from the available analytical solution. For uniaxial extension along the x-direction, the deformation gradient can be expressed in terms of the principal stretches as follows:

$$\mathit{F}=\left[\begin{array}{ccc}{\lambda }_x& 0& 0\\ 0& {\lambda }_y& 0\\ 0& 0& {\lambda }_z\end{array}\right],$$

(47)

where ${\lambda }_y={\lambda }_z=1/\sqrt{{\lambda }_x}$ for the case of an incompressible material. For a given deformation gradient, the Cauchy stress and the corresponding traction on the cube boundaries is calculated based on the selected constitutive relation. Simulations are performed for the principal stretch ${\lambda }_x$ varying in the range 1–$1.8$.

Figure 9 shows the deformed configuration of the cube for ${\lambda }_x=1.8$, calculated using the incompressible Neo-Hookean material. For the same material, Figure 10 presents a comparison between the numerical and analytical solutions for the stretches ${\lambda }_y={\lambda }_z$, the hydrostatic pressure, and the $xx$-component of the Cauchy stress. The numerical results closely follow the analytical solution. An analogous comparison for the incompressible Guccione material is shown in Figure 11.

4.2.2. Simple Shear Test

The simple shear test involves a deformation in which parallel planes of the material slide past each other. The deformation gradient for simple shear in the $xy$-plane is defined as

$$\mathit{F}=\left[\begin{array}{ccc}1& \gamma & 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right],$$

(48)

where $\gamma$ is the shear strain. In this study, a simple shear deformation is applied to the unit cube shown in Figure 12. A zero-displacement boundary condition is imposed on the bottom surface, while fixed traction boundary conditions are applied on the remaining surfaces, with the traction calculated using the selected constitutive relation (Neo-Hookean or Guccione) and the deformation gradient defined in Equation 48. The shear strain is varied in the range 0–$0.6$.

The reference and final (deformed) configurations of the cube are shown in Figure 13 for the Neo-Hookean material model and for $\gamma =0.6$. The current configuration is coloured by the displacement field, where one can note parallel and horizontal contour lines as is expected for the case considered. Figure 14 presents a comparison between the exact and numerically calculated displacement of the top surface of the cube for the Neo-Hookean and Guccione material models as a function of the specified shear strain $\gamma$, where the exact value of the displacement for the unit cube is $u_{x,a}=\gamma \times 1\mathrm{m}$.

4.3. Inflation of a Thick-Wall Cylinder

An infinitely long thick-walled cylinder with inner radius $r_i=8\mathrm{m}$ and outer radius $r_o=16\mathrm{m}$, is subjected to pressure $p=10\mathrm{MPa}$. This case was also solved in [4]. A plane strain deformation is assumed. The material of the cylinder is assumed to be incompressible $(\nu =0.5)$ and hyperelastic with Young’s modulus $E=1000\mathrm{MPa}$. The behaviour of the incompressible hyperelastic material is described by the Neo-Hookean constitutive relation.

Exploiting the geometric and loading symmetry, only one quarter of the domain is modeled, as shown in Figure 15, with symmetry boundary conditions applied on the corresponding planes. The outer wall is modeled as traction-free. The analytical solution for this problem is provided in Appendix A.

Figure 16 shows a hexahedral mesh with 65 cells and a polyhedral mesh with 1,338 cells. These meshes, together with two finer hexahedral meshes consisting of 260 and 1,040 cells, were used to compute the solution for an internal pressure of $p=100\mathrm{MPa}$. The predicted radial displacement, radial stress, hoop stress, and pressure along the radius are compared to the analytical solution in Figure 17. The solutions obtained on the hexahedral meshes converge consistently to the analytical solution with increasing mesh resolution, and good agreement is also observed for the solution obtained on the polyhedral mesh. Temporal accuracy analysis, presented in continuation, is performed for results obtained on the finer hexahedral mesh.

To assess temporal accuracy, the case was solved using different time-step sizes: ${10}^{-2}\mathrm{s}$, $5\times {10}^{-3}\mathrm{s}$, $2.5\times {10}^{-3}\mathrm{s}$, $1.25\times {10}^{-3}\mathrm{s}$, $6.25\times {10}^{-4}\mathrm{s}$ and ${10}^{-5}\mathrm{s}$. The prescribed pressure at the inner surface is increased linearly from $p=0\mathrm{MPa}$ to $p=100\mathrm{MPa}$ over a time period $0.01\mathrm{s}$. The accuracy analysis is carried out at the time instance $t=0.01\mathrm{s}$. The cell-centre displacement error is computed for solution obtained with the five largest time-step sizes with respect to the reference solution obtained with the smallest time-step size ${10}^{-5}\mathrm{s}$. The magnitude of the cell-centre displacement error was quantified using the $L^{\infty }$ norm.

The temporal accuracy is assessed for both linear and non-linear deformation regimes employing four temporal discretisation schemes: Euler, backward, Crank–Nicolson and composite. The corresponding displacement error norms are presented in Figure 18 as functions of the time-step size. For both deformation regimes, the Euler scheme exhibits the expected first-order convergence of the error, whereas the remaining schemes attain the expected second-order convergence.

Finally, the long-term temporal response of the cylinder was analysed for the case in which the internal pressure was first increased linearly from $p=0\mathrm{MPa}$ to $p=100\mathrm{MPa}$ over the interval $t=0\mathrm{s}$ to $t=0.01\mathrm{s}$, and then decreased linearly back to $p=0\mathrm{MPa}$ over the interval $t=0.01\mathrm{s}$ to $t=0.02\mathrm{s}$. For this scenario, the total energy of the isolated system must remain constant for $t>0.02\mathrm{s}$, in the absence of external forcing. Figure 19 shows temporal response of kinetic, potential and total energy obtained using backward temporal discretisation scheme, where the total energy is computed as a sum of kinetic and potential energy and potential energy is calculated according to Neo-Hookean constitutive relation. One can note, expected behaviour of the total energy after time instance $t=0.02\mathrm{s}$. Other second-order temporal discretisation schemes shows similar behaviour.

4.4. Heart Tissue Beam

In this benchmark test case proposed by Land et. al [28], one considers deformation of a thick beam of $10\mathrm{mm}$ length, and $(1\mathrm{mm}\times 1\mathrm{mm})$ quadratic cross section (see Figure 21). The test case is analysed for two hyperelastic constitutive models: incompressible Guccione model as originally proposed in [28] (fibre direction along long axis, constitutive parameters: $C=2000\mathrm{Pa}$, $b_f=8$, $b_t=2$, $b_{fs}=4$) and incompressible Neo-Hookean model with shear modulus $\mu =15000\mathrm{Pa}$. The properties of the incompressible Guccione model approximately correspond to the heart tissue properties.

In first step, a mesh sensitivity study is performed where steady-state FV solution is compared with the corresponding solution obtained using FEBio software suit [25,26], verified on selected test cases proposed in [28]. The beam is discretised by a uniform structured mesh, and it is subjected to a uniform pressure load of $4\mathrm{Pa}$ on the bottom side, while the left side is fixed in all directions. Figure 21 shows the z-coordinate of the tip point A (defined in Figure 20) as a function of cell dimension, where one can note consistent convergence toward the FE benchmark solution. Figure 21 shows the maximal deflection of the beam. Results are given for a $50\times 5\times 5$ computational mesh.

A temporal accuracy study was performed with the Neo-Hookean model, following the procedure outlined in Section 4.3. The results confirm first- and second-order accuracy for the Euler and backward time discretization schemes, respectively, as illustrated in Figure 22, for both linear and nonlinear models. Results obtained with the trapezoidal scheme are not reported, as numerical instabilities were observed for the smallest time-step sizes. To address this issue in nonlinear cases, the composite scheme proposed by Bathe and Noh [21] was employed. For sufficiently small time-step sizes, the composite scheme proposed by Bathe and Noh [21] attains second-order accuracy for both linear and nonlinear model.

Figure 20. Heart tissue beam: reference configuration, with point A indicated.

Figure 20. Heart tissue beam: reference configuration, with point A indicated.

Figure 21. Heart tissue beam: the computational mesh ($50\times 5\times 5$) in deformed configuration. Colors represent z-component of the displacement vector obtained for the Neo-Hookean material law.

Figure 21. Heart tissue beam: the computational mesh ($50\times 5\times 5$) in deformed configuration. Colors represent z-component of the displacement vector obtained for the Neo-Hookean material law.

Figure 22. Heart tissue beam: Time convergence study of the displacements for Euler and backward temporal schemes for (a) non-linear and (b) linear solver.

Figure 22. Heart tissue beam: Time convergence study of the displacements for Euler and backward temporal schemes for (a) non-linear and (b) linear solver.

To check the stability of the solver for the long-term temporal response of the beam. The pressure, which is applied to the bottom side, is linearly increased to the maximum value of $4\mathrm{Pa}$ until $0.02\mathrm{s}$, and then it is decreased back to $0\mathrm{Pa}$ until $0.04\mathrm{s}$. Figure 24 and Figure 25 show the displacement and velocity of the tip of the beam, along with the kinetic, elastic, and total energy of the system, calculated using the Euler, backward, composite, and trapezoidal temporal schemes. Euler discretization scheme does not conserve the energy of the isolated system. In contrast, total energy remains constant after the force is no longer applied when using backward, Crank-Nicolson, and composite temporal discretisation schemes.

Figure 23. Heart tissue beam: Error in maximal z-deflection at beam tip for different mesh densities for Neo-Hookean material law. Results converge to a consensus solution as the number of volumes increases.

Figure 23. Heart tissue beam: Error in maximal z-deflection at beam tip for different mesh densities for Neo-Hookean material law. Results converge to a consensus solution as the number of volumes increases.

Figure 24. Heart tissue beam: Displacements, velocity, and energy of the system with Euler temporal discretization.

Figure 24. Heart tissue beam: Displacements, velocity, and energy of the system with Euler temporal discretization.

Figure 25. Heart tissue beam: Displacements, velocity, and energy of the system with backward temporal discretization.

Figure 25. Heart tissue beam: Displacements, velocity, and energy of the system with backward temporal discretization.

Figure 26. Heart tissue beam: Displacements, velocity, and energy of the system with composite temporal discretization.

Figure 26. Heart tissue beam: Displacements, velocity, and energy of the system with composite temporal discretization.

Figure 27. Heart tissue beam: Displacements, velocity, and energy of the system with trapezoidal temporal discretization.

Figure 27. Heart tissue beam: Displacements, velocity, and energy of the system with trapezoidal temporal discretization.

4.5. Inflation of a Rat Carotid Artery

This test case is based on the experimental study reported in [29], where carotid arteries were harvested from rats and subjected to internal pressurization, during which the outer radius was measured as a function of the applied pressure. Based on these experimental data, the HGO material model [15] was calibrated in [30]. This constitutive model effectively captures the anisotropic mechanical response of arterial walls by decomposing the strain energy density function into two components: an isotropic part, dominant at low pressures and consistent with a Neo-Hookean material behavior, and an anisotropic part, which becomes prominent at higher pressures which models the effect of progressive straightening of helically oriented collagen fibres. The material parameters were defined as follows [30]: Young’s modulus $E=132.69\mathrm{kPa}$, $k_1=206\mathrm{Pa}$, and $k_2=1.465$. The collagen fibers were assumed to be symmetrically oriented at $\pm 39.{76}^{\circ }$ with respect to the circumferential direction.

In this study, numerical simulations of the rat carotid artery inflation, as described above, are performed following a similar FE study by Sun et al. [31]. The computational domain is defined as a cylindrical segment with an inner radius of $r_i=3\mathrm{mm}$, an outer radius of $r_o=4\mathrm{mm}$, and a length of $5\mathrm{mm}$. Axial displacements are constrained. Owing to geometric and loading symmetry, only one quarter of the domain is modeled, as illustrated in Figure 28, with appropriate symmetry boundary conditions applied on the corresponding planes. The outer wall is assumed to be traction-free, while an internal pressure of $p=25\mathrm{kPa}$ is applied to the inner wall through 100 loading steps.

Figure 28 also shows the deformed geometry at the final loading step. The evolution of the inner and outer radii of the cylindrical surfaces as a function of the applied inner pressure, obtained with the FV solver, is presented in Figure 29. Results for the inner surface are compared with experimental data from [29] and FE results from [31], while results for the outer surface are compared against FE results only. Overall, good agreement is observed.

4.6. Inflation of an Idealised Ventricle

Inflation of an idealised ventricle (Figure 30) was proposed by Land et al. [28] as a benchmark problem for cardiac mechanics software. The case is three-dimensional and static, with finite hyperelastic deformation. The initial geometry is defined as a truncated ellipsoid:

$$x=r_{s}sin\left(u\right)cos\left(v\right),y=r_{s}sin\left(u\right)sin\left(v\right),z=r_{l}cos\left(u\right).$$

(49)

The endocardial (inner) surface is defined as $r_s=7\mathrm{mm}$, $r_l=17\mathrm{mm}$, $u\in \left[-\pi ,-arccos{\displaystyle \frac{5}{17}}\right]$, and $v\in [-\pi ,\pi ]$. The epicardial (outer) surface is defined as $r_s=10\mathrm{mm}$, $r_l=20\mathrm{mm}$, $u\in \left[-\pi ,-arccos{\displaystyle \frac{5}{20}}\right]$, and $v\in [-\pi ,\pi ]$. The implicitly defined base plane is positioned at $z=5\mathrm{mm}$. The behavior of hyperelastic material is described by the transversely isotropic incompressible Guccione model [16], where the parameters are $C=10\mathrm{kPa}$, $b_f=b_t=b_{fs}=1$ (the specified parameters produce isotropic behavior).

Although the geometry is axisymmetric, one quarter of the domain is modeled, as shown in Figure 30. The geometry is discretized entirely with hexahedral prisms. The epicardial surface is modelled as traction-free, while a fixed displacement boundary condition is prescribed on the top base plane. Symmetry boundary conditions are applied to the two cutting planes used to extract the quarter section of the ventricle. The endocardial surface is subjected to an internal pressure of $2300\mathrm{Pa}$. The internal pressure value of $10000\mathrm{Pa}$ proposed in [28] could not be reached with the proposed FV solver due to the occurrence of unphysical pressure oscillations.

The load is applied incrementally in steps of $100\mathrm{Pa}$. Simulations are performed on three successively refined meshes containing 2728, 16992, and 38292 cells. The reference solution is obtained using the FEBio software on the finest mesh. Figure 31 shows the magnitude of the apex displacement at the endocardial and epicardial surfaces as a function of the pressure load applied to the endocardial surface. The solution obtained using the proposed FV solver consistently converges with mesh refinement toward the reference FE solution. Figure 32 shows a comparison of the ventricle shape at loading pressure of $2300\mathrm{Pa}$. Figure 32a shows ventricle wall midline in its undeformed and deformed configurations obtained using FV and FEBio solvers. Similarly, Figure 32b shows undeformed and deformed configuration for a line on the endocardial surface. In both cases one can observe very good agreement between FV and FE results.

The simulations described above (performed using incompressible Guccione constitutive model) were repeated with the incompressible Neo-Hookean constitutive model, with the shear modulus set to $\mu =4000\mathrm{kPa}$, corresponding to the initial shear modulus of the Guccione model. The maximum achievable pressure load on the endocardial surface was $1700\mathrm{Pa}$ for both the proposed FV solver and the FE solver. The remaining boundary conditions were identical to those described previously.

Figure 33 illustrates the apex displacement magnitude at the endocardial and epicardial surfaces as a function of the applied internal pressure load, showing convergence toward the reference results with mesh refinement. Figure 34 compares the overall ventricle shape at an internal pressure of $1700\mathrm{Pa}$, as computed with the proposed FV solver and the FEBio tool. Figure 34a presents the deformed configuration of the ventricle wall midline, while Figure 34b examines the deformed shape of a line on the endocardial surface. In both cases, the FV and FE results demonstrate very good agreement.

Appendix D.1. HGO Model

The Holzapfel–Gasser–Ogden (HGO) constitutive model is widely used to describe the nonlinear mechanical response of arterial tissue. In this implementation, the model is formulated with three material parameters and uses two families of collagen fibers. The fibers are assumed to be arranged symmetrically around the circumferential axis at angles $\pm \phi$. They reinforce the isotropic matrix and make the tissue stiffer and highly anisotropic.

The strain energy density function is defined as

$$\mathsf{\Psi }={\mathsf{\Psi }}_{\mathrm{iso}}+{\mathsf{\Psi }}_{\mathrm{aniso}}={\displaystyle \frac{\mu }{2}}(I_1-3)+{\displaystyle \frac{k_1}{2k_2}}\sum _{i=4,6}\left(exp\left[k_2{(I_i-1)}^2\right]-1\right),$$

(A12)

where $\mu$ is the shear modulus of the isotropic ground matrix, $k_1$ and $k_2$ are material parameters associated with the collagen fibers, and $I_1$, $I_4$, and $I_6$ are strain invariants.

The invariants are computed as

$$\begin{array}{cccc}\hfill I_1& =\mathrm{tr}\left(\mathbf{C}\right),I_4\hfill & \hfill =\mathbf{C}:({\mathbf{a}}_0\otimes {\mathbf{a}}_0),I_6& =\mathbf{C}:({\mathbf{b}}_0\otimes {\mathbf{b}}_0),\hfill \end{array}$$

(A13)

where $\mathbf{C}={\mathbf{F}}^{\mathrm{T}}\mathbf{F}$ is the right Cauchy–Green deformation tensor, and ${\mathbf{a}}_0$, ${\mathit{b}}_0$ are unit vectors in the reference configuration defining the preferred fiber directions.

The corresponding Cauchy stress tensor is obtained as:

$$\begin{array}{cc}\hfill \sigma & =-p\mathbf{I}+\mu (\mathbf{B}-\mathbf{I})\hfill \\ \hfill & +2k_1(I_4-1)exp\left[k_2{(I_4-1)}^2\right](\mathbf{a}\otimes \mathbf{a})\hfill \\ \hfill & +2k_1(I_6-1)exp\left[k_2{(I_6-1)}^2\right](\mathbf{b}\otimes \mathbf{b}).\hfill \end{array}$$

(A14)

Appendix D.2. Guccione Model

To describe the mechanical behavior of cardiac tissue, the Guccione constitutive model was implemented. This model takes into account the presence of collagen fibers, which reinforce the tissue matrix and contribute to its stiffness and strength. As a result, the tissue exhibits anisotropic behavior, specifically transverse isotropy, due to the directional alignment of fibers.

The strain energy density function $\mathsf{\Psi }$ is defined as:

$$\mathsf{\Psi }={\displaystyle \frac{k}{2}}(e^Q-1)$$

(A15)

where k is a material constant and Q is a scalar function representing the anisotropic response of the tissue, given by:

$$Q=b_t{I}_1^2-2b_t{I}_2+(b_f-2b_{fs}+b_t)I_4^2+2(b_{fs}-b_t)I_5$$

(A16)

Here, $b_t$, $b_f$, and $b_{fs}$ are material parameters, and $I_1$, $I_2$, $I_4$, and $I_5$ are invariants of the Green-Lagrange strain tensor $\mathbf{C}$. The Green-Lagrange strain is computed from the deformation gradient $\mathbf{F}$ as:

$$\mathbf{C}={\displaystyle \frac{1}{2}}({\mathbf{F}}^T·\mathbf{F}-\mathbf{I})$$

(A17)

with the deformation gradient defined as:

$$\mathbf{F}=\mathbf{I}+{\displaystyle \frac{\partial \mathbf{u}}{\partial \mathbf{X}}}$$

(A18)

Where $\mathbf{u}$ is the displacement vector and $\mathbf{X}$ is the reference configuration. The strain invariants are computed as follows:

$$\begin{array}{c}\hfill I_1=\mathrm{tr}\left(\mathbf{C}\right),I_2={\displaystyle \frac{1}{2}}\left[{\left(\mathrm{tr}\left(\mathbf{C}\right)\right)}^2-\mathrm{tr}(\mathbf{C}·\mathbf{C})\right],I_4=\mathbf{C}:({\mathbf{f}}_0\otimes {\mathbf{f}}_0),I_5={\mathbf{C}}^{\mathbf{2}}:({\mathbf{f}}_0\otimes {\mathbf{f}}_0),\end{array}$$

(A19)

where $f_0$ is a unit vector of reference fiber direction.

The second Piola–Kirchhoff stress tensor $\mathbf{S}$ is then obtained by differentiating the strain energy density concerning $\mathbf{C}$:

$$S={\displaystyle \frac{\partial \mathsf{\Psi }}{\partial \mathbf{C}}}={\displaystyle \frac{\partial Q}{\partial \mathbf{C}}}{\displaystyle \frac{k}{2}}{e}^Q$$

(A20)

The derivative of Q with respect to the Green-Lagrange strain ${\displaystyle \frac{\partial Q}{\partial \mathbf{C}}}$ is:

$${\displaystyle \frac{\partial Q}{\partial \mathbf{C}}}=2c_t\mathbf{C}+2(b_f-2b_{fs}+b_t)I_4(f_0\otimes f_0)+2(b_{fs}-b_t)(\mathbf{C}(·f_0\otimes f_0)+(f_0\otimes f_0)·\mathbf{C}))$$

(A21)