## Abstract

We propose a mathematical model for relaxation modulus and its numerical solution. The model formula is extended from sigmoidal function considering nonlinear strain hardening. Its physical meaning can be interpreted by a macroscale elastic network-viscous medium model with only five model parameters in a simpler format than the molecular-chain-based polymer models to represent general solid materials. We also developed a finite-element (FE) framework and robust numerical algorithm to implement this model for simulating responses under both static and dynamic loadings. We validated the model through both experimental data and numerical simulations on a variety of materials including asphalt concrete, polymer, spider silk, hydrogel, agar and bone. By satisfying the second law of thermodynamics in the form of Calusius–Duhem inequality, the model is able to simulate creep and sinusoidal deformation as well as energy dissipation. Compared to the Prony series, the widely used model with a large number of model parameters, the proposed model has improved accuracy in fitting experimental data and prediction stability outside of the experimental range with competitive numerical stability and computation speed. We also present simulation results of nonlinear stress–strain relationships of spider silk and hydrogels, and dynamic responses of a multilayer structure.

## 1. Introduction

Almost all materials exhibit viscoelastic behaviour more or less. Viscoelastic material exhibits both elastic and viscous characteristics, and its deformation is temperature dependent. Relaxation modulus *E*(*t*) is a characteristic of material viscoelasticity as used to describe the stress relaxation of materials with time (*t*). It is important to accurately simulate the stress relaxation and viscoelastic deformation of subjects in order to evaluate and design materials. Thermal transitions of viscoelastic materials can be described either in terms of free volume changes or relaxation time [1]. The Crankshaft model [2] simulates polymer molecules as a series of jointed segments involving in a few stages of thermal transition with elevated time or temperature. For example, with the elevation of temperature, the first or *γ* transition may start with the local motion of molecules, and then the *β* transition may appear when *E*(*t*) drops slightly due to the bend and stretch of molecules with elevated temperatures. Consequently, the glass or *α* transition occurs when *E*(*t*) significantly reduces until reaching the rubbery stage, and lastly, the terminal transition may exit when the polymer melts into liquids due to the slippage of molecular chains. Different materials may present different thermal transition behaviours, but generally they all include the glass transition at the solid state, which is also the focus of this study. Different theories and physical models for describing *E*(*t*) were proposed primarily in the field of polymer science. The majority of these models focused on linear viscoelasticity. Among these, the Rouse model simulated the single chain diffusion of polymers as Brownian motion of a beads–harmonic-springs system based on the theory of molecular dynamics (figure 1*a*) [3], and the Kremer–Grest model used hundreds of chains and beads in its simulation framework [4] The well-known (entangled) tube model has defined entangled polymer chains that are confined in tubes with permanent topological interactions and move along tubes (figure 1*b*) [5]. This model especially concerned the phenomena associated with polymers of complex topology or long chain branching [6]. The stress relaxation of each chain is calculated as the fraction of the tube that has not been vacated, where the relaxation time is related to the molecular mass of the cube [7]. The arm-retraction model [8] described the entangled monomers (unit of polymer macromolecular) as retracted by arms (figure 1*c*), which was validated for the rheology of entangled polymer liquids.

Maxwell model proposed by Maxwell and Wiechert in late nineteenth century consists of a linear spring and a dashpot in series. With an additional spring in parallel, it is known as the standard linear solid or Zener model. For the dashpot itself stress has a linear relationship with strain rate as *η* is viscosity. However, the Maxwell model can only capture the relaxation behaviour of polymers in a very limited time range. Accordingly, a series of such models have been assembled in parallel, delivering a physical system referred to as the generalized Maxwell (GM) model (figure 1*d*) to improve accuracy in fitting experimental data. The combination of elastic springs and viscous dashpots in this physical system could be interpreted by variable molecular-chain-segment lengths under different time distributions for polymer structures [9]. The GM model expresses *E*(*t*) as follows:
*E*_{∞} is the model's elastic modulus at *t* = ∞, *E _{i}* is the elastic modulus of the spring and

*η*is the viscosity of the dashpot for the

_{i}*i*series, respectively,

^{th}*n*is the number of spring–dashpot series, and

*E*to 0 for a single spring–dashpot series.

