Abstract

In this erratum, we correct a mistake in a subcomponent of the numerical algorithm proposed in our recent study for modeling anisotropic reactive nonlinear viscoelasticity (doi:10.1115/1.4054983), for the special case where multiple weak bond families may be recruited with loading. This correction overcomes a nonphysical response noted under uni-axial cyclical loading.

1 Introduction

In our recently published study [1], we reported on a numerical scheme for the framework of reactive viscoelasticity initially presented in our prior theoretical report [2]. We also extended that original theoretical framework to incorporate multiple families of weak bonds, which break and reform into a stress-free state under loading, to be optionally recruited with increasing amount of loading (Sec. 2.3 of the paper [1]). We explained that this weak-bond recruitment could be easily incorporated into the existing framework via Eq. (2.31) of the paper. We provided a pseudo-code that incorporated Eq. (2.31) in Algorithm 2 of the paper. Finally, we illustrated the application of weak-bond recruitment in the curve-fitting results of Fig. 9 of the paper, which analyzed the response of the tissue to consecutive ramp-and-hold displacement profiles.

Figure 9
Fig. 9 Curve-fitting of nucleus pulposus unconfined compression stress-relaxation response to data from Ref. [28]. (a) Experimental data and reactive viscoelasticity curvefit of the compressive engineering stress, for five consecutive increases of 5% in the compressive engineering strain. (b) Dependence of weak bond recruitment F and relaxation time τ2 on distortional strain K2. Symbols represent the values of K2 corresponding to each engineering strain increment in the stress-relaxation response.
Figure 9
Fig. 9 Curve-fitting of nucleus pulposus unconfined compression stress-relaxation response to data from Ref. [28]. (a) Experimental data and reactive viscoelasticity curvefit of the compressive engineering stress, for five consecutive increases of 5% in the compressive engineering strain. (b) Dependence of weak bond recruitment F and relaxation time τ2 on distortional strain K2. Symbols represent the values of K2 corresponding to each engineering strain increment in the stress-relaxation response.
Close modal

Upon release of the code in the open-source finite element software FEBio [3], we got feedback from a user (see Acknowledgments below) that this algorithm produced nonphysical responses when a tissue is subjected to uni-axial prescribed cyclical displacement with constant amplitude, such as load amplification along the loading direction, instead of the expected load attenuation. Here, we present a correction to Eq. (2.31), and an update to Algorithm 2 and Fig. 9, which overcome this original limitation.

2 Correction

The strain energy density of a reactive viscoelastic material was given in Eq. (2.1) of the paper. To account for weak-bond recruitment, we replace the original Eq. (2.31), which had proposed a modification to the function fu(tv) used in the calculation of the bond mass fraction wu(t) in Eq. (2.2), with a modification of Eq. (2.1),
Ψr(F(t))=Ψre(F(t))+uwu(t)F(Ξ(t))Ψ0a(Fu(t))
(2.31)

where Ξ is the strain measure described in Eq. (2.29) of the paper, and F(Ξ) represents the enhancement in weak bond mass fraction relative to the zero strain state, as illustrated in the example constitutive model of Eq. (2.30). The main difference between this corrected model and the original model is that the new model accounts for the enhancement F(Ξ(t)) in the weak bond mass fraction based on the strain measure Ξ evaluated at the current time t. The original formulation employed the maximum strain measure achieved in the strain history of the material, up until the first time-step tv when weak bonds from generation u started breaking.

This formulation adjusts the availability of weak bonds as the strain varies; in particular, it accommodates cyclical strains.

Table

Algorithm 2. The Function Update is called at each iteration of the nonlinear solver until convergence is achieved at the current time tn+1. The Function Ψr is called when evaluating the mixture free energy density at each iteration of the nonlinear solver. This function serves as a template for the similar calculations of the stress σ and elasticity C as described in Eq. (2.20) of Sec. 2.1.3. When u =0 in this function, u − 1 is equivalent to the superscripted in the text, see Eq. (2.23). These functions depend on those presented in Algorithm 1.

Function update
  Letm=number of generations
  If (m =0) Or (tn+1>tm1)
  If NewGeneration is True
   Pushtn+1tmontotu+1stack
   PushU(tn+1)ontoUu+1stack
   //Implement Sec. 2.3
   EvaluateΞ(U(tn+1))using Eq. (2.29)
   PushF(Ξ(U(tn+1)))
   ontoFu+1stack
   Pushfm1(tm)=
   ReformingBondMassFraction
   ontofu+1stack
   Call CullGenerations
