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.
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:where is the nearest integer, and is the binomial coefficient. If the selected calculation step is small enough, then the limit operation in (1) can be ignored.
The Riemann–Liouville definition is extended from integral order integral to noninteger order, and it is expressed aswhere 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 , its derivative expression under the Riemann–Liouville definition is expressed as follows:where 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
The fractional derivative expression of the sine (cosine) chord function is obtained by replacing its order n with any real number
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.


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 aswhere and are the normal stress and normal strain, respectively, and is the elastic modulus. Dashpot is the simplest viscous fluid model, which follows Newton’s law of viscosity: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:where and . 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:where and .


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 (Figure 5). It is also known as the standard linear solid model, and its constitutive equation is expressed asor normal formwhere , , and . The Maxwell model and Kelvin model are connected in series to form a four-element Burgers model (Figure 6). Its normative constitutive equation iswhere , , , and .


From (13), we can further generalize to the general form of the differential viscoelastic constitutive equation:where P and Q are differential operators.where and are material parameters. 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 iswhere is the viscosity coefficient, is the fractional derivative operator, and . Referring to the general form (14) of integral order constitutive equation, the general form of the constitutive equation expressed by the fractional derivative iswhere and are the material parameters, and both are greater than 0 and less than 1. When and are positive integers, equation (17) becomes the integral derivative constitutive equation.
The fractional Maxwell model is shown in Figure 7. 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.

Figure 8 shows the fractional Kelvin model, and the corresponding constitutive equation can also be obtained as

The fractional three-element model is shown in Figure 9, and its constitutive equation is

The fractional Burgers model contains two different Abel dashpots, and its constitutive model is Figure 10

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 , and is the excitation amplitude, and then equation (19) can be rewritten as follows:
Then, the storage modulus and dissipation modulus of the model are the coefficients in front of and , respectively,
The loss angle is
The above expressions are the same as the results obtained by Fourier transformation and separating the real and imaginary parts. In particular, when , the resulting complex modulus and loss angle are
The above formula is completely consistent with the integer-order model. The fractional Maxwell model is slightly more complicated. Suppose that under harmonic excitation , is the stress response, and is the phase angle. Then, operate on the left and right sides of equation (18)where
Move A and B to the right side of the equation and rewrite it in the form of equation (22)
The storage modulus, energy dissipation modulus, and loss angle of the fractional Maxwell model are further obtained as follows:
When ,
By substituting (30) into (29), the complex modulus and phase angle of the Maxwell model are expressed as follows:
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 Section 2.2 is very easy. The dynamic model shown in Figure 11 is studied [38], where represents the base excitation; and represent the damping coefficients of conventional linear damping and fractional derivative damping, respectively; is the fractional order; and is the stiffness coefficient of the linear spring. Using Newton’s second law, the equation of motion in Figure 11 is

When the basic excitation is sinusoidal excitation , where is the excitation amplitude and is the excitation frequency. Then, we assume that the response has the form
Transform equation (32) according to equation (6)
Further, transform this equation to obtainwhere
Move and on the left side of equation (35) to the right to obtain the response expression
Then, the magnitude and phase angle of the response are
When the system parameters are , , , and , 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 , , and , 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 with the increase in . Especially, when , the obtained transmissibility curve is completely consistent with the result obtained by Laplace transform method, which also shows the correctness of this method.

(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 of , , and , 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 . Here, the formula is derived based on simple harmonic excitation. The quadratic nonlinear term is expressed as follows:where
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 , 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:where
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:
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.
, , , and in the above formula can be expanded as
Substitute equation (45) into equation (44) to further organizewhere
5.3. Test Verification
This paper takes the test in reference [39] as an example, as shown in Figure 13. In the figure, is the shear displacement and 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.

Substitute the fitting data in Table 1 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 Figure 14. 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.

(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.
Data Availability
The data used to support the findings of this study are cited at relevant places within the text as references.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was funded by the National Natural Science Foundation of China (Grant no. 11472133).