_{i}The *E*(*t*) value of this physical system is the sum of *E _{i}* distributed at variable time spectrums. This formula with accumulated exponential terms in a discrete spectrum with a finite number of series is also referred to as Prony series (PS). The model parameters are determined by fitting on experimental data using mathematical optimization skills. Usually a higher

*n*value may achieve a higher fitting accuracy, but also involves more complexity and variability as will be discussed later. The GM model or the PS has been widely used to describe the linear viscoelastic behaviour of a broad range of solid materials. Examples of these materials include polymers [10], waterproof membrane [11], glasses [12], poroelastic medium [13], asphalt concrete (AC) [14], ligament [15] and skin [16], etc. The numerical implementation of the PS is computationally convenient due to its exponential format for a closed-form solution of time integration [17]. It has been implemented as a standard viscoelastic material input in the popular finite-element (FE) software (e.g. ANSYS and ABAQUS). Dynamic modulus

*E**(

*ω*) is also used to describe the viscoelastic behaviour based on frequency domain [18], which is not covered in this paper.

Many other viscoelastic models have also been proposed to simulate specific materials. Iyo and Sasaki *et al*. proposed a Kohlraush–Williams–Watts function based model for biomaterials as follows [19]:
*τ*_{1} and *τ*_{2} are characteristic time of relaxation process, and ϑ is a model parameter.

The fractional derivative model was proposed to express the stress–strain relationship of viscosity as a fractional order of time as follows [20]:
*θ* is a non-integer parameter. The stress–strain relationship of the standard viscoelastic system using a fractional derivative order can be expressed as follows:
*D* is a fractional derivative, and *E*_{0} and *E*_{1} are two modulus parameters. The fractional derivative model can reduce the number of model parameters while attaining a desired accuracy. It is especially appropriate for predicting the dynamic behaviour of polymers as the non-integer parameter, *θ*, offers nonlinearity and flexibility to capture a wide range of relaxation time. In addition, it may simulate the laboratory-measured asymmetric peak of loss modulus by using different fractional terms for stress and strain, respectively [21]. However, the fractional derivative model involves more complexity in fitting experimental data and numerical implementation. The Huet–Sayegh model expresses *E**(*ω*) in fractional order [22], but it is not designed to simulate stress relaxation with time domain [23]. Nonlinearity of viscoelastic models have been considered for solids with large deformation or properties that change with deformation [24]. Schapery's model considers the spring's elastic modulus as a nonlinear function of time [25]. Other models instead present the dashpot's stress as a nonlinear function of strain rate such as following the power-law [26] or exponential function [27]. However, computation instability may arise for nonlinear models [13,28].

Based on literature review, a few research questions may arise for existing models and approaches. The (multi-molecular) chain-based models and theories discussed above (e.g. Rouse model, tube model and arm-retraction model) were primarily derived from and validated for monomers and polymers. These models may also involve difficulty in characterizing model parameters and numerical solutions at macroscale. More importantly, the widely used PS and its modified ones can raise a few research questions. Firstly, the PS only describes linear viscoelasticity within a small strain range. Secondly, it is known that the PS produces unsatisfied curve fitting on experimental data [23,29] due to its discrete spectrum as a finite set of exponential functions. This is especially true when a relatively small number of terms is used as a common practice (e.g. *n* = 7 with 15 model parameters or even less [30,31]). With a larger *n* number data fitting may get more accurate. However, it becomes difficult to fit a large number of model parameters mathematically with multiple solutions because the optimization intrinsically involves more variability [29]. It also becomes difficult to provide a clear physical interpretation of this extensive spring–dashpot system, although it can improve data fitting accuracy. In addition, the PS is primarily used for fitting experimental data rather than predicting *E*(*t*) outside of the experimental range, while prediction may be warranted when experimental tests have limitations to capture *E*(*t*) at the low and high relaxation time ranges. Lastly, strain hardening or softening may appear in some special materials [32]. This hardening or softening feature has not been acknowledged for its potential effect in the PS and many other existing (nonlinear) viscoelastic models [25–27].

Accordingly, the objective of this research is to develop an alternative material model in a relatively simple format with fewer model parameters for describing *E*(*t*) and simulating viscoelastic responses of solid materials to improve accuracy and efficiency. We extended our study from an early version published on the arXiv website [33] by conducting more comprehensive theoretical and experimental validations, including model derivation, detailed experimental validations and simulations of the stress–strain relationships of broad materials. We tested the model through experimental data on a variety of materials including asphalt concrete (AC), shape-memory polymer, hydrogels, spider silk, agar and bone. We also developed a FE scheme to implement this model for simulating responses. The completion of the model also included a check for numerical accuracy, stability and convergence of the proposed model versus that of the PS.

