# Foundations of viscoelasticity and application to soft tissues mechanics

Michele Righi and Valentina Balbi

###### Abstract

Soft tissues are complex media, they display a wide range of mechanical
properties such as anisotropy and non-linear stress-strain behaviour. They
undergo large deformations and they exhibit a time-dependent mechanical
behaviour, i.e. they are viscoelastic. In this chapter we review the
foundations of the linear viscoelastic theory and the theory of Quasi-Linear
Viscoelasticity (QLV) in view of developing new methods to estimate the
viscoelastic properties of soft tissues through model fitting. To this aim, we
consider the simple torsion of a viscoelastic Mooney-Rivlin material in two
different testing scenarios: step-strain and ramp tests. These tests are
commonly performed to characterise the time-dependent properties of soft
tissues and allow to investigate their stress relaxation behaviour. Moreover,
commercial torsional rheometers measure both the torque and the normal force,
giving access to two sets of data. We show that for a step test, the linear
and the QLV models predict the same relaxation curves for the torque. However,
when the strain history is in the form of a ramp function, the non-linear
terms appearing in the QLV model affect the relaxation curve of the torque
depending on the final strain level and on the rising time of the ramp.
Furthermore, our results show that the relaxation curve of the normal force
predicted by the QLV theory depends on the level of strain both for a step and
a ramp tests. To quantify the effect of the non-linear terms, we evaluate the
maximum and the equilibrium (as $t\\!\rightarrow\\!\infty$) values of the
relaxation curves. Our results provide useful guidelines to accurately fit QLV
models in view of estimating the viscoelastic properties of soft tissues.

##  1 Introduction

Soft tissues, such as the brain, the skin, tendons and ligaments are
viscoelastic materials, their mechanical behaviour is therefore time-
dependent. Two typical experiments that show the time-dependent nature of soft
tissues consists in stress relaxation and creep tests. In a stress relaxation
test the tissue is suddenly stretched and then held in position for a certain
time while the resulting stress is measured. Conversely, in a creep test, the
load is applied to the tissue and the resulting deformation is measured. For
many soft tissues the stress-relaxation curve has a decaying exponential form.
Stress relaxation has been observed in the brain [2, 3], in ligaments and
tendons [5, 6] and in the skin [15]. At the microscale, the physical
mechanisms behind stress relaxation differ from tissue to tissue. In tendons,
for example, it has been observed that crimping and un-crimping of the
hierarchical structures that build up the tissues, i.e. the individual
collagen fibrils, are responsible for the stress relaxation of the tissue [9].
In the skin, the interaction between collagen and elastic fibers plays a
crucial role in determining the time-dependent behaviour of the tissue. When
the tissue is deformed, the cross-links maintain the structure and allow the
elastic fibers to stretch and relax [10].
However, in practice there is no machine that can instantaneously deform a
tissue. A more realistic test is indeed a ramp test, where the tissue is
deformed in a finite time and then held in that position. The duration of the
ramp phase is called rising time $t^{*}$. When the rising time of the ramp is
nearly zero, the ramp test can be well approximated by a step-strain test.
However, if $t^{*}$ is not small (compared to the characteristic time
constants of the material) modelling the ramp-test as a step-test can
introduce errors in the estimation of the viscoelastic parameters.
From the modelling viewpoint, the simplest constitutive theory that can be
used to describe the time-dependent behaviour of soft tissues is the linear
viscoelastic theory, where the stress is related to the strain by a time-
dependent function which in turn depends on the tissue’s viscoelastic
parameters. Linear models are based on three main assumptions:

  1. 1.

the tissue remembers the past deformation history through a fading memory, so
that contributions to recent strain increments are more important than past
contributions. A typical form of the time-dependent parameters that satisfies
this assumption is a decaying exponential form;

  2. 2.

according to the Boltzmann superposition principle, the total stress at the
current time $t$ is given by the sum of all past stress contribution;

  3. 3.

the deformation applied to the tissue is small.

In early times, linear models have been employed to predict the viscoelastic
behaviour of soft tissues. However, soon scientists have realised that these
models do not provide accurate predictions, mainly because in reality soft
tissues undergo large deformations. To overcome this limitation, Fung proposed
what is now called the Quasi Linear Viscoelastic (QLV) theory, which is the
simplest extension of the linear theory to large deformations [8]. QLV models
can capture stress-relaxation and creep, the strain-rate dependent response
and account for large deformations. Moreover, the governing equations of the
viscoelastic problem can be solved analytically for the most common modes of
deformations used in experiments (e.g. tension, compression, equi-biaxial
tension, simple shear and torsion). Therefore, the constitutive parameters can
be directly estimated through fitting of the experimental data, by
implementing a minimisation algorithm. Linear and QLV models have a common
limitation: being based on the linear superposition principle, they cannot
account for the coupling between different time-scales, which is a limitation,
especially for tissues with hierarchical structures. Although more complex
non-linear models that account for this coupling have been proposed, they are
numerically costly when it comes to model fitting and material parameters
estimation [16, 17]. Another class of non-linear models goes under the
umbrella of internal variable or rate-type models which have recently gained
popularity among the biomechanical community [12, 13, 14]. These models are
based on thermodynamics foundations. According to the multiplicative
decomposition, the gradient of the deformation is split into an elastic and a
viscous part. The resulting stress is then split into the sum of an elastic
and a viscous term. The elastic stress is generally written with respect to an
elastic strain energy function. The viscous stress is written with respect to
a number of internal variables, whose evolution laws are dictated by the
second law of thermodynamics and motivated by the linear theory [18]. This
approach has the advantage of allowing an easy implementation of the
constitutive model into finite element codes. However, when it comes to model
fitting, the resulting equations are in implicit forms and need to be solved
numerically, even for simple deformation modes.
Finally, differential-type models formulate the time-dependent constitutive
equation in terms of the derivatives of the right stretch tensor evaluated at
the current time [11]. Despite being computationally easy to implement, these
models do not allow for an explicit form with respect to the relaxation
functions, therefore they are less straightforward to fit with experimental
data. We conclude this brief introductory review by noting that
viscoelasticity is not the only time-dependent property of soft tissues. Rate-
type effects, such as stiffening and softening as a results of cyclic loading
and unloading and ageing are other common effects displayed by biological
tissues [29, 30].
In this chapter, we focus on viscoelasticity with the aim of providing useful
guidelines on model fitting and estimation of the viscoelastic parameters for
linear and QLV models. We consider two main experimental scenarios, the step-
and-hold test and the ramp-and-hold test for the torsion of a cylindrical
tissue. These tests are common experimental protocols used to investigate the
viscoelastic properties of soft tissues, in particular stress-relaxation. In
Section 2 we review the standard linear viscoelastic theory and its
rheological interpretation. In Section 3 we review the QLV theory following
the formulation proposed in [19]. In Section 4, we consider the simple torsion
of a cylindrical sample. This deformation can be performed with commercially
available rheometers which measure both the torque and the normal force
required to twist a cylindrical sample, giving access to two independent sets
of data. Torsion has been successfully used to characterise the elastic
properties of the brain in large deformations [20] suggesting that the tissue
behaves as a Mooney-Rivlin material. The same behaviour was previously
observed in simple shear experiments [21]. In view of applications to brain
mechanics, we therefore solve the equilibrium equations for a viscoelastic
material whose elastic stress obeys a Mooney-Rivlin law and we calculate the
expressions for the torque and the normal force. In Section 5 we compare the
predictions of the QLV model in the scenario of a step-strain test and of a
ramp test. We conclude the chapter by discussing our results and by
summarising the main findings.

##  2 Linear viscoelastic models

Linear viscoelastic constitutive models are formulated by introducing the time
dependency in the material parameters, a sort of fading memory which remembers
the strain history of the material up to the current configuration.
Accordingly, Hooke’s law $\bm{\sigma}=\mathbb{K}:\bm{\varepsilon}$ rewrites as
follows:

$\bm{\sigma}(t)=\int_{-\infty}^{t}\mathbb{K}(t-\tau):\derivative{\bm{\varepsilon}(\tau)}{\tau}\differential{\tau},$
(1)

where $\bm{\sigma}$ is the stress tensor, $\bm{\varepsilon}$ is the
infinitesimal strain tensor, and $\mathbb{K}(t)$ is called the tensorial
relaxation function and is a fourth-order tensor whose entries are the time-
dependent material parameters. The symbol $:$ denotes the double contraction
between a fourth-order tensor $\mathbb{Y}$ and a second-order tensor Z such
that $(\mathbb{Y}\\!:\\!\textbf{Z})_{ab}=Y_{abcd}Z_{cd}$.
Eq. (1) is based on the Boltzmann superposition principle [22, 23].
Accordingly, the total stress at the current time $t$ can be written as sum of
past stress contributions up to the time $t$. In a one-dimensional setting,
for the generic component $\sigma$ we can then write
$\sigma=\sum_{i}\Delta\sigma_{i}$, as sketched in Figure 1.

Figure 1: Boltzmann superposition principle: the strain history $\varepsilon$
is approximated as sum of steps $\Delta\varepsilon_{i}$ and the resulting
total stress response $\sigma$ is the sum of the stress responses
$\Delta\sigma_{i}$ to each step increment.

Each $\Delta\sigma_{i}$ is the stress response to the step increment
$\Delta\varepsilon_{i}=\derivative{\varepsilon}{t}\Delta t_{i}$ and is
governed by the relaxation function $k_{\text{step}}(t-t_{i})$. Therefore,
each stress increment can be written as
$\Delta\sigma_{i}=k_{\text{step}}(t-t_{i})\Delta\varepsilon_{i}$. By assuming
that the strain history is continuous over time, the sum can be converted into
an integral over time. The stress component $\sigma$ at time $t$ is then given
by the following convolution integral:

$\sigma(t)=\int_{-\infty}^{t}k_{\text{step}}(t-\tau)\derivative{\varepsilon(\tau)}{\tau}\differential{\tau},$
(2)

for a given deformation history $\varepsilon(t)$.
Eq. (1) is the tensorial version of Eq. (2). The components of the tensor
$\mathbb{K}(t)$ are the different relaxation functions of the tissue, i.e. the
time-dependent mechanical parameters. We will show in Section 2.1 that the
tensor $\mathbb{K}(t)$ can be split into its components according to a set of
fourth-order bases. Thus, the resulting constitutive equation can be written
with respect to different sets of relaxation functions according to the choice
of bases.
From an experimental viewpoint, the choice of the bases might be dictated by
which mechanical properties we want to estimate. For instance, if we are
interested in estimating the time-dependent shear modulus we will perform a
simple shear test or a torsion test, whereas if we want to estimate the time-
dependent Young’s modulus we will perform a tensile test. Moreover, in order
to be able to estimate the components of $\mathbb{K}(t)$ we first have to
specify their functional form with respect to time. In Section 2.2 we will
review a common form used in the biomechanics community, i.e. the Prony series
form. Furthermore, we will discuss two common experimental protocols that are
performed to estimate the components of $\mathbb{K}(t)$, namely the step-
strain and the ramp tests.

###  2.1 Bases decomposition for the tensor $\mathbb{K}(t)$

In this section, we focus on the tensorial nature of the relaxation function
$\mathbb{K}(t)$ and we show that the constitutive equation (1) can be written
with respect to different sets of components of $\mathbb{K}(t)$, i.e. the
time-dependent mechanical properties of the tissue, according to different
choices of fourth-order tensorial bases. To simplify the analysis, in this
chapter we restrict out attention to homogeneous isotropic tissues. The
mechanical behaviour of such tissues is fully described by two independent
mechanical parameters, e.g. the bulk and the shear modulus, $\kappa(t)$ and
$\mu(t)$ respectively, or the first Lamé parameter $\lambda(t)$ and the shear
modulus. According to the elasticity theory, the elasticity tensor
$\mathbb{C}$ for a homogeneous isotropic material depends only on two
independent elastic constants. Similarly, the tensorial relaxation function
$\mathbb{K}(t)$ for an isotropic material has two independent components
$K_{1}(t)$ and $K_{2}(t)$ with respect to two bases $\mathbb{I}_{1}$ and
$\mathbb{I}_{2}$, respectively, such that:

$\mathbb{K}(t)=\sum_{n=1,2}K_{n}(t)\mathbb{I}_{n},$ (3)

A well-known form of the constitutive equation for a homogeneous isotropic and
compressible material follows by splitting the infinitesimal strain tensor
$\bm{\varepsilon}$ into its hydrostatic and deviatoric parts. The hydrostatic
part is associated to volume changes, whereas the deviatoric part is
associated to the volume-preserving part of the deformation. The following set
of bases splits the strain tensor into its hydrostatic and deviatoric parts:

