Applying the fractional derivative Zener model to fitting the time‐dependent material viscoelasticity tested by nanoindentation

Abstract

Material viscoelasticity is of particular importance for polymers, biotissues and many others, due to the significant influence not only on the elastic response of structures, but also on their failure mechanism in actual applications. In this study, the experimental tests by nanoindentations are categorized according to the objective parameter combined with the indenter shape into four distinct cases, and the Fractional Derivative Zener model (FDZM) is employed to characterize this time‐dependent property. The explicit data‐fitting formulas are derived to determine the relaxation modulus and creep compliance tested with the conical and spherical indenters. Besides, a fitting scheme is provided to facilitate the identification of these viscoelastic parameters. In comparison with the general methods adopting the Prony series or empirical models, two silent merits present: (1) fewer model coefficients necessitate being determined based on the experimental data; (2) the relaxation modulus G(t) and the creep compliance J(t) can be obtained through one indentation test, and their exact interconversion is guaranteed as well. Several sets of experimental data are applied to examine the validity of the developed approach, and three more examples are adopted to demonstrate the exact interconversion between G(t) and J(t). A few concluding remarks are drawn eventually.

1

Introduction

Material viscoelasticity has been a research focus of wide interest since Boltzmann [] first established the linear governing formulation, and is of particular importance not only for the homogeneous/heterogeneous materials themselves but also for the structures fabricated ascribed to its time-dependent essence. Two practical examples are the polymer bonded explosives for rocket propellants and the bionic materials for tooth repairs. In this study, our focus is mainly on determining the viscoelastic model coefficients accurately according to the experimental data by nanoindentation and achieving the exact interconversion of the relevant elastic moduli.

Different from the general elastomers whose constitutive formula and elastic moduli are time independent, structures and/or materials involving viscoelasticity are characterised by the time-related governing formulations. The two commonly employed quantities are the relaxation modulus G(t) and the creep compliance J(t). By far, many conventional methods and instruments [, ] have been developed to determine these moduli experimentally. However, these measurements encounter severe challenges due to the limited size of the experimental samples when applied to the nanocomposites, hierarchical biomaterials and microelectromechanical devices. In recent decades, the nanoindentation technique developed based on the theory in [] has become a popular and effective method to determinate the viscoelastic moduli of such objects []. Lu et al. [] proposed and validated a nanoindentation technique to measure the creep function of linearly viscoelastic materials of very small amounts, and further established a differentiation as well as a curve-fitting method to measure the Young's relaxation function of rigid conical indenters in []. Hereafter, Oyen [] applied the elastic–viscoelastic correspondence analysis to generate solutions for indentation creep under the fixed loadings, and presented a simple experimental technique for analysing the viscoelastic behaviours tested by the load-controlled spherical indenter. Peng et al. [] further established a revised step-hold method by replacing the measured load–depth curve in [] with the revised load–depth curve to determine the shear creep compliance of linear viscoelastic–plastic solids. To be noted, most of these works applied the Prony series to charactering the material viscoelasticity in view of its simplicity and convenience for mathematical operations. However, an important shortcoming (i.e. the presence of negative physically meaningless model coefficients) of these methods should be checked carefully and addressed artificially when a wide scatter and unsmoothed set of experimental data intervenes. Commonly, interactive adjustments of the initial fitting parameters or other optimisation techniques are adopted to cope with this issue. Some other deficiencies of the methods adopting the Prony series include: (i) few data-fitting formulas have been proposed by far for spherical indenters [, ]; (ii) the exact interconversion between G(t) and J(t) can be hardly achieved when the number of the exponential terms in Prony series is larger than 3.

However, in practice, only one set of these two modulus functions may be tested due to the instruments available or only the stress-controlled creep test can be conducted due to the stiffness of the viscoelastic materials. Under these situations, it is very important to realise the accurate transformation between G(t) and J(t) []. Thus far, many approximated as well as numerical attempts [] have been made to properly handle this issue. However, the reliability and accuracy of the transformed results are largely affected by the fitted Prony series constants as well as the interconversion terms selected.

To alleviate the aforementioned difficulties, the fractional derivative viscoelastic models (FDVMs) attract more and more scholars around the world []. Recently, Xiao et al. [] demonstrated that the generalised Maxwell model and the fractional derivative Zener model (FDZM) can equivalently describe the material viscoelastic behaviours in the time and frequency domains. Further, the experimental data of viscoelastic materials can be better and more accurately approximated with FDVMs [], as a small number of model parameters necessitate being determined. More importantly, the exact interconversion between G(t) and J(t) can be guaranteed as well. To the best of the authors’ knowledge, most works regarding FDVMs in the open literature (see [] and the references therein) are mainly concerned with the experimental viscoelastic data tested by the conventional measures, while few studies address the problem involving nanoindentation, which motivates the generation of this paper. In what follows, we shall apply the FDZM to fitting the experimental data tested by two commonly applied shaped indenters, and derive the relevant data fitting formulas for the first time to determine the relaxation modulus and creep compliance by using the Laplace transformation and convolution theorem.

This paper is organised as follows: In the section below, the experimental tests by nanoindentation are first categorised into four distinct cases and the relevant viscoelastic governing formulations are then briefly recalled for this paper self-contained. Hereafter, the explicit formulas derived to determinate the linear viscoelastic moduli together with a fitting procedure are presented in Section 3. Section 4 is dedicated to examining the validity of our elaborated approach with the help of the experimental data as well as the fitting results in the open literature, and Section 5 is devoted to demonstrating the exact conversion of the two sets of viscoelastic moduli G(t) and J(t). Finally, a few concluding remarks are given.

2

Preliminaries

In this section, experimental tests by the nanoindentation are first categorised to facilitate the afterwards derivation and discussions, and some preliminary governing formulations are then briefly recalled for this paper self-contained.

2.1

Classification of the experimental tests

As already mentioned in the last section, both the relaxation modulus function G(t) and the creep compliance function J(t) can individually characterise the viscoelastic response of these elastomers. Hence, two feasible experimental strategies may be selected based on the practical requirements and the equipment available. Accordingly, we divide all these experimental tests on the basis of the objective function chosen into two main groups I and II, see Table . In addition, differently shaped indenters are commonly applied to carry out the practical tests. In this study, we primarily focus on the conical and spherical indenters. Therefore, each of the main groups is further categorised into two cases, see Table . We will arrange the subsequent derivations and discussions in the same sequence as the cases in Table .