## 2. Model development

### (a) Model formula

Before deriving the model, we reviewed a few continuous-time-spectrum based mathematical models. Nutting [34] proposed a general stress–strain relationship as relating to a power function of time that *ε* = *σt ^{A}*.

*E*(

*t*) derived as the reciprocal of a power function is called the power law model [35]:

*δ*are model parameters. The power law model has the advantage of simplicity by using only two model parameters. However, it is unable to capture

*E*(

*t*) at the high and low time ranges (thresholds) accurately [23]. The S-shape sigmoidal function may capture these two thresholds more accurately as:

By including additional model parameters into the sigmoidal function, the modified model can more accurately control the model shape and trend. For example, the modified model is proposed to fit the absolute value of dynamic modulus for asphalt materials [36,37]:
*δ*, *α*, *β*, and *γ* are fitted model parameters, and *f* = *ω*/2*π* is frequency. The modified sigmoidal function model has been successfully implemented for fitting experimental data of materials including bitumen and asphalt. However, as a mathematical function it does not describe material viscoelasticity based on a continuum-mechanics framework. For example, even the unit of the left side of equation (2.3) does not match the right side physically. As a result, it has not been used to simulate viscoelastic responses of structures other than mathematical fitting on experimental data.

By extending the continuous-time-spectrum-based models (i.e. the power-law and sigmoidal function), we introduced a new mathematical model with physical mechanism considered to describe the material viscoelasticity at both experimental and numerical scales, for improving accuracy of existing models. With validation of the model formula on experimental data of a variety of materials, we proposed the model as follows, where the log–log scale is applied for improving accuracy:
*E*_{∞} is the minimum modulus at the rubbery stage, *E*_{0} is the instantaneous modulus at the glassy stage, *μ* and *α* are unitless model parameters, log = log_{10} for our implementation and *τ*_{0} is a characteristic time of relaxation process (in seconds) as related to the time after loading required for the strain to reach its equilibrium, and the value. *t*/*τ*_{0} becomes unitless such that the left and right part of the equation equals physically. One shall note that in this model the exponential term *t*/*τ*_{0})* ^{α}*. However, the logarithmic scale expression (the same as equation (2.3)) is adopted here for the purposes of: (1) using logarithmical expression for both

*E*(

*t*) and

*t*to be continuous; (2) the logarithmical scale of the model improves accuracy in fitting experimental data; and (3) attaining computation convenience for numerical implementation with the exponential format.

In addition, some viscoelastic materials may exhibit strain hardening behaviour [38] that represents the strengthening of materials with strain [39]. Accordingly, the proposed model can be extended to consider the nonlinear strain hardening or softening, in which *E*(*t*) decreases with strain for strain softening as opposite to that of strain hardening, as follows:
*β* is a strain hardening or softening coefficient. By setting *τ*_{0} as a constant value such as 1 (s) for the simplicity purpose, the model parameters are reduced to five. To consider the temperature effect, a shift factor *α _{T}* is included to convert time to a reduced time

*t*at the reference temperature according to the temperature–time superposition principle as:

_{r}Here we adopted Arrhenius law to describe the temperature–time superposition [40] for *α _{T}* as:

*T*is temperature (K),

*T*

_{0}is the reference temperature (K),

*E*is a fitted model parameter that may be referred as an ‘energy’ (J Mol

_{a}^{−1}) required for thermal transition,

*R*is a constant (like the gas constant in J K

^{−1}Mol

^{−1}), and

*E*/

_{a}*R*is a model parameter index (K). Other formulae for

*α*have also been proposed including the WLF (Williams–Landel–Ferry) model [41] and the polynomial function [37]. The WLF equation might be a better description for temperatures above the glass transition, while the Arrhenius equation could be more suitable at temperatures below glass transition [42]. The formula of

_{T}*α*shall be properly chosen as dependent on the specific material and temperature range. In this study, we have not evaluated the effect of

_{T}*α*formula.

_{T}### (b) Physical interpretation

