Fractional Burgers models in creep and stress relaxation tests
FFractional Burgers models increep and stress relaxation tests
Aleksandar S. Okuka ∗ , Duˇsan Zorica † March 28, 2019
Abstract
Classical and thermodynamically consistent fractional Burgers models are examined in creep and stressrelaxation tests. Using the Laplace transform method, the creep compliance and relaxation modulus areobtained in integral form, that yielded, when compared to the thermodynamical requirements, the narrowerrange of model parameters in which the creep compliance is a Bernstein function while the relaxation modulusis completely monotonic. Moreover, the relaxation modulus may even be oscillatory function with decreasingamplitude. The asymptotic analysis of the creep compliance and relaxation modulus is performed near theinitial time-instant and for large time as well.
Key words : thermodynamically consistent fractional Burgers models, creep compliance and relaxationmodulus, Bernstein and completely monotonic functions
The classical Burgers model (cid:18) a dd t + a d d t (cid:19) σ ( t ) = (cid:18) b dd t + b d d t (cid:19) ε ( t ) , (1)where σ and ε denote stress and strain, that are functions of time t > , while a , a , b , b > a a ≤ b b ≤ a (2)on model parameters appearing in the classical Burgers model (1) are derived in [16] by requiring non-negativityof the storage and loss modulus that is obtained as a consequence of the dissipativity inequality in the steadystate regime, see [5]. Figure 1: Rheological representation of the classical Burgers model.The fractional generalization of Burgers model is derived in [16] by considering the Scott-Blair (fractional)element instead of the dash-pot element in the rheological representation from Figure 1, with the orders offractional differentiation corresponding to the fractional elements and their sums being replaced by the arbitraryorders of fractional derivatives α, β, µ, γ, and ν. Such obtained fractional Burgers model takes the form (cid:16) a D αt + a D βt + a D γt (cid:17) σ ( t ) = ( b D µt + b D νt ) ε ( t ) , (3) ∗ Department of Mechanics, Faculty of Technical Sciences, University of Novi Sad, Trg D. Obradovi´ca 6, 21000 Novi Sad, Serbia,[email protected] † Mathematical Institute, Serbian Academy of Arts and Sciences, Kneza Mihaila 36, 11000 Belgrade, Serbia, du-san [email protected] and Department of Physics, Faculty of Sciences, University of Novi Sad, Trg D. Obradovi´ca 4, 21000Novi Sad, Serbia a r X i v : . [ phy s i c s . c l a ss - ph ] M a r here the model parameters are denoted by a , a , a , b , b > , α, β, µ ∈ [0 , , with α ≤ β, and γ, ν ∈ [1 , , while D ξt denotes the operator of Riemann-Liouville fractional derivative of order ξ ∈ [ n, n + 1] , n ∈ N , definedby D ξt y ( t ) = d n +1 d t n +1 (cid:18) t − ( ξ − n ) Γ (1 − ( ξ − n )) ∗ y ( t ) (cid:19) , t > , see [11], with ∗ denoting the convolution in time: f ( t ) ∗ g ( t ) = (cid:82) t f ( u ) g ( t − u ) d u, t > . Thermodynamical consistency analysis of the fractional Burgers model (3), conducted in [16] by the useof storage and loss modulus non-negativity requirement, implied that the orders of fractional derivatives frominterval [1 ,
2] cannot be independent of the orders of fractional derivatives from interval [0 , , and this led toformulation of eight thermodynamically consistent fractional Burgers models.Two classes of thermodynamically consistent fractional Burgers models are distinguished according to theorders of fractional derivatives acting on stress and strain. The first class of fractional Burgers models consistsof models having different orders of fractional differentiation of stress and strain from both intervals [0 ,
1] and[1 , . Namely, in the case of Model I the highest order of fractional differentiation of stress is γ ∈ [0 , , whilethe highest order of fractional differentiation of strain is ν = µ + κ ∈ [1 , , with 0 ≤ α ≤ β ≤ γ ≤ µ ≤ κ ∈ { α, β, γ } , while in the case of Models II - V one has 0 ≤ α ≤ β ≤ µ ≤ ≤ γ ≤ ν ≤ , with ( γ, ν ) ∈ { (2 α, µ + α ) , ( α + β, µ + α ) , ( α + β, µ + β ) , (2 β, µ + β ) } , see the unified constitutive equation(33) below. Likewise the classical Burgers model (1), Models VI - VIII contain the same orders of fractionalderivatives acting on stress and strain from both intervals [0 ,
1] and [1 , , such that for Models VI and VII infractional Burgers model (3) one has µ = β ∈ [0 ,
1] and γ = ν = β + η ∈ [1 , , with η ∈ { α, β } , while ModelVIII is obtained from (3) for β = α, ¯ a = a + a , ¯ a = a , µ = α, γ = ν = 2 α, see the unified constitutiveequation (34) below.Models I - VIII, along with corresponding thermodynamical constraints, are listed below. Model I : (cid:16) a D αt + a D βt + a D γt (cid:17) σ ( t ) = (cid:0) b D µt + b D µ + κ t (cid:1) ε ( t ) , (4)0 ≤ α ≤ β ≤ γ ≤ µ ≤ , ≤ µ + κ ≤ α, b b ≤ a i cos ( µ − κ ) π (cid:12)(cid:12)(cid:12) cos ( µ + κ ) π (cid:12)(cid:12)(cid:12) , (5)with ( κ , i ) ∈ { ( α, , ( β, , ( γ, } ; Model II : (cid:16) a D αt + a D βt + a D αt (cid:17) σ ( t ) = (cid:0) b D µt + b D µ + αt (cid:1) ε ( t ) , (6)12 ≤ α ≤ β ≤ µ ≤ , a a (cid:12)(cid:12)(cid:12) sin ( µ − α ) π (cid:12)(cid:12)(cid:12) sin µπ ≤ b b ≤ a cos ( µ − α ) π (cid:12)(cid:12)(cid:12) cos ( µ + α ) π (cid:12)(cid:12)(cid:12) ; (7) Model III : (cid:16) a D αt + a D βt + a D α + βt (cid:17) σ ( t ) = (cid:0) b D µt + b D µ + αt (cid:1) ε ( t ) , (8)0 ≤ α ≤ β ≤ µ ≤ , α + β ≥ , a a (cid:12)(cid:12)(cid:12) sin ( µ − β − α ) π (cid:12)(cid:12)(cid:12) sin ( µ − β + α ) π ≤ b b ≤ a cos ( µ − α ) π (cid:12)(cid:12)(cid:12) cos ( µ + α ) π (cid:12)(cid:12)(cid:12) ; (9) Model IV : (cid:16) a D αt + a D βt + a D α + βt (cid:17) σ ( t ) = (cid:16) b D µt + b D µ + βt (cid:17) ε ( t ) , (10)0 ≤ α ≤ β ≤ µ ≤ , − α ≤ β ≤ − ( µ − α ) , a a (cid:12)(cid:12)(cid:12) sin ( µ − α − β ) π (cid:12)(cid:12)(cid:12) sin ( µ − α + β ) π ≤ b b ≤ a cos ( µ − β ) π (cid:12)(cid:12)(cid:12) cos ( µ + β ) π (cid:12)(cid:12)(cid:12) ; (11) Model V : (cid:16) a D αt + a D βt + a D βt (cid:17) σ ( t ) = (cid:16) b D µt + b D µ + βt (cid:17) ε ( t ) , (12)0 ≤ α ≤ β ≤ µ ≤ , ≤ β ≤ − ( µ − α ) , a a (cid:12)(cid:12)(cid:12) sin ( µ − β ) π (cid:12)(cid:12)(cid:12) sin µπ ≤ b b ≤ a cos ( µ − β ) π (cid:12)(cid:12)(cid:12) cos ( µ + β ) π (cid:12)(cid:12)(cid:12) . (13)2 odel VI : (cid:16) a D αt + a D βt + a D α + βt (cid:17) σ ( t ) = (cid:16) b D βt + b D α + βt (cid:17) ε ( t ) , (14)0 ≤ α ≤ β ≤ , α + β ≥ , a a ≤ b b ≤ a cos ( β − α ) π (cid:12)(cid:12)(cid:12) cos ( β + α ) π (cid:12)(cid:12)(cid:12) ; (15) Model VII : (cid:16) a D αt + a D βt + a D βt (cid:17) σ ( t ) = (cid:16) b D βt + b D βt (cid:17) ε ( t ) , (16)0 ≤ α ≤ β ≤ , ≤ β ≤ α , a a ≤ b b ≤ a | cos ( βπ ) | ; (17) Model VIII : (cid:0) a D αt + ¯ a D αt (cid:1) σ ( t ) = (cid:0) b D αt + b D αt (cid:1) ε ( t ) , (18)12 ≤ α ≤ , ¯ a ¯ a ≤ b b ≤ ¯ a | cos ( απ ) | . (19)Models I - VIII describe different mechanical behavior, that will be illustrated by examining the responsesin creep and stress relaxation tests. Recall, creep compliance (relaxation modulus) is the strain (stress) historyfunction obtained as a response to the imposed sudden and time-constant stress (strain), i.e., stress (strain)assumed as the Heaviside step function. The difference between models belonging to the first and second classreflects in the behavior they describe near the initial time-instant, since models of the first class predict zeroglass compliance and thus infinite glass modulus, while the glass compliance is non-zero implying the non-zeroglass modulus in the case of models belonging to the second class, see (41) and (57) below. For both modelclasses equilibrium compliance is infinite implying the zero equilibrium modulus.For both classes of fractional Burgers models as well as for the classical Burgers model, thermodynamicalrequirements will prove to be less restrictive than the conditions guaranteeing that the creep compliance is aBernstein function, while the relaxation modulus is a completely monotonic function. Therefore, if the modelparameters fulfill the thermodynamical requirements but not the restrictive ones, the creep compliance andrelaxation modulus may even not be a monotonic function. Moreover, conditions on model parameters guar-anteeing oscillatory behavior of the relaxation modulus having amplitudes decreasing in time will be obtained.Conditions guaranteeing that the creep compliance (relaxation modulus) is a Bernstein (completely monotonic)function and conditions guaranteeing the oscillatory behavior of the relaxation modulus are independent. Re-call, completely monotonic function is a positive, monotonically decreasing convex function, or more preciselyfunction f satisfying ( − n f ( n ) ( t ) ≥ , n ∈ N , while Bernstein function is a positive, monotonically increasing,concave function, or more precisely non-negative function having its first derivative completely monotonic.The properties of creep compliance as a Bernstein function and relaxation modulus as a completely mono-tonic function are discussed in [12], while [6] deal with the complete monotonicity of the relaxation modulicorresponding to distributed-order fractional Zener model. The review of creep compliances in the frequencydomain corresponding to the integer-order models of viscoelasticity is presented in [14].Classical (1), fractional and generalized fractional Burgers models in the form (cid:0) a D αt + a D αt (cid:1) σ ( t ) = η (cid:16) b D βt (cid:17) ˙ ε ( t ) , (cid:0) a D αt + a D αt (cid:1) σ ( t ) = η (cid:16) b D βt + b D βt (cid:17) ˙ ε ( t ) , were extensively used in modeling various fluid flow problems, see references in [16]. The thermodynamicalconsistency of the generalized fractional Burgers model for viscoelastic fluid (cid:0) a D αt + a D αt (cid:1) σ ( t ) = η (cid:16) b D βt + b D βt (cid:17) ˙ ε ( t ) , is examined in [7].Some form of the fractional Burgers constitutive equation is used in [8, 15, 20] for modeling asphalt con-crete mixtures according to the experimental data obtained in creep and creep-recovery experiments, while thefractional Burgers model (cid:0) a D νt + a D νt (cid:1) σ ( t ) = (cid:0) b D νt + b D νt (cid:1) ε ( t ) , is examined in creep and stress relaxation tests in [12, 13].Thermodynamical constraints on model parameters appearing in the fractional-order models of viscoelasticbody in the case when orders of fractional differentiation do not exceed the first order, along with the materialresponses in cases of damped oscillations and wave propagation are considered in [1, 2, 17, 18, 19]. The behaviorin creep and stress relaxation tests of a distributed-order fractional viscoelastic material with the inertial effectstaken into account is considered in [3, 4]. 3 Classical Burgers model: creep and stress relaxation tests
The aim is to investigate the behavior of thermodynamically consistent classical Burgers model (1) in creep andstress relaxation tests. It will be shown that the requirement for creep compliance to be the Bernstein func-tion and for relaxation modulus to be the completely monotonic function narrows down the thermodynamicalrestriction (2) to a a ≤ a (cid:32) − (cid:115) − a a (cid:33) < b b < a (cid:32) (cid:115) − a a (cid:33) ≤ a , provided that a a ≥ . Otherwise, still having the thermodynamical requirements fulfilled, the creep compliancewill prove to be a non-negative, monotonically increasing, but convex function, lying above its oblique asymptote,contrary to the case when it is a Bernstein function when the creep compliance is a concave function lying belowthe asymptote. If b b = a (cid:16) ± (cid:113) − a a (cid:17) , then the creep compliance is a linear, increasing function havingthe same form as the asymptote. In any case it starts from a finite value of strain and tends to infinity.Apart from being completely monotonic, the relaxation modulus will prove to be either non-monotonic func-tion having a negative minimum if a (cid:16) (cid:113) − a a (cid:17) < b b ≤ a , or a non-negative, monotonically decreasingfunction that may change its convexity if a a ≤ b b < a (cid:16) − (cid:113) − a a (cid:17) . Moreover, the relaxation modulus canbe an oscillatory function having decreasing amplitude if 1 ≤ a a < . In any case it starts from a finite value ofstress and tends to zero.The creep compliance in the form ε cr ( t ) = a − b b b (cid:16) − e − b b t (cid:17) + 1 b t + a b e − b b t , t ≥ , (20)having the glass compliance finite and equilibrium compliance infinite, i.e., ε ( g ) cr = ε cr (cid:0) + (cid:1) = a b and ε ( e ) cr = lim t →∞ ε cr ( t ) = ∞ , is obtained by: assuming σ = H ( H is the Heaviside function); applying the Laplace transform˜ f ( s ) = L [ f ( t )] ( s ) = (cid:90) ∞ f ( t ) e − st d t, Re s > , to the classical Burgers model (1) yielding˜ ε cr ( s ) = 1 s a s + a s b s + b s = a − b b b s + 1 b s + (cid:32) a b − a − b b b (cid:33) s + b b ; (21)and by the subsequent inversion of the Laplace transform in (21). Due to the thermodynamical restriction (2),the creep compliance (20) is a non-negative function and since˙ ε cr ( t ) = 1 b (cid:16) − e − b b t (cid:17) + a b b (cid:18) b b − a a (cid:19) e − b b t > , t ≥ , where ( · ) · = dd t ( · ) , it is also a monotonically increasing function, again due to (2). The second derivative ofcreep compliance (20) ¨ ε cr ( t ) = b b (cid:32) a b − a − b b b (cid:33) e − b b t , t ≥ , (22)is either non-negative or negative function for all t ≥ a b < b (cid:16) a − b b (cid:17) , since ε cr ( t ) > , ˙ ε cr ( t ) > , ¨ ε cr ( t ) < a b > b (cid:16) a − b b (cid:17) , then thecreep compliance is a non-negative, monotonically increasing, convex function. In the case a b = b (cid:16) a − b b (cid:17) , the creep compliance (20) becomes ε cr ( t ) = a − b b b + 1 b t = a b + 1 b t, t ≥ . (23)4he oblique asymptote of creep compliance (20) takes the form g ( t ) = a − b b b + 1 b t, as t → ∞ . (24)If ε ( g ) cr < g (0) , i.e., a b < b (cid:16) a − b b (cid:17) , ( ε ( g ) cr > g (0) , i.e., a b > b (cid:16) a − b b (cid:17) ), then the creep complianceis below (above) the oblique asymptote for all t ≥ , since it is a monotonically increasing concave (convex)function for all t ≥ . Note that if ε ( g ) cr = g (0) , i.e. , a b = b (cid:16) a − b b (cid:17) , then the creep compliance coincideswith the asymptote for all t ≥ . The relaxation modulus, corresponding to the classical Burgers model (1), is obtained by the Laplace trans-form method either as a non-oscillatory σ sr ( t ) = b a ν (cid:18)(cid:18) − ( µ − ν ) b b (cid:19) e − ( µ − ν ) t + (cid:18) ( µ + ν ) b b − (cid:19) e − ( µ + ν ) t (cid:19) , i.e. , (25) σ sr ( t ) = b a (cid:18) cosh ( νt ) + (cid:18) b b − µ (cid:19) sinh ( νt ) ν (cid:19) e − µt , t ≥ , (26)or as an oscillatory function having decreasing amplitude σ sr ( t ) = b a (cid:18) cos ( ωt ) + (cid:18) b b − µ (cid:19) sin ( ωt ) ω (cid:19) e − µt , t ≥ , (27)with µ = a a , ν = (cid:114) µ − a , and ω = ν i = (cid:114) a − µ . (28)The relaxation moduli (26) and (27) yield finite glass modulus and zero equilibrium modulus, i.e., σ ( g ) sr = σ sr (0) = 1 ε ( g ) cr = b a and σ ( e ) sr = lim t →∞ σ sr ( t ) = 1 ε ( e ) cr = 0 . The Laplace transform applied to the classical Burgers model (1), with the assumption ε = H, yields˜ σ sr ( s ) = 1 a b + b ss + a a s + a = 1 a b + b s ( s + µ − ν ) ( s + µ + ν )= 1 a (cid:18)(cid:18) b b − µb ν (cid:19) s + µ − ν + (cid:18) b − b − µb ν (cid:19) s + µ + ν (cid:19) , (29)so that (25) is obtained by the Laplace transform inversion in (29) if µ > a , while (27) follows from (26) if µ < a . If µ = a , then (29) yields˜ σ sr ( s ) = 1 a (cid:32) b s + µ + b − µb ( s + µ ) (cid:33) , i.e., σ sr ( t ) = b a (cid:18) (cid:18) b b − µ (cid:19) t (cid:19) e − µt , t ≥ . (30)Thermodynamical restriction (2), rewritten as a a ≥ , allows for both non-oscillatory and oscillatory formsof the relaxation moduli (26) and (27), since in (28), the condition ν = a (cid:16) a a − (cid:17) ≥ a a ≥
1, implies a a ≥ , while the condition ν = a (cid:16) a a − (cid:17) < ≤ a a < . The non-oscillatory relaxation modulus (25), transformed into σ sr ( t ) = b a ν (cid:18) ( µ − ν ) (cid:18) a ( µ + ν ) − b b (cid:19) e − ( µ − ν ) t + ( µ + ν ) (cid:18) b b − a ( µ − ν ) (cid:19) e − ( µ + ν ) t (cid:19) , (31)using a (cid:0) µ − ν (cid:1) = 1 from (28), is a completely monotonic function if the creep compliance (20) is a Bernsteinfunction. Namely, the condition a b < b (cid:16) a − b b (cid:17) solved with respect to b b yields a ( µ − ν ) < b b < a ( µ + ν ) , t ≥ σ sr ( t ) > σ sr have the alternatingsign due to the exponential function. The condition a a ≤ a ( µ − ν ) < b b < a ( µ + ν ) ≤ a for creep compliance (20) to be the Bernstein function and relaxation modulus (26) to be the completelymonotonic function is narrower than the thermodynamical restriction (2) since a ( µ + ν ) = a (cid:32) (cid:115) − a a (cid:33) ≤ a and a ( µ − ν ) = a (cid:32) − (cid:115) − a a (cid:33) = a (cid:32)
12 4 a a + ∞ (cid:88) k =2 (2 k − k k ! (cid:18) a a (cid:19) k (cid:33) ≥ a a , where the binomial formula √ − x = 1 − x − ∞ (cid:88) k =2 (2 k − k k ! x k , for 0 < x < , (32)is used.If the creep compliance is a linear function (23), i.e., if a b = b (cid:16) a − b b (cid:17) , i.e., b b = a ( µ − ν ) , or b b = a ( µ + ν ) , then the relaxation modulus (25) reduces to σ sr ( t ) = b a e − b a b t . Assume a ( µ + ν ) < b b ≤ a . The first derivative of relaxation modulus (31) reads˙ σ sr ( t ) = − b a ν ( µ + ν ) (cid:18) b b − a ( µ − ν ) (cid:19) (cid:32) e − νt − (cid:18) µ − νµ + ν (cid:19) b b − a ( µ + ν ) b b − a ( µ − ν ) (cid:33) e − ( µ − ν ) t , with 0 < (cid:16) µ − νµ + ν (cid:17) b b − a ( µ + ν ) b b − a ( µ − ν ) ≤ , due to the thermodynamical restriction a a ≤ b b , see (2). Since the exponentialfunction in the previous expression decreases from one to zero, at the time-instant t ∗ = 12 ν ln (cid:32)(cid:18) µ + νµ − ν (cid:19) b b − a ( µ − ν ) b b − a ( µ + ν ) (cid:33) the relaxation modulus has a minimum, since ˙ σ sr changes the sign from non-negative to negative at t ∗ . Thisfact, along with the finite glass and zero equilibrium modulus, implies that the relaxation modulus decreasesfrom σ ( g ) sr = b a to a negative minimum and further, being negative, increases to σ ( e ) sr = 0 . Assume a a ≤ b b < a ( µ − ν ) . The first derivative of relaxation modulus (31) is˙ σ sr ( t ) = − b a ν ( µ − ν ) (cid:18) a ( µ + ν ) − b b (cid:19) (cid:32) e νt − (cid:18) µ + νµ − ν (cid:19) a ( µ − ν ) − b b a ( µ + ν ) − b b (cid:33) e − ( µ + ν ) t < , since, for t > , e νt > < (cid:16) µ + νµ − ν (cid:17) a ( µ − ν ) − b b a ( µ + ν ) − b b ≤ , due to the thermodynamical restriction a a ≤ b b , see(2). Thus, the relaxation modulus, being a non-negative function, monotonically decreases from σ ( g ) sr = b a to σ ( e ) sr = 0 . However, the relaxation modulus may change its convexity, due to (cid:18) µ + νµ − ν (cid:19) a ( µ − ν ) − b b a ( µ + ν ) − b b ≥ b b ≤ a a − a a , that appears in the second derivative of the relaxation modulus¨ σ sr ( t ) = b a ν ( µ − ν ) (cid:18) a ( µ + ν ) − b b (cid:19) (cid:32) e νt − (cid:18) µ + νµ − ν (cid:19) a ( µ − ν ) − b b a ( µ + ν ) − b b (cid:33) e − ( µ + ν ) t . Fractional Burgers models: creep and stress relaxation tests
The fractional Burgers models I - VIII will be examined in creep and stress relaxation tests. Models I - V,respectively given by (4), (6), (8), (10), and (12), have zero glass compliance and thus infinite glass modulus,while Models VI - VIII, respectively given by (14), (16), and (18), behave similarly as the classical Burgersmodel (1), having non-zero glass compliance and thus glass modulus as well. In the case of both classical andfractional Burgers models, the equilibrium compliance is infinite and therefore equilibrium modulus is zero.The Laplace transform method will be used in calculating the creep compliance and relaxation modulus, sothat both will be obtained in an integral form using the definition of inverse Laplace transform, while the creepcompliance will additionally be expressed in terms of the Mittag-Leffler function. In Section 4, the integral formswill prove to be useful in showing that the thermodynamical requirements (5), (7), (9), (11), (13), (15), (17),and (19) allow wider range of model parameters than the range in which the creep compliance is a Bernsteinfunction and the relaxation modulus is a completely monotonic function, see (73), (80), (86), (92), (98), (104),(109), and (115). Hence, still being non-oscillatory, the creep compliance does not have to be a monotonicfunction, while the relaxation modulus may even be oscillatory function having decreasing amplitude if themodel parameters still fulfill the thermodynamical requirements.Models I - V, i.e., models having zero glass compliance, can be written in an unified manner as (cid:16) a D αt + a D βt + a D γt (cid:17) σ ( t ) = (cid:0) b D µt + b D µ + ηt (cid:1) ε ( t ) , (33)while Models VI - VIII, i.e., models having non-zero glass compliance, take the following unified form (cid:16) a D αt + a D βt + a D β + ηt (cid:17) σ ( t ) = (cid:16) b D βt + b D β + ηt (cid:17) ε ( t ) , (34)where in (33) the highest order of fractional differentiation of strain is µ + η ∈ [1 ,
2] with η ∈ { α, β } , whilethe highest order of fractional differentiation of stress is either γ ∈ [0 ,
1] in the case of Model I (4), with0 ≤ α ≤ β ≤ γ ≤ µ ≤ η = κ ∈ { α, β, γ } , or γ ∈ [1 ,
2] in the case of Models II - V, with 0 ≤ α ≤ β ≤ µ ≤ η, γ ) ∈ { ( α, α ) , ( α, α + β ) , ( β, α + β ) , ( β, β ) } , see (6), (8), (10), and (12), while in (34) one has0 ≤ α ≤ β ≤ β + η ∈ [1 , , with η = α in the case of Model VI (14) and η = β in the case of Model VII(16), while Model VIII (18) is obtained for η = β = α, ¯ a = a + a , and ¯ a = a . The creep compliances in complex domain˜ ε cr ( s ) = 1 s a s α + a s β + a s γ b s µ + b s µ + η = 1 s µ Ψ ( s ) b + b s η , (35)˜ ε cr ( s ) = 1 s a s α + a s β + a s β + η b s β + b s β + η = a b s (cid:32) s β ψ ( s ) b b + s η (cid:33) , (36)˜ ε cr ( s ) = 1 s a s α + ¯ a s α b s α + b s α = ¯ a b s (cid:32) s α ¯ ψ ( s ) b b + s α (cid:33) , (37)respectively corresponding to Models I - V, Models VI and VII, and Model VIII, with functions Ψ , ψ, and ¯ ψ, defined for s ∈ C by respective expressionsΨ ( s ) = 1 + a s α + a s β + a s γ , (38) ψ ( s ) = 1 a + a a s α + (cid:18) a a − b b (cid:19) s β , (39)¯ ψ ( s ) = 1¯ a + (cid:18) ¯ a ¯ a − b b (cid:19) s α , (40)are obtained by applying the Laplace transform to the unified constitutive equations (33) and (34) havingassumed that σ = H. Note that a a − b b ≥ ¯ a ¯ a − b b ≥ ε ( g ) cr = ε cr (0) = lim s →∞ ( s ˜ ε cr ( s )) = , for Models I - V, a b , for Models VI and VII, ¯ a b , for Model VIII, (41) ε ( e ) cr = lim t →∞ ε cr ( t ) = lim s → ( s ˜ ε cr ( s )) = ∞ , ε cr being the creep compliance in complex domain (35), or (36), or(37).Rewriting the creep compliances in complex domain (35), (36), and (37) respectively as˜ ε cr ( s ) = 1 b s η − (1+ µ + η ) b b + s η + a b s η − (1+ µ + η − α ) b b + s η + a b s η − (1+ µ + η − β ) b b + s η + a b s η − (1+ µ + η − γ ) b b + s η , ˜ ε cr ( s ) = a b s + 1 b s η − (1+ η + β ) b b + s η + a b s η − (1+ η + β − α ) b b + s η + a b (cid:18) a a − b b (cid:19) s η − (1+ η ) b b + s η , ˜ ε cr ( s ) = ¯ a b s + 1 b s α − (1+2 α ) b b + s α + ¯ a b (cid:18) ¯ a ¯ a − b b (cid:19) s α − (1+ α ) b b + s α , the creep compliance is expressed in terms of the two-parameter Mittag-Leffler function e η,ζ,λ ( t ) = t ζ − E η,ζ ( − λt η ) = L − (cid:20) s η − ζ s η + λ (cid:21) ( t ) , with E η,ζ ( z ) = ∞ (cid:88) k =0 z k Γ ( ηk + ζ ) , (42)see [10], as ε cr ( t ) = 1 b e η, µ + η, b b ( t ) + a b e η, µ + η − α, b b ( t ) + a b e η, µ + η − β, b b ( t ) + a b e η, µ + η − γ, b b ( t ) , (43) ε cr ( t ) = a b + 1 b e η, η + β, b b ( t ) + a b e η, η + β − α, b b ( t ) + a b (cid:18) a a − b b (cid:19) e η, η, b b ( t ) , (44) ε cr ( t ) = ¯ a b + 1 b e α, α, b b ( t ) + ¯ a b (cid:18) ¯ a ¯ a − b b (cid:19) e α, α, b b ( t ) , (45)respectively corresponding to Models I - V, Models VI and VII, and Model VIII.In the case of Models I - V, the integral form of creep compliance takes the form ε cr ( t ) = 1 π (cid:90) ∞ K ( ρ ) | b + b ρ η e i ηπ | − e − ρt ρ µ d ρ, (46)where K ( ρ ) = b K ( ρ ) + b ρ η K ( ρ ) , with (47) K ( ρ ) = sin ( µπ ) + a ρ α sin (( µ − α ) π ) + a ρ β sin (( µ − β ) π ) + a ρ γ sin (( µ − γ ) π ) ,K ( ρ ) = sin (( µ + η ) π ) + a ρ α sin (( µ + η − α ) π ) + a ρ β sin (( µ + η − β ) π ) + a ρ γ sin (( µ + η − γ ) π ) , while the creep compliance corresponding to Models VI and VII is given by ε cr ( t ) = a b + a b (cid:90) t f cr ( τ ) d τ , with f cr ( t ) = 1 π (cid:90) ∞ Q ( ρ ) (cid:12)(cid:12)(cid:12) b b + ρ η e i ηπ (cid:12)(cid:12)(cid:12) e − ρt ρ β d ρ, (48)where ε ( g ) cr = a b is the glass compliance (41), and where Q ( ρ ) = b b Q ( ρ ) + ρ η Q ( ρ ) , with (49) Q ( ρ ) = 1 a sin ( βπ ) + a a ρ α sin (( β − α ) π ) ,Q ( ρ ) = 1 a sin (( β + η ) π ) + a a ρ α sin (( β + η − α ) π ) + (cid:18) a a − b b (cid:19) ρ β sin ( ηπ ) . In the case of Model VIII, the creep compliance in integral form is obtained as ε cr ( t ) = ¯ a b + ¯ a b (cid:90) t ¯ f cr ( τ ) d τ , with ¯ f cr ( t ) = 1 π (cid:90) ∞ ¯ Q ( ρ ) (cid:12)(cid:12)(cid:12) b b + ρ α e i απ (cid:12)(cid:12)(cid:12) e − ρt ρ α d ρ, (50) If ζ = 1 in (42), the two-parameter Mittag-Leffler function reduces to the one-parameter Mittag-Leffler function e η,λ ( t ) = e η, ,λ ( t ) = E η ( − λt η ) = L − (cid:20) s η − s η + λ (cid:21) ( t ) , with E η ( z ) = ∞ (cid:88) k =0 z k Γ ( ηk + 1) . ε ( g ) cr = ¯ a b is the glass compliance (41), and where¯ Q ( ρ ) = b b ¯ Q ( ρ ) + ρ α ¯ Q ( ρ ) , with (51)¯ Q ( ρ ) = 1¯ a sin ( απ ) , ¯ Q ( ρ ) = 1¯ a sin (2 απ ) + (cid:18) ¯ a ¯ a − b b (cid:19) ρ α sin ( απ ) . The the creep compliances in integral form (46), (48), and (50), respectively corresponding to Models I -V, Models VI and VII, and Model VIII will be calculated in Section B.1 by the definition of inverse Laplacetransform using the creep compliance in complex domain (35), (36), and (37).In order to prove that the thermodynamical restrictions are less restrictive than the conditions for creepcompliance to be the Bernstein function, the creep compliances in integral form (46), (48), and (50) will beused. Namely, in the case of Models I - V, by requiring the non-negativity of kernel K (47), appearing in (46),one has ε cr ( t ) ≥ − k d k d t k ˙ ε cr ( t ) ≥ , k ∈ N , t > , (52)i.e., the creep compliance (46) is a Bernstein function. The conditions for non-negativity of the kernel K will be derived in Sections 4.1 - 4.5 and it will be proved that these conditions are more restrictive than thecorresponding thermodynamical restrictions. In the case of Models VI and VII from (48) and for Model VIIIfrom (50) one respectively has ˙ ε cr ( t ) = a b f cr ( t ) and ˙ ε cr ( t ) = ¯ a b ¯ f cr ( t ) . If functions f cr (48) and ¯ f cr (50) are completely monotonic, then (52) holds, i.e., the creep compliances (48)and (50) are Bernstein functions. The conditions for completely monotonicity of functions f cr and ¯ f cr will bederived in Sections 4.6 - 4.8 by requiring that the kernels Q (49) and ¯ Q (51) are non-negative and it will beproved that these conditions are narrower than the thermodynamical restrictions. Assuming ε = H and by applying the Laplace transform to the unified constitutive equation (33) in the case ofModels I - V, as well as to the unified constitutive equation (34) in the case of Models VI and VII, and ModelVII, the relaxation moduli in complex domain take the respective forms˜ σ sr ( s ) = 1 s b s µ + b s µ + η a s α + a s β + a s γ = 1 s − µ b + b s η Ψ ( s ) , (53)˜ σ sr ( s ) = 1 s b s β + b s β + η a s α + a s β + a s β + η = b a s − ψ ( s ) ψ ( s ) + s β (cid:16) b b + s η (cid:17) , (54)˜ σ sr ( s ) = 1 s b s α + b s α a s α + ¯ a s α = b ¯ a s − ¯ ψ ( s )¯ ψ ( s ) + s α (cid:16) b b + s α (cid:17) , (55)with functions Ψ , ψ, and ¯ ψ given by (38), (39), (40), respectively. Note that functions in the denominator ofrelaxation moduli in complex domain (54) and (55) are (up to multiplication with constant a or ¯ a ) eitherfunction Ψ , with γ = β + η, given by (38), or the following function of complex variableΦ ( s ) = 1 + ¯ a s α + ¯ a s α . (56)The relaxation modulus at initial time-instant either tends to infinity, or has a jump, while it tends to zerofor large time, since by the Tauberian theorem glass and equilibrium moduli are σ ( g ) sr = σ sr (0) = lim s →∞ ( s ˜ σ sr ( s )) = ∞ for Models I - V, b a , for Models VI and VII, b ¯ a , for Model VIII, (57) σ ( e ) sr = lim t →∞ σ sr ( t ) = lim s → ( s ˜ σ sr ( s )) = 0 , (58)where ˜ σ sr is the relaxation modulus in complex domain, given by (53), or (54), or (55).9n the case of Models I - V, by inverting the Laplace transform in (53), the relaxation modulus is obtainedin the integral form as σ sr ( t ) = 1 π (cid:90) ∞ K ( ρ ) | Ψ ( ρ e i π ) | e − ρt ρ − µ d ρ + g sr ( t ) , with (59) g ( t ) = , in Case 1, f ∗ sr ( ρ ∗ ) e − ρ ∗ t , in Case 2, f ( r ) sr ( t ) , in Case 3, (60)where functions K and Ψ are given by (47) and (38), ρ ∗ is determined from the equation a sin ( απ ) a | sin ( γπ ) | + a sin ( βπ ) a | sin ( γπ ) | ( ρ ∗ ) β − α = ( ρ ∗ ) γ − α , (61)while functions f ∗ sr and f ( r ) sr are given by f ∗ sr ( ρ ∗ ) = Re (cid:32) b + b ( ρ ∗ ) η e i ηπ αa ( ρ ∗ ) α e i απ + βa ( ρ ∗ ) β e i βπ + γa ( ρ ∗ ) γ e i γπ e i µπ (cid:33) ( ρ ∗ ) µ , (62) f ( r ) sr ( t ) = 2e ρ t cos ϕ Re (cid:32) b + b s η αa s α − + βa s β − + γa s γ − (cid:12)(cid:12)(cid:12)(cid:12) s = ρ e i ϕ e − i(1 − µ ) ϕ ρ − µ e i ρ t sin ϕ (cid:33) , (63)with s = ρ e i ϕ , ϕ ∈ (cid:0) π , π (cid:1) , being one of the complex conjugated zeros of function Ψ , given by (38), whilethe relaxation modulus in integral form corresponding to Models VI and VII is obtained using the relaxationmodulus in complex domain (54) as σ sr ( t ) = b a − b a (cid:90) t f sr ( τ ) d τ − b a g sr ( t ) , with (64) f sr ( t ) = 1 π (cid:90) ∞ Q ( ρ ) (cid:12)(cid:12)(cid:12) ψ ( ρ e i π ) + ρ β e i βπ (cid:16) b b + ρ η e i ηπ (cid:17)(cid:12)(cid:12)(cid:12) ρ β e − ρt d ρ and (65) g sr ( t ) = , in Case 1, − f ∗ sr ( ρ ∗ ) (cid:0) − e − ρ ∗ t (cid:1) , in Case 2, (cid:82) t f ( r ) sr ( τ ) d τ , in Case 3, (66)where σ ( g ) sr = ε ( g ) cr = b a is the glass modulus (57), function Q is given by (49), ρ ∗ is determined from the equation(61) with γ = β + η, while functions f ∗ sr and f ( r ) sr are given by f ∗ sr ( ρ ∗ ) = Re a ( ρ ∗ ) α e i απ + a (cid:16) a a − b b (cid:17) ( ρ ∗ ) β e i βπ αa ( ρ ∗ ) α e i απ + βa ( ρ ∗ ) β e i βπ + ( β + η ) a ( ρ ∗ ) β + η e i( β + η ) π , (67) f ( r ) sr ( t ) = 2e ρ t cos ϕ Re a s α + a (cid:16) a a − b b (cid:17) s β αa s α − + βa s β − + ( β + η ) a s β + η − (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) s = ρ e i ϕ e i ρ t sin ϕ , (68)with s = ρ e i ϕ , ϕ ∈ (cid:0) π , π (cid:1) , being one of the complex conjugated zeros of function Ψ , given by (38). In thecase of Model VIII, the relaxation modulus is also obtained by the Laplace transform inversion of the relaxationmodulus in complex domain (55), and it takes the following form σ sr ( t ) = b ¯ a − b ¯ a (cid:90) t ¯ f sr ( τ ) d τ − b a ¯ g sr ( t ) , with (69)¯ f sr ( t ) = 1 π (cid:90) ∞ ¯ Q ( ρ ) (cid:12)(cid:12)(cid:12) ¯ ψ ( ρ e i π ) + ρ α e i απ (cid:16) b b + ρ α e i απ (cid:17)(cid:12)(cid:12)(cid:12) ρ α e − ρt d ρ and (70)¯ g sr ( t ) = , in Case 1, − ¯ f ∗ sr ( ρ ∗ ) (cid:0) − e − ρ ∗ t (cid:1) , in Case 2, (cid:82) t ¯ f ( r ) sr ( τ ) d τ , in Case 3, (71)10here σ ( g ) sr = ε ( g ) cr = b ¯ a is the glass modulus (57), function ¯ Q is given by (51), ρ ∗ is determined by ρ ∗ = απ ) (cid:115) a − (cid:18) ¯ a a (cid:19) α , while functions ¯ f ∗ sr and ¯ f ( r ) sr are given by¯ f ∗ sr ( ρ ∗ ) = 1 α ( ρ ∗ ) α Re a (cid:16) ¯ a ¯ a − b b (cid:17) ( ρ ∗ ) α e i απ ¯ a + 2¯ a ( ρ ∗ ) α e i απ e − i απ , ¯ f ( r ) sr ( t ) = 2e ρ t cos ϕ Re a (cid:16) ¯ a ¯ a − b b (cid:17) s α αs α − (¯ a + 2¯ a s α ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) s = ρ e i ϕ e i ρ t sin ϕ , (72)with s = ρ e i ϕ , ϕ ∈ (cid:0) π , π (cid:1) , being one of the complex conjugated zeros of function Φ , given by (56). Thecalculation of the relaxation moduli (59), (64), and (69) by the definition of the inverse Laplace transform andintegration in the complex plane will be given in Section B.2.In Cases 1 and 2, the relaxation moduli (59), (64), and (69) may have a non-monotonic behavior, due topossibly non-monotonic behavior of functions f sr (65) and ¯ f sr (70), even though functions g sr in (60), g sr in (66),and ¯ g sr in (71) are monotonic. However, for Models I - V the relaxation modulus (59) in Case 1 is a completelymonotonic function in the range of model parameters narrower than the thermodynamical requirements, whichis the same range as for the creep compliance (46) to be a Bernstein function, since σ sr ( t ) ≥ − k d k d t k ˙ σ sr ( t ) ≤ , k ∈ N , t > , provided that kernel K, given by (47), is non-negative. Also, for each of Models VI - VIII the relaxation moduli(64) and (69) in Case 1 are completely monotonic functions in the same, more restrictive domain of modelparameters when the creep compliance is the Bernstein function. Namely, in the case of Models VI and VIIfrom (64) and in the case of Model VIII from (69), one respectively has˙ σ sr ( t ) = − b a f sr ( t ) and ˙ σ sr ( t ) = − b ¯ a ¯ f sr ( t ) . If functions f sr (65) and ¯ f cr (70) are completely monotonic, then ( − k d k d t k ˙ σ sr ( t ) ≤ , k ∈ N , t > , and also σ sr ( t ) ≥ , t > , since σ sr monotonically decreases from σ ( g ) sr = b a , in the case of Models VI and VII, or from σ ( g ) sr = b ¯ a in the case of Model VIII, to σ ( e ) sr = 0 , see (57) and (58), implying that the relaxation moduli (64)and (69) in Case 1 are completely monotonic functions. The conditions for complete monotonicity of functions f sr and ¯ f sr are the same as for the functions f cr and ¯ f cr , since they depend on the same kernels, so that inthe same range of model parameters, narrower than the thermodynamical restrictions, the relaxation modulusis a completely monotonic function and the creep compliance is a Bernstein function. In Case 3, the relaxationmoduli (59), (64), and (69) have damped oscillatory behavior, since functions g sr in (60), g sr in (66), and ¯ g sr in (71) are oscillatory with decreasing amplitude, due to functions f ( r ) sr (63), f ( r ) sr (68), and ¯ f ( r ) sr (72).Cases in functions g sr (60), g sr (66), and ¯ g sr (71) that appear in the relaxation moduli (59), (64), and (69)are determined according to the number and position of zeros of function Ψ (38) in cases of Models I - VII andof function Φ (56) in the case of Model VIII. Namely, in Case 1 function Ψ (Φ) has no zeros in the complexplane, while in Case 2 it has a negative real zero − ρ ∗ , and in Case 3 function Ψ (Φ) has a pair of complexconjugated zeros s and ¯ s having negative real part. In the case of Model I, since γ ∈ [0 , , function Ψ has nozeros in the complex plane, implying only non-oscillatory behavior of the relaxation modulus in time domain,i.e., there exists only Case 1 in (60). In cases of Models II - V, since γ ∈ [1 , , and in cases of Models VI andVII, since γ = β + η ∈ [1 , , the number and position of zeros of function Ψ is as followsCase 1: if Re Ψ ( ρ ∗ ) < , then Ψ has no zeros in the complex plane;Case 2: if Re Ψ ( ρ ∗ ) = 0 , then Ψ has one negative real zero − ρ ∗ ;Case 3: if Re Ψ ( ρ ∗ ) > , then Ψ has a pair of complex conjugatedzeros s and ¯ s having negative real part;where Re Ψ ( ρ ∗ ) = 1 + a ( ρ ∗ ) α cos ( απ ) + a ( ρ ∗ ) β cos ( βπ ) + a ( ρ ∗ ) γ cos ( γπ ) , ρ ∗ determined from (61), allowing for both non-oscillatory and oscillatory behavior of the relaxationmodulus. The analysis of the number and position of zeros of function Ψ (38) using the argument principle isperformed in Section A.1. In the case of Model VIII it will be shown in Section A.2 that the following holdstrue for function Φ , given by (56):Case 1: if (cid:16) ¯ a a (cid:17) ≥ a , orif (cid:16) ¯ a a (cid:17) < a and a < b | cos( απ ) | sin( απ ) , then Φ has no zeros in the complex plane;Case 2: if (cid:16) ¯ a a (cid:17) < a and a = b | cos( απ ) | sin( απ ) , then Φ has one negative real zero − ρ ∗ determined by ρ ∗ = (cid:16) b sin( απ ) (cid:17) α ;Case 3: if (cid:16) ¯ a a (cid:17) < a and a > b | cos( απ ) | sin( απ ) , then Φ has a pair of complex conjugatedzeros s and ¯ s having negative real part;where a = ¯ a a , and b = (cid:114) a − (cid:16) ¯ a a (cid:17) . Kernels
K, Q, and ¯ Q, respectively given by (47), (49), and (51), that appear in creep compliances (46), (48),and (50) and relaxation moduli (59), (64), and (69), respectively corresponding to Models I - V, Models VIand VII, and Model VIII, will be examined. Namely, by requiring kernels’ non-negativity, the range of modelparameters in which the creep compliance is a Bernstein function, while the relaxation modulus is completelymonotonic, will be explicitly obtained. Moreover, it will be proved that such obtained range is narrower thanthe corresponding thermodynamical restriction.The asymptotic analysis will reveal that in the vicinity of initial time-instant creep compliance starts fromthe zero value of deformation and increases: proportionally to t µ − γ + κ , with κ ∈ { α, β, γ } , in the case of ModelI, see (75); proportionally to t µ − α in the case of Models II and IV, see (82) and (94); and proportionally to t µ − β in the case of Models III and V, see (88) and (100), while the creep compliance starts from non-zero valueof deformation and increases: proportionally to t α in the case of Models VI and VIII, see (105) and (117); andproportionally to t β in the case of Model VII, see (111).The relaxation modulus for small time either decreases from infinity: proportionally to t − ( µ − γ ) − κ , with κ ∈ { α, β, γ } , in the case of Model I, see (78); proportionally to t − ( µ − α ) in the case of Models II and IV, see(84) and (96); and proportionally to t − ( µ − β ) in the case of Models III and V, see (90) and (102), or decreasesfrom finite glass modulus: proportionally to t α in the case of Models VI and VIII, see (107) and (119); andproportionally to t β in the case of Model VII, see (113).The growth of creep compliance for large time is governed: by t µ in the case of all Models I - V, see (77),(83), (89), (95), and (101); by t β in the case of Models VI and VII, see (106) and (112); and by t α in the caseof Model VIII, see (118). For all fractional Burgers models, the growth of creep compliance for large time isslower than in the case of classical Burgers model when the growth in infinity is linear, see (24).The relaxation modulus for large time tends to zero: proportionally to t − µ in the case of all Models I - V,see (79), (85), (91), (97), and (103); proportionally to t − β in the case of Models VI and VII, see (108) and (114);and proportionally to t − α in the case of Model VIII, see (120).The asymptotic analysis will be performed using the property of Laplace transform that if ˜ f ( s ) ∼ ˜ g ( s ) as s → ∞ ( s → f ( t ) ∼ g ( t ) as t → t → ∞ ), where ˜ f = L [ f ] and ˜ g = L [ g ] , i.e., if the function ˜ g is asymptotic expansion of Laplace image ˜ f , then its inverse Laplace transform g is asymptotic expansion oforiginal f. Model I, given by (4) and subject to thermodynamical restrictions (5), is obtained from the unified model (33)for η = κ ∈ { α, β, γ } . The requirement that the creep compliance (46) is a Bernstein function, while the relaxation modulus (59) iscompletely monotonic narrows down the thermodynamical restriction (5) to b b ≤ a i sin ( µ − κ ) π sin ( µ + κ ) π cos ( µ − κ ) π (cid:12)(cid:12)(cid:12) cos ( µ + κ ) π (cid:12)(cid:12)(cid:12) , (73)12ith ( κ , i ) ∈ { ( α, , ( β, , ( γ, } , since sin ( µ − κ ) π sin ( µ + κ ) π ≤ . (74)The requirement (73) is obtained by insuring the non-negativity of the terms appearing in brackets infunction K, given by (47) and having for κ ∈ { α, β, γ } the respective forms K ( ρ ) = b sin ( µπ )+ b ρ α | sin (( µ + α ) π ) | (cid:16) a µ − α ) π ) | sin(( µ + α ) π ) | − b b (cid:17) + a b ρ β sin (( µ − β ) π ) + a b ρ γ sin (( µ − γ ) π )+ a b ρ α sin ( µπ ) + a b ρ α + β sin (( µ − β + α ) π ) + a b ρ α + γ sin (( µ − γ + α ) π ) ,a b ρ α sin (( µ − α ) π ) + b ρ β | sin (( µ + β ) π ) | (cid:16) a µ − β ) π ) | sin(( µ + β ) π ) | − b b (cid:17) + a b ρ γ sin (( µ − γ ) π )+ a b ρ α + β sin (( µ − α + β ) π ) + a b ρ β sin ( µπ ) + a b ρ β + γ sin (( µ − γ + β ) π ) ,a b ρ α sin (( µ − α ) π ) + a b ρ β sin (( µ − β ) π ) + b ρ γ | sin (( µ + γ ) π ) | (cid:16) a µ − γ ) π ) | sin(( µ + γ ) π ) | − b b (cid:17) + a b ρ α + γ sin (( µ − α + γ ) π ) + a b ρ β + γ sin (( µ − β + γ ) π ) + a b ρ γ sin ( µπ ) , since all other terms in K are non-negative.Rewriting the left-hand-side of (74), one obtainssin ( µ − κ ) π sin ( µ + κ ) π = tan µπ − tan κ π tan µπ + tan κ π ≤ , with κ ∈ { α, β, γ } . According to thermodynamical restriction (5), one has 0 ≤ κ π ≤ µπ ≤ π implying0 ≤ tan κ π ≤ tan µπ , that yields the inequality (74). The asymptotic expansion of creep compliance (46) for Model I near initial time-instant is obtained in the form ε cr ( t ) = a b t µ − γ + α Γ(1+ µ − γ + α ) + O ( t µ + α − ρ ) , for κ = α, a b t µ − γ + β Γ(1+ µ − γ + β ) + O (cid:0) t µ + β − ρ (cid:1) , for κ = β, a b t µ Γ(1+ µ ) + a b t µ + γ − β Γ(1+ µ + γ − β ) + a b t µ + γ − α Γ(1+ µ + γ − α ) − b b (cid:16) a − b b (cid:17) t µ + γ Γ(1+ µ + γ ) + O (cid:0) t µ +2 γ − β (cid:1) , for κ = γ, as t → , (75)with ρ = max { β, γ − α } and ρ = max { β, γ − β } , as the inverse Laplace transform of the creep compliancein complex domain (35) rewritten as˜ ε cr ( s ) = b s µ + α ( a s γ + O ( s ρ )) , for κ = α, b s µ + β ( a s γ + O ( s ρ )) , for κ = β, b s µ + γ (cid:16) a s γ + a s β + a s α + 1 − a b b + O (cid:0) s − γ + β (cid:1)(cid:17) , for κ = γ, as s → ∞ , using the binomial formula 11 + x = 1 − x + x + O (cid:0) x (cid:1) , as x → , (76)while the asymptotic expansion of creep compliance (46) for large time takes the form ε cr ( t ) = b t µ Γ(1+ µ ) + b (cid:16) a − b b (cid:17) t µ − α Γ(1+ µ − α ) + O ( t µ − ρ ) , for κ = α, b t µ Γ(1+ µ ) + a b t µ − α Γ(1+ µ − α ) + b (cid:16) a − b b (cid:17) t µ − β Γ(1+ µ − β ) + O ( t µ − ρ ) , for κ = β, b t µ Γ(1+ µ ) + a b t µ − α Γ(1+ µ − α ) + a b t µ − β Γ(1+ µ − β ) + b (cid:16) a − b b (cid:17) t µ − γ Γ(1+ µ − γ ) + O ( t µ − γ − α ) , for κ = γ, as t → ∞ , (77)with ρ = min { α, β } and ρ = min { α + β, γ } , and it is obtained as the inverse Laplace transform of (46):˜ ε cr ( s ) = b s µ (cid:16) a b − b b s α + O ( s ρ ) (cid:17) , for κ = α, b s µ (cid:16) a s α + a b − b b s β + O ( s ρ ) (cid:17) , for κ = β, b s µ (cid:16) a s α + a s β + a b − b b s γ + O ( s γ + α ) (cid:17) , for κ = γ, as s → , σ sr ( t ) = b a t − ( µ − γ ) − α Γ(1 − ( µ − γ ) − α ) + O ( t − µ + γ + ρ ) , for κ = α, b a t − ( µ − γ ) − β Γ(1 − ( µ − γ ) − β ) + O ( t − µ + γ + ρ ) , for κ = β, b a t − µ Γ(1 − µ ) − a b a t − ( µ − γ ) − β Γ(1 − ( µ − γ ) − β ) + O ( t − µ + γ − ρ ) , for κ = γ, as t → , (78)with ρ = min { , γ − β − α } , ρ = min { , γ − β } , and ρ = max { α, γ − β } , as the inverse Laplace trans-form of the relaxation modulus in complex domain (53) rewritten as˜ σ sr ( s ) = a s − µ + γ (cid:0) b s α + O (cid:0) s ρ (cid:1)(cid:1) , for κ = α, a s − µ + γ (cid:0) b s β + O (cid:0) s ρ (cid:1)(cid:1) , for κ = β, a s − µ + γ (cid:16) b s γ − a b a s β + O ( s ρ ) (cid:17) , for κ = γ, as s → ∞ , using the binomial formula (76), while the asymptotic expansion of the relaxation modulus (59) for large timetakes the form σ sr ( t ) = b t − µ Γ(1 − µ ) − b (cid:16) a − b b (cid:17) t − µ − α Γ(1 − µ − α ) + O ( t − µ − ρ ) , for κ = α,b t − µ Γ(1 − µ ) − a b t − µ − α Γ(1 − µ − α ) + O ( t − µ − ρ ) , for κ = { β, γ } , as t → ∞ , (79)with ρ = min { α, β } , and it is obtained as the inverse Laplace transform of (53):˜ σ sr ( s ) = (cid:40) s − µ ( b + ( b − a b ) s α + O ( s ρ )) , for κ = α, s − µ ( b − a b s α + O ( s ρ )) , for κ = { β, γ } , as s → , where the binomial formula (76) is used. Model II, given by (6) and subject to thermodynamical restrictions (7), is obtained from the unified model (33)for η = α and γ = 2 α. The requirement that the creep compliance (46) is a Bernstein function, while the relaxation modulus (59) iscompletely monotonic narrows down the thermodynamical restriction (7) to a a (cid:12)(cid:12)(cid:12) sin ( µ − α ) π (cid:12)(cid:12)(cid:12) sin µπ cos ( µ − α ) π cos µπ ≤ b b ≤ a sin ( µ − α ) π sin ( µ + α ) π cos ( µ − α ) π (cid:12)(cid:12)(cid:12) cos ( µ + α ) π (cid:12)(cid:12)(cid:12) , (80)since cos ( µ − α ) π cos µπ ≥ ( µ − α ) π sin ( µ + α ) π ≤ . (81)The requirement (80) is obtained by insuring the non-negativity of the terms appearing in brackets infunction K, given by (47) and having the form K ( ρ ) = b sin ( µπ ) + b ρ α | sin (( µ + α ) π ) | (cid:18) a sin (( µ − α ) π ) | sin (( µ + α ) π ) | − b b (cid:19) + a b ρ β sin (( µ − β ) π ) + a b ρ α sin ( µπ ) (cid:18) b b − a a | sin (( µ − α ) π ) | sin ( µπ ) (cid:19) + a b ρ α + β sin (( µ − β + α ) π ) + a b ρ α sin (( µ − α ) π ) , since all other terms in K are non-negative.According to thermodynamical restriction (7), one has 0 ≤ α − µ ≤ α ≤ µ ≤ , so that 0 ≤ (2 α − µ ) π ≤ µπ ≤ π implies the first inequality in (81), while the second inequality holds sincesin ( µ − α ) π sin ( µ + α ) π = tan µπ − tan απ tan µπ + tan απ ≤ , due to 0 ≤ απ ≤ µπ ≤ π implying 0 ≤ tan απ ≤ tan µπ . .2.2 Asymptotic expansions The asymptotic expansion of creep compliance (46) for Model II near initial time-instant is obtained in the form ε cr ( t ) = a b t µ − α Γ (1 + µ − α ) + a b t µ − ( β − α ) Γ (1 + µ − ( β − α )) + a b b (cid:18) b b − a a (cid:19) t µ Γ (1 + µ ) − a b b t µ +2 α − β Γ (1 + µ + 2 α − β ) + O (cid:0) t µ + α (cid:1) , as t → , (82)as the inverse Laplace transform of (35) rewritten as˜ ε cr ( s ) = 1 b s µ + α a s α + a s β + a s α b b s α ˜ ε cr ( s ) = 1 b s µ + α (cid:18) a s α + a s β + a b − a b b s α − a b b s β − α + O (1) (cid:19) , as s → ∞ , using the binomial formula (76), while the asymptotic expansion of creep compliance (46) for large time takesthe form ε cr ( t ) = 1 b t µ Γ (1 + µ ) + 1 b (cid:18) a − b b (cid:19) t µ − α Γ (1 + µ − α ) + a b t µ − β Γ (1 + µ − β ) + O (cid:0) t µ − α (cid:1) , as t → ∞ , (83)and it is obtained as the inverse Laplace transform of (46):˜ ε cr ( s ) = 1 b s µ a s α + a s β + a s α b b s α ˜ ε cr ( s ) = 1 b s µ (cid:18) a b − b b s α + a s β + O (cid:0) s α (cid:1)(cid:19) , as s → , where the binomial formula (76) is used.The asymptotic expansion of relaxation modulus (59) for Model II near initial time-instant is obtained inthe form σ sr ( t ) = b a t − ( µ − α ) Γ (1 − ( µ − α )) + O (cid:0) t α − µ − ρ (cid:1) , as t → , (84)with ρ = max { β − α, α − β } , as the inverse Laplace transform of the relaxation modulus in complex domain(53) rewritten as ˜ σ sr ( s ) = 1 a s α − µ b + b s α a a s α + a a s α − β + a s α ˜ σ sr ( s ) = 1 a s α − µ ( b s α + O ( s ρ )) , as s → ∞ , using the binomial formula (76), while the asymptotic expansion of the relaxation modulus (59) for large timetakes the form σ sr ( t ) = b t − µ Γ (1 − µ ) − b (cid:18) a − b b (cid:19) t − µ − α Γ (1 − µ − α ) − a b t − µ − β Γ (1 − µ − β ) + O (cid:0) t − µ − α (cid:1) , as t → ∞ , (85)and it is obtained as the inverse Laplace transform of (53):˜ σ sr ( s ) = 1 s − µ b + b s α a s α + a s β + a s α ˜ σ sr ( s ) = 1 s − µ (cid:0) b + ( b − a b ) s α − a b s β + O (cid:0) s α (cid:1)(cid:1) , as s → , where the binomial formula (76) is used. Model III, given by (8) and subject to thermodynamical restrictions (9), is obtained from the unified model(33) for η = α and γ = α + β. .3.1 Restrictions on range of model parameters The requirement that the creep compliance (46) is a Bernstein function, while the relaxation modulus (59) iscompletely monotonic narrows down the thermodynamical restriction (9) to a a (cid:12)(cid:12)(cid:12) sin ( µ − β − α ) π (cid:12)(cid:12)(cid:12) sin ( µ − β + α ) π cos ( µ − β − α ) π cos ( µ − β + α ) π ≤ b b ≤ a sin ( µ − α ) π sin ( µ + α ) π cos ( µ − α ) π (cid:12)(cid:12)(cid:12) cos ( µ + α ) π (cid:12)(cid:12)(cid:12) , (86)since cos ( µ − β − α ) π cos ( µ − β + α ) π ≥ ( µ − α ) π sin ( µ + α ) π ≤ . (87)The requirement (86) is obtained by insuring the non-negativity of the terms appearing in brackets infunction K, given by (47) and having the form K ( ρ ) = b sin ( µπ ) + b ρ α | sin (( µ + α ) π ) | (cid:18) a sin (( µ − α ) π ) | sin (( µ + α ) π ) | − b b (cid:19) + a b ρ β sin (( µ − β ) π )+ a b ρ α sin ( µπ ) + a b ρ α + β sin (( µ − β + α ) π ) (cid:18) b b − a a | sin (( µ − β − α ) π ) | sin (( µ − β + α ) π ) (cid:19) + a b ρ α + β sin (( µ − β ) π ) , since all other terms in K are non-negative.Rewriting the left-hand-sides of (87), one obtainscos ( µ − β − α ) π cos ( µ − β + α ) π = 1 + tan ( µ − β ) π tan απ − tan ( µ − β ) π tan απ ≥ ( µ − α ) π sin ( µ + α ) π = tan µπ − tan απ tan µπ + tan απ ≤ . According to thermodynamical restriction (9), one has 0 ≤ µ − β + α ≤ , so that 0 ≤ ( µ − β + α ) π ≤ π implies cos ( µ − β + α ) π ≥ , thus yielding the first inequality in (87), while the second inequality holds since0 ≤ απ ≤ µπ ≤ π implies 0 ≤ tan απ ≤ tan µπ . The asymptotic expansion of creep compliance (46) for Model III near initial time-instant is obtained in theform ε cr ( t ) = a b t µ − β Γ (1 + µ − β ) + b a b (cid:18) b b − a a (cid:19) t µ − ( β − α ) Γ (1 + µ − ( β − α )) + O (cid:0) t µ + α − ρ (cid:1) , as t → , (88)with ρ = max { α, β − α } , as the inverse Laplace transform of (35) rewritten as˜ ε cr ( s ) = 1 b s µ + α a s α + a s β + a s α + β b b s α ˜ ε cr ( s ) = 1 b s µ + α (cid:18) a s α + β + a b − a b b s β + O ( s ρ ) (cid:19) , as s → ∞ , using the binomial formula (76), while the asymptotic expansion of creep compliance (46) for large time takesthe form ε cr ( t ) = 1 b t µ Γ (1 + µ ) + 1 b (cid:18) a − b b (cid:19) t µ − α Γ (1 + µ − α ) + O (cid:0) t µ − ρ (cid:1) , as t → ∞ , (89)with ρ = min { α, β } , and it is obtained as the inverse Laplace transform of (46):˜ ε cr ( s ) = 1 b s µ a s α + a s β + a s α + β b b s α ˜ ε cr ( s ) = 1 b s µ (cid:18) a b − b b s α + O ( s ρ ) (cid:19) , as s → , where the binomial formula (76) is used.The asymptotic expansion of relaxation modulus (59) for Model III near initial time-instant is obtained inthe form σ sr ( t ) = b a t − ( µ − β ) Γ (1 − ( µ − β )) − a b a (cid:18) b b − a a (cid:19) t α + β − µ Γ (1 + α + β − µ ) + O (cid:0) t α + β + ρ − µ (cid:1) , as t → , (90)16ith ρ = min { α, β − α } , as the inverse Laplace transform of the relaxation modulus in complex domain (53)rewritten as ˜ σ sr ( s ) = 1 a s α + β − µ b + b s α a a s β + a a s α + a s α + β ˜ σ sr ( s ) = 1 a s α + β − µ (cid:18) b s α + a b − a b a + O (cid:18) s ρ (cid:19)(cid:19) , as s → ∞ , using the binomial formula (76), while the asymptotic expansion of the relaxation modulus (59) for large timetakes the form σ sr ( t ) = b t − µ Γ (1 − µ ) − b (cid:18) a − b b (cid:19) t − µ − α Γ (1 − µ − α ) + O (cid:0) t − µ − ρ (cid:1) , as t → ∞ , (91)with ρ = min { α, β } , and it is obtained as the inverse Laplace transform of (53):˜ σ sr ( s ) = 1 s − µ b + b s α a s α + a s β + a s α + β ˜ σ sr ( s ) = 1 s − µ ( b + ( b − a b ) s α + O ( s ρ )) , as s → , where the binomial formula (76) is used. Model IV, given by (10) and subject to thermodynamical restrictions (11), is obtained from the unified model(33) for η = β and γ = α + β. The requirement that the creep compliance (46) is a Bernstein function, while the relaxation modulus (59) iscompletely monotonic narrows down the thermodynamical restriction (11) to a a (cid:12)(cid:12)(cid:12) sin ( µ − α − β ) π (cid:12)(cid:12)(cid:12) sin ( µ − α + β ) π cos ( µ − α − β ) π cos ( µ − α + β ) π ≤ b b ≤ a sin ( µ − β ) π sin ( µ + β ) π cos ( µ − β ) π (cid:12)(cid:12)(cid:12) cos ( µ + β ) π (cid:12)(cid:12)(cid:12) , (92)since cos ( µ − α − β ) π cos ( µ − α + β ) π ≥ ( µ − β ) π sin ( µ + β ) π ≤ . (93)The requirement (92) is obtained by insuring the non-negativity of the terms appearing in brackets infunction K, given by (47) and having the form K ( ρ ) = b sin ( µπ ) + a b ρ α sin (( µ − α ) π ) + b ρ β | sin (( µ + β ) π ) | (cid:18) a sin (( µ − β ) π ) | sin (( µ + β ) π ) | − b b (cid:19) + a b ρ α + β sin (( µ + β − α ) π ) (cid:18) b b − a a | sin (( µ − α − β ) π ) | sin (( µ − α + β ) π ) (cid:19) + a b ρ β sin ( µπ ) + a b ρ α +2 β sin (( µ − α ) π ) , since all other terms in K are non-negative.Rewriting the left-hand-sides of (93), one obtainscos ( µ − α − β ) π cos ( µ − α + β ) π = 1 + tan ( µ − α ) π tan βπ − tan ( µ − α ) π tan βπ ≥ ( µ − β ) π sin ( µ + β ) π = tan µπ − tan βπ tan µπ + tan βπ ≤ . According to thermodynamical restriction (11), one has 0 ≤ µ − α + β ≤ , so that 0 ≤ ( µ − α + β ) π ≤ π implies cos ( µ − α + β ) π ≥ , thus yielding the first inequality in (93), while the second inequality holds since0 ≤ βπ ≤ µπ ≤ π implies 0 ≤ tan βπ ≤ tan µπ . .4.2 Asymptotic expansions The asymptotic expansion of creep compliance (46) for Model IV near initial time-instant is obtained in theform ε cr ( t ) = a b t µ − α Γ (1 + µ − α ) + a b t µ Γ (1 + µ ) + a b b (cid:18) b b − a a (cid:19) t µ + β − α Γ (1 + µ + β − α ) − b b (cid:18) a − b b (cid:19) t µ + β Γ (1 + µ + β ) + O (cid:0) t µ +2 β − α (cid:1) , as t → , (94)as the inverse Laplace transform of (35) rewritten as˜ ε cr ( s ) = 1 b s µ + β a s α + a s β + a s α + β b b s β ˜ ε cr ( s ) = 1 b s µ + β (cid:18) a s α + β + a s β + a b − a b b s α + 1 − a b b + O (cid:0) s − β + α (cid:1)(cid:19) , as s → ∞ , using the binomial formula (76), while the asymptotic expansion of creep compliance (46) for large time takesthe form ε cr ( t ) = 1 b t µ Γ (1 + µ ) + a b t µ − α Γ (1 + µ − α ) + 1 b (cid:18) a − b b (cid:19) t µ − β Γ (1 + µ − β ) − a b (cid:18) b b − a a (cid:19) t − ( α + β − µ ) Γ (1 − ( α + β − µ )) + O (cid:0) t µ − β (cid:1) , as t → ∞ , (95)and it is obtained as the inverse Laplace transform of (46):˜ ε cr ( s ) = 1 b s µ a s α + a s β + a s α + β b b s β ˜ ε cr ( s ) = 1 b s µ (cid:18) a s α + a b − b b s β + a b − a b b s α + β + O (cid:0) s β (cid:1)(cid:19) , as s → , where the binomial formula (76) is used.The asymptotic expansion of relaxation modulus (59) for Model IV near initial time-instant is obtained inthe form σ sr ( t ) = b a t − ( µ − α ) Γ (1 − ( µ − α )) − a b a t α − µ Γ (1 + 2 α − µ ) − a b a (cid:18) b b − a a (cid:19) t α + β − µ Γ (1 + α + β − µ ) + O (cid:0) t α − µ (cid:1) , as t → , (96)as the inverse Laplace transform of the relaxation modulus in complex domain (53) rewritten as˜ σ sr ( s ) = 1 a s α + β − µ b + b s β a a s β + a a s α + a s α + β ˜ σ sr ( s ) = 1 a s α + β − µ (cid:18) b s β − a b a s β − α + a b − a b a + O (cid:18) s α − β (cid:19)(cid:19) , as s → ∞ , using the binomial formula (76), while the asymptotic expansion of the relaxation modulus (59) for large timetakes the form σ sr ( t ) = b t − µ Γ (1 − µ ) − a b t − µ − α Γ (1 − µ − α ) + O (cid:0) t − µ − ρ (cid:1) , as t → ∞ , (97)with ρ = min { α, β } , and it is obtained as the inverse Laplace transform of (53):˜ σ sr ( s ) = 1 s − µ b + b s β a s α + a s β + a s α + β ˜ σ sr ( s ) = 1 s − µ ( b − a b s α + O ( s ρ )) , as s → , where the binomial formula (76) is used. Model V, given by (12) and subject to thermodynamical restrictions (13), is obtained from the unified model(33) for η = β and γ = 2 β. .5.1 Restrictions on range of model parameters The requirement that the creep compliance (46) is a Bernstein function, while the relaxation modulus (59) iscompletely monotonic narrows down the thermodynamical restriction (13) to a a (cid:12)(cid:12)(cid:12) sin ( µ − β ) π (cid:12)(cid:12)(cid:12) sin µπ cos ( µ − β ) π cos µπ ≤ b b ≤ a sin ( µ − β ) π sin ( µ + β ) π cos ( µ − β ) π (cid:12)(cid:12)(cid:12) cos ( µ + β ) π (cid:12)(cid:12)(cid:12) , (98)since cos ( µ − β ) π cos µπ ≥ ( µ − β ) π sin ( µ + β ) π ≤ . (99)The requirement (98) is obtained by insuring the non-negativity of the terms appearing in brackets infunction K, given by (47) and having the form K ( ρ ) = b sin ( µπ ) + a b ρ α sin (( µ − α ) π )+ b ρ β | sin (( µ + β ) π ) | (cid:18) a sin (( µ − β ) π ) | sin (( µ + β ) π ) | − b b (cid:19) + a b ρ α + β sin (( µ + β − α ) π )+ a b ρ β sin ( µπ ) (cid:18) b b − a a | sin (( µ − β ) π ) | sin ( µπ ) (cid:19) + a b ρ β sin (( µ − β ) π ) , since all other terms in K are non-negative.According to thermodynamical restriction (13), one has 0 ≤ β − µ ≤ β ≤ µ ≤ , so that 0 ≤ (2 β − µ ) π ≤ µπ ≤ π implies the first inequality in (99), while the second inequality holds sincesin ( µ − β ) π sin ( µ + β ) π = tan µπ − tan βπ tan µπ + tan βπ ≤ , due to 0 ≤ βπ ≤ µπ ≤ π implying 0 ≤ tan βπ ≤ tan µπ . The asymptotic expansion of creep compliance (46) for Model V near initial time-instant is obtained in the form ε cr ( t ) = a b t µ − β Γ (1 + µ − β ) + a b b (cid:18) b b − a a (cid:19) t µ Γ (1 + µ ) + a b t µ + β − α Γ (1 + µ + β − α ) + O (cid:0) t µ + β (cid:1) , as t → , (100)as the inverse Laplace transform of (35) rewritten as˜ ε cr ( s ) = 1 b s µ + β a s α + a s β + a s β b b s β ˜ ε cr ( s ) = 1 b s µ + β (cid:18) a s β + a b − a b b s β + a s α + O (1) (cid:19) , as s → ∞ , using the binomial formula (76), while the asymptotic expansion of creep compliance (46) for large time takesthe form ε cr ( t ) = 1 b t µ Γ (1 + µ ) + a b t µ − α Γ (1 + µ − α ) + 1 b (cid:18) a − b b (cid:19) t µ − β Γ (1 + µ − β ) − a b b t − ( α + β − µ ) Γ (1 − ( α + β − µ )) + O (cid:0) t µ − β (cid:1) , as t → ∞ , (101)and it is obtained as the inverse Laplace transform of (46):˜ ε cr ( s ) = 1 b s µ a s α + a s β + a s β b b s β ˜ ε cr ( s ) = 1 b s µ (cid:18) a s α + a b − b b s β − a b b s α + β + O (cid:0) s β (cid:1)(cid:19) , as s → , where the binomial formula (76) is used. 19he asymptotic expansion of relaxation modulus (59) for Model V near initial time-instant is obtained inthe form σ sr ( t ) = b a t − ( µ − β ) Γ (1 − ( µ − β )) − a b a (cid:18) b b − a a (cid:19) t β − µ Γ (1 + 2 β − µ ) − a b a t β + α − µ Γ (1 + 3 β + α − µ ) + O (cid:0) t β − µ (cid:1) , as t → , (102)as the inverse Laplace transform of the relaxation modulus in complex domain (53) rewritten as˜ σ sr ( s ) = 1 a s β − µ b + b s β a a s β − α + a a s β + a s β ˜ σ sr ( s ) = 1 a s β − µ (cid:18) b s β + a b − a b a − a b a s β − α + O (cid:18) s β (cid:19)(cid:19) , as s → ∞ , using the binomial formula (76), while the asymptotic expansion of the relaxation modulus (59) for large timetakes the form σ sr ( t ) = b t − µ Γ (1 − µ ) − a b t − µ − α Γ (1 − µ − α ) + O (cid:0) t − µ − ρ (cid:1) , as t → ∞ , (103)with ρ = min { α, β } , and it is obtained as the inverse Laplace transform of (53):˜ σ sr ( s ) = 1 s − µ b + b s β a s α + a s β + a s β ˜ σ sr ( s ) = 1 s − µ ( b − a b s α + O ( s ρ )) , as s → , where the binomial formula (76) is used. Model VI, given by (14) and subject to thermodynamical restrictions (15), is obtained from the unified consti-tutive model (34) for η = α. The requirement that the creep compliance (48) is a Bernstein function, while the relaxation modulus (64) iscompletely monotonic narrows down the thermodynamical restriction (15) to a a ≤ b b ≤ a sin ( β − α ) π sin ( β + α ) π cos ( β − α ) π (cid:12)(cid:12)(cid:12) cos ( β + α ) π (cid:12)(cid:12)(cid:12) ≤ a cos ( β − α ) π (cid:12)(cid:12)(cid:12) cos ( β + α ) π (cid:12)(cid:12)(cid:12) . (104)The requirement (104) is obtained by insuring the non-negativity of function Q, given by (49) and having theform Q ( ρ ) = 1 a b b sin ( βπ ) + 1 a b b ρ α | sin (( β + α ) π ) | (cid:18) a sin (( β − α ) π ) | sin (( β + α ) π ) | − b b (cid:19) + a a ρ α sin ( βπ ) + (cid:18) a a − b b (cid:19) ρ α + β sin ( απ )in the case of Model VI. Due to thermodynamical requirement (15), one has a a − b b ≥ , while the non-negativityof Q is guaranteed if the term in brackets is non-negative, yielding b b ≤ a sin ( β − α ) π sin ( β + α ) π cos ( β − α ) π (cid:12)(cid:12)(cid:12) cos ( β + α ) π (cid:12)(cid:12)(cid:12) ≤ a cos ( β − α ) π (cid:12)(cid:12)(cid:12) cos ( β + α ) π (cid:12)(cid:12)(cid:12) , since sin ( β − α ) π sin ( β + α ) π = tan βπ − tan απ tan βπ + tan απ ≤ , due to 0 ≤ απ ≤ βπ ≤ π implying 0 ≤ tan απ ≤ tan βπ . .6.2 Asymptotic expansions The asymptotic expansion of creep compliance (48) for Model VI near initial time-instant is obtained in theform ε cr ( t ) = a b + a b b (cid:18) b b − a a (cid:19) t α Γ (1 + α ) + O (cid:0) t α + β − ρ (cid:1) , as t → , (105)with ρ = max { α, β − α } , as the inverse Laplace transform of the creep compliance in complex domain (36)rewritten as ˜ ε cr ( s ) = 1 b s α + β a s α + a s β + a s α + β b b s α ˜ ε cr ( s ) = 1 b s α + β (cid:18) a s α + β + a b − a b b s β + O ( s ρ ) (cid:19) , as s → ∞ , using the binomial formula (76), while the asymptotic expansion of creep compliance (48) for large time takesthe form ε cr ( t ) = 1 b t β Γ (1 + β ) + 1 b (cid:18) a − b b (cid:19) t β − α Γ (1 + β − α ) + O (cid:0) t β − ρ (cid:1) , as t → ∞ , (106)with ρ = min { α, β } , and it is obtained as the inverse Laplace transform of (36):˜ ε cr ( s ) = 1 b s β a s α + a s β + a s α + β b b s α ˜ ε cr ( s ) = 1 b s β (cid:18) a b − b b s α + O ( s ρ ) (cid:19) , as s → , where the binomial formula (76) is used.The asymptotic expansion of relaxation modulus (64) for Model VI near initial time-instant is obtained inthe form σ sr ( t ) = b a − a b a (cid:18) b b − a a (cid:19) t α Γ (1 + α ) + O (cid:0) t α + ρ (cid:1) , as t → , (107)with ρ = min { α, β − α } , as the inverse Laplace transform of the relaxation modulus in complex domain (54)rewritten as ˜ σ sr ( s ) = 1 a s α b + b s α a a s β + a a s α + a s α + β ˜ σ sr ( s ) = 1 a s α (cid:18) b s α + a b − a b a + O (cid:18) s ρ (cid:19)(cid:19) , as s → ∞ , using the binomial formula (76), while the asymptotic expansion of the relaxation modulus (64) for large timetakes the form σ sr ( t ) = b t − β Γ (1 − β ) − b (cid:18) a − b b (cid:19) t − β − α Γ (1 − β − α ) + O (cid:0) t − β − ρ (cid:1) , as t → ∞ , (108)with ρ = min { α, β } , and it is obtained as the inverse Laplace transform of (54):˜ σ sr ( s ) = 1 s − β b + b s α a s α + a s β + a s α + β ˜ σ sr ( s ) = 1 s − β ( b + ( b − a b ) s α + O ( s ρ )) , as s → , where the binomial formula (76) is used.Note that the term b b − a a is non-negative due to the thermodynamical requirements (15). Model VII, given by (16) and subject to thermodynamical restrictions (17), is obtained from the unified model(34) for η = β. .7.1 Restrictions on range of model parameters The requirement that the creep compliance (48) is a Bernstein function, while the relaxation modulus (64) iscompletely monotonic narrows down the thermodynamical restriction (17) to a a ≤ a ( βπ ) (cid:32) − (cid:115) − a cos ( βπ ) a (cid:33) ≤ b b ≤ a | cos ( βπ ) | , (109)provided that 2 √ a ≤ a | cos ( βπ ) | , guaranteeing the non-negativity of term under the square root. The requirement (109) is obtained by insuringthe non-negativity of function Q, given by (49) and having the form Q ( ρ ) = a a b b ρ α sin (( β − α ) π ) + a a ρ α + β sin ((2 β − α ) π )+ a a b b (cid:18)(cid:18) b b − a a (cid:19) ρ β − ρ β | cos ( βπ ) | a b b + 1 a (cid:19) sin ( βπ )in the case of Model VII. Due to thermodynamical requirement (17), the first two terms in function Q arenon-negative, as well as b b − a a ≥ , so the non-negativity of Q is guaranteed if the quadratic function in ρ β isnon-negative, i.e., if its discriminant is non-positive, yielding (cid:18) | cos ( βπ ) | a (cid:19) (cid:18) b b (cid:19) − a b b + a a ≤ , that solved with respect to b b gives a ( βπ ) (cid:32) − (cid:115) − a cos ( βπ ) a (cid:33) ≤ b b ≤ a ( βπ ) (cid:32) (cid:115) − a cos ( βπ ) a (cid:33) . (110)In order to prove the left-hand-side of (109), one uses the binomial formula (32) in (110) and obtains a ( βπ ) (cid:32) − (cid:115) − a cos ( βπ ) a (cid:33) = a a + a ( βπ ) ∞ (cid:88) k =2 (2 k − k k ! (cid:18) a cos ( βπ ) a (cid:19) k ≤ a a . On the other hand, from (110) one has b b ≤ a ( βπ ) (cid:32) (cid:115) − a cos ( βπ ) a (cid:33) ≤ a cos ( βπ )and since | cos( βπ ) | ≤ ( βπ ) , the thermodynamical requirement remains on the right-hand-side of (109). The asymptotic expansion of creep compliance (48) for Model VII near initial time-instant is obtained in theform ε cr ( t ) = a b + a b b (cid:18) b b − a a (cid:19) t β Γ (1 + β ) + a b t β − α Γ (1 + 2 β − α )+ 1 b (cid:18) − a b b + a b b (cid:19) t β Γ (1 + 2 β ) + O (cid:0) t β (cid:1) , as t → , (111)as the inverse Laplace transform of the creep compliance in complex domain (36) rewritten as˜ ε cr ( s ) = 1 b s β a s α + a s β + a s β b b s β ˜ ε cr ( s ) = 1 b s β (cid:18) a s β + (cid:18) a − a b b (cid:19) s β + a s α + 1 − a b b + a b b + O (cid:0) s − β (cid:1)(cid:19) , as s → ∞ , ε cr ( t ) = 1 b t β Γ (1 + β ) + a b t β − α Γ (1 + β − α ) + 1 b (cid:18) a − b b (cid:19) − a b b t − α Γ (1 − α ) + O (cid:0) t − β (cid:1) , as t → ∞ , (112)and it is obtained as the inverse Laplace transform of (36):˜ ε cr ( s ) = 1 b s β a s α + a s β + a s β b b s β ˜ ε cr ( s ) = 1 b s β (cid:18) a s α + (cid:18) a − b b (cid:19) s β − a b b s α + β + O (cid:0) s β (cid:1)(cid:19) , as s → , where the binomial formula (76) is used.The asymptotic expansion of relaxation modulus (64) for Model VII near initial time-instant is obtained inthe form σ sr ( t ) = b a − a b a (cid:18) b b − a a (cid:19) t β Γ (1 + β ) − a b a t β − α Γ (1 + 2 β − α ) + O (cid:0) t β (cid:1) , as t → , (113)as the inverse Laplace transform of the relaxation modulus in complex domain (54) rewritten as˜ σ sr ( s ) = 1 a s β b + b s β a a s β − α + a a s β + a s β ˜ σ sr ( s ) = 1 a s β (cid:18) b s β + a b − a b a − a b a s α + O (1) (cid:19) , as s → ∞ , using the binomial formula (76), while the asymptotic expansion of the relaxation modulus (64) for large timetakes the form σ sr ( t ) = b t − β Γ (1 − β ) − a b t − β − α Γ (1 − β − α ) − b (cid:18) a − b b (cid:19) t − β Γ (1 − β )+ a b t − β − α Γ (1 − β − α ) + O (cid:0) t − β − α (cid:1) , as t → ∞ , (114)and it is obtained as the inverse Laplace transform of (54):˜ σ sr ( s ) = 1 s − β b + b s β a s α + a s β + a s β ˜ σ sr ( s ) = 1 s − β (cid:0) b − a b s α + ( b − a b ) s β + a b s α + O (cid:0) s α + β (cid:1)(cid:1) , as s → , where the binomial formula (76) is used.Note that the term b b − a a is non-negative due to the thermodynamical requirements (17). Model VIII, given by (18) and subject to thermodynamical restrictions (19), is obtained from the unified model(34) for η = β = α, ¯ a = a + a and ¯ a = a . The requirement that the creep compliance (50) is a Bernstein function, while the relaxation modulus (69) iscompletely monotonic narrows down the thermodynamical restriction (19) to¯ a ¯ a ≤ ¯ a ( απ ) (cid:32) − (cid:115) − a cos ( απ )¯ a (cid:33) ≤ b b ≤ ¯ a | cos ( απ ) | , (115)provided that 2 √ ¯ a ≤ ¯ a | cos ( απ ) | , Q, given by (51) and having the form¯ Q ( ρ ) = ¯ a ¯ a b b (cid:18)(cid:18) b b − ¯ a ¯ a (cid:19) ρ α − b b | cos ( απ ) | ¯ a ρ α + 1¯ a (cid:19) sin ( απ )in the case of Model VIII. Due to thermodynamical requirement (19), one has b b − ¯ a ¯ a ≥ , while the non-negativity of ¯ Q is guaranteed if the quadratic function in ρ α is non-negative, i.e., if its discriminant is non-positive, yielding (cid:18) | cos ( απ ) | ¯ a (cid:19) (cid:18) b b (cid:19) − a b b + ¯ a ¯ a ≤ , that solved with respect to b b gives¯ a ( απ ) (cid:32) − (cid:115) − a cos ( απ )¯ a (cid:33) ≤ b b ≤ ¯ a ( απ ) (cid:32) (cid:115) − a cos ( απ )¯ a (cid:33) . (116)In order to prove the left-hand-side of (115), one uses the binomial formula (32) in (116) and obtains¯ a ( απ ) (cid:32) − (cid:115) − a cos ( απ )¯ a (cid:33) = ¯ a ¯ a + ¯ a ( απ ) ∞ (cid:88) k =2 (2 k − k k ! (cid:18) a cos ( απ )¯ a (cid:19) k ≤ ¯ a ¯ a . On the other hand, from (116) one has b b ≤ ¯ a ( απ ) (cid:32) (cid:115) − a cos ( απ )¯ a (cid:33) ≤ ¯ a cos ( απ )and since | cos( απ ) | ≤ ( απ ) , the thermodynamical requirement remains on the right-hand-side of (115). The asymptotic expansion of creep compliance (50) for Model VIII near initial time-instant is obtained in theform ε cr ( t ) = ¯ a b + ¯ a b b (cid:18) b b − ¯ a ¯ a (cid:19) t α Γ (1 + α ) + O (cid:0) t α (cid:1) , as t → , (117)as the inverse Laplace transform the creep compliance in complex domain (37) rewritten as˜ ε cr ( s ) = 1 b s α a s α + ¯ a s α b b s α ˜ ε cr ( s ) = 1 b s α (cid:18) ¯ a s α + ¯ a b − ¯ a b b s α + O (1) (cid:19) , as s → ∞ , using the binomial formula (76), while the asymptotic expansion of creep compliance (50) for large time takesthe form ε cr ( t ) = 1 b t α Γ (1 + α ) + ¯ a − b b + O (cid:0) t − α (cid:1) , as t → ∞ , (118)and it is obtained as the inverse Laplace transform of (37):˜ ε cr ( s ) = 1 b s α a s α + ¯ a s α b b s α ˜ ε cr ( s ) = 1 b s α (cid:18) a b − b b s α + O (cid:0) s α (cid:1)(cid:19) , as s → , where the binomial formula (76) is used.The asymptotic expansion of relaxation modulus (69) for Model VIII near initial time-instant is obtained inthe form σ sr ( t ) = b ¯ a − ¯ a b ¯ a (cid:18) b b − ¯ a ¯ a (cid:19) t α Γ (1 + α ) + O (cid:0) t α (cid:1) , as t → , (119)24s the inverse Laplace transform of the relaxation modulus in complex domain (55) rewritten as˜ σ sr ( s ) = 1¯ a s α b + b s α ¯ a ¯ a s α + a s α ˜ σ sr ( s ) = 1¯ a s α (cid:18) b s α + ¯ a b − ¯ a b ¯ a + O (cid:18) s α (cid:19)(cid:19) , as s → ∞ , using the binomial formula (76), while the asymptotic expansion of relaxation modulus (69) for large time takesthe form σ sr ( t ) = b t − α Γ (1 − α ) − b (cid:18) ¯ a − b b (cid:19) t − α Γ (1 − α ) + O (cid:0) t − α (cid:1) , as t → ∞ , (120)and it is obtained as the inverse Laplace transform of (55):˜ σ sr ( s ) = 1 s − α b + b s α a s α + ¯ a s α ˜ σ sr ( s ) = 1 s − α (cid:0) b + ( b − ¯ a b ) s α + O (cid:0) s α (cid:1)(cid:1) , as s → , where the binomial formula (76) is used.Note that the term b b − ¯ a ¯ a is non-negative due to the thermodynamical requirements (19). Models I - VIII display similar behavior in creep and stress relation tests, except for the behavior near the initialtime-instant when the creep compliance either starts at zero deformation (Models I - V), or has a jump (ModelsVI - VIII), while the relaxation modulus decreases either from infinite (Models I - V), or from a finite value ofstress (Models VI - VIII), see (41) and (57). The growth of creep compliances to infinity is governed by (43)or (46) in the case of Models I - V and by (44) or (48) in the case of Models VI - VII, i.e., by (45) or (50) forModel VIII, while the relaxation modulus either monotonically, or non-monotonically tends to zero, governedby (59) for Models I - V, or by (64) for Models VI - VII, i.e., by (69) in the case of Model VII. The materialresponses in creep and stress relaxation tests will be illustrated taking the example of Model VII.Figure 2 displays the comparison of creep compliances obtained through the representation via two-parameterMittag-Leffler function (44), presented by dots and the integral representation (64), presented by the solid line.The agreement between two different representations of the creep compliance is evidently very good. The modelparameters: α = 0 . , β = 0 . , a = 0 . , a = 0 . , a = 0 . , b = 0 . , and b = 0 . , taken in order toproduce plots from Figure 2 guarantee not only that the thermodynamical restrictions (17) are fulfilled, butalso that the requirement (109) is met, implying that the creep compliance is a Bernstein function. The sameset of model parameters is used for plots shown in Figures 3a, 3b, and 4a. t ε cr Figure 2: Creep compliance as Bernstein function expressed through: integral form - solid line, two-parameterMittag-Leffler function - dots.The plots of asymptotic expansions of the creep compliance near initial time-instant (111) and for largetime (112) are presented by the dashed lines in Figures 3a and 3b respectively, along with the creep compliance25alculated by (64) and presented by the solid line. One notices the good agreement between creep compliancecurves and its asymptotic expansions. t ε cr (a) Asymptotics for small time. t - ε cr (b) Asymptotics for large time. Figure 3: Comparison of creep compliance (solid line) and its asymptotic expansions (dashed line).Figures 4a and 4b display the relaxation modulus, calculated according to (64) in Case 1 and 2, respectively.Model parameters used to produce plots from Figure 4a are the same as for the previous plots guaranteeingthat the function Ψ (38) has no zeros in the complex plane (thus corresponding to Case 1) and moreover thatthe relaxation modulus is a completely monotonic function. Curve from Figure 4b is produced for the modelparameters: α = 0 . , β = 0 . , a = 1 , a = 1 . , a = 1 . , b = 0 . , and b = 0 . , guaranteeing that thefunction Ψ has one negative real zero. t σ sr (a) Case 1: completely monotonic function. t σ sr (b) Case 2. Figure 4: Relaxation modulus in Cases 1 and 2.Figure 5 displays the relaxation modulus, calculated for model parameters: α = 0 . , β = 0 . , a = 1 . ,a = 1 . , a = 1 . , b = 0 . , and b = 0 .
25 according to (64) in Case 3, i.e., in the case when function Ψ (38)has a pair of complex conjugated having negative real part. One notices the non-monotonic behavior of therelaxation modulus.Figures 6a and 6b display the comparison of relaxation modulus calculated according to (64) in Case 3,presented by the solid line, and its asymptotic expansions near initial time-instant (113) and for large time(114), presented by the dashed line, for the previously mentioned model parameters. Regardless of the non-monotonic behavior of relaxation modulus, the agreement between curves is good.
Thermodynamically consistent classical and fractional Burgers models I - VIII are examined in creep and stressrelaxation tests. Using the Laplace transform method, explicit forms of creep compliance in representationvia Mittag-Leffler function and in integral representation, as well as the integral representation of relaxationmodulus are obtained taking into account the material behavior at initial time-instant, since Models I - V26 t σ sr Figure 5: Relaxation modulus in Case 3: oscillatory function having decreasing amplitude. t σ sr (a) Asymptotics for small time. t - σ sr (b) Asymptotics for large time. Figure 6: Comparison of relaxation modulus (solid line) and its asymptotic expansions (dashed line).Table 1: Summary of the asymptotic analysis.Asymptotics as t →
0. Asymptotics as t → ∞ .Model creep compliance relaxation modulus creep compliance relaxation modulusI ∼ t µ − γ + κ ∼ t − ( µ − γ ) − κ ∼ t µ ∼ t − µ II ∼ t µ − α ∼ t − ( µ − α ) III ∼ t µ − β ∼ t − ( µ − β ) IV ∼ t µ − α ∼ t − ( µ − α ) V ∼ t µ − β ∼ t − ( µ − β ) VI incr. ∼ t α from a b decr. ∼ t α from b a ∼ t β ∼ t − β VII incr. ∼ t β from a b decr. ∼ t β from b a VIII incr. ∼ t α from ¯ a b decr. ∼ t α from b ¯ a ∼ t α ∼ t − α A Determination of position and number of zeros of functions Ψ and Φ A.1 Case of function Ψ Introducing the substitution s = ρ e i ϕ into function Ψ , given by (38), and by separating real and imaginaryparts in Ψ one obtainsRe Ψ ( ρ, ϕ ) = 1 + a ρ α cos ( αϕ ) + a ρ β cos ( βϕ ) + a ρ γ cos ( γϕ ) , (121)Im Ψ ( ρ, ϕ ) = a ρ α sin ( αϕ ) + a ρ β sin ( βϕ ) + a ρ γ sin ( γϕ ) , (122)where 0 < α ≤ β < , γ ∈ (0 , , and a , a , a > . Note that if s = ρ e i ϕ , ϕ ∈ (0 , π ) , is zero of function Ψ , then its complex conjugate, ¯ s = ρ e − i ϕ , is also a zero of Ψ , since Im Ψ ( ρ, − ϕ ) = − Im Ψ ( ρ, ϕ ) . Therefore, itis sufficient to consider the upper complex half-plane only.Considering the imaginary part of Ψ , (122), one concludes that function Ψ does not have zeroes in the rightcomplex half-plane. Namely, function Ψ cannot have positive real zeroes, since for ϕ = 0 , although Im Ψ = 0,it holds that Re Ψ >
0. Also, since α, β ∈ (0 ,
1) and γ ∈ (0 , , for ϕ ∈ (cid:0) , π (cid:3) all sine functions in (122) arepositive, implying that Im Ψ >
0. Moreover, in the case when γ ∈ (0 , , function Ψ does not have zeroes inthe whole complex plane, since for ϕ ∈ (0 , π ] all sine functions in (122) are positive implying that Im Ψ > . Therefore, zeros of function Ψ may exist only if γ ∈ (1 , , and, if zeros exist, they are complex conjugated andlocated in the left complex half-plane.In the case when γ ∈ (1 , ∪ Γ ∪ Γ ∪ Γ placed in the upper left complex quarter-plane, shown in Figure 7.Figure 7: Contour Γ = Γ ∪ Γ ∪ Γ ∪ Γ .Parametrizing the contour Γ by s = ρ e i π , ρ ∈ ( r, R ) , with r → R → ∞ , the real and imaginary partsof Ψ , (121) and (122), becomeRe Ψ ( ρ ) = 1 + a ρ α cos απ a ρ β cos βπ a ρ γ cos γπ , Im Ψ ( ρ ) = a ρ α sin απ a ρ β sin βπ a ρ γ sin γπ > , which in the limiting cases yieldRe Ψ ( ρ ) ∼ , Im Ψ ( ρ ) ∼ a ρ α sin απ → + , as ρ = r → , Re Ψ ( ρ ) ∼ a ρ γ cos γπ → −∞ , Im Ψ ( ρ ) ∼ a ρ γ sin γπ → ∞ , as ρ = R → ∞ . , (121) and (122), along the contour Γ , parametrized by s = R e i ϕ ,ϕ ∈ (cid:2) π , π (cid:3) , with R → ∞ , take the formRe Ψ ( ϕ ) ∼ a R γ cos( γϕ ) , Im Ψ ( ϕ ) ∼ a R γ sin( γϕ ) , as R → ∞ . (123)Note, when γ ∈ (cid:0) , (cid:1) it holds that Re Ψ < , since cos( γϕ ) < . If ϕ = π , then (123) becomesRe Ψ (cid:16) π (cid:17) → −∞ , Im Ψ (cid:16) π (cid:17) → ∞ , as R → ∞ , while if ϕ = π, then, as R → ∞ , (123) becomes eitherRe Ψ ( π ) → −∞ , Im Ψ ( π ) → −∞ , for γ ∈ (cid:18) , (cid:19) , orRe Ψ ( π ) → ∞ , Im Ψ ( π ) → −∞ , for γ ∈ (cid:20) , (cid:19) . Using the parametrization s = ρ e i π , ρ ∈ ( r, R ) , with r → R → ∞ , of contour Γ , the real andimaginary parts of Ψ , (121) and (122), becomeRe Ψ ( ρ ) = 1 + a ρ α cos ( απ ) + a ρ β cos ( βπ ) + a ρ γ cos ( γπ ) , (124)Im Ψ ( ρ ) = a ρ α sin ( απ ) + a ρ β sin ( βπ ) + a ρ γ sin ( γπ ) , (125)which, as ρ = r → , yield Re Ψ ( ρ ) ∼ , Im Ψ ( ρ ) ∼ a ρ α sin ( απ ) → + , and, as ρ = R → ∞ , eitherRe Ψ ( ρ ) ∼ a ρ γ cos ( γπ ) → −∞ , Im Ψ ( ρ ) ∼ a ρ γ sin ( γπ ) → ∞ , for γ ∈ (cid:18) , (cid:19) , orRe Ψ ( ρ ) ∼ a ρ γ cos ( γπ ) → ∞ , Im Ψ ( ρ ) ∼ a ρ γ sin ( γπ ) → ∞ , for γ ∈ (cid:20) , (cid:19) , Since α, β ∈ (0 ,
1) and γ ∈ (1 , , the first two terms in (125) are positive, while the third term is negative, sothat there exists at least one ρ ∗ (cid:54) = 0 such thatIm Ψ ( ρ ∗ ) = a ( ρ ∗ ) α | sin ( γπ ) | (cid:18) a sin ( απ ) a | sin ( γπ ) | + a sin ( βπ ) a | sin ( γπ ) | ( ρ ∗ ) β − α − ( ρ ∗ ) γ − α (cid:19) = 0 . The equation a sin ( απ ) a | sin ( γπ ) | + a sin ( βπ ) a | sin ( γπ ) | ( ρ ∗ ) β − α = ( ρ ∗ ) γ − α (126)has a single solution ρ ∗ . Namely, since 0 < α ≤ β < γ ∈ (1 , , the function ( ρ ∗ ) β − α is concave, dueto β − α ∈ (0 ,
1) (or even constant if β = α ), and ( ρ ∗ ) γ − α is either convex, if γ − α ∈ (1 , γ − α ∈ (0 ,
1) (or even linear if γ − α = 1), but always increases faster than ( ρ ∗ ) β − α . Depending on the valueof ρ ∗ , the real part of Ψ , (124), can be either negative, zero, or positive, which will determine the change ofargument of function Ψ along the contour Γ . Along the contour Γ , parametrized by s = r e i ϕ , ϕ ∈ (cid:2) π , π (cid:3) , with r → , the real and imaginary parts ofΨ , (121) and (122), take the formRe Ψ ( ϕ ) ∼ , Im Ψ ( ϕ ) ∼ a r γ sin( αϕ ) > , as r → , so that Im Ψ ( ϕ ) → + , for ϕ ∈ (cid:2) π , π (cid:3) , as r → . By the argument principle, function Ψ has no zeros (has one complex zero) in the upper left complex quarter-plane if Re Ψ ( ρ ∗ ) < ρ ∗ ) >
0) in (124), since the change of argument of function Ψ is ∆ arg Ψ = 0(∆ arg Ψ = 2 π ) along the contour Γ , while if Re Ψ ( ρ ∗ ) = 0 , then function Ψ has a negative real zero − ρ ∗ , with ρ ∗ determined from (126). This holds true for the whole complex plane as well if γ ∈ (1 , , while if γ ∈ (0 , , then function Ψ does not have zeros in the complex plane.29 .2 Case of function Φ Function Φ , given by (56), being a quadratic function in terms of s α , is decomposed asΦ ( s ) = ¯ a ( s α + a ) ( s α + a ) , with a , = ¯ a a ± (cid:115)(cid:18) ¯ a a (cid:19) − a > , for (cid:16) ¯ a a (cid:17) ≥ a , implying that function Φ has no zeros in the principal Riemann branch, i.e., for arg s ∈ ( − π, π ) , since it is well-known that equation s α + a , = 0 , s ∈ C , has no solutions for arg s ∈ ( − π, π ) . For (cid:16) ¯ a a (cid:17) < a , function Φ is decomposed asΦ ( s ) = ¯ a ( s α + a − i b ) ( s α + a + i b ) , with a = ¯ a a and b = (cid:115) a − (cid:18) ¯ a a (cid:19) . Define a function φ by φ ( s ) = s α + a − i b, s ∈ C , (127)where α ∈ (0 , , The real and imaginary parts of function φ are obtained asRe φ ( ρ, ϕ ) = ρ α cos ( αϕ ) + a, (128)Im φ ( ρ, ϕ ) = ρ α sin ( αϕ ) − b, (129)using s = ρ e i ϕ in (127). If function φ has zero s = ρ e i ϕ in the upper complex half-plane (since sin ( αϕ ) > ϕ ∈ (0 , π )), then the complex conjugated zero ¯ s = ρ e − i ϕ is obtained as a solution of equation s α + a + i b = 0 . Thus, s and ¯ s are complex conjugated zeros of function Φ as well. Function φ does not have zeroes in theupper right complex quarter-plane, since for ϕ ∈ (cid:2) , π (cid:3) its real part (128) is positive.Therefore, the zeros of φ will be sought for in the upper left complex quarter-plane within the contour Γ , shown in Figure 7, by using the argument principle. Parametrizing the contour Γ = Γ ∪ Γ ∪ Γ ∪ Γ as inSection A.1, along Γ , by (128) and (129), one obtainsRe φ ( ρ ) = ρ α cos απ a > , Im φ ( ρ ) = ρ α sin απ − b, Re φ ( ρ ) ∼ a, Im φ ( ρ ) ∼ − b, as ρ = r → , Re φ ( ρ ) → ∞ , Im φ ( ρ ) → ∞ , as ρ = R → ∞ , as real and imaginary parts of φ , and their limiting cases. Along Γ , as R → ∞ , for arbitrary value of ϕ , aswell as for ϕ = π and ϕ = π, (128) and (129) take the following formsRe φ ( ϕ ) ∼ R α cos( αϕ ) , Im φ ( ϕ ) ∼ R α sin( αϕ ) > , Re φ (cid:16) π (cid:17) → ∞ , Im φ (cid:16) π (cid:17) → ∞ , Re φ ( π ) → −∞ , Im φ ( π ) → ∞ . In the case of contour Γ , (128), (129) and their limiting cases areRe φ ( ρ ) = − ρ α | cos ( απ ) | + a, Im φ ( ρ ) = ρ α sin ( απ ) − b, (130)Re φ ( ρ ) ∼ a, Im φ ( ρ ) ∼ − b, as ρ = r → , Re φ ( ρ ) → −∞ , Im φ ( ρ ) → ∞ , as ρ = R → ∞ . Depending on the value of solution ρ ∗ = (cid:16) b sin( απ ) (cid:17) α of equation Im φ ( ρ ∗ ) = 0 , (130) , the real part of φ, (130) ,can be either negative, zero, or positive, which will determine the change of argument of function φ along thewhole contour Γ . Along the contour Γ , one hasRe φ ( ϕ ) ∼ a, Im φ ( ϕ ) ∼ − b + , as r → . By the argument principle, function φ has no zeros (has one complex zero) in the upper left complex quarter-plane and therefore in the whole upper complex half-plane as well, if Re φ ( ρ ∗ ) = − b | cos( απ ) | sin( απ ) + a < φ ( ρ ∗ ) = − b | cos( απ ) | sin( απ ) + a >
0) in (130), since the change of argument of function φ is ∆ arg φ = 0 (∆ arg φ = 2 π ) alongthe contour Γ. If Re φ ( ρ ∗ ) = 0 , then function φ has a negative real zero − ρ ∗ . Calculation of creep compliance and relaxation modulus
The creep compliance and relaxation modulus corresponding to Models I - V, respectively to Models VI andVII, will be obtained in the integral forms (46) and (59), respectively as (48) and (64), by inverting the creepcompliance and relaxation modulus in complex domain, given by (35) and (53), respectively by (36) and (54),using the definition of inverse Laplace transform f ( t ) = L − (cid:104) ˜ f ( s ) (cid:105) ( t ) = 12 π i (cid:90) p +i ∞ p − i ∞ ˜ f ( s ) e st d s. (131)In the case of creep compliance in complex domain, as well as for the relaxation modulus in complex domainwhen function Ψ , given by (38), has either no zeros in complex plane, or has one negative real zero, the Laplacetransform inversion will be performed by using the Cauchy integral theorem (cid:72) Γ f ( z )d z = 0 . In the case ofrelaxation modulus in complex domain when function Ψ (38) has a pair of complex conjugated zeros, theCauchy residue theorem (cid:72) Γ f ( z )d z = 2 π i (cid:80) k Res ( f ( z ) , z k ) will be employed.In the case of Model VIII, the creep compliance and relaxation modulus, given by (50) and (69), are alsocalculated by the Laplace transform inversion of the creep compliance and relaxation modulus in complexdomain, given by (37) and (55). Being analogous to the case of Models VI and VII the calculation is omitted. B.1 Calculation of creep compliance
B.1.1 Creep compliance in the case of models having zero glass compliance
The rate of creep is obtained as˙ ε cr ( t ) = L − [ s ˜ ε cr ( s )] ( t ) = 12 π i lim R →∞ ,r → (cid:90) Γ s ˜ ε cr ( s ) e st d s, (132)using the inverse Laplace transform (131) of the creep compliance in complex domain (35) and zero value ofthe glass compliance (41). The Cauchy integral theorem (cid:73) Γ (I) s ˜ ε cr ( s ) e st d s = 0 , (133)where the contour Γ (I) = Γ ∪ Γ ∪ Γ ∪ Γ ∪ Γ ∪ Γ ∪ Γ ∪ Γ is chosen as in Figure 8, yields the rate of creepFigure 8: Contour Γ (I) . in the form ˙ ε cr ( t ) = 12 π i (cid:90) ∞ (cid:32) Ψ (cid:0) ρ e − i π (cid:1) b + b ρ η e − i ηπ e i µπ − Ψ (cid:0) ρ e i π (cid:1) b + b ρ η e i ηπ e − i µπ (cid:33) e − ρt ρ µ d ρ (134)= 1 π (cid:90) ∞ K ( ρ ) | b + b ρ η e i ηπ | e − ρt ρ µ d ρ, (135)and thus the creep compliance becomes ε cr ( t ) = 1 π (cid:90) ∞ K ( ρ ) | b + b ρ η e i ηπ | − e − ρt ρ µ d ρ, K ( ρ ) = b K ( ρ ) + b ρ η K ( ρ ) , with K ( ρ ) = sin ( µπ ) Re Ψ (cid:0) ρ e i π (cid:1) − cos ( µπ ) Im Ψ (cid:0) ρ e i π (cid:1) = sin ( µπ ) + a ρ α sin (( µ − α ) π ) + a ρ β sin (( µ − β ) π ) + a ρ γ sin (( µ − γ ) π ) ,K ( ρ ) = sin (( µ + η ) π ) Re Ψ (cid:0) ρ e i π (cid:1) − cos (( µ + η ) π ) Im Ψ (cid:0) ρ e i π (cid:1) = sin (( µ + η ) π ) + a ρ α sin (( µ + η − α ) π ) + a ρ β sin (( µ + η − β ) π ) + a ρ γ sin (( µ + η − γ ) π ) , since Ψ (¯ s ) = − ¯Ψ ( s ) . The integrals along contours Γ (parametrized by s = ρ e i π , ρ ∈ ( r, R )) and Γ (parametrized by s = ρ e − i π ,ρ ∈ ( r, R )) readlim R →∞ ,r → (cid:90) Γ s ˜ ε cr ( s ) e st d s = (cid:90) ∞ ρ µ e i µπ Ψ (cid:0) ρ e i π (cid:1) b + b ρ η e i ηπ e − ρt e i π d ρ = (cid:90) ∞ Ψ (cid:0) ρ e i π (cid:1) b + b ρ η e i ηπ e − i µπ e − ρt ρ µ d ρ, lim R →∞ ,r → (cid:90) Γ s ˜ ε cr ( s ) e st d s = (cid:90) ∞ ρ µ e − i µπ Ψ (cid:0) ρ e − i π (cid:1) b + b ρ η e − i ηπ e − ρt e − i π d ρ = − (cid:90) ∞ Ψ (cid:0) ρ e − i π (cid:1) b + b ρ η e − i ηπ e i µπ e − ρt ρ µ d ρ, yielding the rate of creep (134) according to the Cauchy integral theorem (133), since the inverse Laplacetransform of the rate of creep in complex domain is (132) and integrals along Γ , Γ , Γ , Γ , Γ tend to zero as R → ∞ and r → . The contour Γ is parametrized by s = p + i R, p ∈ (0 , p ) , with R → ∞ , so that the integral I Γ = (cid:90) Γ s ˜ ε cr ( s ) e st d s = (cid:90) p p + i R ) µ Ψ ( p + i R ) b + b ( p + i R ) η e ( p +i R ) t d p is estimated as | I Γ | ≤ (cid:90) p | p + i R | µ | Ψ ( p + i R ) || b + b ( p + i R ) η | e pt d p. (136)Assuming s = ρ e i ϕ , since R → ∞ , one obtains ρ = (cid:112) p + R ∼ R and ϕ = arctan Rp ∼ π , so that (136)becomes lim R →∞ | I Γ | ≤ lim R →∞ (cid:90) p R µ (cid:12)(cid:12) Ψ (cid:0) R e i π (cid:1)(cid:12)(cid:12)(cid:12)(cid:12) b + b R η e i ηπ (cid:12)(cid:12) e pt d p ≤ a b lim R →∞ (cid:90) p R µ + η − γ e pt d p = 0 , since, by (38) one has (cid:12)(cid:12) Ψ (cid:0) R e i π (cid:1)(cid:12)(cid:12) ∼ a R γ , as R → ∞ , (137)as well as µ + η − γ = µ, or µ − ( γ − β ) , or µ − ( γ − α ) , for Model I, µ − α, for Models II and IV, µ − β, for Models III and V, ∈ (0 , , (138)due to η = κ ∈ { α, β, γ } , with thermodynamical restriction (5) for Model I, and due to ( η, γ ) ∈ { ( α, α ) , ( α, α + β ) , ( β, α + β ) , ( β, β ) } , with thermodynamical restrictions (7), (9), (11), and (13) for Models II - V.Analogously, it can be proved that lim R →∞ | I Γ | = 0 . The integral along contour Γ , parametrized by s = R e i ϕ , ϕ ∈ (cid:0) π , π (cid:1) , with R → ∞ , is I Γ = (cid:90) Γ s ˜ ε cr ( s ) e st d s = (cid:90) π π R µ e i µϕ Ψ (cid:0) R e i ϕ (cid:1) b + b R η e i ηϕ e Rt e i ϕ i R e i ϕ d ϕ, yielding the estimate | I Γ | ≤ (cid:90) π π R − µ (cid:12)(cid:12) Ψ (cid:0) R e i ϕ (cid:1)(cid:12)(cid:12) | b + b R η e i ηϕ | e Rt cos ϕ d ϕ. (139)Similarly as in (137), one has (cid:12)(cid:12) Ψ (cid:0) R e i ϕ (cid:1)(cid:12)(cid:12) ∼ a R γ , as R → ∞ , which, along with cos ϕ < ϕ ∈ (cid:0) π , π (cid:1) , implies that (139) in the limit when R → ∞ becomeslim R →∞ | I Γ | ≤ a b lim R →∞ (cid:90) π π R − ( µ + η − γ ) e Rt cos ϕ d ϕ = 0 , although, by (138), 1 − ( µ + η − γ ) > . By the similar arguments, lim R →∞ | I Γ | = 0 , as well.32arametrization of the contour Γ is s = r e i ϕ , ϕ ∈ ( − π, π ) , with r → , so that I Γ = (cid:90) Γ s ˜ ε cr ( s ) e st d s = (cid:90) π − π r µ e i µϕ Ψ (cid:0) r e i ϕ (cid:1) b + b r η e i ηϕ e rt e i ϕ i r e i ϕ d ϕ can be estimated by lim r → | I Γ | ≤ lim r → (cid:90) π − π r − µ (cid:12)(cid:12) Ψ (cid:0) r e i ϕ (cid:1)(cid:12)(cid:12) | b + b r η e i ηϕ | e rt cos ϕ d ϕ ≤ πb lim r → r − µ = 0 , since (cid:12)(cid:12) Ψ (cid:0) r e i ϕ (cid:1)(cid:12)(cid:12) ∼ r → . B.1.2 Creep compliance in the case of models having non-zero glass compliance
Inverting the Laplace transform in the creep compliance in complex domain (36), the creep compliance isobtained in the form ε cr ( t ) = a b + a b (cid:90) t f cr ( τ ) d τ , where ε ( g ) cr = a b is the glass compliance (41) and function f cr is defined by its Laplace transform as˜ f cr ( s ) = L [ f cr ( t )] ( s ) = 1 s β ψ ( s ) b b + s η , with function ψ defined by (39). Function f cr , having the form f cr ( t ) = 12 π i (cid:90) ∞ (cid:32) ψ (cid:0) ρ e − i π (cid:1) b b + ρ η e − i ηπ e i βπ − ψ (cid:0) ρ e i π (cid:1) b b + ρ η e i ηπ e − i βπ (cid:33) e − ρt ρ β d ρ = 1 π (cid:90) ∞ Q ( ρ ) (cid:12)(cid:12)(cid:12) b b + ρ η e i ηπ (cid:12)(cid:12)(cid:12) e − ρt ρ β d ρ, where, due to ψ (¯ s ) = − ¯ ψ ( s ) ,Q ( ρ ) = b b Q ( ρ ) + ρ η Q ( ρ ) , with Q ( ρ ) = sin ( βπ ) Re ψ (cid:0) ρ e i π (cid:1) − cos ( βπ ) Im ψ (cid:0) ρ e i π (cid:1) = 1 a sin ( βπ ) + a a ρ α sin (( β − α ) π ) ,Q ( ρ ) = sin (( β + η ) π ) Re ψ (cid:0) ρ e i π (cid:1) − cos (( β + η ) π ) Im ψ (cid:0) ρ e i π (cid:1) = 1 a sin (( β + η ) π ) + a a ρ α sin (( β + η − α ) π ) + (cid:18) a a − b b (cid:19) ρ β sin ( ηπ ) , is calculated by the inverse Laplace transform f cr ( t ) = 12 π i lim R →∞ ,r → (cid:90) Γ ˜ f cr ( s ) e st d s, see (131), using the Cauchy integral theorem (cid:73) Γ (I) ˜ f cr ( s ) e st d s = 0 , where the contour Γ (I) = Γ ∪ Γ ∪ Γ ∪ Γ ∪ Γ ∪ Γ ∪ Γ ∪ Γ is chosen as in Figure 8, since the integrals alongcontours Γ (parametrized by s = ρ e i π , ρ ∈ ( r, R )) and Γ (parametrized by s = ρ e − i π , ρ ∈ ( r, R )) readlim R →∞ ,r → (cid:90) Γ ˜ f cr ( s ) e st d s = (cid:90) ∞ ρ β e i βπ ψ (cid:0) ρ e i π (cid:1) b b + ρ η e i ηπ e − ρt e i π d ρ = (cid:90) ∞ ψ (cid:0) ρ e i π (cid:1) b b + ρ η e i ηπ e − i βπ e − ρt ρ β d ρ, lim R →∞ ,r → (cid:90) Γ ˜ f cr ( s ) e st d s = (cid:90) ∞ ρ β e − i βπ ψ (cid:0) ρ e − i π (cid:1) b b + ρ η e − i ηπ e − ρt e − i π d ρ = − (cid:90) ∞ ψ (cid:0) ρ e − i π (cid:1) b b + ρ η e − i ηπ e i βπ e − ρt ρ β d ρ, , Γ , Γ , Γ , Γ tend to zero as R → ∞ and r → . The contour Γ is parametrized by s = p + i R, p ∈ (0 , p ) , with R → ∞ , so that the integral I Γ = (cid:90) Γ ˜ f cr ( s ) e st d s = (cid:90) p p + i R ) β ψ ( p + i R ) b b + ( p + i R ) η e ( p +i R ) t d p is estimated as | I Γ | ≤ (cid:90) p | p + i R | β | ψ ( p + i R ) | (cid:12)(cid:12)(cid:12) b b + ( p + i R ) η (cid:12)(cid:12)(cid:12) e pt d p. (140)Assuming s = ρ e i ϕ , since R → ∞ , one obtains ρ = (cid:112) p + R ∼ R and ϕ = arctan Rp ∼ π , so that (140)becomes lim R →∞ | I Γ | ≤ lim R →∞ (cid:90) p R β (cid:12)(cid:12) ψ (cid:0) R e i π (cid:1)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) b b + R η e i ηπ (cid:12)(cid:12)(cid:12) e pt d p ≤ (cid:18) a a − b b (cid:19) lim R →∞ (cid:90) p R η e pt d p = 0 , since, by (39), one has (cid:12)(cid:12) ψ (cid:0) R e i π (cid:1)(cid:12)(cid:12) ∼ (cid:18) a a − b b (cid:19) R β , as R → ∞ . (141)Analogously, it can be proved that lim R →∞ | I Γ | = 0 . The integral along contour Γ , parametrized by s = R e i ϕ , ϕ ∈ (cid:0) π , π (cid:1) , with R → ∞ , is I Γ = (cid:90) Γ ˜ f cr ( s ) e st d s = (cid:90) π π R β e i βϕ ψ (cid:0) R e i ϕ (cid:1) b b + R η e i ηϕ e Rt e i ϕ i R e i ϕ d ϕ, yielding the estimate | I Γ | ≤ (cid:90) π π R − β (cid:12)(cid:12) ψ (cid:0) R e i ϕ (cid:1)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) b b + R η e i ηϕ (cid:12)(cid:12)(cid:12) e Rt cos ϕ d ϕ. (142)Similarly as in (141), (cid:12)(cid:12) ψ (cid:0) R e i ϕ (cid:1)(cid:12)(cid:12) ∼ (cid:16) a a − b b (cid:17) R β , as R → ∞ , which along with cos ϕ < ϕ ∈ (cid:0) π , π (cid:1) impliesthat (142) in the limit when R → ∞ becomeslim R →∞ | I Γ | ≤ (cid:18) a a − b b (cid:19) lim R →∞ (cid:90) π π R − η e Rt cos ϕ d ϕ = 0 . By the similar arguments, lim R →∞ | I Γ | = 0 , as well.Parametrization of the contour Γ is s = r e i ϕ , ϕ ∈ ( − π, π ) , with r → , so that I Γ = (cid:90) Γ ˜ f cr ( s ) e st d s = (cid:90) π − π r β e i βϕ ψ (cid:0) r e i ϕ (cid:1) b b + r η e i ηϕ e rt e i ϕ i r e i ϕ d ϕ can be estimated bylim r → | I Γ | ≤ lim r → (cid:90) π − π r − β (cid:12)(cid:12) ψ (cid:0) r e i ϕ (cid:1)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) b b + r η e i ηϕ (cid:12)(cid:12)(cid:12) e rt cos ϕ d ϕ ≤ π b a b lim r → r − β = 0 , since (cid:12)(cid:12) ψ (cid:0) r e i ϕ (cid:1)(cid:12)(cid:12) ∼ a as r → . B.2 Calculation of relaxation modulus
B.2.1 Case when function Ψ has no zeros in complex planeRelaxation modulus in the case of models having zero glass compliance. The Cauchy integral theoremwith the relaxation modulus in complex domain (53) as an integrand becomes (cid:73) Γ (I) ˜ σ sr ( s ) e st d s = 0 , (143)where the contour Γ (I) = Γ ∪ Γ ∪ Γ ∪ Γ ∪ Γ ∪ Γ ∪ Γ ∪ Γ is chosen as in Figure 8. The relaxation modulusis obtained in the form σ sr ( t ) = 12 π i (cid:90) ∞ (cid:18) b + b ρ η e i ηπ Ψ ( ρ e i π ) e i µπ − b + b ρ η e − i ηπ Ψ ( ρ e − i π ) e − i µπ (cid:19) e − ρt ρ − µ d ρ = 1 π (cid:90) ∞ K ( ρ ) | Ψ ( ρ e i π ) | e − ρt ρ − µ d ρ, (144)34ith function K given by (47), as a consequence of the Cauchy integral theorem (143), since the integrals alongcontours Γ (parametrized by s = ρ e i π , ρ ∈ ( r, R )) and Γ (parametrized by s = ρ e − i π , ρ ∈ ( r, R )) readlim R →∞ ,r → (cid:90) Γ ˜ σ sr ( s ) e st d s = (cid:90) ∞ ρ e i π ) − µ b + b ρ η e i ηπ Ψ ( ρ e i π ) e − ρt e i π d ρ = − (cid:90) ∞ b + b ρ η e i ηπ Ψ ( ρ e i π ) e i µπ e − ρt ρ − µ d ρ, (145)lim R →∞ ,r → (cid:90) Γ ˜ σ sr ( s ) e st d s = (cid:90) ∞ ρ e − i π ) − µ b + b ρ η e − i ηπ Ψ ( ρ e − i π ) e − ρt e − i π d ρ = (cid:90) ∞ b + b ρ η e − i ηπ Ψ ( ρ e − i π ) e − i µπ e − ρt ρ − µ d ρ, (146)respectively, with σ sr ( t ) = 12 π i lim R →∞ ,r → (cid:90) Γ ˜ σ sr ( s ) e st d s, (147)as the inverse Laplace transform, given by (131), while the integrals along Γ , Γ , Γ , Γ , Γ tend to zero as R → ∞ and r → . The contour Γ is parametrized by s = p + i R, p ∈ (0 , p ) , with R → ∞ , so that the integral I Γ = (cid:90) Γ ˜ σ sr ( s ) e st d s = (cid:90) p p + i R ) − µ b + b ( p + i R ) η Ψ ( p + i R ) e ( p +i R ) t d p is estimated as | I Γ | ≤ (cid:90) p | p + i R | − µ | b + b ( p + i R ) η || Ψ ( p + i R ) | e pt d p. (148)Assuming s = ρ e i ϕ , since R → ∞ , one obtains ρ = (cid:112) p + R ∼ R and ϕ = arctan Rp ∼ π , so that (148)becomes lim R →∞ | I Γ | ≤ lim R →∞ (cid:90) p R − µ (cid:12)(cid:12) b + b R η e i ηπ (cid:12)(cid:12)(cid:12)(cid:12) Ψ (cid:0) R e i π (cid:1)(cid:12)(cid:12) e pt d p ≤ b a lim R →∞ (cid:90) p R − µ − η + γ e pt d p = 0 , due to (137) and inequality 1 − ( µ + η − γ ) > R →∞ | I Γ | = 0 . The integral along contour Γ , parametrized by s = R e i ϕ , ϕ ∈ (cid:0) π , π (cid:1) , with R → ∞ , is I Γ = (cid:90) Γ ˜ σ sr ( s ) e st d s = (cid:90) π π R − µ e i(1 − µ ) ϕ b + b R η e i ηϕ Ψ ( R e i ϕ ) e Rt e i ϕ i R e i ϕ d ϕ, yielding the estimate | I Γ | ≤ (cid:90) π π R µ (cid:12)(cid:12) b + b R η e i ηϕ (cid:12)(cid:12) | Ψ ( R e i ϕ ) | e Rt cos ϕ d ϕ. (149)Similarly as in (137), (cid:12)(cid:12) Ψ (cid:0) R e i ϕ (cid:1)(cid:12)(cid:12) ∼ a R γ , as R → ∞ , which along with cos ϕ < ϕ ∈ (cid:0) π , π (cid:1) implies that(149) in the limit when R → ∞ becomeslim R →∞ | I Γ | ≤ b a lim R →∞ (cid:90) π π R µ + η − γ e Rt cos ϕ d ϕ = 0 , although, by (138) µ + η − γ > . By the similar arguments, lim R →∞ | I Γ | = 0 , as well.Parametrization of the contour Γ is s = r e i ϕ , ϕ ∈ ( − π, π ) , with r → , so that I Γ = (cid:90) Γ ˜ σ sr ( s ) e st d s = (cid:90) π − π r − µ e i(1 − µ ) ϕ b + b r η e i ηϕ Ψ ( r e i ϕ ) e rt e i ϕ i r e i ϕ d ϕ, can be estimated by lim r → | I Γ | ≤ lim r → (cid:90) π − π r µ (cid:12)(cid:12) b + b r η e i ηϕ (cid:12)(cid:12) | Ψ ( r e i ϕ ) | e rt cos ϕ d ϕ ≤ πb lim r → r µ = 0 , since (cid:12)(cid:12) Ψ (cid:0) r e i ϕ (cid:1)(cid:12)(cid:12) ∼ r → . Relaxation modulus in the case of models having non-zero glass compliance.
Inverting theLaplace transform in the relaxation modulus in complex domain (54), the relaxation modulus is obtained in theform σ sr ( t ) = b a − b a (cid:90) t L − (cid:104) ˜ f sr ( s ) (cid:105) ( τ ) d τ = b a − b a (cid:90) t f sr ( τ ) d τ , σ ( g ) sr = a b is the glass modulus (57), and function f sr is defined by its Laplace transform as˜ f sr ( s ) = ψ ( s ) ψ ( s ) + s β (cid:16) b b + s η (cid:17) , (150)with function ψ defined by (39). Function f sr , having the form f sr ( t ) = 12 π i (cid:90) ∞ ψ (cid:0) ρ e − i π (cid:1) ψ ( ρ e − i π ) + ρ β e − i βπ (cid:16) b b + ρ η e − i ηπ (cid:17) − ψ (cid:0) ρ e i π (cid:1) ψ ( ρ e i π ) + ρ β e i βπ (cid:16) b b + ρ η e i ηπ (cid:17) e − ρt d ρ = 1 π (cid:90) ∞ Q ( ρ ) (cid:12)(cid:12)(cid:12) ψ ( ρ e i π ) + ρ β e i βπ (cid:16) b b + ρ η e i ηπ (cid:17)(cid:12)(cid:12)(cid:12) ρ β e − ρt d ρ, with function Q given by (49), is calculated as the inverse Laplace transform L − (cid:104) ˜ f sr ( s ) (cid:105) = 12 π i lim R →∞ ,r → (cid:90) Γ ˜ f sr ( s ) e st d s, (151)see (131), using the Cauchy integral theorem (cid:73) Γ (I) ˜ f sr ( s ) e st d s = 0 , where the contour Γ (I) = Γ ∪ Γ ∪ Γ ∪ Γ ∪ Γ ∪ Γ ∪ Γ ∪ Γ is chosen as in Figure 8, since the integrals alongcontours Γ (parametrized by s = ρ e i π , ρ ∈ ( r, R )) and Γ (parametrized by s = ρ e − i π , ρ ∈ ( r, R )) readlim R →∞ ,r → (cid:90) Γ ˜ f sr ( s ) e st d s = (cid:90) ∞ ψ (cid:0) ρ e i π (cid:1) ψ ( ρ e i π ) + ρ β e i βπ (cid:16) b b + ρ η e i ηπ (cid:17) e − ρt e i π d ρ = (cid:90) ∞ ψ (cid:0) ρ e i π (cid:1) ψ ( ρ e i π ) + ρ β e i βπ (cid:16) b b + ρ η e i ηπ (cid:17) e − ρt d ρ, (152)lim R →∞ ,r → (cid:90) Γ ˜ f sr ( s ) e st d s = (cid:90) ∞ ψ (cid:0) ρ e − i π (cid:1) ψ ( ρ e − i π ) + ρ β e − i βπ (cid:16) b b + ρ η e − i ηπ (cid:17) e − ρt e − i π d ρ = − (cid:90) ∞ ψ (cid:0) ρ e − i π (cid:1) ψ ( ρ e − i π ) + ρ β e − i βπ (cid:16) b b + ρ η e − i ηπ (cid:17) e − ρt d ρ, (153)respectively, while the integrals along Γ , Γ , Γ , Γ , Γ tend to zero as R → ∞ and r → . The contour Γ is parametrized by s = p + i R, p ∈ (0 , p ) , with R → ∞ , so that the integral I Γ = (cid:90) Γ ˜ f sr ( s ) e st d s = (cid:90) p ψ ( p + i R ) ψ ( p + i R ) + ( p + i R ) β (cid:16) b b + ( p + i R ) η (cid:17) e ( p +i R ) t d p is estimated as | I Γ | ≤ (cid:90) p | ψ ( p + i R ) | (cid:12)(cid:12)(cid:12) ψ ( p + i R ) + ( p + i R ) β (cid:16) b b + ( p + i R ) η (cid:17)(cid:12)(cid:12)(cid:12) e pt d p. (154)Assuming s = ρ e i ϕ , since R → ∞ , one obtains ρ = (cid:112) p + R ∼ R and ϕ = arctan Rp ∼ π , so that (154)becomes lim R →∞ | I Γ | ≤ (cid:18) a a − b b (cid:19) lim R →∞ (cid:90) p R η e pt d p = 0 , due to (141). Analogously, it can be proved that lim R →∞ | I Γ | = 0 . The integral along contour Γ , parametrized by s = R e i ϕ , ϕ ∈ (cid:0) π , π (cid:1) , with R → ∞ , is I Γ = (cid:90) Γ ˜ f sr ( s ) e st d s = (cid:90) π π ψ (cid:0) R e i ϕ (cid:1) ψ ( R e i ϕ ) + R β e i βϕ (cid:16) b b + R η e i ηϕ (cid:17) e Rt e i ϕ i R e i ϕ d ϕ, | I Γ | ≤ (cid:90) π π (cid:12)(cid:12) ψ (cid:0) R e i ϕ (cid:1)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ψ ( R e i ϕ ) + R β e i βϕ (cid:16) b b + R η e i ηϕ (cid:17)(cid:12)(cid:12)(cid:12) R e Rt cos ϕ d ϕ. (155)Similarly as in (141), (cid:12)(cid:12) ψ (cid:0) R e i ϕ (cid:1)(cid:12)(cid:12) ∼ (cid:16) a a − b b (cid:17) R β , as R → ∞ , which along with cos ϕ < ϕ ∈ (cid:0) π , π (cid:1) impliesthat (155) in the limit when R → ∞ becomeslim R →∞ | I Γ | ≤ (cid:18) a a − b b (cid:19) lim R →∞ (cid:90) π π R − η e Rt cos ϕ d ϕ = 0 . By the similar arguments, lim R →∞ | I Γ | = 0 , as well.Parametrization of the contour Γ is s = r e i ϕ , ϕ ∈ ( − π, π ) , with r → , so that I Γ = (cid:90) Γ ˜ f sr ( s ) e st d s = (cid:90) π − π ψ (cid:0) r e i ϕ (cid:1) ψ ( r e i ϕ ) + r β e i βϕ (cid:16) b b + r η e i ηϕ (cid:17) e rt e i ϕ i r e i ϕ d ϕ can be estimated bylim r → | I Γ | ≤ lim r → (cid:90) π − π (cid:12)(cid:12) ψ (cid:0) r e i ϕ (cid:1)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ψ ( r e i ϕ ) + r β e i βϕ (cid:16) b b + r η e i ηϕ (cid:17)(cid:12)(cid:12)(cid:12) r e rt cos ϕ d ϕ ≤ π lim r → r = 0 , since (cid:12)(cid:12) ψ (cid:0) r e i ϕ (cid:1)(cid:12)(cid:12) ∼ a as r → . B.2.2 Case when function Ψ has a negative real zeroRelaxation modulus in the case of models having zero glass compliance. The Cauchy integral theoremwith the relaxation modulus in complex domain (53) as an integrand becomes (cid:73) Γ (II) ˜ σ sr ( s ) e st d s = 0 , (156)where the contour Γ (II) = Γ ∪ Γ ∪ Γ ∪ Γ a ∪ Γ ∗ ∪ Γ b ∪ Γ ∪ Γ a ∪ Γ ∗ ∪ Γ b ∪ Γ ∪ Γ is chosen as in Figure 9due to the existence of negative real zero − ρ ∗ , with ρ ∗ determined from (61), of function Ψ , given by (38).Figure 9: Contour Γ (II) . The relaxation modulus is obtained in the form σ sr ( t ) = 1 π (cid:90) ∞ K ( ρ ) | Ψ ( ρ e i π ) | e − ρt ρ − µ d ρ + f ∗ sr ( ρ ∗ ) e − ρ ∗ t , (157)with functions K and f ∗ sr given by (47) and (62), using the Cauchy integral theorem (156). Namely, the contourΓ a is parametrized by s = ρ e i π , ρ ∈ ( ρ ∗ + r ∗ , R ) , while Γ b has the same parametrization with ρ ∈ ( ρ ∗ − r ∗ , r ),so that in the limit when R → ∞ , r → , and r ∗ → a ∪ Γ b , as in (145), readslim R →∞ ,r → ,r ∗ → (cid:90) Γ a ∪ Γ b ˜ σ sr ( s ) e st d s = − (cid:90) ∞ b + b ρ η e i ηπ Ψ ( ρ e i π ) e i µπ e − ρt ρ − µ d ρ, (158)37nd similarly the integral along contour Γ a ∪ Γ b (parametrized by s = ρ e − i π , ρ ∈ ( ρ ∗ − r ∗ , r ) for Γ a and ρ ∈ ( ρ ∗ + r ∗ , R ) for Γ b ), as in (146), islim R →∞ ,r → ,r ∗ → (cid:90) Γ a ∪ Γ b ˜ σ sr ( s ) e st d s = (cid:90) ∞ b + b ρ η e − i ηπ Ψ ( ρ e − i π ) e − i µπ e − ρt ρ − µ d ρ, (159)so that, along with the inverse Laplace transform (147), the integrals in (158) and (159) yield the first term inthe relaxation modulus (157), while, as it will be calculated, the integrals along contours Γ ∗ and Γ ∗ yield thesecond term in (157), as the integrals along Γ , Γ , Γ , Γ , Γ tend to zero as R → ∞ and r → , as alreadyproved in Section B.2.1.The contour Γ ∗ is parametrized by s − ρ ∗ e i π = r ∗ e i ϕ , ϕ ∈ (0 , π ) , with r ∗ → , and the corresponding integralreads I Γ ∗ = (cid:90) Γ ∗ ˜ σ sr ( s ) e st d s = (cid:90) π ρ ∗ e i π + r ∗ e i ϕ ) − µ b + b (cid:0) ρ ∗ e i π + r ∗ e i ϕ (cid:1) η Ψ ( ρ ∗ e i π + r ∗ e i ϕ ) e( ρ ∗ e i π + r ∗ e i ϕ ) t i r ∗ e i ϕ d ϕ, so that by letting r ∗ → r ∗ → I Γ ∗ = − i b + b ( ρ ∗ ) η e i ηπ ( ρ ∗ ) − µ e i(1 − µ ) π e − ρ ∗ t lim r ∗ → (cid:90) π r ∗ e i ϕ Ψ ( ρ ∗ e i π + r ∗ e i ϕ ) d ϕ, = − i π ( ρ ∗ ) µ e i µπ b + b ( ρ ∗ ) η e i ηπ αa ( ρ ∗ ) α e i απ + βa ( ρ ∗ ) β e i βπ + γa ( ρ ∗ ) γ e i γπ e − ρ ∗ t . (160)In calculating (160), function Ψ (38) for r ∗ → (cid:0) ρ ∗ e i π + r ∗ e i ϕ (cid:1) = 1 + a ( ρ ∗ ) α e i απ (cid:18) − r ∗ e i ϕ ρ ∗ (cid:19) α + a ( ρ ∗ ) β e i βπ (cid:18) − r ∗ e i ϕ ρ ∗ (cid:19) β + a ( ρ ∗ ) γ e i γπ (cid:18) − r ∗ e i ϕ ρ ∗ (cid:19) γ ≈ a ( ρ ∗ ) α e i απ (cid:18) − α r ∗ e i ϕ ρ ∗ (cid:19) + a ( ρ ∗ ) β e i βπ (cid:18) − β r ∗ e i ϕ ρ ∗ (cid:19) + a ( ρ ∗ ) γ e i γπ (cid:18) − γ r ∗ e i ϕ ρ ∗ (cid:19) ≈ − r ∗ e i ϕ ρ ∗ (cid:16) αa ( ρ ∗ ) α e i απ + βa ( ρ ∗ ) β e i βπ + γa ( ρ ∗ ) γ e i γπ (cid:17) , (161)where the approximation (1 + x ) ξ ≈ ξx, for | x | (cid:28) , and fact that Ψ (cid:0) ρ ∗ e i π (cid:1) = 0 (since ρ ∗ is negative realzero of Ψ) are used, implyinglim r ∗ → r ∗ e i ϕ Ψ ( ρ ∗ e i π + r ∗ e i ϕ ) = − ρ ∗ αa ( ρ ∗ ) α e i απ + βa ( ρ ∗ ) β e i βπ + γa ( ρ ∗ ) γ e i γπ . (162)The integral corresponding to contour Γ ∗ , parametrized by s − ρ ∗ e − i π = r ∗ e i ϕ , ϕ ∈ ( − π, , in the limit when r ∗ → r ∗ → I Γ ∗ = lim r ∗ → (cid:90) Γ ∗ ˜ σ sr ( s ) e st d s = i b + b ( ρ ∗ ) η e − i ηπ ( ρ ∗ ) − µ e − i(1 − µ ) π e − ρ ∗ t lim r ∗ → (cid:90) − π r ∗ e i ϕ Ψ ( ρ ∗ e − i π + r ∗ e i ϕ ) d ϕ = − i π ( ρ ∗ ) µ e − i µπ b + b ( ρ ∗ ) η e − i ηπ αa ( ρ ∗ ) α e − i απ + βa ( ρ ∗ ) β e − i βπ + γa ( ρ ∗ ) γ e − i γπ e − ρ ∗ t , using the similar procedure as in calculating lim r ∗ → I Γ ∗ .Therefore, one haslim r ∗ → ( I Γ ∗ + I Γ ∗ ) = − i π (cid:32) b + b ( ρ ∗ ) η e i ηπ αa ( ρ ∗ ) α e i απ + βa ( ρ ∗ ) β e i βπ + γa ( ρ ∗ ) γ e i γπ e i µπ + b + b ( ρ ∗ ) η e − i ηπ αa ( ρ ∗ ) α e − i απ + βa ( ρ ∗ ) β e − i βπ + γa ( ρ ∗ ) γ e − i γπ e − i µπ (cid:33) ( ρ ∗ ) µ e − ρ ∗ t , yielding function f ∗ sr in the form (62) and representing the second term in (157).38 elaxation modulus in the case of models having non-zero glass compliance. Inverting theLaplace transform in the relaxation modulus in complex domain (54), one obtains σ sr ( t ) = b a − b a (cid:90) t L − (cid:104) ˜ f sr ( s ) (cid:105) ( τ ) d τ = b a − b a (cid:90) t f sr ( τ ) d τ + b a f ∗ sr ( ρ ∗ ) (cid:16) − e − ρ ∗ t (cid:17) , (163)where functions ˜ f sr , f sr , and f ∗ sr are respectively given by (150), (65), and (67), due to the existence of negativereal zero − ρ ∗ of function Ψ (38), with ρ ∗ determined from (61) with γ = β + η. Namely, using function ˜ f sr inthe form ˜ f sr ( s ) = ψ ( s ) ψ ( s ) + s β (cid:16) b b + s η (cid:17) = 1 + a s α + a (cid:16) a a − b b (cid:17) s β Ψ ( s ) , (164)see (39) and (38), as an integrand in the Cauchy integral theorem (cid:73) Γ (II) ˜ f sr ( s ) e st d s = 0 , (165)where the contour Γ (II) = Γ ∪ Γ ∪ Γ ∪ Γ a ∪ Γ ∗ ∪ Γ b ∪ Γ ∪ Γ a ∪ Γ ∗ ∪ Γ b ∪ Γ ∪ Γ is chosen as in Figure 9,one obtains the second and third term in (163).The contour Γ a is parametrized by s = ρ e i π , ρ ∈ ( ρ ∗ + r ∗ , R ) , while Γ b has the same parametrization with ρ ∈ ( ρ ∗ − r ∗ , r ), so that in the limit when R → ∞ , r → r ∗ → a ∪ Γ b , asin (152), reads lim R →∞ ,r → ,r ∗ → (cid:90) Γ a ∪ Γ b ˜ f sr ( s ) e st d s = (cid:90) ∞ ψ (cid:0) ρ e i π (cid:1) ψ ( ρ e i π ) + ρ β e i βπ (cid:16) b b + ρ η e i ηπ (cid:17) e − ρt d ρ. (166)Similarly, the integral along contour Γ a ∪ Γ b (parametrized by s = ρ e − i π , ρ ∈ ( ρ ∗ − r ∗ , r ) for Γ a and ρ ∈ ( ρ ∗ + r ∗ , R ) for Γ b ), as in (153), islim R →∞ ,r → ,r ∗ → (cid:90) Γ a ∪ Γ b ˜ f sr ( s ) e st d s = − (cid:90) ∞ ψ (cid:0) ρ e − i π (cid:1) ψ ( ρ e − i π ) + ρ β e − i βπ (cid:16) b b + ρ η e − i ηπ (cid:17) e − ρt d ρ. (167)The contour Γ ∗ is parametrized by s − ρ ∗ e i π = r ∗ e i ϕ , ϕ ∈ (0 , π ) , with r ∗ → , so that the correspondingintegral reads I Γ ∗ = (cid:90) Γ ∗ ˜ f sr ( s ) e st d s = (cid:90) π a (cid:0) ρ ∗ e i π + r ∗ e i ϕ (cid:1) α + a (cid:16) a a − b b (cid:17) (cid:0) ρ ∗ e i π + r ∗ e i ϕ (cid:1) β Ψ ( ρ ∗ e i π + r ∗ e i ϕ ) e( ρ ∗ e i π + r ∗ e i ϕ ) t i r ∗ e i ϕ d ϕ, so that by letting r ∗ → r ∗ → I Γ ∗ = − i (cid:18) a ( ρ ∗ ) α e i απ + a (cid:18) a a − b b (cid:19) ( ρ ∗ ) β e i βπ (cid:19) e − ρ ∗ t lim r ∗ → (cid:90) π r ∗ e i ϕ Ψ ( ρ ∗ e i π + r ∗ e i ϕ ) d ϕ, = i π a ( ρ ∗ ) α e i απ + a (cid:16) a a − b b (cid:17) ( ρ ∗ ) β e i βπ αa ( ρ ∗ ) α e i απ + βa ( ρ ∗ ) β e i βπ + ( β + η ) a ( ρ ∗ ) β + η e i( β + η ) π ρ ∗ e − ρ ∗ t . (168)In calculating (168), the expression (162), with γ = β + η, is used. The integral corresponding to contour Γ ∗ ,parametrized by s − ρ ∗ e − i π = r ∗ e i ϕ , ϕ ∈ ( − π, , in the limit when r ∗ → r ∗ → I Γ ∗ = lim r ∗ → (cid:90) Γ ∗ ˜ f sr ( s ) e st d s = i (cid:18) a ( ρ ∗ ) α e − i απ + a (cid:18) a a − b b (cid:19) ( ρ ∗ ) β e − i βπ (cid:19) e − ρ ∗ t lim r ∗ → (cid:90) − π r ∗ e i ϕ Ψ ( ρ ∗ e − i π + r ∗ e i ϕ ) d ϕ = i π a ( ρ ∗ ) α e − i απ + a (cid:16) a a − b b (cid:17) ( ρ ∗ ) β e − i βπ αa ( ρ ∗ ) α e − i απ + βa ( ρ ∗ ) β e − i βπ + ( β + η ) a ( ρ ∗ ) β + η e − i( β + η ) π ρ ∗ e − ρ ∗ t , (169)39sing the similar procedure as in calculating lim r ∗ → I Γ ∗ .In the limit when R → ∞ , r → , and r ∗ → , the inverse Laplace transform L − (cid:104) ˜ f sr ( s ) (cid:105) (151), i.e., thesecond and the third term in the relaxation modulus (163), is obtained from the Cauchy integral theorem (165)as the sum of integrals along contours Γ a ∪ Γ b , Γ a ∪ Γ b , Γ ∗ , and Γ ∗ , respectively given by (166), (167), (168),and (169), since the integrals along Γ , Γ , Γ , Γ , Γ tend to zero as R → ∞ and r → , as already proved inSection B.2.1. Therefore, one has L − (cid:104) ˜ f sr ( s ) (cid:105) ( t ) = f sr ( t ) − ρ ∗ f ∗ sr ( ρ ∗ ) e − ρ ∗ t , with f sr given by (65) and f ∗ sr ( ρ ∗ ) = 12 a ( ρ ∗ ) α e i απ + a (cid:16) a a − b b (cid:17) ( ρ ∗ ) β e i βπ αa ( ρ ∗ ) α e i απ + βa ( ρ ∗ ) β e i βπ + ( β + η ) a ( ρ ∗ ) β + η e i( β + η ) π + 1 + a ( ρ ∗ ) α e − i απ + a (cid:16) a a − b b (cid:17) ( ρ ∗ ) β e − i βπ αa ( ρ ∗ ) α e − i απ + βa ( ρ ∗ ) β e − i βπ + ( β + η ) a ( ρ ∗ ) β + η e − i( β + η ) π , yielding (67). B.2.3 Case when function Ψ has a pair of complex conjugated zerosRelaxation modulus in the case of models having zero glass compliance. The Cauchy residue theoremwith the relaxation modulus in complex domain (53) as an integrand takes the form (cid:73) Γ (I) ˜ σ sr ( s ) e st d s = 2 π i (cid:0) Res (cid:0) ˜ σ sr ( s ) e st , s (cid:1) + Res (cid:0) ˜ σ sr ( s ) e st , ¯ s (cid:1)(cid:1) , (170)due to the existence of complex conjugated zeroes s and ¯ s having negative real part of function Ψ , given by(38), with the contour Γ (I) = Γ ∪ Γ ∪ Γ ∪ Γ ∪ Γ ∪ Γ ∪ Γ ∪ Γ chosen as in Figure 8.The relaxation modulus is obtained in the form σ sr ( t ) = 1 π (cid:90) ∞ K ( ρ ) | Ψ ( ρ e i π ) | e − ρt ρ − µ d ρ + f ( r ) sr ( t ) , (171)with functions K and f ( r ) sr given by (47) and (63), using the Cauchy residue theorem (170). Namely, theintegration along the contour Γ (I) yields the inverse Laplace transform (147) and the first term in the relaxationmodulus (171), that is already obtained in Section B.2.1, while the second term in (171) consists of the residuesof function ˜ σ sr ( s ) e st , since s = ρ e i ϕ and ¯ s = ρ e − i ϕ , ϕ ∈ (cid:0) π , π (cid:1) , are its poles of the first order, i.e., thefirst order zeros of function Ψ (38). Therefore, one has f ( r ) sr ( t ) = (cid:34) s − µ b + b s η dd s Ψ ( s ) e st (cid:35) s = s + (cid:34) s − µ b + b s η dd s Ψ ( s ) e st (cid:35) s =¯ s = (cid:32) b + b s η αa s α − + βa s β − + γa s γ − (cid:12)(cid:12)(cid:12)(cid:12) s = ρ e i ϕ e − i(1 − µ ) ϕ ρ − µ e i ρ t sin ϕ + b + b s η αa s α − + βa s β − + γa s γ − (cid:12)(cid:12)(cid:12)(cid:12) s = ρ e − i ϕ e i(1 − µ ) ϕ ρ − µ e i ρ t sin ϕ (cid:33) e ρ t cos ϕ , yielding (63). Relaxation modulus in the case of models having non-zero glass compliance.
Inverting theLaplace transform in the relaxation modulus in complex domain (54), one obtains σ sr ( t ) = b a − b a (cid:90) t L − (cid:104) ˜ f sr ( s ) (cid:105) ( τ ) d τ = b a − b a (cid:90) t f sr ( τ ) d τ − b a (cid:90) t f ( r ) sr ( τ ) d τ , (172)where functions ˜ f sr , f sr , and f ( r ) sr are respectively given by (150), (65), and (68), due to the existence of complexconjugated zeroes s and ¯ s of function Ψ (38) having negative real part.40sing function ˜ f sr in the form (164) as an integrand in the Cauchy residue theorem (cid:73) Γ (I) ˜ f sr ( s ) e st d s = 2 π i (cid:16) Res (cid:16) ˜ f sr ( s ) e st , s (cid:17) + Res (cid:16) ˜ f sr ( s ) e st , ¯ s (cid:17)(cid:17) , with the contour Γ (I) = Γ ∪ Γ ∪ Γ ∪ Γ ∪ Γ ∪ Γ ∪ Γ ∪ Γ chosen as in Figure 8, one obtains the inverseLaplace transform L − (cid:104) ˜ f sr ( s ) (cid:105) (151), i.e., the second and the third term in (172), in the form L − (cid:104) ˜ f sr ( s ) (cid:105) ( t ) = f sr ( t ) + f ( r ) sr ( t ) . (173)The first term in (173), being a consequence of the integration along the contour Γ (I) , is already obtainedin Section B.2.1, while the second one consists of the residues of function ˜ f sr ( s ) e st , since s = ρ e i ϕ and¯ s = ρ e − i ϕ , ϕ ∈ (cid:0) π , π (cid:1) , are its poles of the first order, i.e., the first order zeros of function Ψ , as proved inSection A.1. Therefore in (173), one has f sr given by (65) and f ( r ) sr ( t ) = a s α + a (cid:16) a a − b b (cid:17) s β dd s Ψ ( s ) e st s = s + a s α + a (cid:16) a a − b b (cid:17) s β dd s Ψ ( s ) e st s =¯ s = a s α + a (cid:16) a a − b b (cid:17) s β αa s α − + βa s β − + ( β + η ) a s β + η − (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) s = ρ e i ϕ e i ρ t sin ϕ + 1 + a s α + a (cid:16) a a − b b (cid:17) s β αa s α − + βa s β − + ( β + η ) a s β + η − (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) s = ρ e − i ϕ e − i ρ t sin ϕ e ρ t cos ϕ , yielding (68). Acknowledgment
This work is supported by the Serbian Ministry of Education, Science and Technological Development undergrant 174005, as well as by the Provincial Secretariat for Higher Education and Scientific Research under grant142 − − / References [1] T. M. Atanackovic, S. Pilipovic, B. Stankovic, and D. Zorica.
Fractional Calculus with Applications inMechanics: Vibrations and Diffusion Processes . Wiley-ISTE, London, 2014.[2] T. M. Atanackovic, S. Pilipovic, B. Stankovic, and D. Zorica.
Fractional Calculus with Applications inMechanics: Wave Propagation, Impact and Variational Principles . Wiley-ISTE, London, 2014.[3] T. M. Atanackovic, S. Pilipovic, and D. Zorica. Distributed-order fractional wave equation on a finitedomain: creep and forced oscillations of a rod.
Continuum Mechanics and Thermodynamics , 23:305–318,2011.[4] T. M. Atanackovic, S. Pilipovic, and D. Zorica. Distributed-order fractional wave equation on a finitedomain. Stress relaxation in a rod.
International Journal of Engineering Science , 49:175–190, 2011.[5] R. L. Bagley and P. J. Torvik. On the fractional calculus model of viscoelastic behavior.
Journal ofRheology , 30:133–155, 1986.[6] E. Bazhlekova and I. Bazhlekov. Complete monotonicity of the relaxation moduli of distributed-orderfractional Zener model. In
Proceedings of the 44 th International Conference on Applications of Mathematicsin Engineering and Economics, AIP Conference Proceedings 2048 , pages 050008–1–8, 2018.[7] E. Bazhlekova and K. Tsocheva. Fractional Burgers’ model: thermodynamic constraints and completelymonotonic relaxation function.
Comptes rendus de l’Acad´emie bulgare des Sciences , 69:825–834, 2016.[8] C. Celauro, C. Fecarotti, A. Pirrotta, and A. C. Collop. Experimental validation of a fractional model forcreep/recovery testing of asphalt mixtures.
Construction and Building Materials , 36:458–466, 2012.419] W. N. Findley, J. S. Lai, and K. Onaran.
Creep and Relaxation of Nonlinear Viscoelastic Materials - Withan Introduction to Linear Viscoelasticity . Dover Publications, New York, 1976.[10] R. Gorenflo and F. Mainardi.
Fractional Calculus: Integral and Differential Equations of Fractional Order.In: Fractals and Fractional Calculus in Continuum Mechanics (eds A. Carpinteri, F. Mainardi) , volume378 of
CISM Courses and Lecture Notes . Springer Verlag, Wien and New York, 1997.[11] A. A. Kilbas, H. M. Srivastava, and J. J. Trujillo.
Theory and Applications of Fractional DifferentialEquations . Elsevier B.V., Amsterdam, 2006.[12] F. Mainardi.
Fractional Calculus and Waves in Linear Viscoelasticity . Imperial College Press, London,2010.[13] F. Mainardi and G. Spada. Creep, relaxation and viscosity properties for basic fractional models in rheology.
European Physical Journal Special Topics , 193:133–160, 2011.[14] N. Makris. The frequency response function of the creep compliance.
Meccanica , 2018.[15] M. Oeser, T. Pellinen, T. Scarpas, and C. Kasbergen. Studies on creep and recovery of rheological bodiesbased upon conventional and fractional formulations and their application on asphalt mixture.
InternationalJournal of Pavement Engineering , 9:373–386, 2008.[16] A. S. Okuka and D. Zorica. Formulation of thermodynamically consistent fractional Burgers models.
ActaMechanica , 229:3557–3570, 2018.[17] Yu. A. Rossikhin and M. V. Shitikova. Analysis of rheological equations involving more than one fractionalparameters by the use of the simplest mechanical systems based on these equations.
Mechanics of Time-Dependent Materials , 5:131–175, 2001.[18] Yu. A. Rossikhin and M. V. Shitikova. Analysis of the viscoelastic rod dynamics via models involvingfractional derivatives or operators of two different orders.
Shock and Vibration Digest , 36:3–26, 2004.[19] Yu. A. Rossikhin and M. V. Shitikova. Free damped vibrations of a viscoelastic oscillator based on Rabot-nov’s model.
Mechanics of Time-Dependent Materials , 12:129–149, 2008.[20] A. Zbiciak. Mathematical description of rheological properties of asphalt-aggregate mixes.