General Finite-Element Framework of the Virtual Fields Method in Nonlinear Elasticity

Abstract

This paper presents a method to derive the virtual fields for identifying constitutive model parameters using the Virtual Fields Method (VFM). The VFM is an approach to identify unknown constitutive parameters using deformation fields measured across a given volume of interest. The general principle for solving identification problems with the VFM is first to derive parametric stress field, where the stress components at any point depend on the unknown constitutive parameters, across the volume of interest from the measured deformation fields. Applying the principle of virtual work to the parametric stress fields, one can write scalar equations of the unknown parameters and solve the obtained system of equations to deduce the values of unknown parameters. However, no rules have been proposed to select the virtual fields in identification problems related to nonlinear elasticity and there are multiple strategies possible that can yield different results. In this work, we propose a systematic, robust and automatic approach to reconstruct the systems of scalar equations with the VFM. This approach is well suited to finite-element implementation and can be applied to any problem provided that full-field deformation data are available across a volume of interest. We also successfully demonstrate the feasibility of the novel approach by multiple numerical examples. Potential applications of the proposed approach are numerous in biomedical engineering where imaging techniques are commonly used to observe soft tissues and where alterations of material properties are markers of diseased states.

1 Introduction

One of the objectives in computational biomechanics is to make predictions: given for instance a complete and patient-specific reconstruction of the aorta and surrounding tissues, we can predict the deformations induced by a deployed stent graft during surgical repair of an aneurysm [1–3]. More generally, we could predict the deformations of the aorta under the action of any radial forces applied in an experiment. This problem of predicting the result of measurements is called the simulation problem, or the forward problem. The inverse problem consists of using the experimental measurements, such as the displacement field, to infer the values of the parameters that characterize the system, such as unknown material parameters [4–11], unknown boundary conditions, or even sometimes the unknown initial geometry of the solid before the application of mechanical loading (load-free configuration in finite deformations [12, 13]).

A subcategory of inverse problems is made by identification problems [6, 14], where a finite number of unknown material parameters have to be recovered from experimental measurements. Such inverse problems can be solved by defining a cost function that estimates the absolute error between the model predictions and the measurements. The cost function is minimized in the least-squares sense. In general situations, the model is solved numerically using a finite-element model updating technique (FEMU) [6, 14–19].

In the more and more common situations when full-field measurements are available [7, 20–25], an alternative to FEMU is possible: the Virtual Fields Method (VFM), which has been shown to be more direct and robust with respect to unknown boundary conditions in these situations [26].

The general principle for solving identification problems with the VFM is first to derive parametric stress field across the volume of interest from the measured deformation field, stress components at any position for a given deformation depending on a number of unknown constitutive parameters. Applying the principle of virtual work to the parametric stress field gives scalar equations for the unknown parameters. The system of equations can be solved for the values of unknown parameters.

There are currently multiple strategies possible for the selection of virtual fields when the VFM is applied in nonlinear elasticity. A simple and efficient approach, well suited for tensile tests, is to define virtual displacements varying linearly from one boundary to the other. This approach was recently applied successfully in biaxial tension by Kazerooni et al [27] to identify the parameters of a Holzapfel model [28] for the skin. For other types of loading, for instance, Zhang et al. [23] introduced analytical expressions of virtual fields which were written in the cylindrical coordinate system using the $arctan$ function. For inflation experiments, Bersi et al. [29] used virtual fields that permitted to derive local equilibrium equations relating the intraluminal pressure and the wall tensile stresses. These virtual fields were valid for the very specific problem of inflation of thin-wall cylinder-like structures. The extension to thick-wall was solved in further papers [25, 30] by weighting locally the virtual fields with a Gaussian function, which was successfully centered at the middle of square patches in order to identify the distribution of material properties across the whole cylinder.

Although the approaches discussed above were all successful in solving specific identification problems, a general and unified VFM approach is still needed for the following reasons:

• 1. the expression of virtual fields used in most of previous studies [27, 30, 31] was well adapted to the 2D geometries such as membranes or shells but there is a range of samples or materials in which boundary conditions or geometries or both cannot be straightforwardly designed or chosen to be 2D. With the acquisition of fully volumetric deformation fields which is expanding for soft tissues [23–25, 32–35], there is a need for a method that can provide virtual fields for any 3D geometry and virtual field.

• 2. In previous studies cited above, the VFM was used to identify the parameters of incompressible hyperelastic materials. The compressibility modulus was usually disregarded. To identify simultaneously the compressibility and shear moduli of a hyperelastic material, two independent virtual fields are needed in order to separate the hydrostatic and deviatoric contributions of the stress. This problem was solved a decade ago for linear elasticity [36, 37] but it remained an open question for nonlinear elasticity.

In this paper, we introduce an original approach addressing these essential questions. An application of this method to identifying material properties of the lamina cribrosa (LC) in the optic nerve head (ONH), using optical coherence tomography (OCT) imaging data, demonstrates the viability of the technique. The paper is organized as follows: We first elaborate the theoretical aspects of the proposed VFM approach in Sect. 2. Subsequently, several numerical examples and an application are presented to test the feasibility of the method in Sect. 3. We then discuss the results and the proposed approach in Sect. 4 and end with a conclusion in Sect. 5.

2 Materials and Methods

2.1 Definition of the Problem

A solid-like body, or part of it, denoted ℬ, is considered in a reference configuration $\kappa _{R}(\mathscr{B})$. Under the application of loads, it is assumed that the solid has a nonlinear elastic response which is governed by a finite number $N$ of material properties (defining local constitutive equations and their possible spatial variations).

These are denoted as a vector $\boldsymbol{\beta }=(\beta _{1},\beta _{2},\dots ,\beta _{q},\dots , \beta _{n},\dots ,\beta _{N})$.

At time $t$, after the application of traction ${\mathbf{h}}(t)$ on a part of the boundary of ℬ, denoted by $\Gamma _{h}$, ℬ undergoes a motion described by a mapping $\chi _{\boldsymbol{\beta }}$, depending on material properties $\boldsymbol{\beta }$, from a reference configuration $\kappa _{R}(\mathscr{B})$ to a current configuration $\kappa _{\boldsymbol{\beta }}(\mathscr{B})$, such as

$$ \textbf{x} = \chi _{\boldsymbol{\beta }}(\textbf{X},t) = \textbf{X} + \textbf{u}(\boldsymbol{\beta },X,t)~, $$

(1)

where $\textbf{X}$ and $\textbf{x}$ are position vectors relative to reference and current configurations and $\textbf{u}$ is the displacement vector field. In the following, we introduce the deformation gradient

$$ \textbf{F} = \frac{\partial \chi _{\boldsymbol{\beta }}(\textbf{X},t)}{\partial \textbf{X}}= \textbf{I}+ \frac{\partial \textbf{u}(\boldsymbol{\beta },\textbf{X},t)}{\partial \textbf{X}}, $$

(2)

the right Cauchy–Green tensor

$$ \textbf{C}=\textbf{F}^{T} \textbf{F}, $$

(3)

and the Jacobian

$$ J=\textrm{det}(\textbf{F}). $$

(4)

We also introduce the isochoric right Cauchy–Green tensor

$$ \hat{\textbf{C}}=\hat{\textbf{F}}^{T} \hat{\textbf{F}}, $$

(5)

where we denote the isochoric part of the deformation gradient as

$$ \hat{\textbf{F}}=J^{-1/3}\textbf{F}. $$

(6)

We assume that an experimental measurement of the displacement field $\textbf{u}(\textbf{X},t)$ is available across the solid at time $t$. The measured displacement field $\tilde{\textbf{u}}(t)$ may be different than the actual displacement $\textbf{u}(\textbf{X},t)$ because of measurement noise.

The inverse / identification problem consists of finding the values of the approximate parameters $\tilde{\boldsymbol{\beta }}$ that minimizes the absolute error between $\textbf{u}(\tilde{\boldsymbol{\beta }},\textbf{X},t)$ and $\tilde{\textbf{u}}(\textbf{X},t)$.

2.2 Constitutive and Equilibrium Equations

Equilibrium equations in $\kappa _{\boldsymbol{\beta }}(\mathscr{B})$ must be satisfied by the Cauchy stress $\textbf{T}$, which may be written, in absence of accelerations and body forces, as

$$ \textstyle\begin{array}{c@{\quad }c@{\quad }c@{\quad }c@{\quad }c} \boldsymbol{\nabla }_{\kappa _{\boldsymbol{\beta }}}.\textbf{T} & = & \textbf{0} & \mathrm{on} & \chi _{\boldsymbol{\beta }}(\mathscr{B},t), \\ \textbf{T}.{\mathbf{n}} & = & {\mathbf{h}}(t) & \mathrm{on} & \chi _{ \boldsymbol{\beta }}(\Gamma _{h},t). \end{array} $$

(7)

In Eq. (7), the equilibrium on the boundaries is only written where known tractions are applied as we will use these equilibrium equations to solve the identification problem and we do not want to introduce supplemental unknowns such as the reaction tractions on the boundaries where Dirichlet boundary conditions would be applied.

In this paper, we focus on compressible hyperelastic materials. It is assumed that their strain energy density function $\Phi$ can be written in an uncoupled form as

$$ \Phi = U(J,\boldsymbol{\beta }) + \hat{W}(\hat{\textbf{C}}, \boldsymbol{\beta }), $$

(8)

where $U(J,\boldsymbol{\beta })$ describes the volumetric response and $\hat{W}(\hat{\textbf{C}},\boldsymbol{\beta })$ describes the deviatoric response.

The second Piola-Kirchhoff stress $\textbf{S}$ can be written as

$$ \textbf{S}=2 \frac{\mathrm{d}\Phi }{\mathrm{d}\textbf{C}} = \frac{\mathrm{d}U}{\mathrm{d}J} J \textbf{C}^{-1} + 2 J^{-2/3} ~ \textbf{Dev} \left ( \frac{\partial \hat{W}}{\partial \hat{\textbf{C}}} \right ), $$

(9)

where

$$ \textbf{Dev} \left ( \frac{\partial \hat{W}}{\partial \hat{\textbf{C}}} \right ) = \left ( \frac{\partial \hat{W}}{\partial \hat{\textbf{C}}} \right ) - \frac{1}{3}\left [ \left ( \frac{\partial \hat{W}}{\partial \hat{\textbf{C}}} \right ): \hat{\textbf{C}} \right ]\hat{\textbf{C}}^{-1}. $$

(10)

In the following, we introduce $\hat{\textbf{S}}=2\frac{\partial \hat{W}}{\partial \hat{\textbf{C}}}$. The Cauchy stress $\textbf{T}$ can be written as

$$ \textbf{T}= J^{-1}\textbf{F}^{T}\textbf{S}\textbf{F} = \frac{\mathrm{d}U}{\mathrm{d}J}\textbf{I}+ \textbf{dev} \left ( \hat{\textbf{T}} \right ), $$

(11)

where $\hat{\textbf{T}}=J^{-1} \hat{\textbf{F}}\hat{\textbf{S}} \hat{\textbf{F}}^{T}$, and

$$ \textbf{dev} \left ( \hat{\textbf{T}} \right ) = \hat{\textbf{T}} - \frac{1}{3}\textrm{Tr}\left ( \hat{\textbf{T}} \right )\textbf{I}~. $$

(12)

Moreover, we introduce

$$ \textbf{S}_{,\beta _{q}} = \frac{\partial \textbf{S} }{\partial \beta _{q}} = \frac{\mathrm{d} U_{,\beta _{q}}}{\mathrm{d}J}J \textbf{C}^{-1} + J^{-2/3} \textbf{Dev} \left ( \hat{\textbf{S}}_{,\beta _{q}} \right ), $$

(13)

and