Network model-based physical mechanisms have been often used for polymers, rubber and elastomers, although the dynamics of polymer networks is complex and not fully understood yet in polymer physics. Examples include the coarse-grain polymer network model, in which the mesh-like networks are stretched in the viscous medium [43], the symmetrically growing treelike micro-network model built on sub-chains and the non-affine microspheres [44]. However, molecular and micromechanical models are more appropriate for specific type of materials, e.g. the multi-molecular-chain-network models specifically for polymers as discussed above. We studied a few different materials as plotted in figure 2: (1) agar material (e.g. agar cake, figure 2*a*), spider silk (figure 2*b*) and polymer which belong to the multi-chain cross-linked network (figure 2*d*). Spider silk is a protein fibre consisting of short polypeptide stretches of about 10 to 50 amino acids. These patterns can be repeated more than 100 times within one individual protein, and each polypeptide duplication results in extremely durable spider silk threads [45]; and (2) other materials that cannot be represented by the multi-chain network including organic asphalt (figure 2*c*) and bone (figure 2*f*). Bone, a rigid and lightweight composite, primarily consists of hydroxyapatite-like mineral particles (e.g. calcium) embedded in a matrix of collagen fibres. The cortical or compact bone, which is the dense outside layer has noticeable viscoelasticity as contributed by the collagen fibres and non-fibrous proteins in the bone matrix [46–48]. Asphalt is heavy organic presented in petroleum and performs as a viscous fluid at high temperature and viscoelastic solid at normal or low temperature. Asphalt is amorphous with complex molecular structure and low molecular weight. Organic asphalt cannot be separated into individual components or narrow fractions [49]. Therefore, its viscoelastic behaviour may not be properly interpreted by the molecular-chain-based models, which is also true for many other materials with different morphologies than polymers and rubber-like materials.

We interpret the physical meaning of the proposed model to represent general materials with various morphologies based on a relatively simple and network-based model. As illustrated in figure 2*e*, the model consists of an elastic network with viscous medium filled between the network meshes. The elastic network carries *E*(*t*) ∈ [*E*_{∞}, *E*_{0}], which transits between the glassy stage and the rubber stage (*T* − *t*) superposition principle, the temperature effect is equivalent to the reduced relaxation time *t _{r}* (i.e. a higher

*T*corresponds to a higher

*t*and vice versa). At higher

_{r}*T*or

*t*the model domain becomes ‘softer’ with lower

_{r}*E*(

*T*) value until reaching

*E*

_{∞}at the infinite time (rubbery stage). At the zero time or infinitely low

*T*, the network carries the highest modulus of

*E*

_{0}(glassy stage). The

*β*coefficient characterizes the strain hardening or softening behaviour. The

*α*coefficient controls the slope or rate of

*E*(

*t*) change such that ∂

*E*(

*t*)/∂

*t*∝

*α*. A higher

*α*value indicates a higher thermal sensitivity, resulting in more

*E*(

*t*) reduction (figure 2

*g*). The

*μ*coefficient is related to friction, and a higher 1/

*μ*value indicates a higher friction force and viscosity leading to a higher

*E*(

*t*) value (figure 2

*h*). The mathematical model formula (on the non-logarithmical scale) may also be represented by a three-element of nonlinear spring–dashpot system with the constitutive stress–strain relation satisfying the Calusius–Duhem inequality [50] for thermodynamic consistency as proofed in the electronic supplementary material a). By extending the model to a three-dimensional (3D) domain, the shear and bulk relaxation moduli

*G*(

*t*) and

*K*(

*t*) can be attained using the same model formula as that of

*E*(

*t*) [51]:

*G*

_{∞}is shear modulus at

*t*= ∞, and

*G*

_{0}is the instantaneous shear modulus,

*K*

_{∞}is bulk modulus at

*t*= ∞,

*K*

_{0}is the instantaneous bulk modulus and

*α*,

_{G}*α*,

_{K}*β*,

_{G}*β*,

_{K}*μ*and

_{G}*μ*are model parameters. When compared to the modified sigmoidal-function model in equation (2.3), the advantages of the proposed model include: (1) it includes a physical mechanism to interpret material viscoelastic behaviours, (2) we developed a FE scheme to implement the model for simulating responses based on a continuum-mechanics perspective; and (3) it includes the nonlinear strain hardening to model broader materials including spider silk. Other referred models [43,44] were designed primarily for polymers at microscale. In comparison, our model is designed to model a broad range of general materials for experimental data fit and prediction, and numerical simulation at macroscale. However, we have not evaluated the model for multiple relaxation processes.

