# Viscoelastic Fractional Model Based on Harmonic Excitation

Authors: Ruiqun Ma; Min Ni; Quanlong Chen; Yinjia Zhou; Jinglong Han
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

This paper analyzes the viscoelastic fractional constitutive model based on the analytical form of the fractional derivative of the sine (cosine) function. First, the polynomials of the linear model are combined into a fractional derivative term, then the complex modulus can be easily obtained, and the response characteristics of the model can be analyzed. Second, the fractional nonlinear model is expressed in harmonic form, which can easily identify nonlinear parameters. Finally, fractional-order nonlinear models can also be transformed into integer-order nonlinear forms. In this paper, the transformation model is verified by literature test data, which show the correctness and accuracy of the transformation form.


---

## Body

## 1. Introduction

Fractional calculus has a long history, dating back to 1695 after the birth of integer calculus. However, in the next nearly 300 years, fractional calculus developed slowly. Until the end of the 19th century, the theoretical system of fractional calculus was basically established and gradually developed and improved. In recent decades, fractional calculus has been used increasingly in various fields, such as materials science [1–3], biology [4, 5], chemistry [6, 7], and finance [8].Viscoelastic materials have the property of degrading memory, and the fractional-order model can describe this property well. Hence, the fractional-order model is more in line with the nature of the material’s performance. When describing an ideal linear elastic material, its load is expressed as a linear function of displacement; when describing Newtonian fluid materials, the load is expressed as a linear function of velocity. Load originally represents the viscoelastic properties of a material, usually in the form of linear springs combined with Newtonian dashpot. The use of this model often requires several springs and sticking pot elements to accurately describe the mechanical properties of the material, such as the generalized Kelvin model and the generalized Maxwell model. From another point of view, displacement is used to describe elasticity, and the first derivative of displacement is used to describe viscosity. Thus, the viscoelasticity between elasticity and viscosity can be described by the fractional derivative of displacement, and the order is between 0 and 1. Many scholars have studied the viscoelastic fractional model. References [9, 10] studied the mechanical properties of viscoelastic materials by using the fractional derivative model. Reference [11] used the fractional Kelvin–Voigt model and the high-order fractional derivative FVMP model to predict the frequency response curve of a rubber material and evaluated the dynamic mechanical properties of a rubber material. Amabili et al. also made significant contributions to the viscoelastic fractional model and its application [12–16], in which literature [13] studied the change of damping with amplitude under harmonic excitation by introducing geometric nonlinearity into the fractional Kelvin model. Reference [17] mentioned that Prony series could be applied to approximate numerical solutions of fractional-order models to any desired accuracy. Reference [18] proposed a fractional-order model that considered the geometric factors of viscoelastic damping system, and a numerical method was derived on the basis of matrix function theory and the Grumwald–Letnikov discrete form of fractional derivative to study the vibration response of viscoelastic suspension installed on heavy tracked vehicles. Reference [19] studied the fractional Maxwell model based on the generalized Caputo definition of fractional derivative.When calculating the complex modulus of fractional-order linear model, the frequency domain transformation method of integer-order linear model is usually used. However, after the frequency domain transformation of the fractional-order model, the fractional power of the imaginary number unit appeared on the denominator, wherein separating the real and imaginary parts is cumbersome [11, 20–23]. In the aspect of parameter identification, the identification method is not as mature as the identification method of integer-order system parameters because the unknown parameters contain derivative orders. With the application of fractional calculus to viscoelasticity, different methods, such as genetic algorithm method [24, 25], spectral parameter estimation method [26], error variable method [27], ant colony optimization method [28], curve fitting method [29–31], nonlinear optimization method [32], and error function method based on minimization [33–35], have been used for the identification of fractional model parameters. Reference [36] proposed a method to identify parameters using hysteresis loop, which was verified by fractional Kelvin–Voigt model and fractional Maxwell model. Reference [37] proposed a new identification method that determined unknown parameters by transforming from the time domain to the frequency domain combined with the least square method.In this paper, based on harmonic excitation, the complex modulus form of the fractional linear model is obtained conveniently by using the fractional derivative formula of the sine (cosine) function, which is completely consistent with the form obtained by the frequency domain transformation method. At the same time, analyzing the response characteristics of linear systems under harmonic excitation by applying the fractional derivative formula of the sine (cosine) function is very easy. The fractional nonlinear model is converted into harmonic form, and the differential evolution algorithm is used for parameter identification. The results show that the fractional nonlinear model is more accurate than the inelastic displacement field model. Also, the form of the model is simpler. Viscoelastic materials are often subjected to periodic loads in the damping mechanism. To facilitate the integration of fractional-order models and integer-order system equations, the harmonic form is further transformed into an integer-order form, and the effect is the same as that of the harmonic form.


## 2. Fractional Calculus

### 2.1. Definitions of Fractional Calculus

Fractional calculus has many definitions in different forms. Currently, the definitions of Grünwald–Letnikov, Riemann–Liouville, and Caputo are widely used. This paper mainly introduces the definitions of Grünwald–Letnikov and Riemann–Liouville. The Grünwald–Letnikov definition is extended from the formula of higher-order derivatives of integer order, and its form is expressed as follows:(1)GLt0Dtαft=limh⟶01hα∑j=0t−t0/h−1jαjft−jh,where ⋅ is the nearest integer, and −1jαj is the binomial coefficient. If the selected calculation step is small enough, then the limit operation in (1) can be ignored.(2)GLt0Dtαft≈1hα∑j=0t−t0/h−1jαjft−jh.The Riemann–Liouville definition is extended from integral order integral to noninteger order, and it is expressed as(3)RLt0Dt−αft=1Γα∫t0tt−τα−1fτdτ,where Γ represents the gamma function. When expressing fractional differential, we cannot simply replace α in the integral formula with −α. Instead, it should be given based on the definition. Given a function ft, its derivative expression under the Riemann–Liouville definition is expressed as follows:(4)RLt0Dtβft=dndtnRLt0Dt−n−βft=1Γn−βdndtn∫t0tt−τn−β−1fτdτ,where β>0 and n is the smallest integer greater than β.