$$ \textbf{T}_{,\beta _{q}} = \frac{\partial \textbf{T} }{\partial \beta _{q}} = \frac{\mathrm{d} U_{,\beta _{q}}}{\mathrm{d}J}\textbf{I}+ J^{-1} \textbf{dev} \left ( \hat{\textbf{F}}~\hat{\textbf{S}}_{,\beta _{q}}~ \hat{\textbf{F}}^{T} \right ), $$

(14)

where

$$ \hat{\textbf{S}}_{,\beta _{q}} = \frac{\partial \hat{\textbf{S}}}{\partial \beta _{q}} = 2 \frac{\partial \frac{\partial \hat{W}}{\partial \hat{\textbf{C}}}}{\partial \beta _{q}} = 2 \frac{\partial \frac{\partial \hat{W}}{\partial \beta _{q}}}{\partial \hat{\textbf{C}}} = 2 \frac{\partial \hat{W}_{,\beta _{q}}}{\partial \hat{\textbf{C}}}, $$

(15)

$$ \hat{W}_{,\beta _{q}} = \frac{\partial \hat{W}}{\partial \beta _{q}} ~~~~ \textrm{and}~~~~ U_{,\beta _{q}} = \frac{\partial U}{\partial \beta _{q}}. $$

(16)

2.3 Definition of an Intermediate Configuration

2.3.1 Decomposition of the Deformation Gradient

Let us start with a set of parameters $\boldsymbol{\beta }^{o}$ which corresponds to an initial estimation of the unknown constitutive parameters based on existing literature on the same tissue for instance.

Unstressed “intermediate configuration” were traditionally used in the literature for the multiplicative split of the deformation gradient in finite strain plasticity [38]. We define here a stressed intermediate configuration $\kappa _{\boldsymbol{\beta }^{o}}(\mathscr{B})$ for which the solid with the vector of material properties $\boldsymbol{\beta }^{o}$ is at equilibrium (Fig. 1). We emphasize that this is a different use of the “intermediate configuration” terminology than in finite strain plasticity [38].

Fig. 1

The relationship of the reference, intermediate and current configurations

Then, let the position in the intermediate (stressed) configuration be denoted by $\textbf{x}^{o} = \chi _{\boldsymbol{\beta }^{o}}(\textbf{X},t)$, corresponding to displacement $\textbf{u}^{o}(\textbf{X},t) = \textbf{x}^{o} - \textbf{X}$.

We assume that a small displacement $\boldsymbol{\delta }\textbf{x}^{o}$ superimposed upon the large deformation $\textbf{u}^{o}$, yields the current position $\textbf{x}$ at time $t$ for which the solid with a vector of material properties $\boldsymbol{\beta }=\boldsymbol{\beta }^{o} + \boldsymbol{\delta }\boldsymbol{\beta }^{o}$ is at equilibrium. We assume that $\boldsymbol{\delta }\boldsymbol{\beta }^{o}$ is a small variation of $\boldsymbol{\beta }^{o}$ such as

$$ \textstyle\begin{array}{c@{\quad }c@{\quad }c} ||\boldsymbol{\delta }\textbf{x}^{o}|| & \ll & ||\textbf{u}^{o} ||, \\ ||\boldsymbol{\delta }\boldsymbol{\beta }^{o}|| & \ll & || \boldsymbol{\beta }^{o}||. \end{array} $$

(17)

The small displacement $\boldsymbol{\delta }\textbf{x}^{o}$ corresponds to the deformation between the intermediate configuration (at which the solid with $\boldsymbol{\beta }^{o}$ material properties is at equilibrium) and the current configuration (at which the solid with $\boldsymbol{\beta }=\boldsymbol{\beta }^{o} + \boldsymbol{\delta }\boldsymbol{\beta }^{o}$ material properties is at equilibrium).

The current position can thus be written as

$$ \textbf{x} = \textbf{x}^{o} + \boldsymbol{\delta }\textbf{x}^{o}. $$

(18)

The deformation gradient associated with mappings from the reference to the intermediate is thus given by

$$ \textbf{F}^{o} = \frac{\partial \chi _{\boldsymbol{\beta }^{o}}(\textbf{X},t)}{\partial \textbf{X}}~. $$

(19)

The deformation gradient representing a mapping from the intermediate configuration to the current configuration (corresponding to the variations of the configuration for a variation of material properties $\boldsymbol{\delta }\boldsymbol{\beta }^{o}$) may be expressed as $\textbf{I} + \boldsymbol{\delta }\textbf{F}^{o}$ where

$$ \boldsymbol{\delta }\textbf{F}^{o} = \frac{\partial \boldsymbol{\delta }\textbf{x}^{o}}{\partial \textbf{x}^{o}}~. $$

(20)

Then, gradients of the successive motions are obtained with the chain rule, such as

$$ \textbf{F}= \frac{\partial \textbf{x}}{\partial \textbf{x}^{o}} \frac{\partial \textbf{x}^{o}}{\partial \textbf{X}} = \frac{\partial \textbf{x}}{\partial \textbf{x}^{o}} \textbf{F}^{o} = \left [ \frac{\partial (\textbf{x}^{o}+\boldsymbol{\delta }\textbf{x}^{o})}{\partial \textbf{x}^{o}} \right ] \textbf{F}^{o} = \frac{\partial \textbf{x}^{o}}{\partial \textbf{x}^{o}}\textbf{F}^{o} + \frac{\partial \boldsymbol{\delta }\textbf{x}^{o}}{\partial \textbf{x}^{o}} \textbf{F}^{o} = \textbf{F}^{o} + \boldsymbol{\delta }\textbf{F}^{o} ~ \textbf{F}^{o} ~. $$

(21)

Then, the identification problem can be formulated as, find the values of $\boldsymbol{\delta }\boldsymbol{\beta }^{o}$ that minimizes the error between $\boldsymbol{\delta }\textbf{x}^{o}$ and $\tilde{\textbf{u}}-\textbf{u}^{o}$.

2.3.2 Cauchy Stress Tensor

Using Eq. (9) and Eq. (11) specifically with $\textbf{F}^{o}$ deformation gradient, the second Piola-Kirchhoff stress $\textbf{S}^{o}$ and Cauchy stress $\textbf{T}^{o}$ in the intermediate configuration may be written respectively

$$ \textstyle\begin{array}{c@{\quad }c@{\quad }c} \textbf{S}^{o} & = & \displaystyle \frac{\mathrm{d}U^{o}}{\mathrm{d}J} J^{o} \left (\textbf{C}^{o}\right )^{-1} + 2 \left (J^{o}\right )^{-2/3} \textbf{Dev} \left ( \frac{\partial \hat{W}^{o}}{\partial \hat{\textbf{C}}} \right )~, \\ \textbf{T}^{o} & = & \displaystyle \frac{\mathrm{d}U^{o}}{\mathrm{d}J} \textbf{I}+ 2\left (J^{o}\right )^{-1} \textbf{dev} \left ( \hat{\textbf{F}^{o}} \frac{\partial \hat{W}^{o}}{\partial \hat{\textbf{C}}} \hat{\textbf{F}^{o}}^{T} \right )~, \end{array} $$

(22)

where $U^{o}=U(J^{o},\boldsymbol{\beta }^{o})$ and $\hat{W}^{o}=\hat{W}(\hat{\textbf{F}^{o}},\boldsymbol{\beta }^{o})$.

Note that the Cauchy stress $\textbf{T}^{o}$ must satisfy the equilibrium equations on the $\kappa _{\boldsymbol{\beta }^{o}}$ configuration (further denoted $\kappa _{o}$) which may be rewritten

$$ \textstyle\begin{array}{c@{\quad }c@{\quad }c@{\quad }c@{\quad }c} \boldsymbol{\nabla }.\textbf{T}^{o} & = & \textbf{0} & \mathrm{on} & \chi _{\boldsymbol{\beta }^{o}}(\mathscr{B},t)~, \\ \textbf{T}^{o}.{\mathbf{n}} & = & {\mathbf{h}}(t) & \mathrm{on} & \chi _{ \boldsymbol{\beta }^{o}}(\Gamma _{h},t)~. \end{array} $$

(23)

The second Piola-Kirchhoff stress $\textbf{S}$ can be related to $\textbf{S}^{o}$ using a Taylor expansion of first order. For that we substitute into Eq. (9), $\boldsymbol{\beta }=\boldsymbol{\beta }^{o}+\boldsymbol{\delta }\boldsymbol{\beta }^{o}$ and $\textbf{C}=\textbf{C}^{o}+\boldsymbol{\delta }\textbf{C}^{o}$, which yields

$$ \mathbf{S} \simeq \textbf{S}^{o} + \boldsymbol{\delta }\textbf{S}^{o} \simeq \textbf{S}^{o} + \frac{\partial \textbf{S}}{\partial \textbf{C}} : \boldsymbol{\delta }\textbf{C}^{o} + \sum _{q=1}^{N} \delta \beta _{q}^{o}~\textbf{S}^{o}_{, \beta _{q}} ~, $$

(24)

where

$$ \textbf{S}^{o}_{,\beta _{q}} = \frac{\mathrm{d}U_{,\beta _{q}}^{o}}{\mathrm{d}J}J^{o} \left ( \textbf{C}^{o}\right )^{-1} + 2 \left (J^{o}\right )^{-2/3} \textbf{Dev} \left ( \frac{\partial \hat{W}_{,\beta _{q}}^{o}}{\partial \hat{\textbf{C}}} \right ) ~. $$

(25)

Neglecting the second order terms in the deformations from the intermediate to the current configuration, it may be written

$$ \textstyle\begin{array}{l@{\quad }c@{\quad }l} \boldsymbol{\delta }\textbf{C}^{o} & = & \left (\left ( \boldsymbol{\delta }\textbf{F}^{o}+\textbf{I} \right )\textbf{F}^{o} \right )^{T} \left (\left (\boldsymbol{\delta }\textbf{F}^{o}+\textbf{I} \right )\textbf{F}^{o} \right ) -\textbf{F}^{oT} \textbf{F}^{o} \\ & = & {\textbf{F}^{o}}^{T} \left ( \left (\boldsymbol{\delta }\textbf{F}^{o}+ \textbf{I} \right )^{T} \left (\boldsymbol{\delta }\textbf{F}^{o}+ \textbf{I} \right ) - \textbf{I} \right ) \textbf{F}^{o} \\ & = & {\textbf{F}^{o}}^{T} \left (\boldsymbol{\delta }\textbf{F}^{o}+ \left (\boldsymbol{\delta }\textbf{F}^{o} \right )^{T} + \left ( \boldsymbol{\delta }\textbf{F}^{o} \right )^{T} \boldsymbol{\delta }\textbf{F}^{o} \right ) \textbf{F}^{o} \\ & \simeq & 2 {\textbf{F}^{o}}^{T} \boldsymbol{\delta }\textbf{E}^{o} ~ \textbf{F}^{o}~, \end{array} $$

(26)

where $\boldsymbol{\delta }\textbf{E}^{o}=\frac{1}{2}\left (\boldsymbol{\delta }\textbf{F}^{o}+\left (\boldsymbol{\delta }\textbf{F}^{o} \right )^{T} + \left (\boldsymbol{\delta }\textbf{F}^{o} \right )^{T} \boldsymbol{\delta }\textbf{F}^{o} \right )$ signifies the infinitesimal strain induced by a slight variation of material properties when $\boldsymbol{\delta }\textbf{F}^{o}$ is small.

Finally,

$$ \boldsymbol{\delta }\textbf{S}^{o} = \pmb{\mathbb{K}}^{o} : \left ( { \textbf{F}^{o}}^{T} \boldsymbol{\delta }\textbf{E}^{o} ~\textbf{F}^{o} \right ) + \sum _{q=1}^{N} \delta \beta _{q}^{o}~\textbf{S}_{,\beta _{q}}^{o} ~, $$