$\mathbb{I}_{1abcd}=\dfrac{1}{3}\delta_{ab}\delta_{cd}\qquad\text{and}\qquad\mathbb{I}_{2abcd}=\dfrac{1}{2}\left(\delta_{ac}\delta_{bd}+\delta_{ad}\delta_{bc}\right)-\dfrac{1}{3}\delta_{ab}\delta_{cd}.$
(4)

Accordingly, Eq. (1) takes the following form:

$\begin{split}\bm{\sigma}(t)&=\int_{-\infty}^{t}\sum_{n=1,2}K_{n}(t-\tau)\mathbb{I}_{n}:\derivative{\bm{\varepsilon}(\tau)}{\tau}\differential{\tau}\\\
&=\int_{-\infty}^{t}K_{1}(t-\tau)\derivative{\mathbb{I}_{1}:\bm{\varepsilon}(\tau)}{\tau}\differential{\tau}+\int_{-\infty}^{t}K_{2}(t-\tau)\derivative{\mathbb{I}_{2}:\bm{\varepsilon}(\tau)}{\tau}\differential{\tau}\\\
&=\int_{-\infty}^{t}K_{1}(t-\tau)\derivative{\tau}\left(\dfrac{1}{3}\tr\left(\bm{\varepsilon}(\tau)\right)\textbf{I}\right)\differential{\tau}+\int_{-\infty}^{t}K_{2}(t-\tau)\derivative{\tau}\left(\bm{\varepsilon}(\tau)-\dfrac{1}{3}\tr\left(\bm{\varepsilon}(\tau)\right)\textbf{I}\right)\differential{\tau}\\\
&=\int_{-\infty}^{t}\kappa(t-\tau)\derivative{\tau}\left(\tr\left(\bm{\varepsilon}(\tau)\right)\textbf{I}\right)\differential{\tau}+2\int_{-\infty}^{t}\mu(t-\tau)\derivative{\tau}\left(\operatorname{dev}\left(\bm{\varepsilon}(\tau)\right)\right)\differential{\tau},\end{split}$
(5)

where $\delta_{ab}$ is the Kronecker delta ($\delta_{ab}\\!=\\!1$ if
$a\\!=\\!b$ and $\delta_{ab}\\!=\\!0$ if $a\\!\neq\\!b$), I is the second-
order identity tensor and
$\operatorname{dev}\bm{\varepsilon}=\bm{\varepsilon}-\frac{1}{3}\tr\left(\bm{\varepsilon}\right)\textbf{I}$
is the deviatoric part of the second-order tensor $\bm{\varepsilon}$. The
bases $\mathbb{I}_{1}$ and $\mathbb{I}_{2}$ defined in Eq. (5) act on a
second-order tensor by splitting the tensor into its spherical and deviatoric
parts, respectively. The associated material parameters $\kappa(t)$ and
$\mu(t)$ are the time-dependent bulk and shear modulus, respectively.
For incompressible materials, i.e. materials that deform by keeping their
volume constant, the bulk modulus is much greater than the shear modulus
($\kappa(t)\gg\mu(t)$ for $\forall t$). Moreover, the following assumptions
are true: $\tr\bm{\varepsilon}(t)\rightarrow 0$,
$\kappa(t)\rightarrow\infty,\forall t$. In these limits, Eq. (5) reduces to:

$\bm{\sigma}=-p(t)\textbf{I}+2\int_{-\infty}^{t}\mu(t-\tau)\derivative{\tau}\left(\operatorname{dev}\bm{\varepsilon}(\tau)\right)\differential{\tau},$
(6)

where we have introduced the Lagrange multiplier $p(t)$:

$-p(t)=\lim_{\tr\bm{\varepsilon}(t)\rightarrow
0}\lim_{\kappa(t)\rightarrow\infty}\int_{-\infty}^{t}\kappa(t-\tau)\derivative{\tau}\left(\tr\bm{\varepsilon}(\tau)\right)\differential{\tau}.$
(7)

The scalar $p(t)$ can be interpreted as a hydrostatic pressure and can be
calculated by solving the governing equations of motions for a continuum body,
upon imposing the boundary conditions. Note that the stress component
$-p(t)\textbf{I}$ represents a workless reaction with respect to the kinematic
constraint of the deformation field. No dissipation is involved in the
isochoric deformation of the body. Hence, for materials that can be treated as
incompressible, only the deviatoric part of the stress exhibits a viscoelastic
nature.

Similarly, by choosing the following bases

$\mathbb{J}_{1abcd}=\delta_{ab}\delta_{cd}\qquad\text{and}\qquad\mathbb{J}_{2abcd}=\dfrac{1}{2}\left(\delta_{ac}\delta_{bd}+\delta_{ad}\delta_{bc}\right),$
(8)

$\mathbb{K}(t)=\sum_{n=1,2}A_{n}(t)\mathbb{J}_{n}$ and Eq. (1) writes as
follows:

$\begin{split}\bm{\sigma}(t)&=\int_{-\infty}^{t}\sum_{n=1,2}A_{n}(t-\tau)\mathbb{J}_{n}:\derivative{\bm{\varepsilon}(\tau)}{\tau}\differential{\tau}=\cdots=\\\
&=\int_{-\infty}^{t}A_{1}(t-\tau)\derivative{\tau}\left(\tr\bm{\varepsilon}(\tau)\textbf{I}\right)\differential{\tau}+\int_{-\infty}^{t}A_{2}(t-\tau)\derivative{\tau}\left(\bm{\varepsilon}(\tau)\right)\differential{\tau}\\\
&=\int_{-\infty}^{t}\lambda(t-\tau)\derivative{\tau}\left(\tr\bm{\varepsilon}(\tau)\textbf{I}\right)\differential{\tau}+2\int_{-\infty}^{t}\mu(t-\tau)\derivative{\tau}\left(\bm{\varepsilon}(\tau)\right)\differential{\tau}.\end{split}$
(9)

Now, we have $A_{1}(t)=\lambda(t)$, which is the time-dependent first Lamé
parameter and $A_{2}(t)=2\mu(t)$. Moreover, note that
$\mathbb{J}_{1}=3\mathbb{I}_{1}$ and
$\mathbb{J}_{2}=\mathbb{I}_{2}-\mathbb{I}_{1}$ and the following link is true
$\lambda(t)=\kappa(t)-\frac{2}{3}\mu(t)$, $\forall t$. Clearly, Eqs. (5) and
(9) are equivalent forms of the constitutive equation (1) and predict the same
stress response to a general strain input $\bm{\varepsilon}(t)$. The choice of
the bases and therefore the final form of the constitutive model is usually
dictated by what type of material properties we want to determine and the type
of experimental devices available for testing (e.g. tensile machines,
rheometers, bi-axial devices, etc.). In Section 3 we will use the bases
decomposition of Eq. (4) to write the constitutive equation for the QLV model.
In the next section, we focus on the mathematical form of the components of
the tensor $\mathbb{K}(t)$.

###  2.2 Rheological models for the relaxation function

In order to use the constitutive equation (1) for model fitting and parameter
estimations, a mathematical form for the components of the tensorial function
$\mathbb{K}(t)$ has to be chosen. On the one hand, the form of the relaxation
functions is restricted by the following physical principles: positive strain
energy and satisfaction of the second law of thermodynamics. Imposing the
energy density to be non-negative during the tissue relaxation results in
requiring the relaxation function to be positive $\forall\,t$ [7].
Furthermore, the second law of thermodynamics (dissipation inequality)
requires that the relaxation function decreases monotonically with time [7].
On the other hand, any function that satisfies the physical constraints and
replicates the shape of the observed stress relaxation curve can be used. From
the experimental viewpoint, a classical experiment that can be done to
determine the relaxation curve is to apply a displacement to the tissue, then
hold the tissue in position (i.e. maintain a constant level of strain) for a
certain time and measure the resulting stress curve. For most soft tissues,
the measured stress relaxation curve has a decaying exponential behaviour [4,
5, 6].

The simplest form for the relaxation function that captures the exponential
decaying behaviour and satisfies the physical constraints is the so-called
Prony series. The Prony-series has its origin in one-dimensional rheological
models [7, 8]. Such models are represented by an arrangement of linear springs
and linear dash-pots. The layout of such arrangements of elements provides the
qualitative behaviour of the system, e.g. solid-like or fluid-like behaviour,
while the values of the constants characterise the quantitative behaviour. In
Figure 2 we sketch the generalized Maxwell scheme, which is used to model the
viscoelastic response of solid materials. The isolated spring $k_{\infty}$
represents the residual (long-term) elasticity of the tissue.

Figure 2: Rheological models: springs and dash-pots arrangement of the
generalized Maxwell model. The parameters $k_{i}$ are spring constants
($\mathrm{P}\mathrm{a}$) and $\eta_{i}$ are dash-pot constants
($\mathrm{Pa}\text{\,}\mathrm{s}$), or viscosities. $\sigma$ is the stress
response of the system to the applied strain $\varepsilon$.

One-dimensional rheological models can be described by a linear ordinary
differential equation in the variables $\sigma$ (the stress) and $\varepsilon$
(the strain). The stress response $\sigma$ of the system to a step-strain
input $\varepsilon$ provides the form of the relaxation function
$k_{\text{step}}$ in (2). A convenient way to derive the response of
rheological models, especially when a large number of elements is involved, is
to employ the Laplace transform. The approach involves the following steps:

  1. 1.

Write the constitutive equations for all the elements in the system:
$\sigma_{i}=k_{i}\varepsilon_{i}$ for springs and
$\sigma_{i}=\eta_{i}\dot{\varepsilon}_{i}$ for dash-pots. Then write the
equilibrium equations $\sigma=\sum_{i}\sigma_{i}$ for elements in parallel
that experience the same strain and $\varepsilon=\sum_{i}\varepsilon_{i}$ for
elements in series that experience the same stress.

  2. 2.

By applying the Laplace transform to each equations at point 1, convert the
system of mixed differential and algebraic equations into a system of only
algebraic equations.

  3. 3.

Apply a variable elimination procedure to the system derived at point 2. The
system has $3n+2$ equations in $3n+3$ unknowns, i.e. $\varepsilon_{i}$ and
$\sigma_{i}$. This reduces the system to the single equation:
$\bar{\sigma}(s)=F(s)\bar{\varepsilon}(s)$, where $\bar{\sigma}$ and
$\bar{\varepsilon}$ are the Laplace transforms of $\sigma(t)$ and
$\varepsilon(t)$, respectively. $F(s)$ is the transfer function in the complex
domain represented by the complex variable $s$.

  4. 4.

By applying the inverse Laplace transform to the equation
$\bar{\sigma}(s)=F(s)\bar{\varepsilon}(s)$, obtain the constitutive equation
in the time domain.

From point 3, the transfer function $F(s)$ of the generalized Maxwell model in
Figure 2 is given by:

$F(s)=\frac{\bar{\sigma}(s)}{\bar{\varepsilon}(s)}=k_{\infty}+\sum_{i=1}^{n}\frac{s\eta_{i}}{1+s\tau_{i}},$
(10)

where the constants $\tau_{i}=\frac{\eta_{i}}{k_{i}}$ are called relaxation
times of the model.

To calculate the response of a generalised Maxwell system to a step input, we
write the strain $\varepsilon(t)=\varepsilon_{0}H(t)$, where $H(t)$ is the
Heaviside function: $H(t)=1$, $\forall t\geq 0$ and $H(t)=0$, $\forall t<0$
and $\varepsilon_{0}$ is the amplitude of the step. Then we calculate the
Laplace transform $\bar{\varepsilon}$ and substitute the result into Eq. (10)
obtaining:

$\bar{\sigma}(s)=\left(\frac{k_{\infty}}{s}+\sum_{i=1}^{n}\frac{\eta_{i}}{1+s\tau_{i}}\right)\varepsilon_{0}.$
(11)

According to point 4, by transforming back into the time domain, we obtain the
stress response to a step-strain input with amplitude $\varepsilon_{0}$:

$\sigma(t)=\left(k_{\infty}+\sum_{i=1}^{n}k_{i}e^{-\frac{t}{\tau_{i}}}\right)\varepsilon_{0}.$
(12)

From Eq. (12) we can then calculate the stress response to a step-strain
input, i.e. the relaxation function $k_{\text{step}}(t)$, by dividing
$\sigma(t)$ by the amplitude of the step, as follows:

$k_{\text{step}}(t)=\frac{\sigma(t)}{\varepsilon_{0}}=k_{\infty}+\sum_{i=1}^{n}k_{i}e^{-\frac{t}{\tau_{i}}}.$
(13)

The function (13) is called Prony series and is a sum of exponential terms,
each corresponding to a branch of the generalised Maxwell in Figure 2. Note
that the relaxation function $k_{\text{step}}(t)$ does not depend on the
strain $\varepsilon$.

