2.1 Theoretical Framework

Let

be the

*reference configuration* of a body in a 3D Euclidean space

at the time

*t*_{r} its geometry was captured experimentally. In

, the body is supposed to be subjected to an external load (e.g. pressure) that is in equilibrium with the internal stress state. From

the body deforms at a certain time

*t* to its

*actual configuration*
due to e.g. changes in the external loading conditions (). At the time

*t*_{0} the body was in its

*unloaded configuration*
in which the external load (pressure) was zero and the body was unstressed. Residual stresses are not considered in the present work.

Let

**x**_{r} be the position vector of a point

in

. The motion of this point to its unloaded

**X** or actual position

**x** may be expressed as a function of time and the reference position

The two deformation gradients that are related to these motions are defined as

The total deformation gradient related to the motion from the unloaded to the actual configuration can be decomposed into these two tensors (

de Putter et al. 2007;

Gee et al. 2010)

Our target is to compute the motion

**x**(

**x**_{r},

*t*) =

**x**_{r} +

**u**_{t}(

**x**_{r},

*t*) of the body form the loaded reference configuration

to its actual configuration

(1)

_{2} by using the finite-element method. This can be done by solving the weak form of the balance equation according to the principal of virtual work. In the absence of body forces, the virtual work reads with respect to the unloaded material configuration

with the external load vector

**t**_{0} and the variation of the right Cauchy-Green tensor

*δ*C =

*δ*(F

^{T}F). Let us assume that the mapping from the unloaded to the reference configuration is known such that

*δ*F

_{r} = 0 and

*δ***x**_{r} =

**0**. In this case and by considering the multiplicative split (

3), the variation of the right Cauchy-Green tensor can be restated to

Equation (5) together with the identities

**t**_{0}d

*A*_{0} =

**t**_{r}d

*A*_{r} and detF

_{r} = d

*V*_{r}/d

*V*_{0} can be used to reformulate the virtual work (

4) with respect to the reference configuration to

where Σ represents the reference stress tensor at

.

Assuming a hyperelastic constitutive response throughout the deformation history of the body and according to (

6), the reference stress can be expressed as a function of total F and the pre-existing deformation gradient F

_{r}Consequently, the deformation

**u**_{t} can only be computed from (

6) if the pre-existing deformation gradient F

_{r} is known. One can also see that the explicit calculation of the unloaded configuration

and the related motion

**X**(

**x**_{r},

*t*_{0}) (1)

_{1} are not needed to solve (

6). However, most of the prestressing approaches available are based on the backward solution of the motion

**X**(

**x**_{r},

*t*_{0}). Note that the Cauchy stress relates to the reference stress through

.

2.2 Forward Incremental Prestressing Method

The presented method is based on the

*modified updated Lagrangian formulation* proposed by

Gee et al. (2010). The basic idea of

Gee et al. (2010) was to estimate the pre-existing stress/strain state by computing the pre-existing deformation gradient F

_{r} (2)

_{1} with a forward calculation and not by calculating the backward motion

**X**(

**x**_{r},

*t*_{0}) (1)

_{1}.

Let

be a configuration in the close neighborhood of the experimentally obtained reference configuration

at the increment

*i* (). In the following we will call

the

*neighboring reference configuration.* At

, the total deformation gradient is split into pre-existing deformations F

_{r} representing the mapping from the unloaded to the true reference configuration and subsequent deformations

to the neighboring reference configuration

In

, the stress state Σ (F

^{i}, F

_{r}) ≠

**0** is required to be in equilibrium with the known external load

**t**_{r}. Following the assumptions made in the previous subsection, the weak form of the balance equation at an increment

*i* can be derived similar to (

6)

where

.

The aim of the forward incremental prestressing method is to incrementally reduce the subsequent deformations

until the neighboring reference configuration

coincides with the true reference configuration

. Let us assume that a limit state

*i* → ∞ exists, where the subsequent deformations vanish

In this case,

equation (9) can only be true if a non-trivial deformation gradient F

_{r} ≠ I exists that results in a non-trivial stress tensor Σ(F

_{r}, F

_{r}) ≠ 0. At this limit state the total deformation gradient (

8) and the stress tensor, consequently, represent the pre-existing strain/stress state

To achieve the limit state (

10) we use relation (

11)

_{1} to motivate the incremental update of the pre-existing deformation gradient in accordance to

Gee et al. (2010)If

and

