Generalized Fractional Derivative Anisotropic Viscoelastic Characterization

Abstract

Isotropic linear and nonlinear fractional derivative constitutive relations are formulated and examined in terms of many parameter generalized Kelvin models and are analytically extended to cover general anisotropic homogeneous or non-homogeneous as well as functionally graded viscoelastic material behavior. Equivalent integral constitutive relations, which are computationally more powerful, are derived from fractional differential ones and the associated anisotropic temperature-moisture-degree-of-cure shift functions and reduced times are established. Approximate Fourier transform inversions for fractional derivative relations are formulated and their accuracy is evaluated. The efficacy of integer and fractional derivative constitutive relations is compared and the preferential use of either characterization in analyzing isotropic and anisotropic real materials must be examined on a case-by-case basis. Approximate protocols for curve fitting analytical fractional derivative results to experimental data are formulated and evaluated.

1. Introduction

A vast body of literature indicates that the predominant choice in viscoelastic research, analysis and design is for constitutive relations consisting of integral relations, which are the Green’s functions associated with integral order differential equations [1,2,3,4,5]. However, an alternative representation for isotropic viscoelasticity in terms of fractional derivative calculus has been developed in [6,7,8,9,10,11,12,13,14] and generalized [15]. This concept, which is reflected in a substantial number of publications, is based on the premise that the derivatives in viscoelastic constitutive relations need not necessarily be of integer order, but may be of some fractional order with similar translations to integral representations. The mathematical foundation of fractional calculus and of the attending Laplace and Fourier transforms is well developed and may be found in [16,17,18,19,20,21,22,23], to mention only a few. In [24] a comprehensive historical review of linear isotropic viscoelastic characterizations based on fractional derivatives has been presented.

Volterra [25] was the first to extend Boltzmann’s isotropic superposition principle [26] to anisotropic viscoelasticity. Subsequently, in [27,28] anisotropic viscoelastic stress-strain relations were formulated. The mathematical and physical foundations of anisotropic viscoelasticity firmly rest on the well established isotropic viscoelastic concepts as shown in [29,30]. The elastic-viscoelastic integral transform correspondence principle (analogy) was extended to anisotropic viscoelasticity and latter’s proper mechanical modeling was formulated in [31]. In [31] it was shown that the powerful tool of separation of variables available under proper conditions in isotropic viscoelasticity is inadmissible in anisotropic viscoelasticity, leaving only Laplace (LT) and Fourier (FT) transforms, finite element methods (FEM) and a limited number of analytical solutions such as asymptotic approaches capable of producing analytical anisotropic solutions. Presently, the main applications of anisotropic viscoelasticity are in the areas of fiber/high polymer composites, fiber reinforced concrete, filament wound vessels, textiles and tires.

In the present paper the concept of single fractional derivative isotropic viscoelastic constitutive relations is formalized for anisotropic viscoelasticity in terms of many parameter generalized Kelvin models. Both differential and integral stress-strain equations and expressions for the associated complex moduli, compliances and creep functions are derived and their applicability and usefulness are evaluated. Viscoelastic integral stress-strain laws are always preferable to differential ones, because they are relatively easier to handle and computationally yield more accurate results.

2. Analysis

2.1. Material Property Dependence on Temperature

Viscoelastic media, be they metals at elevated temperatures [32,33] or high polymers [34], are well known to exhibit material properties that are extremely temperature sensitive. In general, their coefficients of viscosity vary one order of magnitude for roughly every 20 °C (Figure 1), which is described by the Arrhenius relation

$$\eta ={\eta }_{0}exp\left(-\frac{Q_a}{R\Theta }\right)$$

(1)

where $Q_a$ is the activation energy, R the universal gas constant and Θ the material temperature in degrees Kelvin.

Figure 1. Coefficients of viscosity for metals and composites.

Figure 1. Coefficients of viscosity for metals and composites.

Consider a Cartesian system $x=\{x_1,x_2,x_3\}$. The linear viscoelastic constitutive relations may be written in the form:

$$\left\{\sigma (x,t)\right\}=\underset{-\infty }{\overset{t}{\int }}\left[E(x,t,t^{\prime })\right]\left\{\frac{\partial ϵ(x,t^{\prime })}{\partial t^{\prime }}\right\}d{t}^{\prime }-\underset{-\infty }{\overset{t}{\int }}\left[E^T(x,t,t^{\prime })\right]\left\{\frac{\partial \left[\alpha T(x,t^{\prime })\right]}{\partial t^{\prime }}\right\}d{t}^{\prime }$$

(2)

These relations may be then classified into categories as seen in Table 1. Consequently, the first three cases are relatively “simple” integral relations, whereas the third one, as will be seen subsequently, presents considerable difficulties. These issues will be analyzed in detail in the next Sections.

Table 1. Temperature influence on constitutive relations.

Temperature	Modulus	Type of Constitutive Relations, Equation (2)	Type of Material
$T_0$	$E(t-t,{}^{\prime }{T}_0)$	convolution	homogeneous
$T_0$	$E(x,t-t,{}^{\prime }{T}_0)$	convolution	non-homogeneous
$T\left(x\right)$	$E[x,t-t,{}^{\prime }T\left(x\right)]$	convolution	non-homogeneous
$T(x,t)$	$E[x,t,t,{}^{\prime }T(x,t^{\prime })]$	non-convolution	non-homogeneous

Table 1. Temperature influence on constitutive relations.

Temperature	Modulus	Type of Constitutive Relations, Equation (2)	Type of Material
$T_0$	$E(t-t,{}^{\prime }{T}_0)$	convolution	homogeneous
$T_0$	$E(x,t-t,{}^{\prime }{T}_0)$	convolution	non-homogeneous
$T\left(x\right)$	$E[x,t-t,{}^{\prime }T\left(x\right)]$	convolution	non-homogeneous
$T(x,t)$	$E[x,t,t,{}^{\prime }T(x,t^{\prime })]$	non-convolution	non-homogeneous

2.2. General Concepts—Isotropic Materials

In order to systematically formulate anisotropic viscoelastic fractional derivative constitutive relations, it is necessary to first make a few general comments about isotropic relations with constant $\left(T_0\right)$ and spatially varying $\left(T(x\right))$ temperatures (see Table 1). In [15] generalized isotropic fractional derivative (FD) Maxwell models (GFDMM) are examined, while in [1,2,3] integer derivative (ID) differential and integral equations (DE and IE) for isotropic viscoelasticity in terms of generalized Maxwell and Kelvin models (GMM and GKM). A parallel development for integer order derivative DE and IE anisotropic viscoelastic responses has been formulated in [35]. In general, even though such modeling representations are interrelated and, therefore, interchangeable, GMM are better suited for numerical analyses, while GKM are more useful for achieving mathematical and physical understanding of material behavior. The use of GKM or GMM and GFDKM or GFDMM guarantees that any additive combination automatically satisfies thermodynamic principles.