_{K}## 3. Experimental validation

For some materials like AC the *E**(*ω*) is more often measured than *E*(*t*) for practical purpose, and *E*(*t*) can be derived from *E**(*ω*) via the interconversion. However, for other materials like polymer and biomaterials *E*(*t*) are more often measured and then fitted by the PS. In addition, it is important to predict *E*(*t*) values outside of the experimental range to evaluate material property and simulate responses for a wider range of relaxation time. Therefore, it is useful to compare the proposed model and the PS in both fitting and predicting the experimental data of *E*(*t*). When part of the experimental data within a middle time range is used as training data to fit model parameters, the rest of the experimental data at the lower and higher time ranges is predicted using the fitted model parameters. We used the nonlinear reduced gradient method for both the proposed model and the PS to determine model parameters as illustrated in electronic supplementary material b). In addition to the physical-ground based evaluation, we used the reduced *χ*^{2} statistics to quantify the goodness of fit [52] as follows:
*L* is degree of freedom as *L* = *m* − *N* − 1 (*m* is the total number of data points, *N* is the number of model parameters). A lower

In figure 3 we evaluate the fitting accuracy of the PS versus the proposed model for experimental data of *E*(*t*) (converted from the measured *E**(*ω*)). *E*(*t*) data are fitted by the PS with variable term *n* ranging from 1 to 30 (for a total of 3 to 61 model parameters). Results have shown that with a higher *n* the PS intends to capture *E*(*t*) at a wider range of time more accurately. However, as *n* is relatively large (i.e. *n* ≥ 14) the *n*=14 and 30 attain almost identical *E*(*t*) values). Different seed values may result in different *E*(*t*) shapes although with similar *E*(*t*) values at the low time range using relatively large seed values of *n* = 17 are significantly higher than that fitted using smaller seed values for *n* = 14, but both have almost the same *E*_{0} and *E*_{∞} especially when *n* is relatively small (figure 3). With a larger *n* the *E*(*t*) curve may become smoother when using proper seed values. However, it becomes more difficult to properly estimate plenty of model parameters without unique solution. In comparison, the proposed model can achieve a smoother, more accurate and unique curve fitting using much less model parameters.

In the following we present important tests to validate the model on another three materials: polymer, agar and bone. For the PS, a term number of *n* = 14 is used for all materials in the following sections with satisfied accuracy (the model formula and data fitting/prediction are illustrated in electronic supplementary material c) Excel spreadsheet). Figure 4 presents the model fitting and prediction of a polyurethane (PU) polymer material. *E*(*t*) was measured at six temperatures: −30, −22.5, −20, −17.5, −15 and −12.5°C (data from [53]) and ‘shifted’ to that at the reference temperature of −17.5°C through the temperature–time superposition. In comparison, for the PS, the fitted *E*(*t*) curve has a ‘sudden’ jump at the low time range of 10^{−5} ≤ *t* ≤ 10^{−4} s outside of the experimental data range (figure 4). It intends to converge to the glassy and rubbery stages sharply once it is out of the experimental data range. However, it still achieves a high fitting accuracy (for the experimental data range) with a *E*(*t*) at the high time range. This result may imply that the PS has lower stability and higher dependency on seed values than the proposed model especially for the prediction outside of the experimental data range. Figure 5 presents both the fitted and predicted *E*(*t*) of agar materials. The PS and proposed model have relatively close values, which may be due to the relatively high linear relationship of *E*(*t*) versus *t* for this specific experimental dataset. The PS has slightly overpredicted *E*(*t*) at the high time range. Figure 6 presents fitted and prediction of *E*(*t*) of the cortical bone material. We used two groups of seed values for fitting the PS parameters, and the one with higher seed values (*b* versus 6*a*) as indicated by their *E*(*t*) at the high time range as it converged to *E*_{∞} sharply outside of the fitting data range. In comparison, the proposed model produced a smoother and more accurate prediction, and it predicted *E*_{0} as 16.5 GPa which falls within the range of 11–21 GPa as measured by other researchers [55].

Table 1 lists the fitted model parameters for different materials. Bone is the stiffest material with the highest *E*_{0} and *E*_{∞} values and poses the longest relaxation time spectrum, while agar is the softest one with the lowest modulus among these materials. Agar poses a higher *α* value than the others, indicating its relatively higher thermal sensitivity for viscosity, followed by the PU and then the AC, and lastly the bone material which has the elasticity dominant.