(27)

where

$$ \pmb{\mathbb{K}}^{o} = 2 \frac{\partial \textbf{S}^{o}}{\partial \textbf{C}}. $$

(28)

Next, we can derive the Cauchy stress $\textbf{T}$ for the current configuration from $\textbf{T}^{o}$ as

$$ \textstyle\begin{array}{c@{\quad }c@{\quad }l} \mathbf{T} & = & \textrm{det}(\textbf{F}^{o}+\boldsymbol{\delta }\textbf{F}^{o})^{-1}\left (J^{o}\right )^{-1}\left (\textbf{I} + \boldsymbol{\delta }\textbf{F}^{o} \right ) \textbf{F}^{o} \left ( \textbf{S}^{o}+ \boldsymbol{\delta }\textbf{S}^{o} \right ) \left ( \textbf{F}^{o}\right )^{T} \left (\textbf{I} + \boldsymbol{\delta }\textbf{F}^{o} \right )^{T} \\ & \simeq & \mathbf{T}^{o} + \boldsymbol{\delta }\textbf{T}^{o}~, \end{array} $$

(29)

where

$$ \textstyle\begin{array}{l} \boldsymbol{\delta }\textbf{T}^{o} = - \textrm{Tr}(\boldsymbol{\delta }\textbf{E}^{o}) \mathbf{T}^{o} + \boldsymbol{\delta }\textbf{F}^{o} \mathbf{T}^{o} + \mathbf{T}^{o} \left (\boldsymbol{\delta }\textbf{F}^{o} \right )^{T} \\ + \displaystyle \sum _{n=1}^{N} \delta \beta _{n}~ \textbf{T}_{,\beta _{n}}^{o} + \left (J^{o}\right )^{-1} \textbf{F}^{o} \left [ \pmb{\mathbb{K}}^{o} : \left ( {\textbf{F}^{o}}^{T} \boldsymbol{\delta }\textbf{F}^{o}~ \textbf{F}^{o} \right ) \right ] \left (\textbf{F}^{o}\right )^{T}~, \end{array} $$

(30)

and

$$ \textbf{T}_{,\beta _{q}}^{o} = \frac{\mathrm{d} U_{,\beta _{q}}^{o}}{\mathrm{d}J} \textbf{I}+ 2\left (J^{o} \right )^{-1} \textbf{dev} \left ( \hat{\textbf{F}^{o}} \frac{\partial \hat{W}_{,\beta _{q}}^{o}}{\partial \hat{\textbf{C}}} \hat{\textbf{F}^{o}}^{T} \right )~. $$

(31)

Equation 30 may be rewritten by introducing a 4$^{\text{th}}$ order tensor $\pmb{\mathbb{L}}^{o}$ such as

$$ \boldsymbol{\delta }\textbf{T}^{o} = \pmb{\mathbb{L}}^{o}: \boldsymbol{\delta }\textbf{F}^{o} + \displaystyle \sum _{q=1}^{N} \delta \beta _{q}~ \textbf{T}_{,\beta _{q}}^{o}~. $$

(32)

The components of the $\pmb{\mathbb{L}}^{o}$ tensor may be written such as

$$ \mathbb{L}^{o}_{ijkl}=\frac{1}{2}(\mathscr{C}^{o}_{ijkl}+\mathscr{D}^{o}_{ijkl}+ \mathscr{C}^{o}_{ijlk}-\mathscr{D}^{o}_{ijlk})~, $$

(33)

where, the summation convention being adopted for the repeated indices,

$$ \textstyle\begin{array}{c@{\quad }c@{\quad }l} \mathscr{C}^{o}_{ijkl} & = & - \delta _{kl} T^{o}_{ij} + \delta _{ik}T^{o}_{lj} + T^{o}_{ik}\delta _{jl} + \left (J^{o}\right )^{-1} F^{o}_{iA}F^{o}_{jB}F^{o}_{kC}F^{o}_{lD} {\mathbb{K}}^{o}_{ABCD}~, \\ \mathscr{D}^{o}_{ijkl} & = & \delta _{ik}T^{o}_{lj} + T^{o}_{ik} \delta _{jl}~. \end{array} $$

(34)

Equilibrium must be satisfied in the current configuration. Given that Eq. (23) is satisfied by $\textbf{T}^{o}$ on the intermediate configuration, and assuming that the intermediate configuration is infinitesimally close to the current configuration (owing to $\boldsymbol{\delta }\boldsymbol{\beta }^{o}/\boldsymbol{\beta }^{o} \ll 1$), we can assume that the following equations are satisfied by $\boldsymbol{\delta }\textbf{T}^{o}$,

$$ \textstyle\begin{array}{c@{\quad }c@{\quad }c@{\quad }c@{\quad }c} \boldsymbol{\nabla }_{\kappa _{o}}.\boldsymbol{\delta }\textbf{T}^{o} & = & \textbf{0} & \mathrm{on} & \chi _{\boldsymbol{\beta }^{o}}(\mathscr{B},t) ~, \\ \boldsymbol{\delta }\textbf{T}^{o}.{\mathbf{n}} & = & \textbf{0} & \mathrm{on} & \chi _{\boldsymbol{\beta }^{o}}(\Gamma _{h},t)~. \end{array} $$

(35)

The equilibrium equations may be rewritten such as

$$ \textstyle\begin{array}{c@{\quad }c@{\quad }c@{\quad }c@{\quad }c} \displaystyle \sum _{q=1}^{N} \delta \beta _{q}~ \boldsymbol{\nabla }. \textbf{T}_{,\beta _{q}}^{o} & = & -\boldsymbol{\nabla }_{\kappa _{o}}. \left (\pmb{\mathbb{L}}^{o}: \boldsymbol{\delta }\textbf{F}^{o} \right ) & ~~~~\mathrm{on} & \chi _{\boldsymbol{\beta }^{o}}(\mathscr{B},t) ~, \\ \displaystyle \sum _{q=1}^{N} \delta \beta _{q}~\textbf{T}_{,\beta _{n}}^{o}.{ \mathbf{n}} & = & -\left (\pmb{\mathbb{L}}^{o}: \boldsymbol{\delta }\textbf{F}^{o} \right ).{\mathbf{n}} & ~~~~\mathrm{on} & \chi _{ \boldsymbol{\beta }^{o}}(\Gamma _{h},t)~. \end{array} $$

(36)

In the following, for the sake of simplification, the following abuse of notation will be used for the gradient and divergence in the $\kappa _{o}$ configuration: $\boldsymbol{\nabla }\equiv \boldsymbol{\nabla }_{\kappa _{o}}$.

2.4 The Virtual Fields Method

$\boldsymbol{\delta }\textbf{T}^{o}$ should satisfy the equilibrium equations (Eq. (36)) written just above. These equations can be written in their weak form such as

$$ \int _{\chi _{\boldsymbol{\beta }^{o}}(\mathscr{B},t)} \boldsymbol{\delta }\textbf{T}^{o} : \boldsymbol{\nabla }\boldsymbol{\delta }\mathfrak{u}^{o(n)} dV^{o} = 0~, $$

(37)

where $\boldsymbol{\delta }\mathfrak{u}^{o(n)}$ is a virtual displacement field which equals zero on $\chi _{\boldsymbol{\beta }^{o}}(\partial \mathscr{B} \backslash \Gamma _{h})$ (boundary where traction are not applied). In $\boldsymbol{\delta }\mathfrak{u}^{o(n)}$, the index $(n)$ indicates that at least $N$ virtual fields are necessary to establish a system of $N$ equations of the $N$ unknown material properties $\beta _{n}$.

As full-field measurements are available, $\boldsymbol{\delta }\textbf{F}^{o}$ and $\boldsymbol{\delta }\textbf{E}^{o}$ from Eq. (30) can be replaced by their measures denoted respectively $\tilde{\textbf{H}}$ and $\tilde{\boldsymbol{\epsilon }}$, such as

$$ \tilde{\textbf{H}}=\boldsymbol{\nabla }(\tilde{\textbf{u}}-\textbf{u}^{o})~,~ \tilde{\boldsymbol{\epsilon }}=\frac{1}{2}(\tilde{\textbf{H}}+ \tilde{\textbf{H}}^{T})~, $$

(38)

yielding

$$ \textstyle\begin{array}{l} \displaystyle \sum _{q=1}^{N} \delta \beta _{q}^{o}~\displaystyle \int _{\chi _{\boldsymbol{\beta }^{o}}(\mathscr{B})} \textbf{T}_{,\beta _{q}}^{o} : \boldsymbol{\nabla }\boldsymbol{\delta }\mathfrak{u}^{o(n)} dV^{o} \\ = - \displaystyle \int _{\chi _{\boldsymbol{\beta }^{o}}(\mathscr{B})} \left ( - \textrm{Tr}(\tilde{\boldsymbol{\epsilon }}) \mathbf{T}^{o} + \tilde{\textbf{H}}\mathbf{T}^{o} + \mathbf{T}^{o} \tilde{\textbf{H}}^{T} \right ): \boldsymbol{\nabla }\boldsymbol{\delta }\mathfrak{u}^{o(n)} dV^{o} \\ - \displaystyle \int _{\chi _{\boldsymbol{\beta }^{o}}(\mathscr{B})} \left ( \left (J^{o}\right )^{-1} \textbf{F}^{o} \left [ \pmb{\mathbb{K}}^{o} : \left ( {\textbf{F}^{o}}^{T} \tilde{\boldsymbol{\epsilon }} \textbf{F}^{o} \right ) \right ] \left ( \textbf{F}^{o}\right )^{T} \right ) : \boldsymbol{\nabla }\boldsymbol{\delta }\mathfrak{u}^{o(n)} dV^{o} ~. \end{array} $$

(39)

Equation 39 may be rewritten by introducing the $\pmb{\mathbb{L}}^{o}$ tensor such as

$$ \textstyle\begin{array}{l} \displaystyle \sum _{q=1}^{N} \delta \beta _{q}^{o}~\displaystyle \int _{\chi _{\boldsymbol{\beta }^{o}}(\mathscr{B})} \textbf{T}_{,\beta _{q}}^{o} : \boldsymbol{\nabla }\boldsymbol{\delta }\mathfrak{u}^{o(n)} dV^{o} \\ = - \displaystyle \int _{\chi _{\boldsymbol{\beta }^{o}}(\mathscr{B})} \left ( \pmb{\mathbb{L}}^{o}: \tilde{\textbf{H}} \right ): \boldsymbol{\nabla }\boldsymbol{\delta }\mathfrak{u}^{o(n)} dV^{o} \\ = - \displaystyle \int _{\chi _{\boldsymbol{\beta }^{o}}(\mathscr{B})} \tilde{\textbf{H}} : \left ( \boldsymbol{\pmb{\mathbb{L}}^{o}}^{T}: \nabla \boldsymbol{\delta }\mathfrak{u}^{o(n)} \right ) dV^{o} ~. \end{array} $$

(40)

The previous equation can be written $N$ times with $N$ virtual fields $\boldsymbol{\delta }\mathfrak{u}^{o(n)}$. The obtained system of equations may be written such as

(41)

In the next section we provide a methodology to choose the set of $N$ virtual fields.

2.5 Derivation of the Virtual Fields for Parameter Identification

In the previous subsection, we showed that the unknown constitutive parameters can be identified by solving a linear system of equations. To establish this system of equations, one needs to define $N$ virtual fields denoted $\boldsymbol{\delta }\mathfrak{u}^{o(n)}$, such that $[\mathbf{A}^{o}]$ is invertible. The invertibility is ensured by relating each virtual field $\boldsymbol{\delta }\mathfrak{u}^{o(n)}$ to the sensitivity of the Cauchy stress to each unknown parameter, denoted as $\textbf{T}_{,\beta _{n}}^{o}$ in Eq. (31) and introduced in Eq. (30), since it is assumed that the different $\textbf{T}_{,\beta _{n}}^{o}$ constitute a set of $N$ linearly independent tensorial functions.