Figure 2. Generalized Kelvin Model.

Figure 2. Generalized Kelvin Model.

Consider a GKM or GFDKM as shown in Figure 2. The single spring represents instantaneous elastic responses associated with “solids” and the single dashpot represents long time behavior tending to large or infinite displacements for long time periods. For isotropic cases, it is most useful to write separate constitutive relations for changes in shape and in volume, since real materials behave distinctly during these changes. Using Einstein’s Cartesian tensor notation and summation convention, all relations to follow are formulated in terms of non-dimensionalized material parameters and state variables. The mean stresses σ and strains ϵ are defined as

$$\sigma ={\sigma }_{ii}/3\mathrm{and}ϵ={ϵ}_{ii}/3$$

(3)

with stress and strain deviators

$$S_{kl}={\sigma }_{kl}-{\delta }_{kl}\sigma \mathrm{and}{E}_{kl}={ϵ}_{kl}-{\delta }_{kl}ϵ$$

(4)

The various elements of a GFDKM (Figure 2) have the following constitutive relations

$$\begin{array}{cc}\hfill \text{elastic}\text{spring}⟹S_{kl}& =2G_0{E}_{kl}^{\left(o\right)}\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill \text{each}\text{FDKM}⟹S_{kl}& =\underset{{\mathrm{n}}^{\mathrm{th}}\mathrm{dashpot}}{\underbrace{2{\eta }_n{D}^{{\delta }_n}{E}_{kl}^{\left(n\right)}}}+\underset{{\mathrm{n}}^{\mathrm{th}}\mathrm{spring}}{\underbrace{2G_n{E}_{kl}^{\left(n\right)}}}(1\le n\le N)\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill \text{single}\text{dashpot}⟹S_{kl}& =2{\eta }_{N+1}{D}^{{\delta }_{N+1}}{E}_{kl}^{(N+1)}\hfill \end{array}$$

(7)

where ${\delta }_n$ are distinct nonintegers and the derivatives $D^{{\delta }_n}=d^{{\delta }_n}/d{t}^{{\delta }_n}$ are defined in the Fourier sense [36]. For ${\delta }_n=1$ with $1\le n\le N+1$, these fractional derivative (FD) representations immediately reduce to the corresponding integer derivative (ID) ones. Alternatively, without changing any of the relations to follow, the FDs can also be represented in the Riemann–Liouville format [37] as

$$D^{{\delta }_n}\left\{F\left(t\right)\right\}={\displaystyle \frac{d^{{\delta }_n}F\left(t\right)}{d{t}^{{\delta }_n}}}={\displaystyle \frac{1}{\Gamma (1-{\delta }_n)}}{\displaystyle \frac{d}{dt}}\underset{0}{\overset{t}{\int }}{\displaystyle \frac{F\left(t^{\prime }\right)}{{(t-t^{\prime })}^{{\delta }_n}}}d{t}^{\prime }(0<{\delta }_n<1)$$

(8)

but with the limitations ${\delta }_n<1<n$ it is not possible to recover the corresponding integer derivatives $d^n/d{t}^n$. However, the original definitions in [38,39,40] carried no such restrictions, requiring only that $\Re \left({\delta }_n\right)>0$. In the present paper, the original more permissive condition for ${\delta }_n$ will be enforced.

The total strain deviator for the GFDKM and also for the GKM is the sum of the individual contributions, or

$$E_{kl}=\sum _{n=0}^{N+1}{E}_{kl}^{\left(n\right)}$$

(9)

while this sum may be formed by successive fractional differentiations of the individual parts [3], it is much more expeditious to apply FT to Equations (5)–(7), such that

$$2{\overline{\overline{E}}}_{kl}=\left[\frac{1}{G_0}+\frac{1}{{\eta }_{N+1}{\left(ı\omega \right)}^{{\delta }_{N+1}}}+\sum _{n=1}^N\frac{1}{{\eta }_n}\frac{1}{{\left(ı\omega \right)}^{{\delta }_n}+1/{\tau }_n}\right]{\overline{\overline{S}}}_{kl}=\frac{{\overline{\overline{S}}}_{kl}}{{\overline{\overline{G}}}_{FD}}$$

(10)

where $ı=\sqrt{-1}$ and the double bars indicate FT. From Equation (10), the complex shear modulus then can be defined as

$$\frac{1}{{\overline{\overline{G}}}_{FD}}=\frac{{\displaystyle \prod _{n=1}^{N+1}}\left[A_n{\left(ı\omega \right)}^{{\gamma }_n}+B_n\right]}{G_0{\left(ı\omega \right)}^{{\delta }_{N+1}}{\displaystyle \prod _{n=1}^N}\left[{\left(ı\omega \right)}^{{\delta }_n}+1/{\tau }_n\right]}={\overline{\overline{J}}}_{FD}$$

(11)

with ${\overline{\overline{J}}}_{FD}$ as the complex shear compliance. The complex moduli and compliances, of course, are at the heart of the isotropic and anisotropic elastic-viscoelastic integral transform analogies [3,35].

Substitution of Equation (11) into (10) and formally carrying out the multiplications results in

$$2G_0\sum _{n=0}^{N+1}{b}_n^{\prime }{\left(ı\omega \right)}^{{\beta }_n}{\overline{\overline{E}}}_{kl}=\sum _{n=0}^{N+1}{a}_n^{\prime }{\left(ı\omega \right)}^{{\alpha }_n}{\overline{\overline{S}}}_{kl}$$

(12)

and upon inverting the FT yields

$$2G_0{Q}_{FD}\left\{E_{kl}\right\}=P_{FD}\left\{E_{kl}\right\}$$

(13)

where

$$P_{FD}=\sum _{n=0}^r{a}_n^{\prime }{D}^{{\alpha }_n}\text{and}{Q}_{FD}=\sum _{n=0}^s{b}_n^{\prime }{D}^{{\beta }_n}$$

(14)

It should be noted that when $G_0<\infty$ then the single elastic spring is present in the system. (See [3] for other combinations and for the relations between r, s and N.) Some other observations are noteworthy, namely that for $G_0<\infty$ and $r=s=N+1$

$${\alpha }_N={\beta }_N=\sum _{n=1}^{N+1}{\delta }_n$$

(15)