The norm-based method was often used to evaluate the accuracy of model prediction, and a lower norm value indicates higher prediction accuracy [56,57]. We use the *L*^{2} norm to evaluate the model prediction as *E _{i}* is the predicted modulus and

*L*

^{2}values of the PS/proposed models are 0.42/0.24, 4.20/0.89 and 0.69/0.15 for agar, PU and bone, respectively, indicating higher prediction accuracy of the proposed model with lower

*L*

^{2}value. The improved prediction accuracy and stability of the proposed model may be explained as follows. The proposed model captures

*E*(

*t*) in a smooth

*S*curve constrained within the [

*E*

_{∞},

*E*

_{0}] range, with the variation slope and trend controlled by

*α*—the thermal sensitivity parameter and

*u*—the friction coefficient. As a result, the model predicts

*E*(

*t*) growth outside of the experimental range naturally to capture the glass transition. In comparison, the PS intends to attain mostly accurate fitting only within the time range of the available data and it may sharply converge to the more likely false

*E*

_{0}and

*E*

_{∞}values once it is outside the data range (figures 4–6). Given

*E*

_{0}and

*E*

_{∞}values as constraints (then it is not a true prediction), the PS may achieves higher prediction accuracy, but a sharp converge may still be seen due to its discrete finite set of time spectrums. Nevertheless, this comparison is fair because both models use the same mathematical optimization scheme and constraints for determining model parameters. One may argue that given a very large term number for the PS (e.g.

*n*> 30) the fitting and prediction accuracy may be improved, and the pre-smooth method may also be used to improve the accuracy of the PS fitting [29]. However, the proposed model could be superior to the PS for experimental data fitting and prediction regarding the following aspects: (1) a model with less parameters and satisfied accuracy is always preferred for simplicity; (2) a large number of model parameters of the PS produces more instability for data fitting and prediction; (3) the large spring–dashpot systems of the PS may represent multi-chains of polymer and other macromolecular materials physically, but not for many others with different morphology.

## 4. Numerical solution method

### (a) Finite-element solution

After the theoretical demonstration and experimental validation, we implemented the proposed model for numerical simulation under both static and dynamic loadings. The PS has the closed-form solution for the time integration, while for the proposed model numerical discretization is necessary for both the time and space domains. Different numerical methods with different advantage(s) may exist for implementing the same model including the one based on frequency-domain [58]. The objective of our numerical solution is to achieve competitive computation speed and numerical stability when compared to the PS. We proposed a FE solution based on time-domain to integrate *E*(*t*) without domain transformation, which is fast and more straightforward than the frequency-domain method. We proposed a combined time discretization scheme with the Houbolt, forward and backward finite difference methods for shortening time-step length, as extended from a previous work for the PS model [59]. We also conducted error analysis to compare the simulation accuracy of the proposed model versus the PS.

For solids in a continuum state, the strong form of the governing state equation for the dynamic problem can be formed as follows:
*σ* is stress tensor, *u* is displacement, *b* is body force, *ρ* is material density, *t*_{0}, *t _{d}*] is a time domain. The stress–displacement relationship of the viscoelastic solid can be expressed as follows [59]:

*p*(

*t*) on both sides of the strong form (equation (4.1)), and then integrate with the space and time domains to form a weak form as follows:

According to divergence theory, the divergence of stress can be decomposed into two parts:
*f*(*t*). Substitute equation (4.4) into equation (4.3), and the weak form can be re-expressed as follows:

Substitute equation (4.2) into equation (4.5) and then re-arrange to reach the following weak form:

*R*_{0} (*t* − *τ*) is a relaxation modulus matrix for linear viscoelasticity in the 3D domain, which is defined as follows [59]:
*u*(*t*) and *p*(*t*) are then discretized to that at the FE nodes and presented in a vector format through the shape function. The discretized *p*(*t*) value is arbitrary and can be dismissed on both sides of the equation. Thus, the following weak form shall satisfy

The time domain *t* ∈ [*t*_{0}, *t _{d}*] is discretized to

*N*time steps for

*k*= 1, 2…

*N*, and for each of them the sub-time domain

*τ*∈ [0,

*t*] includes

*k*time steps for

*j*= 1, 2 …