### 2.2. Fractional Derivative of the Sine (Cosine) Function

The integer derivative of the sine (cosine) function [15] is(5)dndtnsinωt=ωnsinωt+nπ2,dndtncosωt=ωncosωt+nπ2.The fractional derivative expression of the sine (cosine) chord function is obtained by replacing its ordern with any real number α(6)Dtαsinωt=ωαsinωt+απ2,Dtαcosωt=ωαcosωt+απ2.Equation (6) is verified by the numerical solution of Grünwald–Letnikov definition and the analytical solution of Riemann–Liouville definition, as shown in Figure 1. The 0.5 derivative of sin(2t) is shown in the figure. The overall agreement is very good. Figure 2 shows the difference curve between the two definitions and the extended formula. The figure shows only minimal difference in the initial stage, which is caused by the memory in the definition of fractional calculus. This finding establishes that the fractional derivative of the sine (cosine) function can be extended from the form of the integral derivative. In addition, Figure 2 shows that the error curve of the Riemann–Liouville definition tends to 0 smoothly, whereas the error curve of the Grünwald–Letnikov definition tends to 0 in fluctuations but does not affect the accuracy of the definition.Figure 1
Numerical curve of the extended formula and two definitions.Figure 2
Difference curves between two definitions and extended formula.


## 2.1. Definitions of Fractional Calculus

Fractional calculus has many definitions in different forms. Currently, the definitions of Grünwald–Letnikov, Riemann–Liouville, and Caputo are widely used. This paper mainly introduces the definitions of Grünwald–Letnikov and Riemann–Liouville. The Grünwald–Letnikov definition is extended from the formula of higher-order derivatives of integer order, and its form is expressed as follows:(1)GLt0Dtαft=limh⟶01hα∑j=0t−t0/h−1jαjft−jh,where ⋅ is the nearest integer, and −1jαj is the binomial coefficient. If the selected calculation step is small enough, then the limit operation in (1) can be ignored.(2)GLt0Dtαft≈1hα∑j=0t−t0/h−1jαjft−jh.The Riemann–Liouville definition is extended from integral order integral to noninteger order, and it is expressed as(3)RLt0Dt−αft=1Γα∫t0tt−τα−1fτdτ,where Γ represents the gamma function. When expressing fractional differential, we cannot simply replace α in the integral formula with −α. Instead, it should be given based on the definition. Given a function ft, its derivative expression under the Riemann–Liouville definition is expressed as follows:(4)RLt0Dtβft=dndtnRLt0Dt−n−βft=1Γn−βdndtn∫t0tt−τn−β−1fτdτ,where β>0 and n is the smallest integer greater than β.


## 2.2. Fractional Derivative of the Sine (Cosine) Function

The integer derivative of the sine (cosine) function [15] is(5)dndtnsinωt=ωnsinωt+nπ2,dndtncosωt=ωncosωt+nπ2.The fractional derivative expression of the sine (cosine) chord function is obtained by replacing its ordern with any real number α(6)Dtαsinωt=ωαsinωt+απ2,Dtαcosωt=ωαcosωt+απ2.Equation (6) is verified by the numerical solution of Grünwald–Letnikov definition and the analytical solution of Riemann–Liouville definition, as shown in Figure 1. The 0.5 derivative of sin(2t) is shown in the figure. The overall agreement is very good. Figure 2 shows the difference curve between the two definitions and the extended formula. The figure shows only minimal difference in the initial stage, which is caused by the memory in the definition of fractional calculus. This finding establishes that the fractional derivative of the sine (cosine) function can be extended from the form of the integral derivative. In addition, Figure 2 shows that the error curve of the Riemann–Liouville definition tends to 0 smoothly, whereas the error curve of the Grünwald–Letnikov definition tends to 0 in fluctuations but does not affect the accuracy of the definition.Figure 1
Numerical curve of the extended formula and two definitions.Figure 2
Difference curves between two definitions and extended formula.


## 3. Viscoelastic Model

### 3.1. Integer Order Model

The viscoelastic constitutive models studied in the early stage are mostly integer-order models, which are relatively intuitive; that is, they are obtained by connecting springs and dashpot in series and parallel. The advantage of this model is that it has a very regular form [38] when establishing the equations of motion of the system. The spring is the simplest elastic model, and the stressstrain relationship is expressed as(7)σ=Eε,where σ and ε are the normal stress and normal strain, respectively, and E is the elastic modulus. Dashpot is the simplest viscous fluid model, which follows Newton’s law of viscosity:(8)σ=ηε˙,where η is the viscosity coefficient, and ε˙ represents the derivative of ε with respect to time, that is, the strain rate. A spring and a dashpot are connected in series to form a Maxwell model, as shown in Figure 3. From (7) and (8) and the strain in series is equal to the sum of the strains of each part, the following equation of the Maxwell model can be obtained:(9)ε˙=σ˙E+σηorσ+p1σ˙=q1ε˙,where p1=η/E and q1=η. A spring and a dashpot are connected in parallel to form the Kelvin model (Figure 4). From (7) and (8) and the stress in parallel is equal to the sum of the stress of each part, the constitutive equation of the Kelvin model can be obtained:(10)σ=Eε+ηε˙orσ=q0ε+q1ε˙,where q0=E and q1=η.Figure 3
Maxwell model.Figure 4
Kelvin model.The Maxwell model is essentially closer to fluid, which is called Maxwell fluid. The Kelvin model can represent the creep of solid, but it cannot represent the instantaneous elasticity and relaxation of solid. Thus, it is not a viscoelastic solid in a complete sense. Neither the Maxwell model nor the Kelvin model alone can express the basic characteristics of viscoelastic solids: transient elasticity, creep, and relaxation. To express these characteristics, a larger number of elements is required. The simplest is the three-element solid model. The three-element model connects a spring in series with the Kelvin model but cannot represent instantaneous elasticity (Figure5). It is also known as the standard linear solid model, and its constitutive equation is expressed as(11)E1+E2σ+η1σ˙=E1E2ε+η1E2ε˙,or normal form(12)σ+p1σ˙=q0ε+q1ε˙,where p1=η1/E1+E2, q0=E1E2/E1+E2, and q1=η1E2/E1+E2. The Maxwell model and Kelvin model are connected in series to form a four-element Burgers model (Figure 6). Its normative constitutive equation is(13)σ+p1σ˙+p2σ¨=q1ε˙+q2ε¨,where p1=η2/E1+η2+η3/E3, p2=η2η3/E1E3, q1=η2, and q2=η2η3/E3.Figure 5
Three-element model.Figure 6
Burgers model.From (13), we can further generalize to the general form of the differential viscoelastic constitutive equation:(14)Pσ=Qε,where P and Q are differential operators.(15)P=∑m=0Mpmdmdtm,Q=∑n=0Nqndndtn,where pm and qn are material parameters. dm/dtm represents the mth derivative of the stress or strain function with respect to time.