Table 1 Classification of the nanoindentation tests

Group number	Objective function	Type of indenter	Case number
I	J(t)	conical indenter	①
spherical indenter	②
II	G(t)	conical indenter	③
spherical indenter	④

2.2

Mathematical preliminaries

Mathematical preliminaries about the interconversion of G(t) and J(t) and the FDZM are introduced in what follows.

2.2.1

Interconversion of the relaxation modulus and the creep compliance functions

The foregoing relaxation modulus function G(t) and the creep compliance function J(t) are related to each other through the well-known convolution integrals []: 1a ${\int }_0^{t}G\left(t-\tau \right)\frac{\mathrm{d}J\left(\tau \right)}{\mathrm{d}\tau }\mathrm{d}\tau =1,t>0$ 1b ${\int }_0^{t}J\left(t-\tau \right)\frac{\mathrm{d}G\left(\tau \right)}{\mathrm{d}\tau }\mathrm{d}\tau =1,t>0$where τ is a dummy variable.

Applying the Laplace transformation and convolution theorem to the integral formulas above, we obtain the following equivalent relation in the Laplace domain: 2 $\overline{G}\left(s\right)\overline{J}\left(s\right)s^2=1$Above, $\overline{G}\left(s\right)$ and $\overline{J}\left(s\right)$, respectively, represent the Laplace transform formula of G(t) and J(t), s refers to the transform parameter.

Theoretically, once any one of the modulus functions in () is ascertained, the other one may be retrieved by () without conducting additional experiments. Many pioneering contributions [, , ] have been made to realise the exact interconversion. However, only the approximated solutions can be obtained so far with these methods due to the ill-posed relationship of the Prony series. In 1999, Welch et al. [] introduced the fractional calculus analysis techniques into the quasi-static viscoelastic response of polymeric materials and achieved the exact interconversion of these two modulus functions analytically.

2.2.2

Fractional derivative viscoelastic model

Different from the Prony series-based viscoelastic models which apply the Newton dashpots to characterising the viscous material behaviour, FDVMs employ the new element, namely springpot, to model this time-dependent performance. Interested readers can consult [] for more information about FDVMs. Let us begin the introduction with the famous single parameter Mittag–Leffler function presented in [] 3a $E_{\alpha }\left(x\right)=\sum _{k=0}^{\mathrm{\infty }}\frac{x^k}{\mathrm{\Gamma }\left(\alpha k+1\right)}$where Г(x) is the Gamma function (see [, ]), α is the model coefficient with its value larger than 0, and x is a real or complex number. Equation () was further formulated in a generalised form (see []) as 3b $E_{\alpha ,\beta }\left(x\right)=\sum _{k=0}^{\mathrm{\infty }}\frac{x^k}{\mathrm{\Gamma }\left(\alpha k+\beta \right)}$with α and β being the two positive model coefficients.

It is now readily verify that the relation below holds for the generalised Mittag–Leffler function []: 4 $L^{-1}\left(\frac{s^{\alpha -\beta }}{s^{\alpha }+\lambda }\right)=t^{\beta -1}{E}_{\alpha ,\beta }\left(-\lambda t^{\alpha }\right)$Above, L⁻¹[f(s)] denotes the inverse Laplace transform of f(s). With the aid of (3) and (), the foregoing modulus function G(t) (or J(t)) of FDVM and $\overline{G}\left(s\right)$ (or $\overline{J}\left(s\right)$) can be transformed conveniently.

In this study, the commonly applied FDZM [] is taken as an illustration to elaborate the fitting strategy, and its differential constitutive equation is governed as 5 $\sigma \left(t\right)+a_1\frac{{\mathrm{d}}^{\alpha }\sigma }{\mathrm{d}{t}^{\alpha }}=b_1\frac{{\mathrm{d}}^{\alpha }\epsilon }{\mathrm{d}{t}^{\alpha }}+m\epsilon \left(t\right)$where α is a positive constant varying in the range (0, 1), and a₁, b₁ and m are positive real numbers. σ(t) and ɛ(t) separately represent the time-dependent stress and strain functions.

The creep compliance function J(t) as well as the relaxation modulus function G(t) of the system characterised by () can be derived by performing some transformations as 6a $J\left(t\right)=\frac{1}{m}+\left(\frac{a_1}{b_1}-\frac{1}{m}\right)E_{\alpha }\left[-{\left(\frac{m}{b_1}\right)}^{\alpha }{t}^{\alpha }\right]$ 6b $G\left(t\right)=m+\left(\frac{b_1}{a_1}-m\right)E_{\alpha }\left[-\frac{1}{a_1^{\alpha }}{t}^{\alpha }\right]$It is obvious that the elastic modulus functions above contain fewer coefficients compared with the ones involving the Prony series.

2.3

Some data-fitting formulas in the open literature

In this subsection, the indentation data-fitting formulas in [, , ] are briefly recalled to facilitate the subsequent discussions and comparisons.

2.3.1

Formulas regarding the creep compliance function

Let us first consider the two cases ① and ② in Table . The governing formula of the experimental tests for these two cases by Lu et al. [] can be uniformly expressed as 7 $h^n\left(t\right)=c{\int }_0^{t}J\left(t-\tau \right)\mathrm{d}P\left(\tau \right)$with h(t) being the indentation depth function obtained from the experimental test, and P(τ) the loading function defined as P(τ) = kτ. k refers to the loading rate. n and c are two indenter-shape-related coefficients, and, respectively, take the value n = 2 and c = 0.5π(1 − v²)tan θ for case ①, and n = 1.5 and c = 3(1 − v²)/($4\sqrt{R}$) for case ②. θ represents the angle between the cone generator and the substrate plane; R is the radius of the spherical indenter tip. v is the Poisson's ratio of the material under consideration and is assumed to be constant during the experiment time [, , ]. Note that (i) the contact circle radius of the conical indenter does not affect the parameter c, and interested readers are referred to [] for more details. (ii) The creep compliance in shear J_s(t) in [] equals 2(1 + ν)J(t), and the material Poisson's ratio is specified as a constant.