The constants $k_{\infty},k_{i}$ and $\eta_{i}$ can be determined by fitting
Eq. (13) to the stress relaxation curve experimentally measured from a step-
strain test. The sample is suddenly deformed up to the strain
$\varepsilon_{0}$ and held in position for a certain amount of time. In the
limit $t\rightarrow\infty$ Eq. (12) recovers the elastic equilibrium stress
($\sigma_{\infty}\\!=\\!k_{\infty}\varepsilon_{0}$) and the relaxation
function in Eq. (13) reduces to the long term elastic modulus $k_{\infty}$:

$k_{\infty}=\lim_{t\rightarrow\infty}k_{\text{step}}(t).$ (14)

Experimentally, this limit is equivalent to a very slow ramp test, i.e. a
quasi-static test, where the final value of strain $\varepsilon_{0}$ is
attained as $t\\!\rightarrow\\!\infty$.
On the other hand, in the limit $t\rightarrow 0$, Eq. (13) reduces to the
instantaneous elastic modulus $k_{0}$:

$k_{0}=k_{\text{step}}(0)=k_{\infty}+\sum_{i=1}^{n}k_{i},$ (15)

The value in Eq. (15) corresponds to the maximum of the relaxation function,
see Figure 3. Experimentally, this limit corresponds to the application of an
instantaneous strain, which is practically impossible to perform. The
constants $k_{0}$ and $k_{\infty}$ are the elastic parameters of the
constitutive model and describe the instantaneous and long-term elastic
behaviours of the tissue, respectively.
The viscous behaviour of the tissue is associated to the parameters $\eta_{i}$
and $\tau_{i}$. To get some insights on the viscous parameters, we define the
function:

$\tilde{k}_{\text{step}}(t)=k_{\text{step}}(t)-k_{\infty}=\sum_{i=1}^{n}k_{i}e^{-\frac{t}{\tau_{i}}}$
(16)

and we integrate over the whole time spectrum:

$\eta_{0}=\int_{0}^{\infty}\tilde{k}_{\text{step}}(t)\differential{t}=\sum_{i=1}^{n}\eta_{i}.$
(17)

The value $\eta_{0}$ has an important geometrical interpretation since it
represents the area between the curve $k_{\text{step}}(t)$ and the asymptotic
line $k_{\infty}$, see Figure 3. The bigger the area, the more viscous the
material.
Moreover, we can compute the mean relaxation time $T_{c}$ of the tissue as
follows:

$T_{c}=\frac{\int_{0}^{\infty}t\tilde{k}_{\text{step}}(t)\differential{t}}{\int_{0}^{\infty}\tilde{k}_{\text{step}}(t)\differential{t}}=\frac{\sum_{i=1}^{n}\eta_{i}\tau_{i}}{\sum_{i=1}^{n}\eta_{i}}.$
(18)

Geometrically, $T_{c}$ represents the centroid of the shaded area below the
relaxation function $k_{\text{step}}(t)$ and the asymptotic value
$k_{\infty}$, and can be interpreted as the average relaxation time.

Figure 3: Stress response of the generalised Maxwell model to a step-strain
input.

#### Ramp tests

Now, we recall that Eqs. (13-18) are valid for a strain input in the form of a
step function (i.e. $\varepsilon(t)=H(t)\varepsilon_{0}$), where the strain
value $\varepsilon_{0}$ is attained instantaneously. However, such experiment
is not feasible in laboratory since that would require a testing machine able
to reach an infinite rate of deformation. Real tests are much closer to a ramp
test, where the strain $\varepsilon_{0}$ is reached after a finite rising time
$t^{*}>0$, as shown in Figure 4.

Figure 4: Strain history for a ramp and hold experiment. The constant strain
value $\varepsilon_{0}$ is reached at the end of the loading phase
($t=t^{*}$), where the strain increases at a constant rate.

In view of providing an analytical expression for the relaxation curve to be
fitted with the experimental data, in this section we derive the stress
response of a generalised Maxwell system to a ramp input. The strain input for
a ramp test takes the following form:

$\varepsilon(t)=\frac{\varepsilon_{0}}{t^{*}}t-\frac{\varepsilon_{0}}{t^{*}}(t-t^{*})H(t-t^{*}),$
(19)

where $H(t)$ is the Heaviside function. We can calculate the stress response
to the input in Eq. (19) by substituting Eq. (13) and Eq. (19) into Eq. (2).
We then obtain:

$\sigma(t)=k_{\text{ramp}}(t)\varepsilon_{0},\qquad\text{for}\quad t>t^{*},$
(20)

where:

$k_{\text{ramp}}(t)=\left(k_{\infty}+\frac{1}{t^{*}}\sum_{i=1}^{n}\tau_{i}k_{i}e^{-\frac{t}{\tau_{i}}}\left(e^{\frac{t^{*}}{\tau_{i}}}-1\right)\right),\qquad\text{for}\quad
t>t^{*}.$ (21)

Note that by taking the limit for $t\rightarrow\infty$ in Eq. (21), we recover
the long term elastic modulus $k_{\infty}$. Similarly, the equilibrium stress
as $t\rightarrow\infty$ is given by
$\sigma_{\infty}\\!=\\!k_{\infty}\varepsilon_{0}$ as for the step test.
On the other hand, the stress response at $t\\!=\\!t^{*}$ is now affected by
the previous deformation history at $t\\!<\\!t^{*}$. To quantify the effect of
the rising time of the ramp on the instantaneous elastic and viscous response,
we calculate the function $\tilde{k}_{\text{ramp}}(t)$:

$\tilde{k}_{\text{ramp}}(t)=k_{\text{ramp}}(t)-k_{\infty}=\sum_{i=1}^{n}\frac{k_{i}}{\nu_{i}}e^{-\frac{t}{\tau_{i}}}\left(e^{\nu_{i}}-1\right),$
(22)

where $\nu_{i}=t^{*}/\tau_{i}$ are the ratio between the rise time and the
characteristic time constants of the tissue $\tau_{i}$. Note that by taking
the limit $t^{*}\rightarrow 0$ (i.e. a step input), Eq. (22) recovers Eq.
(16):

$\lim_{t^{*}\rightarrow
0}\tilde{k}_{\text{ramp}}(t)=\tilde{k}_{\text{step}}(t)=\sum_{i=1}^{n}k_{i}e^{-\frac{t}{\tau_{i}}}.$
(23)

Moreover, we define the modified instantaneous elastic modulus
$k_{0\text{ramp}}$ as the function (22) evaluated at the end of the ramp
phase, i.e. at $t=t^{*}$:

$k_{0\text{ramp}}=\tilde{k}_{\text{ramp}}(t^{*})=\sum_{i=1}^{n}k_{i}\zeta_{i}=\bm{k}\cdot\bm{\zeta},\quad\text{with}\quad\zeta_{i}=\nu_{i}^{-1}\left(1-e^{-\nu_{i}}\right).$
(24)

$\bm{k}$ and $\bm{\zeta}$ are vectors with components $k_{i}$ and $\zeta_{i}$,
$i\\!=\\!\\{1,\dots,n\\}$, respectively. Since $0\leq\nu_{i}<\infty$, then the
parameters $\zeta_{i}$ range between $0<\zeta_{i}\leq 1$. When $\nu_{i}=0$,
$\zeta_{i}=1$ and the ramp recovers the perfect step-strain input and the
elastic modulus in Eq. (24) reduces to the instantaneous modulus in Eq. (15).

We note that the elastic constants $k_{i}$ are intrinsic properties of the
tissue, therefore they do not depend on the testing procedure nor on the form
of the strain history or on the strain-rate at which the test is performed. On
the other hand it is well-known that the response of a viscoelastic material
strongly depends on the strain-rate. The coefficients $\zeta_{i}$ account for
the strain rate of the deformation process, i.e. for the fact that the strain
is applied to the tissue in a finite time $t^{*}$. The vector $\bm{\zeta}$
allows to isolate the effect of the strain-rate from the constant elastic
moduli $\bm{k}$ that describe the material. Therefore, fitting stress
relaxation data by using the relaxation function of a step instead of the ramp
will result in underestimating the elastic moduli. The smaller the strain-rate
(i.e. the greater $t^{*}$), the smaller $k_{0\text{ramp}}$. In the limit of a
quasi-static deformation, only the infinite modulus $k_{\infty}$ is
recoverable.

Similarly to Eq. (17) we can now compute the total viscosity of the material
using the modified function $\tilde{k}_{\text{ramp}}(t)$ in Eq. (22):

$\eta_{0\text{ramp}}=\int_{t^{*}}^{\infty}\tilde{k}_{\text{ramp}}(t)\differential{t}=\sum_{i=1}^{n}\eta_{i}\zeta_{i}=\bm{\eta}\cdot\bm{\zeta}.$
(25)

Eq. (25) highlights that the same set of coefficients $\zeta_{i}$ that link
the elastic moduli also link the viscous constants. Since for cases of
practical interest $0<\zeta_{i}<1$, it follows that
$\eta_{0\text{ramp}}<\eta_{0}$. Therefore, neglecting the influence of the
deformation rate will result in underestimating the total viscosity of the
material. This is in agreement with the result
$\tilde{k}_{0\text{ramp}}<\tilde{k}_{0}$, with
$\tilde{k}_{0}=k_{0}-k_{\infty}$.
In conclusion, by fitting a stress relaxation curve obtained from a ramp test
with Eq. (13) we can obtain a correct estimation of the infinite modulus
$k_{\infty}$, which is strain-rate independent. However, the peak of the
relaxation curve, which is related to the instantaneous response of the tissue
and therefore to its instantaneous elastic modulus $k_{0}$, depends on the
strain rate and it is given by Eq. (24). In particular, the lower the strain
rate, the lower the peak. The area enclosed between the relaxation function
and its horizontal asymptote is related to the total viscosity of the
material. Since $k_{\infty}$ is not affected by the rate of deformation,
$k_{0\text{ramp}}$ decreases with the area represented by
$\eta_{0\text{ramp}}$.
These preliminary synthetic information ($k_{\infty}$, $k_{0\text{ramp}}$,
$\eta_{0\text{ramp}}$) derived from observation of the experimental relaxation
function can be used as a starting point of the fitting procedure to determine
the constitutive parameters of the rheological model ($k_{\infty}$, $k_{i}$,
$\tau_{i}$).

##  3 QLV model

The linear model in Eq. (1) predicts accurate results only in the small
deformation regime, i.e. when $\varepsilon\approx 0$. However it fails to
accurately predict the stress response when a tissue is subjected to a large
deformation. To account for large deformations, Fung originally proposed the
theory of Quasi Linear Viscoelasticity (QLV) [8], which is the extension of
the linear theory we reviewed in the previous section to the large deformation
regime. In this section we review the QLV theory, following [19] to derive the
constitutive equation for isotropic compressible and incompressible soft
tissues.
The QLV theory is based on the same assumptions of the linear theory, i.e. the
Boltzmann superposition principle and the assumption of fading memory.
Moreover, Fung postulated that the total stress is separable into the product
of a function of time, i.e. the relaxation function, and a function of the
deformation, i.e. the elastic stress. The former accounts for the time-
decaying relaxation of the stress and the latter accounts for the non-linear
elastic response of the tissue. In the QLV formulation, the relation between
the elastic stress and the strain is non-linear. To write the constitutive
equation for a QLV model, we start by rewriting the linear model in Eq. (1) in
the following equivalent form:

$\bm{\sigma}(t)=\int_{0}^{t}\mathbb{G}(t-\tau):\dfrac{\differential\bm{\sigma}^{\text{e}}(\tau)}{\differential\tau}\differential{\tau}$
(26)

where the tensor $\mathbb{G}(t)$ is now a fourth-order tensor whose components
are non-dimensional and such that $G_{n}(0)=1$ for $n=\\{1,2\\}$. We call
$\mathbb{G}(t)$ the reduced relaxation tensor. Note that in Eq. (26) we have
assumed that the deformation history starts at $t=0$. The stress term
$\bm{\sigma}^{\text{e}}=\mathbb{K}(0)\bm{\varepsilon}$ is the linear elastic
stress. Fung proposed to replace the linear stress $\bm{\sigma}^{\text{e}}$ by
the corresponding instantaneous elastic stress in large deformation and
rewrite the constitutive equation (26) as follows:

$\bm{\Pi}(t)=\int_{0}^{t}\mathbb{G}(t-\tau):\dfrac{\differential\bm{\Pi}^{\text{e}}(\tau)}{\differential\tau}\differential{\tau}.$
(27)

The tensor $\bm{\Pi}(t)$ is the second Piola-Kirchhoff stress tensor and
$\bm{\Pi}^{\text{e}}(t)$ is the elastic second Piola-Kirchhoff stress tensor
defined as follows:

$\bm{\Pi}^{\text{e}}=J\textbf{F}^{-1}\textbf{T}^{\text{e}}\textbf{F}^{-\text{T}}.$
(28)