are in a close neighborhood, the iterative application of solving (

9) and subsequently updating F

_{r} according to (

12) will lead to the limit state outlined in (

10) and (

11).

The forward incremental prestressing procedure is summarized in . To ease the understanding of the forward incremental prestressing method we illustrate the procedure in the following: First, the pre-existing deformation gradient F

_{r} = I is initialized. In the first increment, the pre-existing load

**t**_{r} is applied to the reference configuration as it would be the unloaded configuration

=

. A usual iterative finite element calculation is performed and the body deforms to the neighboring reference configuration

(see ), where external and internal forces are at equilibrium based on (

9). During the first increment, the total deformation gradient does not include pre-existing deformation F

_{r} = 0 but consists solely of subsequent deformations

.

| **Table 1**Forward incremental prestressing procedure. **t**_{r} the known pre-existing external load vector, which can be applied in multiple load steps *n* = 1,2, … *m; tol* anticipated tolerance of the prestressing procedure; *i*=1,2,… cumulated number of (more ...) |

At the end of the first increment, the pre-existing deformation gradient is explicitly updated according to (

12). Consequently, the total deformation gradient of the second increment

is now partially related to pre-existing deformations F

_{r} and subsequent deformations

. As F

^{2} is utilized to compute the constitutive equation Σ(F

^{2}, F

_{r}) to solve (

9), the pre-existing deformations F

_{r} will contribute to the internal forces. In turn, the subsequent deformations

need to be smaller than

to lower the contribution of

to the internal forces, such that equilibrium can be achieved at

. At the end of the second increment, the pre-existing deformation gradient is again updated according to (

12). Further iterations of solving (

9) followed by the update of the pre-existing deformations (

12) will eventually lead to vanishing subsequent deformations as anticipated in (

10). The prestressing procedure is deemed converged once the subsequent displacements are below a given tolerance

.

One can picture the iterative procedure outlined above as the incremental decrease of the distance between the neighboring and the true reference configuration while the unloaded configuration related to the neighboring reference configuration virtually moves from the true reference configuration

=

to the true unloaded configuration

=

(see ). However, the forward incremental prestressing method does not require the computation of the unloaded configuration

or the backward motion

**X**(

**x**_{r},

*t*_{0}). This is a major advantage of prestressing methods that are based on the forward calculations of F

_{r} such as the present and the MULF method (

Gee et al. 2010) compared to inverse elasostatic approaches.

In the original MULF method (

Gee et al. 2010), the neighboring reference configuration

was not updated during the prestressing process. Instead, the Jacobian matrix was updated to perform an equivalent step. Consequently, equilibrium between the pre-existing external load and the internal stress state was not achieved with respect to the true reference configuration

but with respect to an ‘unconverged’ neighboring reference configuration

. Due to this inconsistency, the results of the original prestressing method depend on the load step size and the estimation of the prestressed state can only be considered as an approximation. (

Gee et al. 2010) noted that the MULF method overestimates the pre-existing stresses. This is not the case for the present method as

is incrementally updated until it coincides with the true reference configuration

with respect to a given tolerance (

). The vector field

is consistent with the neighboring reference configuration

, where equilibrium between external and internal forces is achieved at the end of each prestressing increment

*i*. Accordingly, the converged numerical results of the present prestressing method are independent of the applied load step size. However, the present method requires reasonably small load increments during prestressing to guarantee convergence similar to explicit integration schemes.

Note that the boundary value problem defined in (

9) is practically identical to (

6) except for the incremental update of the pre-existing deformation gradient F

_{r} through (

12). Consequently, the implementation of the presented prestressing approach into a standard finite element code is relative simple. There are basically two adjustments to be done compared to a standard nonlinear finite element calculation: (i) throughout the simulation (during prestressing and afterwards) the total deformation gradient (

3) which includes the pre-existing deformation gradient has to be used to compute the constitutive equation Σ(F, F

_{r}) and (ii) during a prestressing increment the pre-existing deformation gradient F

_{r} has to be incrementally updated according to (

12). The implementation can be performed at the material point level, where the components of F

_{r} can be stored as state variables. In contrast, the original MULF method required for an implementation at the finite element level to perform the update of the Jacobian matrix. However, the present approach requires additional increments to calculate the pre-existing deformation gradient and is therefore computationally more expensive than the original approach by

Gee et al. (2010).

2.3 Prestressing Soft Tissues with Distributed Collagen Fibrils