Then, combining the creep compliance function J(t) of the Prony series form below: 8 $J\left(t\right)=J_0+\sum _{i=1}^N{J}_i\left(1-{\mathrm{e}}^{-t/{\tau }_i}\right)$they obtained the explicit data-fitting relation between h(t) and J(t) as 9 $h^n\left(t\right)=ck\cdot \left[J_{\mathrm{\infty }}t-\sum _{i=1}^N{J}_i{\tau }_i\left(1-{\mathrm{e}}^{-t/{\tau }_i}\right)\right]$with J_∞ = (J₀ + J₁ + ⋯ + J_N).

In (), J₀, J₁, …, J_N are the material compliance parameters to be ascertained, and τ₀, τ₁, …, τ_N are the retardation times with N being the total number of items chosen from the Prony series.

The elastic moduli above can hereafter be determined by introducing the experimental data tested into () together with a fitting scheme. It is obvious that the accuracy of the fitting results is largely dependent on the value of N.

2.3.2

Formulas regarding the relaxation modulus function

When the two cases ③ and ④ of group II in Table intervene, the governing formulas in [] may also be written in a unified form such that 10 $P\left(t\right)=\omega \upsilon {\int }_0^{t}G\left(t-\tau \right){\left[h\left(\tau \right)\right]}^{\upsilon -1}\frac{\mathrm{d}h\left(\tau \right)}{\mathrm{d}\tau }\mathrm{d}\tau$where P(t) is the indentation load function, and h(τ) is the indentation depth function defined through h(τ) = κτ as the test is carried out with a constant displacement rate. κ is the displacement rate of the indenter. Parameters ω = 2/[π(1 − v²)tan θ] and υ = 2 are applied for case ③, and $\omega =4\sqrt{R}/\left[3(1-v^2)\right]$ and υ = 1.5 are for case ④. All other coefficients have the same physical meaning as those in the last subsection.

Invoking the following relaxation modulus function G(t) of the Prony series form: 11 $G\left(t\right)=G_{\mathrm{\infty }}+\sum _{i=1}^N{G}_i{\mathrm{e}}^{-{\lambda }_{i}t}$Huang and Lu [] derived the explicit fitting expression of P(t) for nanoindentation tests with the conical indenter (i.e. case ③) 12 $P\left(t\right)=\frac{4k^2}{\pi \left(1-v^2\right)\tan \theta }\left(\frac{1}{2}{G}_{\mathrm{\infty }}{t}^2+\sum _{i=1}^N\frac{G_i}{{\lambda }_i}\left(t-\frac{1}{{\lambda }_i}+{\mathrm{e}}^{-{\lambda }_{i}t}\right)\right)$where G_∞ and G_i are the relaxation moduli to be determined, and λ_i is the reciprocals of relaxation times.

Besides, the data-fitting formula for case ④ is also obtained in [] and described as 13 $P\left(t\right)=\frac{4k^{3/2}\sqrt{R}}{3\left(1-v^2\right)}\left(G_{\mathrm{\infty }}{t}^{3/2}+\sum _{i=1}^N{\int }_0^t{G}_i{\mathrm{e}}^{-{\lambda }_{i}t}\sqrt{t-\tau }\mathrm{d}\tau \right)$It is noted that P(t) in () is expressed in terms of G_∞ and G_i in an implicit manner. Hence, the integral operation would also affect the resulting fitting accuracy.

In [], Martynova replaced () with the empirical formula below to deduce the explicit expression for case ④: 14 $P\left(t\right)=c_0\left(1-{\mathrm{e}}^{-c_{1}t}\right)-c_0{c}_{1}t+c_2{t}^{3/2}$The coefficients c₀, c₁ and c₂ are first ascertained according to the experimental indention data, and then substituted into the following governing formula to determine the relaxation modulus in shear G_s(t): 15 $G_{\mathrm{S}}\left(t\right)=\frac{k^{-3/2}\left(1-v\right)}{\pi \sqrt{R}}\left[\frac{3}{4}\pi c_2-c_0{c}_1^{2}e\left(-c_{1}t\right){\int }_0^t\frac{e\left(c_1\tau \right)}{\sqrt{\tau }}\mathrm{d}\tau \right]$with G_s(t) = G(t)/(1 + v).

It is noticed that the explicit formula in () is derived based on the empirical relation (). Hence, its applicability for a wide range of practical viscoelastic materials should be further verified. Besides, the exact interconversion between G(t) and J(t) can be hardly achieved with this formulation due to its complexity.

3

Data-fitting formulas with FDZM together with the implementation

The key step for applying the FDZM to actual experimental data fittings is to establish the explicit relation between the indentation depth h(t) (or the indentation load P(t)) and the model coefficients of the creep compliance function J(t) (or the relaxation modulus G(t)). Hereafter, a general fitting scheme can be employed to determine the value of these coefficients.

3.1

Explicit fitting formulas with FDZM of all four cases in Table

Let us begin the expatiation with the Stieltjes convolution [] such that: f(t) and g(t) are two continuous functions defined over the range [0,∞) of the time t. When t is negative, both functions g(t) and f(t) equal to 0. Accordingly, the Stieltjes convolution of f(t) and dg(t) can be described by using Riemann–Stieltjes integral [] as 16 $f\left(t\right)\times \mathrm{d}g\left(t\right)=f(t)g(0)+{\int }_0^{t}f\left(t-\tau \right)\mathrm{d}g\left(\tau \right)$where g(0) denotes the function value of g(t) at the time instant t = 0.

The Laplace transform of the Stieltjes convolution in () is 17 $L\left[f\left(t\right)\times \mathrm{d}g\left(t\right)\right]=s\cdot \overline{f}(s)\cdot \overline{g}\left(s\right)$where L is the Laplace transform operator.

It is observed that the aforementioned governing equations () and () can be deemed as the special cases of the Stieltjes convolution. First, let us consider the governing formula () of the cases ① and ② in Table , and set the loading function as P(τ) = kτ. Thus, the relation P(0) = 0 holds, and () is just the Stieltjes convolution of J(t) and dP(τ). Applying the Laplace transform () to (), we obtain 18 $L\left[h^n\left(t\right)\right]=c\cdot L\left[J\left(t\right)\ast \mathrm{d}P\left(t\right)\right]=c\cdot s\cdot \overline{J}\left(s\right)\cdot \overline{P}(s)$The explicit expression of $\overline{J}\left(s\right)$ and $\overline{P}\left(s\right)$ can be written by invoking () and the definition of P(τ) as 19 $\overline{P}\left(s\right)=k\frac{1}{s^2},\overline{J}\left(s\right)=\frac{p_2}{s}+p_1\frac{s^{\alpha -1}}{s^{\alpha }+p_3}$with p₁ = (ma₁ − b₁)/(mb₁), p₂ = 1/m, and p₃ = (m/b₁)^α.