*k*. Thus, equation (4.10) can be rewritten as follows after discretizing

*τ*into

*k*sub-time steps:

*t*and

_{j}*t*

_{j}_{−1}are the

*j*and (

^{th}*j*− 1)

^{th}time.

The gradient term *τ* ∈ [*t _{j}*

_{−1},

*t*] can be discretized

_{j}*τ*∈ [

*t*

_{j}_{−1},

*t*], and thus equation (4.11) is re-expressed as follows:

_{j}*u*(

*j*) is displacement at the

*j*

^{th}sub-time step for

*j*= 1, 2, 3 …

*k*, and

*J*(

*j*) is a viscoelastic stiffness matrix defined as follows:

Poisson's ratio can be time dependent for different test modalities [60]. For viscoelastic materials the temperature effect can be converted to a relaxation time through the temperature–time superposition rule. Given a constant temperature and strain rate within a short loading time period, it can be reasonable to consider Poisson's ratio as constant for the simplicity purpose in numerical implementation. Therefore, *G*(*t*) and *K*(*t*) values can be directly derived from *E*(*t*) following the relationship among the elastic Young's, shear and bulk moduli [59]. Therefore, *J*(*t*) can be rederived as follows:
_{e} is a discretized elastic matrix at FE nodes of the fourth-order elasticity tensor. For the proposed model, substitute equation (2.4) into equation (4.14), the viscoelastic stiffness matrix *J _{P}*(

*j*) can be calculated as follows:

Following the trapezoidal rule, *J _{P}*(

*j*) can be discretized as:

*t*=

*t*−

_{j}*t*

_{j}_{−1}is the time-step length. While for the PS, substitute equation (1.1) into equation (4.14),

*J*(

*j*) can be attained as follows:

The Houbolt method is adopted for time discretization of acceleration due to its low time-step dependency and high stability [61]:
*u*(*k*) is displacement at the *k*^{th} time step for *k* = 1, 2, 3 … *n*. A relatively great time-step length (i.e. 0.001 second) was used to improve computation speed satisfying numerical accuracy. For solving the nonlinear viscoelastic displacement, substitute equation (4.18) into equation (4.12) to attain the discretized:

We employed the Cholesky factorization with the iterative refinement method to solve the global linear system. At each iteration step Newton's method can be applied to solve *u*(*k*) by calculating the tangent of the function *t* = 0 *u*(0) = 0 and

At *t* = 0 the governing state equation satisfies:

Given *ε* (0) = 0, *u*( − 1) and *u*(1) can be solved from equation (4.20) and (4.21), then *u*(2) is solved following the Houbolt time discretization, and so on for *u*(3), *u*(4) … *u*(*N*). With displacements calculated at all FE nodes, strain can be solved as

We developed the VBA/FORTRAN coding to implement the numerical solution (see computer code in electronic supplementary material d) for one-dimensional stress solution and the subroutine of stiffness matrix for two-dimensional solution).

### (b) Error analysis

In this section, we analyse the simulation error for the proposed model versus the PS, in which strain hardening or softening is not considered. The major numerical error of the viscoelastic solution lies in the computation of the viscoelastic stiffness matrix *J*(*j*). Using the trapezoidal rule [63] the truncation error of the proposed model for the discretization of *N* is the number of time grids. Substitute equation (4.16) into equation (4.23):

The first derivative of *E*(*t* − *ξ*) with respect to *ξ* can be derived as:
*A* = *μ ^{α}*, and

*J*(

*j*) has a closed-form solution without needing time discretization. However, the PS produces computational errors due to its less accurate data fitting, for which the computation error can be calculated as follows:

*E*

_{ps}is the PS fitted modulus,

*E*is the true value and C

_{e}is elastic tensor.

Here we present the error analysis results of agar material as a numerical example because the *E*(*t*) values of agar fitted by the PS are closer to that of the proposed model than other materials (figure 5). Figure 7 plots the *J*(*t*) errors of the proposed model versus that of the PS during 10 s. Results indicate that the *J*(*t*) errors of the proposed model are negligible (i.e. average 0.016% and max 0.46%) and smaller than that of the PS.

## 5. Numerical analysis and implementation

With the numerical algorithm implemented and laboratory-measured material property as inputs, we simulate the displacement, stress–strain relationship and hysteresis considering strain hardening, and energy dissipation of a variety of materials. We evaluate the computation accuracy, speed and stability for convergence of the proposed model versus the PS.