$\textbf{T}^{\text{e}}$ is the elastic Cauchy stress,
$\textbf{F}=\textbf{I}-\nabla\textbf{u}=\partial\textbf{x}/\partial\textbf{X}$
is the deformation gradient associated to the large deformation
$\textbf{x}=\chi(\textbf{X})$, and $J=\det\textbf{F}$. x and X are the
position vectors in the undeformed and deformed configurations, respectively.
We use the notation T to avoid confusion with the linear stress tensor
$\bm{\sigma}$. In the small deformation regime the undeformed and deformed
configurations coincide since $\nabla\textbf{u}\approx 0$ and
$\textbf{F}\approx\textbf{I}$ and therefore the second Piola-Kirchhoff and the
Cauchy stress tensors also reduce to the same stress tensor.
Now, we can use the bases in Eq. (9) to split the tensor $\mathbb{G}$ of Eq.
(27). We call the associated components $\mathcal{H}(t)$ and $\mathcal{D}(t)$,
respectively, and we rewrite Eq. (27) in the following form:

$\begin{split}\bm{\Pi}(t)&=\int_{0}^{t}\mathcal{H}(t-\tau)\dfrac{\differential\bm{\Pi}^{\text{e}}_{\text{H}}(\tau)}{\differential\tau}\differential{\tau}+\int_{0}^{t}\mathcal{D}(t-\tau)\dfrac{\differential\bm{\Pi}^{\text{e}}_{\text{D}}(\tau)}{\differential\tau}\differential{\tau}.\end{split}$
(29)

By following [19], we define:

$\bm{\Pi}^{\text{e}}_{\text{H}}=J\textbf{F}^{-1}\left(\dfrac{1}{3}\text{tr}(\textbf{T}^{\text{e}})\textbf{I}\right)\textbf{F}^{-\text{T}}\qquad\text{and}\qquad\bm{\Pi}^{\text{e}}_{\text{D}}=J\textbf{F}^{-1}\left(\text{dev}(\textbf{T}^{\text{e}})\right)\textbf{F}^{-\text{T}},$
(30)

so that
$\bm{\Pi}^{\text{e}}=\bm{\Pi}^{\text{e}}_{\text{H}}+\bm{\Pi}^{\text{e}}_{\text{D}}$.
Note that the relaxation functions $\mathcal{H}(t)$ and $\mathcal{D}(t)$ are
associated to the Piola transformations of the hydrostatic and deviatoric
parts of the Cauchy stress, $\bm{\Pi}^{\text{e}}_{\text{H}}$ and
$\bm{\Pi}^{\text{e}}_{\text{D}}$, respectively. Moreover, by comparing Eqs.
(26) and (27) we see that they are both written with respect to the same
tensorial relaxation function $\mathbb{G}(t)$. Therefore, the components
$\mathcal{H}(t)$ and $\mathcal{D}(t)$ can be determined by performing step-
strain tests in the linear regime. Upon a closer inspection of Eq. (5) we can
also note that $\mathcal{H}(t)=\kappa(t)/\kappa_{0}$ and
$\mathcal{D}(t)=\mu(t)/\mu_{0}$ are the non-dimensional version of the
relaxation functions $\kappa(t)$ and $\mu(t)$, where $\mu_{0}=\mu(0)$ and
$\kappa_{0}=\kappa(0)$ are the instantaneous elastic bulk and shear modulus,
respectively.
Finally, the Cauchy stress tensor follows from applying the transformation
$\textbf{T}=\\!J^{-1}\textbf{F}\bm{\Pi}\textbf{F}^{\text{T}}$ to Eq. (29) and
is given by:

$\textbf{T}(t)=J^{-1}(t)\textbf{F}(t)\left(\int_{0}^{t}\mathcal{H}(t-\tau)\dfrac{\differential\bm{\Pi}^{\text{e}}_{\text{H}}(\tau)}{\differential\tau}\differential{\tau}+\int_{0}^{t}\mathcal{D}(t-\tau)\dfrac{\differential\bm{\Pi}^{\text{e}}_{\text{D}}(\tau)}{\differential\tau}\differential{\tau}\right)\textbf{F}^{\text{T}}(t).$
(31)

Eq. (29) is the QLV form of the constitutive equation for an isotropic
compressible viscoelastic material. In the incompressible limit
$J\\!\rightarrow\\!1$ and
$\kappa(t)\\!\rightarrow\\!\kappa_{0}\\!\rightarrow\\!\infty$, $\forall t$,
therefore Eq. (31) reduces to the following form:

$\textbf{T}(t)=\textbf{F}(t)\left(\int_{0}^{t}\mathcal{D}(t-\tau)\dfrac{\differential\bm{\Pi}^{\text{e}}_{\text{D}}(\tau)}{\differential\tau}\differential{\tau}\right)\textbf{F}^{\text{T}}(t)-p(t)\textbf{I}$
(32)

where the Lagrange multiplier $p(t)$ is given by:

$p(t)=\lim_{k(t)\rightarrow\infty}\lim_{J\rightarrow
1}\left(J^{-1}(t)\textbf{F}(t)\left(\int_{0}^{t}\dfrac{\kappa(t-\tau)}{\kappa_{0}}\dfrac{\differential\bm{\Pi}^{\text{e}}_{\text{H}}(\tau)}{\differential\tau}\differential{\tau}\right)\textbf{F}^{\text{T}}(t)\right).$
(33)

Eq. (32) is the QLV form of the constitutive equation for an isotropic
incompressible material.
In the next section we consider the simple torsion of a solid cylinder. We
derive the analytical expressions of the torque and the normal force required
to twist the cylinder, both in the linear and the large deformations regime.
We then derive the analytical expression for the relaxation curves of the
torque and the normal force in two experimental scenarios: the step-strain
test and the ramp test.

##  4 Simple Torsion

In this section we consider the problem of simple torsion of a solid cylinder.
We start by defining the coordinates of the cylinder in the reference
configuration $\mathcal{B}_{0}$ and in the deformed configuration
$\mathcal{B}(t)$ as $\\{R,\Theta,Z\\}$ and $\\{r(t),\theta(t),z(t)\\}$,
respectively. We assume that the deformation starts at time $t\\!=\\!0$ and
take the reference configuration as the initial configuration
$\mathcal{B}_{0}\\!=\\!\mathcal{B}(0)$. The displacement vectors X and
$\textbf{x}(t)$ in $\mathcal{B}_{0}$ and $\mathcal{B}(t)$, respectively, are
defined with respect to the bases
$\\{\textbf{E}_{R},\textbf{E}_{\Theta},\textbf{E}_{Z}\\}$ and
$\\{\textbf{e}_{r},\textbf{e}_{\theta},\textbf{e}_{z}\\}$, so that
$\textbf{X}=R\textbf{E}_{R}+\Theta\textbf{E}_{\Theta}+Z\textbf{E}_{Z}$ and
$\textbf{x}(t)=r(t)\textbf{e}_{r}+\theta(t)\textbf{e}_{\theta}+z(t)\textbf{e}_{z}$.
The deformation can then be written as follows:

$r(t)=R,\qquad\theta(t)=\Theta+\phi(t)Z,\qquad z(t)=Z,$ (34)

where $\phi(t)=\alpha(t)/l$ is the amount of twist experienced by the cylinder
at time $t$, defined as the angle of rotation $\alpha(t)$ per unit length. $l$
is the length of the cylinder which remains constant at all times. The strain
$\gamma(r,t)$, a non-dimensional measure of the deformation is:

$\gamma(r,t)=\dfrac{r\alpha(t)}{l}=r\phi(t).$ (35)

The deformation gradient
$\textbf{F}(r,t)=\dfrac{\partial\textbf{x}(t)}{\partial\textbf{X}}$ is given
by:

$\textbf{F}(r,t)=\left(\begin{array}[]{ccc}1&0&0\\\ 0&1&r\,\phi(t)\\\ 0&0&1\\\
\end{array}\right)$ (36)

and the left Cauchy-Green tensor
$\textbf{B}(r,t)=\textbf{F}(r,t)\textbf{F}(r,t)^{\text{T}}$ and its inverse
are given by:

$\textbf{B}(r,t)=\left(\begin{array}[]{ccc}1&0&0\\\
0&1+r^{2}\phi^{2}(t)&r\phi(t)\\\ 0&r\phi(t)&1\\\
\end{array}\right)\qquad\text{and}\qquad\textbf{B}(r,t)^{-1}=\left(\begin{array}[]{ccc}1&0&0\\\
0&1&-r\phi(t)\\\ 0&-r\phi(t)&1+r^{2}\phi^{2}(t)\\\ \end{array}\right).$ (37)

Note that the deformation gradient depends on the spatial variable $r$, i.e.
the deformation is non-homogeneous and the stress distribution will depend on
the radial position as well.
The principal stretches and the principal directions associated to the torsion
deformation are the eigenvalues and the eigenvectors of the tensor B,
respectively. Upon diagonalising B, we find that the principal stretches are
given by:

$\lambda_{1}=1,\quad\lambda_{2,3}(r,t)=\sqrt{1+\frac{\gamma(r,t)}{2}\left(\gamma(r,t)\pm\sqrt{\gamma^{2}(r,t)+4}\right)}.$
(38)

$\lambda_{2}$ and $\lambda_{3}$ are the greatest and the smallest stretch
respectively and the associated eigenvectors are the directions where
$\lambda_{2}$ and $\lambda_{3}$ occur. Note that both $\lambda_{2}$ and
$\lambda_{3}$ depend on the spatial variable $r$. Moreover, $\lambda_{2}$ is
maximum at the outer surface $r\\!=\\!r_{o}$. It is useful to define the
strain $\gamma_{o}(t)$ as the strain at the outer surface of the cylinder at
time $t$:

$\gamma_{o}(t)=\gamma(r_{o},t)=r_{o}\phi(t)$ (39)

so that:

$\max_{r}\lambda_{2}(r,t)=\lambda_{2}(r_{o},t)=\sqrt{1+\frac{\gamma_{o}(t)}{2}\left(\gamma_{o}(t)+\sqrt{\gamma_{o}^{2}(t)+4}\right)}$
(40)

(a) a (b) b

Figure 5: Sketch of the lateral surface of a cylinder with length $l$ and
radius $r_{o}$ twisted with $\gamma_{o}=0.02$ (a) and $\gamma_{o}=0.8$ (b).
The green arrow shows the principal direction associated to the maximum
stretch $\lambda_{2}$ in Eq. (38). The dashed and the solid lines represent
the undeformed and the deformed cylinder respectively.

In Figure 5 we show the principal direction (green arrow) associated to the
maximum stretch (i.e. $\lambda_{2}(r_{o},t)$) on the external surface of the
cylinder, when the cylinder experiences a strain $\gamma_{o}(t)=0.02$ (Figure
5a) and $\gamma_{o}(t)=0.8$ (Figure 5b). When $\gamma_{o}(t)\ll 1$, i.e. in
the small deformation regime, the principal direction is aligned with the
diagonal of the rectangle and the maximum stretch $\lambda_{2}(r_{o},t)$ can
be approximated by the following expansion:

$\lambda_{2}(r_{o},t)=1+\dfrac{\gamma_{o}(t)}{2}+\mathcal{O}(\gamma_{o}^{2}(t)),$
(41)

which recovers the relation between stretch $\lambda$ used for large
deformations and the infinitesimal strain $\varepsilon$ used in small
deformations $\lambda=\varepsilon+1$, with $\varepsilon=\gamma/2$.
However, Figure 5b shows that in the large deformation regime, the principal
direction associated to the maximum stretch $\lambda_{2}(r_{o},t)$ is not
aligned with the diagonal and $\lambda_{2}(r_{o},t)$ is given by Eq. (40).

We can now write the governing equations for the simple torsion of a solid
cylinder. Assuming that the inertia is negligible, the equilibrium equations
at any time $t>0$ are given by:

