# Thermodynamical Restrictions and Wave Propagation for a Class of

Authors: Teodor M. Atanacković; Sanja Konjik; Ljubica Oparnica; Dušan Zorica
Journal: Abstract and Applied Analysis
 (2011)

---

## Abstract

We discuss thermodynamical restrictions for a linear constitutive equation containing fractional derivatives of stress and strain of different orders. Such an equation generalizes several known models. The restrictions on coefficients are derived by using entropy inequality for isothermal processes. In addition, we study waves in a rod of finite length modelled by a linear fractional constitutive equation. In particular, we examine stress relaxation and creep and compare results with the quasistatic analysis.


---

## Body

## 1. Introduction

Fractional calculus is intensively used to describe various phenomena in physics and engineering. We mention just few: diffusion and heat conduction processes [1, 2], damping in inelastic bodies [3], thermoelasticity [4], dissipative bending of rods [5], and concrete behavior of rods [6].Generalizations of classical equations of mathematical physics in order to include fractional derivatives, can be conducted in several ways. In the first approach, one changes the ordinary and partial integer order derivatives with the fractional ones in the relevant system of equations. Often, in this approach, the physical meaning of the newly defined terms with fractional derivatives might become unclear. In the second approach, one defines Lagrangian density function in which the integer order derivatives are replaced with fractional ones. In the next step, such a modified Lagrangian is used in the Hamilton principle (minimization of action integral) to obtain a system of the Euler-Lagrange equations, that describes the process (cf. [7, 8]). Although this variational approach have more sound physical meaning, it is not frequently used in applications since it, in principle, leads to equations with both left and right fractional derivatives. In both methods, the obtained equations have to satisfy requirements following from the Second Law of Thermodynamics, often expressed in the form of the Clausius-Duhamel inequality (see e.g., [9, 10]).In this work we propose a new model for linear viscoelastic body (obtained from the one-dimensional model) with an arbitrary number of springs and dashpots (see e.g., [11, page 28]) by replacing integer order derivatives with the Riemann-Liouville derivatives of real order (the first approach). Thus, we consider the class of constitutive equations of the form(1.1)∑n=0NanDtαn0σ=∑m=0MbmDtβm0ɛ,t>0.
Here σ denotes the Cauchy stress, ɛ the strain at the time instant t, while 0Dtη, η∈[0,1], denotes the operator of the left Riemann-Liouville fractional derivation. Recall, for η∈[0,1], the left ηth order Riemann-Liouville fractional derivative 0Dtηy of a function y∈AC([0,T]), T>0, is defined as (1.2)D0tηy(t)∶=1Γ(1-η)ddt∫0ty(τ)(t-τ)ηdτ,t>0,
where Γ is the Euler gamma function and AC([0,T]), T>0, denotes the space of absolutely continuous functions (for a detailed account on fractional calculus see e.g., [12]). For technical purposes, the orders of the fractional derivatives in (1.1) are assumed to satisfy  (1.3)0≤α0<α1<⋯<αn<⋯<αN≤1,0≤β0<β1<⋯<βm<⋯<βM≤1.
Also, in (1.1), the coefficients {an}n=0,…,N, {bm}m=0,…,M have the physical meaning of relaxation times and are assumed to be given. Many known constitutive equations of one-dimensional viscoelasticity are included as special cases of (1.1). For example, the choice α0=β0=0, α1=β1=1, a0=b0=1, a1=a and b1=b, and all other coefficients being equal to zero, leads to the classical Zener model (see [13]). If a0=b0=1, α0=β0=0, a1=a, b1=b, α1=β1=α, where α∈[0,1], and all other coefficients vanishing, then one has the generalized Zener model (see [9]). It is known that the Clausius-Duhamel inequality in both cases restricts constants a and b so that a≤b. However, we are not aware of restrictions on coefficients and orders of derivatives in generalized constitutive equation given by (1.1). Thus, in Section 2 of this work we shall derive these restrictions.In the second part of this work, we chose one specific equation of the type (1.1), proposed in [14]. Namely, the constitutive equation takes the form(1.4)(1+abDtα-β0)σ(x,t)=E(aDtα0+cDtγ0+acbDtα+γ-β0)ɛ(x,t),x∈[0,L],t>0,
where E is the generalized Young modulus (positive constant having dimension of stress), a, b, and c are given positive constants, while 0<β<α<γ<1/2. It should be stressed that in [14] there is no discussion concerning restrictions on the parameters a, b, c, α, β, and γ. Instead, only a,b,c≥0, 0≤β<α≤1 and 0≤γ≤1 were assumed. We will show, by the use of the method presented in [10], that conditions a,b,c≥0, and 0<β<α<γ<1/2 guarantee that the dissipation inequality is satisfied. The constitutive equation of the form (1.4) is obtained in [14] via the rheological model that generalizes the classical Zener rheological model by substituting spring and dashpot elements by fractional elements. It is assumed that the stress-strain relation for the fractional element is given by σ(t)=Dtηɛ0, t>0, where η∈[0,1]. We refer to [14] for more details of the derivation, creep compliance, and relaxation modulus of (1.4).The constitutive equation (1.4) actually describes a viscoelastic fluid body. For such a body of finite dimension we will analyze the wave propagation as well as two characteristic properties of the material: stress relaxation and creep. In this respect, the present work is a continuation of previous investigations (cf. [15–18]), where waves in viscoelastic solid body have been studied. We will consider the equation of motion of one-dimensional continuous body as(1.5)∂∂xσ(x,t)=ρ∂2∂t2u(x,t),x∈[0,L],t>0,
where ρ, σ, and u denote density, stress, and displacement of material at a point positioned at x and time t, respectively, as well as the strain measure, defined by(1.6)ɛ(x,t)=∂∂xu(x,t),x∈[0,L],t>0.
These two equations are coupled with a constitutive equation (1.4). We will impose initial and two types of boundary conditions to system (1.4), (1.5), and (1.6). The first type of boundary conditions describes the rod fixed at one end (displacement is zero during the time), while the other end is subject to a prescribed displacement. The second type describes the same rod but its other end is subject to the prescribed stress. We will obtain solutions in both cases in the convolution form of boundary conditions and the kernel of certain type. This setting is appropriate for examining stress relaxation and creep processes. If the displacement of the body's free end is prescribed as the Heaviside function then one can examine stress relaxation, while prescribing the stress of rod's free end by the Heaviside function will enable the study of creep.Similar problems have already been studied by several authors. System (1.5) and (1.6) coupled with the constitutive law of distributed-order type(1.7)∫01ϕσ(η)Dtησ(x,t)dη0=E∫01ϕɛ(η)Dtηɛ(x,t)dη0,x∈[0,L],t>0,
with ϕσ(η):=aη and ϕɛ(η):=bη, a≤b, was considered in [15] for the case of the stress relaxation in a viscoelastic material described by (1.7), and in [16] for the case of the creep and forced vibrations in a viscoelastic material of the same type. Special cases of (1.7) were studied in [19] (with ϕσ(η):=δ(η)+τεαδ(η-α) and ϕɛ(η):=E∞τεβδ(η-β), where δ is the Dirac distribution), and in [20] (with ϕσ(η):=δ(η)+τεαδ(η-α) and ϕɛ(η):=E0(δ(η)+τσαδ(η-α)+τσβδ(η-β)). Also our constitutive law (1.4) follows from (1.7) by choosing ϕσ(η):=δ(η)+(a/b)δ(η-(α-β)) and ϕɛ(η):=aδ(η-α)+cδ(η-γ)+(ac/b)δ(η-(α+γ-β)).We stress here that there is a strong connection between thermodynamical restrictions on coefficients in (1.1) and conditions for the existence of the inverse Laplace transform of equation of motion (see [17, 18]). It has been proved that the thermodynamical restrictions guarantee the existence of solutions.


## 2. Thermodynamical Restrictions

In this section we consider generalized linear fractional model and distributed-order fractional model of a viscoelastic body and give thermodynamical restrictions on coefficients and orders of fractional derivatives that appear in those models.


### 2.1. Generalized Linear Fractional Model

Our aim is to find restrictions on parameters of the model (1.1), that is, on {αn}n=0,…,N, {βm}m=0,…,M, {an}n=0,…,N, {bm}m=0,…,M, so that the generalized linear fractional model of a viscoelastic body satisfies the requirements of the Second Law of Thermodynamics.Following the procedure analogous to the one presented in [13], we apply the Fourier transform to (1.1) and obtain(2.1)σ̂(ω)∑n=0Nan(iω)αn=ɛ̂(ω)∑m=0Mbm(iω)βm,ω∈R,
where f̂(ω)=ℱ[f(t)](ω)=∫-∞∞f(t)e-iωtdt, ω∈ℝ, and ℱ[Dtηf0](ω)=(iω)ηf̂(ω). Writing (2.1) in the form (2.2)σ̂(ω)=Ê(ω)ɛ̂(ω),ω∈R,
with Ê being the complex modulus defined as (2.3)Ê(ω)∶=∑m=0Mbm(iω)βm∑n=0Nan(iω)αn,ω∈R,
one uses the conditions (cf. [10, 13])
(2.4)ReÊ(ω)≥0,∀ω>0,(2.5)Im⁡Ê(ω)≥0,∀ω>0,
which follow from the Second Law of Thermodynamics in case of the isothermal process, in order to obtain restrictions on parameters αn, βm, an, and bm.A straightforward calculation yields(2.6)Ê(ω)=∑m=0Mbmωβmcos⁡(βmπ/2)+i∑m=0Mbmωβmsin(βmπ/2)∑n=0Nanωαncos⁡(αnπ/2)+i∑n=0Nanωαnsin(αnπ/2).
Introducing Ê′ as (2.7)Ê′(ω)=Ê(ω)|∑n=0Nan(iω)αn|2,
one obtains(2.8)ReÊ′(ω)=(∑n=0Nanωαncos⁡αnπ2)(∑m=0Mbmωβmcos⁡βmπ2)+(∑n=0Nanωαnsinαnπ2)(∑m=0Mbmωβmsinβmπ2),ImÊ′(ω)=(∑n=0Nanωαncos⁡αnπ2)(∑m=0Mbmωβmsinβmπ2)-(∑n=0Nanωαnsinαnπ2)(∑m=0Mbmωβmcos⁡βmπ2).Sinceαn,βm∈[0,1], n=0,1,…,N, m=0,1,…,M, we have that αnπ/2,βmπ/2∈[0,π/2], and consequently, sine and cosine of those angles are positive. Therefore, assuming an,bm≥0 we obtain that ReÊ′(ω)≥0, and hence (2.4) is satisfied. In the sequel we will restrict our attention to this case (i.e., an,bm≥0).After a straightforward calculation of (2.8) one concludes that (2.5) holds if and only if  (2.9)ImÊ′(ω)=-∑n∈{0,1,…,N}m∈{0,1,…,M}ωαn+βmsin(αn-βm)π2anbm≥0,∀ω>0.Lemma 2.1.
Let (1.3) hold and an,bm≥0, n=0,1,…,N, m=0,1,…,M. Suppose that αN≠βM. Then a necessary condition for (2.9) is that αN<βM.
In other words, the highest order of fractional derivatives of stress in (1.1) could not be greater than the highest order of fractional derivatives of strain.Proof.
We observe that for largeω the sign of ImÊ′(ω) coincides with the sign of the term in the sum on the right-hand side of (2.9) with the largest power of ω. It follows from (1.3) that the latter is achieved for αN and βM. Therefore, in order that ImÊ′(ω)>0, αN has to be less than βM, as claimed.Remark 2.2.
(i) Similarly as in Lemma2.1 one can prove that if αN=βM then a necessary condition becomes αn<βm, for the largest αn and βm which do not coincide.
(ii) A similar conclusion for a particular problem with the constitutive equation of the form (1.1) has been obtained in [21].Example 2.3.
Equation of the forma0σ+aαDtασ0=b0ɛ, 0<α<1, cannot be a constitutive equation of a viscoelastic body. It does not obey the Second Law of Thermodynamics. Indeed, (2.9) reduces to -ωαsin(απ/2)aαb0, which is strictly less than zero for all ω>0.For practical purposes it can happen that the condition (2.9) is too general and thus hardly applicable to concrete problems. Therefore, it will be useful to extract from (2.9) particular conditions on parameters an,bm,αn and βm which guarantee that (2.5) is satisfied. We will separately consider several possibilities:(1)αn≠βm, for all n,m, that is, there are N+1 and M+1 terms of different order in (1.1).Since we assumed thatan,bm≥0, we can choose the orders of fractional derivatives as αn≤βm, foralln,m, so that all sin((αn-βm)π/2) are negative, and consequently ImÊ(ω)≥0. This, together with (1.3), further leads to the following condition:
(2.10)0≤α0<α1<⋯<αn<⋯<αN<β0<β1<⋯<βm<⋯<βM≤1.
In other words, an,bm≥0 and (2.10) are sufficient conditions for (2.5).Remark 2.4.
Note that the constitutive equation (1.4) belongs to this class. Indeed, by setting N=1, a0=1, α0=0, a1=a/b, α1=α-β, M=2, b0=Ea, β0=α, b1=Ec, β1=γ, b2=Eac/b, β2=α+γ-β in (1.1) we obtain (1.4). Thermodynamical restrictions are satisfied since constants E, a, b, c are positive, while (2.10) in this case reads
(2.11)0≤α-β<α<γ<α+γ-β≤1,
and is satisfied according to assumption 0<β<α<γ<1/2.(2)M>N and αi=βi, i=0,1,…,N, that is, there are N+1 first terms of the same order and M-N terms left in (1.1).Then (2.9) reduces to
(2.12)ImÊ′(ω)=∑i,j∈{0,1,…,N},i<jωαi+αjsin(αj-αi)π2(aibj-ajbi)-∑i∈{0,1,…,N},m>Nωαi+βmsin(αi-βm)π2aibm≥0∀ω>0.
We are interested in a special case when each term in the above sum is nonnegative. Then the following should hold:
(2.13)aibi≥ajbj,∀i,j∈{0,1,…,N}:i<j,αi≤βm,∀i∈{0,1,…,N},m>N,
which implies that
(2.14)a0b0≥a1b1≥⋯≥aNbN≥0,0≤α0<α1<⋯<αN<βN+1<⋯<βM≤1.(3)N>M and αN-M+i=βi, i=0,1,…,M, that is, there are M+1 last terms of the same order and N-M terms left in (1.1).Then (2.9) reduces to
(2.15)ImÊ′(ω)=-∑n<N-M,i∈{0,1,…,M}ωαn+αN-M+isin(αn-αN-M+i)π2anbi+∑i,j∈{0,1,…,M},i<jωβi+βjsin(βj-βi)π2(aN-M+ibj-aN-M+jbi)≥0,∀ω>0.
Again, we concentrate on a special case when each term in the above sum is nonnegative. By assumption (1.3) it follows that all terms in the first part of the sum are positive. Thus,
(2.16)aN-M+ibi≥aN-M+jbj,∀i,j∈{0,1,…,M}:i<j,which implies(2.17)aN-Mb0≥aN-M+1b1≥⋯≥aNbM≥0.(4)N=M and αi=βi, i=0,1,…,N, that is, there are N+1 terms of the same order on both sides of (1.1).In this case (2.9) becomes
(2.18)ImÊ′(ω)=∑i,j∈{0,1,…,N},i<jωαi+αjsin(αj-αi)π2(aibj-ajbi)≥0,∀ω>0.
As in the previous cases, we look for those parameters for which all terms in the above sum are nonnegative. Hence, we have to restrict our attention only to those parameters an,bm≥0 which in addition satisfy
(2.19)aibi≥ajbj,∀i,j:i<j.
This implies that
(2.20)a0b0≥a1b1≥⋯≥aNbN≥0.Remark 2.5.
It can be shown that in other cases which are not listed above (e.g., when there areK<min{N,M} terms of the same order in (1.1)) it is not possible to make all terms in (2.9) nonnegative simultaneously, thus it is more difficult to find some particular conditions on parameters an, bm, αn, and βm which implies (2.5).


### 2.2. Distributed-Order Fractional Model

The constitutive equation (1.1) can be generalized by (1.7), and further, if the integrals in (1.7) are interpreted in the distributional setting as
(2.21)〈∫supp⁡ϕϕ(α)Dtαhdα,φ(t)0〉∶=〈ϕ(α),〈D0tαh,φ(t)〉〉,φ∈D(R),
then (1.1) is generalized by
(2.22)〈ϕσ(α),〈D0tασ,φ(t)〉〉=〈ϕɛ(α),〈D0tαɛ,φ(t)〉〉,φ∈D(R)
(cf. [22]). Here ϕσ and ϕɛ are positive integrable functions of α, α∈[0,1], or distributions with compact support in [0,1], and φ denotes a test function belonging to the space 𝒟(ℝ) of compactly supported smooth functions on ℝ.Notice that by settingϕσ(α):=∑n=0Nanδ(α-αn), and ϕɛ(α):=∑m=0Mbmδ(α-βm) one obtains (1.1).In the sequel, we will considerϕσ(α):=aα and ϕɛ(α):=bα and look for the restrictions on a and b. Then constitutive equation (1.7) becomes ∫01aα0Dtασdα=∫01bα0Dtαɛdα, t>0. Interpreting these integrals as the Riemann sums, one obtains
(2.23)∑n=0NaαnDtαnσΔαn0=∑n=0NbαnDtαnɛΔαn0,
where N→∞ and Δαn→0. Putting an:=aαnΔαn and bn=bαnΔαn, one obtains (1.1) with N=M and all terms of the same order.Takinga,b>0 we have that (2.4) holds, while, as in the case (4) above (with equal number of terms of the same order on both sides of (1.1)), condition (2.24)(ab)αi≥(ab)αj,∀i,j:i<j,
provides validity of (2.5). Since αi≤αj, for all i,j, i<j, the latter is equivalent to
(2.25)ab<1,thatis,a<b.


## 2.1. Generalized Linear Fractional Model

Our aim is to find restrictions on parameters of the model (1.1), that is, on {αn}n=0,…,N, {βm}m=0,…,M, {an}n=0,…,N, {bm}m=0,…,M, so that the generalized linear fractional model of a viscoelastic body satisfies the requirements of the Second Law of Thermodynamics.Following the procedure analogous to the one presented in [13], we apply the Fourier transform to (1.1) and obtain(2.1)σ̂(ω)∑n=0Nan(iω)αn=ɛ̂(ω)∑m=0Mbm(iω)βm,ω∈R,
where f̂(ω)=ℱ[f(t)](ω)=∫-∞∞f(t)e-iωtdt, ω∈ℝ, and ℱ[Dtηf0](ω)=(iω)ηf̂(ω). Writing (2.1) in the form (2.2)σ̂(ω)=Ê(ω)ɛ̂(ω),ω∈R,
with Ê being the complex modulus defined as (2.3)Ê(ω)∶=∑m=0Mbm(iω)βm∑n=0Nan(iω)αn,ω∈R,
one uses the conditions (cf. [10, 13])
(2.4)ReÊ(ω)≥0,∀ω>0,(2.5)Im⁡Ê(ω)≥0,∀ω>0,
which follow from the Second Law of Thermodynamics in case of the isothermal process, in order to obtain restrictions on parameters αn, βm, an, and bm.A straightforward calculation yields(2.6)Ê(ω)=∑m=0Mbmωβmcos⁡(βmπ/2)+i∑m=0Mbmωβmsin(βmπ/2)∑n=0Nanωαncos⁡(αnπ/2)+i∑n=0Nanωαnsin(αnπ/2).
Introducing Ê′ as (2.7)Ê′(ω)=Ê(ω)|∑n=0Nan(iω)αn|2,
one obtains(2.8)ReÊ′(ω)=(∑n=0Nanωαncos⁡αnπ2)(∑m=0Mbmωβmcos⁡βmπ2)+(∑n=0Nanωαnsinαnπ2)(∑m=0Mbmωβmsinβmπ2),ImÊ′(ω)=(∑n=0Nanωαncos⁡αnπ2)(∑m=0Mbmωβmsinβmπ2)-(∑n=0Nanωαnsinαnπ2)(∑m=0Mbmωβmcos⁡βmπ2).Sinceαn,βm∈[0,1], n=0,1,…,N, m=0,1,…,M, we have that αnπ/2,βmπ/2∈[0,π/2], and consequently, sine and cosine of those angles are positive. Therefore, assuming an,bm≥0 we obtain that ReÊ′(ω)≥0, and hence (2.4) is satisfied. In the sequel we will restrict our attention to this case (i.e., an,bm≥0).After a straightforward calculation of (2.8) one concludes that (2.5) holds if and only if  (2.9)ImÊ′(ω)=-∑n∈{0,1,…,N}m∈{0,1,…,M}ωαn+βmsin(αn-βm)π2anbm≥0,∀ω>0.Lemma 2.1.
Let (1.3) hold and an,bm≥0, n=0,1,…,N, m=0,1,…,M. Suppose that αN≠βM. Then a necessary condition for (2.9) is that αN<βM.
In other words, the highest order of fractional derivatives of stress in (1.1) could not be greater than the highest order of fractional derivatives of strain.Proof.
We observe that for largeω the sign of ImÊ′(ω) coincides with the sign of the term in the sum on the right-hand side of (2.9) with the largest power of ω. It follows from (1.3) that the latter is achieved for αN and βM. Therefore, in order that ImÊ′(ω)>0, αN has to be less than βM, as claimed.Remark 2.2.
(i) Similarly as in Lemma2.1 one can prove that if αN=βM then a necessary condition becomes αn<βm, for the largest αn and βm which do not coincide.
(ii) A similar conclusion for a particular problem with the constitutive equation of the form (1.1) has been obtained in [21].Example 2.3.
Equation of the forma0σ+aαDtασ0=b0ɛ, 0<α<1, cannot be a constitutive equation of a viscoelastic body. It does not obey the Second Law of Thermodynamics. Indeed, (2.9) reduces to -ωαsin(απ/2)aαb0, which is strictly less than zero for all ω>0.For practical purposes it can happen that the condition (2.9) is too general and thus hardly applicable to concrete problems. Therefore, it will be useful to extract from (2.9) particular conditions on parameters an,bm,αn and βm which guarantee that (2.5) is satisfied. We will separately consider several possibilities:(1)αn≠βm, for all n,m, that is, there are N+1 and M+1 terms of different order in (1.1).Since we assumed thatan,bm≥0, we can choose the orders of fractional derivatives as αn≤βm, foralln,m, so that all sin((αn-βm)π/2) are negative, and consequently ImÊ(ω)≥0. This, together with (1.3), further leads to the following condition:
(2.10)0≤α0<α1<⋯<αn<⋯<αN<β0<β1<⋯<βm<⋯<βM≤1.
In other words, an,bm≥0 and (2.10) are sufficient conditions for (2.5).Remark 2.4.
Note that the constitutive equation (1.4) belongs to this class. Indeed, by setting N=1, a0=1, α0=0, a1=a/b, α1=α-β, M=2, b0=Ea, β0=α, b1=Ec, β1=γ, b2=Eac/b, β2=α+γ-β in (1.1) we obtain (1.4). Thermodynamical restrictions are satisfied since constants E, a, b, c are positive, while (2.10) in this case reads
(2.11)0≤α-β<α<γ<α+γ-β≤1,
and is satisfied according to assumption 0<β<α<γ<1/2.(2)M>N and αi=βi, i=0,1,…,N, that is, there are N+1 first terms of the same order and M-N terms left in (1.1).Then (2.9) reduces to
(2.12)ImÊ′(ω)=∑i,j∈{0,1,…,N},i<jωαi+αjsin(αj-αi)π2(aibj-ajbi)-∑i∈{0,1,…,N},m>Nωαi+βmsin(αi-βm)π2aibm≥0∀ω>0.
We are interested in a special case when each term in the above sum is nonnegative. Then the following should hold:
(2.13)aibi≥ajbj,∀i,j∈{0,1,…,N}:i<j,αi≤βm,∀i∈{0,1,…,N},m>N,
which implies that
(2.14)a0b0≥a1b1≥⋯≥aNbN≥0,0≤α0<α1<⋯<αN<βN+1<⋯<βM≤1.(3)N>M and αN-M+i=βi, i=0,1,…,M, that is, there are M+1 last terms of the same order and N-M terms left in (1.1).Then (2.9) reduces to
(2.15)ImÊ′(ω)=-∑n<N-M,i∈{0,1,…,M}ωαn+αN-M+isin(αn-αN-M+i)π2anbi+∑i,j∈{0,1,…,M},i<jωβi+βjsin(βj-βi)π2(aN-M+ibj-aN-M+jbi)≥0,∀ω>0.
Again, we concentrate on a special case when each term in the above sum is nonnegative. By assumption (1.3) it follows that all terms in the first part of the sum are positive. Thus,
(2.16)aN-M+ibi≥aN-M+jbj,∀i,j∈{0,1,…,M}:i<j,which implies(2.17)aN-Mb0≥aN-M+1b1≥⋯≥aNbM≥0.(4)N=M and αi=βi, i=0,1,…,N, that is, there are N+1 terms of the same order on both sides of (1.1).In this case (2.9) becomes
(2.18)ImÊ′(ω)=∑i,j∈{0,1,…,N},i<jωαi+αjsin(αj-αi)π2(aibj-ajbi)≥0,∀ω>0.
As in the previous cases, we look for those parameters for which all terms in the above sum are nonnegative. Hence, we have to restrict our attention only to those parameters an,bm≥0 which in addition satisfy
(2.19)aibi≥ajbj,∀i,j:i<j.
This implies that
(2.20)a0b0≥a1b1≥⋯≥aNbN≥0.Remark 2.5.
It can be shown that in other cases which are not listed above (e.g., when there areK<min{N,M} terms of the same order in (1.1)) it is not possible to make all terms in (2.9) nonnegative simultaneously, thus it is more difficult to find some particular conditions on parameters an, bm, αn, and βm which implies (2.5).


## 2.2. Distributed-Order Fractional Model

The constitutive equation (1.1) can be generalized by (1.7), and further, if the integrals in (1.7) are interpreted in the distributional setting as
(2.21)〈∫supp⁡ϕϕ(α)Dtαhdα,φ(t)0〉∶=〈ϕ(α),〈D0tαh,φ(t)〉〉,φ∈D(R),
then (1.1) is generalized by
(2.22)〈ϕσ(α),〈D0tασ,φ(t)〉〉=〈ϕɛ(α),〈D0tαɛ,φ(t)〉〉,φ∈D(R)
(cf. [22]). Here ϕσ and ϕɛ are positive integrable functions of α, α∈[0,1], or distributions with compact support in [0,1], and φ denotes a test function belonging to the space 𝒟(ℝ) of compactly supported smooth functions on ℝ.Notice that by settingϕσ(α):=∑n=0Nanδ(α-αn), and ϕɛ(α):=∑m=0Mbmδ(α-βm) one obtains (1.1).In the sequel, we will considerϕσ(α):=aα and ϕɛ(α):=bα and look for the restrictions on a and b. Then constitutive equation (1.7) becomes ∫01aα0Dtασdα=∫01bα0Dtαɛdα, t>0. Interpreting these integrals as the Riemann sums, one obtains
(2.23)∑n=0NaαnDtαnσΔαn0=∑n=0NbαnDtαnɛΔαn0,
where N→∞ and Δαn→0. Putting an:=aαnΔαn and bn=bαnΔαn, one obtains (1.1) with N=M and all terms of the same order.Takinga,b>0 we have that (2.4) holds, while, as in the case (4) above (with equal number of terms of the same order on both sides of (1.1)), condition (2.24)(ab)αi≥(ab)αj,∀i,j:i<j,
provides validity of (2.5). Since αi≤αj, for all i,j, i<j, the latter is equivalent to
(2.25)ab<1,thatis,a<b.


## 3. A Model of Viscoelastic Body of Finite Length

In this section we analyze wave propagation in a viscoelastic rod of finite length. The rod is made of a viscoelastic material described by a fractional-type constitutive equation (1.4), which is a special case of generalized linear fractional model (1.1). In fact, we will study an initial boundary value problem for system (1.4), (1.5), and (1.6).


### 3.1. Convolution Form of Solutions

Consider system (1.4), (1.5), and (1.6) supplied with initial conditions
(3.1)u(x,0)=0,∂∂tu(x,0)=0,σ(x,0)=0,ɛ(x,0)=0,x∈[0,L],
as well as with two types of boundary conditions:(3.2)u(0,t)=0,u(L,t)=Υ(t),t∈R,u(0,t)=0,σ(L,t)=Σ(t),t∈R.
Functions Υ and Σ are locally integrable functions supported in [0,∞). Note that if Υ=Υ0H we have the case of stress relaxation, while if Σ=Σ0H we have the case of creep, where H is the Heaviside function.Introducing dimensionless quantities(3.3)x¯=xL,t¯=tLρ/E,u¯=uL,σ¯=σE,Υ¯=ΥL,Σ¯=ΣE,a¯=a(Lρ/E)α,b¯=b(Lρ/E)β,c¯=c(Lρ/E)γ,
and using the fact that fractional derivatives are transformed as 0Dt¯ηu(t¯)=(Lρ/E)η0Dtηu(t) (cf. [17, 23]), after omitting the bar over dimensionless quantities, we obtain the following system: x∈[0,1], t>0,
(3.4)∂∂xσ(x,t)=∂2∂t2u(x,t),ɛ(x,t)=∂∂xu(x,t),(1+abDtα-β0)σ(x,t)=E(aDtα0+cDtγ0+acbDtα+γ-β0)ɛ(x,t).
System (3.4) is subject to initial(3.5)u(x,0)=0,∂∂tu(x,0)=0,σ(x,0)=0,ɛ(x,0)=0,x∈[0,1],
and two types of boundary conditions
(3.6)u(0,t)=0,u(1,t)=Υ(t),t∈R,(3.7)u(0,t)=0,σ(1,t)=Σ(t),t∈R
as described above.Solutions of the above system will be determined by the Laplace transform method. Recall, the Laplace transform off∈Lloc1(ℝ), f≡0 in (-∞,0), and |f(t)|≤Meat, for some M>0, a∈ℝ and t≥0, is defined by (3.8)f̃(s)=L[f(t)](s)∶=∫0∞e-stf(t)dt,Res>a.
Thus, applying the Laplace transform to (3.4) and (3.5), one obtains  (3.9)∂∂xσ̃(x,s)=s2ũ(x,s),ɛ̃(x,s)=∂∂xũ(x,s),(1+absα-β)σ̃(x,s)=(asα+csγ+acbsα+γ-β)ɛ̃(x,s).
System (3.9) reduces to (3.10)∂2∂x2ũ(x,s)-(sM(s))2ũ(x,s)=0,x∈[0,1],s∈D,
where(3.11)M(s)∶=1+(a/b)sα-βasα+csγ+(ac/b)sα+γ-β=1asα1+(a/b)sα-β1+(c/a)sγ-α+(c/b)sγ-β.
It has a solution(3.12)ũ(x,s)=C1(s)exsM(s)+C2(s)e-xsM(s),x∈[0,1],s∈D,
where C1 and C2 are functions of s which will be determined from the boundary conditions. Since the power function sη is analytic on the complex plane except the branch cut along the negative axis (including the origin) we take D:=ℂ∖(-∞,0] to be the domain for variable s in (3.9). Applying either (3.6)1 or (3.7)1, we obtain C1=-C2=:C, and thus(3.13)ũ(x,s)=C(s)(exsM(s)-e-xsM(s)),x∈[0,1],s∈C∖(-∞,0].
From (3.9) and (3.11) it follows that(3.14)σ̃(x,s)=1M2(s)∂∂xũ(x,s),x∈[0,1],s∈C∖(-∞,0].As announced in Section1 we will separately seek solutions in different cases: displacement u and stress σ in the case of prescribed displacement (such as e.g., stress relaxation), and displacement u in the case of prescribed stress (e.g., creep). For the former we supply to the system boundary conditions (3.6), while for the latter we assume (3.7). In the sequel we derive convolution forms of solutions in all these cases.In the case of prescribed displacementΥ, substituting (3.13) into (3.12) and using (3.6), one obtains(3.15)ũ(x,s)=Υ̃(s)P̃(x,s),x∈[0,1],s∈C∖(-∞,0],
where(3.16)P̃(x,s)∶=sinh⁡(xsM(s))sinh⁡(sM(s)),x∈[0,1],s∈C∖(-∞,0].
Clearly, P(1,t)=δ(t), t∈ℝ.SinceΥ and P are supported in [0,∞), displacement u is given by  (3.17)u(x,t)=Υ(t)*P(x,t),x∈[0,1],t∈R,(3.18)u(x,t)=0,x∈[0,1],t<0,
where * denotes the convolution with respect to t. Recall, if f,g∈Lloc1(ℝ), suppf,g⊂[0,∞), then (f*g)(t):=∫0tf(τ)g(t-τ)dτ, t∈ℝ. Explicit calculation of (3.17) will be done by the use of the Laplace inversion formula applied to (3.15).Further, from (3.14), (3.15), and (3.16) it follows that(3.19)σ̃(x,s)=sΥ̃(s)T̃(x,s),x∈[0,1],s∈C∖(-∞,0],
where(3.20)T̃(x,s)=cosh⁡(xsM(s))M(s)sinh⁡(sM(s)),x∈[0,1],s∈C∖(-∞,0].
Applying the Laplace inversion formula to (3.19) we obtain(3.21)σ(x,t)=ddt(Υ(t)*T(x,t)),x∈[0,1],t∈R,
where the derivative is understood in the sense of distributions. Again, σ(x,t)=0 for x∈[0,1], t<0.In the case of prescribed stressΣ, using (3.14) at x=1 and (3.7), we obtain(3.22)∂∂xũ(1,s)=Σ̃(s)M2(s),s∈C∖(-∞,0].
This combined with (3.13) gives(3.23)ũ(x,s)=Σ̃(s)Q̃(x,s),x∈[0,1],s∈C∖(-∞,0],
where(3.24)Q̃(x,s)=1sM(s)sinh⁡(xsM(s))cosh⁡(sM(s)),x∈[0,1],s∈C∖(-∞,0].
Again, applying the Laplace inversion formula to (3.23), the displacement reads  (3.25)u(x,t)=Σ(t)*Q(x,t),x∈[0,1],t>0,(3.26)u(x,t)=0,x∈[0,1],t<0.Remark 3.1.
As in [15] one can show that P, T, and Q are real-valued, locally integrable functions on ℝ, supported in [0,∞) and smooth for t>0.


### 3.2. Explicit Forms of Solutions

This subsection is devoted to the calculation of inverse Laplace transforms of certain distributions and functions introduced above. We begin with examining some basic properties ofM given in (3.11), which will be important for further investigations. In the sequel we shall write A(x)~B(x) if limx→∞(A(x)/B(x))=1.Proposition 3.2.
LetM be the function defined by (3.11). Then  (i)
M is an analytic function on ℂ∖(-∞,0] if  0<β<α<γ<1.(ii)
Fors∈ℂ∖(-∞,0], lim|s|→0M(s)=∞, lim|s|→0sM(s)=0, lim|s|→∞M(s)=0, and lim|s|→∞sM(s)=∞.Proof.
(i) Since1+(a/b)sα-β≠0 and 1+(c/a)sγ-α+(c/b)sγ-β≠0 if args∈(-π,π) and 0<β<α<γ<1, it follows that the function M given in (3.11) is analytic on the complex plane except the branch cut along the negative axis, that is, on ℂ∖(-∞,0].
(ii) Limits in (ii) can easily be calculated.


#### 3.2.1. Determination of Displacementu in the Case of Prescribed Displacement Υ

To begin with, we examine properties ofP̃ given by (3.16). P̃ has isolated singularities at Psn(±), n∈ℕ, where Psn(±) denotes solutions of the equation
(3.27)sinh⁡(sM(s))=0,that is,sM(s)=±niπ.Let us examine their position and multiplicity.Proposition 3.3.
There are infinitely many complex conjugated solutionsPsn(±), n∈ℕ, of (3.27), which all lie in the left complex half plane. Moreover, each Psn(±), n∈ℕ, is a simple pole.Proof.
Let us square (3.27) and define
(3.28)Φ(s,n)∶=(sM(s))2+(nπ)2=s21+(a/b)sα-βasα+csγ+(ac/b)sα+γ-β+(nπ)2.
Writing s=Reiφ it follows that
(3.29)Φ(Reiφ,n)=R2e2iφ1+(a/b)Rα-βei(α-β)φaRαeiαφ+cRγeiγφ+(ac/b)Rα+γ-βei(α+γ-β)φ+(nπ)2=R2(cos⁡(2φ)+isin(2φ))A+iBC+iD+(nπ)2,
where
(3.30)A∶=1+abRα-βcos⁡((α-β)φ),B∶=abRα-βsin((α-β)φ),C∶=aRαcos⁡(αφ)+cRγcos⁡(γφ)+acbRα+γ-βcos⁡((α+γ-β)φ),D∶=aRαsin(αφ)+cRγsin(γφ)+acbRα+γ-βsin((α+γ-β)φ).
Real and imaginary parts of Φ are then given by
(3.31)ReΦ(Reiφ,n)=R2C2+D2[(AC+BD)cos⁡(2φ)+(AD-BC)sin(2φ)]+(nπ)2,ImΦ(Reiφ,n)=R2C2+D2[(AC+BD)sin(2φ)-(AD-BC)cos(2φ)],(3.32)AC+BD=aRαcos⁡(αφ)+cRγcos⁡(γφ)+a2bR2α-βcos⁡(βφ)+a2cb2R2(α-β)+γcos⁡(γφ)+2acbRα+γ-βcos⁡((α-β)φ)cos⁡(γφ),AD-BC=aRαsin(αφ)+cRγsin(γφ)+a2bR2α-βsin(βφ)+a2cb2R2(α-β)+γsin(γφ)+2acbRα+γ-βcos⁡((α-β)φ)sin(γφ).
Let(R,φ) be a solution to Φ(s,n)=0 (or equivalently, ReΦ=ImΦ=0). Then changing φ→-φ, we again obtain that ReΦ=ImΦ=0, which implies that solutions of (3.27) are complex conjugated.
Further, we will prove that the functionΦ has no zeros in the half plane args∈[0,π/2]. For that purpose we will use the argument principle. Recall, if Φ is an analytic function inside and on a regular closed curve c and nonzero on c, then the number of zeros of Φ is given by NZ=(1/2π)ΔargΦ(s). Let γ=γa∪γb∪γc be parametrized as
(3.33)γa:s=x,x∈[0,R],γb:s=Reiφ,φ∈[0,π2],γc:s=xeiπ/2,x∈[0,R],
and let R→∞. Along γa function Φ becomes a real-valued function, hence ΔargΦ(s,n)=0. Along γb we have that AC+BD,AD-BC≥0 for φ∈[0,π], since sin(ηφ),cos(ηφ)>0, for η∈{α,β,γ,α-β}⊂(0,1/2), and φ∈[0,π]. Therefore, (3.31) imply that
(3.34)ReΦ(Reiφ,n)>0,φ∈[0,π4],ImΦ(Reiφ,n)>0,φ∈[π4,π2]
(inequalities at boundary points are easily checked by inserting them into (3.31)). Along γc, we obtain
(3.35)ReΦ(xeiπ/2,n)=-R2C2+D2|R=x,φ=π/2(AC+BD)|R=x,φ=π/2+(nπ)2,(3.36)ImΦ(xeiπ/2,n)=R2C2+D2|R=x,φ=π/2(AD-BC)|R=x,φ=π/2>0.
From (3.34) and (3.36), we may now conclude that along γb and γc,
(3.37)ΔargΦ(s,n)=0.
Indeed, this follows from the following. Along γb, in the case φ∈[0,π/4], ImΦ(Reiφ,n) can change its sign, but ReΦ(Reiφ,n)>0, while for φ∈[π/4,π/2], ImΦ(Reiφ,n)>0. Along γc, ImΦ(xeiπ/2,n)>0. Therefore, there is no change in the argument of Φ along the whole γ, which implies that Φ has no zeros for φ∈[0,π/2]. This further implies that (3.27) has no solutions in the right complex half plane, since its solutions are complex conjugated.
In order to prove that for fixedn∈ℕ (3.27) has one solution (and its complex conjugate), we again use function Φ and the argument principle. Consider now the contour Γ=ΓA∪ΓB∪ΓC parametrized by
(3.38)ΓA:s=xeiπ/2,x∈[0,R],ΓB:s=Reiφ,φ∈[π2,π],ΓC:s=xeiπ,x∈[0,R],
and let R→∞. Along ΓA real and imaginary parts of Φ are given by (3.35) and (3.36). Along ΓB, using (3.31), we conclude that
(3.39)ReΦ(Reiφ,n)<0,φ∈[π2,3π4],forR⟶∞,andnfixed,ImΦ(Reiφ,n)<0,φ∈[3π4,π].
Along ΓC we have
(3.40)ReΦ(Reiφ,n)=R2C2+D2|R=x,φ=π(AC+BD)|R=x,φ=π+(nπ)2>0,ImΦ(Reiφ,n)=-R2C2+D2|R=x,φ=π(AD-BC)|R=x,φ=π<0.
Along ΓA we have that ImΦ>0, while ReΦ changes its sign (since ReΦ(0,n)=(nπ)2 and limx→∞ReΦ(xeiπ/2,n)=-∞, for fixed n∈ℕ). Along the part of ΓB where φ∈[3π/4,π], and along ΓC, ImΦ<0. Also, limR→∞ReΦ(Rei3π/4,n)=-∞ and ImΦ(Rei3π/4,n)<0. This implies that the argument of Φ changes from 0 to 2π.
As a conclusion one obtains that alongΓ(3.41)ΔargΦ(s,n)=2π,
which implies, by the argument principle, that function Φ has exactly one zero in the upper left complex plane, for each fixed n∈ℕ. Since the zeros of Φ are complex conjugated, it follows that Φ also has one zero in the lower left complex plane, for each fixed n∈ℕ.In the next proposition we examine behavior of simple polesPsn(±), n∈ℕ.Proposition 3.4.
SolutionsPsn(±), n∈ℕ, of (3.27), are such that
(3.42)Re(Psn(±))=Rcos⁡φ~c(nπ)22-γcos⁡(π2-γ)<0,Im(Psn(±))=Rsinφ~±c(nπ)22-γsin(π2-γ),
as n→∞.Proof.
Let us square (3.27) and insert Psn(±)=Reiφ, φ∈(-π,π). Then, after separation of real and imaginary parts, we obtain
(3.43)R2cos⁡(2φ)Re(M2(Reiφ))-R2sin(2φ)Im⁡(M2(Reiφ))=-(nπ)2,(3.44)R2sin(2φ)Re(M2(Reiφ))+R2cos⁡(2φ)Im⁡(M2(Reiφ))=0.
Using notation from the proof of Proposition 3.3 we can write
(3.45)Re(M2(Reiφ))=AC+BDC2+D2,Im⁡(M2(Reiφ))=BC-ADC2+D2.
Letting R→∞, one has
(3.46)Re(M2(Reiφ))~(a2c/b2)R2(α-β)+γcos⁡(γφ)(a2c2/b2)R2(α-β)+2γ=1cRγcos⁡(γφ),Im⁡(M2(Reiφ))~-(a2c/b2)R2(α-β)+γsin(γφ)(a2c2/b2)R2(α-β)+2γ=-1cRγsin(γφ).
It now follows from (3.44) and (3.46) that
(3.47)tg(2φ)=-Im(M2(Reiφ))Re(M2(Reiφ))~tg(γφ)⟹sin((2-γ)φ)cos⁡(2φ)cos⁡(γφ)~0⟹φ~±π2-γ.
Inserting (3.47) into (3.46), and subsequently into (3.43), we obtain
(3.48)R2-γccos⁡(2π2-γ)cos⁡(γπ2-γ)+R2-γcsin(2π2-γ)sin(γπ2-γ)~-(nπ)2,R~c(nπ)22-γ.
Thus, real and imaginary parts of Psn(±), as R→∞, are as claimed.Proposition 3.5.
Letp∈(0,s0), s0>0. Then
(3.49)M(p±iR)~1cRγe∓iγπ/4
as R→∞.Proof.
Setμ:=p2+R2 and ν:=arctan(±R/p). Then μ~R and ν~±π/2, as R→∞. By (3.46), we have
(3.50)M(μeiν)~1cμγe∓iγν/4,μ⟶∞,
as claimed.In order to obtain the explicit form of solutionu to initial-boundary value problem (3.4), (3.5), and (3.6), it remains to calculate function P.Theorem 3.6.
The solutionu to initial-boundary value problem (3.4), (3.5), and (3.6) is given by (3.17), that is, u(x,t)=Υ(t)*P(x,t), where P takes the form
(3.51)P(x,t)=12πi∫0∞(sinh⁡(xqM(qe-iπ))sinh⁡(qM(qe-iπ))-sinh⁡(xqM(qeiπ))sinh⁡(qM(qeiπ)))e-qtdq+∑n=1∞[Res(P̃(x,s)est,Psn(+))+Res(P̃(x,s)est,Psn(-))],t>0.
The residues at simple poles Psn(±), n∈ℕ, are given by
(3.52)Res(P̃(x,s)est,Psn(±))=[sinh⁡(xsM(s))(d/ds)[sinh⁡(sM(s))]est]s=Psn(±).Proof.
FunctionP(x,t), x∈[0,1], t>0, will be calculated by integration over a suitable contour. The Cauchy residues theorem yields
(3.53)∮ΓP̃(x,s)estds=2πi∑n=1∞[Res(P̃(x,s)est,Psn(+))+Res(P̃(x,s)est,Psn(-))],
where Γ=Γ1∪Γ2∪Γ3∪Γɛ∪Γ4∪Γ5∪Γ6∪γ0 is such a contour that all poles lie inside the contour Γ (see Figure 1).    First we show that the series of residues in (3.51) is convergent. By Proposition 3.4, Psn(±) are simple poles of P̃, and therefore also simple poles of estP̃. The residues in (3.53) can be calculated as it is given in (3.52), so
(3.54)Res(P̃(x,s)est,Psn(±))=[1(d/ds)(sM(s))sinh⁡(xsM(s))cosh⁡(sM(s))est]s=Psn(±).
Using (3.27) we have that sinh(xsM(s))=±isin(xnπ), while cosh(sM(s))=(-1)n. Also, a calculation gives that (d/ds)(sM(s))=M(s)+M(s)·A(s), where
(3.55)A(s):=1+(a/b)(α-β)sα-β2(1+(a/b)sα-β)-aαsα+cγsγ+(ac/b)(α+γ-β)sα+γ-β2(asα+csγ+(ac/b)sα+γ-β).
Take now that Psn(±)=Re±iφ. Then
(3.56)Res(P̃(x,s)est,Psn(±))=(-1)nsin(nπx)nπReRtcosφe±i(φ+Rtsinφ)A(Re±iφ).
Let n→∞. Then also |Psn(±)|→∞, that is, R→∞, and
(3.57)|A(Re±iφ)|⟶1-γ2.
Further, it follows from Proposition 3.4 that
(3.58)Re(sPn(±))=Rcos⁡φ~c(nπ)22-γcos⁡(π2-γ)≤-Cn,forsomeC>0,
and by (3.48), R/n~cπ22-γ·nγ/(2-γ). Therefore, as n→∞,
(3.59)|Res(P̃(x,s)est,Psn(+))+Res(P̃(x,s)est,Psn(-))|≤|sin(nπx)nπReRtcos⁡φA(Re+iφ)|+|sin(nπx)nπReRtcos⁡φA(Re-iφ)|≤1πRne-Cnt(1|A(Re+iφ)|+1|A(Re-iφ)|)~4π(2-γ)cπ22-γ⋅nγ/(2-γ)e-Cnt,
which implies the convergence of the sum of residues in (3.53).
It remains to calculate the integral overΓ in (3.53). Consider the integral along contour Γ1:s=p+iR, s0>p>0. Then
(3.60)|∫Γ1P̃(x,s)estds|≤∫0s0|P̃(x,p+iR)||e(p+iR)t|dp.
Let R→∞. In order to estimate |P̃(x,p±iR)|, using Proposition 3.5, we write
(3.61)M(p±iR)~v±iw,v=1cRγcos⁡(γπ4),w=-1cRγsin(γπ4).
Then
(3.62)|P̃(x,p±iR)|~|sinh⁡[x(pv-Rw)±ix(pw+Rv)]sinh⁡[(pv-Rw)±i(pw+Rv)]|≤ex(pv-Rw)+e-x(pv-Rw)|epv-Rw-e-(pv-Rw)|=e-(1-x)(pv-Rw)1+e-2x(pv-Rw)|1-e-2(pv-Rw)|⟶0,asR⟶∞.
The above convergence is valid since
(3.63)pv-Rw=p1cRγcos⁡(γπ4)+R1cRγsin(γπ4)⟶∞,asR⟶∞.
Therefore, according to (3.62), we have
(3.64)lim⁡R→∞|∫Γ1P̃(x,s)estds|=0.
The similar argument is valid for the integral along Γ6, thus
(3.65)lim⁡R→∞|∫Γ6P̃(x,s)estds|=0.
Next, we consider the integral along contour Γ2:s=Reiφ, π/2<φ<π:
(3.66)|∫Γ2P̃(x,s)estds|≤∫π/2πR|e-R(1-x)eiφM(Reiφ)||1-e-2xReiφM(Reiφ)1-e-2ReiφM(Reiφ)|eRtcos⁡φdφ.
Since sM(s)→∞ as |s|→∞ (see Proposition 3.2 (ii)) and cosφ≤0 for φ∈[π/2,π], we have
(3.67)lim⁡R→∞|∫Γ2P̃(x,s)estds|≤lim⁡R→∞∫π/2πR|e-R(1-x)eiφM(Reiφ)|eRtcos⁡φdφ=0.
The similar argument is valid for the integral along Γ5, thus
(3.68)lim⁡R→∞|∫Γ5P̃(x,s)estds|=0.
Since sM(s)→0 as |s|→0 (see Proposition 3.2 (ii)), the integration along contour Γɛ:s=ɛeiφ, π>φ>-π, gives
(3.69)lim⁡ɛ→0|∫ΓɛP̃(x,s)estds|≤lim⁡ɛ→0∫π-πɛ|sinh⁡(xɛeiφM(ɛeiφ))sinh⁡(ɛeiφM(ɛeiφ))|eɛtcos⁡φdφ=0.
Integrals along parts of contour Γ3:s=qeiπ, R>q>ɛ, Γ4:s=qe-iπ, ɛ<q<R, and γ0:s=s0+ir, -R<r<R, give
(3.70)lim⁡R→∞ɛ→0∫Γ3P̃(x,s)estds=∫0∞sinh⁡(xqM(qeiπ))sinh⁡(qM(qeiπ))e-qtdq,lim⁡R→∞ɛ→0∫Γ4P̃(x,s)estds=-∫0∞sinh⁡(xqM(qe-iπ))sinh⁡(qM(qe-iπ))e-qtdq,lim⁡R→∞∫γ0P̃(x,s)estds=2πiP(x,t).
Equation (3.51) now follows from (3.53).Figure 1
Integration contourΓ.Corollary 3.7.
In the case of stress relaxation, that is, whenΥ(t)=Υ0H(t), Υ0>0, t∈ℝ, the solution takes the form
(3.71)uH(x,t)=Υ0H(t)*P(x,t),x∈[0,1],t∈R.
We will numerically examine it in the sequel.


#### 3.2.2. Determination of Stressσ in the Case of Prescribed Displacement Υ

In Section3.1 we determined stress σ (cf. (3.21)) that is a solution to (3.4), (3.5), and (3.6). In order to obtain an explicit form of σ we need to calculate function T. As in previous Section 3.2.1 it will be done by inversion of the Laplace transform of T̃.FunctionT̃, which is given by (3.20), is analytic on the complex plane except the branch cut (-∞,0], and has simple poles at the same points as P̃, that is, Psn(±), n∈ℕ.Using the similar arguments as in the proof of Theorem3.6, one can prove the following theorem.Theorem 3.8.
The solutionσ to initial-boundary value problem (3.4), (3.5), and (3.6) is given by (3.21), that is, σ(x,t)=d/dt(Υ(t)*T(x,t)), where T takes the form
(3.72)T(x,t)=12πi∫0∞(cosh⁡(xqM(qeiπ))M(qeiπ)sinh⁡(qM(qeiπ))-cosh⁡(xqM(qe-iπ))M(qe-iπ)sinh⁡(qM(qe-iπ)))e-qtdq+∑n=1∞[Res(T̃(x,s)est,Psn(+))+Res(T̃(x,s)est,Psn(-))],t>0.
The residues at simple poles Psn(±), n∈ℕ, are given by
(3.73)Res(T̃(x,s)est,Psn(±))=[cosh⁡(xsM(s))M(s)(d/ds)[sinh⁡(sM(s))]est]s=Psn(±).Corollary 3.9.
Similarly, in the case of stress relaxationΥ=Υ0H we obtain the solution
(3.74)σH(x,t)=Υ0T(x,t),x∈[0,1],t>0.In order to check our results for large times, we compare them with the quasistatic case. In the quasistatic case one uses only the constitutive equation (3.4)3, that is, the dynamics of the process is neglected. Taking the Laplace transform of the constitutive equation we obtain (3.9)2, and define the relaxation modulus G via its Laplace transform, as follows:
(3.75)G(t)∶=L-1[G̃(s)](t),t>0,G̃(s)∶=σ̃(QS)(s)ɛ̃(QS)(s)∶=asα+csγ+(ac/b)sα+γ-β1+(a/b)sα-β,s∈C∖(-∞,0].
Then the stress in the quasistatic case is
(3.76)σ(QS)=G*ɛ(QS).
Following the proof of Theorem 3.6, we obtain that (3.77)G(t)=12πi∫0∞[G̃(qe-iπ)-G̃(qeiπ)]e-qtdt.In the quasistatic case it holds that(3.78)u(x,t)=x⋅u(1,t),x∈[0,1]t>0,
and consequently by (3.4), (3.79)ɛ(x,t)=u(1,t)=:ɛ(QS)(t),x∈[0,1]t>0.
Since according to boundary condition (3.6), u(1,t)=Υ(t), it follows from (3.76) that (3.80)σ(QS)=Υ*G,
which, in the case of stress relaxation, that is, when Υ=Υ0H, Υ0>0, becomes(3.81)σH(QS)=Υ0H*G.


#### 3.2.3. Determination of Displacementu in the Case of Prescribed Stress Σ

In Section3.1 we determined displacement u (cf. (3.25)) that is a solution to (3.4), (3.5), and (3.7). As above, we now want to find u explicitly, by calculating the inverse Laplace transform of Q̃.FunctionQ̃ given by (3.24) is analytic on the complex plane except the branch cut (-∞,0] and has isolated singularities at solutions Qsn(±) of the equation
(3.82)cosh⁡(sM(s))=0thatis,sM(s)=±2n+12iπ,n∈N0.We state a proposition that is analogous to Propositions3.3 and 3.4. The proof is omitted since it follows the same lines as those of Propositions 3.3 and 3.4.Proposition 3.10.
(i) There are infinitely many complex conjugated solutionsQsn(±), n∈ℕ0, of (3.82), which all lie in the left complex half plane. Moreover, each Qsn(±), n∈ℕ0, is a simple pole.
(ii) SolutionsQsn(±), n∈ℕ0, of (3.82), are such that
(3.83)Re(sQn(±))=Rcos⁡φ~c(2n+12π)22-γcos⁡(π2-γ)<0,Im(sQn(±))=Rsinφ~±c(2n+12π)22-γsin(π2-γ),
as n→∞.In the following theorem we calculate explicitly displacementu.Theorem 3.11.
The solutionu to initial-boundary value problem (3.4), (3.5), and (3.7) is given by (3.25), that is, u(x,t)=Σ(t)*Q(x,t), where Q takes the form
(3.84)Q(x,t)=12πi∫0∞(M(qe-iπ)sinh⁡(xqM(qe-iπ))cosh⁡(qM(qe-iπ))-M(qeiπ)sinh⁡(xqM(qeiπ))cosh⁡(qM(qeiπ)))e-qtqdq+∑n=0∞[Res(Q̃(x,s)est,Qsn(+))+Res(Q̃(x,s)est,Qsn(-))],t>0.
The residues at simple poles Qsn(±), n∈ℕ0, are given by
(3.85)Res(Q̃(x,s)est,Qsn(±))=[1sM(s)sinh⁡(xsM(s))(d/ds)[cosh⁡(sM(s))]est]s=Qsn(±).Corollary 3.12.
The case of creep is described by the boundary conditionΣ(t)=Σ0H(t), Σ0>0, t∈ℝ, in which displacement u, given by (3.25), reads
(3.86)u(t)=Σ0H(t)*Q(x,t),x∈[0,1],t∈R.Similarly as in Section3.2.2, we examine the quasistatic case that corresponds to displacement u. Again, using only the constitutive equation (3.4)3 and its Laplace transform (3.9)2, we define the creep compliance J via its Laplace transform as
(3.87)J(t)∶=L-1[J̃(s)](t),t>0,J̃(s)∶=ɛ̃(QS)(s)σ̃(QS)(s)∶=1+(a/b)sα-βasα+csγ+(ac/b)sα+γ-β,s∈C∖(-∞,0].
Strain measure in the quasistatic case now equals(3.88)ɛ(QS)=J*σ(QS).
Following the proof of Theorem 3.6 one obtains (3.89)J(t)=12πi∫0∞[J̃(qe-iπ)-J̃(qeiπ)]e-qtdt.
In the quasistatic case it holds that (3.90)u(x,t)=x⋅u(1,t)=x⋅u(QS)(t),x∈[0,1],t>0,
and consequently by (3.4),
(3.91)ɛ(QS)=u(QS).
Also, σ(QS)(t):=σ(1,t)=Σ(t), t>0, which is the boundary condition (3.7), hence by (3.88) and (3.91), (3.92)u(QS)=Σ*J.
In the case of creep (cf. Corollary 3.12) we have(3.93)uH(QS)=Σ0H*J.


### 3.3. Numerical Examples

In this subsection we give several numerical examples of displacementuH and stress σH, given by (3.71) and (3.74), respectively, which correspond to the case of stress relaxation, and examine solutions (3.86), which correspond to displacement u in the case of creep. In addition, we investigate solutions uH, σH, and u for different orders of fractional derivatives.Figure2 presents displacements in a stress relaxation experiment, determined according to (3.71), for three different positions. Parameters in (3.71) are chosen as follows: Υ0=1, a=0.2, b=0.6, c=0.45, α=0.3, β=0.1, and γ=0.4. From Figure 2, one sees that the displacements in the case of stress relaxation show damped oscillatory character and that they tend to a constant value for large times, namely, limt→∞uH(x,t)=x, x∈[0,1]. Figure 3 presents the same displacements as Figure 2, but close to initial time instant. It is evident that there is a delay in displacement that increases as the point is further from the end where the prescribed displacement is applied. This is a consequence of the finite wave propagation speed.Figure 2
DisplacementsuH(x,t) in a stress relaxation experiment as functions of time t at x∈{0.25,0.5,0.75} for t∈(0,10).Figure 3
DisplacementsuH(x,t) in a stress relaxation experiment as functions of time t at x∈{0.25,0.5,0.75} for t∈(0,5).Figure4 presents the displacements uH of the points close to the end of the rod where the sudden but afterwards constant displacement is applied (i.e., u(1,t)=H(t), t>0). One sees that the amplitudes of these points deform in shape so that they do not exceed the prescribed value of the displacement of the rod's free end.Figure 4
DisplacementsuH(x,t) in a stress relaxation experiment as functions of time t at x∈{0.75,0.8,0.85,0.9,0.95} for t∈(0,5).In order to examine the influence of the orders of fractional derivatives in the constitutive equation (1.4) on the displacement uH and stress σH in a stress relaxation experiment, we plot the displacement uH and stress σH obtained by (3.71) and (3.74) for the following sets of parameters (α,β,γ)∈{(0.1,0.05,0.15),(0.3,0.1,0.4),(0.45,0.4,0.49)}, while we fix x=0.5 and leave other parameters as before. In the case of the first set, the constitutive equation (1.4) describes a body in which the elastic properties are dominant, since the orders of the fractional derivatives of stress and strain are close to zero, that is, the fractional derivatives of stress and strain almost coincide with the stress and strain. This is also evident from Figure 5, since the oscillations of the point x=0.5 for the first set of parameters vanish quite slowly comparing to the second set and in particular comparing with the third set of parameters. Note that the third set of parameters describes a body in which the fluid properties of the fractional type dominate, since in the constitutive equation (1.4), we have the low-order derivative of stress (almost the stress itself), while almost all the derivatives of strain are of order 0.5. Figure 5 also shows that the dissipative properties of a material grow as the orders of the fractional derivatives increase. Figure 6 shows that the delay in displacement depends on the order of the fractional derivative, so that a material which has dominant elastic properties (the first set) has the longest delay, compared to a material with the dominant fluid properties (the third set).Figure 5
DisplacementsuH(x,t) in a stress relaxation experiment as functions of time t at x=0.5 for t∈(0,20).Figure 6
DisplacementsuH(x,t) in a stress relaxation experiment as functions of time t at x=0.5 for t∈(0,8).Figure7 presents stresses σH in the case of stress relaxation, determined according to (3.74), for different points of the rod. Parameters are the same as in the previous case. Also, Figure 7 presents the quasistatic curve σH(QS) obtained by (3.81).Figure 7
StressesσH(x,t) and σH(QS)(t) in a stress relaxation experiment as functions of time t at x∈{0.1,0.2,0.3,0.4,0.5,0.6,0.7} for t∈(0,5).Stresses, as it can be seen from Figure7, show damped oscillatory character and for large times, in each point x∈[0,1], tend to the quasistatic curve, that is, to the same value. Eventually, the stresses in all points of the rod tend to zero, namely, limt→∞σH(x,t)=0, x∈[0,1]. From Figure 7 it is evident that as further the point is from the rod-free end the greater is the delay. This is again the consequence of the finite wave speed.Figures8 and 9 present the stresses of the points close to the free end. One notices from Figures 7 and 8 that as the point is closer to the end whose displacement is prescribed, the peak of stress is higher and the peak's width is smaller. Figure 9 presents the compressive phase in the stress relaxation process and the quasistatic curve as well.Figure 8
StressesσH(x,t) in a stress relaxation experiment as functions of time t at x∈{0.7,0.8,0.9} for t∈(0,0.5).Figure 9
StressesσH(x,t) and σH(QS)(t) in a stress relaxation experiment as functions of time t at x∈{0.7,0.8,0.9} for t∈(0,6).Again, we fix the midpoint of the rod and investigate the influence of the change of orders of the fractional derivatives (for the same three sets as before) on the stressσH in a stress relaxation experiment. From Figure 10 one notices that there is a stress relaxation in a material regardless of the order of fractional derivatives. However, the relaxed stress depends on the order of the derivatives, as expected from the analysis of the stress σH(QS) in the quasistatic case, given by (3.81). The material with the dominant elastic properties (the first set) relaxes to the highest stress, and as the fluid properties of the material become more and more dominant, the relaxed stress decreases. Again, the oscillations of the value of the stress for the first set (material with dominant elastic properties) are the least damped compared to the second and third set. Figure 11 shows that the conclusion about the dependence of delay on the orders of fractional derivatives drown earlier holds.Figure 10
StressesσH(x,t) in a stress relaxation experiment as functions of time t at x=0.5 for t∈(0,10).Figure 11
StressesσH(x,t) in a stress relaxation experiment as functions of time t at x=0.5 for t∈(0,4.5).Figure12 presents displacements u in the creep experiment, determined according to (3.86), for four different points. Parameters are the same as in the previous cases (except that in this case instead of Υ0=1 we have Σ0=1). For large times, as it can be seen from Figure 12, the displacement curves are monotonically increasing. This indicates that we deal with the viscoelastic fluid. Figure 12 also shows good agreement between the displacements obtained by the dynamic model (displacement is given by (3.86)) and the quasistatic model (displacement is given by (3.93)). Figure 13 presents the same displacements as Figure 12, but close to initial time instant and it is evident that, again, there is a delay, due to the finite speed of wave propagation.Figure 12
Displacementsu(x,t) and uH(QS)(t) in a creep experiment as functions of time t at x∈{0.25,0.5,0.75,1} for t∈(0,40).Figure 13
Displacementsu(x,t) in creep experiment as functions of time t at x∈{0.25,0.5,0.75,1} for t∈(0,10).In order to examine the dependence of the displacementu in the case of creep experiment, given by (3.86), on the orders of fractional derivatives in the constitutive equation (1.4), we again fix the point x=0.5 of the rod and plot the displacement u for the same set of values of α, β, and γ as before. Figure 14 clearly shows that all of the materials exhibits creep, but the displacements do not tend to a constant value. However, the material which has the dominant elastic properties (the first set) creeps slower comparing to the material with the dominant fluid properties (the third set). Figure 15 again supports the conclusion that the more elastic the material is, the greater the delay is.Figure 14
Displacementsu(x,t) in creep experiment as functions of time t at x=0.5 for t∈(0,40).Figure 15
Displacementsu(x,t) in creep experiment as functions of time t at x=0.5 for t∈(0,10).


## 3.1. Convolution Form of Solutions

Consider system (1.4), (1.5), and (1.6) supplied with initial conditions
(3.1)u(x,0)=0,∂∂tu(x,0)=0,σ(x,0)=0,ɛ(x,0)=0,x∈[0,L],
as well as with two types of boundary conditions:(3.2)u(0,t)=0,u(L,t)=Υ(t),t∈R,u(0,t)=0,σ(L,t)=Σ(t),t∈R.
Functions Υ and Σ are locally integrable functions supported in [0,∞). Note that if Υ=Υ0H we have the case of stress relaxation, while if Σ=Σ0H we have the case of creep, where H is the Heaviside function.Introducing dimensionless quantities(3.3)x¯=xL,t¯=tLρ/E,u¯=uL,σ¯=σE,Υ¯=ΥL,Σ¯=ΣE,a¯=a(Lρ/E)α,b¯=b(Lρ/E)β,c¯=c(Lρ/E)γ,
and using the fact that fractional derivatives are transformed as 0Dt¯ηu(t¯)=(Lρ/E)η0Dtηu(t) (cf. [17, 23]), after omitting the bar over dimensionless quantities, we obtain the following system: x∈[0,1], t>0,
(3.4)∂∂xσ(x,t)=∂2∂t2u(x,t),ɛ(x,t)=∂∂xu(x,t),(1+abDtα-β0)σ(x,t)=E(aDtα0+cDtγ0+acbDtα+γ-β0)ɛ(x,t).
System (3.4) is subject to initial(3.5)u(x,0)=0,∂∂tu(x,0)=0,σ(x,0)=0,ɛ(x,0)=0,x∈[0,1],
and two types of boundary conditions
(3.6)u(0,t)=0,u(1,t)=Υ(t),t∈R,(3.7)u(0,t)=0,σ(1,t)=Σ(t),t∈R
as described above.Solutions of the above system will be determined by the Laplace transform method. Recall, the Laplace transform off∈Lloc1(ℝ), f≡0 in (-∞,0), and |f(t)|≤Meat, for some M>0, a∈ℝ and t≥0, is defined by (3.8)f̃(s)=L[f(t)](s)∶=∫0∞e-stf(t)dt,Res>a.
Thus, applying the Laplace transform to (3.4) and (3.5), one obtains  (3.9)∂∂xσ̃(x,s)=s2ũ(x,s),ɛ̃(x,s)=∂∂xũ(x,s),(1+absα-β)σ̃(x,s)=(asα+csγ+acbsα+γ-β)ɛ̃(x,s).
System (3.9) reduces to (3.10)∂2∂x2ũ(x,s)-(sM(s))2ũ(x,s)=0,x∈[0,1],s∈D,
where(3.11)M(s)∶=1+(a/b)sα-βasα+csγ+(ac/b)sα+γ-β=1asα1+(a/b)sα-β1+(c/a)sγ-α+(c/b)sγ-β.
It has a solution(3.12)ũ(x,s)=C1(s)exsM(s)+C2(s)e-xsM(s),x∈[0,1],s∈D,
where C1 and C2 are functions of s which will be determined from the boundary conditions. Since the power function sη is analytic on the complex plane except the branch cut along the negative axis (including the origin) we take D:=ℂ∖(-∞,0] to be the domain for variable s in (3.9). Applying either (3.6)1 or (3.7)1, we obtain C1=-C2=:C, and thus(3.13)ũ(x,s)=C(s)(exsM(s)-e-xsM(s)),x∈[0,1],s∈C∖(-∞,0].
From (3.9) and (3.11) it follows that(3.14)σ̃(x,s)=1M2(s)∂∂xũ(x,s),x∈[0,1],s∈C∖(-∞,0].As announced in Section1 we will separately seek solutions in different cases: displacement u and stress σ in the case of prescribed displacement (such as e.g., stress relaxation), and displacement u in the case of prescribed stress (e.g., creep). For the former we supply to the system boundary conditions (3.6), while for the latter we assume (3.7). In the sequel we derive convolution forms of solutions in all these cases.In the case of prescribed displacementΥ, substituting (3.13) into (3.12) and using (3.6), one obtains(3.15)ũ(x,s)=Υ̃(s)P̃(x,s),x∈[0,1],s∈C∖(-∞,0],
where(3.16)P̃(x,s)∶=sinh⁡(xsM(s))sinh⁡(sM(s)),x∈[0,1],s∈C∖(-∞,0].
Clearly, P(1,t)=δ(t), t∈ℝ.SinceΥ and P are supported in [0,∞), displacement u is given by  (3.17)u(x,t)=Υ(t)*P(x,t),x∈[0,1],t∈R,(3.18)u(x,t)=0,x∈[0,1],t<0,
where * denotes the convolution with respect to t. Recall, if f,g∈Lloc1(ℝ), suppf,g⊂[0,∞), then (f*g)(t):=∫0tf(τ)g(t-τ)dτ, t∈ℝ. Explicit calculation of (3.17) will be done by the use of the Laplace inversion formula applied to (3.15).Further, from (3.14), (3.15), and (3.16) it follows that(3.19)σ̃(x,s)=sΥ̃(s)T̃(x,s),x∈[0,1],s∈C∖(-∞,0],
where(3.20)T̃(x,s)=cosh⁡(xsM(s))M(s)sinh⁡(sM(s)),x∈[0,1],s∈C∖(-∞,0].
Applying the Laplace inversion formula to (3.19) we obtain(3.21)σ(x,t)=ddt(Υ(t)*T(x,t)),x∈[0,1],t∈R,
where the derivative is understood in the sense of distributions. Again, σ(x,t)=0 for x∈[0,1], t<0.In the case of prescribed stressΣ, using (3.14) at x=1 and (3.7), we obtain(3.22)∂∂xũ(1,s)=Σ̃(s)M2(s),s∈C∖(-∞,0].
This combined with (3.13) gives(3.23)ũ(x,s)=Σ̃(s)Q̃(x,s),x∈[0,1],s∈C∖(-∞,0],
where(3.24)Q̃(x,s)=1sM(s)sinh⁡(xsM(s))cosh⁡(sM(s)),x∈[0,1],s∈C∖(-∞,0].
Again, applying the Laplace inversion formula to (3.23), the displacement reads  (3.25)u(x,t)=Σ(t)*Q(x,t),x∈[0,1],t>0,(3.26)u(x,t)=0,x∈[0,1],t<0.Remark 3.1.
As in [15] one can show that P, T, and Q are real-valued, locally integrable functions on ℝ, supported in [0,∞) and smooth for t>0.


## 3.2. Explicit Forms of Solutions

This subsection is devoted to the calculation of inverse Laplace transforms of certain distributions and functions introduced above. We begin with examining some basic properties ofM given in (3.11), which will be important for further investigations. In the sequel we shall write A(x)~B(x) if limx→∞(A(x)/B(x))=1.Proposition 3.2.
LetM be the function defined by (3.11). Then  (i)
M is an analytic function on ℂ∖(-∞,0] if  0<β<α<γ<1.(ii)
Fors∈ℂ∖(-∞,0], lim|s|→0M(s)=∞, lim|s|→0sM(s)=0, lim|s|→∞M(s)=0, and lim|s|→∞sM(s)=∞.Proof.
(i) Since1+(a/b)sα-β≠0 and 1+(c/a)sγ-α+(c/b)sγ-β≠0 if args∈(-π,π) and 0<β<α<γ<1, it follows that the function M given in (3.11) is analytic on the complex plane except the branch cut along the negative axis, that is, on ℂ∖(-∞,0].
(ii) Limits in (ii) can easily be calculated.


### 3.2.1. Determination of Displacementu in the Case of Prescribed Displacement Υ

To begin with, we examine properties ofP̃ given by (3.16). P̃ has isolated singularities at Psn(±), n∈ℕ, where Psn(±) denotes solutions of the equation
(3.27)sinh⁡(sM(s))=0,that is,sM(s)=±niπ.Let us examine their position and multiplicity.Proposition 3.3.
There are infinitely many complex conjugated solutionsPsn(±), n∈ℕ, of (3.27), which all lie in the left complex half plane. Moreover, each Psn(±), n∈ℕ, is a simple pole.Proof.
Let us square (3.27) and define
(3.28)Φ(s,n)∶=(sM(s))2+(nπ)2=s21+(a/b)sα-βasα+csγ+(ac/b)sα+γ-β+(nπ)2.
Writing s=Reiφ it follows that
(3.29)Φ(Reiφ,n)=R2e2iφ1+(a/b)Rα-βei(α-β)φaRαeiαφ+cRγeiγφ+(ac/b)Rα+γ-βei(α+γ-β)φ+(nπ)2=R2(cos⁡(2φ)+isin(2φ))A+iBC+iD+(nπ)2,
where
(3.30)A∶=1+abRα-βcos⁡((α-β)φ),B∶=abRα-βsin((α-β)φ),C∶=aRαcos⁡(αφ)+cRγcos⁡(γφ)+acbRα+γ-βcos⁡((α+γ-β)φ),D∶=aRαsin(αφ)+cRγsin(γφ)+acbRα+γ-βsin((α+γ-β)φ).
Real and imaginary parts of Φ are then given by
(3.31)ReΦ(Reiφ,n)=R2C2+D2[(AC+BD)cos⁡(2φ)+(AD-BC)sin(2φ)]+(nπ)2,ImΦ(Reiφ,n)=R2C2+D2[(AC+BD)sin(2φ)-(AD-BC)cos(2φ)],(3.32)AC+BD=aRαcos⁡(αφ)+cRγcos⁡(γφ)+a2bR2α-βcos⁡(βφ)+a2cb2R2(α-β)+γcos⁡(γφ)+2acbRα+γ-βcos⁡((α-β)φ)cos⁡(γφ),AD-BC=aRαsin(αφ)+cRγsin(γφ)+a2bR2α-βsin(βφ)+a2cb2R2(α-β)+γsin(γφ)+2acbRα+γ-βcos⁡((α-β)φ)sin(γφ).
Let(R,φ) be a solution to Φ(s,n)=0 (or equivalently, ReΦ=ImΦ=0). Then changing φ→-φ, we again obtain that ReΦ=ImΦ=0, which implies that solutions of (3.27) are complex conjugated.
Further, we will prove that the functionΦ has no zeros in the half plane args∈[0,π/2]. For that purpose we will use the argument principle. Recall, if Φ is an analytic function inside and on a regular closed curve c and nonzero on c, then the number of zeros of Φ is given by NZ=(1/2π)ΔargΦ(s). Let γ=γa∪γb∪γc be parametrized as
(3.33)γa:s=x,x∈[0,R],γb:s=Reiφ,φ∈[0,π2],γc:s=xeiπ/2,x∈[0,R],
and let R→∞. Along γa function Φ becomes a real-valued function, hence ΔargΦ(s,n)=0. Along γb we have that AC+BD,AD-BC≥0 for φ∈[0,π], since sin(ηφ),cos(ηφ)>0, for η∈{α,β,γ,α-β}⊂(0,1/2), and φ∈[0,π]. Therefore, (3.31) imply that
(3.34)ReΦ(Reiφ,n)>0,φ∈[0,π4],ImΦ(Reiφ,n)>0,φ∈[π4,π2]
(inequalities at boundary points are easily checked by inserting them into (3.31)). Along γc, we obtain
(3.35)ReΦ(xeiπ/2,n)=-R2C2+D2|R=x,φ=π/2(AC+BD)|R=x,φ=π/2+(nπ)2,(3.36)ImΦ(xeiπ/2,n)=R2C2+D2|R=x,φ=π/2(AD-BC)|R=x,φ=π/2>0.
From (3.34) and (3.36), we may now conclude that along γb and γc,
(3.37)ΔargΦ(s,n)=0.
Indeed, this follows from the following. Along γb, in the case φ∈[0,π/4], ImΦ(Reiφ,n) can change its sign, but ReΦ(Reiφ,n)>0, while for φ∈[π/4,π/2], ImΦ(Reiφ,n)>0. Along γc, ImΦ(xeiπ/2,n)>0. Therefore, there is no change in the argument of Φ along the whole γ, which implies that Φ has no zeros for φ∈[0,π/2]. This further implies that (3.27) has no solutions in the right complex half plane, since its solutions are complex conjugated.
In order to prove that for fixedn∈ℕ (3.27) has one solution (and its complex conjugate), we again use function Φ and the argument principle. Consider now the contour Γ=ΓA∪ΓB∪ΓC parametrized by
(3.38)ΓA:s=xeiπ/2,x∈[0,R],ΓB:s=Reiφ,φ∈[π2,π],ΓC:s=xeiπ,x∈[0,R],
and let R→∞. Along ΓA real and imaginary parts of Φ are given by (3.35) and (3.36). Along ΓB, using (3.31), we conclude that
(3.39)ReΦ(Reiφ,n)<0,φ∈[π2,3π4],forR⟶∞,andnfixed,ImΦ(Reiφ,n)<0,φ∈[3π4,π].
Along ΓC we have
(3.40)ReΦ(Reiφ,n)=R2C2+D2|R=x,φ=π(AC+BD)|R=x,φ=π+(nπ)2>0,ImΦ(Reiφ,n)=-R2C2+D2|R=x,φ=π(AD-BC)|R=x,φ=π<0.
Along ΓA we have that ImΦ>0, while ReΦ changes its sign (since ReΦ(0,n)=(nπ)2 and limx→∞ReΦ(xeiπ/2,n)=-∞, for fixed n∈ℕ). Along the part of ΓB where φ∈[3π/4,π], and along ΓC, ImΦ<0. Also, limR→∞ReΦ(Rei3π/4,n)=-∞ and ImΦ(Rei3π/4,n)<0. This implies that the argument of Φ changes from 0 to 2π.
As a conclusion one obtains that alongΓ(3.41)ΔargΦ(s,n)=2π,
which implies, by the argument principle, that function Φ has exactly one zero in the upper left complex plane, for each fixed n∈ℕ. Since the zeros of Φ are complex conjugated, it follows that Φ also has one zero in the lower left complex plane, for each fixed n∈ℕ.In the next proposition we examine behavior of simple polesPsn(±), n∈ℕ.Proposition 3.4.
SolutionsPsn(±), n∈ℕ, of (3.27), are such that
(3.42)Re(Psn(±))=Rcos⁡φ~c(nπ)22-γcos⁡(π2-γ)<0,Im(Psn(±))=Rsinφ~±c(nπ)22-γsin(π2-γ),
as n→∞.Proof.
Let us square (3.27) and insert Psn(±)=Reiφ, φ∈(-π,π). Then, after separation of real and imaginary parts, we obtain
(3.43)R2cos⁡(2φ)Re(M2(Reiφ))-R2sin(2φ)Im⁡(M2(Reiφ))=-(nπ)2,(3.44)R2sin(2φ)Re(M2(Reiφ))+R2cos⁡(2φ)Im⁡(M2(Reiφ))=0.
Using notation from the proof of Proposition 3.3 we can write
(3.45)Re(M2(Reiφ))=AC+BDC2+D2,Im⁡(M2(Reiφ))=BC-ADC2+D2.
Letting R→∞, one has
(3.46)Re(M2(Reiφ))~(a2c/b2)R2(α-β)+γcos⁡(γφ)(a2c2/b2)R2(α-β)+2γ=1cRγcos⁡(γφ),Im⁡(M2(Reiφ))~-(a2c/b2)R2(α-β)+γsin(γφ)(a2c2/b2)R2(α-β)+2γ=-1cRγsin(γφ).
It now follows from (3.44) and (3.46) that
(3.47)tg(2φ)=-Im(M2(Reiφ))Re(M2(Reiφ))~tg(γφ)⟹sin((2-γ)φ)cos⁡(2φ)cos⁡(γφ)~0⟹φ~±π2-γ.
Inserting (3.47) into (3.46), and subsequently into (3.43), we obtain
(3.48)R2-γccos⁡(2π2-γ)cos⁡(γπ2-γ)+R2-γcsin(2π2-γ)sin(γπ2-γ)~-(nπ)2,R~c(nπ)22-γ.
Thus, real and imaginary parts of Psn(±), as R→∞, are as claimed.Proposition 3.5.
Letp∈(0,s0), s0>0. Then
(3.49)M(p±iR)~1cRγe∓iγπ/4
as R→∞.Proof.
Setμ:=p2+R2 and ν:=arctan(±R/p). Then μ~R and ν~±π/2, as R→∞. By (3.46), we have
(3.50)M(μeiν)~1cμγe∓iγν/4,μ⟶∞,
as claimed.In order to obtain the explicit form of solutionu to initial-boundary value problem (3.4), (3.5), and (3.6), it remains to calculate function P.Theorem 3.6.
The solutionu to initial-boundary value problem (3.4), (3.5), and (3.6) is given by (3.17), that is, u(x,t)=Υ(t)*P(x,t), where P takes the form
(3.51)P(x,t)=12πi∫0∞(sinh⁡(xqM(qe-iπ))sinh⁡(qM(qe-iπ))-sinh⁡(xqM(qeiπ))sinh⁡(qM(qeiπ)))e-qtdq+∑n=1∞[Res(P̃(x,s)est,Psn(+))+Res(P̃(x,s)est,Psn(-))],t>0.
The residues at simple poles Psn(±), n∈ℕ, are given by
(3.52)Res(P̃(x,s)est,Psn(±))=[sinh⁡(xsM(s))(d/ds)[sinh⁡(sM(s))]est]s=Psn(±).Proof.
FunctionP(x,t), x∈[0,1], t>0, will be calculated by integration over a suitable contour. The Cauchy residues theorem yields
(3.53)∮ΓP̃(x,s)estds=2πi∑n=1∞[Res(P̃(x,s)est,Psn(+))+Res(P̃(x,s)est,Psn(-))],
where Γ=Γ1∪Γ2∪Γ3∪Γɛ∪Γ4∪Γ5∪Γ6∪γ0 is such a contour that all poles lie inside the contour Γ (see Figure 1).    First we show that the series of residues in (3.51) is convergent. By Proposition 3.4, Psn(±) are simple poles of P̃, and therefore also simple poles of estP̃. The residues in (3.53) can be calculated as it is given in (3.52), so
(3.54)Res(P̃(x,s)est,Psn(±))=[1(d/ds)(sM(s))sinh⁡(xsM(s))cosh⁡(sM(s))est]s=Psn(±).
Using (3.27) we have that sinh(xsM(s))=±isin(xnπ), while cosh(sM(s))=(-1)n. Also, a calculation gives that (d/ds)(sM(s))=M(s)+M(s)·A(s), where
(3.55)A(s):=1+(a/b)(α-β)sα-β2(1+(a/b)sα-β)-aαsα+cγsγ+(ac/b)(α+γ-β)sα+γ-β2(asα+csγ+(ac/b)sα+γ-β).
Take now that Psn(±)=Re±iφ. Then
(3.56)Res(P̃(x,s)est,Psn(±))=(-1)nsin(nπx)nπReRtcosφe±i(φ+Rtsinφ)A(Re±iφ).
Let n→∞. Then also |Psn(±)|→∞, that is, R→∞, and
(3.57)|A(Re±iφ)|⟶1-γ2.
Further, it follows from Proposition 3.4 that
(3.58)Re(sPn(±))=Rcos⁡φ~c(nπ)22-γcos⁡(π2-γ)≤-Cn,forsomeC>0,
and by (3.48), R/n~cπ22-γ·nγ/(2-γ). Therefore, as n→∞,
(3.59)|Res(P̃(x,s)est,Psn(+))+Res(P̃(x,s)est,Psn(-))|≤|sin(nπx)nπReRtcos⁡φA(Re+iφ)|+|sin(nπx)nπReRtcos⁡φA(Re-iφ)|≤1πRne-Cnt(1|A(Re+iφ)|+1|A(Re-iφ)|)~4π(2-γ)cπ22-γ⋅nγ/(2-γ)e-Cnt,
which implies the convergence of the sum of residues in (3.53).
It remains to calculate the integral overΓ in (3.53). Consider the integral along contour Γ1:s=p+iR, s0>p>0. Then
(3.60)|∫Γ1P̃(x,s)estds|≤∫0s0|P̃(x,p+iR)||e(p+iR)t|dp.
Let R→∞. In order to estimate |P̃(x,p±iR)|, using Proposition 3.5, we write
(3.61)M(p±iR)~v±iw,v=1cRγcos⁡(γπ4),w=-1cRγsin(γπ4).
Then
(3.62)|P̃(x,p±iR)|~|sinh⁡[x(pv-Rw)±ix(pw+Rv)]sinh⁡[(pv-Rw)±i(pw+Rv)]|≤ex(pv-Rw)+e-x(pv-Rw)|epv-Rw-e-(pv-Rw)|=e-(1-x)(pv-Rw)1+e-2x(pv-Rw)|1-e-2(pv-Rw)|⟶0,asR⟶∞.
The above convergence is valid since
(3.63)pv-Rw=p1cRγcos⁡(γπ4)+R1cRγsin(γπ4)⟶∞,asR⟶∞.
Therefore, according to (3.62), we have
(3.64)lim⁡R→∞|∫Γ1P̃(x,s)estds|=0.
The similar argument is valid for the integral along Γ6, thus
(3.65)lim⁡R→∞|∫Γ6P̃(x,s)estds|=0.
Next, we consider the integral along contour Γ2:s=Reiφ, π/2<φ<π:
(3.66)|∫Γ2P̃(x,s)estds|≤∫π/2πR|e-R(1-x)eiφM(Reiφ)||1-e-2xReiφM(Reiφ)1-e-2ReiφM(Reiφ)|eRtcos⁡φdφ.
Since sM(s)→∞ as |s|→∞ (see Proposition 3.2 (ii)) and cosφ≤0 for φ∈[π/2,π], we have
(3.67)lim⁡R→∞|∫Γ2P̃(x,s)estds|≤lim⁡R→∞∫π/2πR|e-R(1-x)eiφM(Reiφ)|eRtcos⁡φdφ=0.
The similar argument is valid for the integral along Γ5, thus
(3.68)lim⁡R→∞|∫Γ5P̃(x,s)estds|=0.
Since sM(s)→0 as |s|→0 (see Proposition 3.2 (ii)), the integration along contour Γɛ:s=ɛeiφ, π>φ>-π, gives
(3.69)lim⁡ɛ→0|∫ΓɛP̃(x,s)estds|≤lim⁡ɛ→0∫π-πɛ|sinh⁡(xɛeiφM(ɛeiφ))sinh⁡(ɛeiφM(ɛeiφ))|eɛtcos⁡φdφ=0.
Integrals along parts of contour Γ3:s=qeiπ, R>q>ɛ, Γ4:s=qe-iπ, ɛ<q<R, and γ0:s=s0+ir, -R<r<R, give
(3.70)lim⁡R→∞ɛ→0∫Γ3P̃(x,s)estds=∫0∞sinh⁡(xqM(qeiπ))sinh⁡(qM(qeiπ))e-qtdq,lim⁡R→∞ɛ→0∫Γ4P̃(x,s)estds=-∫0∞sinh⁡(xqM(qe-iπ))sinh⁡(qM(qe-iπ))e-qtdq,lim⁡R→∞∫γ0P̃(x,s)estds=2πiP(x,t).
Equation (3.51) now follows from (3.53).Figure 1
Integration contourΓ.Corollary 3.7.
In the case of stress relaxation, that is, whenΥ(t)=Υ0H(t), Υ0>0, t∈ℝ, the solution takes the form
(3.71)uH(x,t)=Υ0H(t)*P(x,t),x∈[0,1],t∈R.
We will numerically examine it in the sequel.


### 3.2.2. Determination of Stressσ in the Case of Prescribed Displacement Υ

In Section3.1 we determined stress σ (cf. (3.21)) that is a solution to (3.4), (3.5), and (3.6). In order to obtain an explicit form of σ we need to calculate function T. As in previous Section 3.2.1 it will be done by inversion of the Laplace transform of T̃.FunctionT̃, which is given by (3.20), is analytic on the complex plane except the branch cut (-∞,0], and has simple poles at the same points as P̃, that is, Psn(±), n∈ℕ.Using the similar arguments as in the proof of Theorem3.6, one can prove the following theorem.Theorem 3.8.
The solutionσ to initial-boundary value problem (3.4), (3.5), and (3.6) is given by (3.21), that is, σ(x,t)=d/dt(Υ(t)*T(x,t)), where T takes the form
(3.72)T(x,t)=12πi∫0∞(cosh⁡(xqM(qeiπ))M(qeiπ)sinh⁡(qM(qeiπ))-cosh⁡(xqM(qe-iπ))M(qe-iπ)sinh⁡(qM(qe-iπ)))e-qtdq+∑n=1∞[Res(T̃(x,s)est,Psn(+))+Res(T̃(x,s)est,Psn(-))],t>0.
The residues at simple poles Psn(±), n∈ℕ, are given by
(3.73)Res(T̃(x,s)est,Psn(±))=[cosh⁡(xsM(s))M(s)(d/ds)[sinh⁡(sM(s))]est]s=Psn(±).Corollary 3.9.
Similarly, in the case of stress relaxationΥ=Υ0H we obtain the solution
(3.74)σH(x,t)=Υ0T(x,t),x∈[0,1],t>0.In order to check our results for large times, we compare them with the quasistatic case. In the quasistatic case one uses only the constitutive equation (3.4)3, that is, the dynamics of the process is neglected. Taking the Laplace transform of the constitutive equation we obtain (3.9)2, and define the relaxation modulus G via its Laplace transform, as follows:
(3.75)G(t)∶=L-1[G̃(s)](t),t>0,G̃(s)∶=σ̃(QS)(s)ɛ̃(QS)(s)∶=asα+csγ+(ac/b)sα+γ-β1+(a/b)sα-β,s∈C∖(-∞,0].
Then the stress in the quasistatic case is
(3.76)σ(QS)=G*ɛ(QS).
Following the proof of Theorem 3.6, we obtain that (3.77)G(t)=12πi∫0∞[G̃(qe-iπ)-G̃(qeiπ)]e-qtdt.In the quasistatic case it holds that(3.78)u(x,t)=x⋅u(1,t),x∈[0,1]t>0,
and consequently by (3.4), (3.79)ɛ(x,t)=u(1,t)=:ɛ(QS)(t),x∈[0,1]t>0.
Since according to boundary condition (3.6), u(1,t)=Υ(t), it follows from (3.76) that (3.80)σ(QS)=Υ*G,
which, in the case of stress relaxation, that is, when Υ=Υ0H, Υ0>0, becomes(3.81)σH(QS)=Υ0H*G.


### 3.2.3. Determination of Displacementu in the Case of Prescribed Stress Σ

In Section3.1 we determined displacement u (cf. (3.25)) that is a solution to (3.4), (3.5), and (3.7). As above, we now want to find u explicitly, by calculating the inverse Laplace transform of Q̃.FunctionQ̃ given by (3.24) is analytic on the complex plane except the branch cut (-∞,0] and has isolated singularities at solutions Qsn(±) of the equation
(3.82)cosh⁡(sM(s))=0thatis,sM(s)=±2n+12iπ,n∈N0.We state a proposition that is analogous to Propositions3.3 and 3.4. The proof is omitted since it follows the same lines as those of Propositions 3.3 and 3.4.Proposition 3.10.
(i) There are infinitely many complex conjugated solutionsQsn(±), n∈ℕ0, of (3.82), which all lie in the left complex half plane. Moreover, each Qsn(±), n∈ℕ0, is a simple pole.
(ii) SolutionsQsn(±), n∈ℕ0, of (3.82), are such that
(3.83)Re(sQn(±))=Rcos⁡φ~c(2n+12π)22-γcos⁡(π2-γ)<0,Im(sQn(±))=Rsinφ~±c(2n+12π)22-γsin(π2-γ),
as n→∞.In the following theorem we calculate explicitly displacementu.Theorem 3.11.
The solutionu to initial-boundary value problem (3.4), (3.5), and (3.7) is given by (3.25), that is, u(x,t)=Σ(t)*Q(x,t), where Q takes the form
(3.84)Q(x,t)=12πi∫0∞(M(qe-iπ)sinh⁡(xqM(qe-iπ))cosh⁡(qM(qe-iπ))-M(qeiπ)sinh⁡(xqM(qeiπ))cosh⁡(qM(qeiπ)))e-qtqdq+∑n=0∞[Res(Q̃(x,s)est,Qsn(+))+Res(Q̃(x,s)est,Qsn(-))],t>0.
The residues at simple poles Qsn(±), n∈ℕ0, are given by
(3.85)Res(Q̃(x,s)est,Qsn(±))=[1sM(s)sinh⁡(xsM(s))(d/ds)[cosh⁡(sM(s))]est]s=Qsn(±).Corollary 3.12.
The case of creep is described by the boundary conditionΣ(t)=Σ0H(t), Σ0>0, t∈ℝ, in which displacement u, given by (3.25), reads
(3.86)u(t)=Σ0H(t)*Q(x,t),x∈[0,1],t∈R.Similarly as in Section3.2.2, we examine the quasistatic case that corresponds to displacement u. Again, using only the constitutive equation (3.4)3 and its Laplace transform (3.9)2, we define the creep compliance J via its Laplace transform as
(3.87)J(t)∶=L-1[J̃(s)](t),t>0,J̃(s)∶=ɛ̃(QS)(s)σ̃(QS)(s)∶=1+(a/b)sα-βasα+csγ+(ac/b)sα+γ-β,s∈C∖(-∞,0].
Strain measure in the quasistatic case now equals(3.88)ɛ(QS)=J*σ(QS).
Following the proof of Theorem 3.6 one obtains (3.89)J(t)=12πi∫0∞[J̃(qe-iπ)-J̃(qeiπ)]e-qtdt.
In the quasistatic case it holds that (3.90)u(x,t)=x⋅u(1,t)=x⋅u(QS)(t),x∈[0,1],t>0,
and consequently by (3.4),
(3.91)ɛ(QS)=u(QS).
Also, σ(QS)(t):=σ(1,t)=Σ(t), t>0, which is the boundary condition (3.7), hence by (3.88) and (3.91), (3.92)u(QS)=Σ*J.
In the case of creep (cf. Corollary 3.12) we have(3.93)uH(QS)=Σ0H*J.


## 3.2.1. Determination of Displacementu in the Case of Prescribed Displacement Υ

To begin with, we examine properties ofP̃ given by (3.16). P̃ has isolated singularities at Psn(±), n∈ℕ, where Psn(±) denotes solutions of the equation
(3.27)sinh⁡(sM(s))=0,that is,sM(s)=±niπ.Let us examine their position and multiplicity.Proposition 3.3.
There are infinitely many complex conjugated solutionsPsn(±), n∈ℕ, of (3.27), which all lie in the left complex half plane. Moreover, each Psn(±), n∈ℕ, is a simple pole.Proof.
Let us square (3.27) and define
(3.28)Φ(s,n)∶=(sM(s))2+(nπ)2=s21+(a/b)sα-βasα+csγ+(ac/b)sα+γ-β+(nπ)2.
Writing s=Reiφ it follows that
(3.29)Φ(Reiφ,n)=R2e2iφ1+(a/b)Rα-βei(α-β)φaRαeiαφ+cRγeiγφ+(ac/b)Rα+γ-βei(α+γ-β)φ+(nπ)2=R2(cos⁡(2φ)+isin(2φ))A+iBC+iD+(nπ)2,
where
(3.30)A∶=1+abRα-βcos⁡((α-β)φ),B∶=abRα-βsin((α-β)φ),C∶=aRαcos⁡(αφ)+cRγcos⁡(γφ)+acbRα+γ-βcos⁡((α+γ-β)φ),D∶=aRαsin(αφ)+cRγsin(γφ)+acbRα+γ-βsin((α+γ-β)φ).
Real and imaginary parts of Φ are then given by
(3.31)ReΦ(Reiφ,n)=R2C2+D2[(AC+BD)cos⁡(2φ)+(AD-BC)sin(2φ)]+(nπ)2,ImΦ(Reiφ,n)=R2C2+D2[(AC+BD)sin(2φ)-(AD-BC)cos(2φ)],(3.32)AC+BD=aRαcos⁡(αφ)+cRγcos⁡(γφ)+a2bR2α-βcos⁡(βφ)+a2cb2R2(α-β)+γcos⁡(γφ)+2acbRα+γ-βcos⁡((α-β)φ)cos⁡(γφ),AD-BC=aRαsin(αφ)+cRγsin(γφ)+a2bR2α-βsin(βφ)+a2cb2R2(α-β)+γsin(γφ)+2acbRα+γ-βcos⁡((α-β)φ)sin(γφ).
Let(R,φ) be a solution to Φ(s,n)=0 (or equivalently, ReΦ=ImΦ=0). Then changing φ→-φ, we again obtain that ReΦ=ImΦ=0, which implies that solutions of (3.27) are complex conjugated.
Further, we will prove that the functionΦ has no zeros in the half plane args∈[0,π/2]. For that purpose we will use the argument principle. Recall, if Φ is an analytic function inside and on a regular closed curve c and nonzero on c, then the number of zeros of Φ is given by NZ=(1/2π)ΔargΦ(s). Let γ=γa∪γb∪γc be parametrized as
(3.33)γa:s=x,x∈[0,R],γb:s=Reiφ,φ∈[0,π2],γc:s=xeiπ/2,x∈[0,R],
and let R→∞. Along γa function Φ becomes a real-valued function, hence ΔargΦ(s,n)=0. Along γb we have that AC+BD,AD-BC≥0 for φ∈[0,π], since sin(ηφ),cos(ηφ)>0, for η∈{α,β,γ,α-β}⊂(0,1/2), and φ∈[0,π]. Therefore, (3.31) imply that
(3.34)ReΦ(Reiφ,n)>0,φ∈[0,π4],ImΦ(Reiφ,n)>0,φ∈[π4,π2]
(inequalities at boundary points are easily checked by inserting them into (3.31)). Along γc, we obtain
(3.35)ReΦ(xeiπ/2,n)=-R2C2+D2|R=x,φ=π/2(AC+BD)|R=x,φ=π/2+(nπ)2,(3.36)ImΦ(xeiπ/2,n)=R2C2+D2|R=x,φ=π/2(AD-BC)|R=x,φ=π/2>0.
From (3.34) and (3.36), we may now conclude that along γb and γc,
(3.37)ΔargΦ(s,n)=0.
Indeed, this follows from the following. Along γb, in the case φ∈[0,π/4], ImΦ(Reiφ,n) can change its sign, but ReΦ(Reiφ,n)>0, while for φ∈[π/4,π/2], ImΦ(Reiφ,n)>0. Along γc, ImΦ(xeiπ/2,n)>0. Therefore, there is no change in the argument of Φ along the whole γ, which implies that Φ has no zeros for φ∈[0,π/2]. This further implies that (3.27) has no solutions in the right complex half plane, since its solutions are complex conjugated.
In order to prove that for fixedn∈ℕ (3.27) has one solution (and its complex conjugate), we again use function Φ and the argument principle. Consider now the contour Γ=ΓA∪ΓB∪ΓC parametrized by
(3.38)ΓA:s=xeiπ/2,x∈[0,R],ΓB:s=Reiφ,φ∈[π2,π],ΓC:s=xeiπ,x∈[0,R],
and let R→∞. Along ΓA real and imaginary parts of Φ are given by (3.35) and (3.36). Along ΓB, using (3.31), we conclude that
(3.39)ReΦ(Reiφ,n)<0,φ∈[π2,3π4],forR⟶∞,andnfixed,ImΦ(Reiφ,n)<0,φ∈[3π4,π].
Along ΓC we have
(3.40)ReΦ(Reiφ,n)=R2C2+D2|R=x,φ=π(AC+BD)|R=x,φ=π+(nπ)2>0,ImΦ(Reiφ,n)=-R2C2+D2|R=x,φ=π(AD-BC)|R=x,φ=π<0.
Along ΓA we have that ImΦ>0, while ReΦ changes its sign (since ReΦ(0,n)=(nπ)2 and limx→∞ReΦ(xeiπ/2,n)=-∞, for fixed n∈ℕ). Along the part of ΓB where φ∈[3π/4,π], and along ΓC, ImΦ<0. Also, limR→∞ReΦ(Rei3π/4,n)=-∞ and ImΦ(Rei3π/4,n)<0. This implies that the argument of Φ changes from 0 to 2π.
As a conclusion one obtains that alongΓ(3.41)ΔargΦ(s,n)=2π,
which implies, by the argument principle, that function Φ has exactly one zero in the upper left complex plane, for each fixed n∈ℕ. Since the zeros of Φ are complex conjugated, it follows that Φ also has one zero in the lower left complex plane, for each fixed n∈ℕ.In the next proposition we examine behavior of simple polesPsn(±), n∈ℕ.Proposition 3.4.
SolutionsPsn(±), n∈ℕ, of (3.27), are such that
(3.42)Re(Psn(±))=Rcos⁡φ~c(nπ)22-γcos⁡(π2-γ)<0,Im(Psn(±))=Rsinφ~±c(nπ)22-γsin(π2-γ),
as n→∞.Proof.
Let us square (3.27) and insert Psn(±)=Reiφ, φ∈(-π,π). Then, after separation of real and imaginary parts, we obtain
(3.43)R2cos⁡(2φ)Re(M2(Reiφ))-R2sin(2φ)Im⁡(M2(Reiφ))=-(nπ)2,(3.44)R2sin(2φ)Re(M2(Reiφ))+R2cos⁡(2φ)Im⁡(M2(Reiφ))=0.
Using notation from the proof of Proposition 3.3 we can write
(3.45)Re(M2(Reiφ))=AC+BDC2+D2,Im⁡(M2(Reiφ))=BC-ADC2+D2.
Letting R→∞, one has
(3.46)Re(M2(Reiφ))~(a2c/b2)R2(α-β)+γcos⁡(γφ)(a2c2/b2)R2(α-β)+2γ=1cRγcos⁡(γφ),Im⁡(M2(Reiφ))~-(a2c/b2)R2(α-β)+γsin(γφ)(a2c2/b2)R2(α-β)+2γ=-1cRγsin(γφ).
It now follows from (3.44) and (3.46) that
(3.47)tg(2φ)=-Im(M2(Reiφ))Re(M2(Reiφ))~tg(γφ)⟹sin((2-γ)φ)cos⁡(2φ)cos⁡(γφ)~0⟹φ~±π2-γ.
Inserting (3.47) into (3.46), and subsequently into (3.43), we obtain
(3.48)R2-γccos⁡(2π2-γ)cos⁡(γπ2-γ)+R2-γcsin(2π2-γ)sin(γπ2-γ)~-(nπ)2,R~c(nπ)22-γ.
Thus, real and imaginary parts of Psn(±), as R→∞, are as claimed.Proposition 3.5.
Letp∈(0,s0), s0>0. Then
(3.49)M(p±iR)~1cRγe∓iγπ/4
as R→∞.Proof.
Setμ:=p2+R2 and ν:=arctan(±R/p). Then μ~R and ν~±π/2, as R→∞. By (3.46), we have
(3.50)M(μeiν)~1cμγe∓iγν/4,μ⟶∞,
as claimed.In order to obtain the explicit form of solutionu to initial-boundary value problem (3.4), (3.5), and (3.6), it remains to calculate function P.Theorem 3.6.
The solutionu to initial-boundary value problem (3.4), (3.5), and (3.6) is given by (3.17), that is, u(x,t)=Υ(t)*P(x,t), where P takes the form
(3.51)P(x,t)=12πi∫0∞(sinh⁡(xqM(qe-iπ))sinh⁡(qM(qe-iπ))-sinh⁡(xqM(qeiπ))sinh⁡(qM(qeiπ)))e-qtdq+∑n=1∞[Res(P̃(x,s)est,Psn(+))+Res(P̃(x,s)est,Psn(-))],t>0.
The residues at simple poles Psn(±), n∈ℕ, are given by
(3.52)Res(P̃(x,s)est,Psn(±))=[sinh⁡(xsM(s))(d/ds)[sinh⁡(sM(s))]est]s=Psn(±).Proof.
FunctionP(x,t), x∈[0,1], t>0, will be calculated by integration over a suitable contour. The Cauchy residues theorem yields
(3.53)∮ΓP̃(x,s)estds=2πi∑n=1∞[Res(P̃(x,s)est,Psn(+))+Res(P̃(x,s)est,Psn(-))],
where Γ=Γ1∪Γ2∪Γ3∪Γɛ∪Γ4∪Γ5∪Γ6∪γ0 is such a contour that all poles lie inside the contour Γ (see Figure 1).    First we show that the series of residues in (3.51) is convergent. By Proposition 3.4, Psn(±) are simple poles of P̃, and therefore also simple poles of estP̃. The residues in (3.53) can be calculated as it is given in (3.52), so
(3.54)Res(P̃(x,s)est,Psn(±))=[1(d/ds)(sM(s))sinh⁡(xsM(s))cosh⁡(sM(s))est]s=Psn(±).
Using (3.27) we have that sinh(xsM(s))=±isin(xnπ), while cosh(sM(s))=(-1)n. Also, a calculation gives that (d/ds)(sM(s))=M(s)+M(s)·A(s), where
(3.55)A(s):=1+(a/b)(α-β)sα-β2(1+(a/b)sα-β)-aαsα+cγsγ+(ac/b)(α+γ-β)sα+γ-β2(asα+csγ+(ac/b)sα+γ-β).
Take now that Psn(±)=Re±iφ. Then
(3.56)Res(P̃(x,s)est,Psn(±))=(-1)nsin(nπx)nπReRtcosφe±i(φ+Rtsinφ)A(Re±iφ).
Let n→∞. Then also |Psn(±)|→∞, that is, R→∞, and
(3.57)|A(Re±iφ)|⟶1-γ2.
Further, it follows from Proposition 3.4 that
(3.58)Re(sPn(±))=Rcos⁡φ~c(nπ)22-γcos⁡(π2-γ)≤-Cn,forsomeC>0,
and by (3.48), R/n~cπ22-γ·nγ/(2-γ). Therefore, as n→∞,
(3.59)|Res(P̃(x,s)est,Psn(+))+Res(P̃(x,s)est,Psn(-))|≤|sin(nπx)nπReRtcos⁡φA(Re+iφ)|+|sin(nπx)nπReRtcos⁡φA(Re-iφ)|≤1πRne-Cnt(1|A(Re+iφ)|+1|A(Re-iφ)|)~4π(2-γ)cπ22-γ⋅nγ/(2-γ)e-Cnt,
which implies the convergence of the sum of residues in (3.53).
It remains to calculate the integral overΓ in (3.53). Consider the integral along contour Γ1:s=p+iR, s0>p>0. Then
(3.60)|∫Γ1P̃(x,s)estds|≤∫0s0|P̃(x,p+iR)||e(p+iR)t|dp.
Let R→∞. In order to estimate |P̃(x,p±iR)|, using Proposition 3.5, we write
(3.61)M(p±iR)~v±iw,v=1cRγcos⁡(γπ4),w=-1cRγsin(γπ4).
Then
(3.62)|P̃(x,p±iR)|~|sinh⁡[x(pv-Rw)±ix(pw+Rv)]sinh⁡[(pv-Rw)±i(pw+Rv)]|≤ex(pv-Rw)+e-x(pv-Rw)|epv-Rw-e-(pv-Rw)|=e-(1-x)(pv-Rw)1+e-2x(pv-Rw)|1-e-2(pv-Rw)|⟶0,asR⟶∞.
The above convergence is valid since
(3.63)pv-Rw=p1cRγcos⁡(γπ4)+R1cRγsin(γπ4)⟶∞,asR⟶∞.
Therefore, according to (3.62), we have
(3.64)lim⁡R→∞|∫Γ1P̃(x,s)estds|=0.
The similar argument is valid for the integral along Γ6, thus
(3.65)lim⁡R→∞|∫Γ6P̃(x,s)estds|=0.
Next, we consider the integral along contour Γ2:s=Reiφ, π/2<φ<π:
(3.66)|∫Γ2P̃(x,s)estds|≤∫π/2πR|e-R(1-x)eiφM(Reiφ)||1-e-2xReiφM(Reiφ)1-e-2ReiφM(Reiφ)|eRtcos⁡φdφ.
Since sM(s)→∞ as |s|→∞ (see Proposition 3.2 (ii)) and cosφ≤0 for φ∈[π/2,π], we have
(3.67)lim⁡R→∞|∫Γ2P̃(x,s)estds|≤lim⁡R→∞∫π/2πR|e-R(1-x)eiφM(Reiφ)|eRtcos⁡φdφ=0.
The similar argument is valid for the integral along Γ5, thus
(3.68)lim⁡R→∞|∫Γ5P̃(x,s)estds|=0.
Since sM(s)→0 as |s|→0 (see Proposition 3.2 (ii)), the integration along contour Γɛ:s=ɛeiφ, π>φ>-π, gives
(3.69)lim⁡ɛ→0|∫ΓɛP̃(x,s)estds|≤lim⁡ɛ→0∫π-πɛ|sinh⁡(xɛeiφM(ɛeiφ))sinh⁡(ɛeiφM(ɛeiφ))|eɛtcos⁡φdφ=0.
Integrals along parts of contour Γ3:s=qeiπ, R>q>ɛ, Γ4:s=qe-iπ, ɛ<q<R, and γ0:s=s0+ir, -R<r<R, give
(3.70)lim⁡R→∞ɛ→0∫Γ3P̃(x,s)estds=∫0∞sinh⁡(xqM(qeiπ))sinh⁡(qM(qeiπ))e-qtdq,lim⁡R→∞ɛ→0∫Γ4P̃(x,s)estds=-∫0∞sinh⁡(xqM(qe-iπ))sinh⁡(qM(qe-iπ))e-qtdq,lim⁡R→∞∫γ0P̃(x,s)estds=2πiP(x,t).
Equation (3.51) now follows from (3.53).Figure 1
Integration contourΓ.Corollary 3.7.
In the case of stress relaxation, that is, whenΥ(t)=Υ0H(t), Υ0>0, t∈ℝ, the solution takes the form
(3.71)uH(x,t)=Υ0H(t)*P(x,t),x∈[0,1],t∈R.
We will numerically examine it in the sequel.


## 3.2.2. Determination of Stressσ in the Case of Prescribed Displacement Υ

In Section3.1 we determined stress σ (cf. (3.21)) that is a solution to (3.4), (3.5), and (3.6). In order to obtain an explicit form of σ we need to calculate function T. As in previous Section 3.2.1 it will be done by inversion of the Laplace transform of T̃.FunctionT̃, which is given by (3.20), is analytic on the complex plane except the branch cut (-∞,0], and has simple poles at the same points as P̃, that is, Psn(±), n∈ℕ.Using the similar arguments as in the proof of Theorem3.6, one can prove the following theorem.Theorem 3.8.
The solutionσ to initial-boundary value problem (3.4), (3.5), and (3.6) is given by (3.21), that is, σ(x,t)=d/dt(Υ(t)*T(x,t)), where T takes the form
(3.72)T(x,t)=12πi∫0∞(cosh⁡(xqM(qeiπ))M(qeiπ)sinh⁡(qM(qeiπ))-cosh⁡(xqM(qe-iπ))M(qe-iπ)sinh⁡(qM(qe-iπ)))e-qtdq+∑n=1∞[Res(T̃(x,s)est,Psn(+))+Res(T̃(x,s)est,Psn(-))],t>0.
The residues at simple poles Psn(±), n∈ℕ, are given by
(3.73)Res(T̃(x,s)est,Psn(±))=[cosh⁡(xsM(s))M(s)(d/ds)[sinh⁡(sM(s))]est]s=Psn(±).Corollary 3.9.
Similarly, in the case of stress relaxationΥ=Υ0H we obtain the solution
(3.74)σH(x,t)=Υ0T(x,t),x∈[0,1],t>0.In order to check our results for large times, we compare them with the quasistatic case. In the quasistatic case one uses only the constitutive equation (3.4)3, that is, the dynamics of the process is neglected. Taking the Laplace transform of the constitutive equation we obtain (3.9)2, and define the relaxation modulus G via its Laplace transform, as follows:
(3.75)G(t)∶=L-1[G̃(s)](t),t>0,G̃(s)∶=σ̃(QS)(s)ɛ̃(QS)(s)∶=asα+csγ+(ac/b)sα+γ-β1+(a/b)sα-β,s∈C∖(-∞,0].
Then the stress in the quasistatic case is
(3.76)σ(QS)=G*ɛ(QS).
Following the proof of Theorem 3.6, we obtain that (3.77)G(t)=12πi∫0∞[G̃(qe-iπ)-G̃(qeiπ)]e-qtdt.In the quasistatic case it holds that(3.78)u(x,t)=x⋅u(1,t),x∈[0,1]t>0,
and consequently by (3.4), (3.79)ɛ(x,t)=u(1,t)=:ɛ(QS)(t),x∈[0,1]t>0.
Since according to boundary condition (3.6), u(1,t)=Υ(t), it follows from (3.76) that (3.80)σ(QS)=Υ*G,
which, in the case of stress relaxation, that is, when Υ=Υ0H, Υ0>0, becomes(3.81)σH(QS)=Υ0H*G.


## 3.2.3. Determination of Displacementu in the Case of Prescribed Stress Σ

In Section3.1 we determined displacement u (cf. (3.25)) that is a solution to (3.4), (3.5), and (3.7). As above, we now want to find u explicitly, by calculating the inverse Laplace transform of Q̃.FunctionQ̃ given by (3.24) is analytic on the complex plane except the branch cut (-∞,0] and has isolated singularities at solutions Qsn(±) of the equation
(3.82)cosh⁡(sM(s))=0thatis,sM(s)=±2n+12iπ,n∈N0.We state a proposition that is analogous to Propositions3.3 and 3.4. The proof is omitted since it follows the same lines as those of Propositions 3.3 and 3.4.Proposition 3.10.
(i) There are infinitely many complex conjugated solutionsQsn(±), n∈ℕ0, of (3.82), which all lie in the left complex half plane. Moreover, each Qsn(±), n∈ℕ0, is a simple pole.
(ii) SolutionsQsn(±), n∈ℕ0, of (3.82), are such that
(3.83)Re(sQn(±))=Rcos⁡φ~c(2n+12π)22-γcos⁡(π2-γ)<0,Im(sQn(±))=Rsinφ~±c(2n+12π)22-γsin(π2-γ),
as n→∞.In the following theorem we calculate explicitly displacementu.Theorem 3.11.
The solutionu to initial-boundary value problem (3.4), (3.5), and (3.7) is given by (3.25), that is, u(x,t)=Σ(t)*Q(x,t), where Q takes the form
(3.84)Q(x,t)=12πi∫0∞(M(qe-iπ)sinh⁡(xqM(qe-iπ))cosh⁡(qM(qe-iπ))-M(qeiπ)sinh⁡(xqM(qeiπ))cosh⁡(qM(qeiπ)))e-qtqdq+∑n=0∞[Res(Q̃(x,s)est,Qsn(+))+Res(Q̃(x,s)est,Qsn(-))],t>0.
The residues at simple poles Qsn(±), n∈ℕ0, are given by
(3.85)Res(Q̃(x,s)est,Qsn(±))=[1sM(s)sinh⁡(xsM(s))(d/ds)[cosh⁡(sM(s))]est]s=Qsn(±).Corollary 3.12.
The case of creep is described by the boundary conditionΣ(t)=Σ0H(t), Σ0>0, t∈ℝ, in which displacement u, given by (3.25), reads
(3.86)u(t)=Σ0H(t)*Q(x,t),x∈[0,1],t∈R.Similarly as in Section3.2.2, we examine the quasistatic case that corresponds to displacement u. Again, using only the constitutive equation (3.4)3 and its Laplace transform (3.9)2, we define the creep compliance J via its Laplace transform as
(3.87)J(t)∶=L-1[J̃(s)](t),t>0,J̃(s)∶=ɛ̃(QS)(s)σ̃(QS)(s)∶=1+(a/b)sα-βasα+csγ+(ac/b)sα+γ-β,s∈C∖(-∞,0].
Strain measure in the quasistatic case now equals(3.88)ɛ(QS)=J*σ(QS).
Following the proof of Theorem 3.6 one obtains (3.89)J(t)=12πi∫0∞[J̃(qe-iπ)-J̃(qeiπ)]e-qtdt.
In the quasistatic case it holds that (3.90)u(x,t)=x⋅u(1,t)=x⋅u(QS)(t),x∈[0,1],t>0,
and consequently by (3.4),
(3.91)ɛ(QS)=u(QS).
Also, σ(QS)(t):=σ(1,t)=Σ(t), t>0, which is the boundary condition (3.7), hence by (3.88) and (3.91), (3.92)u(QS)=Σ*J.
In the case of creep (cf. Corollary 3.12) we have(3.93)uH(QS)=Σ0H*J.


## 3.3. Numerical Examples

In this subsection we give several numerical examples of displacementuH and stress σH, given by (3.71) and (3.74), respectively, which correspond to the case of stress relaxation, and examine solutions (3.86), which correspond to displacement u in the case of creep. In addition, we investigate solutions uH, σH, and u for different orders of fractional derivatives.Figure2 presents displacements in a stress relaxation experiment, determined according to (3.71), for three different positions. Parameters in (3.71) are chosen as follows: Υ0=1, a=0.2, b=0.6, c=0.45, α=0.3, β=0.1, and γ=0.4. From Figure 2, one sees that the displacements in the case of stress relaxation show damped oscillatory character and that they tend to a constant value for large times, namely, limt→∞uH(x,t)=x, x∈[0,1]. Figure 3 presents the same displacements as Figure 2, but close to initial time instant. It is evident that there is a delay in displacement that increases as the point is further from the end where the prescribed displacement is applied. This is a consequence of the finite wave propagation speed.Figure 2
DisplacementsuH(x,t) in a stress relaxation experiment as functions of time t at x∈{0.25,0.5,0.75} for t∈(0,10).Figure 3
DisplacementsuH(x,t) in a stress relaxation experiment as functions of time t at x∈{0.25,0.5,0.75} for t∈(0,5).Figure4 presents the displacements uH of the points close to the end of the rod where the sudden but afterwards constant displacement is applied (i.e., u(1,t)=H(t), t>0). One sees that the amplitudes of these points deform in shape so that they do not exceed the prescribed value of the displacement of the rod's free end.Figure 4
DisplacementsuH(x,t) in a stress relaxation experiment as functions of time t at x∈{0.75,0.8,0.85,0.9,0.95} for t∈(0,5).In order to examine the influence of the orders of fractional derivatives in the constitutive equation (1.4) on the displacement uH and stress σH in a stress relaxation experiment, we plot the displacement uH and stress σH obtained by (3.71) and (3.74) for the following sets of parameters (α,β,γ)∈{(0.1,0.05,0.15),(0.3,0.1,0.4),(0.45,0.4,0.49)}, while we fix x=0.5 and leave other parameters as before. In the case of the first set, the constitutive equation (1.4) describes a body in which the elastic properties are dominant, since the orders of the fractional derivatives of stress and strain are close to zero, that is, the fractional derivatives of stress and strain almost coincide with the stress and strain. This is also evident from Figure 5, since the oscillations of the point x=0.5 for the first set of parameters vanish quite slowly comparing to the second set and in particular comparing with the third set of parameters. Note that the third set of parameters describes a body in which the fluid properties of the fractional type dominate, since in the constitutive equation (1.4), we have the low-order derivative of stress (almost the stress itself), while almost all the derivatives of strain are of order 0.5. Figure 5 also shows that the dissipative properties of a material grow as the orders of the fractional derivatives increase. Figure 6 shows that the delay in displacement depends on the order of the fractional derivative, so that a material which has dominant elastic properties (the first set) has the longest delay, compared to a material with the dominant fluid properties (the third set).Figure 5
DisplacementsuH(x,t) in a stress relaxation experiment as functions of time t at x=0.5 for t∈(0,20).Figure 6
DisplacementsuH(x,t) in a stress relaxation experiment as functions of time t at x=0.5 for t∈(0,8).Figure7 presents stresses σH in the case of stress relaxation, determined according to (3.74), for different points of the rod. Parameters are the same as in the previous case. Also, Figure 7 presents the quasistatic curve σH(QS) obtained by (3.81).Figure 7
StressesσH(x,t) and σH(QS)(t) in a stress relaxation experiment as functions of time t at x∈{0.1,0.2,0.3,0.4,0.5,0.6,0.7} for t∈(0,5).Stresses, as it can be seen from Figure7, show damped oscillatory character and for large times, in each point x∈[0,1], tend to the quasistatic curve, that is, to the same value. Eventually, the stresses in all points of the rod tend to zero, namely, limt→∞σH(x,t)=0, x∈[0,1]. From Figure 7 it is evident that as further the point is from the rod-free end the greater is the delay. This is again the consequence of the finite wave speed.Figures8 and 9 present the stresses of the points close to the free end. One notices from Figures 7 and 8 that as the point is closer to the end whose displacement is prescribed, the peak of stress is higher and the peak's width is smaller. Figure 9 presents the compressive phase in the stress relaxation process and the quasistatic curve as well.Figure 8
StressesσH(x,t) in a stress relaxation experiment as functions of time t at x∈{0.7,0.8,0.9} for t∈(0,0.5).Figure 9
StressesσH(x,t) and σH(QS)(t) in a stress relaxation experiment as functions of time t at x∈{0.7,0.8,0.9} for t∈(0,6).Again, we fix the midpoint of the rod and investigate the influence of the change of orders of the fractional derivatives (for the same three sets as before) on the stressσH in a stress relaxation experiment. From Figure 10 one notices that there is a stress relaxation in a material regardless of the order of fractional derivatives. However, the relaxed stress depends on the order of the derivatives, as expected from the analysis of the stress σH(QS) in the quasistatic case, given by (3.81). The material with the dominant elastic properties (the first set) relaxes to the highest stress, and as the fluid properties of the material become more and more dominant, the relaxed stress decreases. Again, the oscillations of the value of the stress for the first set (material with dominant elastic properties) are the least damped compared to the second and third set. Figure 11 shows that the conclusion about the dependence of delay on the orders of fractional derivatives drown earlier holds.Figure 10
StressesσH(x,t) in a stress relaxation experiment as functions of time t at x=0.5 for t∈(0,10).Figure 11
StressesσH(x,t) in a stress relaxation experiment as functions of time t at x=0.5 for t∈(0,4.5).Figure12 presents displacements u in the creep experiment, determined according to (3.86), for four different points. Parameters are the same as in the previous cases (except that in this case instead of Υ0=1 we have Σ0=1). For large times, as it can be seen from Figure 12, the displacement curves are monotonically increasing. This indicates that we deal with the viscoelastic fluid. Figure 12 also shows good agreement between the displacements obtained by the dynamic model (displacement is given by (3.86)) and the quasistatic model (displacement is given by (3.93)). Figure 13 presents the same displacements as Figure 12, but close to initial time instant and it is evident that, again, there is a delay, due to the finite speed of wave propagation.Figure 12
Displacementsu(x,t) and uH(QS)(t) in a creep experiment as functions of time t at x∈{0.25,0.5,0.75,1} for t∈(0,40).Figure 13
Displacementsu(x,t) in creep experiment as functions of time t at x∈{0.25,0.5,0.75,1} for t∈(0,10).In order to examine the dependence of the displacementu in the case of creep experiment, given by (3.86), on the orders of fractional derivatives in the constitutive equation (1.4), we again fix the point x=0.5 of the rod and plot the displacement u for the same set of values of α, β, and γ as before. Figure 14 clearly shows that all of the materials exhibits creep, but the displacements do not tend to a constant value. However, the material which has the dominant elastic properties (the first set) creeps slower comparing to the material with the dominant fluid properties (the third set). Figure 15 again supports the conclusion that the more elastic the material is, the greater the delay is.Figure 14
Displacementsu(x,t) in creep experiment as functions of time t at x=0.5 for t∈(0,40).Figure 15
Displacementsu(x,t) in creep experiment as functions of time t at x=0.5 for t∈(0,10).



---

*Source: 975694-2011-06-02.xml*