Rheological representation of fractional order viscoelastic material models

Abstract

We develop rheological representations, i.e., discrete spectrum models, for the fractional derivative viscoelastic element (fractional dashpot or springpot). Our representations are generalized Maxwell models or series of Kelvin-Voigt units, which, however, maintain the number of parameters of the corresponding fractional order model. Accordingly, the number of parameters of the rheological representation is independent of the number of rheological units. We prove that the representations converge to the corresponding fractional model in the limit as the number of units tends to infinity. The representations extend to compound fractional derivative models such as the fractional Maxwell model, fractional Kelvin-Voigt model, and fractional standard linear solid. Computational experiments show that the rheological representations are accurate approximations of the fractional order models even for a small number of units.

Appendices

Appendix A: Frequency representation

A review of classical material is provided here for the sake of the reader and to clarify notation. For a strain loading of $\epsilon(t)=\epsilon_0e^{i\omega t}$, where ε ₀ is a nonzero real coefficient, the stress as determined from (2) can be written in the form $\sigma(t)=\sigma_0e^{i\omega t}$, where $G^{*}=\sigma_0/\epsilon_0$ is the complex shear modulus, then

$$\begin{array}{lll} G^{*}(\omega)=G(0)+\int_{0}^{\infty} \dot{G}(t) e^{-i\omega t}\,{\rm d}t. \label{nishu-4} \end{array}$$

(75)

For example, if G is a constant function, then G* is the same constant. If σ ₀ is written in polar form as $s_0e^{i\delta}$, then δ is the angle by which the strain lags behind the stress. The complex modulus can be written as G*(ω) = G′(ω) + iG′′(ω) where G′(ω) and G′′(ω) are known as the storage and loss moduli, respectively.

Similarly the complex compliance is obtained by the following equation

$$\begin{array}{lll} J^{*}(\omega)=J'(\omega)-iJ''(\omega)=J(0)+\int_{0}^{\infty} \dot{J}(t) e^{-i\omega t}\,{\rm d}t. \label{nishu-5} \end{array}$$

(76)

J′(ω) and J′′(ω) are called the storage and loss compliance respectively. The loss tangent, tanδ, is defined as G′′(ω)/G′(ω) or equivalently as J′′(ω)/J′(ω).

Alternatively, the complex modulus and complex compliance can be derived from the differential equation describing the viscoelastic model. Assume the equation takes the form P(D)[σ] = Q(D)[ε], where P and Q are polynomials and D is the time-differentiation operator. In the case of fractional models, P and Q involve fractional powers of D. Upon the application of a harmonic strain $\epsilon^*=\epsilon_0e^{i\omega t}$, a steady-state response of the form $\sigma^*=\sigma_0e^{i\omega t}$ is obtained as above. The complex modulus is then given by

$$ G^*=\frac{\sigma^*}{\epsilon^*}=\frac{\sigma_0}{\epsilon_0}=\frac{Q(i\omega)}{P(i\omega)} \label{eq:Gstar_from diffeq} $$

(77)

and the complex compliance J* is the reciprocal.

Appendix B: Data-fitting methodology

In our computational experiments reported in the “Introduction” and in more detail in “Computational results,” we fitted various models (FMM, FKM, FSLS, and the generalized Maxwell model) to experimental data. Although the fitting of fractional derivative models to experimental data is not the principal concern of this paper, we briefly describe the methodology used to find the optimal model parameters in this appendix.

In order to fit the three parameters J ₀,η,α of the FMM to experimentally measured storage and loss compliances, we first observe that the fitting problem becomes linear in J ₀ and η once α is known. Thus, our fitting algorithm loops over 200 evenly-spaced choices of α between 0.0001 and 0.9999; for each α, it seeks the optimal choice for the other two parameters. Then the parameter set that yields the smallest overall residual is kept. The linear problem for η and J ₀ is given as follows. Let the experimentally measured storage and loss compliances be denoted by triples (ω _i ,J′ _i ,J′′ _i ) for i = 1,...,M. Then parameters J ₀ and η are chosen to minimize

$$ \sum\limits_{i=1}^M\left(\frac{J'(\omega_i)}{J'_i}-1\right)^2 +\sum\limits_{i=1}^M\left(\frac{J'(\omega_i)}{J'_i}-\frac{J''(\omega_i)}{J''_i} \right)^2 ,\label{fmmlls1} $$

(78)

where J′(·) and J′′(·) are given by (55) and (56). This particular linear least-squares formulation was chosen over the perhaps more obvious formula

$$ \sum\limits_{i=1}^M\left(J'(\omega_i)-J'_i\right)^2 +\sum\limits_{i=1}^M\left(J''(\omega_i)-J''_i\right)^2 \label{fmmlls2} $$

(79)