and with the single dashpot absent, ${\delta }_{N+1}=0$ or ${\eta }_{N+1}=\infty$ or both. For ${\alpha }_n={\beta }_n=n$, the fractional derivative operators $P_{FD}$ and $Q_{FD}$ reduce to integer order derivative ones denoted by P and Q. On the other hand, the fractional derivative orders ${\alpha }_n$ are not necessarily equal to the corresponding ${\beta }_n$ for each n. One of the advantages of the fractional derivative models is that for the same number N of individual Kelvin or Maxwell bodies a third more material parameters, i.e., the ${\delta }_n$, are available for fitting actual creep or relaxation data, thus in principle, reducing the needed value of N for the same goodness of fit when compared to ID characterizations.

The FT of the creep and relaxation functions are given by

$$ı\omega {\overline{\overline{\varphi }}}_{FD}=\frac{1}{ı\omega {\overline{\overline{\psi }}}_{FD}}={\overline{\overline{G}}}_{FD}=\frac{1}{{\overline{\overline{J}}}_{FD}}$$

(16)

which relates the FT of the creep function ${\overline{\overline{\psi }}}_{FD}$ directly to the expression in square bracket in Equation (10). The FT inversion of the first two terms poses no difficulty, but the individual terms in the summation have no tabulated transform. Consequently, their inversion can be accomplished only by formal integration, fast Fourier transforms (FFT) [41], asymptotic expansions or some approximation. In [15] one of these for a single $\beta = .5$ and large times are presented. In [11] the approximation applicable for long times, i.e., small ωs,

$$\frac{1}{1+{\left(ı\omega \right)}^{\beta }\tau }\simeq 1-{\left(ı\omega \right)}^{\beta }\tau \mathrm{for}{\left(ı\omega \right)}^{\beta }\tau \ll 1\mathrm{and}t\gg 1$$

(17)

has been used. Alternatively, a different approximate inversion scheme is suggested here for each fraction in the sum of Equation (10) by noting the near similarity of the functions

$$\frac{1}{{\left(ı\omega \right)}^{{\delta }_n}+1/{\tau }_n}\simeq \frac{1}{{(ı\omega +c_n)}^{{\mu }_n}}⟹F\left(t\right)={\displaystyle \frac{e^{-c_{n}t}{t}^{{\mu }_n-1}}{\Gamma \left(c_n\right)}}$$

(18)

as seen in Figure 3 with $p=ı\omega$ for corresponding Laplace transforms. While the agreement is poorest in the range of small values of p, it must be remembered that this is the region of very long time responses where solutions may not be required and where asymptotic expressions can be used to greater benefit. Furthermore, while the approximation between individual terms may be less than desired, all parameters $c_n$ and ${\mu }_n$ can be adjusted by appropriate curve fitting of creep and relaxation data so that the sums match each other collectively with much greater accuracy than that of each individually paired term. This is in keeping with the practical requirements of actually modeling the behavior of real materials in terms of the entire complex modulus representations without the need to necessarily produce precise values for each of its constituent parameters. The advantage of the Equation (18) is that the approximation on the r. h. s. has a known inverse and an approximate analytical expression can then be obtained by inverting Equations (10) and (16) to yield

$${\psi }_{FD}\left(t\right)\simeq \underset{\text{initial spring}}{\underbrace{\frac{1}{G_0}}}+\underset{\text{final dashpot}}{\underbrace{\frac{t^{{\beta }_{N+1}+1}}{{\eta }_{N+1}\Gamma ({\beta }_{N+1}+1)}}}+\sum _{n=1}^N\underset{\text{Kelvin models}}{\underbrace{\left[\frac{1}{{\eta }_n\Gamma \left({\mu }_n\right)}\underset{0}{\overset{t}{\int }}{p}^{{\mu }_n-1}exp\left(-\frac{p}{{\tau }_n}\right)dp\right]}}$$

(19)

where the first two expressions are exact, but the last is approximate for ${\psi }_{FD}\left(t\right)$ even though it is the exact inverse of the r. h. s. of Equation (18). The exact expression for the Kelvin model integrals is given in Equation (44).

Figure 3. Comparison of exact and approximate relaxation moduli Laplace transforms of $1/(p^{\beta }+1/\tau )$ and $1/{(p+c)}^{\mu }$.

Figure 3. Comparison of exact and approximate relaxation moduli Laplace transforms of $1/(p^{\beta }+1/\tau )$ and $1/{(p+c)}^{\mu }$.

The isothermal integral fractional derivative constitutive relations are

$$2E_{kl}\left(t\right)=\underset{0}{\overset{t}{\int }}{\psi }_{FD}(t-t^{\prime })\frac{\partial S_{kl}\left(t^{\prime }\right)}{\partial t^{\prime }}d{t}^{\prime }=\underset{0}{\overset{t}{\int }}{J}_{FD}(t-t^{\prime })S_{kl}\left(t^{\prime }\right)d{t}^{\prime }$$

(20)

and they differ from the integer derivative ones only in the definition of the creep and relaxation functions, i.e., by replacing ${\psi }_{FD}$ with ψ in Equation (20).

Expressions similar to Equations (5) through Equation (20) can also be formulated for volumetric constitutive relations for σ and ϵ resulting in

$$ϵ\left(t\right)=\underset{0}{\overset{t}{\int }}{\psi }_{FD}^V(t-t^{\prime })\frac{\partial \sigma \left(t^{\prime }\right)}{\partial t^{\prime }}d{t}^{\prime }+\underset{0}{\overset{t}{\int }}{\psi }_{FD}^{VT}(t-t^{\prime })\frac{\partial \left[\mathcal{A}T\left(t^{\prime }\right)\right]}{\partial t^{\prime }}d{t}^{\prime }$$

(21)

where $\mathcal{A}$ is the coefficient of thermal expansion. Using Equations (3) and (4) and substituting Equation (21) into Equation (20) produces a single set of isotropic constitutive relations between ${ϵ}_{kl}$ and ${\sigma }_{kl}$ such that

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill 2{ϵ}_{kl}\left(t\right)=& \underset{0}{\overset{t}{\int }}{\psi }_{FD}(t-t^{\prime })\frac{\partial [{\sigma }_{kl}\left(t^{\prime }\right)-{\delta }_{kl}\sigma \left(t^{\prime }\right)]}{\partial t^{\prime }}d{t}^{\prime }\hfill \\ & +2{\delta }_{kl}\left\{\underset{0}{\overset{t}{\int }}{\psi }_{FD}^V(t-t^{\prime })\frac{\partial \sigma \left(t^{\prime }\right)}{\partial t^{\prime }}d{t}^{\prime }+\underset{0}{\overset{t}{\int }}{\psi }_{FD}^{VT}(t-t^{\prime })\frac{\partial \left[\mathcal{A}T\left(t^{\prime }\right)\right]}{\partial t^{\prime }}d{t}^{\prime }\right\}\hfill \end{array}\end{array}$$

(22)

As can be seen from these isotropic stress-strain relations, the normal strains are related only to normal stresses and shear strains are expressed in terms of only a single corresponding shear stress component, such that