Therefore, a possible choice of $N$ linearly independent virtual fields could simply be: $\boldsymbol{\delta }\mathfrak{u}^{o(n)}=\boldsymbol{\nabla }.\textbf{T}_{, \beta _{n}}^{o}$. However, the virtual field must equal zero on $\chi _{\boldsymbol{\beta }^{o}}(\partial \mathscr{B} \backslash \Gamma _{h})$ and remain continuous. To meet these requirements and benefit from the linear independence between the $\textbf{T}_{,\beta _{n}}^{o}$ fields, we defined $\boldsymbol{\delta }\mathfrak{u}^{o(n)}$ as the vectorial fields satisfying

$$ \textstyle\begin{array}{c@{\quad }c@{\quad }c@{\quad }c@{\quad }c} \boldsymbol{\nabla }.\left ({\pmb{\mathbb{L}}^{o}}^{T}: \boldsymbol{\nabla }\boldsymbol{\delta }\mathfrak{u}^{o(n)} \right ) & = & \boldsymbol{\nabla }.\left (\textbf{T}_{,\beta _{n}}^{o}\right ) & \mathrm{on} & \chi _{\boldsymbol{\beta }^{o}}(\mathscr{B},t), \\ \left ({\pmb{\mathbb{L}}^{o}}^{T}:\boldsymbol{\nabla }\boldsymbol{\delta }\mathfrak{u}^{o(n)} \right ).\textbf{n} & = & \left ( \textbf{T}_{,\beta _{n}}^{o}\right ) .\textbf{n}~~ & \mathrm{on} & \chi _{\boldsymbol{\beta }^{o}}(\Gamma _{h}), \\ \boldsymbol{\delta }\mathfrak{u}^{o(n)} & = & \textbf{0} & \mathrm{on} & \chi _{\boldsymbol{\beta }^{o}}(\partial \mathscr{B} \backslash \Gamma _{h}). \end{array} $$

(42)

The virtual fields $\boldsymbol{\delta }\mathfrak{u}^{o(n)}$ are eventually obtained by solving the linear elastic problems defined in Eq. (42) using the finite-element method (Sect. 2.6). Let us introduce the $\mathfrak{L}$ linear operator such as $\boldsymbol{\delta }\mathfrak{u}^{o(n)}=\mathfrak{L}(\textbf{T}_{, \beta _{n}}^{o})$ is the solution of Eq. (42).

Moreover using the integration by parts, we can write

$$ \textstyle\begin{array}{l} \displaystyle \int _{\chi _{\boldsymbol{\beta }^{o}}(\mathscr{B})} \tilde{\textbf{H}} : \left ( {\pmb{\mathbb{L}}^{o}}^{T} : \boldsymbol{\nabla }\boldsymbol{\delta }\mathfrak{u}^{n} \right ) dV^{o} = \displaystyle \oint _{\chi _{\boldsymbol{\beta }^{o}}(\Gamma _{h})} ( \tilde{\textbf{u}}-\textbf{u}^{o}) . \left ( {\pmb{\mathbb{L}}^{o}}^{T} : \boldsymbol{\nabla }\boldsymbol{\delta }\mathfrak{u}^{n} \right ). \textbf{n}~dS^{o} \\ - \displaystyle \int _{\chi _{\boldsymbol{\beta }^{o}}(\mathscr{B})} ( \tilde{\textbf{u}}-\textbf{u}^{o}) . \boldsymbol{\nabla }. \left ( { \pmb{\mathbb{L}}^{o}}^{T} : \boldsymbol{\nabla }\boldsymbol{\delta }\mathfrak{u}^{n} \right ) dV^{o}~, \\ \\ \displaystyle \int _{\chi _{\boldsymbol{\beta }^{o}}(\mathscr{B})} \tilde{\textbf{H}} : \left ( {\pmb{\mathbb{L}}^{o}}^{T} : \boldsymbol{\nabla }\boldsymbol{\delta }\mathfrak{u}^{n} \right ) dV^{o} = \displaystyle \oint _{\chi _{\boldsymbol{\beta }^{o}}(\Gamma _{h})} ( \tilde{\textbf{u}}-\textbf{u}^{o}) . \left ( \textbf{T}_{,\beta _{n}}^{o} \right ).\textbf{n}~dS^{o} \\ - \displaystyle \int _{\chi _{\boldsymbol{\beta }^{o}}(\mathscr{B})} ( \tilde{\textbf{u}}-\textbf{u}^{o}) . \boldsymbol{\nabla }.\left ( \textbf{T}_{,\beta _{n}}^{o}\right ) dV^{o}~, \\ \\ \displaystyle \int _{\chi _{\boldsymbol{\beta }^{o}}(\mathscr{B})} \tilde{\textbf{H}} : \left ( {\pmb{\mathbb{L}}^{o}}^{T} : \boldsymbol{\nabla }\boldsymbol{\delta }\mathfrak{u}^{n} \right ) dV^{o} = \displaystyle \int _{\chi _{\boldsymbol{\beta }^{o}}(\mathscr{B})} \tilde{\textbf{H}} : \textbf{T}_{,\beta _{n}}^{o} dV^{o}. \end{array} $$

(43)

Then we can build the system of equations of Eq. (41) by replacing each $\boldsymbol{\delta }\mathfrak{u}^{o(n)}$ by their expression coming from Eq. (42), and use also the simplifications derived in Eq. (43), yielding

$$\begin{array}{c}\left[\begin{array}{ccc}{\int }_{{\chi }_{{\mathit{\beta }}^o}(\mathcal{B})}{\mathbf{\text{T}}}_{,{\beta }_1}^o:\mathbf{\nabla }\mathfrak{L}({\mathbf{\text{T}}}_{,{\beta }_1}^o)d{V}^o& \dots & {\int }_{{\chi }_{{\mathit{\beta }}^o}(\mathcal{B})}{\mathbf{\text{T}}}_{,{\beta }_N}^o:\mathbf{\nabla }\mathfrak{L}({\mathbf{\text{T}}}_{,{\beta }_1}^o)d{V}^o\\ & & \\ ⋮& \ddots & ⋮\\ & & \\ {\int }_{{\chi }_{{\mathit{\beta }}^o}(\mathcal{B})}{\mathbf{\text{T}}}_{,{\beta }_1}^o:\mathbf{\nabla }\mathfrak{L}({\mathbf{\text{T}}}_{,{\beta }_N}^o)d{V}^o& \dots & {\int }_{{\chi }_{{\mathit{\beta }}^o}(\mathcal{B})}{\mathbf{\text{T}}}_{,{\beta }_N}^o:\mathbf{\nabla }\mathfrak{L}({\mathbf{\text{T}}}_{,{\beta }_N}^o)d{V}^o\end{array}\right]\hfill \\ \hfill \times \left(\begin{array}{c}\delta {\beta }_1^o\\ \\ \delta {\beta }_2^o\\ \\ ⋮\\ \\ \delta {\beta }_N^o\end{array}\right)=-\left(\begin{array}{c}{\int }_{{\chi }_{{\mathit{\beta }}^o}(\mathcal{B})}\widetilde{\mathbf{\text{H}}}:{\mathbf{\text{T}}}_{,{\beta }_1}^{o}d{V}^o\\ \\ ⋮\\ \\ {\int }_{{\chi }_{{\mathit{\beta }}^o}(\mathcal{B})}\widetilde{\mathbf{\text{H}}}:{\mathbf{\text{T}}}_{,{\beta }_N}^{o}d{V}^o\end{array}\right).\end{array}$$

(44)

This system of equations can be rewritten

$$ [\mathbf{A}^{o}] \{ \boldsymbol{\delta }\boldsymbol{\beta }^{o} \} = \{ \mathbf{b}^{o} \} ~. $$

(45)

Then, the unknown constitutive parameters can be obtained as

$$ \boldsymbol{\beta }= \boldsymbol{\beta }^{o} + [\mathbf{A}^{o}]^{-1}\{ \mathbf{b}^{o} \} ~. $$

(46)

2.6 Finite-Element Implementation

The numerical implementation of the proposed method requires the use of finite-element analyses at two different levels:

• 1. find the intermediate configuration and the deformation gradient $\textbf{F}^{o}$ by solving the partial differential equations (PDEs) of the forward problem (Eq. (23)) with a set of parameters $\boldsymbol{\beta }^{o}$,

• 2. compute integrals in Eq. (41) or Eq. (44) to identify the unknown constitutive parameters.

In this section, we propose a possible numerical implementation. As standard in finite-elements [39–42]), the computational domain ℬ is discretised into a finite number of non-overlapping elements $e \in \mathbb{E}$ such that

$$ \textstyle\begin{array}{c@{\quad }c@{\quad }c@{\quad }c@{\quad }c} \mathscr{B} & \approx & \mathscr{B}^{h} & = & \bigcup \limits _{e \in \mathbb{E}}\mathscr{B}^{e}~, \\ \chi _{\boldsymbol{\beta }^{o}}(\mathscr{B}) & \approx & \chi _{ \boldsymbol{\beta }^{o}}(\mathscr{B}^{h}) & = & \bigcup \limits _{e \in \mathbb{E}} \chi _{\boldsymbol{\beta }^{o}}(\mathscr{B}^{e}) ~.\end{array} $$

(47)

The fields are discretised using the following standard vectorial shape functions $\pmb{\mathscr{N}}(\mathbf{X})$ as

$$ \textstyle\begin{array}{c@{\quad }c@{\quad }c} \mathbf{u}^{o}(\mathbf{X}) & = & \displaystyle \sum _{a=1}^{N_{w}} {u}^{o}_{a} \pmb{\mathscr{N}}^{a}(\mathbf{X})~, \\ \tilde{\mathbf{u}}(\mathbf{X}) & = & \displaystyle \sum _{a=1}^{N_{w}} \tilde{{u}}_{a} \pmb{\mathscr{N}}^{a}(\mathbf{X})~, \\ \boldsymbol{\delta }\mathfrak{u}^{o(n)}(\mathbf{X}) & = & \displaystyle \sum _{a=1}^{N_{w}} \delta {u}^{o(n)}_{a} \pmb{\mathscr{N}}^{a}(\mathbf{X})~, \end{array} $$

(48)

where ${u}^{o}_{a}$ are nodal variables for the $\textbf{u}^{o}$ field, $\tilde{{u}}_{a}$ are nodal variables for the $\textbf{u}^{o}$ field and $\delta {u}^{o(n)}_{a}$ are nodal variables for the $\boldsymbol{\delta }\mathfrak{u}^{o(n)}$ field.

A standard finite-element discretisation enables the following discrete form for the components of Eq. (44),

$$ \textstyle\begin{array}{c@{\quad }c@{\quad }l} \{\textbf{b}^{o}\}_{n} & = & \displaystyle \int _{\chi _{ \boldsymbol{\beta }^{o}}(\mathscr{B})} \tilde{\textbf{H}} : \textbf{T}_{, \beta _{n}}^{o} dV^{o} \\ & & \\ & = & \displaystyle \sum _{a=1}^{N_{w}} \left (\tilde{{u}}_{a}-{u}_{a}^{o} \right ) {\mathscr{F}}^{o(n)}_{a} \\ & & \\ & = & \left ( \pmb{\mathscr{F}}^{o(n)} \right )^{T} \left ( \tilde{\textbf{u}}-\textbf{u}^{o}\right )~, \end{array} $$