Substituting () into () gives 20 $L\left[h^n\left(t\right)\right]=ck\left[\frac{p_2}{s^2}+p_1\frac{s^{\alpha -2}}{s^{\alpha }+p_3}\right]$Introducing () into () and performing some mathematical simplifications finally yield 21 $h^n\left(t\right)=ck\cdot t\left[p_2+p_1{E}_{\alpha ,2}\left(-p_3{t}^{\alpha }\right)\right]$Up to now, the explicit expression for fitting the indentation data of the foregoing group I has been obtained and presented in (). It is clear that only four parameters, α, p₁, p₂ and p₃, necessitate being determined according to the experimental data tested.

Next, we recall the governing formula () of the group II in Table . This time, the indentation depth function is specified as h(τ) = κτ. It is easy to verify that () is the Stieltjes convolution of G(t) and $\mathrm{d}{\left[h(\tau )\right]}^{\upsilon }$, since ${\left[h(\tau )\right]}^{\upsilon }$ equals 0 at the instant τ = 0. Taking () into the Laplace transform (), we have 22 $L\left[P(t)\right]=\omega \cdot L\left[G\left(t\right)\ast \mathrm{d}{\left[h\left(t\right)\right]}^{\upsilon }\right]=\omega \cdot s\cdot \overline{G}\left(s\right)\cdot \overline{h}{\left(s\right)}^{\upsilon }$$\overline{G}\left(s\right)$ and $\overline{h}{\left(s\right)}^{\upsilon }$ are derived by taking into account of () and the definition of h(τ) as 23 $\overline{h}{\left(s\right)}^{\upsilon }={\kappa }^{\upsilon }\frac{\mathrm{\Gamma }\left(\upsilon +1\right)}{s^{\upsilon +1}},\overline{G}\left(s\right)=\frac{q_2}{s}+q_1\frac{s^{\alpha -1}}{s^{\alpha }+q_3}.$where q₁ = (b₁ − ma₁)/a₁, q₂ = m, and q₃ = (1/a₁)^α.

Introducing () into () delivers 24 $L\left[P\left(t\right)\right]=\eta \left[\frac{q_2}{s^{\upsilon +1}}+q_1\frac{s^{\alpha -1-\upsilon }}{s^{\alpha }+q_3}\right]$with η = Γ(υ + 1)ωκ^υ.

Conducting the inverse Laplace transformation of (), we eventually get the formula to fit the indentation data of the aforementioned group II 25 $P\left(t\right)=\eta t^{\upsilon }\left[\frac{q_2}{\mathrm{\Gamma }\left(\upsilon +1\right)}+q_1{E}_{\alpha ,\upsilon +1}\left(-q_3{t}^{\alpha }\right)\right]$Similar to (), () also contains four unknown coefficients α, q₁, q₂ and q₃ merely. Therefore, the fitting formulas derived in this study have less model coefficients compared with those in [, , ], and are expressed in an explicit manner for all indentation tests involving differently shaped indenters.

3.2

Indentation data fitting implementation

With the help of the explicit data-fitting formulas () and (), we can conveniently determine the material viscoelasticity characterised by the FDZM and tested by the nanoindentation. The whole fitting process is summarised as the flowchart shown in Fig. . At first, the nanoindentation tests are conducted, and the relevant experimental data are extracted. According to the data (h(t) or P(t)) obtained, the viscoelastic modulus function to be fitted as well as the group number in Table is ascertained. Next, the shape of the indenter applied is taken into account, and the relevant constant parameters in the governing formula () for Group I or () for Group II are selected. Afterwards, a general fitting method is employed to calculate the model coefficients' value. In our study, the least-square method in MATLAB is used to perform the data fitting process. Once the J(t) (or G(t)) expression is determined, its counterpart modulus function G(t) (or J(t)) can be computed by using the relation in () and ().

[IMAGE OMITTED. SEE PDF]

4

Validations and discussions

In what follows, several sets of indentation data together with the relevant predicted modulus functions in the open literature [, , ] are applied to demonstrate the validity of our new formulation and the data-fitting scheme in the last section.

4.1

Validation and discussion of the predictions with FDZM of group I

The indentation depth data h(t) together with the creep compliance in shear J_s(t) predicted in [] are applied to examine the validity of the elaborated approach. Two kinds of viscoelastic materials, i.e. polycarbonate (PC) and polymethyl methacrylate (PMMA), are employed for versatility. The PC material is taken as the first example, and the experimental data required for the curve fitting with () are given as: θ = 19.7°, v = 0.3 and k = 0.038 mN/s. Interested readers are referred to [] for more other experimental conditions. The mathematical expression of () fitted by our approach is given in (), and the relevant load-indentation depth curve together with the test data employed is plotted in Fig. . The creep compliance in shear J_s(t) determined is presented in () and visualised in Fig. 26a $h^2\left(t\right)=0.007481t\left[1.575-0.2752E_{0.8657,2}\left(-0.02022t^{0.8657}\right)\right]$ 26b $J_s\left(t\right)=1.575-0.2752E_{0.8657}\left(-0.02022t^{0.8657}\right)$

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

It is observed that the predicted load-indentation displacement results in Fig. match the test data by the conical indenter very well, though only a limited number of sample points measured from the graph in [] are applied to the curve fitting. The relevant viscoelastic modulus function J_s(t) (see Fig. ) determined by our approach is pretty close to the predictions in [] at the beginning, and gradually approaches the results by the conventional test after the instant t = 40 s. Further, both J_s(t) functions by our approach and in [] agree with the conventional evaluations with reasonably good accuracy. Hence, our approach is valid to fit the indentation data tested with the conical indenters and can correctly predict the associated creep compliance function. The difference between the predicted J_s(t) functions separately with the nanoindentation and the conventional tests is considered to be the influence of the indentation depth h(t). Nanoindentation tests are mainly applied to determine the local surface viscoelasticity of a certain material, hence, the plane stress state dominates the deformation of the experimental sample when h(t) is small and its effect is gradually depressed with the increment of h(t). Meanwhile, the general 3D governing formula becomes valid as h(t) increases to a certain critical value. However, the specimens applied in the conventional tensile or shear tests are commonly finite sized, hence, they are governed by the 3D elastic equations and the influence of the plane stress state is quite small or even ignorable.