### 3.2. Fractional Order Model

The fractional derivative constitutive model replaces the dashpot element in the combined spring-dashpot model with an Abel dashpot. The constitutive relationship of Abel dashpot is(16)σ=ηDγε,where η is the viscosity coefficient, Dγ is the fractional derivative operator, and 0<γ<1. Referring to the general form (14) of integral order constitutive equation, the general form of the constitutive equation expressed by the fractional derivative is(17)∑m=0MamDαmσ=∑n=0NbnDβnε,where αm and βn are the material parameters, and both are greater than 0 and less than 1. When αm and βn are positive integers, equation (17) becomes the integral derivative constitutive equation.The fractional Maxwell model is shown in Figure7. According to the strain of the series being equal to the sum of the strains of each part, the constitutive equation of the fractional Maxwell model can be obtained.(18)σ+ηEDγσ=ηDγε.Figure 7
Fractional Maxwell model.Figure8 shows the fractional Kelvin model, and the corresponding constitutive equation can also be obtained as(19)σ=Eε+ηDγε.Figure 8
Fractional Kelvin model.The fractional three-element model is shown in Figure9, and its constitutive equation is(20)E1+E2σ+η1Dγσ=E1E2ε+E2η1Dγε.Figure 9
Fractional three-element model.The fractional Burgers model contains two different Abel dashpots, and its constitutive model is Figure10(21)E1E3σ+E1+E3η2Dασ+E1η3Dβσ+η2η3Dα+βσ=E1E3η2Dαε+E1η2η3Dα+βε.Figure 10
Fractional Burgers model.


## 3.1. Integer Order Model

The viscoelastic constitutive models studied in the early stage are mostly integer-order models, which are relatively intuitive; that is, they are obtained by connecting springs and dashpot in series and parallel. The advantage of this model is that it has a very regular form [38] when establishing the equations of motion of the system. The spring is the simplest elastic model, and the stressstrain relationship is expressed as(7)σ=Eε,where σ and ε are the normal stress and normal strain, respectively, and E is the elastic modulus. Dashpot is the simplest viscous fluid model, which follows Newton’s law of viscosity:(8)σ=ηε˙,where η is the viscosity coefficient, and ε˙ represents the derivative of ε with respect to time, that is, the strain rate. A spring and a dashpot are connected in series to form a Maxwell model, as shown in Figure 3. From (7) and (8) and the strain in series is equal to the sum of the strains of each part, the following equation of the Maxwell model can be obtained:(9)ε˙=σ˙E+σηorσ+p1σ˙=q1ε˙,where p1=η/E and q1=η. A spring and a dashpot are connected in parallel to form the Kelvin model (Figure 4). From (7) and (8) and the stress in parallel is equal to the sum of the stress of each part, the constitutive equation of the Kelvin model can be obtained:(10)σ=Eε+ηε˙orσ=q0ε+q1ε˙,where q0=E and q1=η.Figure 3
Maxwell model.Figure 4
Kelvin model.The Maxwell model is essentially closer to fluid, which is called Maxwell fluid. The Kelvin model can represent the creep of solid, but it cannot represent the instantaneous elasticity and relaxation of solid. Thus, it is not a viscoelastic solid in a complete sense. Neither the Maxwell model nor the Kelvin model alone can express the basic characteristics of viscoelastic solids: transient elasticity, creep, and relaxation. To express these characteristics, a larger number of elements is required. The simplest is the three-element solid model. The three-element model connects a spring in series with the Kelvin model but cannot represent instantaneous elasticity (Figure5). It is also known as the standard linear solid model, and its constitutive equation is expressed as(11)E1+E2σ+η1σ˙=E1E2ε+η1E2ε˙,or normal form(12)σ+p1σ˙=q0ε+q1ε˙,where p1=η1/E1+E2, q0=E1E2/E1+E2, and q1=η1E2/E1+E2. The Maxwell model and Kelvin model are connected in series to form a four-element Burgers model (Figure 6). Its normative constitutive equation is(13)σ+p1σ˙+p2σ¨=q1ε˙+q2ε¨,where p1=η2/E1+η2+η3/E3, p2=η2η3/E1E3, q1=η2, and q2=η2η3/E3.Figure 5
Three-element model.Figure 6
Burgers model.From (13), we can further generalize to the general form of the differential viscoelastic constitutive equation:(14)Pσ=Qε,where P and Q are differential operators.(15)P=∑m=0Mpmdmdtm,Q=∑n=0Nqndndtn,where pm and qn are material parameters. dm/dtm represents the mth derivative of the stress or strain function with respect to time.


## 3.2. Fractional Order Model

The fractional derivative constitutive model replaces the dashpot element in the combined spring-dashpot model with an Abel dashpot. The constitutive relationship of Abel dashpot is(16)σ=ηDγε,where η is the viscosity coefficient, Dγ is the fractional derivative operator, and 0<γ<1. Referring to the general form (14) of integral order constitutive equation, the general form of the constitutive equation expressed by the fractional derivative is(17)∑m=0MamDαmσ=∑n=0NbnDβnε,where αm and βn are the material parameters, and both are greater than 0 and less than 1. When αm and βn are positive integers, equation (17) becomes the integral derivative constitutive equation.The fractional Maxwell model is shown in Figure7. According to the strain of the series being equal to the sum of the strains of each part, the constitutive equation of the fractional Maxwell model can be obtained.(18)σ+ηEDγσ=ηDγε.Figure 7
Fractional Maxwell model.Figure8 shows the fractional Kelvin model, and the corresponding constitutive equation can also be obtained as(19)σ=Eε+ηDγε.Figure 8
Fractional Kelvin model.The fractional three-element model is shown in Figure9, and its constitutive equation is(20)E1+E2σ+η1Dγσ=E1E2ε+E2η1Dγε.Figure 9
Fractional three-element model.The fractional Burgers model contains two different Abel dashpots, and its constitutive model is Figure10(21)E1E3σ+E1+E3η2Dασ+E1η3Dβσ+η2η3Dα+βσ=E1E3η2Dαε+E1η2η3Dα+βε.Figure 10
Fractional Burgers model.