(49)

where $\pmb{\mathscr{F}}^{o(n)}$ is obtained such as

$$ {\mathscr{F}}^{o(n)}_{a} = \displaystyle \int _{\chi _{ \boldsymbol{\beta }^{o}}(\mathscr{B})} \boldsymbol{\nabla }\pmb{\mathscr{N}}(\textbf{x}^{o}) : \textbf{T}_{,\beta _{n}}^{o} dV^{o}~, $$

(50)

and

$$ \textstyle\begin{array}{c@{\quad }c@{\quad }l} [\textbf{A}^{o}]_{qn} & = & \displaystyle \int _{\chi _{ \boldsymbol{\beta }^{o}}(\mathscr{B})} \textbf{T}_{,\beta _{q}}^{o} : \boldsymbol{\nabla }\mathfrak{L}(\textbf{T}_{,\beta _{n}}^{o}) dV^{o} \\ & & \\ & = & \displaystyle \sum _{a=1}^{N_{w}}\sum _{b=1}^{N_{w}} { \mathscr{F}}^{o(n)}_{a} [{\mathscr{K}}^{o}]^{-1}_{ab} {\mathscr{F}}^{o(q)}_{b} \\ & & \\ & = & \left ( \pmb{\mathscr{F}}^{o(n)} \right )^{T} [{ \pmb{\mathscr{K}}^{o}}]^{-1} \pmb{\mathscr{F}}^{o(q)}~, \end{array} $$

(51)

where $\pmb{\mathscr{K}}^{o}$, which is the stiffness matrix needed to solve the elastic problem of Eq. (42) to derive the virtual fields, and which satisfies $\pmb{\mathscr{K}}^{o} \boldsymbol{\delta }\mathfrak{u}^{o(n)} = \pmb{\mathscr{F}}^{o(n)}$, is obtained such as

$$ {\mathscr{K}}^{o}_{ab} = \int _{\chi _{\boldsymbol{\beta }^{o}}( \mathscr{B})} \left (\boldsymbol{\nabla }\pmb{\mathscr{N}}(\textbf{x}^{o}) \right ):\left ({\pmb{\mathbb{L}}^{o}}^{T}(\textbf{x}^{o})\right ): \left (\boldsymbol{\nabla }\pmb{\mathscr{N}}(\textbf{x}^{o})\right )~dV^{o} ~. $$

(52)

Finally, the expression of $[\textbf{A}^{o}]$ and $\{\textbf{b}^{o}\}$ with this finite-element discretisation are

$$ \{\textbf{b}^{o}\}_{n} =- \left (\pmb{\mathscr{F}}^{o(n)}\right )^{T} \left (\tilde{\textbf{u}}-\textbf{u}^{o}\right )~, $$

(53)

and

$$ [\textbf{A}^{o}]_{qn} = \left (\pmb{\mathscr{F}}^{o(n)}\right )^{T} \left [\pmb{\mathscr{K}}^{o}\right ]^{-1} \pmb{\mathscr{F}}^{o(q)}~. $$

(54)

2.7 Convergence of Parameter Identification

Based on the details discussed in the previous subsection, it is possible to derive $\boldsymbol{\delta }\boldsymbol{\beta }^{o}$ from the choice of an initial guess of $\boldsymbol{\beta }^{o}$ by solving Eq. (45). The question that remains to be solved is how to choose $\boldsymbol{\beta }^{o}$. In practice, an initial estimation of the unknown constitutive parameters based on existing literature on the same tissue can be used to initiate the resolution. However, this does not guarantee the criterion in Eq. (17) is satisfied. Then, a non infinitesimal deviation between the intermediate configuration and the reference configuration would cause Eq. (36) to not be satisfied on $\chi _{\boldsymbol{\beta }^{o}}(\mathscr{B},t)$. Therefore, the concept of “intermediate configuration” has to be iterative. From the choice of a first set of parameters $\boldsymbol{\beta }^{o}$, an intermediate configuration can be found by solving the forward problem in Eq. (23). Evaluating $[\textbf{A}^{o}]$ and $\{\textbf{b}^{o}\}$ (from equations 53 and 54) using the obtained $\pmb{\mathbb{L}^{o}}$ and $\textbf{T}_{,\beta _{n}}^{o}$ expressions (from equations 33 and 31) provides an update for $\boldsymbol{\beta }^{o}$. The process is repeated until the deviation between $\mathbf{u}^{o}$ and $\tilde{\mathbf{u}}$ becomes small enough (Fig. 2).

Fig. 2

Flowchart of the proposed inverse method

The convergence criterion of the inverse algorithm is that the relative difference between the current estimation of material properties $\boldsymbol{\beta }^{o}$ and its update $[\mathbf{A}^{o}]^{-1}\{ \mathbf{b}^{o} \}$ is less than the tolerance delta

$$ \frac{\left \| [\mathbf{A}^{o}]^{-1}\{ \mathbf{b}^{o} \} \right \| }{\left \| \boldsymbol{\beta }^{o} \right \| } < 10^{-6}~. $$

(55)

Although deriving analytically the convergence rate is not possible, we verified numerically in the following results that a quadratic convergence was obtained for cases where existence and uniqueness of the solution are guaranteed.

2.8 Examples of Hyperelastic Constitutive Models

In the following sections, we applied the VFM method and algorithm in Fig. 2 to determine the parameters of 3 hyperelastic strain energy potentials commonly used to describe collagenous tissues.

2.8.1 Neo-Hookean Model

For the Neo-Hookean constitutive model, the strain energy density function is

$$ \Phi = \frac{1}{2} \left [ \mu (\hat{I}_{1} - 3) + \kappa ( \mathrm{ln} J)^{2} \right ] , $$

(56)

where $\hat{I}_{1}=\textrm{Tr}\left ( \hat{\textbf{C}} \right )=J^{-2/3} \textrm{Tr}\left ( \textbf{C} \right )$.

The strain energy density function depends linearly on two unknown constitutive parameters denoted $\mu$ and $\kappa$. Then in an identification problem, the vector of unknown parameters would be: $\beta =(\mu ,\kappa )$.

The second Piola-Kirchhoff stress is written

$$ \textbf{S} = \kappa (\mathrm{ln} J) \mathbf{C}^{-1} + \mu \left ( J^{-2/3} \mathbf{I} -\frac{1}{3}\hat{I}_{1} \mathbf{C}^{-1} \right ) . $$

(57)

The associated Cauchy stress is written as

$$ \textbf{T} = J^{-1} \left [ \kappa (\mathrm{ln} J) \mathbf{I} + \mu ( \hat{\mathbf{B}} -\frac{1}{3}\hat{I}_{1} \mathbf{I} ) \right ] . $$

(58)

The sensitivity of the second Piola-Kirchhoff stress to each parameter is written

$$ \textstyle\begin{array}{c@{\quad }c@{\quad }c@{\quad }c@{\quad }l} \textbf{S}_{,\beta _{1}} & = & \frac{\partial \textbf{S} }{\partial \beta _{1}} & = & J^{-2/3} \mathbf{I} -\frac{1}{3}\hat{I}_{1} \mathbf{C}^{-1}, \\ \textbf{S}_{,\beta _{2}} & = & \frac{\partial \textbf{S} }{\partial \beta _{2}} & = & (\mathrm{ln} J) \mathbf{C}^{-1}. \end{array} $$

(59)

The sensitivity of the Cauchy stress to each parameter is written

$$ \textstyle\begin{array}{c@{\quad }c@{\quad }c@{\quad }c@{\quad }l} \textbf{T}_{,\beta _{1}} & = & \frac{\partial \textbf{T} }{\partial \beta _{1}} & = & J^{-1} ( \hat{\mathbf{B}} -\frac{1}{3}\hat{I}_{1} \mathbf{I} ) , \\ \textbf{T}_{,\beta _{2}} & = & \frac{\partial \textbf{T} }{\partial \beta _{2}} & = & J^{-1} ( \mathrm{ln} J) \mathbf{I} . \end{array} $$

(60)

Therefore,

$$ \mathbf{T}^{o} = \mu ^{o} \mathbf{T}_{,\beta _{1}}^{o} + \kappa ^{o} \mathbf{T}_{,\beta _{2}}^{o} , $$

(61)

with

$$ \textstyle\begin{array}{c@{\quad }c@{\quad }l} \textbf{T}_{,\beta _{1}}^{o} & = & (J^{o})^{-1} (\hat{\mathbf{B}}^{o} - \frac{1}{3}\hat{I}_{1}^{o} \mathbf{I} ) , \\ \textbf{T}_{,\beta _{2}}^{o} & = & (J^{o})^{-1} [\mathrm{ln} (J^{o}) ] \mathbf{I} . \end{array} $$

(62)

Moreover, we have

$$ \pmb{\pmb{\mathbb{K}}}^{o} = \mu ^{o} \pmb{\mathbb{K}}_{,\beta _{1}}^{o} + \kappa ^{o} \pmb{\mathbb{K}}_{,\beta _{2}}^{o} , $$

(63)

with

$$ \textstyle\begin{array}{c@{\quad }c@{\quad }l} \pmb{\mathbb{K}}_{,\beta _{1}}^{o} & = & -\frac{2}{3} J^{-2/3} \left ( \mathbf{I} \otimes (\mathbf{C}^{o})^{-1} + (\mathbf{C}^{o})^{-1} \otimes \mathbf{I} \right ) + \frac{2}{9} \hat{I}_{1} (\mathbf{C}^{o})^{-1} \otimes (\mathbf{C}^{o})^{-1} \\ & & + \frac{2}{3} \hat{I}_{1} (\mathbf{C}^{o})^{-1} \odot (\mathbf{C}^{o})^{-1} , \\ & & \\ \pmb{\mathbb{K}}_{,\beta _{2}}^{o} & = & (J^{o})^{-1} (\mathbf{C}^{o})^{-1} \otimes (\mathbf{C}^{o})^{-1} - (\mathrm{ln} J^{o}) (\mathbf{C}^{o})^{-1} \odot (\mathbf{C}^{o})^{-1} , \end{array} $$

(64)

where $\mathbf{C}^{-1} \odot \mathbf{C}^{-1}=- \frac{\partial \mathbf{C}^{-1}}{\partial \mathbf{C}}$ and ⊗ represents the dyadic multiplication symbol.

2.8.2 Mooney-Rivlin Model

For the Mooney-Rivlin constitutive model, the strain energy density function is

$$ \Phi = \frac{1}{2} \left [ \mu (\hat{I}_{1} - 3) + \alpha (\hat{I}_{2} - 3) + \frac{1}{2} \kappa (\mathrm{ln} J)^{2} \right ] , $$

(65)

where $\hat{I}_{2}=J^{-4/3} I_{2} = J^{-4/3}\textrm{Tr}\left ( \hat{\textbf{C}} \right )=J^{-2/3}\textrm{Tr}\left ( \textbf{C} \right )$.

The strain energy density function depends linearly on three unknown constitutive parameters denoted $\mu$, $\alpha$ and $\kappa$. Then in an identification problem, the vector of unknown parameters would be: $\beta =(\mu , \alpha , \kappa )$.

The second Piola-Kirchhoff stress is written

$$ \textbf{S} = \kappa (\mathrm{ln} J) \mathbf{C}^{-1} + \mu \left ( J^{-2/3} \mathbf{I} -\frac{1}{3}\hat{I}_{1} \mathbf{C}^{-1} \right ) + \alpha \left ( J^{-2/3} \hat{I}_{1} \mathbf{I} -\frac{2}{3}\hat{I}_{2} \mathbf{C}^{-1} - J^{-4/3} \mathbf{C} \right ) . $$

(66)

The associated Cauchy stress is written