In the second example, the nanoindentation data of PMMA in [] is considered. The experimental parameters applied in () are θ = 19.7°, v = 0.3 and k = 0.267 mN/s. Interested reader may consult [] for more other information. The explicit expression of () is provided in (), and Fig. further shows the predicted indentation depth curve by our approach against the test data in []. Fig. plots the creep compliance in shear J_s(t) in () associated with the curve predicted 27a $h^2\left(t\right)=0.05255t\left[0.9909-0.1409E_{0.5367,2}\left(-0.4982t^{0.5367}\right)\right]$ 27b $J_s\left(t\right)=0.9909-0.1409E_{0.5367}\left(-0.4982t^{0.5367}\right)$

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

Again, a good agreement is observed from Fig. between the load-indentation displacement curve obtained by our approach and the data tested by the nanoindentation. The J_s(t) functions predicted by our approach and in [] based on the nanoindentation data are in a reasonably good agreement with those by the conventional test as well. Slight difference between the evaluations of J_s(t) by the nanoindentation and the conventional test attributes to the influence of the indentation depth h(t), which is already explained in the first example. In addition, the J_s(t) evaluations with our approach get more closer to the conventional results with the increment of h(t), which is also consistent with the foregoing explanations.

Further, our approach is also readily applied to fit the nanoindentation data with spherical indenters and predict the relevant creep compliance function (i.e. case ②). We have searched the experimental data of this type in the open literature extensively, but quite few sets of test data with enough experimental information are available. Although some indentation results with the detailed experimental parameters are provided in [], we cannot retrieve their results with the relevant formula and the parameters provided. Interested readers are suggested to test the validity of the established approach with their own indentation data. It is believed that our approach is also valid for these nanoindentation tests, since the load-indentation displacement data tested with the spherical indenter are similar to those in Figs. and .

4.2

Validation and discussion of the predictions with FDZM of group II

In what follows, the nanoindentation results in [, ] are employed to further examine the validity of our elaborated approach in predicting the relaxation modulus G(t) of viscoelastic materials.

Let us begin the validation with the experimental data of the PC material tested with the conical indenter in [] (case ③ in Table ). The experiment parameters in () are θ = 19.7°, v = 0.3 and κ = 5 nm/s. The model coefficients (α, q₁, q₂ and q₃) are first determined with the help of the least-square method and the experimental results in Fig. , and given in (). Afterwards, they are used to calculate the relaxation modulus function G(t), and the resulting formula is provided in () and visualised in Fig. . The G(t) function by the conventional test in [] is also provided in Fig. for comparison 28a $\begin{array}{r}P\left(t\right)=9.769\times {10}^{-5}{t}^2\left[0.7785+0.3428E_{0.2634,3}\left(-0.8353t^{0.2634}\right)\right]\end{array}$ 28b $G\left(t\right)=1.557+0.3428E_{0.2634}\left(-0.8353t^{0.2634}\right)$

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

It is seen that the predicted load-indentation displacement curve in Fig. perfectly matches the test data by the nanoindentation. The relaxation modulus function G(t) determined by our approach in Fig. is very close to but slightly higher than that in [], and shows similar variation as the conventional result at the very beginning. Both G(t) functions predicted based on the nanoindentation data agree with the conventional solution with reasonably good accuracy. The difference between the evaluations by the nanoindentation and conventional test is due to the influence of the indentation depth h(t).

The viscoelastic response of the PMMA material tested with the conical indenter in [] is taken as the second example. The experiment parameters in () are given as θ = 19.7°, v = 0.3 and κ = 5 nm/s. The mathematical expression of () fitted by our approach is given in (), and the indentation depth curve obtained by our approach together with the experimental data adopted is plotted in Fig. . The relevant relaxation modulus function G(t) is presented in () and shown in Fig. 29a $P\left(t\right)=9.769\times {10}^{-5}{t}^2\left[1.249+0.4998E_{0.6,3}\left(-0.5259t^{0.6}\right)\right]$ 29b $G\left(t\right)=2.498+0.4998E_{0.6}\left(-0.5259t^{0.6}\right)$

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

Similar to the former example, an excellent agreement is observed again between the predicted load-indentation displacement curve by our approach and the experimental data in Fig. . Meanwhile, the relaxation modulus function G(t) determined by our approach shows perfect agreement with the prediction in [] at the beginning, and presents slight deviation as the time t increases and gets much closer to the solution by the conventional test. Both G(t) functions on the basis of the nanoindentation data are in good agreement with reasonable accuracy compared to the conventional solution.

Hereafter, our attention is mainly on examining the validity of our approach in predicting the relaxation modulus function G(t) with the experimental data tested by the spherical indenter, i.e. case ④ in Table . Two examples are considered here as well. The experimental test on the PMMA material in [] is taken as the first example. The radius of the indenter applied is R = 10 μm, the Poisson's ratio of PMMA is specified as v = 0.3, the displacement rate of the indenter is given as κ = 4.9 nm/s during the test. This time, () with υ = 1.5 and ω = 4R^1/2/[3(1 − v²)], is applied to implement the data fitting. The resulting formula is presented in (), and the curve predicted together with the experiment data is plotted in Fig. . The relaxation modulus in shear G_s(t) determined by the fitted curve is written in () and shown in Fig. 30a $P\left(t\right)=0.002746t^{1.5}\left[1.134+1.792E_{0.7216,2.5}\left(-0.06518t^{0.7216}\right)\right]$ 30b $G_s\left(t\right)=1.508+1.792E_{0.7216}\left(-0.06518t^{0.7216}\right)$

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

An excellent agreement is clearly observed from Fig. between the load-indentation displacement curve by our approach and the experimental data by the nanoindentation. The predicted G_s(t) functions by our approach and in [] based on the indentation data in Fig. are in a reasonably good agreement with that by the conventional test, and our result is slightly closer to the conventional solution compared with that in []. The difference between the evaluations of G_s(t) by the nanoindentation and the conventional test is ascribed to the influence of the indentation depth h(t), which is explained in the first example of the last subsection.

