Stop using Prony series - Try KWW model

1. Introduction

To start with, I’m not saying that using the Prony series is wrong. In some cases, there's simply a better model. Currently, commercial finite element software models linear viscoelastic materials using the Generalized Maxwell model. This approach connects springs and dashpots in parallel, as shown in Fig. 1a. The number of parallel networks is increased until the stress relaxation curve, measured in experiments, is adequately reproduced. This method can be effectively used when the stress relaxation curve from uniaxial tension tests is available.

However, obtaining such curves may be challenging, as in the case of soft tissue, where indentation tests are often used. Similarly, in some cases, only stress-strain curves under different strain rates may be available without a stress relaxation curve. In this cases, Prony series coefficients cannot be obtained directly from stress relaxation curve, requiring inverse FE analysis, where material parameters are adjusted to match experimental load-displacement curves.

There is a limit to how many viscoelastic networks can be added for inverse FE analysis, as each new network adds two additional material parameters, leading to more iteration. For instance, three networks require six parameters. Therefore, a more flexible model with fewer parameters is needed.

The KWW (Kohlrausch-Williams-Watts) model, also known as the Stretched Exponential model, provides such flexibility. The KWW model introduces a power term within the exponential function, offering more adaptable stress relaxation curves, as shown in Fig. 1b.

Since the KWW model only requires three parameters, it is advantageous compared to continuously adding Prony series terms. This post discusses two ways to implement the KWW model in Abaqus: using a user subroutine and the Prony Spectrum approach.

2. Implementation

2.1 UCREEPNETWORK Subroutine

Abaqus allows the implementation of user defined viscoelastic network through the UCREEPNETWORK subroutine. To implement the KWW model, the creep strain increment of the viscoelastic network must be defined based on stress relaxation behavior (Eq. 1). Fig. 2 illustrates an Elastic-Viscoelastic network used to derive the creep strain increment. In Fig. 2, Ee and Ev represent the elastic moduli of the springs in the elastic and viscoelastic networks, respectively, and η is the viscosity of the dashpot. If this value is constant, the network is equivalent to a single-term Prony series. In the KWW model, however, the viscosity is assumed to vary. ε represents strain, with the subscript s denoting the spring and cr representing creep. In the viscoelastic network, the spring and dashpot are subjected to the same stress, which leads to the following equation.

In a stress relaxation test, the total strain of the network remains constant, and can be expressed as follows.

Differentiating both sides with respect to time gives the following

Now, we differentiate Eq. 2 with respect to time.

After rearranging both sides and integrating, the equation is expressed as follows.

Eq. 6 represents the stress relaxation behavior and should be expressed in the form of the KWW model in Eq 1.

The expression inside the exponential function must be the same in the above equation, and using this, we can derive the following equations for the elastic modulus and viscosity of the viscoelastic network

By substituting Eq. 8 into Eq. 2, the creep strain can be expressed as follows.

Here, Ev cannot be directly input into Abaqus. Instead, the stress ratio, r, between the elastic and viscoelastic networks is defined. Expressing Ev in terms of the stress ratio, we get the following.

Eq. 9 is a one-dimensional relation under uniaxial stress. However, since real materials are in a multi-axial stress state, it must be expressed using the effective creep strain and the effective deviatoric Kirchhoff stress. The following equation is substituted into Eq. 9.

Here, J is the determinant of the deformation gradient tensor. The expression for the effective creep strain is given as follows.

In the UCREEPNETWORK subroutine, the effective creep strain increment and the partial derivatives for the Newton iteration are defined as follows.

2.2 Prony spectrum

In polymers, cross-links with various lengths contribute to stress relaxation, with shorter links relaxing faster and longer links slower. The frequency of certain cross-link lengths leads to a higher contribution from corresponding relaxation modes. Prony Spectrum quantifies this, expressing the Prony series coefficients, gi as functions of relaxation times, Ï„i. Here, gi can be thought of as the occurrence frequency of cross-links with a specific length, and Ï„i corresponds to the relaxation time of those cross-links.

The Prony spectrum of the KWW model was derived by Berberan-Santos et al., 2005, as shown below. In this post, we will not cover the derivation.

Here, g0, τ0 and β is KWW model parameter in Eq. 1. In Eq. 18, f(τ) can be regarded as a probability distribution function. Although it appears complex due to its representation as an infinite series, a numerical approximation can be obtained by summing a sufficient number of terms. Once f(τ) is known, gi can be calculated using the following relation.

Fig. 3 shows the probability distribution f(Ï„). The distribution shows a skewed bell-shaped curve.

Using the above equation, the KWW model can be represented as a Prony series, and the coefficients can be entered into the Abaqus input file in table form.

3. Validation

In Section 2, two methods for implementing the KWW model in Abaqus were introduced: the UCREEPNETWORK subroutine and the Prony spectrum approach. To validate these methods, a single-element model was subjected to uniaxial tension, and the displacement was fixed to predict the stress relaxation curve. Fig. 4 shows that the stress relaxation curves predicted by both methods shows good agreement with the analytical curve.

4. Conclusion

In this post, I introduced the KWW model, a linear viscoelastic model capable of simulating highly flexible stress relaxation curves using just three parameters. Two methods for implementing the model in Abaqus were presented: a user subroutine and the Prony spectrum approach. The predicted stress relaxation curves using both methods were also validated. The KWW model can significantly reduce the number of iterations required in the optimization process for material property evaluation through inverse FE method and decrease the data needed when building an training database for AI model. Since there is no open resource in the literature on how to implement the KWW model in Abaqus, I hope this post will be helpful to engineers working in this field. The input files and code used in this post can be downloaded from the Github repository linked below.

https://github.com/YB-LIM/KWW_Viscoelastic

Reference

[1] Abaqus Documentation, 2024, Dassault Systèmes Simulia Corp., Providence, RI, USA.

[2] Berberan-Santos MN, Bodunov EN, Valeur B, 2005, Mathematical functions for the analysis of luminescence decays with underlying distributions 1. Kohlrausch decay function, Chemical Physics 315, 171-182.

[3] Mauro JC, Mauro YZ, 2018, On the prony series representation of stretched exponential relaxation, Physica A 506, 1181-1190.