## 4. Fractional Linear Model

### 4.1. Complex Modulus of Linear Model

The usual way to find the complex modulus is to perform Fourier transform or Laplace transform on the constitutive equation. When the fractional constitutive equation uses the two aforementioned transformations to find the complex modulus, the fractional power of the imaginary unit will be substituted, and the formula is more complicated when the real part and the imaginary part are separated. In this paper, the fractional derivative (6) of the sine (cosine) function is used to transform the linear model in the previous section, and the complex modulus of the model is obtained in a very convenient way. First, the fractional Kelvin model is taken as an example, assuming that the strain is a harmonic excitation ε=ε0sinωt, and ε0 is the excitation amplitude, and then equation (19) can be rewritten as follows:(22)σ=Eε0sinωt+ηωγε0sinωt+γπ2=E+ηωγcosγπ2ε0sinωt+ηωγsinγπ2ε0cosωt.Then, the storage modulus and dissipation modulus of the model are the coefficients in front ofε0sinωt and ε0cosωt, respectively,(23)Ys=E+ηωγcosγπ2Yl=ηωγsinγπ2.The loss angle is(24)δ=Arctanηωγsinγπ/2E+ηωγcosγπ/2.The above expressions are the same as the results obtained by Fourier transformation and separating the real and imaginary parts. In particular, whenγ=1, the resulting complex modulus and loss angle are(25)Ys=EYl=ηωδ=ArctanηωE.The above formula is completely consistent with the integer-order model. The fractional Maxwell model is slightly more complicated. Suppose that under harmonic excitationε=ε0sinωt, σ=σ0sinωt+φ is the stress response, and φ is the phase angle. Then, operate on the left and right sides of equation (18)(26)σ+ηEDγσ=σ01+ωγηEcosγπ2sinωt+φ+σ0ωγηEsinγπ2cosωt+φ=Aσ0sinωt+φ+BηDγε=ηωγε0sinωt+γπ2,where(27)A=1+ωγηEcosγπ22+ωγηEsinγπ22B=Arctanωγη/Esinγπ/21+ωγη/Ecosγπ/2.MoveA and B to the right side of the equation and rewrite it in the form of equation (22)(28)σ=ηωγAcosγπ2−Bε0sinωt+ηωγAsinγπ2−Bε0cosωt.The storage modulus, energy dissipation modulus, and loss angle of the fractional Maxwell model are further obtained as follows:(29)Ys=ηωγAcosγπ2−BYl=ηωγAsinγπ2−Bδ=γπ2−B.Whenγ=1,(30)A=1+ωηE2B=ArctanωηE.By substituting (30) into (29), the complex modulus and phase angle of the Maxwell model are expressed as follows:(31)Ys=Eη2ω2E2+η2ω2Yl=E2ηωE2+η2ω2δ=ArctanEωη.The obtained results are completely consistent with the results obtained by the integer-order Maxwell model by the Laplace transform method. In the same way, the expressions of storage modulus, dissipation modulus, and loss angle of the two other models can be obtained by the same method, which will not be repeated in this article.


### 4.2. Response Characteristics of System Based on the Linear Model

Analyzing the response characteristics of linear system under simple harmonic excitation by using the fractional derivative formula of sine (cosine) function in Section2.2 is very easy. The dynamic model shown in Figure 11 is studied [38], where rt represents the base excitation; c and η represent the damping coefficients of conventional linear damping and fractional derivative damping, respectively; α is the fractional order; and k is the stiffness coefficient of the linear spring. Using Newton’s second law, the equation of motion in Figure 11 is(32)mx¨+cx˙+μDαx+kx=cr˙+μDαr+kr.Figure 11
Dynamic model of the fractional-order linear system under basic excitation.When the basic excitation is sinusoidal excitationrt=Rsinωt, where R is the excitation amplitude and ω is the excitation frequency. Then, we assume that the response has the form(33)xt=Xsinωt+ϕ.Transform equation (32) according to equation (6)(34)−mω2+μωαcosαπ2+kXsinωt+φ+cω+μωαsinαπ2Xcosωt+φ=μωαcosαπ2+kRsinωt+cω+μωαsinαπ2Rcosωt.Further, transform this equation to obtain(35)A1Xsinωt+φ+B1=A2Rsinωt+B2,where(36)A1=−mω2+μωαcosαπ2+k2+cω+μωαsinαπ22B1=Arctancω+μωαsinαπ/2−mω2+μωαcosαπ/2+kA2=μωαcosαπ2+k2+cω+μωαsinαπ22B2=Arctancω+μωαsinαπ/2μωαcosαπ/2+k.MoveA1 and B1 on the left side of equation (35) to the right to obtain the response expression(37)xt=Xsinωt+φ=A2A1Rsinωt+B2−B1.Then, the magnitude and phase angle of the response are(38)X=A2A1Rφ=B2−B1.When the system parameters arem=5, k=10, c=0.5, and μ=1.5, transmissibility of the system is shown in Figure 12. Figure 12(a) is the three-dimensional transmissibility diagram; Figures 12(b)–12(d) are the transmissibility diagrams of α=0.5, α=1, and α=1.5, respectively. The three-dimensional diagram of the transmissibility in Figure 12(a) shows that when α is far from 1, the transmissibility of the system response will increase. In addition, in combination with Figures 12(b)–12(d), the peak point of the transmissibility is closer and closer to the position of ω=1 with the increase in α. Especially, when α=1, the obtained transmissibility curve is completely consistent with the result obtained by Laplace transform method, which also shows the correctness of this method.Figure 12
Transmissibility of system: (a) three-dimensional diagram of transmissibility; (b)α=0.5 transmissibility curve; (c) α=1 transmissibility curve; (d) α=1.5 transmissibility curve.
(a)(b)(c)(d)