In this last example, the virtual experimental data given in [] is adopted for the validation. The shape of the indenter is also spherical, and its radius equals R = 10 μm. The material Poisson's ratio v is set as 0.3, and the displacement rate κ equals 4.9 nm/s. The model coefficients (α, q₁, q₂ and q₃) are first determined based on the data in Fig. and given in (), and then applied to compute the relaxation modulus function G_s(t) (see () and Fig. ) 31a $P\left(t\right)=0.002746t^{1.5}\left[1.484+0.2179E_{0.6058,2.5}\left(-0.0524t^{0.6058}\right)\right]$ 31b $G_s\left(t\right)=1.974+0.2179E_{0.6058}\left(-0.0524t^{0.6058}\right)$

[IMAGE OMITTED. SEE PDF]

[IMAGE OMITTED. SEE PDF]

It is seen from Fig. that the predicted load-indentation displacement curve in Fig. perfectly matches the virtual experimental data in [], and the relaxation modulus function G_s(t) evaluated by our approach is also in an excellent agreement with the conventional solution as well as the prediction in []. These results highlight the correctness and validity of our elaborated approach.

In summary, all foregoing examples indicate that our approach developed based on the FDZM is applicable in identifying the linear viscoelastic material properties (J(t) or G(t)) tested by nanoindentation with differently shaped indenters at the constant loading or displacement rates. Further, the data-fitting formulas () and () derived contain less unknown coefficients to be ascertained, and are thus more convenient compared with the methods adopting the Prony series and the empirical formula.

5

Exact interconversion between G(t) and J(t)

As mentioned earlier, the viscoelastic modulus function G(t) (or J(t)) described with Prony series can convert into its exact counterpart J(t) (or G(t)) only when the exponential terms chosen from Prony series are <3, while the data-fitting approach elaborated in this study can realise the analytical interconversion of G(t) and J(t) without this limitation. In what follows, three examples are presented to show this fact and demonstrate the importance of the exact interconversion between G(t) and J(t).

The relaxation modulus function G(t) presented in [] is taken as our first example 32 $G\left(t\right)=1+0.6{\mathrm{e}}^{-t/5}+0.4{\mathrm{e}}^{-t/20}$It is clear that there are only two exponential terms in (). Hence, the exact counterpart J(t) of G(t) can be derived by resorting to the analytic solution in [] as 33 $J\left(t\right)=1-0.34{\mathrm{e}}^{-t/29.155}-0.16{\mathrm{e}}^{-t/6.864}$Relation G(t) in () is approximated by the FDZM in () as 34 $G\left(t\right)=0.9751+1.0461E_{0.8482}\left(-0.1590t^{0.8482}\right)$The corresponding J(t) with respect to () is obtained by taking account of () and () as 35 $J\left(t\right)=0.0050+0.0471E_{0.7193}\left(-0.0462t^{0.7193}\right)$Fig. shows the predicted viscoelastic moduli by ()–(). It is observed that: (i) both the relaxation modulus G(t) and the converted creep compliance J(t) predicted by the FDZM are in good agreement with the relevant evaluations in []; (ii) the maximal relative error of either G(t) or J(t) over the whole time period under consideration is smaller than 1%. These results highlight the correctness and effectiveness of our established approach.

[IMAGE OMITTED. SEE PDF]

In the second example, the creep compliance function J(t) obtained by Xu et al. [] and shown below is employed 36 $\begin{array}{rl}J\left(t\right)& =0.7104-0.1872{\mathrm{e}}^{-(t/0.01)}-0.1823{\mathrm{e}}^{-(t/0.1)}\\ & -0.1361{\mathrm{e}}^{-t}-0.01664{\mathrm{e}}^{-(t/10)}-0.0003143{\mathrm{e}}^{-(t/100)}\end{array}$Since five exponential terms are applied in (), the exact expression of G(t) cannot be obtained by directly solving the interconversion integral formula. Xu et al. [] adopted the approximate method in [] and got the following approximated solution of G(t): 37 $\begin{array}{rl}G\left(t\right)& =3.345+6.917{\mathrm{e}}^{-205.4t}+1.898{\mathrm{e}}^{-14.63t}+0.8917{\mathrm{e}}^{-1.273t}\\ & +4.585\times {10}^{-31}{\mathrm{e}}^{-t}+0.08929{\mathrm{e}}^{-0.1027t}\end{array}$Then, J(t) characterised by the FDZM in () is applied to the description of () and delivers 38 $J\left(t\right)=0.7201+0.5324E_{0.6038}\left(-3.778t^{0.6038}\right)$Accordingly, the relaxation modulus function G(t) associated with J(t) in () is derived by using the relations () and () as 39 $G\left(t\right)=3.388+9.611E_{0.6038}\left(-14.49t^{0.6038}\right)$The resulting predictions of J(t) and G(t) are plotted in Fig. . As expected, all J(t) values predicted by our approach agree well with the results in [], which indicates that the FDZM model adopted can properly depict the viscoelastic deformation of Example 2 as well. To be specified, the relaxation modulus function G(t) by our approach is in good agreement with that in [] at first, but presents an obvious deviation after the time instant 1 s. The reason is considered to be the fact that G(t) by our approach is the exact transformation of J(t) in (), while G(t) in [] is only an approximated conversion of J(t) in ().

[IMAGE OMITTED. SEE PDF]

In this last example, we adopted the creep compliance function J(t) in [] as an illustration 40 $\begin{array}{rl}J\left(t\right)& =0.0149-0.000526{\mathrm{e}}^{-t}-0.000167{\mathrm{e}}^{-(t/10)}\\ & -0.001841{\mathrm{e}}^{-(t/100)}-0.005101{\mathrm{e}}^{-(t/1000)}\\ & -0.007227{\mathrm{e}}^{-(t/10000)}\end{array}$In [], Ebrahimi et al. employed an approximation method based on the power-law interrelation to compute the G(t) value from the relation in (), and no explicit expression of G(t) is provided there. Interested readers can consult [] for the detailed implementation.

Applying the FDZM in ()–(), we have 41 $J\left(t\right)=0.02006-0.02E_{0.5106}\left(-0.01135t^{0.5106}\right)$The explicit expression of G(t) corresponding to J(t) in () is then obtained by using () and () 42 $G\left(t\right)=49.95+1.628\times {10}^4{E}_{0.5106}\left(-3.712t^{0.5106}\right)$Fig. shows the resulting evaluations of the viscoelastic moduli in ()–() together with G(t) by the approximation method in []. Similar to the former example, the J(t) curve with our approach and that in [] are overlapped with each other, which implies that both functions can correctly describe the variation of the creep compliance J(t). However, the G(t) function predicted with the power-law interrelation in [] gets small much quicker than our results with the increment of time t ascribed to its approximation essence, and a huge difference is observed between these two kinds of predictions after 200 s. It is thus concluded that the exact interconversion between J(t) and G(t) is of vital importance for practical failure analysis, structure design and other applications.