$$\begin{array}{cccc}\hfill {ϵ}_{\underset{¯}{ii}}& =f_{\underset{¯}{ii}}[{\sigma }_{11},{\sigma }_{22},{\sigma }_{33},t]\hfill & \hfill \text{or}{\sigma }_{\underset{¯}{ii}}& =F_{\underset{¯}{ii}}[{ϵ}_{11},{ϵ}_{22},{ϵ}_{33},t]\hfill \end{array}$$

(23)

$$\begin{array}{cccc}\hfill {ϵ}_{ij}& =f_{\underset{¯}{ij}}[{\sigma }_{\underset{¯}{ij}},t]\hfill & \hfill \text{or}{\sigma }_{ij}& =F_{\underset{¯}{ij}}[{ϵ}_{\underset{¯}{ij}},t]i\ne j\hfill \end{array}$$

(24)

These isotropic constitutive relations—and in the next section their anisotropic cousins—are cast in the form of relaxation and/or creep functions, rather than combinations of the latter and viscoelastic Poisson’s ratios (PR). This is due to the fact viscoelastic PRs have been shown to be intrinsic functions of time [42,43,44], stresses and loading conditions and not “pure” material properties such as relaxation functions or moduli [45,46,47,48].

From a practical analysis point of view, the integral constitutive Equations (20) to (22) are always preferable to those expressed as differential operators such as Equation (14) even if the derivatives are of integer order. While the numerical evaluations of convolution integrals are CPU and memory intensive, they are by far more accurate than numerical solutions of high order PDEs. Furthermore, the discretized recurrence relations for such convolution integrals based on only two previous time steps developed in [49] are computationally economical and accurate.

For distinct temperature conditions viscoelastic materials admit the presence of temperature shift functions $a_T\left(T\right)$, which for thermorheologically simple materials (TSM) are material property functions of temperature, moisture, degree of cure and aging effects [34,50,51]. An associated reduced time ξ can then be defined as

$$\xi (x,t)=\underset{0}{\overset{t}{\int }}{a}_T(x,t^{\prime })d{t}^{\prime }$$

(25)

which transforms Equation (20) to a single relation in ξ instead of t for all service temperatures called the master relaxation or creep curve. The volumetric constitutive Equation (21) will have two additional reduced times, one for stresses or strains and the other for thermal expansions. Thus isotropic viscoelastic materials may have three distinct reduced times, although for convenience they are often approximated as all equal to each other, which makes characterization in terms of deviatoric and volumetric master curves possible. More complex nonisothermal behavior has been described in [24,52,53]. During the curing processes the material properties are also dependent on the degree of cure as well as on temperature and moisture content, thus further complicating the modeling of viscoelastic responses [54], but not conceptually altering any of the above results except that the shift function, as well as other parameters, become $a_T(T,M,\alpha )$, $\xi (T,M,\alpha )$, $G_n(T,M,\alpha )$ and ${\eta }_n(T,M,\alpha )$ where α is the degree of cure and M the moisture content, and with similar expressions for the volumetric components.

Nonlinear viscoelastic behavior necessitates that creep and relaxation functions are strain dependent, i.e., $\psi =\psi [x,\xi (x,t),I_1(x,\xi ),I_2(x,\xi ),I_3(x,\xi )]$, where $I_m$ are the first three invariants of the strain tensor ${ϵ}_{kl}$ [55]. The nonlinearities introduce additional complications into the analysis and make analytical solutions difficult if not impossible, necessitating numerical finite element formulations. While FT cannot be successfully applied to nonlinear constitutive relations, they can be used in the generalized linear formulation and then nonlinearities can be introduced a posteriori through proper interpretation when the entire linear modeling has been formulated and completed.

2.3. Anisotropic Relations

For anisotropic viscoelastic materials it is no longer possible to readily separate the stresses and strains into deviatoric and volumetric components, because volumetric stresses are now material property dependent [35]. Material parameters associated with mechanical models, such as $E_{ijkln}$ and ${\eta }_{ijkln}$ instead of the isotropic ones $E_n$ and ${\eta }_n$, must be identified with specific directions. Anisotropic states can then be achieved in three separate ways:

(1)

different parametric values $E_{ijkln}$ and ${\eta }_{ijkln}$ with equal numbers N in all directions

(2)

same parameters $E_{ijkln}$ and ${\eta }_{ijkln}$ in all directions but with distinct $N_{ijkl}$, thus generating different numbers of GKM parameters in each direction

(3)

combinations of (1) and (2) above

Inherent in such anisotropic behavior is the need to prescribe distinct integral and differential operator equations with property functions which are directionally dependent and which are associated with ${\varphi }_{kl}$, ${\psi }_{kl}$, $P_{kl}$, $Q_{kl}$, etc. In matrix form this can be written as

$$\left\{ϵ\right\}=\left[\mathbf{C}\right]\left\{\sigma \right\}+\left[{\mathbf{C}}_{\mathbf{T}}\right]\left\{\mathcal{A}\mathbf{T}\right\}$$

(26)

where the compliance matrices $\left[\mathbf{C}\right]$ and $\left[{\mathbf{C}}_{\mathbf{T}}\right]$ are either a set of integral or differential operators, resulting in FT expressions

$$\left\{\overline{\overline{ϵ}}\right\}=\left[\overline{\overline{\mathbf{C}}}\right]\left\{\overline{\overline{\sigma }}\right\}+\left[\overline{\overline{{\mathbf{C}}_{\mathbf{T}}}}\right]\left\{\overline{\overline{\mathcal{A}\mathbf{T}}}\right\}$$

(27)

Consequently, the anisotropic integral constitutive relations incorporating the above features and following the pattern of the isotropic relations of Equations (20) and (21) now read

$${ϵ}_{ij}(x,t)=\underset{-\infty }{\overset{t}{\int }}{\psi }_{ijkl}(x,t,t^{\prime })\frac{\partial {\sigma }_{kl}(x,t^{\prime })}{\partial t^{\prime }}d{t}^{\prime }+\underset{-\infty }{\overset{t}{\int }}{\psi }_{ijklT}(x,t,t^{\prime })\frac{\partial \left[{\mathcal{A}}_{kl}T\right](x,t^{\prime })}{\partial t^{\prime }}d{t}^{\prime }$$

(28)