$\left\\{\begin{split}\text{div}\textbf{T}(r,t)&=\bm{0}\\\
T_{rr}(r,t)&=0\qquad\text{at}\qquad r=r_{o}\end{split}\right.$ (42)

where $r_{o}$ is the outer radius of the cylinder and the last equation
follows from imposing that the lateral surface of the cylinder is free of
traction at any time $t$. The operator div is the divergence operator (in
cylindrical coordinates), see [25] for details. In writing Eq. (42) we have
implicitly neglected inertial forces. For a discussion on the validity of this
assumption we refer to [26, 27].

The Cauchy stress $\textbf{T}(r,t)$ is given by the constitutive equation
(32). Given the form of the tensors $\textbf{F}(r,t)$ and $\textbf{B}(r,t)$
from Eqs. (36) and (37), the Cauchy stress tensor will have components
$T_{r\theta}(r,t)\\!=\\!T_{\theta
r}(r,t)\\!=\\!T_{rz}(r,t)\\!=\\!T_{zr}(r,t)\\!=\\!0$, $\forall t$. The
remaining non-zero components are:

$\begin{split}T_{rr}(r,t)&=\int_{0}^{t}\mathcal{D}(t-\tau)\dfrac{\partial}{\partial\tau}\Pi^{\text{e}}_{\text{D}11}(r,\tau)\differential{\tau}-p(r,\theta,z,t)\\\
T_{\theta\theta}(r,t)&=\int_{0}^{t}\mathcal{D}(t-\tau)\dfrac{\partial}{\partial\tau}\Pi^{\text{e}}_{\text{D}22}(r,\tau)\differential{\tau}+2r\phi(t)\int_{0}^{t}\mathcal{D}(t-\tau)\dfrac{\partial}{\partial\tau}\Pi^{\text{e}}_{\text{D}23}(r,\tau)\differential{\tau}\\\
&+r^{2}\phi(t)^{2}\int_{0}^{t}\mathcal{D}(t-\tau)\dfrac{\partial}{\partial\tau}\Pi^{\text{e}}_{\text{D}33}(r,\tau)\differential{\tau}-p(r,\theta,z,t)\\\
T_{zz}(r,t)&=\int_{0}^{t}\mathcal{D}(t-\tau)\dfrac{\partial}{\partial\tau}\Pi^{\text{e}}_{\text{D}33}(r,\tau)\differential{\tau}-p(r,\theta,z,t)\\\
T_{\theta
z}(r,t)&=\int_{0}^{t}\mathcal{D}(t-\tau)\dfrac{\partial}{\partial\tau}\Pi^{\text{e}}_{\text{D}23}(r,\tau)\differential{\tau}+r\phi(t)\int_{0}^{t}\mathcal{D}(t-\tau)\dfrac{\differential}{\partial\tau}\Pi^{\text{e}}_{\text{D}33}(r,\tau)\differential{\tau}\\\
\end{split}$ (43)

Therefore, the governing equations reduce to:

$\left\\{\begin{split}&\dfrac{\partial T_{rr}(r,t)}{\partial
r}+\dfrac{T_{rr}(r,t)-T_{\theta\theta}(r,t)}{r}=0\\\ &\dfrac{\partial
T_{\theta\theta}(r,t)}{\partial\theta}=0\\\ &\dfrac{\partial
T_{zz}(r,t)}{\partial z}=0\end{split}\right.$ (44)

From the last two equations in Eqs. (44) we can conclude that the Lagrange
multiplier $p$ only depends on the spatial variable $r$ and at any time $t$
the governing problem reduces to a single Ordinary Differential Equation (ODE)
in the argument $r$:

$\left\\{\begin{split}\dfrac{\differential T_{rr}}{\differential
r}+\dfrac{T_{rr}-T_{\theta\theta}}{r}&=0\\\ T_{rr}&=0\qquad\text{at}\qquad
r=r_{o}\end{split}\right.$ (45)

Now, we restrict our attention to soft tissues whose elastic behaviour can be
considered hyperelastic. For such tissues a strain energy function $W$ can be
defined. Here, we choose $W$ in the form of the Mooney-Rivlin model:

$W(I_{1},I_{2})=(\frac{\mu_{0}}{2}-c_{2})(I_{1}-3)+c_{2}(I_{2}-3),$ (46)

where $I_{1}\\!=\\!\tr\textbf{B}$ and
$I_{2}=\nicefrac{{1}}{{2}}\left((\tr\textbf{B})^{2}-\tr\textbf{B}^{2}\right)$
are the first and second invariants of the tensor B, respectively, $\mu_{0}$
is the infinitesimal shear modulus and $c_{2}$ is the second Mooney-Rivlin
parameter. The Neo-Hookean model is recovered by setting $c_{2}=0$. The choice
for $W$ is motivated by many experimental observations on soft tissues. In
particular, it has been observed that the brain behaves as a Mooney-Rivlin
material in torsion [20] and in simple shear [21]. Moreover, the Mooney-Rivlin
model has the key feature of predicting a linear elastic response in torsion,
i.e. the torque required to twist the cylinder depends linearly on the strain.
The elastic Cauchy stress $\textbf{T}^{\text{e}}$ for an incompressible
material is given by the following relation:

$\textbf{T}^{\text{e}}=2W_{1}\textbf{B}-2W_{2}\textbf{B}^{-1}-p^{\text{e}}\textbf{I},$
(47)

where $W_{i}=\partial W/\partial I_{i}$, $i=\\{1,2\\}$ and $p^{\text{e}}$ is
the elastic Lagrange multiplier [25]. By combining Eqs. (46), (47) and (30) we
can calculate the components of the tensor $\bm{\Pi}^{\text{e}}_{\text{D}}$:

$\begin{split}\Pi^{\text{e}}_{\text{D}11}(r,t)&=\nicefrac{{1}}{{3}}(4c_{2}-\mu_{0})r^{2}\phi^{2}(t),\\\
\Pi^{\text{e}}_{\text{D}22}(r,t)&=-\nicefrac{{2}}{{3}}(c_{2}+2\mu_{0})r^{2}\phi^{2}(t)-\nicefrac{{1}}{{3}}(2c_{2}+\mu_{0})r^{4}\phi^{4}(t),\\\
\Pi^{\text{e}}_{\text{D}33}(r,t)&=-\nicefrac{{1}}{{3}}(2c_{2}+\mu_{0})r^{2}\phi^{2}(t),\\\
\Pi^{\text{e}}_{\text{D}23}(r,t)&=\mu_{0}r\phi(t)+\nicefrac{{1}}{{3}}(2c_{2}+\mu_{0})r^{3}\phi^{3}(t).\end{split}$
(48)

By combining Eqs. (43) and Eqs. (48) and substituting into Eq. (45) we obtain
an ODE for the variable $p$:

$\begin{split}\dfrac{\differential p}{\differential
r}&=\dfrac{r}{3}(14c_{2}+\mu_{0})\int^{t}_{0}\mathcal{D}(t-\tau)\dfrac{\differential}{\differential\tau}\phi^{2}(\tau)\differential\tau+\dfrac{r^{3}}{3}(2c_{2}+\mu_{0})\int^{t}_{0}\mathcal{D}(t-\tau)\dfrac{\differential}{\differential\tau}\phi^{4}(\tau)\differential{\tau}\\\
&-\phi(t)\left(2r\mu\int^{t}_{0}\mathcal{D}(t-\tau)\dfrac{\differential}{\differential\tau}\phi(\tau)\differential{\tau}+\dfrac{2}{3}r^{3}(2c_{2}+\mu)\int^{t}_{0}\mathcal{D}(t-\tau)\dfrac{\differential}{\differential\tau}\phi^{3}(\tau)\differential{\tau}\right)\\\
&+\phi^{2}(t)\dfrac{r^{3}}{3}(2c_{2}+\mu_{0})\int^{t}_{0}\mathcal{D}(t-\tau)\dfrac{\differential}{\differential\tau}\phi^{2}(\tau)\differential{\tau}\end{split}$
(49)

with the initial condition:

$p=\dfrac{r_{o}^{2}}{3}(4c_{2}-\mu_{0})\int_{0}^{t}\mathcal{D}(t-\tau)\dfrac{\differential}{\differential\tau}\phi^{2}(\tau)\differential{\tau}\qquad\text{at}\,\,r=r_{o}$
(50)

whose solution is:

$\begin{split}p&=\dfrac{1}{6}\left((14c_{2}+\mu_{0})r^{2}-3(2c_{2}+\mu_{0})r_{o}^{2}\right)\int^{t}_{0}\mathcal{D}(t-\tau)\dfrac{\differential}{\differential\tau}\phi^{2}(\tau)\differential{\tau}\\\
&+\dfrac{1}{12}(2c_{2}+\mu_{0})(r^{4}-r_{o}^{4})\int^{t}_{0}\mathcal{D}(t-\tau)\dfrac{\differential}{\differential\tau}\phi^{4}(\tau)\differential{\tau}\\\
&-\left(\mu_{0}(r^{2}-r_{o}^{2})\int^{t}_{0}\mathcal{D}(t-\tau)\dfrac{\differential}{\differential\tau}\phi(\tau)\differential{\tau}\right.\\\
&\left.+\dfrac{1}{6}(2c_{2}+\mu_{0})(r^{4}-r_{o}^{4})\int^{t}_{0}\mathcal{D}(t-\tau)\dfrac{\differential}{\differential\tau}\phi^{3}(\tau)\differential{\tau}\right)\phi(t)\\\
&+\left(\dfrac{1}{12}(2c_{2}+\mu_{0})(r^{4}-r_{o}^{4})\int^{t}_{0}\mathcal{D}(t-\tau)\dfrac{\differential}{\differential\tau}\phi^{2}(\tau)\differential{\tau}\right)\phi^{2}(t)\end{split}$
(51)

Finally, the components of the stress $\textbf{T}(r,t)$ can be obtained by
substituting Eq. (51) into Eqs. (43) to fully determine the final stress
distribution in the cylinder.

The torque $T(t)$ required to twist the cylinder can be computed as:

$T(t)=\int_{0}^{2\pi}\int_{0}^{r_{o}}T_{\theta
z}(r,t)r^{2}\differential{r}\differential{\theta},$ (52)

and the normal force $N$ necessary to keep to cylinder length constant reads:

$N(t)=\int_{0}^{2\pi}\int_{0}^{r_{o}}T_{zz}(r,t)r\differential{r}\differential{\theta}.$
(53)

The components $T_{\theta z}(r,t)$ and $T_{zz}(r,t)$ are given by Eqs. (43)
upon substituting Eqs. (48) and (51). The final expressions for the torque and
the normal force read:

$\begin{split}T(t)&=\dfrac{\pi}{2}\mu_{0}r_{o}^{4}\int^{t}_{0}\mathcal{D}(t-\tau)\dfrac{\differential}{\differential\tau}\phi(\tau)\differential{\tau}\\\
&+\dfrac{\pi}{9}(2c_{2}+\mu_{0})r_{o}^{6}\left(\int^{t}_{0}\mathcal{D}(t-\tau)\dfrac{\differential}{\differential\tau}\phi^{3}(\tau)\differential{\tau}-\phi(t)\int^{t}_{0}\mathcal{D}(t-\tau)\dfrac{\differential}{\differential\tau}\phi^{2}(\tau)\differential{\tau}\right)\end{split}$
(54)

and

$\begin{split}N(t)&=-\dfrac{\pi}{2}\mu_{0}r_{o}^{4}\phi(t)\int^{t}_{0}\mathcal{D}(t-\tau)\dfrac{\differential}{\differential\tau}\phi(\tau)\differential{\tau}-\dfrac{\pi}{4}(2c_{2}-\mu_{0})r_{o}^{4}\int^{t}_{0}\mathcal{D}(t-\tau)\dfrac{\differential}{\differential\tau}\phi^{2}(\tau)\differential{\tau}\\\
&+\dfrac{\pi}{18}(2c_{2}+\mu_{0})r_{o}^{6}\left(\int^{t}_{0}\mathcal{D}(t-\tau)\dfrac{\differential}{\differential\tau}\phi^{2}(\tau)\differential{\tau}\right)\phi^{2}(t)\\\
&-\dfrac{\pi}{9}(2c_{2}+\mu_{0})r_{o}^{6}\left(\int^{t}_{0}\mathcal{D}(t-\tau)\dfrac{\differential}{\differential\tau}\phi^{3}(\tau)\differential{\tau}\right)\phi(t)\\\
&+\dfrac{\pi}{18}(2c_{2}+\mu_{0})r_{o}^{6}\int^{t}_{0}\mathcal{D}(t-\tau)\dfrac{\differential}{\differential\tau}\phi^{4}(\tau)\differential{\tau},\end{split}$
(55)

respectively.
Note that Eqs. (54) and (55) are written with respect to the twist $\phi(t)$
which is a dimensional measure of the deformation. In view of comparing the
predictions of the QLV theory with those of the linear theory, it is useful to
rewrite (54) and (55) in terms of the strain $\gamma_{o}(t)$, defined in Eq.
(39), which is the strain at the outer surface of the cylinder and is a non-
dimensional measure of the deformation. Then, Eqs. (54) and (55) rewrite as
follows:

$\begin{split}T(t)&=\dfrac{\pi}{2}r_{o}^{3}\int^{t}_{0}\mu(t-\tau)\dfrac{\differential}{\differential\tau}\gamma_{o}(\tau)\differential{\tau}\\\
&+\dfrac{\pi}{9}(\dfrac{2c_{2}}{\mu_{0}}+1)r_{o}^{3}\left(\int^{t}_{0}\mu(t-\tau)\dfrac{\differential}{\differential\tau}\gamma_{o}^{3}(\tau)\differential{\tau}-\gamma_{o}(t)\int^{t}_{0}\mu(t-\tau)\dfrac{\differential}{\differential\tau}\gamma_{o}^{2}(\tau)\differential{\tau}\right)\end{split}$
(56)

and

$\begin{split}N(t)&=-\dfrac{\pi}{2}r_{o}^{2}\gamma_{o}(t)\int^{t}_{0}\mu(t-\tau)\dfrac{\differential}{\differential\tau}\gamma_{o}(\tau)\differential{\tau}-\dfrac{\pi}{4}(\dfrac{2c_{2}}{\mu_{0}}-1)r_{o}^{2}\int^{t}_{0}\mu(t-\tau)\dfrac{\differential}{\differential\tau}\gamma_{o}^{2}(\tau)\differential{\tau}\\\
&+\dfrac{\pi}{18}(\dfrac{2c_{2}}{\mu_{0}}+1)r_{o}^{2}\left(\int^{t}_{0}\mu(t-\tau)\dfrac{\differential}{\differential\tau}\gamma_{o}^{2}(\tau)\differential{\tau}\right)\gamma_{o}^{2}(t)\\\
&-\dfrac{\pi}{9}(\dfrac{2c_{2}}{\mu_{0}}+1)r_{o}^{2}\left(\int^{t}_{0}\mu(t-\tau)\dfrac{\differential}{\differential\tau}\gamma_{o}^{3}(\tau)\differential{\tau}\right)\gamma_{o}(t)\\\
&+\dfrac{\pi}{18}(\dfrac{2c_{2}}{\mu_{0}}+1)r_{o}^{2}\int^{t}_{0}\mu(t-\tau)\dfrac{\differential}{\differential\tau}\gamma_{o}^{4}(\tau)\differential{\tau}\end{split}$
(57)

where we have used the connection $\mathcal{D}(t)\\!=\\!\mu(t)/\mu_{0}$.
In the next section we will use Eqs. (56) and (57) to calculate the
predictions of the QLV model in the experimental scenarios of a step-strain
test and a ramp test.

##  5 Results

In view of predicting the relaxation behaviour of a tissue in simple torsion,
we consider two scenarios which are important from the experimental viewpoint:
the step-strain test and the ramp test. By using Eqs. (56) and (57), we then
derive the analytical expressions of the relaxation curves for the torque and
the normal force.

###  5.1 Small Deformations

We start by considering the torsion of a cylindrical tissue in the small
deformation regime and we calculate the torque and the normal force predicted
by the linear viscoelastic theory presented in Section 2. We assume that the
viscoelastic response of the tissue can be modelled as that of a generalised
Maxwell model, i.e. a system of spring and dash-pots arranged as shown in
Figure 2. For small deformations, the only non-zero component of the
infinitesimal strain tensor $\bm{\varepsilon}$ is $\varepsilon_{\theta
z}=\gamma(r,t)/2$, thus it follows that
$\text{dev}\bm{\varepsilon}=\bm{\varepsilon}$ and the constitutive equation
(6) reduces to the following equation:

$\sigma_{\theta
z}(r,t)=\int_{-\infty}^{t}\mu(t-\tau)\dfrac{\differential}{\differential\tau}\gamma(r,\tau)\differential{\tau}.$
(58)

The governing equations (42) are automatically satisfied. Eq. (58) is a one-
dimensional equation in the same form of Eq. (2), where $\mu(t)$ is now the
relaxation function. According to the formula
$T_{\text{lin}}=2\pi\int_{0}^{r_{o}}\sigma_{\theta z}r^{2}\differential r$,
the torque is then given by:

$T_{\text{lin}}(t)=\dfrac{\pi}{2}r_{o}^{3}\int_{-\infty}^{t}\mu(t-\tau)\dfrac{\differential}{\differential\tau}\gamma_{o}(\tau)\differential{\tau},$
(59)

where $\gamma_{o}(t)$ is the shear evaluated at the outer radius $r_{o}$, see
Eq. (39). Note that Eq. (59) is the linearised version of Eq. (56) for
$\gamma_{o}\ll 1$.

Following Section 2.2, we can calculate the analytical expressions of the
relaxation curve for the torque in response to a step and a ramp input,
respectively.
For the step-strain test, we consider the following form for the strain
$\gamma_{o}(t)$:

$\gamma_{o}(t)=\frac{r_{o}}{l}\alpha_{0}H(t)=\gamma_{o,0}H(t),$ (60)

where $\gamma_{o,0}=\frac{r_{o}}{l}\alpha_{0}$ is the amplitude of the step
and $H(t)$ is the Heaviside function as defined in Section 2.2. By
substituting Eq. (60) into Eq. (59) and upon integrating, we obtain:

$T^{\text{lin}}(t)=\frac{1}{2}\pi r_{o}^{3}\mu_{\text{step}}(t)\gamma_{o,0}$
(61)

where the relaxation function $\mu_{\text{step}}(t)$ is the following Prony
series:

$\mu_{\text{step}}(t)=\dfrac{2T^{\text{lin}}(t)}{\pi
r_{o}^{3}\gamma_{o,0}}=\mu_{\infty}+\sum_{i=1}^{n}\mu_{i}e^{-\frac{t}{\tau_{i}}}.$
(62)

The parameters $\mu_{i}$ in Eq. (62) are the shear moduli of the $n$ branches
of the generalized Maxwell model, see Figure 2, and
$\tau_{i}=\frac{\eta_{i}}{\mu_{i}}$ are the associated relaxation times. The
instantaneous shear modulus is obtained by evaluating the maximum of the
relaxation function $\mu_{\text{step}}(t)$ in $t=0$, as follows:

$\mu_{0}=\mu_{\text{ramp}}(0)=\mu_{\infty}+\sum_{i=1}^{n}\mu_{i}.$ (63)

We now consider a strain history in the form of a ramp, as follows:

$\gamma_{o}(t)=\frac{\gamma_{o,0}}{t^{*}}t-\frac{\gamma_{o,0}}{t^{*}}(t-t^{*})H(t-t^{*}).$
(64)

where $t^{*}$ is the rising time of the ramp. Substituting Eq. (64) into Eq.
(59) and computing the integral provides an analogous expression for the
torque as Eq. (61). The relaxation curve of the torque is now given by:

$\mu_{\text{ramp}}(t)=\frac{2T^{\text{lin}}(t)}{\pi
r_{o}^{3}\gamma_{o,0}}=\mu_{\infty}+\frac{1}{t^{*}}\sum_{i=1}^{n}\tau_{i}\mu_{i}e^{-\frac{t}{\tau_{i}}}\left(e^{\frac{t^{*}}{\tau_{i}}}-1\right),\qquad\text{for}\quad
t>t^{*}.$ (65)

Eq. (65) describes the relaxation curve of the torque in response to the ramp
function in Eq. (64). The right-hand side term in Eq. (65) is the response of
the generalized Maxwell system to a ramp input. Eqs. (62) and (65) can then be
used to estimate the viscoelastic parameters $\mu_{i},\mu_{\infty}$ and
$\tau_{i}$ by fitting the data from a step-test and a ramp-test, respectively.
The left-hand side of Eq. (65) can be computed from the experimental data,
i.e. the measured torque $T_{\text{lin}}$ and the imposed strain
$\gamma_{o,0}$ and the radius of the sample $r_{o}$. The right-hand side is
the analytical expression to be fitted in order to estimate the viscoelastic
parameters.
Note that the expression on the right-hand side of Eq. (65) does not depend on
the amount of shear $\gamma_{o,0}$, therefore Eq. (65) will predict accurate
results for those tissues that display the same relaxation response when
subjected to different levels of strain. Moreover, we define the following
non-dimensional function:

$\tilde{\mu}_{0\text{ramp}}=\dfrac{\mu_{\text{ramp}}(t^{*})-\mu_{\infty}}{\mu_{0}}=\frac{1}{t^{*}}\sum_{i=1}^{n}\dfrac{\tau_{i}\mu_{i}}{\mu_{0}}\left(1-e^{-\frac{t^{*}}{\tau_{i}}}\right)=\sum_{i=1}^{n}\dfrac{\nu_{i}^{-1}\mu_{i}}{\mu_{0}}\left(1-e^{-\nu_{i}}\right),$
(66)

which provides an estimate of the maximum value of the relaxation curve in Eq.
(65).

In the next section, we derive the corresponding expressions for the
relaxation curves of the torque and the normal force for the QLV model.

###  5.2 Large Deformations

In this section we derive the relaxation curves for the torque and the normal
force for a step test and a ramp test in the large deformation regime. We
first address the ramp test scenario and then by taking the limit case
$t^{*}\rightarrow 0$ we consider the step test scenario.

####  5.2.1 Torque

To compute the analytical expression of the relaxation curve for the torque,
we substitute the strain history (64) and the relaxation function (62) into
Eq. (56). Upon integrating we obtain:

$T(t,\gamma_{o,0})=\frac{1}{2}\pi
r_{o}^{3}\,\mu_{\text{ramp}}^{\text{QLV}}(t,\gamma_{o,0})\,\gamma_{o,0},\qquad\text{for}\quad
t>t^{*}$ (67)

where:

$\begin{split}\mu_{\text{ramp}}^{\text{QLV}}(t,\gamma_{o,0})=\mu_{\text{ramp}}(t)+\frac{2(1+\nicefrac{{2c_{2}}}{{\mu_{0}}})}{9t^{*3}}\sum_{i=1}^{n}\tau_{i}\mu_{i}e^{-\frac{t}{\tau_{i}}}&\Big{(}-2\tau_{i}(t^{*}+3\tau_{i})\\\
&+e^{\frac{t^{*}}{\tau_{i}}}(t^{*2}-4t^{*}\tau_{i}+6\tau_{i}^{2})\Big{)}\gamma_{o,0}^{2},\end{split}$
(68)

which is valid for $t>t^{*}$. Eq. (68) shows that the relaxation curve for the
torque predicted by the QLV model is the sum of the relaxation function
$\mu_{\text{ramp}}(t)$ in Eq. (65) and a term proportional to the square of
the imposed strain $\gamma_{o,0}$. Therefore, in the limit $\gamma_{o,0}\ll
1$, i.e. small deformation regime, Eq. (68) recovers the result in (65), which
is the relaxation function predicted by the linear theory for a ramp test.
A less straightforward comment is that the quadratic term in Eq. (68) arises
as a consequence of the time-dependent nature of the strain history. In other
words, if the tissue is deformed infinitely fast (as it is the case for a
perfect step-strain), Eq. (68) recovers the function $\mu_{\text{step}}(t)$ in
(62), without making the assumption of small deformations! By taking the limit
for $t^{*}\rightarrow 0$ in Eq. (68) we indeed obtain:

$\lim_{t^{*}\rightarrow
0}\mu_{\text{ramp}}^{\text{QLV}}(t,\gamma_{o,0})=\mu_{\infty}+\sum_{i=1}^{n}\mu_{i}e^{-\frac{t}{\tau_{i}}}=\mu_{\text{step}}(t),$
(69)

which is the relaxation curve predicted by the linear theory for a step test.
Moreover, opposite to the linear viscoelastic theory, the relaxation curve of
the torque predicted by the QLV model, i.e.
$\mu_{\text{ramp}}^{\text{QLV}}(t,\gamma_{o,0})$, depends on both the time and
the level of strain $\gamma_{o,0}$. This is in agreement with the original
assumption made by Fung when he first formulated the QLV theory [8].
With the aim of plotting the relaxation curve $\mu_{\text{ramp}}^{\text{QLV}}$
in Eq. (68), we now consider a simplified version of the generalised Maxwell
model with only one branch. This layout is also called Standard Linear Solid
model, and it is the simplest arrangement of elements which is able to
describe the behaviour of a viscoelastic solid. Note that the QLV theory,
despite taking into account large deformations, obeys the superposition
principle in time. Therefore, adding more branches to the generalized Maxwell
model will increase the accuracy of the model without producing mixed higher
order terms as, for example, in the multiple integral formulation [28]. We can
then set $n=1$ in Eq. (68) without loss of generality. Furthermore, we set
$c_{2}=c_{1}/2$ according to the observed values of $c_{1}$ and $c_{2}$ for
brain tissues [20].

Figure 6: Effect of the rising time $t^{*}$ on the relaxation curve
$\mu_{\text{ramp}}^{\text{QLV}}(t,\gamma_{o,0})/\mu_{0}$ in Eq. (68). The
parameter $\nu_{1}=t^{*}/\tau_{1}$ spans over $\\{0,0.2,0.4,0.6,0.8,1\\}$. The
following parameters have been fixed: $n=1$, $\tau_{1}=1$,
$\mu_{1}=\mu_{\infty}=1$, $c_{2}/\mu_{0}=2/3$ and the strain level
$\gamma_{o,0}=0.02$.

In Figure 6 we plot the relaxation curve
$\mu_{\text{ramp}}^{\text{QLV}}(t,\gamma_{o,0})$ in Eq. (68) for six different
ramp histories. To quantify the effect of the rising time on the profile of
the relaxation curve we vary the parameter $\nu_{1}=t^{*}/\tau_{1}$. We note
that Eq. (68) is valid for $t>t^{*}$. Experimentally, the ramp phase
($t<t^{*}$) is the noisy part of the data and cannot be used to perform the
model fitting. Therefore, we plot the ramp phase of the curves in Figure 6
with dashed lines. The dashed lines can be obtained by integrating Eq. (56)
with the strain history (64) in the time interval $0\leq t<t^{*}$.
From Figure 6 we can conclude that the faster the ramp phase, the higher the
peak of the relaxation curve. However, the limiting value of the curves as $t$
approaches $\infty$ is $\mu_{\infty}/\mu_{0}$ and is not affected by the
rising time $t^{*}$.

Figure 7: Effect of the final strain level $\gamma_{o,0}$ on the relaxation
curve $\mu_{\text{ramp}}^{\text{QLV}}(t,\gamma_{o,0})/\mu_{0}$ in Eq. (68).
The strain $\gamma_{o,0}$ spans over $\\{0.02,\pi/5,\pi/4,\pi/3,\pi/2,\pi\\}$.
The following parameters have been fixed: $n=1$, $\tau_{1}=1$,
$\mu_{1}=\mu_{\infty}=1$, $c_{2}/\mu_{0}=2/3$ and the parameter
$\nu_{1}=t^{*}/\tau_{1}=0.5$.

To quantify the effect of the strain level reached at the end of the ramp
phase, in Figure 7 we plot Eq. (68) for different values of $\gamma_{o,0}$ at
fixed $\nu_{1}=0.5$.

We now look at two limiting values of Eq. (68), namely the long-term
equilibrium and the instantaneous values of
$\mu_{\text{ramp}}^{\text{QLV}}(t,\gamma_{o,0})$. First, by taking the limit
for $t\rightarrow\infty$ we obtain:

$\lim_{t\rightarrow\infty}\mu_{\text{ramp}}^{\text{QLV}}(t,\gamma_{o,0})=\mu_{\infty},$
(70)

From Eq. (70) we observe that the equilibrium value of the relaxation curve of
the torque predicted by the QLV model is not affected by the quadratic terms
in $\gamma_{o,0}$. In other words, after an infinite time, i.e. at the
equilibrium, the response of the tissue is dominated by the long term shear
modulus $\mu_{\infty}$. Therefore, time and its non-linear effect through the
quadratic terms in Eq. (68) have no influence on the long-term modulus. This
effect is also observable from the horizontal asymptotes of Figures 6 and 7.

Then, we look at the instantaneous value of the function
$\mu_{\text{ramp}}^{\text{QLV}}(t,\gamma_{o,0})$, defined as
$\mu_{\text{ramp}}^{\text{QLV}}(t^{*},\gamma_{o,0})$. The point
$\mu_{\text{ramp}}^{\text{QLV}}(t^{*},\gamma_{o,0})$ is the maximum of the
relaxation curve. In the linear model of Eq. (62), the maximum of the
relaxation curve (obtained for a step-strain test) is equal to the
instantaneous shear modulus $\mu_{0}$ in Eq. (63). For a ramp test, the
instantaneous response is related to the instantaneous modulus through the
formula in Eq. (24), where the elastic constants $k_{i}$ are replaced by
$\mu_{i}$ for $i=\\{1,\dots,n\\}$, respectively.

We can now investigate how the non-linear terms in Eq. (68) affect the
instantaneous response of the torque. To quantify this effect, we compute the
following function:

$\begin{split}\tilde{\mu}_{0\text{ramp}}^{\text{QLV}}&=\dfrac{\mu_{\text{ramp}}^{\text{QLV}}(t^{*},\gamma_{o,0})-\mu_{\infty}}{\mu_{0}}\\\
&=\nu_{1}^{-1}\left(1-e^{-\nu_{1}}\right)+\frac{2}{9}(1+\nicefrac{{2c_{2}}}{{\mu_{0}}})\nu_{1}^{-3}e^{-\nu_{1}}\left(e^{\nu_{1}}\left(\nu_{1}^{2}-4\nu_{1}+6\right)-2(\nu_{1}+3)\right)\gamma_{o,0}^{2}\end{split}$
(71)

where $\nu_{1}=\nicefrac{{t^{*}}}{{\tau_{1}}}$.

In Figure 8, we plot Eq. (71) with respect to
$\nu_{1}=\nicefrac{{t^{*}}}{{\tau_{1}}}$ and for different values of
$\gamma_{o,0}$. The parameter $\nu_{1}$ spans from 0 to 1, where the value $0$
corresponds to perfect step input ($t^{*}=0$) and the value $1$ corresponds to
a ramp input when the rising time of the ramp is equal to the characteristic
relaxation time of the tissue ($t^{*}=\tau_{1}$).

Figure 8: Effect of the rising time $t^{*}$ and of the strain level
$\gamma_{o,0}$ on the maximum of the relaxation curve of the torque. The
function in Eq. (71) is plotted with respect to $\nu_{1}=t^{*}/\tau_{1}$ for
different values of strain
$\gamma_{o,0}=\\{0.02,\pi/5,\pi/4,\pi/3,\pi/2,\pi\\}$. The solid blue line is
obtained from the linear model, see Eq. (24) and does not depend on the level
of strain.

The curves in Figure 8 show that the higher $\gamma_{o,0}$, the higher is the
peak of the relaxation curve. For $\gamma_{o,0}=0.02$ (solid blue line), Eq.
(71) recovers the relaxation curve in Eq. (66). In the small deformation
regime $\gamma_{o,0}\ll 1$, the non-linear terms in Eq. (71) are very small
and as expected the predictions of the QLV model recover those of the linear
model.
When $t^{*}\rightarrow 0$ (i.e. $\nu_{1}\rightarrow 0$), the ramp is
infinitely fast and is very close to be a step. In this limit Eq. (68) reduces
to Eq. (62) and the linear and quasi-linear viscoelastic theories predict the
same results for any value of $\gamma_{o,0}$. This limit corresponds to the
point $(0,1)$ in Figure 8.

The effect of the rising time $t^{*}$ (and therefore of the strain-rate if
$\gamma_{o,0}$ is fixed) on the relaxation curve is accounted by the parameter
$\nu_{1}$: the slower the test, the lower the peak of the experimental
relaxation curve in Eq. (68) and the bigger is the difference between the
predictions of the linear and the QLV model.
When $t^{*}\rightarrow\infty$, i.e. for very slow (quasi-static) tests,
$\tilde{\mu}_{0\text{ramp}}^{\text{QLV}}\rightarrow 0$, which means that only
the quasi-static elastic properties of the material, i.e. $\mu_{\infty}$, can
be estimated. This would correspond to a quasi-static experiment.

Moreover, from Figure 8 we can quantify the influence of the quadratic term in
Eq. (71). For small values of $\nu_{1}$, the quadratic contribution is minimal
even for large values of strain. Its effect increases with both $\nu_{1}$ and
$\gamma_{o,0}$, as expected.

Note that $\tilde{\mu}_{0\text{ramp}}^{\text{QLV}}$ decreases as $\nu_{1}$
increases, but increases with $\gamma_{o,0}$. These two opposite effects can
lead to an apparent compensation, i.e. constant value of
$\tilde{\mu}_{0\text{ramp}}^{\text{QLV}}$, when slow experiments are performed
in the large deformation range.

Finally, we note that the second Mooney-Rivlin parameter $c_{2}$ enters Eq.
(68) through the quadratic term $\gamma_{o,0}^{2}$ of the function
$\mu_{\text{ramp}}^{\text{QLV}}(t,\gamma_{o,0})$. As we highlighted above, the
quadratic term vanishes in both the limits $t^{*}\rightarrow 0$ (step-strain
test) and $t\rightarrow\infty$ (elastic equilibrium). Therefore, it is
difficult to determine the parameter $c_{2}$ from the torque data. The
identification of $c_{2}$ requires information on the normal force.
In conclusion, the relaxation curve of the torque allows us to estimate the
instantaneous and the long-term moduli $\mu_{0}$ and $\mu_{\infty}$,
respectively, the moduli $\mu_{i}$ and the relaxation times $\tau_{i}$ from a
step or a ramp test.

####  5.2.2 Normal Force

In this section we derive the analytical expression of the relaxation curve
for the normal force predicted by the QLV model. As in the previous section,
we consider both the ramp test and the step-strain test scenarios.
Furthermore, we will show how to use the normal force data measured from the
two tests to estimate the second Mooney-Rivlin parameter $c_{2}$.
We recall that the only non-zero component of the infinitesimal strain tensor
is the shear component $\varepsilon_{\theta z}$. It follows that all three
normal components of the Cauchy stress tensor are zero, particularly
$\sigma_{zz}=0,\forall t$. Therefore, in the small deformation regime the
linear theory predicts a zero normal force response for any strain history
input.
On the other hand, in the large deformation regime the QLV theory predicts a
non-zero normal force for an incompressible Mooney-Rivlin viscoelastic
material under torsion, see Eq. (57). By substituting the deformation history
(64) into Eq. (57), we calculate the relaxation curve for the normal force in
response to a ramp input. Upon integrating, we obtain:

$\begin{split}N(t,\gamma_{o,0})=&-\frac{\pi}{4}r_{o}^{2}\left(\left(1+2c_{2}/\mu_{0}\right)\mu_{\infty}+2\sum_{i=1}^{n}\mu_{i}\nu_{i}^{-2}e^{-\frac{t}{\tau_{i}}}\left(2c_{2}/\mu_{0}\left(e^{\nu_{i}}\left(\nu_{i}-1\right)+1\right)+e^{\nu_{i}}-\nu_{i}-1\right)\right)\gamma_{o,0}^{2}\\\
&-\frac{\pi}{9}\left(1+2c_{2}/\mu_{0}\right)r_{o}^{2}\sum_{i=1}^{n}\mu_{i}\nu_{i}^{-4}e^{-\frac{t}{\tau_{i}}}\left(e^{\nu_{i}}\left(\nu_{i}^{2}-6\nu_{i}+12\right)-\nu_{i}^{2}-6\nu_{i}-12\right)\gamma_{o,0}^{4},\end{split}$
(72)

which is valid for $t>t^{*}$ and where $\nu_{i}=t^{*}/\tau_{i}$. From Eq. (72)
we observe that the normal force $N$ is given the sum of two terms
proportional to $\gamma_{o,0}^{2}$ and $\gamma_{o,0}^{4}$, respectively.
Therefore, in the small deformation limit ($\gamma_{o,0}\\!\rightarrow 0$) Eq.
(72) recovers the predictions of the linear theory, i.e. $N\\!=\\!0,\forall
t$. According to the definition in (13), the relaxation function associated to
the normal force is $\sigma_{zz}(t,\gamma_{o,0})/\varepsilon_{0}$, where
$\sigma_{zz}(t,\gamma_{o,0})=-N(t,\gamma_{o,0})/(\pi r_{o}^{2})$ and
$\varepsilon_{0}=\gamma_{o,0}/2$, namely:

$\mu_{N}(t,\gamma_{o,0})=-\frac{2N(t,\gamma_{o,0})}{\pi
r_{o}^{2}\gamma_{o,0}}.$ (73)

In the limit $t^{*}\rightarrow 0$ (step strain input), Eq. (73) reduces to:

$\lim_{t^{*}\rightarrow
0}\mu_{N}(t,\gamma_{o,0})=\frac{\gamma_{o,0}}{2}\,\left(1+2\frac{c_{2}}{\mu_{0}}\right)\mu(t),\quad\text{with}\quad\mu(t)=\mu_{\infty}+\sum_{i=1}^{n}\mu_{i}e^{-\frac{t}{\tau_{i}}}.$
(74)

The factor $\gamma_{o,0}/2$ represents the dependence on the deformation,
while the factor $(1+2c_{2}/\mu_{0})\mu(t)$ represents the dependence on time,
which is indeed strain-independent. Thus, it follows that different relaxation
curves obtained for different values of the strain $\gamma_{o,0}$ do not
overlap. The functional dependence of $\mu_{N}$ on the strain $\gamma_{o,0}$
is dictated by the form of the elastic constitutive model. In this particular
case, the choice of a Mooney-Rivlin model yields a linear dependence on
$\gamma_{o,0}$ in Eq. (74). However, the relaxation curves display the same
exponential decay in time, which is dictated by the choice of the rheological
model. For a generalised Maxwell model, the decay is exponential according to
the Prony series in Eq. (74). If the experimental curves do not display the
same decay in time, then the QLV model will not fit the data accurately and
more advanced/non-linear models should be considered [1].

Since $N(t)$ depends on the deformation through quadratic terms in
$\gamma_{o,0}$ and the aim here is to illustrate how the relaxation curve of
the normal force is affected by the level of deformation and the ramp phase,
we choose to introduce the following function:

$\begin{split}f_{N}(t,\gamma_{o,0})=&-\frac{2N(t)}{\pi
r_{o}^{2}\gamma_{o,0}^{2}}\\\
=&\frac{1+2c_{2}/\mu_{0}}{2}\mu_{\infty}+\sum_{i=1}^{n}\mu_{i}\nu_{i}^{-2}e^{-\frac{t}{\tau_{i}}}\left(2c_{2}/\mu_{0}\left(e^{\nu_{i}}\left(\nu_{i}-1\right)+1\right)+e^{\nu_{i}}-\nu_{i}-1\right)\\\
&+\frac{1+2c_{2}/\mu_{0}}{9}\sum_{i=1}^{n}\mu_{i}\nu_{i}^{-4}e^{-\frac{t}{\tau_{i}}}\left(e^{\nu_{i}}\left(\nu_{i}^{2}-6\nu_{i}+12\right)-\nu_{i}^{2}-6\nu_{i}-12\right)\gamma_{o,0}^{2},\end{split}$
(75)

which is valid for $t>t^{*}$. By taking the limit $t\rightarrow\infty$, Eq.
(75) reduces to:

$f_{N_{\infty}}=\lim_{t\rightarrow\infty}f_{N}(t,\gamma_{o,0})=\left(1/2+c_{2}/\mu_{0}\right)\mu_{\infty}.$
(76)

The limiting value in Eq. (76) explicitly depends on the second Mooney-Rivlin
coefficient $c_{2}$. Therefore, experimentally we can determine $c_{2}$ by
performing a ramp (or a step) test and by measuring the asymptotic value of
the normal force curve as $t\rightarrow\infty$. We point out that the
parameters $c_{2}$ appears explicitly in the limiting value (as
$t\rightarrow\infty$) of the normal force only (it does not appear in the
asymptotic value on the torque). This is a consequence of the fact that the
vertical force is associated to the change in area of the section of the
cylinder. Furthermore, this is consistent with the fact that $c_{2}$ is the
parameter associated to the second invariant $I_{2}$ of the strain tensor B,
which indeed accounts for the changes in the area of the material due to the
deformation.
Finally, to show the influence of the rising time $t^{*}$ and of the level of
strain at the end of the ramp $\gamma_{o,0}$ on the function $f_{N}$, we plot
$f_{N}/\mu_{0}$ for different values of $\nu_{1}$ and $\gamma_{o,0}$ in Figure
9a and Figure 9b, respectively.

Figure 9: Effect of the rising time (a) and of the strain level (b) on the
relaxation curve of the normal force in Eq. (75). The following parameters are
fixed: $n=1$, $\tau_{1}=1$, $\mu_{1}=\mu_{\infty}=1$, $c_{2}/\mu_{0}=2/3$. In
(a) we fix $\gamma_{o,0}=0.02$ and let $\nu_{1}$ spanning over
$\\{0,0.2,0.4,0.6,0.8,1\\}$. $\nu_{1}=0$ represents a perfect step test and
$\nu_{1}=1$ represents a ramp test with rising time $t^{*}=\tau_{1}$. In (b)
we fix $\nu_{1}=1/2$ and let
$\gamma_{o,0}=\\{0.02,\pi/5,\pi/4,\pi/3,\pi/2,\pi\\}$.

Figure 9b shows that the final level of strain reached at the end of the ramp
has a negligible effect on the relaxation curves of the normal force. However,
for a fixed strain level, the slower is the ramp, the lower is the peak of the
relaxation curve, similarly as we observed for the relaxation curves of the
torque (Figure 6).

In conclusion to fully characterise the viscoelastic behaviour of a soft
tissue in torsion (that obeys a Mooney-Rivlin hyperelastic law), according to
the QLV theory, we only need to perform a single step-strain test. This test
can be performed by using a rheometer that gives access to two sets of data,
the torque and the normal force. Moreover, as we showed in Figure 8 and Figure
9b, the level of strain does not affect the relaxation curves of the torque
and the normal force for a step test. Therefore, if the rheometer can achieve
high strain rates and if we are only interested in estimating
$\mu_{0},\mu_{\infty},\mu_{i}$ and $\tau_{i}$ (for instance if the material is
Neo-Hookean, i.e. $c_{2}=0$), the test can be performed in the small
deformation regime and the parameters can be estimated by fitting the torque
data with Eq. (62). Otherwise, the test should be modelled as a ramp-test in
the large deformation regime (see Section 5.2). In this case, from the
relaxation curve of the torque in Eq. (68) we can estimate the long-term shear
modulus $\mu_{\infty}$ from the value of the data as $t\rightarrow\infty$,
according to Eq. (70). We can obtain all the other moduli $\mu_{i}$ and the
relaxation times $\tau_{i}$ by fitting the relaxation curve in Eq. (68).
Finally, we can estimate the parameter $c_{2}$ from the value of the normal
force data as $t\rightarrow\infty$, according to Eq. (76).

##  6 Conclusions

In this Chapter we reviewed the foundations of linear viscoelasticity and the
theory of Quasi Linear Viscoelasticity (QLV). With the aim of providing a
fitting procedure for the QLV model and estimating the viscoelastic properties
of a soft tissue, we considered the torsion of a soft solid cylinder and we
wrote the governing equations of the viscoelastic problem for a tissue that
elastically behaves as a Mooney-Rivlin material, such as the brain. We derived
the analytical predictions of the relaxation curves for the torque and the
normal force necessary to twist a cylindrical sample. We considered two
experimental scenarios: the step test, where the tissue is instantaneously
deformed and held in position, and the ramp test, where the tissue is deformed
in a finite time and then held in position. These tests are commonly performed
to characterise the time-dependent properties of soft tissues and allow to
investigate their stress relaxation behaviour. We investigated the effect of
the strain level and rising time of the ramp on the relaxation curves of the
torque and the normal force. Our results show that in a step test, the linear
and the QLV models predict the same relaxation curves for the torque. However,
when the strain input is in the form of a ramp function, the non-linear terms
appearing in the QLV model affect the relaxation curve of the torque depending
on the strain level attained at the end of the loading phase (see Figure 8).
In particular, the higher is the strain level, the higher is the maximum of
the relaxation curve, whilst the equilibrium value remains constant and
unchanged.
The linear model predicts a zero normal force $\forall t$, whilst the QLV
model predicts a non-zero normal force that depends on $\gamma_{o,0}^{2}$ and
$\gamma_{o,0}^{4}$. Our results show that the relaxation curve of the normal
force depends on the level of strain both for a step test (see Eq. (74)) and a
ramp test (see Eq. (73)). Although the contributions of the non-linear terms
are negligible for a Mooney-Rivlin material (see Figure 9b), their effect
might be relevant for materials that obey a different elastic law (e.g. Fung,
Ogden, Gent, etc.).
Finally, our results provide useful guidelines to accurately fit QLV models in
view of estimating the viscoelastic properties of soft tissues. We showed how
to use the data from a ramp test in torsion (i.e. the relaxation curves of the
torque and the normal force) to estimate the constitutive parameters of a
Mooney-Rivlin viscoelastic tissue.

## References

  * [1] Christensen R, (2012) Theory of viscoelasticity: an introduction. Elsevier.
  * [2] Budday S, Sommer G, Haybaeck J, Steinmann P, Holzapfel G A, Kuhl E. (2017). Rheological characterization of human brain tissue. Acta biomaterialia, 60, 315-329.
  * [3] Rashid B, Destrade M, Gilchrist M D (2014). Mechanical characterization of brain tissue in tension at dynamic strain rates. Journal of the mechanical behavior of biomedical materials, 33, 43-54.
  * [4] Chatelin S, Constantinesco A, Willinger R (2010) Fifty years of brain tissue mechanical testing: from in vitro to in vivo investigations. Biorheology, 47:255–276.
  * [5] Pack E, Dubik J, Snyder W, Simon A, Clark S, De Vita R (2020) Biaxial Stress Relaxation of Vaginal Tissue in Pubertal Gilts. Journal of Biomechanical Engineering, 142(3).
  * [6] Shearer T, Parnell W J, Lynch B, Screen H R, Abrahams D I (2020) A Recruitment Model of Tendon Viscoelasticity That Incorporates Fibril Creep and Explains Strain-Dependent Relaxation. Journal of Biomechanical Engineering, 142(7).
  * [7] Lakes R, (2009) Viscoelastic Materials. Cambridge University Press.
  * [8] Fung Y C (2013) Biomechanics: mechanical properties of living tissues. Springer Science and Business Media.
  * [9] Rigby B J, Hirai N, Spikes J D, Eyring H (1959) The mechanical properties of rat tail tendon. The Journal of general physiology, 43(2), 265-283.
  * [10] Everett J S, Sommers M S. (2013). Skin viscoelasticity: physiologic mechanisms, measurement issues, and application to nursing science. Biological research for nursing, 15(3), 338-346.
  * [11] Pucci E, Saccomandi G (2012) On the nonlinear theory of viscoelasticity of differential type, Mathematics and mechanics of solids, 17(6), 624-630.
  * [12] Reese, S (2003) A micromechanically motivated material model for the thermo-viscoelastic material behaviour of rubber-like polymers. International Journal of Plasticity, 19(7), 909-940.
  * [13] Reese S, Govindjee S (1998) A theory of finite viscoelasticity and numerical aspects. International journal of solids and structures, 35(26-27), 3455-3482.
  * [14] Lion A (1997) A physically based method to represent the thermo-mechanical behaviour of elastomers. Acta Mechanica, 123(1), 1-25.
  * [15] Yazdi S J M, Cho K S, Kang N (2018). Characterization of the viscoelastic model of in vivo human posterior thigh skin using ramp-relaxation indentation test. Korea-Australia Rheology Journal, 30(4), 293-307.
  * [16] Green A E, Rivlin R S (1957). The mechanics of non-linear materials with memory. Archive for Rational Mechanics and Analysis, 1(1), 1-21.
  * [17] Findley W N, Davis F A. (2013). Creep and relaxation of nonlinear viscoelastic materials. Courier Corporation.
  * [18] Valanis K C. (2014). Irreversible thermodynamics of continuous media: internal variable theory (Vol. 77). Springer.
  * [19] De Pascalis R, Abrahams I D, Parnell W J (2014) On nonlinear viscoelastic deformations: a reappraisal of Fung’s quasi-linear viscoelastic model. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 470(2166), 20140058.
  * [20] Balbi V, Trotta A, Destrade M, Annaidh A N (2019) Poynting effect of brain matter in torsion. Soft matter, 15(25), 5147-5153.
  * [21] Rashid B, Destrade M, Gilchrist M D (2013) Mechanical characterization of brain tissue in simple shear at dynamic strain rates. Journal of the mechanical behavior of biomedical materials, 28, 71-85.
  * [22] Boltzmann L (1874) Zur Theorie der elastischen Nachwirkungen. Sitzungsberichte der Kaiserlichen Akademie der Wissenschaften. Mathematisch-Naturwissenschaftliche Classe, 70 (II) 275.
  * [23] Markovitz H (1977) Boltzmann and the beginnings of linear viscoelasticity. Transactions of the Society of Rheology, 21(3), 381-398.
  * [24] Falcioni M, Vulpiani A (2015) Ludwig Boltzmann: a tribute on his 170th birthday. Lettera Matematica, 2(4), 171-183.
  * [25] Ogden R W (1997) Non-linear elastic deformations. Courier Corporation.
  * [26] Balbi V, Shearer T, Parnell W J (2018) A modified formulation of quasi-linear viscoelasticity for transversely isotropic materials under finite deformation. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 474(2217), 20180231.
  * [27] Gilchrist M D, Rashid B, Murphy J G, Saccomandi G (2013) Quasi-static deformations of biological soft tissue. Mathematics and Mechanics of Solids, 18(6), 622-633.
  * [28] Lockett F J (1972) Nonlinear viscoelastic solids. Academic Press.
  * [29] Afsar-Kazerooni N, Srinivasa A R, Criscione J C (2020) Experimental investigation of the inelastic response of pig and rat skin under uniaxial cyclic mechanical loading. Experimental Mechanics, 1-17.
  * [30] Kazerooni N A, Srinivasa A R, Criscione J C (2020) A multinetwork inelastic model for the hysteretic response during cyclic loading of pig and rat skin. International Journal of Non-Linear Mechanics, 126, 103555.