[IMAGE OMITTED. SEE PDF]

6

Conclusions

In this study, the FDZM is applied to fit the experimental viscoelastic data tested by nanoindentation, and the explicit data-fitting formulas for determining both the relaxation modulus G(t) and the creep compliance J(t) with conical and spherical nanoindenters are derived and given in () and (). Based on these data-fitting formulas together with the scheme established, the material viscoelasticity characterised by the FDZM can be conveniently and correctly identified according to the experimental tests, which has been demonstrated by the predicted results in Section 4 and the comparisons with those in the open literature [, , ]. In addition, the exact interconversion between G(t) and J(t) is also retained with the new approach. Hence, the two modulus functions G(t) and J(t) can be obtained with one nanoindentation test, which is of particular importance for actual applications.

Although only the FDZM is considered in this paper, fractional viscoelastic models of other types may also be adopted to characterise the material viscoelasticity by deriving the corresponding data-fitting formulas and using the data-fitting scheme herein.

7

Acknowledgments

This project was supported by National Natural Science Foundation of China (grant nos. 51290291, 51675447, 51305362, 11772274), the Young Scientific Innovation Team of Science and Technology of Sichuan (grant no. 2017TD0017) and the Fundamental Research Funds for the Central Universities (grant no. 2682016CX024). The financial support from the China Scholar Council is also gratefully acknowledged.

References

[1]

Boltzmann, L.: ‘Zur Theorie der elastischen Nahwirkung’, Ann. Phys., 1878, 241, (11), pp. 430–432 (doi: 10.1002/andp.18782411107)

[2]

Jafarzadeh, T., Erfan, M., Watts, D.: ‘Creep and viscoelastic behaviour of human dentin’, J. Dent. Tehran Univ. Med. Sci., 2004, 1, (1), pp. 5–14

[3]

Syu, Y.J., Lin, R.L., Liou, N.S.: ‘The modeling of time dependent mechanical properties of cervine Enamel’, Mech. Biol. Syst. Mater., 2015, 7, pp. 33–36

[4]

Lee, E.H., Radok, J.R.M.: ‘The contact problem for viscoelastic bodies’, J. Appl. Mech., 1960, 27, (3), pp. 438–444 (doi: 10.1115/1.3644020)

[5]

Zhang, Q., Gao, S., Min, J. et al.: ‘Graded viscoelastic behavior of human enamel by nanoindentation’, Mater. Lett., 2016, 179, pp. 126–129 (doi: 10.1016/j.matlet.2016.05.034)

[6]

Ashrafi, H., Shariyat, M.: ‘A nano‐indentation identification technique for viscoelastic constitutive characteristics of periodontal ligaments’, J. Biomed. Phys. Eng., 2016, 6, (2), pp. 109–118

[7]

Mencík, J., He, L.H., Swain, M.V.: ‘Determination of viscoelastic–plastic material parameters of biomaterials by instrumented indentation’, J. Mech. Behav. Biomed. Mater., 2009, 2, (4), pp. 318–325 (doi: 10.1016/j.jmbbm.2008.09.002)

[8]

Lucca, D.A., Herrmann, K., Klopfstein, M.J.: ‘Nanoindentation: measuring methods and applications’, CIRP Ann.‐Manuf. Technol., 2010, 59, (2), pp. 803–819 (doi: 10.1016/j.cirp.2010.05.009)

[9]

Lu, H., Wang, B., Ma, J. et al.: ‘Measurement of creep compliance of solid polymers by nanoindentation’, Mech. Time‐Depend. Mater., 2003, 7, (3), pp. 189–207 (doi: 10.1023/B:MTDM.0000007217.07156.9b)

[10]

Huang, G., Lu, H.: ‘Measurement of Young's relaxation modulus using nanoindentation’, Mech. Time‐Depend. Mater., 2006, 10, (3), pp. 229–243 (doi: 10.1007/s11043‐006‐9020‐3)

[11]

Oyen, M.L.: ‘Spherical indentation creep following ramp loading’, J. Mater. Res., 2005, 20, (8), pp. 2094–2100 (doi: 10.1557/JMR.2005.0259)

[12]

Peng, G., Zhang, T., Feng, Y. et al.: ‘Determination of shear creep compliance of linear viscoelastic‐plastic solids by instrumented indentation’, Polym. Test., 2012, 31, (8), pp. 1038–1044 (doi: 10.1016/j.polymertesting.2012.07.007)

[13]

Martynova, E.: ‘Determination of the properties of viscoelastic materials using spherical nanoindentation’, Mech. Time‐Depend. Mater., 2016, 20, (1), pp. 85–93 (doi: 10.1007/s11043‐015‐9285‐5)

[14]

Ebrahimi, G.M., Saleh, M., Gonzalez, M.A.M.: ‘Interconversion between viscoelastic functions using the Tikhonov regularisation method and its comparison with approximate techniques’, Road Mater. Pavement Des., 2014, 15, (4), pp. 820–840 (doi: 10.1080/14680629.2014.924428)

[15]

Loy, R., Anderssen, R.: ‘On the construction of Dirichlet series approximations for completely monotone functions’, Math. Comput., 2014, 83, (286), pp. 835–846 (doi: 10.1090/S0025‐5718‐2013‐02725‐1)

[16]

Park, S.W., Kim, Y.R.: ‘Interconversion between relaxation modulus and creep compliance for viscoelastic solids’, J. Mater. Civ. Eng., 1999, 11, (1), pp. 76–82 (doi: 10.1061/(ASCE)0899‐1561(1999)11:1(76))

[17]

Anderssen, R.S., Davies, A.R., De Hoog, F.R.: ‘On the interconversion integral equation for relaxation and creep’, ANZIAM J., 2007, 48, pp. 346–363 (doi: 10.21914/anziamj.v48i0.64)

[18]

Honerkamp, J., Weese, J.: ‘Tikhonovs regularization method for ill‐posed problems’, Contin. Mech. Thermodyn., 1990, 2, (1), pp. 17–30 (doi: 10.1007/BF01170953)