## 4.1. Complex Modulus of Linear Model

The usual way to find the complex modulus is to perform Fourier transform or Laplace transform on the constitutive equation. When the fractional constitutive equation uses the two aforementioned transformations to find the complex modulus, the fractional power of the imaginary unit will be substituted, and the formula is more complicated when the real part and the imaginary part are separated. In this paper, the fractional derivative (6) of the sine (cosine) function is used to transform the linear model in the previous section, and the complex modulus of the model is obtained in a very convenient way. First, the fractional Kelvin model is taken as an example, assuming that the strain is a harmonic excitation ε=ε0sinωt, and ε0 is the excitation amplitude, and then equation (19) can be rewritten as follows:(22)σ=Eε0sinωt+ηωγε0sinωt+γπ2=E+ηωγcosγπ2ε0sinωt+ηωγsinγπ2ε0cosωt.Then, the storage modulus and dissipation modulus of the model are the coefficients in front ofε0sinωt and ε0cosωt, respectively,(23)Ys=E+ηωγcosγπ2Yl=ηωγsinγπ2.The loss angle is(24)δ=Arctanηωγsinγπ/2E+ηωγcosγπ/2.The above expressions are the same as the results obtained by Fourier transformation and separating the real and imaginary parts. In particular, whenγ=1, the resulting complex modulus and loss angle are(25)Ys=EYl=ηωδ=ArctanηωE.The above formula is completely consistent with the integer-order model. The fractional Maxwell model is slightly more complicated. Suppose that under harmonic excitationε=ε0sinωt, σ=σ0sinωt+φ is the stress response, and φ is the phase angle. Then, operate on the left and right sides of equation (18)(26)σ+ηEDγσ=σ01+ωγηEcosγπ2sinωt+φ+σ0ωγηEsinγπ2cosωt+φ=Aσ0sinωt+φ+BηDγε=ηωγε0sinωt+γπ2,where(27)A=1+ωγηEcosγπ22+ωγηEsinγπ22B=Arctanωγη/Esinγπ/21+ωγη/Ecosγπ/2.MoveA and B to the right side of the equation and rewrite it in the form of equation (22)(28)σ=ηωγAcosγπ2−Bε0sinωt+ηωγAsinγπ2−Bε0cosωt.The storage modulus, energy dissipation modulus, and loss angle of the fractional Maxwell model are further obtained as follows:(29)Ys=ηωγAcosγπ2−BYl=ηωγAsinγπ2−Bδ=γπ2−B.Whenγ=1,(30)A=1+ωηE2B=ArctanωηE.By substituting (30) into (29), the complex modulus and phase angle of the Maxwell model are expressed as follows:(31)Ys=Eη2ω2E2+η2ω2Yl=E2ηωE2+η2ω2δ=ArctanEωη.The obtained results are completely consistent with the results obtained by the integer-order Maxwell model by the Laplace transform method. In the same way, the expressions of storage modulus, dissipation modulus, and loss angle of the two other models can be obtained by the same method, which will not be repeated in this article.


## 4.2. Response Characteristics of System Based on the Linear Model

Analyzing the response characteristics of linear system under simple harmonic excitation by using the fractional derivative formula of sine (cosine) function in Section2.2 is very easy. The dynamic model shown in Figure 11 is studied [38], where rt represents the base excitation; c and η represent the damping coefficients of conventional linear damping and fractional derivative damping, respectively; α is the fractional order; and k is the stiffness coefficient of the linear spring. Using Newton’s second law, the equation of motion in Figure 11 is(32)mx¨+cx˙+μDαx+kx=cr˙+μDαr+kr.Figure 11
Dynamic model of the fractional-order linear system under basic excitation.When the basic excitation is sinusoidal excitationrt=Rsinωt, where R is the excitation amplitude and ω is the excitation frequency. Then, we assume that the response has the form(33)xt=Xsinωt+ϕ.Transform equation (32) according to equation (6)(34)−mω2+μωαcosαπ2+kXsinωt+φ+cω+μωαsinαπ2Xcosωt+φ=μωαcosαπ2+kRsinωt+cω+μωαsinαπ2Rcosωt.Further, transform this equation to obtain(35)A1Xsinωt+φ+B1=A2Rsinωt+B2,where(36)A1=−mω2+μωαcosαπ2+k2+cω+μωαsinαπ22B1=Arctancω+μωαsinαπ/2−mω2+μωαcosαπ/2+kA2=μωαcosαπ2+k2+cω+μωαsinαπ22B2=Arctancω+μωαsinαπ/2μωαcosαπ/2+k.MoveA1 and B1 on the left side of equation (35) to the right to obtain the response expression(37)xt=Xsinωt+φ=A2A1Rsinωt+B2−B1.Then, the magnitude and phase angle of the response are(38)X=A2A1Rφ=B2−B1.When the system parameters arem=5, k=10, c=0.5, and μ=1.5, transmissibility of the system is shown in Figure 12. Figure 12(a) is the three-dimensional transmissibility diagram; Figures 12(b)–12(d) are the transmissibility diagrams of α=0.5, α=1, and α=1.5, respectively. The three-dimensional diagram of the transmissibility in Figure 12(a) shows that when α is far from 1, the transmissibility of the system response will increase. In addition, in combination with Figures 12(b)–12(d), the peak point of the transmissibility is closer and closer to the position of ω=1 with the increase in α. Especially, when α=1, the obtained transmissibility curve is completely consistent with the result obtained by Laplace transform method, which also shows the correctness of this method.Figure 12
Transmissibility of system: (a) three-dimensional diagram of transmissibility; (b)α=0.5 transmissibility curve; (c) α=1 transmissibility curve; (d) α=1.5 transmissibility curve.
(a)(b)(c)(d)


## 5. Fractional Nonlinear Model