### (a) Creep and cyclic deformations

In order to validate the accuracy of the developed FE algorithm, we first compared its simulation results with that of ANSYS—one of the most powerful commercial FE software using the PS model. They attain almost identical results although they used different algorithms (e.g. different time discretization methods). We simulated a unit agar body as illustration to evaluate the numerical stability and convergence. Figure 8*a* presents the simulated deformations under a constant Heaviside loading (1 kPa) for *t* ∈ [0, 10 s]. Note that the initial deformation at zero time is the instantaneous elastic responses. Different time-step lengths *dt* ∈ [0.01, 1 s] are used, and both models converge at *dt* = 0.01 s although the PS is slightly less sensitive to *dt*. Figure 8*b* presents simulated deformations under the sinusoidal loading *f*(*t*) = 0.5sin(*πt* − *π*/2) + 0.5 with an amplitude of 1 kPa. Again simulation results show that both models converge at the same *dt* = 0.1 s. The constant loading requires a shorter *dt* for convergence than the cyclic loading due to the resulting larger creep deformation. Figure 8*c* compares the simulated deformations of the two models with *dt* = 0.01 s. The difference is as high as 4.2 and 8.1% for the constant Heaviside loading and sinusoidal loading, respectively. This discrepancy is primarily introduced by the PS's less accurate *E*(*t*) model fitting on experimental data.

### (b) Stress–strain relationship of materials

We simulated stress–strain curves of AC, PU polymer, agar and bone as plotted in figure 9*a–d*, clearly showing material viscoelastic behaviours. Note that different strain ranges are used for different materials due to their different viscoelastic ranges [64–67]. We also simulated the nonlinear stress–strain relationship of different materials. In the laboratory test, a spider silk thread was glued to a glass substrate with the glycoprotein glue (figure 9*e*) and was subjected to a stretching loading (data reproduced from [68] with permission). Results have shown that the model is able to capture the nonlinear strain hardening behaviour of spider silk (figure 9*e*), and the nonlinear strain softening and then hardening of the polyampholyte-formed hydrogel material (figure 9*f*, laboratory data reproduced from [68] with permission). This hydrogel performed viscoelastic behaviour with stress recovery for strains up to 700% [68]. These materials may perform plastic deformation and fracture at larger strain levels which are studied by our model. Finally, we extended the model to simulate the dynamic responses of a multilayered pavement structure as illustrated in the Supplementary material e).

## 6. Conclusion

We proposed a mathematical model for the data fitting and prediction of *E*(*t*) and numerical solution of viscoelastic responses for a broad range of materials. We interpret its physical mechanism as an elastic network–viscous medium system with five model parameters to represent general materials considering nonlinear strain hardening. We then developed a Galerkin FE method and robust numerical algorithm to implement the model for simulating dynamic responses. We have validated the model on both experimental data and numerical simulation for a variety of materials including asphalt concrete, polymer, spider silk, agar and bone. The model has shown some advantages when compared to the PS—the widely used model including: (a) it improves the fitting and prediction accuracy for experimental data while using much less model parameters; (b) its numerical solution has competitive computation speed and numerical stability for convergence; and (c) it has considered the nonlinear strain hardening behaviour of some special materials. The model is able to simulate the creep and sinusoidal responses reasonably under both static and dynamic loadings. Therefore, the model may serve as an alternative to simulate the viscoelastic behaviours of solids at both experimental and numerical scales with improving accuracy while reducing complexity.

## Data accessibility

Part of the data along with the model formula is presented in the electronic supplementary materials c and d, and the other data can be attained by replacing the model parameters.

## Authors’ contributions

Q.X. developed the model, performed analysis and wrote the paper. B.E. advised the model development and reviewed the paper with revision advice.

## Competing interests

We have no competing interests.

## Funding

Part of the data collection was supported by the Federal Highway Administration project DTFH61-07-C-0032.

## Acknowledgements

We appreciate the courtesy of those researchers or publishers [46,53,55,63,68] for offering part of experimental data.

## Footnotes

Electronic supplementary material is available online at https://doi.org/10.6084/m9.figshare.c.4088969.

- Received August 11, 2017.
- Accepted March 6, 2018.

- © 2018 The Author(s)

Published by the Royal Society. All rights reserved.