$$ \textbf{T} = J^{-1} \left [ \kappa (\mathrm{ln} J) \mathbf{I} + \mu ( \hat{\mathbf{B}} -\frac{1}{3}\hat{I}_{1} \mathbf{I} ) + \alpha ( \hat{I}_{1} \hat{\mathbf{B}} -\frac{2}{3}\hat{I}_{2} \mathbf{I} - \hat{\mathbf{B}}^{2} ) \right ] . $$

(67)

The sensitivity of the second Piola-Kirchhoff stress to each parameter is written

$$ \textstyle\begin{array}{c@{\quad }c@{\quad }c@{\quad }c@{\quad }l} \textbf{S}_{,\beta _{1}} & = & \frac{\partial \textbf{S} }{\partial \beta _{1}} & = & J^{-2/3} \mathbf{I} -\frac{1}{3}\hat{I}_{1} \mathbf{C}^{-1} , \\ \textbf{S}_{,\beta _{2}} & = & \frac{\partial \textbf{S} }{\partial \beta _{2}} & = & J^{-2/3} \hat{I}_{1} \mathbf{I} -\frac{2}{3}\hat{I}_{2} \mathbf{C}^{-1} - J^{-4/3} \mathbf{C} , \\ \textbf{S}_{,\beta _{3}} & = & \frac{\partial \textbf{S} }{\partial \beta _{3}} & = & (\mathrm{ln} J) \mathbf{C}^{-1} . \end{array} $$

(68)

The sensitivity of the Cauchy stress to each parameter is written

$$ \textstyle\begin{array}{c@{\quad }c@{\quad }c@{\quad }c@{\quad }l} \textbf{T}_{,\beta _{1}} & = & \frac{\partial \textbf{T} }{\partial \beta _{1}} & = & J^{-1} ( \hat{\mathbf{B}} -\frac{1}{3}\hat{I}_{1} \mathbf{I} ) , \\ \textbf{T}_{,\beta _{2}} & = & \frac{\partial \textbf{T} }{\partial \beta _{2}} & = & J^{-1} ( \hat{I}_{1} \hat{\mathbf{B}} -\frac{2}{3}\hat{I}_{2} \mathbf{I} - \hat{\mathbf{B}}^{2} ) , \\ \textbf{T}_{,\beta _{3}} & = & \frac{\partial \textbf{T} }{\partial \beta _{3}} & = & J^{-1} ( \mathrm{ln} J) \mathbf{I} . \end{array} $$

(69)

Therefore,

$$ \mathbf{T}^{o} = \mu ^{o} \mathbf{T}_{,\beta _{1}}^{o} + \alpha ^{o} \mathbf{T}_{,\beta _{2}}^{o} + \kappa ^{o} \mathbf{T}_{,\beta _{3}}^{o} , $$

(70)

with

$$ \textstyle\begin{array}{c@{\quad }c@{\quad }l} \textbf{T}_{,\beta _{1}}^{o} & = & (J^{o})^{-1} \left [ \hat{\mathbf{B}}^{o} -\frac{1}{3}\hat{I}_{1}^{o} \mathbf{I} \right ] , \\ \textbf{T}_{,\beta _{2}}^{o} & = & (J^{o})^{-1} \left [ \hat{I}_{1}^{o} \hat{\mathbf{B}}^{o} -\frac{2}{3}\hat{I}_{2}^{o} \mathbf{I} - ( \hat{\mathbf{B}}^{o})^{2} \right ] , \\ \textbf{T}_{,\beta _{3}}^{o} & = & (J^{o})^{-1} [\mathrm{ln} (J^{o}) ] \mathbf{I} . \end{array} $$

(71)

Moreover, we have

$$ \pmb{\mathbb{K}}^{o} = \mu ^{o} \pmb{\mathbb{K}}_{,\beta _{1}}^{o} + \alpha ^{o} \pmb{\mathbb{K}}_{,\beta _{2}}^{o} + \kappa ^{o} \pmb{\mathbb{K}}_{,\beta _{3}}^{o} , $$

(72)

with

$$ \textstyle\begin{array}{rcl} \pmb{\mathbb{K}}_{,\beta _{1}}^{o} & = & -\frac{2}{3} J^{-2/3} \left ( \mathbf{I} \otimes (\mathbf{C}^{o})^{-1} + (\mathbf{C}^{o})^{-1} \otimes \mathbf{I} \right ) + \frac{2}{9} \hat{I}_{1} (\mathbf{C}^{o})^{-1} \otimes (\mathbf{C}^{o})^{-1} \\ & & + \frac{2}{3} \hat{I}_{1} (\mathbf{C}^{o})^{-1} \odot (\mathbf{C}^{o})^{-1} , \\ \pmb{\mathbb{K}}_{,\beta _{2}}^{o} & = & 2 J^{-4/3} \left ( \mathbf{I} \otimes \mathbf{I} - \mathbb{I} \right ) -\frac{4}{3} \hat{I}_{1} J^{-2/3} \left ( \mathbf{I} \otimes (\mathbf{C}^{o})^{-1} + (\mathbf{C}^{o})^{-1} \otimes \mathbf{I} \right ) \\ & & + \frac{8}{9} \hat{I}_{2} (\mathbf{C}^{o})^{-1} \otimes ( \mathbf{C}^{o})^{-1} + \frac{4}{3} \hat{I}_{2} (\mathbf{C}^{o})^{-1} \odot (\mathbf{C}^{o})^{-1} \\ & & + \frac{4}{3} J^{-4/3} \left [ \mathbf{C}^{o} \otimes (\mathbf{C}^{o})^{-1} + (\mathbf{C}^{o})^{-1} \otimes \mathbf{C}^{o} \right ] , \\ \pmb{\mathbb{K}}_{,\beta _{3}}^{o} & = & (J^{o})^{-1} (\mathbf{C}^{o})^{-1} \otimes (\mathbf{C}^{o})^{-1} - (\mathrm{ln} J^{o}) (\mathbf{C}^{o})^{-1} \odot (\mathbf{C}^{o})^{-1} , \end{array} $$

(73)

where $\mathbb{I}$ is the fourth order identity tensor.

2.8.3 Veronda-Westmann Model

In the two previous constitutive models, the strain energy density function depends linearly on the material properties. In this subsection, the Veronda-Westmann model [43] which involves an exponential function of a material property introduced. The strain energy density function is written such as

$$ \Phi = \frac{c_{1}}{2} \left ( e^{c_{2}(\hat{I}_{1} - 3)} -1 \right ) - \frac{c_{1} c_{2}}{2} (\hat{I}_{2} - 3) + \frac{1}{2} \kappa ( \mathrm{ln} J)^{2} , $$

(74)

where $c_{1}$, $c_{2}$ and $\kappa$ are material properties. Note that $c_{2}$ is a parameter controlling the nonlinear behavior of the Veronda-Westmann solid, and $\mu = c_{1} c_{2}$ is the shear modulus controlling the linear behavior of the Veronda-Westmann solid. In the following, we rewrite the model such as

$$ \Phi = \frac{\mu }{2c_{2}} \left ( e^{c_{2}(\hat{I}_{1} - 3)} -1 \right ) -\frac{\mu }{2} (\hat{I}_{2} - 3) + \frac{1}{2} \kappa ( \mathrm{ln} J)^{2} . $$

(75)

The strain energy density function depends linearly on two unknown constitutive parameters denoted $\mu$ and $\kappa$, and non linearly on the unknown parameter $c_{2}$. In the identification problem, the vector of unknown parameters would be written: $\beta =(\mu ,c_{2},\kappa )$.

The second Piola-Kirchhoff stress is written

$$\begin{array}{c}\mathbf{\text{S}}=\kappa (\ln J){\mathbf{C}}^{-1}+\mu [(J^{-2/3}\mathbf{I}-\frac{1}{3}{\hat{I}}_1{\mathbf{C}}^{-1})e^{c_2({\hat{I}}_1-3)}\hfill \\ \hfill +(J^{-2/3}{\hat{I}}_1\mathbf{I}-\frac{2}{3}{\hat{I}}_2{\mathbf{C}}^{-1}-J^{-4/3}\mathbf{C})].\end{array}$$

(76)

The associated Cauchy stress is written

$$ \textbf{T} = J^{-1} \left [ \kappa (\mathrm{ln} J) \mathbf{I} + \mu \left ( (\hat{\mathbf{B}} -\frac{1}{3}\hat{I}_{1} \mathbf{I} )e^{c_{2}( \hat{I}_{1} - 3)} + ( \hat{I}_{1} \hat{\mathbf{B}} -\frac{2}{3} \hat{I}_{2} \mathbf{I} - \hat{\mathbf{B}}^{2} ) \right ) \right ] . $$

(77)

The sensitivity of the second Piola-Kirchhoff stress to each parameter is written

$$ \textstyle\begin{array}{rcl} \textbf{S}_{,\beta _{1}} & = & \frac{\partial \textbf{S} }{\partial \beta _{1}} = \left ( J^{-2/3} \mathbf{I} -\frac{1}{3}\hat{I}_{1} \mathbf{C}^{-1} \right )e^{c_{2}( \hat{I}_{1} - 3)} \\ & & + \left ( J^{-2/3} \hat{I}_{1} \mathbf{I} -\frac{2}{3}\hat{I}_{2} \mathbf{C}^{-1} - J^{-4/3} \mathbf{C} \right ) , \\ \textbf{S}_{,\beta _{2}} & = & \frac{\partial \textbf{S} }{\partial \beta _{2}} = \mu (\hat{I}_{1} - 3)\left ( J^{-2/3} \mathbf{I} -\frac{1}{3}\hat{I}_{1} \mathbf{C}^{-1} \right )e^{c_{2}(\hat{I}_{1} - 3)} , \\ \textbf{S}_{,\beta _{3}} & = & \frac{\partial \textbf{S} }{\partial \beta _{3}} = (\mathrm{ln} J) \mathbf{C}^{-1} . \end{array} $$

(78)

The sensitivity of the Cauchy stress to each parameter is written

$$ \textstyle\begin{array}{rcl} \textbf{T}_{,\beta _{1}} & = & \frac{\partial \textbf{T} }{\partial \beta _{1}} = J^{-1} \left [ ( \hat{\mathbf{B}} -\frac{1}{3}\hat{I}_{1} \mathbf{I} )e^{c_{2}(\hat{I}_{1} - 3)} + ( \hat{I}_{1} \hat{\mathbf{B}} -\frac{2}{3}\hat{I}_{2} \mathbf{I} - \hat{\mathbf{B}}^{2} ) \right ] , \\ \textbf{T}_{,\beta _{2}} & = & \frac{\partial \textbf{T} }{\partial \beta _{2}} = J^{-1} \mu ( \hat{I}_{1} - 3) (\hat{\mathbf{B}} -\frac{1}{3}\hat{I}_{1} \mathbf{I} )e^{c_{2}( \hat{I}_{1} - 3)} , \\ \textbf{T}_{,\beta _{3}} & = & \frac{\partial \textbf{T} }{\partial \beta _{3}} = J^{-1} ( \mathrm{ln} J) \mathbf{I} . \end{array} $$

(79)

Then we approximate $\mathbf{T}^{o}$ such as

$$ \mathbf{T}^{o} = \mu ^{o} \mathbf{T}_{,\beta _{1}}^{o} + c_{2}^{o} \mathbf{T}_{,\beta _{2}}^{o} + \kappa ^{o} \mathbf{T}_{,\beta _{3}}^{o} , $$

(80)

with

