Reactive Viscoelasticity Curvefitting Problem
Overview
The objective of this case study is to provide a representative problem of curvefitting a viscoelastic response to the reactive viscoelasticity material model in FEBio. Reactive viscoelasticity is described in two prior papers, Ateshian (2015) and Ateshian et al. (2022). Representative experimental data for the tensile response of a sample of human glenohumeral joint articular cartilage was obtained from the study of Huang et al. (2005). Model files for this case study can be downloaded from the FEBio Studio Repository, in the Tutorial folder.
Reactive Viscoelasticity in FEBio
The reative viscoelasticity framework in FEBio is based on constrained reactive mixture theory. It accommodates anisotropic nonlinear viscoelasticity and is well-suited for modeling fibrous biological tissues where fibers only sustain tension. The strain energy density of a reactive viscoelastic material is given by
$\Psi_{r}\left(\mathbf{F}\left(t\right)\right)=\Psi_{r}^{e}\left(\mathbf{F}\left(t\right)\right)+\sum_{u}w^{u}\Psi_{0}^{a}\left(\mathbf{F}^{u}\left(t\right)\right)\,.$ (1)
Here, $\mathbf{F}\left(t\right)$ is the observable deformation gradient at the current time $t$, and $t^{u}$ represents the time $t$ when all currently reformed (non-breaking) weak bonds break and start reforming into a stress-free state; we refer to it as bond generation $u$, with $-\infty<t^{u}\le t$. The deformation gradient $\mathbf{F}^{u}\left(t\right)$ of generation $u$ is a suitably chosen deformation gradient relative to the time that $u-$generation bonds are reformed, based on a constitutive modeling assumption. The function $\Psi_{0}^{a}$ in eq.(1) represents the strain energy density of solid matter associated with weak bonds, whereas $\Psi_{r}^{e}$ is that of strong (elastic) bonds. The mass fraction $w^{u}$ of each weak bond generation $u$ increases from generation $u$ to the subsequent generation $v$ as those bonds are reforming in a stress-free state; this increase takes place during the time interval $t^{u}\le t<t^{v}$. In particular, at time $t^{u}$, $\mathbf{F}^{u}\left(t^{u}\right)$ must be a proper orthogonal tensor in order to satisfy the stress-free reformation of weak bonds, $\Psi_{0}^{a}\left(\mathbf{F}^{u}\left(t^{u}\right)\right)=0$.
Loading that takes place at generation $v$ leads to a subsequent decay of $w^{u}$ with time $t\ge t^{v}$ due to the evolving bond-breaking-and-reforming reaction, as governed by a user-defined reduced relaxation function $g\left(\mathbf{F}\left(t^{v}\right),t-t^{v}\right)$ for $t\ge t^{v}$. By definition, the reduced relaxation function must satisfy $\lim_{t\to0}g\left(\mathbf{F},t\right)=1$, whereas $\lim_{t\to\infty}g\left(\mathbf{F},t\right)=0$. To accommodate nonlinear viscoelasticity this relaxation function may depend on some suitable measure of the state of strain at the time of bond breaking, as embodied in the dependence of $g$ on $\mathbf{F}\left(t^{v}\right)$. The general recursive relation for the bond mass fractions $w^{u}$ is
$w^{u}\left(t\right)=f^{u}\left(t^{v}\right)g\left(\mathbf{F}\left(t^{v}\right),t-t^{v}\right)\,,\quad t^{v}\ge t\,,$ (2)
where
$f^{u}\left(t\right)=1-\sum_{\gamma<u}w^{\gamma}\left(t\right)\,.$ (3)
The Cauchy stress tensor $\boldsymbol{\sigma}$ of a reactive viscoelastic material may be evaluated from the strain energy density $\Psi_{r}$ in eq.(1) using the standard hyperelasticity relation
$\boldsymbol{\sigma}=2J^{-1}\mathbf{F}\cdot\frac{∂\Psi_{r}}{∂\mathbf{C}}\cdot\mathbf{F}^{T}\,,$ (4)
thus producing the relation
$\boldsymbol{\sigma}\left(\mathbf{F}\right)=\boldsymbol{\sigma}^{e}\left(\mathbf{F}\right)+\sum_{u}w^{u}J^{-1}\left(t^{u}\right)\boldsymbol{\sigma}_{0}^{a}\left(\mathbf{F}^{u}\right)\,,$ (5)
where $J^{u}=\det\mathbf{F}$. Thus, the complete formulation of a reactive viscoelastic material in FEBio requires the specification of a constitutive model for the strong bond stress $\boldsymbol{\sigma}^{e}$, the weak bond stress $\boldsymbol{\sigma}_{0}^{a}$, and the reduced relaxation function $g\left(\mathbf{F}\left(t^{v}\right),t-t^{v}\right)$. This formulation implies that the weak bond constitutive model need not be the same as that of the strong bonds. It also implies that any desired relaxation function may be used to describe the time-dependent viscoelastic response. This model thus stands in contrast with the standard quasilinear viscoelastic model described in the FEBio Theory Manual, which always employs exponential relaxation (in a Prony series) and assumes that the viscous part of the response is proportional to the elastic part, using scalar proportionality constants (called $\gamma_{i}$ in the FEBio formulation).
Experimental Data
The tensile testing experiments on human articular cartilage described in Huang et al. (2005) consisted of a series of stress relaxation responses to prescribed step strains of 2%, 4%, 6%, 8%, 10%, 12%, 14%, 16%, 20%, and 25%. Based on the earlier study by Huang et al. (2001), it was assumed that the tensile stress-relaxation response of cartilage was due mostly to intrinsic viscoelasticity of its solid matrix rather than flow-dependent viscoelasticity. Therefore, in this case study, we choose to only fit a reactive viscoelastic model to those data, without accounting for the biphasic nature of articular cartilage. The experimental data used in this case study is presented in the figure below. The load response has been divided by the initial cross-sectional area of the sample cartilage strip to produce the engineering stress as a function of time.
FEBio Stress-Relaxation Model
The FEBio model needed to replicate this experimental response uses a unit cube (one hex8 element) to represent the gauge section of the long cartilage strip used in the experimental study (see Figure 3 below). The cube is aligned with the XYZ-coordinate axes, and zero displacements are prescribed normal to each of the coordinate planes to produce a homogeneous state of strain and stress in this model. On the positive Z-face (Z=1), a displacement profile is prescribed which replicates the experimental step strain sequence (Figure 2). The displacement load curve uses a linear function to ramp up the prescribed displacement in one second each time that a step strain is applied, to allow FEBio to take smaller time increments that produce smaller changes in displacement in case the solver fails to achieve convergence following a relatively large step increase in displacement.
Since large strains are prescribed in this analysis (engineering strains up to 25%), the strain is prescribed on a rigid platen (Figure 3) that uses the sliding-elastic contact interface with the tension flag turned on. This configuration allows the lateral sides of the cube to shrink as needed when the sample undergoes axial extension. Then, the engineering stress is evaluated from the force exerted by this rigid platen on the cube. Since the cube has an initial cross-sectional area of unity, the resulting force is also equal to the engineering stress.
Material Models
Choosing the correct material models when curve-fitting experimental data is an iterative process. Here, we assume that we have already settled on the right choice of material models and we seek to find suitable material properties using parameter optimization. Evidently, if one doesn’t know which material model is going to work best, one should try a set of material models, using the methodology described here, then modify this choice until one achieves a satisfactory curvefitting outcome.
For this particular case study, we have opted to model cartilage as a fiber-reinforced ground matrix, using the simplifying assumptions that (a) the fibers are unidirectional and aligned with the loading direction, and (b) only the fibers are viscoelastic, whereas (c) the ground matrix is isotropic and elastic. The properties of the ground matrix are taken from the original study of Huang et al. (2005), who performed compressive testing on cartilage cylindrical plugs adjacent to the cartilage strips used for tensile testing. Based on their study, we use the Holmes-Mow isotropic elastic material model for the ground matrix, with $E=0.141\,\text{MPa}$, $\nu=0$, and $\beta=0.81$.
For the strong bonds and weak bonds of the reactive viscoelastic fiber materials, we choose the fiber-exp-pow-linear material model, which requires up to four material parameters: $E$ which is Young’s modulus in the linear range of this toe-linear model; $\alpha$ and $\beta$ which respectively regulate the exponential and power parts of the toe region’s response, and $\lambda_{0}$ which determines the stretch ratio for the transition between the toe and linear regimes. We also need to select a reduced relaxation function that will suitably fit the transient response of this cartilage material. By trial and error we settled on the continuous relaxation spectrum function proposed by Fung (1972)1Y. C Fung, Nicholas Perrone, M Anliker. Biomechanics, its foundations and objectives. Prentice-Hall, 1972, which is implemented in FEBio as a special case of the Malkin reduced relaxation function. The material parameters for this relaxation function are $\tau_{1}$, $\tau_{2}$ and $\beta$, and setting $\beta=1$ recovers the desired special case.
In principle, we thus need to fit 11 material parameters to produce a suitable material model for the viscoelastic fibers. These parameters are $E^{e}$, $\alpha^{e}$, $\beta^{e}$ and $\lambda_{0}^{e}$ for the elastic bonds, $E^{a}$, $\alpha^{a}$, $\beta^{a}$ and $\lambda_{0}^{a}$ for the weak bonds, and $\tau_{1}$, $\tau_{2}$ and $\beta$ for the reduced relaxation function. However, we have already chosen to set $\beta=1$ in the Malkin relaxation function. We can also fix the value of $\tau_{1}$ on the basis that the experimental time-versus-stress data used in this case study (Figure 1) employs time increments ≥ 1 s. Therefore, there is no opportunity to explore relaxation responses that depend on shorter time scales, hence we can safely set $\tau_{1}=0.1\,\text{s}$ and not worry about finding an optimal value for it. These two choices reduce our list of unknown material parameters to nine.
It is important to understand that not all these decisions were made a priori. Many were the result of trial and error, and finally settling on this decision.
Curvefitting
A brute-force approach, whereby we attempt to fit nine parameters to the transient data of Figure 1, is not practical until we provide reasonable estimates for the values to be used as initial guesses. This can be done as suggested in the next two sections.
Initial Guess for Properties of the Strong Bond Material
We start by extracting the engineering stress versus engineering strain achieved at the last time step of each stress-relaxation response (orange points in Figure 1). We treat this stress-strain response as an approximation to the equilibrium strong bond response of the fibers and ground matrix. Then we create an elastic material model (a solid mixture) that only includes the Holmes-Mow ground matrix (whose properties are fixed) and a single fiber-exp-pow-linear material model representing the strong bond response.
<material id="2" name="tissue" type="solid mixture">
<density>1</density>
<solid type="Holmes-Mow">
<E>0.141</E>
<v>0</v>
<beta>0.81</beta>
</solid>
<solid type="fiber-exp-pow-linear">
<E>10.0</E>
<alpha>3.0</alpha>
<beta>3</beta>
<lam0>1.10</lam0>
<fiber type="vector">0,0,1</fiber>
</solid>
</material>We set this analysis to run from time $t=0$ to $t=1$ in 100 time increments and prescribe the platen displacement to ramp up linearly from $0$ to $0.25$, thus producing up to 25% engineering strain. We perform a parameter optimization on this model using the settings shown below. Here, initial guesses, their ranges and scales are set arbitrarily. If the parameter optimization hits one of the lower (min) or upper (max) values of the range(s) specified in this optimization file, the corresponding range is extended and the analysis is restarted using the last optimization solution as initial guesses for the next optimization analysis. Note that the data-fit model used in this analysis uses engineering stress versus engineering strain (not engineering stress versus time), where engineering stress = platen force and engineering strain = platen displacement.
<Options type="levmar">
<obj_tol>1e-4</obj_tol>
<f_diff_scale>1e-4</f_diff_scale>
<print_level>PRINT_VERBOSE</print_level>
</Options>
<Parameters>
<param name="fem.material('tissue').solid[1].E">10.0,1,20,5</param>
<param name="fem.material('tissue').solid[1].alpha">0.1,0.001,20,0.1</param>
<param name="fem.material('tissue').solid[1].beta">3,2.01,4,2</param>
<param name="fem.material('tissue').solid[1].lam0">1.10,1.06,1.5,1</param>
</Parameters>
<Objective type="data-fit">
<fnc type="parameter">
<ordinate name="fem.rigidbody('platen').position.z"/>
<param name="fem.rigidbody('platen').Fz"/>
</fnc>
<data>
<point>0.02,0.0496649349240781</point>
<point>0.04,0.115826043383948</point>
<point>0.06,0.191439650759219</point>
<point>0.08,0.285956127982646</point>
<point>0.10,0.405676431670282</point>
<point>0.12,0.556903646420824</point>
<point>0.14,0.745934472885032</point>
<point>0.16,0.957021839479393</point>
<point>0.20,1.37919444468547</point>
<point>0.25,1.94314176572668</point>
</data>
</Objective>The results of this parameter optimization ($R^{2}=0.999915$) produce $E^{e}=11.51\,\text{MPa}$, $\alpha^{e}=3.892$, $\beta^{e}=2.187$ and $\lambda_{0}^{e}=1.142$. We then verify that these optimal parameters produce a good fit by running the strong bond model using these properties (Figure 4).
Initial Guess for Properties of the Weak Bond Material
The objective of this analysis is to predict the peak stress at the start of each step-relaxation response (yellow symbols in Figure 1). Since the peak response represents the combination of strong and weak bond contributions, we modify the model file used in the previous section to include one more solid in the mixture, to represent the weak bond material. For the strong bond material we use the material properties obtained from the parameter optimization completed in the previous section. For the weak bond material we will start with some initial guesses.
<material id="2" name="tissue" type="solid mixture">
<density>1</density>
<solid type="Holmes-Mow">
<E>0.141</E>
<v>0</v>
<beta>0.81</beta>
</solid>
<solid type="fiber-exp-pow-linear">
<E>11.51</E>
<alpha>3.892</alpha>
<beta>2.187</beta>
<lam0>1.142</lam0>
<fiber type="vector">0,0,1</fiber>
</solid>
<solid type="fiber-exp-pow-linear">
<E>10</E>
<alpha>3</alpha>
<beta>3</beta>
<lam0>1.1</lam0>
<fiber type="vector">0,0,1</fiber>
</solid>
</material>Now we use this material model to perform a parameter optimization only for the material coefficients of the weak bond material (the third solid in this mixture). The parameter optimization file takes on the following form:
<Options type="levmar">
<obj_tol>1e-4</obj_tol>
<f_diff_scale>1e-4</f_diff_scale>
<print_level>PRINT_VERBOSE</print_level>
</Options>
<Parameters>
<param name="fem.material('tissue').solid[2].E">10.0,1,20,5</param>
<param name="fem.material('tissue').solid[2].alpha">0.1,0.001,20,0.1</param>
<param name="fem.material('tissue').solid[2].beta">3.0,2.01,4,2</param>
<param name="fem.material('tissue').solid[2].lam0">1.10,1.06,1.5,1</param>
</Parameters>
<Objective type="data-fit">
<fnc type="parameter">
<ordinate name="fem.rigidbody('platen').position.z"/>
<param name="fem.rigidbody('platen').Fz"/>
</fnc>
<data>
<point>0.02,0.131579498915401</point>
<point>0.04,0.232396932754881</point>
<point>0.06,0.352117236442516</point>
<point>0.08,0.478138496746204</point>
<point>0.1,0.629365711496746</point>
<point>0.12,0.771140427331887</point>
<point>0.14,1.01688199132321</point>
<point>0.16,1.27837701084599</point>
<point>0.2,1.97464654880694</point>
<point>0.25,2.82844425813449</point>
</data>This parameter optimization ($R^{2}=0.997715$) produces $E^{a}=6.406\,\text{MPa}$, $\alpha^{a}=1.467$, $\beta^{a}=2.01$ and $\lambda_{0}^{a}=1.208$ and the resulting fit is shown in Figure 5 below.
Curvefitting of Entire Transient Response
We are finally ready to perform a parameter optimization for fitting all eight material constants to the transient data of Figure 1. The material model employed in this final step represents the full viscoelastic model for the this tissue.
<material id="2" name="tissue" type="solid mixture">
<density>1</density>
<solid type="Holmes-Mow">
<E>0.141</E>
<v>0</v>
<beta>0.81</beta>
</solid>
<solid type="reactive viscoelastic">
<density>1</density>
<kinetics>1</kinetics>
<trigger>0</trigger>
<wmin>0</wmin>
<emin>0.0002</emin>
<elastic type="fiber-exp-pow-linear">
<E>11.51</E>
<alpha>3.892</alpha>
<beta>2.187</beta>
<lam0>1.142</lam0>
<fiber type="vector">0,0,1</fiber>
</elastic>
<bond type="fiber-exp-pow-linear">
<E>6.406</E>
<alpha>1.467</alpha>
<beta>2.010</beta>
<lam0>1.208</lam0>
<fiber type="vector">0,0,1</fiber>
</bond>
<relaxation type="relaxation-Malkin">
<tau1>0.1</tau1>
<tau2>1350</tau2>
<beta>1.0</beta>
</relaxation>
</solid>
</material>For this transient analysis the list of parameters appearing in the optimization file takes the form
<param name="fem.material('tissue').solid[1].elastic.E">11.51,1,20,5</param>
<param name="fem.material('tissue').solid[1].elastic.alpha">3.9,0.001,20,0.1</param>
<param name="fem.material('tissue').solid[1].elastic.beta">2.2,2.01,4,2</param>
<param name="fem.material('tissue').solid[1].elastic.lam0">1.14,1.06,1.5,1</param>
<param name="fem.material('tissue').solid[1].bond.E">6.41,1,100,10</param>
<param name="fem.material('tissue').solid[1].bond.alpha">1.47,0.02,1000,100</param>
<param name="fem.material('tissue').solid[1].bond.beta">2.01,2.01,4,2</param>
<param name="fem.material('tissue').solid[1].bond.lam0">1.21,1.001,1.5,1</param>
<param name="fem.material('tissue').solid[1].relaxation.tau2">1350,1,4000,1000</param>At the completion of this parameter optimization we find that the regression coefficient is $R^{2}=0.997971$ and the final optimal material parameters are given by $E^{e}=10.93\,\text{MPa}$, $\alpha^{e}=0$, $\beta^{e}=2.890$, $\lambda_{0}^{e}=1.123$, $E^{a}=75.72\,\text{MPa}$, $\alpha^{a}=182.4$, $\beta^{a}=2.051$, $\lambda_{0}^{a}=1.038$, and $\tau_{2}=1141\,\text{s}$, with $\tau_{1}=0.1\,\text{s}$ and $\beta=1$. The resulting fit is shown in Figure 6 below.
Though this fit is not perfect, it still represents a remarkable outcome given the large number and range of strain increments. Importantly, this fit was produced with a reduced relaxation function which is not dependent on strain, implying that the response observed for this cartilage sample represents quasilinear (not nonlinear) viscoelasticity. It is interesting to note that the final optimization solution produced results for the strong bonds that were reasonably close to those obtained from the initial guess analysis (Figure 7), whereas the results produced for the weak bonds (Figure 8) is substantially different. This outcome may be explained by the fact that the experimental data has time increments no smaller than one second, implying that the prescribed strain at each step was not truly instantaneous, whereas the fitted weak bond properties reflect the true response of the weak bonds to an instantaneous jump in the strain.
Summary
The objective of this case study was to provide a representative problem of curvefitting a viscoelastic response to the reactive viscoelasticity material model in FEBio. A sample data set for the response of human articular cartilage to sequential tensile stress relaxations, over a range of strains spanning from 2% to 25%, was used in this case study. The procedure shown in this case study representing a suggestion on how to perform this type of curvefitting on experimental data. Results also illustrate the ability of the reactive viscoelastic framework to fit complex experimental data from biological soft tissues, using the wide range of material models available in FEBio.
Footnotes
- 1Y. C Fung, Nicholas Perrone, M Anliker. Biomechanics, its foundations and objectives. Prentice-Hall, 1972