### 5.1. Harmonic Form of Nonlinear Model

The quadratic nonlinearity of integer-order differential equation means that the equation contains one or more terms ofx2, xx˙, and x˙2, which is also named here; that is, the quadratic fractional-order nonlinear differential equation means that the equation contains one or more terms similar to Dαx⋅Dβx. Here, the formula is derived based on simple harmonic excitation. The quadratic nonlinear term is expressed as follows:(39)Dβ¯ε⋅Dβ˜ε=ε02ωβ¯+β˜sinωt+β¯π2sinωt+β˜π2=−ε02ωβ¯+β˜2cos2ωt+β¯+β˜2π−cosβ¯−β˜2π=−ε02ωβ¯+β˜2cos2ωt+β¯+β˜2π+ε02ωβ¯+β˜2cosβ¯−β˜2π=A2sin2ωt+B2+C2,where(40)A2=−ε02ωβ¯+β˜2B2=β¯+β˜+12πC2=ε02ωβ¯+β˜2cosβ¯−β˜2π.From the above formula, any quadratic nonlinear term can be transformed into a second-order harmonic term and a constant term. The coefficient of the second-order harmonic term is related to the excitation amplitude, frequency, and fractional order, whereas the phase angle is only related to fractional order. The constant term is also related to the excitation amplitude, frequency, and fractional order. If the model contains multiple second-order terms, then all second-order harmonic terms can be combined into one term according to the linear model method in the previous section. The constant term refers to the translation along the longitudinal axis; that is, whenβ¯−β˜≠±1, a DC component will exist.The cubic nonlinear term can also follow the previous operation but without a constant term. The derivation process is presented as follows:(41)Dtγ1ε⋅Dtγ2ε⋅Dtγ3ε=ε03ωγ1sinωt+γ1π2ωγ2sinωt+γ2π2ωγ3sinωt+γ3π2=−ε03ωγ1+γ2+γ32cos2ωt+γ1+γ22π−cosγ1−γ22πsinωt+γ3π2=−ε03ωγ1+γ2+γ34sin3ωt+γ1+γ2+γ32π−ε03ωγ1+γ2+γ34sinωt+γ1+γ2−γ32π+ε03ωγ1+γ2+γ32cosγ1−γ22πsinωt+γ3π2=A1sinωt+B1+A3sin3ωt+B3,where(42)A1=ε03ωγ1+γ2+γ32cosγ1−γ22πcosγ32π−12cosγ1+γ2+γ32π2+cosγ1−γ22πsinγ32π−12sinγ1+γ2+γ32π21/2B1=Arctancosγ1−γ2/2πsinγ3/2π−1/2sinγ1+γ2+γ3/2πcosγ1−γ2/2πcosγ3/2ππ−1/2cosγ1+γ2+γ3/2πA3=−ε03ωγ1+γ2+γ34B3=γ1+γ2+γ32π.Equation (41) shows that any cubic nonlinear term can be expressed in the form of first-order and third-order harmonics. Similarly, the expressions of the fifth-order and the seventh-order nonlinear terms can also be obtained, and the highest order of the transformed harmonic term is equal to the nonlinear order. To avoid the influence of the DC component caused by the even-order nonlinear term, the odd-order nonlinear order is adopted in this paper. The transformation expressions of linear model and cubic, quintic, and seventh-order nonlinear models are given below:(43)σ=A1sinωt+B1σ=A1sinωt+B1+A3sin3ωt+B3σ=A1sinωt+B1+A3sin3ωt+B3+A5sin5ωt+B5σ=A1sinωt+B1+A3sin3ωt+B3+A5sin5ωt+B5+A7sin7ωt+B7.Notably, the coefficients and phase angles of the same order are not the same in each formula. For example, in the fifth-order harmonic form, the first-order term contains the first-order harmonic term transformed from the third-order and fifth-order nonlinear terms, and the third-order term includes the third-order harmonic term transformed from the fifth-order nonlinear term.


### 5.2. Integer-Order Form of Nonlinear Model

Under harmonic excitation, the fractional-order nonlinear model can be transformed into the form of integer-order nonlinearity. Taking the third formula in (43) as an example, we can make further derivation.(44)σ=A1sinωt+B1+A3sin3ωt+B3+A5sin5ωt+B5=A1cosB1sinωt+A1sinB1cosωt+A3cosB3sin3ωt+A3sinB3cos3ωt+A5cosB5sin5ωt+A5sinB5cos5ωt.sin3ωt, cos3ωt, sin5ωt, and cos5ωt in the above formula can be expanded as(45)sin3ωt=sin2ωtcosωt+cos2ωtsinωt=3sinωtcos2ωt−sin3ωtcos3ωt=cos2ωtcosωt−sin2ωtsinωt=cos3ωt−3sin2ωtcosωtsin5ωt=sin2ωtcos3ωt+cos2ωtsin3ωt=5sinωtcos4ωt−10sin3ωtcos2ωt+sin5ωtcos5ωt=cos2ωtcosωt−sin2ωtsinωt=cos5ωt−10sin2ωtcos3ωt+5cosωtsin4ωt.Substitute equation (45) into equation (44) to further organize(46)σ=k1ε+k2ε˙+k3ε3+k4ε2ε˙+k5εε˙2+k6ε˙3+k7ε5+k8ε4ε˙+k9ε3ε˙2+k10ε2ε˙3+k11εε˙4+k12ε˙5,where(47)k1=A1cosB1ε0k2=A1sinB1ωε0k3=−A3cosB3ε03k4=−3A3sinB3ωε03k5=3A3cosB3ω2ε03k6=A3sinB3ω3ε03k7=A5cosB5ε05k8=5A5sinB5ωε05k9=−10A5cosB5ω2ε05k10=−10A5sinB5ω3ε05k11=5A5cosB5ω4ε05k12=A5sinB5ω5ε05.


### 5.3. Test Verification