with similar expressions for the fractional derivative ones obtained by changing ${\psi }_{ijkl}$ and ${\psi }_{ijklT}$ to ${\psi }_{ijkl}^{FD}$ and ${\psi }_{ijklT}^{FD}$ respectively. A modification of Equation (19) incorporating anisotropic representations then leads to

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\psi }_{ijkl}^{FD}\left(t\right)\widetilde{=}& \frac{1}{G_0}+\frac{t^{{\beta }_{N_{\underset{¯}{ijkl}}+1}+1}}{{\eta }_{N_{\underset{¯}{ijkl}}+1}\Gamma ({\beta }_{N_{\underset{¯}{ijkl}}+1}+1)}\hfill \\ & +\sum _{n=1}^{N_{\underset{¯}{ijkl}}}\left[\frac{1}{{\eta }_{\underset{¯}{ijkl}n}\Gamma \left({\mu }_{\underset{¯}{ijkl}n}\right)}\underset{0}{\overset{t}{\int }}{p}^{{\mu }_{\underset{¯}{ijkl}n}-1}exp\left(-\frac{p}{{\tau }_{\underset{¯}{ijkl}n}}\right)dp\right]\hfill \end{array}\end{array}$$

(29)

where the underscore indicates no summations on the marked indices.

Similarly, associated differential operators can be found by properly interpreting Equation (14) as

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {\widehat{P}}_{ij}(x,t)& =\sum _{n=0}^{r_{\underset{¯}{ij}}}{\widehat{a}}_{\underset{¯}{ij}n}(x,t)D^{{\alpha }_{\underset{¯}{ij}n}}\hfill \\ \hfill {\widehat{Q}}_{ijkl}(x,t)& =\sum _{n=0}^{s_{\underset{¯}{ijkl}}}{\widehat{b}}_{\underset{¯}{ijkl}n}(x,t)D^{{\beta }_{\underset{¯}{ijkl}n}}\hfill \\ \hfill {\widehat{Q}}_{ijklT}(x,t)& =\sum _{n=0}^{s_{\underset{¯}{ijkl}}^{\prime }}{\widehat{b}}_{\underset{¯}{ijkl}Tn}(x,t)D^{{\gamma }_{\underset{¯}{ijkl}n}}\hfill \end{array}\end{array}$$

(30)

where the symbol $\widehat{}$ is generic and refers to either integer or fractional derivative representations and for integer derivative characterization ${\alpha }_{ijkln}={\beta }_{ijkln}={\delta }_{ijkln}=n$. The anisotropic DE constitutive relations then become

$${\widehat{P}}_{\underset{¯}{ij}}\{{\sigma }_{\underset{¯}{ij}}\}={\widehat{Q}}_{ijkl}\left\{{ϵ}_{kl}\right\}-{\widehat{Q}}_{ijklT}\left\{{\mathcal{A}}_{kl}T\right\}$$

(31)

Temperature shift functions and reduced times must similarly be directionally identified as

$$\begin{array}{cc}\hfill {\xi }_{ijkl}(x,t)& =\underset{0}{\overset{t}{\int }}{a}_{Tijkl}[x,s,T(x,s),M(x,s),\alpha (x,s)]ds\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill {\xi }_{ijklT}(x,t)& =\underset{0}{\overset{t}{\int }}{a}_{TijklT}[x,s,T(x,s),M(x,s),{\alpha }_{cure}(x,s)]ds\hfill \end{array}$$

(33)

While in isotropic viscoelasticity there are only two ψs, one ${\psi }_T$, two ξs and one ${\xi }_T$, anisotropic viscoelasticity requires due to symmetry a maximum of 21 ${\psi }_{ijkl}$ and ${\xi }_{ijkl}$ each and six ${\psi }_{ijklT}$ and ${\xi }_{ijklT}$. The possible presence of even two or three reduced times in isotropic viscoelasticity, much less 27 for anisotropic responses, makes any solution by analytical means totally unmanageable. For effective analytical anisotropic viscoelastic formulations, it is, therefore, necessary to introduce only a single shift function and a single reduced time. This is analogous to what is done in isotropic viscoelasticity where the shift functions for changes in volume and shape and the one for thermal expansions are all equated to each other. Then

$$a_T=a_{Tijkl}=a_{TijklT}$$

(34)

and

$$\xi (x,t)={\xi }_{ijkl}={\xi }_{ijklT}$$

(35)

This simplifies the constitutive relations to

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill {ϵ}_{ij}(x,\xi )=& \underset{-\infty }{\overset{\xi }{\int }}{\widehat{\psi }}_{ijkl}(x,T_0,M_0,\xi -{\xi }^{\prime })\frac{\partial {\sigma }_{kl}(x,{\xi }^{\prime })}{\partial {\xi }^{\prime }}d{\xi }^{\prime }\hfill \\ & +\underset{-\infty }{\overset{\xi }{\int }}{\widehat{\psi }}_{ijklT}(x,T_0,M_0,\xi -{\xi }^{\prime })\frac{\partial \left[{\mathcal{A}}_{kl}T\right](x,{\xi }^{\prime })}{\partial {\xi }^{\prime }}d{\xi }^{\prime }\hfill \end{array}\end{array}$$

(36)