In general, the forward incremental prestressing method outlined in the previous subsection is independent of the constitutive formulation. However, special care must be exercised when dealing with anisotropic constitutive models in which (generalized) structure tensors account for predominant collagen fibril orientations. These models usually define preferred fibril orientations at the unloaded configuration, which is an unknown configuration in our case.

In this work, we use the microstructure-oriented constitutive formulation initially derived in

Grytz (2008). The model is based on the assumption that, at the micro-scale, collagen fibrils crimp into the shape of a helix and are solely subjected to an axial load in fiber direction

**e** (

Grytz and Meschke 2009). At the meso-scale, collagen fibrils are assumed to form a collagen network which can be characterized by a planar distribution of its fibril orientations

**e** (

Grytz and Meschke 2010). The constitutive model contains two material parameters (the elastic modulus of the fibril per macroscopic tissue volume

*E* and the stiffness parameter of the standard Neo-Hookean model

*c*), two microstructural parameters (the crimp angle

*θ*_{0} and the ratio

*R*_{0}/

*r*_{0} between the radii of the helical wave form and of the fibril cross section) and three meso-structural parameters (the collagen fibril dispersion parameter

*κ* [0; 1/2] and two vectors

**M**_{α}(

*α* = 1, 2) that span the plane in which collagen fibrils are dispersed around the mean fibril orientation

**M**_{1}). The meso-structural parameters were used to formulate a generalized structure tensor similar to the work of

Gasser et al. (2006)
where

**M**_{α} were assumed to form a orthonormal basis

**M**_{i} with

**M**_{3} =

**M**_{1} ×

**M**_{2}. The tensor H

_{2D} links the macro- to the micro-scale in the constitutive model, representing the dispersion of the collagen fibril orientations in an integral sense.

In contrast to the concept of the generalized structure tensor,

Roberts et al. (2009) used the MIL method to calculate a fabric tensor

representing the predominant beam orientations of the regional connective tissue architecture in the LC. Note that the generalized structure tensor (

13) can be represented by the fabric tensor (

14) by setting

*H*_{1} = 1 −

*κ*,

*H*_{2} =

*κ* and

*H*_{3} = 0. Accordingly, we propose that the fabric tensor (

14) can be used as an alternative to the generalized structure tensor (

13) to represent the collagen fibril network in our constitutive formulation.

Together with the macroscopic Cauchy-Green tensor C the generalized structure tensor (

13) or the fabric tensor (

14) is used to compute the average fibril stretch of the collagen fibril network

, which is then used to compute the strain energy contribution of the collagen network from the one-dimensional strain energy density function originally derived for one crimped collagen fibril at the micro-scale

*W*_{fib}(

*λ*_{col}). For the complete derivation and a detailed explanation of the constitutive model, please see

Grytz (2008);

Grytz and Meschke (2009,

2010). Note that in contrast to our previous work the collagen fibril dispersion is represented here by one instead of two generalized structure tensors of the form (

13).

Note that tissue deformation changes the micro-and meso-structure of the collagen network including the collagen fibril orientations

**e** (see ). As the constitutive model was formulated with respect to the unloaded configuration, so were the micro- and meso-structural parameters related to the unloaded configuration of the collagen fibril network. Here, however, the macro-structure of the specimen is assumed to be known from experimental observations only at a preloaded reference configuration

while the unloaded configuration

is unknown. While this presents no problem for the scalar parameters in our constitutive formulation, it is more practical to introduce a second orthonormal basis

to define predominant orientations of the collagen network at the pre-loaded reference configuration. However, to construct the generalized structure tensor (

13) the orthonormal basis

**M**_{i} at the unloaded configuration needs to be estimated from

. This cannot be accomplished by a standard pull-back operation using the pre-existing deformation gradient F

_{r}, where the orthonormality of a pulled-back basis is generally not preserved. Therefore we introduce a orthonormal pull-back operation of the vector basis

as proposed by

Grytz and Meschke (2010)The orthonormal pull-back operation (

15) ensures that vector basis

derived by a standard push-forward application satisfies the following relationship with the orthonormal vector basis

at the pre-loaded reference configuration

: (i) the mean fibril orientations will point in the same direction (

) and (ii) the plane in which the fibrils are dispersed will coincide (

). The fulfillment of these requirements can be easily proven as shown here

Note that together with (

15) and (

13) the stress and elasticity tensor can then be calculated in a standard way using the total deformation gradient F (3) as outlined in detail in

Grytz (2008).