This paper takes the test in reference [39] as an example, as shown in Figure 13. In the figure, ut is the shear displacement and Ft is the shear load. When the excitation frequency is 4 Hz and the dynamic strain amplitude is 10%, 50%, and 100%, respectively, the test data are obtained. In this paper, the harmonic and integer-order models derived from the fractional-order model are used for fitting and compared with the ADF model in the reference. After the fifth-order nonlinear fractional order model is converted into the fifth-order harmonic form, the amplitude and phase angle parameters of each order are shown in Table 1. Table 2 shows the coefficients of the corresponding integer-order model.Figure 13
Specimen configuration.Table 1
Fitting parameters in harmonic form.
10% strain50% strain100% strainA113.1662A143.5214A181.0966B10.394584B10.340361B10.327254A3−0.527916A33.71492A312.2928B30.719419B3−1.0746B3−0.8738A5−0.159523A50.838552A53.52747B50.68787B50.09556B50.14204Table 2
Coefficients in integer-order form.
10% strain amplitude50% strain amplitude100% strain amplitudek11.242970.8313320.795736k20.02059470.01171420.0107477k34.24684e−41.47175e−58.78007e−6k44.4409e−53.24489e−61.25179e−6k5−2.01701e−6−6.99e−8−4.17003e−8k6−2.34353e−8−1.71237e−9−6.60587e−10k7−1.37844e−62.85226e−94.17191e−10k8−2.25371e−75.4382e−111.18381e−11k92.18227e−8−4.51553e−11−6.60473e−12k107.13591e−10−1.72189e−13−3.74829e−14k11−1.72743e−113.57437e−145.22811e−15k12−1.12972e−132.726e−175.93408e−18Substitute the fitting data in Table1 into equation (47) to obtain the coefficients of integer-order form, as shown in Table 2. Table 2 shows that the coefficient gradually decreases with the increase in the nonlinear order. Considering that the coefficient of the integer-order form is derived from the harmonic form, it is completely equivalent to the harmonic form. According to the test data in reference [39], the fitting curve is completely consistent with the fractional harmonic form curve in Figure 13, thereby indicating the correctness of the integer-order form.Since the coefficients of the integer-order form are derived from the harmonic form, they are completely equivalent to the harmonic form, so only one fitting curve of the fractional-order model is shown in Figure14. As can be seen from Figure 14, when the strain is small, the fitting effect of the fractional model and the ADF model is not very good. With the increase of the strain amplitude, the fitting curve of the fractional model is better than that of the ADF model. In terms of the form of the mathematical model, the fractional-order model is also simpler than the ADF model, and the harmonic form is very convenient for parameter identification. The disadvantage is that the fractional order nonlinear model cannot construct a physical model with springs and dashpots in series and parallel like the linear model.Figure 14
Test data, ADF model fitting curve, and fractional model fitting curve. The shear strain amplitudes are (a) 10%, (b) 50%, and (c) 100%, respectively.
(a)(b)(c)


## 5.1. Harmonic Form of Nonlinear Model

The quadratic nonlinearity of integer-order differential equation means that the equation contains one or more terms ofx2, xx˙, and x˙2, which is also named here; that is, the quadratic fractional-order nonlinear differential equation means that the equation contains one or more terms similar to Dαx⋅Dβx. Here, the formula is derived based on simple harmonic excitation. The quadratic nonlinear term is expressed as follows:(39)Dβ¯ε⋅Dβ˜ε=ε02ωβ¯+β˜sinωt+β¯π2sinωt+β˜π2=−ε02ωβ¯+β˜2cos2ωt+β¯+β˜2π−cosβ¯−β˜2π=−ε02ωβ¯+β˜2cos2ωt+β¯+β˜2π+ε02ωβ¯+β˜2cosβ¯−β˜2π=A2sin2ωt+B2+C2,where(40)A2=−ε02ωβ¯+β˜2B2=β¯+β˜+12πC2=ε02ωβ¯+β˜2cosβ¯−β˜2π.From the above formula, any quadratic nonlinear term can be transformed into a second-order harmonic term and a constant term. The coefficient of the second-order harmonic term is related to the excitation amplitude, frequency, and fractional order, whereas the phase angle is only related to fractional order. The constant term is also related to the excitation amplitude, frequency, and fractional order. If the model contains multiple second-order terms, then all second-order harmonic terms can be combined into one term according to the linear model method in the previous section. The constant term refers to the translation along the longitudinal axis; that is, whenβ¯−β˜≠±1, a DC component will exist.The cubic nonlinear term can also follow the previous operation but without a constant term. The derivation process is presented as follows:(41)Dtγ1ε⋅Dtγ2ε⋅Dtγ3ε=ε03ωγ1sinωt+γ1π2ωγ2sinωt+γ2π2ωγ3sinωt+γ3π2=−ε03ωγ1+γ2+γ32cos2ωt+γ1+γ22π−cosγ1−γ22πsinωt+γ3π2=−ε03ωγ1+γ2+γ34sin3ωt+γ1+γ2+γ32π−ε03ωγ1+γ2+γ34sinωt+γ1+γ2−γ32π+ε03ωγ1+γ2+γ32cosγ1−γ22πsinωt+γ3π2=A1sinωt+B1+A3sin3ωt+B3,where(42)A1=ε03ωγ1+γ2+γ32cosγ1−γ22πcosγ32π−12cosγ1+γ2+γ32π2+cosγ1−γ22πsinγ32π−12sinγ1+γ2+γ32π21/2B1=Arctancosγ1−γ2/2πsinγ3/2π−1/2sinγ1+γ2+γ3/2πcosγ1−γ2/2πcosγ3/2ππ−1/2cosγ1+γ2+γ3/2πA3=−ε03ωγ1+γ2+γ34B3=γ1+γ2+γ32π.Equation (41) shows that any cubic nonlinear term can be expressed in the form of first-order and third-order harmonics. Similarly, the expressions of the fifth-order and the seventh-order nonlinear terms can also be obtained, and the highest order of the transformed harmonic term is equal to the nonlinear order. To avoid the influence of the DC component caused by the even-order nonlinear term, the odd-order nonlinear order is adopted in this paper. The transformation expressions of linear model and cubic, quintic, and seventh-order nonlinear models are given below:(43)σ=A1sinωt+B1σ=A1sinωt+B1+A3sin3ωt+B3σ=A1sinωt+B1+A3sin3ωt+B3+A5sin5ωt+B5σ=A1sinωt+B1+A3sin3ωt+B3+A5sin5ωt+B5+A7sin7ωt+B7.Notably, the coefficients and phase angles of the same order are not the same in each formula. For example, in the fifth-order harmonic form, the first-order term contains the first-order harmonic term transformed from the third-order and fifth-order nonlinear terms, and the third-order term includes the third-order harmonic term transformed from the fifth-order nonlinear term.