Else iftn+1=tm
  Update previously pushed
  U(tm), F(Ξ(U(tm)))andfm1(tm)
FunctionΨr
  Letm=number of generations
  Ifm =0 ReturnΨr=Ψre(F(tn+1))+Ψ0a(F(tn+1))
  LetΨr=Ψre(F(tn+1))
  Foru =0 tom − 1
   Evaluatewu1(tn+1)=fu1(tu)g(F(tu),tn+1tu)
   EvaluateFu1(tn+1)=F(tn+1)·U1(tu1)
   EvaluateΨr+= wu1(tn+1)F(Ξ(U(tn+1)))Ψ0a(Fu1(tn+1))
  ReturnΨr
Function update
  Letm=number of generations
  If (m =0) Or (tn+1>tm1)
  If NewGeneration is True
   Pushtn+1tmontotu+1stack
   PushU(tn+1)ontoUu+1stack
   //Implement Sec. 2.3
   EvaluateΞ(U(tn+1))using Eq. (2.29)
   PushF(Ξ(U(tn+1)))
   ontoFu+1stack
   Pushfm1(tm)=
   ReformingBondMassFraction
   ontofu+1stack
   Call CullGenerations
Else iftn+1=tm
  Update previously pushed
  U(tm), F(Ξ(U(tm)))andfm1(tm)
FunctionΨr
  Letm=number of generations
  Ifm =0 ReturnΨr=Ψre(F(tn+1))+Ψ0a(F(tn+1))
  LetΨr=Ψre(F(tn+1))
  Foru =0 tom − 1
   Evaluatewu1(tn+1)=fu1(tu)g(F(tu),tn+1tu)
   EvaluateFu1(tn+1)=F(tn+1)·U1(tu1)
   EvaluateΨr+= wu1(tn+1)F(Ξ(U(tn+1)))Ψ0a(Fu1(tn+1))
  ReturnΨr

3 Results

The results corresponding to Fig. 9 of the original paper become updated when employing this new formulation as follows.

We refitted the experimental data to this model to extract the material parameters μe and μa (shear moduli of the strong and weak bonds); τ10, τ20, τ21, and s2 for the reduced relaxation function; and μ and Ξ0 for the recruitment function. This eight-parameter fit produced μe=0.40kPa, μa=1.79kPa,τ10=0.244s,τ20=205s, τ21=2575s,s2=0.0096,μ=30.0, and Ξ0=0.971. The nonlinear regression coefficient between the fit and data was R2=0.991, and the fit is shown in Fig. 9(a). To better understand the strain-dependence of τ2(K2) and F(K2) for this fit, these functions are shown in Fig. 9(b), calculated using τ20, τ21, and s2 for τ2(K2) and μ, Ξ0 and α for F(K2).

The most significant difference we found in the curve-fit of the experimental data of Fig. 9, when compared to the original publication, was the observation of a smoother decrease with increasing strain, in the relaxation time τ2 of the Malkin-distortion reduced relaxation function (Fig. 9(b)).

Acknowledgment

The authors would like to thank Dr. Lance Frazer of the Southwest Research Institute, San Antonio, Texas, for identifying the nonphysical response when modeling cyclical loading of a reactive viscoelastic material with weak-bond recruitment.

Funding Data

  • National Institute of General Medical Sciences (Award ID: R01 GM083925; Funder ID: 10.13039/100000057).

Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

Disclaimer

The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

References

1.
Ateshian
,
G. A.
,
Petersen
,
C. A.
,
Maas
,
S. A.
, and
Weiss
,
J. A.
,
2023
, “
A Numerical Scheme for Anisotropic Reactive Nonlinear Viscoelasticity
,”
ASME J. Biomech. Eng.
,
145
(
1
), p.
011004
.10.1115/1.4054983
2.
Ateshian
,
G. A.
,
2015
, “
Viscoelasticity Using Reactive Constrained Solid Mixtures
,”
J. Biomech.
,
48
(
6
), pp.
941
947
.10.1016/j.jbiomech.2015.02.019
3.
Maas
,
S. A.
,
Ellis
,
B. J.
,
Ateshian
,
G. A.
, and
Weiss
,
J. A.
,
2012
, “
FEBio: Finite Elements for Biomechanics
,”
ASME J. Biomech. Eng.
,
134
(
1
), p.
011005
.10.1115/1.4005694