$$ \textstyle\begin{array}{rcl} \textbf{T}_{,\beta _{1}}^{o} & = & (J^{o})^{-1} \left [ ( \hat{\mathbf{B}}^{o} -\frac{1}{3}\hat{I}_{1}^{o} \mathbf{I} )e^{c_{2}^{o}( \hat{I}_{1}^{o} - 3)} + ( \hat{I}_{1}^{o} \hat{\mathbf{B}}^{o} - \frac{2}{3}\hat{I}_{2}^{o} \mathbf{I} - (\hat{\mathbf{B}}^{o})^{2} ) \right ] , \\ \textbf{T}_{,\beta _{2}}^{o} & = & (J^{o})^{-1} \mu ^{o} (\hat{I}_{1}^{o} - 3) (\hat{\mathbf{B}}^{o} -\frac{1}{3}\hat{I}_{1}^{o} \mathbf{I} )e^{c_{2}^{o}( \hat{I}_{1}^{o} - 3)} , \\ \textbf{T}_{,\beta _{3}}^{o} & = & (J^{o})^{-1} [\mathrm{ln} (J^{o}) ] \mathbf{I} . \end{array} $$

(81)

Moreover, we also approximate

$$ \pmb{\mathbb{K}}^{o} = \mu ^{o} \pmb{\mathbb{K}}_{,\beta _{1}}^{o} + c_{2}^{o} \pmb{\mathbb{K}}_{,\beta _{2}}^{o} + \kappa ^{o} \pmb{\mathbb{K}}_{, \beta _{3}}^{o} ~, $$

(82)

with

$$ \textstyle\begin{array}{rcl} \pmb{\mathbb{K}}_{,\beta _{1}}^{o} & = & \left [ -\frac{2}{3} J^{-2/3} \left ( \mathbf{I} \otimes (\mathbf{C}^{o})^{-1} + (\mathbf{C}^{o})^{-1} \otimes \mathbf{I} \right ) + \frac{2}{9} \hat{I}_{1} (\mathbf{C}^{o})^{-1} \otimes (\mathbf{C}^{o})^{-1} \right . \\ & & + \left . \frac{2}{3} \hat{I}_{1} (\mathbf{C}^{o})^{-1} \odot ( \mathbf{C}^{o})^{-1} \right ] e^{c_{2}^{o}(\hat{I}_{1}^{o} - 3)} + 2 J^{-4/3} \left ( \mathbf{I} \otimes \mathbf{I} - \mathbb{I} \right ) \\ & & -\frac{4}{3} \hat{I}_{1} J^{-2/3} \left ( \mathbf{I} \otimes ( \mathbf{C}^{o})^{-1} + (\mathbf{C}^{o})^{-1} \otimes \mathbf{I} \right ) + \frac{8}{9} \hat{I}_{2} (\mathbf{C}^{o})^{-1} \otimes ( \mathbf{C}^{o})^{-1} \\ & & + \frac{4}{3} \hat{I}_{2} (\mathbf{C}^{o})^{-1} \odot (\mathbf{C}^{o})^{-1} + \frac{4}{3} J^{-4/3} \left [ \mathbf{C}^{o} \otimes (\mathbf{C}^{o})^{-1} + (\mathbf{C}^{o})^{-1} \otimes \mathbf{C}^{o} \right ] \\ & & + 2 c_{2}^{o} \left [ J^{-4/3} \mathbf{I} \otimes \mathbf{I} - \frac{1}{3} J^{-2/3} \hat{I}_{1}^{o} \left ( (\mathbf{C}^{o})^{-1} \otimes \mathbf{I} + \mathbf{I} \otimes (\mathbf{C}^{o})^{-1} \right ) \right . \\ & & \left . + \frac{1}{9} (\hat{I}_{1}^{o})^{2} (\mathbf{C}^{o})^{-1} \otimes (\mathbf{C}^{o})^{-1} \right ] e^{c_{2}^{o}(\hat{I}_{1}^{o} - 3)} , \\ \pmb{\mathbb{K}}_{,\beta _{2}}^{o} & = & \mu ^{o}(\hat{I}_{1}^{o} - 3) e^{c_{2}^{o}(\hat{I}_{1}^{o} - 3)} \left [ -\frac{2}{3} J^{-2/3} \left ( \mathbf{I} \otimes (\mathbf{C}^{o})^{-1} + (\mathbf{C}^{o})^{-1} \otimes \mathbf{I} \right ) \right . \\ & & \left . + \frac{2}{9} \hat{I}_{1} (\mathbf{C}^{o})^{-1} \otimes ( \mathbf{C}^{o})^{-1} + \frac{2}{3} \hat{I}_{1} (\mathbf{C}^{o})^{-1} \odot (\mathbf{C}^{o})^{-1} \right ] \\ & & + 2 \mu ^{o} \left (1 + c_{2}^{o}(\hat{I}_{1}^{o} - 3) \right ) e^{c_{2}^{o}( \hat{I}_{1}^{o} - 3)} \left [ J^{-4/3} \mathbf{I} \otimes \mathbf{I} \right . \\ & & \left . -\frac{1}{3} J^{-2/3} \hat{I}_{1}^{o} \left ( (\mathbf{C}^{o})^{-1} \otimes \mathbf{I} + \mathbf{I} \otimes (\mathbf{C}^{o})^{-1} \right ) + \frac{1}{9} (\hat{I}_{1}^{o})^{2} (\mathbf{C}^{o})^{-1} \otimes ( \mathbf{C}^{o})^{-1} \right ] , \\ \pmb{\mathbb{K}}_{,\beta _{3}}^{o} & = & (J^{o})^{-1} (\mathbf{C}^{o})^{-1} \otimes (\mathbf{C}^{o})^{-1} - (\mathrm{ln} J^{o}) (\mathbf{C}^{o})^{-1} \odot (\mathbf{C}^{o})^{-1} . \end{array} $$

(83)

2.9 Case Study for Verification

Using the previous formulas derived for the Neo-Hookean, the Mooney-Rivlin and the Veronda-Westmann models, it is straightforward to implement the novel VFM method and apply it to solve identification problems based on full-field data. Note also that the approach also extends naturally to more complex strain energy density functions involving exponential functions with invariants $I_{4}$ and $I_{6}$ [28] for instance. In the next section we show results for the 3 constitutive models, which were previously introduced, and discuss the convergence and the feasibility of the proposed approach in the case of compression and tension tests onto 2 simple geometries (Fig. 3):

• 1. Compression is applied on a cube of 1 cm edge which is fully fixed on its bottom face (Fig. 3a). For the sake of verification of the approach, the cube is simply discretized uniformly by $4 \times 4 \times 4$ elements.

  Fig. 3

  

  The geometry and finite-element meshes used to verify the implementation of the novel VFM approach

• 2. We also consider another geometric model as presented in (Fig. 3b) where we apply tension on the top face and fix the bottom face.

In both cases, full-field displacement measurements were simulated by solving finite-element models using the open source software FEBio [44]. For each case, the full-field displacement measurements were simulated with a set of target material properties, named the target, and the proposed VFM approach was employed to recover the target for the Neo-Hookean, the Mooney-Rivlin and the Veronda-Westmann model, using different initialization values to assess the convergence of the method.

2.10 Case Study for Validation

The proposed VFM method will be used in this case study to determine the material properties of the LC, a connective tissue structure in the ONH of great interest to researchers studying development and progression of glaucoma [45]. The LC in humans is approximately 1.5 mm - 2.0 mm in diameter and 450 microns thick and is located in the back of the eye. The mechanics of the LC is thought [46] to play an important role in mediating the progressive vision loss associated with glaucoma, a degenerative disease that is a leading cause of blindness worldwide [47]. As shown in Fig. 4, the LC is situated beneath the prelaminar neural tissue (PLNT). We further split the LC into an anterior region (ALC), 250 ${\mu }$m posterior to the anterior posterior surface, and a posterior region (PLC), which includes the remainder of the imageable volume of the LC.

Fig. 4

(a) Schematic showing the solid model of the optical nerve head (ONH). The anterior (top) face is subjected to the intra-ocular pressure. The posterior (bottom) face is in direct contact with the optic nerve. (b) Cross sectional view of the ONH solid model showing the 3 regions with different material properties: the prelaminar neural tissue (PLNT), the anterior lamina cribosa (ALC) and the posterior lamina cribosa (PLC). (c) Displacement field caused by reducing the intraocular pressure of 10 mmHg (colorbar in mm)

Spectral domain OCT imaging (Heildelberg Engineering) was applied to acquire 24 radial scans centered about the ONH of the left eye of a glaucoma patient. The OCT imaging was performed at Johns Hopkins University’s Wilmer Eye Institute in the Glaucoma Center of Excellence, and was approved by the appropriate Institutional Review Board. The following structural features were marked in the 24 radial scans to segment the tissue structures of the ONH as described in Midgett et al. [48]: Bruch’s membrane opening, the anterior boundary of the PLNT, the anterior LC surface, and the boundary of the imageable volume below the anterior LC surface. The manual marking were imported into Cubit (Coreform, Orem, UT, USA) to construct surface geometries using closed splines. The PLNT, ALC, and PLC volumes were defined by extruding a cylinder from Bruch’s membrane posteriorly to intersect the anterior PLNT surface and the anterior LC surface (for the PLNT), the anterior ALC surface and a surface positioned 250 microns posterior to the anterior ALC surface (for the ALC), and that posterior surface and a surface marking the end of the imageable volume of the LC (see Fig. 4a). The final solid volume was meshed in Cubit linear 4-node tetrahedral elements and exported into FEBio, where the linear elements were converted to 10-node quadratic tetrahedral elements [49] to avoid mesh locking and improve accuracy. The compressible Neo-Hookean constitutive model (Eq. (56)) was chosen for all three materials.

To validate the present VFM method, we simulated a displacement field by solving the forward problem with target parameter values, which are reported in Table 1. In the forward problem, zero displacement boundary conditions were applied to the posterior and lateral surfaces. To account for the 10 mmHg pressure decrease, a pressure boundary condition was applied to the ONH surface, with a magnitude of −10 mmHg. The induced displacement is shown in Fig. 4c.

Table 1 Results obtained with the novel VFM for identifying compressibility modulus values of the 3 separate regions in the ONH using simulated data

Eventually the proposed VFM approach was employed to process the simulated displacement fields and recover the target compressibility modulus values for each separate material.

3 Results

3.1 Neo-Hookean Model

In this section, we tested the feasibility of the proposed VFM approach for the compressible neo-Hookean model using the 10-mm-edge cubic specimen under compression shown in Fig. 3a. The cube was discretized uniformly by $4 \times 4 \times 4$ elements. We applied a deformation of 2 mm to the cubic specimen. The target values of material properties were: Young’s modulus $\bar{E} = 10~{\mathrm{{MPa}}}$ and Poisson’s ratio $\bar{\nu }=0.3$, which corresponds to $\mu =\frac{E}{2(1+\nu )}=3.85~\mathrm{{MPa}}$ and $\kappa =\frac{E}{3(1-2\nu )}=8.33~\mathrm{{MPa}}$. We choose three different initial guess pairs $({E^{o}} = 40~{\mathrm{{MPa}}}, {\nu ^{o}} = 0.45)$, $({E^{o}} = 5~{\mathrm{{MPa}}}, {\nu ^{o}} = 0.45)$ and $({E^{o}} = 40~{\mathrm{{MPa}}}, {\nu ^{o}} = 0.15)$. The convergence plots for each case are shown in Fig. 5 for the objective function value and in Fig. 6 for the evolution of the obtained parameter values. The material parameters were recovered regardless of the initial guesses. Convergence was reached after 4 to 6 iterations indicating the quadratic convergence of the new VFM approach.

Fig. 5