[19]

Fernández, P., Rodríguez, D., Lamela, M.J. et al.: ‘Study of the interconversion between viscoelastic behaviour functions of PMMA’, Mech. Time‐Depend. Mater., 2011, 15, (2), pp. 169–180 (doi: 10.1007/s11043‐010‐9128‐3)

[20]

Park, S.W., Schapery, R.A.: ‘Methods of interconversion between linear viscoelastic material functions. Part I – a numerical method based on Prony series’, Int. J. Solids Struct., 1999, 36, (11), pp. 1653–1675 (doi: 10.1016/S0020‐7683(98)00055‐9)

[21]

Nikonov, A., Davies, A.R., Emri, I.: ‘The determination of creep and relaxation functions from a single experiment’, J. Rheol., 2005, 49, (6), pp. 1193–1211 (doi: 10.1122/1.2072027)

[22]

Sorvari, J., Malinen, M.: ‘Numerical interconversion between linear viscoelastic material functions with regularization’, Int. J. Solids Struct., 2007, 44, (3), pp. 1291–1303 (doi: 10.1016/j.ijsolstr.2006.06.029)

[23]

Sorvari, J., Malinen, M.: ‘On the direct estimation of creep and relaxation functions’, Mech. Time‐Depend. Mater., 2007, 11, (2), pp. 143–157 (doi: 10.1007/s11043‐007‐9038‐1)

[24]

Loy, R.J., De Hoog, F.R., Anderssen, R.S.: ‘Interconversion of Prony series for relaxation and creep’, J. Rheol., 2015, 59, (5), pp. 1261–1270 (doi: 10.1122/1.4929398)

[25]

Zhuravkov, M.A., Romanova, N.S.: ‘Review of methods and approaches for mechanical problem solutions based on fractional calculus’, Math. Mech. Solids, 2016, 21, (5), pp. 595–620 (doi: 10.1177/1081286514532934)

[26]

Sasso, M., Palmieri, G., Amodio, D.: ‘Application of fractional derivative models in linear viscoelastic problems’, Mech. Time‐Depend. Mater., 2011, 15, (4), pp. 367–387 (doi: 10.1007/s11043‐011‐9153‐x)

[27]

Sapora, A., Cornetti, P., Carpinteri, A. et al.: ‘The use of fractional calculus to model the experimental creep‐recovery behavior of modified bituminous binders’, Mater. Struct., 2016, 49, (1–2), pp. 45–55 (doi: 10.1617/s11527‐014‐0473‐6)

[28]

Zhang, H., Wang, Y., Fatemi, M. et al.: ‘Assessing composition and structure of soft biphasic media from Kelvin–Voigt fractional derivative model parameters’, Meas. Sci. Technol., 2017, 28, (3), p. 035703 (doi: 10.1088/1361‐6501/aa5531)

[29]

Hernández‐Jiménez, A., Hernández‐Santiago, J., Macias‐Garcıa, A. et al.: ‘Relaxation modulus in PMMA and PTFE fitting by fractional Maxwell model’, Polym. Test., 2002, 21, (3), pp. 325–331 (doi: 10.1016/S0142‐9418(01)00092‐7)

[30]

Friedrich, C.H.R.: ‘Relaxation and retardation functions of the Maxwell model with fractional derivatives’, Rheol. Acta, 1991, 30, (2), pp. 151–158 (doi: 10.1007/BF01134604)

[31]

Liu, J.G., Xu, M.Y.: ‘Higher‐order fractional constitutive equations of viscoelastic materials involving three different parameters and their relaxation and creep functions’, Mech. Time‐Depend. Mater., 2006, 10, (4), pp. 263–279 (doi: 10.1007/s11043‐007‐9022‐9)

[32]

Bagley, R.L., Torvik, J.: ‘Fractional calculus – a different approach to the analysis of viscoelastically damped structures’, AIAA J., 1983, 21, (5), pp. 741–748 (doi: 10.2514/3.8142)

[33]

Xiao, R., Sun, H., Chen, W.: ‘An equivalence between generalized Maxwell model and fractional Zener model’, Mech. Mater., 2016, 100, pp. 148–153 (doi: 10.1016/j.mechmat.2016.06.016)

[34]

Ferry, J.D., Rice, S.A.: ‘Viscoelastic properties of polymers’ (Wiley, New York, 1961)

[35]

Welch, S.W.J., Rorrer, R.A.L., Duren, Jr., R.G.: ‘Application of time‐based fractional calculus methods to viscoelastic creep and stress relaxation of materials’, Mech. Time‐Depend. Mater., 1999, 3, (3), pp. 279–303 (doi: 10.1023/A:1009834317545)

[36]

Mainardi, F.: ‘Fractional calculus and waves in linear viscoelasticity: an introduction to mathematical models’ (World Scientific, London, 2010, 1st edn.), pp. 61–63

[37]

Mittag‐Leffler, G.M.: ‘Sur la nouvelle function E_α(x)’, Comptes Rendus, 1903, 137, pp. 554–558

[38]

Kilbas, A.A., Ssigo, M., Saxena, R.K.: ‘Generalized Mittag‐Leffler function and generalized fractional calculus operators’, Integr. Transforms Spec. Functions, 2004, 15, (1), pp. 31–49 (doi: 10.1080/10652460310001600717)

[39]

Koeller, R.C.: ‘Polynomial operators, stieltjes convolution, and fractional calculus in hereditary mechanics’, Acta Mech., 1986, 58, pp. 251–264 (doi: 10.1007/BF01176603)

[40]

Polyanin, A.D., Manzhirov, A.V.: ‘Handbook of integral equations’ (CRC Press, Boca Raton, 2008, 2nd edn.), pp. 149–151

[41]

Xu, T., Yoo, J.H., Babu, S. et al.: ‘Characterization of the mechanical behavior of SU‐8 at microscale by viscoelastic analysis’, J. Micromech. Microeng., 2016, 26, (10), p. 105001 (doi: 10.1088/0960‐1317/26/10/105001)

[42]

Luk‐Cyr, J., Crochon, T., Lévesque, M.: ‘Interconversion of linearly viscoelastic material functions expressed as Prony series: a closure’, Mech. Time‐Depend. Mater., 2013, 17, (1), pp. 53–82 (doi: 10.1007/s11043‐012‐9176‐y)