for the following reasons. The first term of (78) optimizes the relative error of J _i ′ versus J′(ω) whereas the first term of (79) optimizes the absolute difference, and hence attaches much greater significance to the terms in which J _i ′ is large. It seems preferable that a good model should have a low relative error over the whole range of frequencies. The second term of (78) is chosen to optimize the fit to the loss tangent, i.e., to make J′′(ω _i )/J′(ω _i ) close to J′′ _i /J′ _i . This is important because in practice the storage compliance and the loss tangent are the most significant viscoelastic properties in applications of damping and seismic loading.

In the case of fitting the FMM to time-dependent creep data, that is, data of the form (t ₁,J ₁), ..., (t _M ,J _M ), the algorithm is slightly simpler: one loops over values of α and for each choice of α, the optimal η and J ₀ are chosen to minimize

$$ \sum\limits_{i=1}^M \left(\frac{J(t_i)}{J_i}-1\right)^2. $$

We do not report on fitting this kind of data herein.

For fitting the FKM to storage and loss modulus experimental data (ω ₁,G′₁,G′′₁), ..., (ω _M , G′ _M ,G′′ _M ), again we loop over a range of possible α’s and for each α, we choose G _∞ and η to minimize

$$ \sum\limits_{i=1}^M\left(\frac{G'(\omega_i)}{G'_i}-1\right)^2 +\sum\limits_{i=1}^M\left(\frac{G'(\omega_i)}{G'_i}-\frac{G''(\omega_i)}{G''_i} \right)^2 ,\label{fkmlls1} $$

(80)

where G′(·) and G′′(·) are given by (62) and (63).

Finally, for fitting the four-parameter fractional SLS to storage and loss modulus experimental data, we observe that the fractional SLS is linear in the parameters B,C if the parameters A,α are fixed. Therefore, we loop over 200 possible choices of α and 401 possible choices of A (hence 200 × 401 = 80,200 least squares problems are solved in total) to find B,C that minimize (80), in which G′ and G′′ are defined by (68) and (69). We also use a discrete search for the best A after first identifying a reasonable search domain. We adopted the following procedure for determining the domain of A. First, a fractional Kelvin model was fitted to the data using the procedure in the previous paragraph. The fractional Kelvin model corresponds to the special case of fractional standard linear solid in which G ₀ = ∞ (and hence A = 0 by (70)). Once η*, α* and $G_d^*$ are obtained from the fractional Kelvin model fit, it appears from (70) that the maximum value of A should be around $\eta^*/G_d^*$. Therefore, a selection of 400 A values spaced geometrically with ratio of 1.2, with maximum value $100\eta^*/G_d^*$ and with minimum value $100\eta/(1.2^{399}G_d)$ is tried. In addition, A = 0 is tried since the fractional Kelvin model is a special case of the fractional SLS.

The last model fitted to the data was a generalized Maxwell model. This model was fitted for the purpose of the comparison made previously in the “Introduction” between generalized Maxwell model fits and fractional models. The problem of fitting a generalized Maxwell model to experimental data is nonlinear, and there are several published algorithms (see, e.g., Tschoegl 1989) for addressing this problem. Since this topic is not the focus of our paper, we instead employed a fairly brute-force fitting procedure. First, recall that the generalized Maxwell model is a model with one spring and N Maxwell units (spring-dashpot series) in parallel. The relaxation function for N elements is

$$\begin{array}{lll} G(t)= E + \sum\limits_{i=1}^{N}g_{i}e^{-t/\tau_{i}}, \end{array}$$

(81)

where g _i ,τ _i are the discrete elastic moduli and relaxation times, respectively, of the material, and E is the equilibrium modulus. The complex moduli are derived as

$$\begin{array}{lll} && G'(\omega)= E + \sum\limits_{i=1}^{N}\frac{\omega^{2} \tau_{i}^{2}g_{i}}{1+\omega^{2} \tau_{i}^{2}},\label{GMMsto}\end{array}$$

(82)

$$\begin{array}{lll} &&G''(\omega)= \sum\limits_{i=1}^{N}\frac{\omega \tau_{i}g_{i}}{1+\omega^{2} \tau_{i}^{2}}.\label{GMMlos} \end{array}$$

(83)

Observe that G′ and G′′ are linear in E and g _i . Therefore, for fitting the N = 4 case, we looped over a range of possible choices of τ ₁,...,τ ₄ and for each such quadruplet, we solved a linear least squares problem of selecting E,g ₁,...,g ₄ to minimize (80), where G′(·) and G′′(·) are given by (82) and (83). Each τ _i was allowed to take on one of 60 values logarithmically spaced between $1/(\omega_{\min} e^{-2})$ and $1/(\omega_{\max} e^2)$, where ω _min and ω _max are the minimum and maximum values of frequency in the data. Thus, 60×59×58×57/(1×2×3 ×4) = 487635 linear least squares problems were solved. We also checked that increasing this number 60 did not seem to significantly improve the fit. It should be noted that this procedure could in principle yield either positive or (unphysical) negative values of the g _i ’s and E. In the optimal fit from this procedure for the data reported herein, however, all the parameters were positive.