where $T_0$ and $M_0$ are the reference temperature and the moisture content respectively [35,49,54,55,56]. It should be noted that the shift function $a_T$ and $\widehat{\psi }$ itself and consequently the reduced time ξ may change their functional definitions drastically from one temperature and moisture region to another as the material undergoes phase changes, as is for instance the case during cure processes where the material property functions become additionally dependent on the degree of cure α [54, i.e., $\widehat{\psi }(x,T,M,\alpha ,\xi -{\xi }^{\prime })$ or $\widehat{\psi }(x,T,M,\alpha ,\xi ,{\xi }^{\prime })$.

Nonlinear responses may again be introduced into the characterization in the same manner as in the isotropic case, i.e.,

$${\psi }_{ijkl}={\psi }_{ijkl}\left[x,t,T(x,t),M(x,t),\alpha (x,t),I_1(x,t),I_2(x,t),I_3(x,t)\right]$$

(37)

Similar functional expressions can be generated for ${\psi }_{ijklT}$, ${\varphi }_{ijkl}$ and ${\varphi }_{ijklT}$. Such nonlinear representations, of course, turn the material property parameters ${\mu }_{ijkln}$, ${\eta }_{ijkln}$ and ${\tau }_{ijkln}$ of Equations (19) and (29) into functions of the variables contained in the preceding square bracket defining ψ. Specific isotropic examples have been formulated in [57] and extended to anisotropic viscoelasticity in [55].

Practical computational questions need to be raised about the FD characterization involving creep or relaxation functions. Neither the approximate Equations (19) and (29) nor the corresponding exact ones can be directly analytically integrated except by (asymptotic) series expansions for each ${\mu }_{ijkln}$, thereby increasing computational times for FEM solutions. It must be remembered that this solution process is inherently costly since these integrals or series must be repeatedly evaluated at every node point and at each time step in every direction.

One must, therefore, turn to an alternative but equally valid characterization for both isotropic and anisotropic viscoelastic materials in terms of compliances rather than creep functions. For anisotropic materials, Equation (16) can be rewritten as

$$ı\omega \overline{\overline{{\varphi }_{ijkl}^{FD}}}=\frac{1}{ı\omega \overline{\overline{{\psi }_{ijkl}^{FD}}}}=\overline{\overline{G_{ijkl}^{FD}}}=\frac{1}{\overline{\overline{J_{ijkl}^{FD}}}}$$

(38)

The constitutive Equation (28) in terms of compliances are

$${\widehat{ϵ}}_{ij}(x,\xi )=\underset{-\infty }{\overset{\xi }{\int }}{\widehat{J}}_{ijkl}(x,\xi -{\xi }^{\prime })\frac{\partial {\sigma }_{kl}(x,{\xi }^{\prime })}{\partial {\xi }^{\prime }}d{\xi }^{\prime }+\underset{-\infty }{\overset{\xi }{\int }}{\widehat{J}}_{ijklT}(x,\xi -{\xi }^{\prime })\frac{\partial \left[{\mathcal{A}}_{kl}T(x,{\xi }^{\prime })\right]}{\partial {\xi }^{\prime }}d{\xi }^{\prime }$$

(39)

with similar relations in the ξ space replacing Equations (14) and including both isotropic and anisotropic ID and FD representations. The inverse of the FD compliance $\overline{\overline{J_{ijkl}^{FD}}}$ for time dependent properties in the ξ space defined by Equations (10) and (18) is

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill J_{ijkl}^{FD}(x,\xi )\widetilde{=}& \frac{\delta \left(\xi \right)}{G_{0ijkl}(x,\xi )}+\frac{{\xi }^{{\beta }_{N_{\underset{¯}{ijkl}}+1}}}{{\eta }_{N_{\underset{¯}{ijkl}}+1}(x,\xi )\Gamma \left({\beta }_{N_{\underset{¯}{ijkl}}+1}\right)}\hfill \\ & +\sum _{n=1}^{N_{\underset{¯}{ijkl}}}\left\{\frac{{\xi }^{(-1+{\mu }_{\underset{¯}{ijkl}}n)}}{{\eta }_{\underset{¯}{ijkl}n}(x,\xi )\Gamma \left({\mu }_{\underset{¯}{ijkl}n}\right)}exp\left(-\underset{0}{\overset{\xi }{\int }}\frac{d{\xi }^{\prime }}{{\tau }_{\underset{¯}{ijkl}n}^{\prime }(x,{\xi }^{\prime })}\right)\right\}\hfill \end{array}\end{array}$$

(40)

where $1/{\tau }_{ijkln}^{\prime }=c_{ijkln}$ and the ${\tau }_{ijkln}^{\prime }$ do not necessarily represent actual relaxation times. For time independent nonhomogeneous anisotropic properties, the compliances reduce to

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill J_{ijkl}^{FD}(x,\xi )\widetilde{=}& \frac{\delta (\xi )}{G_{0ijkl}(x)}+\frac{{\xi }^{{\beta }_{N_{\underset{¯}{ijkl}}+1}}}{{\eta }_{N_{\underset{¯}{ijkl}}+1}(x)\Gamma ({\beta }_{N_{\underset{¯}{ijkl}}+1})}\hfill \\ & +\sum _{n=1}^{N_{\underset{¯}{ijkl}}}\left\{\frac{{\xi }^{(-1+{\mu }_{\underset{¯}{ijkl}n})}}{{\eta }_{\underset{¯}{ijkl}n}(x)\Gamma ({\mu }_{\underset{¯}{ijkl}n})}exp\left(\frac{-\xi }{{\tau }_{\underset{¯}{ijkl}n}^{\prime }(x)}\right)\right\}\hfill \end{array}\end{array}$$

(41)

The constitutive Equation (39) are considerably more tractable for analytical and computational analyses based on either isotropic or anisotropic materials, since they involve single time integrals, rather than the double ones needed for the ψ or φ relations. The use of moduli or compliances in the FD relations offers minimal computational penalties as it entails almost identical computational time and labor, except for the extra ${\xi }^{\mu -1}$ in the series terms, when compared to the integer derivative characterization.

The integer derivative constitutive Equation (18) are recoverable from Equations (19) and (20) by setting ${\beta }_{N_{ijkl}+1}={\mu }_{ijkln}=1$ for $1\le n\le N_{ijkl}$.

3. Discussion

Three of the important problem areas associated with viscoelastic constitutive relations and stress-deformation analyses are:

(a)

Performing “simple” experiments for which analytic solution can be formulated and evaluated

(b)

Curve fitting of actual creep and/or relaxation data by least square or other methods in order to determine modulus, creep function or compliance parameters

(c)

Inversion of FT or LT expressions for moduli, stresses and deformations

Curve fitting of moduli, etc. expressions can be carried out in either the real time space t or with FT in the frequency domain ω [1,2,3]. Values of the parameters $N_{ijkl}$ are selected a priori in order to achieve a preassigned accuracy of fit. For the ID representation there are only two sets of other parameters to be determined for each n ($1\le n\le N_{ijkl}$), namely ${\eta }_{ijkln}$ and ${\tau }_{ijkln}$, while for the FD relations a third set ${\mu }_{ijkln}$ or ${\beta }_{ijkln}$ must be added.

It can be readily seen from Equations (10), (19) and (41) that these expressions contain nonlinear functional dependencies on τs, βs and μs. (Changing the ηs to their inverses results in their linear contributions and, consequently, they present no computational problems.) For isotropic ID constitutive relations, Schapery [58] has suggested that an approximation be applied to the ${\tau }_n$ by setting each equal to ${10}^{2n}$ and then leaving only the $1/{\eta }_n$ to be determined by the fit. The same scheme is extended here to anisotropic ID relations. Experience shows that by proper selection of $N_{ijkl}$ any reasonable accuracy of fit can be attained, but for assumed τ values the Ns may be roughly twice as large as those corresponding to the correct τ evaluations, i.e., twice as many Prony series [59] terms. Under the Schapery approximation, the τs loose their physical meaning as relaxation times, but the modulus and/or creep functions maintain their accurate representations and physical meaning in toto, while the separate parameters do not. This represents no loss of generality in the physical definitions or actual evaluations by curve fit of experimental data of the entire modulus expressions [60,61,62,63].

Similarly, if one extends this idea to anisotropic FD relations, then approximate values can be generated for the exponents by setting $0<{\delta }_n\simeq |n\pm {\tilde{a}}_n|$ and $0<{\mu }_n\simeq {|n\pm \tilde{b}}_n|$ and where ${\tilde{a}}_n$ and ${\tilde{b}}_n>0$ can be adjusted for a “best” fit of the experimental data. This greatly simplifies the procedure for determining $1/{\eta }_n$ if additionally the same Schapery approximation is applied to the relaxation times, i.e., ${\tau }_n={\tau }_n^{\prime }=O\left({10}^{2n}\right)$.

The ID representations lead to integrable analytical expressions for ψ and φ and even though their values of the Ns may possibly be up to two to three times as large as the comparable ones for the FDs, the general computational time is considerably shorter. However, as was shown in the previous section, in order to avoid the second time integral needed for the ${\psi }_{ijkl}$ inversions in the FD forms, one must use the compliance form of the constitutive relations as shown in Equation (38). This still necessitates FT inversions to determine compliances as time functions. Such FT inversions can be accomplished by

• exact procedures, such as table look up or analytical evaluation of the ω FT inversion integrals
• numerical procedures, such as FFT [41]
• approximate procedures, such as those devised in [58], given by Equation (45) and limited to quasi-static problems where inertia and body forces are negligible
• exact for $k=\infty$, but approximate for $k<\infty$ LT inversions due to Post [64] as

$$f(x,t)=\underset{k\to \infty }{lim}{\left[\frac{{(-1)}^k}{k!}{p}^{k+1}\frac{{\partial }^k\overline{f}(x,p)}{\partial p^k}\right]}_{p=k/t}$$

(42)

where k are positive integers. Note that when the transform integrals exist, the relationship between LTs and FTs is given by

$$\overline{f}(x,p)=\underset{0}{\overset{\infty }{\int }}f(x,t)exp(-pt)dt={\left.\overline{\overline{f}}(x,\omega )\right|}_{p=ı\omega }={\left.\underset{\text{= half interval FT}}{\underbrace{\underset{0}{\overset{\infty }{\int }}f(x,t)exp(-ı\omega t)dt}}\right|}_{p=ı\omega }$$

(43)

Table look up and/or formal inversions through integration are, of course, the approaches of choice. Unfortunately, these work only for limited classes of LT and FT functions. In principle FFT always work but, as has been shown in [65] they are computationally intensive even on supercomputers because in order to obtain accurate FT inversion values one must sweep some 13 or more decades of ω, thus making the process computationally very costly. Finally, the approximate inversion scheme of Schapery based on producing inversions by multiplying the FT by $ı\omega$ and then setting $ı\omega =p=.5/t$ in the transform leads to simple procedures. Recalling that at the reference temperature $\xi =t$, the exact and approximate inverses of the r. h. s. of Equation (18) are given by

$$\begin{array}{cc}\hfill \text{exact}{F}_{1n}\left(\xi \right)& =\frac{{\xi }^{{\mu }_n-1}exp(-\xi /{\tau }_n^{\prime })}{\Gamma \left({\mu }_n\right)}\hfill \end{array}$$

(44)

$$\begin{array}{cc}\hfill \text{approximate}{F}_{2n}\left(\xi \right)& \simeq \frac{p}{{(p+1/{\tau }_n^{\prime })}^{{\mu }_n}}{|}_{p=.5/\xi }\hfill \end{array}$$

(45)

These functions together with the absolute error $F_{1n}\left(\xi \right)-F_{2n}\left(\xi \right)$ are displayed in Figure 4 . (The relative error cannot be shown because the exact solution tends to zero for larger values of ξ.) It is seen that the approximate inverse, while simple, is not very accurate. Equations (44) and (45) deal with only a single term of the Equation (40) representing the FD anisotropic compliances $J_{ijkl}^{FD}\left(\xi \right)$. The cumulative effect of such errors makes the Schapery approximate inversion scheme unacceptable in this instance.

Figure 4. Fractional derivative inverse.

Figure 4. Fractional derivative inverse.

Accuracy limitations of Schapery’s approximate inversion method when applied to dynamic isotropic viscoelastic problems has been previously pointed out in [66,67,68], among others. Although the method has found favor with much justification in quasi-static viscoelastic problems, this appears to be the first time that Schapery’s approximate inversion of a non-Prony series complex compliance itself is seen to also yield undesirable time functions.

In [69] the simple quasi-elastic and more involved self-consistent approximations for evaluating convolution integrals are compared. The reference’s conclusion is that the self-consistent method is more accurate, but yields more complicated results. The quasi-elastic method simply stated is

$${\sigma }_{ij}\left(\xi \right)=\underset{0}{\overset{\xi }{\int }}{E}_{ijkl}(\xi -{\xi }^{\prime }){ϵ}_{kl}\left({\xi }^{\prime }\right)d{\xi }^{\prime }\simeq E_{ijkl}\left(\xi \right){ϵ}_{kl}\left(\xi \right)$$

(46)

still requires knowledge of the inverses of ${\overline{\overline{E}}}_{ijkl}$ and ${\overline{\overline{ϵ}}}_{kl}$, but avoids the need for costly convolution integral evaluations. However, it is applicable only intermittently over relatively short time intervals and unfortunately, it is adequate only for simple classes of functions, such as single power laws, and cannot begin to cope with the complicated set of functions describing real material relaxation moduli [70]. The self-consistent approach leads to complicated results, even when tested against relatively simple functions, and is of limited practical use for real materials with the usual series expansions for their viscoelastic moduli.

More accurate, but also considerably more cumbersome, approximate inversion schemes are available [71]. Other pertinent works on approximate integral transform inversions may be found in [72,73,74,75].

Another possibility is the expansion of the complex moduli or compliances for both ID and FD characterizations in terms of a Lorenz series [3], such that

$${\overline{\overline{J}}}_{ijkl}\left(\omega \right)=\sum _{n=0}^{\infty }\left[\frac{A_{ijkln}}{{\left(ı\omega \right)}^n}+B_{ijkln}{\left(ı\omega \right)}^n\right]$$

(47)

where the $A_{ijkln}$ and $B_{ijkln}$ coefficients are determined by curve fitting experimental data. Even when truncated, the first sum gives excellent results for short times, $\omega \gg 1$, and the second one does the same for long times when $\omega \ll 1$. However, both sums must be used in the entire space $0<t<\infty$, except as noted in the relatively narrow bands at either time end, making this approach very cumbersome when combined with FEM for isotropic or anisotropic viscoelasticity.

The selection between FD and ID characterization then comes down to a choice between (a) fitting and using more complicated FD stress-strain laws with fewer individual terms but three parameters per term, and (b) simpler ID ones with two parameters per term but needing more terms in their series. In the final analysis, the ease with which the ID Prony series parameters can be generated from real material data in either the t or ω domains outweighs any advantages the FD representations may offer with its relatively shorter sums. A case in point is the work by [76] comparing the performance of a single ID Kelvin body with a 3- and 4-parameter FD model, showing marked improvements in accuracy of fit for FD models over ID ones for rubber and acrylic materials. However, as seen in Figure 6 , more complex real materials such as composites require ID characterizations with values for each $N_{ijkl}$ of approximately 30, i.e., about 61 Prony series parameters. For the corresponding ID and FD models in this example, the $N_{ijkl}$ are 29 and 10 with maximum errors of fit of 0.056 % and 10.15 % respectively. The $N_{ijkl}$ for FD representations can, of course, be increased but this defeats the primary purpose of their use. In the final analysis, the choice between ID and FD characterization has to be examined on a case by case basis for specific problem solutions, particularly those involving highly computationally intensive finite element approaches.

Figure 5. Comparison of ID and FT modeling of $1/{(p+c)}^{\mu }$.

Figure 5. Comparison of ID and FT modeling of $1/{(p+c)}^{\mu }$.

Figure 6. Laplace transform modulus.

Figure 6. Laplace transform modulus.

The requirement for relatively large $N_{ijkl}$ values for both ID and FD representations make the FD characterizations less attractive and modeling with the relatively simpler ID integral constitutive relations is less computationally intensive. Furthermore, for large scale engineering problems requiring upward of 10,000 FEA nodes, integral transforms are prohibitive and real time or reduced time ξ integrations need to be undertaken.

It should be further noted that while the area of highest percentage relative error is in the long time domain (small p values), the continuing mismatch between the actual and FD relaxation moduli spans over three initial time decades which are of primary interest. Because of the hereditary nature of the constitutive relation time integrals, these early time errors propagate throughout the entire time regime and adversely influence stress and displacement evaluations at all times.

Figure 6 represents comparisons for only one direction with normalized relaxation moduli and is typical of what is seen in other directions for anisotropic viscoelastic materials, except for adjustments in the fully relaxed modulus values and in the time scales.

In Figure 7 the approximate Kelvin model function of Equation (18) is displayed for a number of $\delta \le 1$ values. As can be readily seen these terms are sensitive to the exponent over the entire range of Laplace transform variable p and therefore correspondingly over the entire time range $-\infty \le t\le \infty$. Figure 8 exhibits the exact curve ${\left(p^{\delta }+1\right)}^{-1}$ and its approximate representation by ${\left(p+1\right)}^{-\mu }$ for specific, but representative, values of $\delta =0.75$ and $\mu =0.88$. The approximate representation mimics the exact one quite well with an error of of less than 0.1 over the entire time range. It must, of course, be remembered that experimental data for viscoelastic material properties is customarily in a comparable scatter range.

Figure 7. The influence of δ.

Figure 7. The influence of δ.

Figure 9 depicts the time functions associated with an ID and an approximate FD $(\delta =0.88)$ single Kelvin mode for $\tau =1$. It is readily seen that agreement between the two models is excellent, i.e., error $\cong 0$, except in the neighborhood of the origin $t=0$. This is due to the fact that for a single Kelvin model

$$\underset{t\to 0}{lim}\left\{exp\left(-{\displaystyle \frac{t}{\tau }}\right)\right\}\to 1\text{and}\underset{t\to 0}{lim}\left\{{\displaystyle \frac{exp(-t/c)t^{\mu -1}}{\Gamma \left(\mu \right)}}\right\}\to \infty 0<\mu <1$$

(48)

For materials characterized by a sum of approximate FD models,

$$\underset{t\to 0}{lim}\left\{\sum _{n=1}^N{J}_n{\displaystyle \frac{exp(-t/c_n)t^{{\mu }_n-1}}{\Gamma \left({\mu }_n\right)}}\right\}=\underset{t\to 0}{lim}\left\{\sum _{n=1}^N{J}_n{F}_n\left(t\right)\right\}\to \infty 0<{\mu }_n<1$$

(49)

the singularity still cannot be removed even through proper choices of the $J_n$, $c_n$ and ${\mu }_n$ are made to minimize the error of fit for ${\sum }_{n=1}^N{J}_n{F}_n\left(t\right)$ elsewhere at $0<t\le \infty$, since no value of $J_n\ne 0$ can correct or ameliorate the discrepancy at the time origin when for any one n ${\mu }_n<1$.

Figure 8. Exact and approximate FD.

Figure 8. Exact and approximate FD.

Figure 9. Comparison of ID and approximate FD.

Figure 9. Comparison of ID and approximate FD.

On the other hand, if the Riemann–Liouville condition of $\Re \left({\delta }_n\right)>0$ is enforced and ${\mu }_n>1$, then

$$\underset{t\to 0}{lim}\left\{J_n{F}_n\left(t\right)\right\}\to 01<n<N$$

(50)

For the Prony series representations in both real time t-space and reduced time ξ space, the constitutive relations can be transformed into ordinary instead convolution integrals, i.e.,

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \underset{-\infty }{\overset{\xi (x,t)}{\int }}exp\left(-\frac{\xi (x,t)-{\xi }^{\prime }(x,t^{\prime })}{{\tau }_{0n}}\right)& \widehat{f}\left[{\xi }^{\prime }(x,t^{\prime })\right]d{\xi }^{\prime }\to \hfill \\ & exp\left(-\frac{\xi (x,t)}{{\tau }_{0n}}\right)\underset{-\infty }{\overset{\xi (x,t)}{\int }}exp\left(\frac{{\xi }^{\prime }(x,t^{\prime })}{{\tau }_{0n}}\right)\widehat{f}\left[{\xi }^{\prime }(x,t^{\prime })\right]d{\xi }^{\prime }\hfill \end{array}\end{array}$$

(51)

so as to avoid CPU intensive computations [77]. Unfortunately, no such transformations have been derived for FD constitutive relations.

It is worthy to note that from a practical point of view either FD or ID are use to characterize the same material behavior and are obtained from the same experimental stress, strain and time data by fitting either or both appropriate analytical expressions to the latter. The accuracy of these fits depends on the number of terms used in each truncated series and on the scatter of the experimental data.

The least square fits used to determine parameters for the truncated FD and ID series lead to nonlinear transcendental relations in the first instance and to linear algebraic ones for the Prony series [60,61,62,63]. Comparisons of the goodness of fit for the two approaches can only be made only in a case by case approach for a specific material and a given set of experimental data.

4. Conclusions

While sufficiently accurate approximate procedures have been formulated to reasonably characterize real anisotropic materials with generalized integral and differential viscoelastic fractional derivative constitutive relations, their use with or without approximate characterization protocols requires functional complexity. Although integer derivative representations require possibly twice or three times as many terms in the expressions defining viscoelastic moduli, the relative simplicity of the resulting Prony series leads to computationally efficient anisotropic viscoelastic constitutive relations in real time as well as in FT frequency domains. In both fractional and integer derivative use, characterizations by integral rather than differential constitutive laws remain operationally preferred models because of their higher computational accuracies.