Evolution of the objective function relative value $||\textbf{u}^{o}-\tilde{ \textbf{u}}||/||\tilde{ \textbf{u}}||$ (log-log representation) for the identification of each unknown parameter of the neo-Hookean model with the novel VFM approach, using initial guesses ${E^{o}} = 40~{\mathrm{{MPa}}}$, ${\nu ^{o}} = 0.45$ in (a), ${E^{o}} = 5~{\mathrm{{MPa}}}$, ${\nu ^{o}} = 0.45$, in (b) and ${E^{o}} = 40~{\mathrm{{MPa}}}$, ${\nu ^{o}} = 0.15$ in (c)

Fig. 6

Convergence plots for the identification of each unknown parameter of the Neo-Hookean model with the novel VFM approach, using initial guesses ${E^{o}} = 40~{\mathrm{{MPa}}}$, ${\nu ^{o}} = 0.45$ in (a) and (b), ${E^{o}} = 5~{\mathrm{{MPa}}}$, ${\nu ^{o}} = 0.45$ in (c) and (d) and ${E^{o}} = 40~{\mathrm{{MPa}}}$, ${\nu ^{o}} = 0.15$ in (e) and (f)

3.2 Mooney-Rivlin Model

Then, we tested the feasibility of the proposed VFM approach for the Mooney-Rivlin model, which includes three unknown material parameters. In this case, the material property vector was denoted $\beta =[\mu ,\alpha ,\kappa ]$. Note that when $\alpha =0$, the Mooney-Rivlin model simplifies into the neo-Hookean model. For the verification of the VFM, we used again the cubic model under compression. The target material properties were set to $\bar{\beta }=[5~\mathrm{{MPa}},10~\mathrm{{MPa}},10~\mathrm{{MPa}]}$. We chose two different sets of initial guess to evaluate the material identification problem. Convergence plots for each material parameters are shown in Fig. 7. The material parameters were recovered regardless of the initial guesses. Convergence was reached after 4 to 6 iterations indicating the quadratic convergence of the new VFM approach with the Mooney-Rivlin model.

Fig. 7

Convergence plots of the novel VFM approach for each parameter of the Mooney-Rivlin model using two sets of initial guesses: (a, b and c): $\beta ^{o}=[4~\mathrm{{MPa}},6~\mathrm{{MPa}},30~\mathrm{{MPa}]}$ and (d, e and f) $\beta ^{o}=[20~\mathrm{{MPa}},20~\mathrm{{MPa}},50~\mathrm{{MPa}]}$

3.3 Veronda-Westmann Model

Usually, to estimate the material parameters of the Veronda-Westmann model, displacement data with multiple loading stages are required due to the exponential term. In this work, we attempted to use only one loading stage to estimate the material parameters in the Veronda-Westmann model using the novel VFM approach. For that, we used the geometric model on the right hand side of Fig. 3. We applied tension on the top face and fix the bottom face. A deformation of 2 mm was applied. The target material parameters were $\bar{c}_{1}=0.1~\mathrm{{MPa}}$, $\bar{c}_{2}=10$, $\bar{\kappa }=10~\mathrm{{MPa}}$. We used two different initial guesses. The initial guesses of the material parameters were $c_{1}^{0}=0.2~\mathrm{{MPa}}$, $c_{2}=30$, $\kappa =20~\mathrm{{MPa}}$. Another initial guess is set to $c_{1}^{0}=0.5~\mathrm{{MPa}}$, $c_{2}=5$, $\kappa =5~\mathrm{{MPa}}$

The convergence performance of each material parameters is shown in Fig. 8. We observed that 12 iterations were needed this time to reach convergence. Although the proposed VFM approach was capable of estimating material properties of the Veronda-Westmann solid using only one loading stage with a quadratic convergence, the convergence rate was slower due to the nonlinearity of the constitutive equations.

Fig. 8

Convergence plots of the novel VFM approach for each parameter of the Veronda-Westmann model with initial guesses $c_{1}^{0}=0.2~\mathrm{{MPa}}$, $c_{2}=30$, $\kappa =20~\mathrm{{MPa}}$ in (a), (b) and (c), and initial guesses $c_{1}^{0}=0.5~\mathrm{{MPa}}$, $c_{2}=5$, $\kappa =5~\mathrm{{MPa}}$ in (d), (e) and (d)

3.4 Sensitivity to the Compressibility Parameters

Most biological tissues are nearly or truly incompressible. Thus, the sensitivity of the proposed method to the compressibility should be investigated. For incompressible materials, fewer material properties have to be identified compared to the associated compressible constitutive model. Thus, we considered the Neo-Hookean model as an example and focused on compressible cases with Poisson’s ratio $\nu =0.45$ and Young’s modulus $\bar{E} = 10~\mathrm{MPa}$, which corresponds to $\mu =\frac{E}{2(1+\nu )}=3.45~\mathrm{MPa}$ and $\kappa =\frac{E}{3(1-2\nu )}=33.33~\mathrm{MPa}$.

We chose two different initial guess pairs ($E^{o} = 5~\mathrm{MPa}$, $\nu ^{o} = 0.499$), ($E^{o} = 20~\mathrm{MPa}$ and $\nu ^{o} = 0.499$). The convergence plots for each case are shown in Fig. 9. We observed that it took 8 iterations for convergence, slightly more than the more compressible cases where $\nu$ was taken equal to 0.3.

Fig. 9

Convergence plots for the identification of each unknown parameter of the neo-Hookean model with Poisson’s ratio $\nu =0.45$, using initial guesses $E^{o} = 5~\mathrm{{MPa}}$, $\nu ^{o} = 0.499$ in (a) and (b), and $E^{o} = 20~\mathrm{{MPa}}$, $\nu ^{o} = 0.499$ in (c) and (d)

3.5 Case Study for Validation

The patient-specific geometry of the ONH, which was reconstructed using OCT data (Fig. 4), was split into 3 parts with different compressibility moduli. It was assumed that each part had the same shear modulus $\mu =0.2$ MPa. Identifying the 3 unknown compressibility moduli with the VFM requires to define 3 independent virtual fields across the volume of interest. As any kinematically admissible virtual displacement field is a possible option, there are an infinite number of choices for these virtual fields. However, only very specific virtual fields ensure that the final system of equations is well conditioned and it may be cumbersome to determine these virtual fields by trial and errors. Here we applied the new VFM approach, where systems of equations are automatically built, without trial and error verifications, on the patient-specific geometry of the lamina cribrosa. As our VFM approach is iterative, values need to be chosen for initializing the compressibility moduli at the first iteration. We drew randomly these initial values within the range [1–5] MPa. The VFM always converged towards the target moduli. The largest error was obtained when the initial values were 5 MPa for the 3 compressibility moduli. Results obtained in that case are reported in Table 1. Errors below 1% can be reached after only 4 iterations of the new VFM approach.

4 Discussion

In this paper, we proposed a novel general framework to solve identification problems in hyperelasticity with the VFM. For that, we employed the finite-element method and generated the virtual fields by solving novel equations derived in this paper. The main advantage of the new approach for generating the virtual fields is that it can be applied to any kind of geometric shapes.

The VFM was previously applied to identify multiaxial hyperelastic properties of murine dissecting aneurysm samples [25, 29, 30] or hyperelastic properties in human ocular tissues [23]. This permitted significant progress especially in the characterization of heterogeneous material properties of lesions exhibiting complex morphologies, with different regions characterized by localized changes in tissue composition, microstructure, and properties. This was rendered possible by combining extension-distension data with full-field multimodality measurements of wall strain and thickness to inform the inverse material characterization using the virtual fields method.

Latest key advances were the use of a digital volume correlation data [23, 25] that allowed for characterization of properties in the bulk of soft tissues. However, despite these progresses, generic rules had never been defined for the choice of virtual fields. For instance, in [23], analytical expressions of virtual fields were written in the cylindrical coordinate system using the $arctan$ function. In [29], virtual fields were also defined analytically, permitting to yield expressions generalizing the Laplace’s law for pressurized membranes.

The main difficulty of these previous studies was to ensure that the choice of virtual fields can provide enough equations to identify all the unknown constitutive parameters. A subsequent difficulty is to build systems of equations which are well conditioned. Moreover, in these previous studies focused on nearly incompressible solids, the identification of compressibility parameters was avoided by selecting virtual fields which would filter out the effects of the hydrostatic pressure in the principle of virtual power. This was generally required as the volume changes could not be measured with enough accuracy by optical techniques, causing too much uncertainty on the local values of the hydrostatic pressure anyhow. Generating virtual fields that can filter out the hydrostatic pressure was never easy with analytic functions. The use of a finite-element implementation of the VFM simplified this task [50] but the novel approach presented in the current paper further generalizes this for any situation. There is a price to pay though as the approach proposed in the current paper is iterative, requiring to choose a first guess of the unknown material parameters and repeating the identification a number of time. Fortunately, we found that the convergence of this iterative approach was always quadratic.

It is important to emphasize that any kinematically admissible virtual displacement with a non-zero volume change is always an option to identify material properties of compressible materials with the VFM. However, only very specific virtual fields can ensure that the final system of equations is well conditioned and it can be cumbersome to determine these virtual fields by trial and errors. For complex geometries and when there are several parameters to identify, as in the ONH problem shown in this paper, it is convenient to use our novel VFM approach, which automatically generates the virtual fields as solutions of Eq. (42). Although these specific virtual fields satisfy the equilibrium conditions written in Eq. (42), virtual fields in general do not have to satisfy such equations.

To test the feasibility of the novel VFM approach, we first utilized the simplest hyperelastic model: the Neo-Hookean model. This model remains very popular for many biological tissues as for instance ocular tissues [23]. The novel VFM approach, which is iterative, showed a fast and quadratic convergence for neo-Hookean identification. Additionally, we employed the proposed approach for more sophisticated model, such as the Mooney-Rivlin model and the Veronda-Westmann model. The estimation results demonstrated that the proposed method is able to identify the hyperelastic parameters of such models, and even that constitutive equations involving an exponential could still be identified with high accuracy even using a single loading stage.

We then applied the proposed VFM method to the clinically-relevant challenge of extracting material properties of ONH tissues from full-field OCT data. Modeling the LC in living subjects is difficult, as this tissue is strongly inhomogeneous and exhibits very large variations in strain. In addition, OCT imaging does not allow the user to definitively probe the tissue itself and often cannot fully resolve tissue boundaries, complicating inverse finite element methods which rely upon strictly defined geometry to delineate tissues and material models. Results obtained on simulated data are very promising, showing that our VFM approach can remove these difficulties by providing material properties in different regions of the ONH. Future use on experimental measurements will allow to examine regional stiffness variations, which otherwise may have been neglected.

Despite these significant improvements of the VFM, the content of the current paper remains mostly theoretical, introducing the concepts, verifying them for simple cases and validating the approach on a first case study. The main limitation of the current approach is that it was developed in Matlab. It needs to be fully incorporated into a finite-element package to optimize computational costs and constitute a systematic alternative to FEMU techniques [10, 11] for all identification problems in hyperelasticity based on full-field displacement measurements.

Another limitation of the VFM is that it requires full-field data to achieve parameter identification. Although it performs very well when these full-field data are available, other methods based on the FEMU approach [10] do not have this limitation. For many years, full-field data were mostly available in two-dimensional bodies (membranes, shells) on which displacement fields could be measured using the DIC technique [22, 27, 31]. However, bulk measurements of displacements fields using digital volume correlation and imaging techniques such as Optical Coherence Tomography or Magnetic Resonance Imaging are becoming commonplace in soft tissues [23–25, 32–35], offering more and more applications for the VFM in the biomedical field.

5 Conclusions

In this paper we presented a novel framework of the virtual fields method based on finite-element implementation and automatic generation of virtual fields. The proposed approach makes a great improvement in the theoretical aspects of the VFM as the choice of virtual fields is usually critical in the VFM. Future work will be focused on applying this finite-element based VFM into more complex cases with potential mechanobiological and clinical applications.