## 5.2. Integer-Order Form of Nonlinear Model

Under harmonic excitation, the fractional-order nonlinear model can be transformed into the form of integer-order nonlinearity. Taking the third formula in (43) as an example, we can make further derivation.(44)σ=A1sinωt+B1+A3sin3ωt+B3+A5sin5ωt+B5=A1cosB1sinωt+A1sinB1cosωt+A3cosB3sin3ωt+A3sinB3cos3ωt+A5cosB5sin5ωt+A5sinB5cos5ωt.sin3ωt, cos3ωt, sin5ωt, and cos5ωt in the above formula can be expanded as(45)sin3ωt=sin2ωtcosωt+cos2ωtsinωt=3sinωtcos2ωt−sin3ωtcos3ωt=cos2ωtcosωt−sin2ωtsinωt=cos3ωt−3sin2ωtcosωtsin5ωt=sin2ωtcos3ωt+cos2ωtsin3ωt=5sinωtcos4ωt−10sin3ωtcos2ωt+sin5ωtcos5ωt=cos2ωtcosωt−sin2ωtsinωt=cos5ωt−10sin2ωtcos3ωt+5cosωtsin4ωt.Substitute equation (45) into equation (44) to further organize(46)σ=k1ε+k2ε˙+k3ε3+k4ε2ε˙+k5εε˙2+k6ε˙3+k7ε5+k8ε4ε˙+k9ε3ε˙2+k10ε2ε˙3+k11εε˙4+k12ε˙5,where(47)k1=A1cosB1ε0k2=A1sinB1ωε0k3=−A3cosB3ε03k4=−3A3sinB3ωε03k5=3A3cosB3ω2ε03k6=A3sinB3ω3ε03k7=A5cosB5ε05k8=5A5sinB5ωε05k9=−10A5cosB5ω2ε05k10=−10A5sinB5ω3ε05k11=5A5cosB5ω4ε05k12=A5sinB5ω5ε05.


## 5.3. Test Verification

This paper takes the test in reference [39] as an example, as shown in Figure 13. In the figure, ut is the shear displacement and Ft is the shear load. When the excitation frequency is 4 Hz and the dynamic strain amplitude is 10%, 50%, and 100%, respectively, the test data are obtained. In this paper, the harmonic and integer-order models derived from the fractional-order model are used for fitting and compared with the ADF model in the reference. After the fifth-order nonlinear fractional order model is converted into the fifth-order harmonic form, the amplitude and phase angle parameters of each order are shown in Table 1. Table 2 shows the coefficients of the corresponding integer-order model.Figure 13
Specimen configuration.Table 1
Fitting parameters in harmonic form.
10% strain50% strain100% strainA113.1662A143.5214A181.0966B10.394584B10.340361B10.327254A3−0.527916A33.71492A312.2928B30.719419B3−1.0746B3−0.8738A5−0.159523A50.838552A53.52747B50.68787B50.09556B50.14204Table 2
Coefficients in integer-order form.
10% strain amplitude50% strain amplitude100% strain amplitudek11.242970.8313320.795736k20.02059470.01171420.0107477k34.24684e−41.47175e−58.78007e−6k44.4409e−53.24489e−61.25179e−6k5−2.01701e−6−6.99e−8−4.17003e−8k6−2.34353e−8−1.71237e−9−6.60587e−10k7−1.37844e−62.85226e−94.17191e−10k8−2.25371e−75.4382e−111.18381e−11k92.18227e−8−4.51553e−11−6.60473e−12k107.13591e−10−1.72189e−13−3.74829e−14k11−1.72743e−113.57437e−145.22811e−15k12−1.12972e−132.726e−175.93408e−18Substitute the fitting data in Table1 into equation (47) to obtain the coefficients of integer-order form, as shown in Table 2. Table 2 shows that the coefficient gradually decreases with the increase in the nonlinear order. Considering that the coefficient of the integer-order form is derived from the harmonic form, it is completely equivalent to the harmonic form. According to the test data in reference [39], the fitting curve is completely consistent with the fractional harmonic form curve in Figure 13, thereby indicating the correctness of the integer-order form.Since the coefficients of the integer-order form are derived from the harmonic form, they are completely equivalent to the harmonic form, so only one fitting curve of the fractional-order model is shown in Figure14. As can be seen from Figure 14, when the strain is small, the fitting effect of the fractional model and the ADF model is not very good. With the increase of the strain amplitude, the fitting curve of the fractional model is better than that of the ADF model. In terms of the form of the mathematical model, the fractional-order model is also simpler than the ADF model, and the harmonic form is very convenient for parameter identification. The disadvantage is that the fractional order nonlinear model cannot construct a physical model with springs and dashpots in series and parallel like the linear model.Figure 14
Test data, ADF model fitting curve, and fractional model fitting curve. The shear strain amplitudes are (a) 10%, (b) 50%, and (c) 100%, respectively.
(a)(b)(c)


## 6. Conclusion

This paper studies the viscoelastic fractional derivative model based on harmonic excitation. First, a numerical example is used to illustrate the correctness of the fractional derivative formula obtained by extending the integer derivative formula of the sine (cosine) function. Then, a calculation method of the complex modulus of the fractional model under harmonic excitation, which is simpler than the frequency domain transform method, is given. Second, the response characteristics of the fractional-order linear system under basic excitation are analyzed by using the fractional derivative expansion formula of sinusoidal function. Finally, the nonlinear model was converted into harmonic form and integer-order form, and the test data in the reference were used as the original data for parameter fitting. Compared with the hysteresis elastic displacement field model, the results showed that the fractional-order nonlinear model has higher accuracy.



---

*Source: 4799387-2022-07-18.xml*