Next Article in Journal
Some Local Fractional Hilbert-Type Inequalities
Next Article in Special Issue
Fractional Order Sequential Minimal Optimization Classification Method
Previous Article in Journal
Forecasting Cryptocurrency Prices Using LSTM, GRU, and Bi-Directional LSTM: A Deep Learning Approach
Previous Article in Special Issue
Inflation and Fractional Quantum Cosmology
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Solving and Numerical Simulations of Fractional-Order Governing Equation for Micro-Beams

1
College of Science, North China University of Science and Technology, Tangshan 063210, China
2
Hebei Key Laboratory of Data Science and Application, North China University of Science and Technology, Tangshan 063210, China
3
The Key Laboratory of Engineering Computing in Tangshan City, North China University of Science and Technology, Tangshan 063210, China
4
Tangshan Intelligent Industry and Image Processing Technology Innovation Center, North China University of Science and Technology, Tangshan 063210, China
5
Hebei Engineering Research Center for the Intelligentization of Iron Ore Optimization and Ironmaking Raw Materials Preparation Processes, North China University of Science and Technology, Tangshan 063210, China
6
School of Science, Yanshan University, Qinhuangdao 066004, China
7
LE STUDIUM, Loire Valley Institute for Advanced Studies Orléans & Tours, 45000 Orléans, France
*
Author to whom correspondence should be addressed.
Submission received: 31 December 2022 / Revised: 15 February 2023 / Accepted: 16 February 2023 / Published: 18 February 2023

Abstract

:
This paper applies a recently proposed numerical algorithm to discuss the deflection of viscoelastic micro-beams in the time domain with direct access. A nonlinear-fractional order model for viscoelastic micro-beams is constructed. Before solving the governing equations, the operator matrices of shifted Chebyshev polynomials are derived first. Shifted Chebyshev polynomials are used to approximate the deflection function, and the nonlinear fractional order governing equation is expressed in the form of operator matrices. Next, the collocation method is used to discretize the equations into the form of algebraic equations for solution. It simplifies the calculation. MATLAB software was used to program this algorithm to simulate the numerical solution of the deflection. The effectiveness and accuracy of the algorithm are verified by the numerical example. Finally, numerical simulations are carried out on the viscoelastic micro-beams. It is found that the viscous damping coefficient is inversely proportional to the deflection, and the length scale parameter of the micro-beam is also inversely proportional to the deflection. In addition, the stress and strain of micro-beam, the change of deflection under different simple harmonic loads, and potential energy of micro-beam are discussed. The results of the study fully demonstrated that the shifted Chebyshev polynomial algorithm is effective for the numerical simulations of viscoelastic micro-beams.

1. Introduction

Fractional calculus is the definition of arbitrary order differential and integral. It is unified with integer order calculus and is the extension of integer order calculus. Integer order calculus has been generally accepted as a mathematical tool to describe theories of classical physics and related subjects. However, when doing research on some complex systems and complex phenomena, the description of these systems by classical integer order calculus equations will encounter some problems. For example, it is necessary to construct nonlinear equations, introduce some empirical parameters and assumptions, construct new models, etc. Based on the above reasons, people urgently hope to have an available mathematical tool and basic principles to model these complex systems. Fractional differential equations are well suited to characterize materials and processes with memory and genetic properties [1,2,3,4,5,6]. It has the advantages of simple modeling, clear physical meaning of parameters, and accurate description for the description of complex systems. Therefore, fractional calculus is known as one of the important tools for mathematical modeling of complex mechanics and physical processes [7,8,9,10]. Nowadays, inspired by many physical phenomena, scholars worldwide have made more and more contributions in different fields by applying fractional-order constitutive relations and fractional-order systems containing fractional-order derivatives and integrals [11,12,13], and these new models are superior to the traditional integer-order ones. Fractal methods are widely used in the modeling of deformation–relaxation processes. Among them, viscoelastic materials are very closely coupled with fractional-order differential and fractional-order integral operators. Viscoelastic materials have four major characteristics, which are relaxation, creep, preconditioning, and hysteresis, and the traditional integer-order vibration equations cannot accurately describe the viscoelastic material properties. However, the fractional order derivative model has fewer parameters and a better fitting effect. Hendy et al. [14] proposed a unified new model for the thermal viscoelasticity theory of fractional heat transfer. Bonfanti et al. [15] proposed that fractional calculus is a promising tool for accurately describing the rheological properties of soft materials and provided an easy-to-understand fractional viscoelastic model description.
Researchers pay more and more attention to the combination of viscoelastic properties and fractional constitutive models. They have established that the fractional-order constitutive model can better describe the viscoelastic behavior, so a series of related studies have been conducted on it. Lewandowski [16], for example, used the fractional constitutive model of the beam to study the nonlinear vibration behavior of the beam under simple harmonic loading. Duan et al. [17] considered relaxation, creep, dissipation, and hysteresis resulting from a six-parameter fractional constitutive model and its particular cases. Usuki [18] investigated the effect of viscoelastic fractional order through drawing of the phase and group rate curves on beams with solid round cross-sectional sections. Catania [19] analyzed the vibration of fractional-order nonuniform beams using the finite element method. The rheological model used is a combination of the fractional Zener model and fractional standard linear solid model. Xu [20] discretized the beam according to the discontinuous characteristic position of the beam during vibration analysis, and then calculated the vibration extent of the beam from the spreading, transfer, and reflection matrices associated with the beam.
The viscoelastic control equations established using the fractional order constitutive model are not easily solved analytically [21]. Therefore, for ideal results, researchers mostly use the discussion of numerical solutions to make approximations. Approaches to the solution of fractional order governing equations also include the Galerkin method [22], the finite element method [23], the finite difference method [24], and the Laplace transform [25]. However, the limitation of these methods is that they cannot be solved directly in the time domain and require complex transformation processes. In other words, a Laplace transform of the frequency domain solution is required to obtain the time domain solution. The complexity of Laplace transform and inverse transform makes the fractional order equations difficult to answer efficiently [26]. The variational iterative method is an influential method for solving some linear and nonlinear problems. Birol et al. [27] studied the basic idea of the variational iterative method and its convergence, and used the variational iterative method to solve some fractional order differential equations. The Adomian decomposition method was used by Duan et al. [28] to solve fractional order differential equations. Muhammad Abbas et al. [29] used cubic B-spline functions to solve third-order time fractional differential equations. In addition, polynomial algorithms also work well when solving approximations to fractional order systems, such as Legendre wavelet method [30], Legendre polynomials method, Bernstein polynomials method [31], Chebyshev wavelet method [32], etc. The method approximates the unknown function with polynomials, uses the properties of fractional order differentiation to infer the matrices of differential operators of the polynomials, in addition to combining the idea of operator matrices to turn it into a problem of solving a system of the algebraic equations.
In recent years, microstructure-scale experiments confirm the existence of size effects on the mechanical properties of microstructures. Traditional mechanical theories do not contain any scale-associated parameters in their constitutive relations, and are not able to relate and explain the phenomenon of dimensional effects on mechanical properties. Therefore, the research on microstructures in the engineering field is gradually increasing. For example, in [33], the size-dependent buckling of variable-thickness FG-CNTRC microplates under inhomogeneous biaxial compression was investigated using third-order shear deformation theory (TSDT) and modified couple stress theory (MCST) combined with the total potential energy principle. In [34], a general PD beam-shell model based on microbeam bonding is proposed. The bond deformation is obtained directly by interpolation method, the micropotential energy of the bond is established, and the micromodulus of the beam-shell model is solved spontaneously.
Micro-beams are one of the most diffusible microstructures and have been analyzed by many researchers. The vibration amplitude of functionally graded (FGM) microbeams has been studied in detail by Al-Basyouni et al. [35]. Farokhi et al. [36] discussed the free vibration theory of crystalline microbeams and solved the equations of motion of PQC microbeams using modified coupled stress theory (MCST) and differential quadrature method (DQM). Zhang et al. [37] have numerically researched the modal properties of a microperforated sandwich beam having a corrugated square honeycomb hybrid core using the finite element method. It can be seen that research on microstructures is gradually increasing in the engineering field.
The chapter distribution is shown below. In this paper, Section 2 introduces the fractional order calculus definition and properties, and establishes the dimensionless fractional order constitutive equations for viscoelastic micro-beams. Section 3 presents the relevant properties of shifted Chebyshev polynomials and computes integer-order and fractional-order differential operator matrices. A number of examples of this algorithm are presented in Section 4. The stress, strain, deflection, and potential energy of viscoelastic micro-beams were obtained in Section 5. The practical application of the algorithm in the study of viscoelastic materials is demonstrated. Section 6 concludes the paper.

2. Pre-Requisite Knowledge

In this section, the fractional order calculus is defined and its properties are described in detail for easier understanding of the applications that follow [38]. Combining the theoretical knowledge of viscoelastic materials, a governing equation for a viscoelastic micro-beam is derived and established using a fractional order model.

2.1. Definition and Related Properties of Fractional Order Calculus

From the existing literature [21], it is clear that the fractional order calculus has several different ways of definition, including the Riemann–Liouville definition [39], the Caputo definition, and the Grünwald–Letnikov definition [40]. This paper applies the Caputo definition of a fractional order differential operator [39].
Definition 1. 
The α > 0 order derivative of Caputo may be expressed by the following equation:
c D x α f ( x ) = d q f ( x ) d x q α = q N + , 1 Γ ( q α ) 0 x f ( q ) ( τ ) ( x τ ) α + 1 q 0 q 1 < α < q .
where f ( x ) is a continuous differentiable function of order q on [ 0 , + ) . When 0 < α < 1 , c D x α f ( x ) is a Caputo fractional order derivative operator. Γ ( ) is Gamma function, Γ ( y ) = 0 + e x x y 1 d x .
From Definition 1, the fractional order derivative of a polynomial can be introduced as the following equation:
c D x α x p = Γ ( p + 1 ) Γ ( p + 1 α ) x p α , p = 1 , 2 , 0 , p = 0
The fractional order derivative of Caputo has the following properties.
The result of finding the fractional order derivative of Caputo for a constant is zero: c D x α C = 0 , where C is a constant.
The Caputo fractional order derivative satisfies the linearity property [41]: c D x α [ λ f ( x ) + μ g ( x ) ] = λ c D x α f ( x ) + μ c D x α g ( x ) , where λ , μ R .
Definition 2. 
Supposing that f ( x ) is a continuous and integrable function on ( 0 , + ) , 0 < β < 1 , we have [21]:
c D x β f ( x ) = 1 Γ ( β ) 0 x ( x τ ) β 1 f ( τ ) d τ
From Definition 2, the fractional order integration of a polynomial can be introduced as the following equation:
c D x β x p = Γ ( p + 1 ) Γ ( p + 1 + β ) x p + β , p = 0 , 1 , 2 ,
Definition 3. 
If f ( x ) and φ ( x ) along with all its derivatives are continuous in [ 0 , + ] , under this condition, the Leibniz rule for fractional differentiation takes the form [21]:
c D x α ( φ ( x ) f ( x ) ) = k = 0 α k φ ( k ) ( x ) c D x α k f ( x )

2.2. Establishment of Viscoelastic Micro-Beams Constitutive Equations

In this part, we will discuss the motion equations for viscoelastic micro-beams. The motion equations are obtained from the Euler–Bernoulli beam theory and the Modified Coupled Stress Theory. For the convenience of the derivation, we will ignore the effects of rotational inertia and shear stress.
The cross section of viscoelastic micro-beams is denoted by A, the density by ρ , and the modulus of elasticity by E. The viscoelastic micro-beam studied in this paper is a beam structure with both ends fixed, which is illustrated in Figure 1.
The kinetic (T) and potential ( U p ) energies for viscoelastic micro-beams can be obtained using Modified Couple Stress Theory (MCST) [22]:
T = 1 2 ρ A 0 L ω ( x , t ) t 2 d x U p = V ( σ : ε + m : γ )
where ω represents the transverse deflection, t represents the time, x represents the axial coordinate, V represents the occupied area volume, σ represents the stress tensor, ε represents the strain tensor, m represents the deviation part of the couple stress tensor, and γ represents the symmetric curvature tensor.
Through the von Kármán stress–strain relationship, the strain expression for the nonlinear Euler–Bernoulli beam can be written as:
ε = 1 2 ω ( x , t ) x 2 z 2 ω ( x , t ) x 2
where z represents the vertical coordinate.
Furthermore, the relationship between the symmetric curvature tensor γ and the rotation vector θ might usefully be described in [42]:
θ = 1 2 curl ( u ) γ i j = 1 2 θ j x i + θ i x j
where u is the displacement vector.
By means of Equation (8), the elements of γ are obtained as:
γ 12 = γ 21 = 1 2 2 ω ( x , t ) x 2
using Equation (9) to find the coupled stress tensor as m i j = 2 μ L 2 . γ i j , where the components are found to be:
m 12 = m 21 = μ L 2 2 ω ( x , t ) x 2
μ is the shear modulus of elasticity, and L is the length scale parameter.
From Equations (6), (7), (9), and (10), the potential energy for the viscoelastic micro-beam can be obtained by calculating [22] as:
U p = 1 2 V σ ε + m x y γ x y + m y x γ y x d V = 0 L 1 8 E A ω ( x , t ) x 4 + 1 2 E I 2 ω ( x , t ) x 2 2 + 1 2 μ A L 2 2 ω ( x , t ) x 2 2 d x
For a viscoelastic micro-beam affected by an external loads f, the amount of change in virtual work carried out with an applied load and the amount of change in virtual work carried out with non-conservative viscous damping loads, considering the viscous damping factor, are [22,43]:
δ W f = 0 L f ( x , t ) δ W d x δ W D = a 0 L ω t δ W d x
From Hamilton’s principle, it is obtained that [22,43]:
t 1 t 2 ( δ T δ U p + δ W ) d t = 0
The kinetic energy expressions of Equation (6) and Equations (11) and (12) are brought into Equation (13) to obtain the equations of motion governing for a nonlinear viscoelastic micro-beam as:
ρ A 2 ω ( x , t ) t 2 + E I + μ A L 2 4 ω ( x , t ) x 4 3 2 E A ω ( x , t ) x 2 2 ω ( x , t ) x 2 + a ω ( x , t ) t = f ( x , t )
In studying the properties of viscoelastic micro-beams, we will invoke the fractional Kelvin–Voigt model. For obtaining a fractional order constitutive equation, the moduli of elasticity E and shear modulus μ are modified to [22]:
E = E + E · η d α t α μ = μ + μ · η d α t α
where η d denotes the viscoelastic coefficient, and α denotes the order of the fractional order derivative ( 0 < α 1 , integer order at α = 1 ). Thus, Equation (14) can be rewritten as:
ρ A 2 ω ( x , t ) t 2 + E I + μ A L 2 4 ω ( x , t ) x 4 + E I + μ A L 2 η d C D t α 4 ω ( x , t ) x 4 3 2 E A ω ( x , t ) x 2 2 ω ( x , t ) x 2 3 2 E A η d C D t α ω ( x , t ) x 2 2 ω ( x , t ) x 2 + a ω ( x , t ) t = f ( x , t )
Next, the calculations are simplified by replacing the physical quantities involved in the above with suitable variables by means of dimensionless methods [22]:
ω ¯ = ω h x ¯ = x L t ¯ = E I / ρ A L 4 · t η 1 = μ A L 2 E I χ = η d E I ρ A L 4 α 2 a ¯ = a L 4 E I E I / ρ A L 4 η 2 = μ A L 2 E I η d E I ρ A L 4 α 2 f ¯ ( x , t ) = f ( x , t ) L 4 E I h β 1 = 3 2 A h 2 I β 2 = 3 2 A h 2 I η d E I ρ A L 4 α 2
where h is the height of the micro-beam. Replacing the dimensionless parameters into the equation, the following fractional order viscoelastic constitutive equation can be obtained:
2 ω ( x , t ) t 2 + 1 + η 1 4 ω ( x , t ) x 4 + χ + η 2 C D t α 4 ω ( x , t ) x 4 β 1 2 ω ( x , t ) x 2 ω ( x , t ) x 2 β 2 C D t α 2 ω ( x , t ) x 2 ω ( x , t ) x 2 + a ω ( x , t ) t = f ( x , t )
In the above equation, the overbar has been removed from the symbol for simplicity, where the boundary conditions of the equation are:
ω ( 0 , t ) = ω ( L , t ) = ω ( 0 , t ) x = ω ( L , t ) x = 0
and the initial conditions are:
ω ( x , 0 ) = ω ( x , 0 ) t = 0

3. Shift Chebyshev Polynomial Numerical Algorithm

In the literature reviewed, we can find that there are many ways to solve the governing equations of fractional viscoelastic micro-beams, such as using the finite difference method (FDM) and the finite element method (FEM) to solve at the same time. In addition, the FDM method is required to be discretized in the time domain, and the FEM method is required to discretize in space. Alternatively, Laplace transform and inverse transform can be used to solve the problem. These methods are more complicated, and it is not easy to directly solve the deflection changes of viscoelastic material micro-beams at different times and different positions in the time domain [44]. Using the algorithm proposed in this paper to solve the fractional differential equation only needs to apply the least squares method and the MATLAB program, and the program takes about 20 s to run, which is simple and fast, and has a good accuracy rate.
The shifted Chebyshev polynomial algorithm can be better applied to the solution of fractional differential equations. Therefore, the properties are described in detail in this section for the shifted Chebyshev polynomials. In addition, the viscoelastic micro-beam governing equation with operator matrix form is derived.

3.1. Correlation Properties of Shifted Chebyshev Polynomials

Chebyshev polynomials satisfy the recurrence relation [45,46,47]:
R i + 1 ( s ) = 2 s R i ( s ) R i 1 ( s ) , i = 1 , 2 ,
where R 0 ( s ) = 1 , R 1 ( s ) = s , s [ 1 , 1 ] , i = 1 , 2 , …
For the Shifted Chebyshev polynomials, the interval of the independent variable s becomes [ 0 , L ] , in which L is a non-negative real number. According to the shifted Chebyshev polynomials, a new independent variable u is introduced such that s = 2 u L 1 . Therefore, the recurrence relation of the new independent variable u with respect to the shifted Chebyshev polynomials is obtained as:
R ¯ i + 1 ( u ) = 2 2 u L 1 R ¯ i ( u ) R ¯ i 1 ( u ) , i = 1 , 2 ,
where R ¯ 0 ( u ) = 1 , R ¯ 1 ( u ) = 2 u L 1 .
On the interval [ 0 , L ] , the general formula for the shift Chebyshev based polynomials has been written as [48]:
R ¯ i ( u ) = i k = 0 i ( 1 ) i k ( i + k 1 ) ! 2 2 k ( i k ) ! ( 2 k ) ! L k u k , i = 1 , 2 ,
where R ¯ i ( 0 ) = ( 1 ) i , R ¯ i ( L ) = 1 . The shifted Chebyshev polynomials with weighted orthogonality are in the closed range [ 0 , L ] . The orthogonality relation is satisfied by:
0 L R ¯ j ( u ) R ¯ k ( u ) ω L ( u ) d u = h k ,
where ω L ( u ) = 1 L u u 2 , h k = b k 2 π , k = j 0 , k j , b 0 = 2 , b k = 1 , k 1 .
Thus, the vector consisting of the shifted Chebyshev polynomials is:
Φ n ( u ) = [ R ¯ 0 ( u ) , R ¯ 1 ( u ) , , R ¯ n ( u ) ] T
Let M n ( u ) = [ 1 , u , u 2 , , u n ] T , and let C n be the coefficient matrix of the shifted Chebyshev polynomials. Then,
Φ n ( u ) = C n M n ( u )
C n = c 00 0 0 c 10 c 11 0 c n 0 c n 1 c n n
where C n = c i j i , j = 0 n , c i j = 1 , i = j = 0 2 2 L c i 1 , j 1 c i 1 , j c i 2 , j , others 0 , i < j or i < 0 or j < 0 . Clearly, C n is a full rank invertible matrix.

3.2. Functional Approximation of ω ( x , t )

Let ω ( x ) be a differentiable and square integrable function on the closed interval [ 0 , L ] . Then, the approximate form of ω ( x ) with respect to the shifted Chebyshev polynomials can be expressed as:
ω ( x ) = i = 0 a i R ¯ i ( x ) , i N
where a i = 1 h i 0 L ω ( x ) R ¯ i ( x ) ω L ( x ) d x , a i is the shifted Chebyshev polynomials coefficient. ω ( x ) can be approximated by the first n + 1 terms as:
ω ( x ) ω n ( x ) = i = 0 n a i R ¯ i ( x ) = A T Φ n ( x )
where A = a i i = 0 n , Φ n ( x ) is the matrix consisting of shifted Chebyshev polynomials. The form of Φ n ( x ) is shown below:
Φ n ( x ) = [ R ¯ 0 ( x ) , R ¯ 1 ( x ) , , R ¯ n ( x ) ] T
Similarly, for any continuous two-dimensional function ω ( x , t ) L 2 ( [ 0 , L ] × [ 0 , T ] ) , ω ( x , t ) is a continuous integrable function on this interval. It can be written in the form of a truncated sequence as follows:
ω ( x , t ) = i = 0 j = 0 U i j R ¯ i ( x ) R ¯ j ( t ) Φ n T ( x ) U Φ n ( t )
where U = U i j i j = 0 n , Φ n ( t ) is to replace the independent variable x in Φ n ( x ) with t.

3.3. Differential Operator Matrices on the Basis of Shifted Chebyshev Polynomials

3.3.1. Integer-Order Differential Operator Matrices

Definition 4. 
There exists a matrix Q x ( 1 ) satisfying condition d d x Φ n ( x ) = Q x ( 1 ) Φ n ( x ) . Then, Q x ( 1 ) is said to be a matrix of a first order differential operator based on shifted Chebyshev polynomials.
From Equation (26), the derivative d d x Φ n ( x ) of the first-order differential operator of Φ n ( x ) can be calculated as:
d d x Φ n ( x ) = d d x C n M n ( x ) = C n d d x M n ( x ) = C n 1 x x 2 x n 1 x n = C n 0 1 2 x ( n 1 ) x n 2 n x n 1 = C n H M n ( x )
where H is a square matrix of order n + 1 . H = 0 0 0 0 0 1 0 0 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 n 0 that is H = h i j i , j = 0 n , h i j = i , i = j + 1 0 , others .
From Φ n ( x ) = C n M n ( x ) , we have:
M n ( x ) = C n 1 Φ n ( x )
Therefore,
d d x Φ n ( x ) = C n H C n 1 Φ n ( x ) = Q x ( 1 ) Φ n ( x )
That is, Q x ( 1 ) = C n H C n 1 . Q x ( 1 ) is the matrix of a first-order differential operator based on shifted Chebyshev polynomials of Φ n ( x ) with respect to x.
Similarly, we can proceed to define the second order differential operator matrices.
Definition 5. 
There exists a matrix Q x ( 2 ) satisfying the condition d 2 d x 2 Φ n ( x ) = Q x ( 2 ) Φ n ( x ) . Then, Q x ( 2 ) is said to be a matrix of a second order differential operator based on shifted Chebyshev polynomials:
d 2 d x 2 Φ n ( x ) = d d x d d x Φ n ( x ) = d d x C n H C n 1 Φ n ( x ) = C n H C n 1 d d x Φ n ( x ) = C n H C n 1 2 Φ n ( x ) = Q x ( 2 ) Φ n ( x )
where Q x ( 2 ) = C n H C n 1 2 is the matrix of a second-order differential operator based on shifted Chebyshev polynomials of Φ n ( x ) with respect to x.
Thus, the m-order differential operator matrix Q x m based on the shifted Chebyshev polynomials can be derived:
d m d x m Φ n ( x ) = d m 1 d x m 1 d d x Φ n ( x ) = d m 1 d x m 1 C n H C n 1 Φ n ( x ) = C n H C n 1 d m 1 d x m 1 Φ n ( x ) = = C n H C n 1 m 1 d d x Φ n ( x ) = C n H C n 1 m Φ n ( x ) = Q x ( m ) Φ n ( x )
That is, Q x ( m ) = C n H C n 1 m = ( Q x ( 1 ) ) m , m N + .
In summary, it is obtained that:
ω ( x , t ) x Φ n T ( x ) U Φ n ( t ) x = Φ n T ( x ) C n H C n 1 T U Φ n ( t ) = Φ n T ( x ) Q x ( 1 ) T U Φ n ( t )
2 ω ( x , t ) x 2 2 Φ n T ( x ) U Φ n ( t ) x 2 = Φ n T ( x ) C n H C n 1 T 2 U Φ n ( t ) = Φ n T ( x ) Q x ( 1 ) T 2 U Φ n ( t )
4 ω ( x , t ) x 4 4 Φ n T ( x ) U Φ n ( t ) x 4 = Φ n T ( x ) C n H C n 1 T 4 U Φ n ( t ) = Φ n T ( x ) Q x ( 1 ) T 4 U Φ n ( t )
ω ( x , t ) t Φ n T ( x ) U Φ n ( t ) t = Φ n T ( x ) U C n H C n 1 Φ n ( t ) = Φ n T ( x ) U Q t ( 1 ) Φ n ( t )
2 ω ( x , t ) t 2 2 Φ n T ( x ) U Φ n ( t ) t 2 = Φ n T ( x ) U C n H C n 1 2 Φ n ( t ) = Φ n T ( x ) U Q t ( 1 ) 2 Φ n ( t )
ω ( x , t ) x 2 Φ n T ( x ) Q x ( 1 ) T U Φ n ( t ) Φ n T ( x ) Q x ( 1 ) T U Φ n ( t ) T = Φ n T ( x ) Q x ( 1 ) T U Φ n ( t ) Φ n T ( t ) U T Q x ( 1 ) Φ n ( x ) = Φ n T ( x ) Q x ( 1 ) T U C n M n ( t ) M n T ( t ) C n T U T Q x ( 1 ) Φ n ( x ) = Φ n T ( x ) Q x ( 1 ) T U C n 1 t t n t t 2 t n + 1 t n t n + 1 t 2 n C n T U T Q x ( 1 ) Φ n ( x )
2 ω ( x , t ) x 2 ω ( x , t ) x 2 Φ n T ( x ) Q x ( 1 ) T 2 U Φ n ( t ) · Φ n T ( x ) Q x ( 1 ) T U C n 1 t t n t t 2 t n + 1 t n t n + 1 t 2 n C n T U T Q x ( 1 ) Φ n ( x )
where Q t ( 1 ) is the matrix of a first-order differential operator based on shifted Chebyshev polynomials for Φ n ( t ) with respect to t.

3.3.2. Fractional Order Differential Operator Matrices

Definition 6. 
There exists a matrix Q t ( α ) ( 0 < α < 1 ) satisfying condition C D t α Φ n ( t ) = Q t ( α ) Φ n ( t ) ; then, Q t ( α ) is said to be a matrix of a differential operator of order α based on the shifted Chebyshev polynomials.
From Equation (26), the derivative C D t α Φ n ( x ) of the α -order differential operator of Φ n ( x ) can be calculated as:
C D t α Φ n ( t ) = C D t α C n M n ( t ) = C n C D t α M n ( t ) = C n 0 Γ ( 2 ) Γ ( 2 α ) t 1 α Γ ( n + 1 ) Γ ( n + 1 α ) t n α = C n H α M n ( t ) = C n H α C n 1 Φ n ( t )
where H α = 0 0 0 0 Γ ( 2 ) Γ ( 2 α ) t α 0 0 0 Γ ( n + 1 ) Γ ( n + 1 α ) t α , that is, H α = h α , i j i , j = 0 n , h α , i j = Γ ( i + 1 ) Γ ( i + 1 α ) t α , i = j 0 0 , others
Thus, Q t ( α ) = C n H α C n 1 , Q t ( α ) is the matrix of fractional order differential operator based on shifted Chebyshev polynomials for Φ n ( t ) with respect to t.
In summary, it is obtained that:
C D t α 4 ω ( x , t ) x 4 C D t α Φ n T ( x ) C n H C n 1 T 4 U Φ n ( t ) = Φ n T ( x ) C n H C n 1 T 4 U C D t α Φ n ( t ) = Φ n T ( x ) C n H C n 1 T 4 U C n H α C n 1 Φ n ( t ) = Φ n T ( x ) Q x ( 1 ) T 4 U Q t ( α ) Φ n ( t )

3.3.3. Handling of Nonlinear Terms

According to Definition 2 and Definition 3, the nonlinear terms of the governing Equation (18) can be dealt with as follows:
c D t α 2 ω ( x , t ) x 2 ω ( x , t ) x 2 = k = 0 + α k 2 ω ( x , t ) x 2 ( k ) · c D t α k ω ( x , t ) x 2 k = 0 n α k 2 ω ( x , t ) x 2 ( k ) · c D t α k ω ( x , t ) x 2 k = 0 n α k Φ n T ( x ) Q x ( 1 ) T 2 U Φ n ( t ) ( k ) · c D t α k Φ n T ( x ) Q x ( 1 ) T U C n 1 t t n t t 2 t n + 1 t n t n + 1 t 2 n C n T U T Q x ( 1 ) Φ n ( x ) = k = 0 n α k Φ n T ( x ) Q x ( 1 ) T 2 U Q t ( 1 ) k Φ n ( t ) · Φ n T ( x ) Q x ( 1 ) T U C n W C n T U T Q x ( 1 ) Φ n ( x )
where W = w i j i , j = 1 n + 1 , w i j = 0 α k > 0 , i = j = 1 Γ i + j 2 + 1 Γ i + j 2 + 1 α k x i + j 2 α k α k > 0 , others Γ i + j 2 + 1 Γ i + j 2 + 1 + k α x i + j 2 + k α α k < 0
Note in Formula (46) that, because d k M n ( t ) d t k = 0 ( n + 1 ) × 1 for k n + 1 , the infinite sum can be cut; just take the first n terms.

3.4. Algebraic Equation Form of the Viscoelastic Micro-Beam

Using Equations (37)–(46), Equation (18) is discretized into the following form:
Φ n T ( x ) U Q t ( 1 ) 2 Φ n ( t ) + 1 + η 1 Φ n T ( x ) ( Q x ( 1 ) ) T 4 U Φ n ( t ) + χ + η 2 Φ n T ( x ) ( Q x ( 1 ) ) T 4 U Q t ( α ) Φ n ( t ) β 1 Φ n T ( x ) Q x ( 1 ) T 2 U Φ n ( t ) Φ n T ( x ) Q x ( 1 ) T U Φ n ( t ) 2 β 2 k = 0 n α k Φ n T ( x ) Q x ( 1 ) T 2 U Q t ( 1 ) k Φ n ( t ) · Φ n T ( x ) Q x ( 1 ) T U C n W C n T U T Q x ( 1 ) Φ n ( x ) + a Φ n T ( x ) U Q t ( 1 ) Φ n ( t ) = f ( x , t )
The boundary conditions are: ω ( 0 , t ) = ω ( L , t ) = ω ( 0 , t ) x = ω ( L , t ) x = 0 .
The initial conditions are: ω ( x , 0 ) = ω ( x , 0 ) t = 0 .
The algorithm is based on the collocation method, which converts the fractional differential equation into a set of algebraic equations, solves the coefficient matrix U, uses the least squares method and MATLAB software to calculate, and substitutes the result into Equation (31) to obtain the numerical solution of the fractional differential equation. Through the introduction of the shifted Chebyshev polynomial algorithm in Section 3, we can summarize the steps required to use the algorithm as follows:
  • First, approximate the function: ω ( x , t ) Φ n T ( x ) U Φ n ( t ) ;
  • Integer-order fractional-order differential operators are derived;
  • Substituting the operator matrix into Equation (18) converts the initial equation into an algebraic equation;
  • Next, it is necessary to use the collocation method to discretize the variables; take the nodes to discretize the variables ( x , t ) into ( x i , y j )
    x i = 2 i 1 2 ( n + 1 ) L , i = 0 , 1 , 2 , , n , t j = 2 j 1 2 ( n + 1 ) T , j = 0 , 1 , 2 , , n ;
  • Finally, the system of algebraic equations is solved using the method of least squares.

4. Convergence Analysis and Numerical Example

4.1. Convergence Analysis

This section will introduce the uniform convergence of shifted Chebyshev polynomials.
Theorem 1. 
Assume that the function ω ( x ) : [ 0 , 1 ] R is m + 1 times continuously differentiable, i.e., ω ( x ) C m + 1 [ 0 , 1 ] and Y = Span B 0 , n , B 1 , n , B 2 , n , B n , n . If A T Φ n ( x ) is the best approximation to ω ( x ) from Y, the mean error bound can be expressed as:
ω ( x ) A T Φ n ( x ) 2 2 M S 2 m + 3 2 ( m + 1 ) ! 2 m + 3
where M = max x [ 0 , 1 ] ω m + 1 ( x ) , S = max 1 x 0 , x 0 , 0 x 0 1
Proof. 
Let the Taylor expansion of ω ( x ) be ω 1 ( x ) ; then,
ω 1 ( x ) = ω x 0 + ω x 0 x x 0 + ω x 0 x x 0 2 2 + + ω ( m ) x 0 x x 0 m m !
In addition, because ω ( x ) ω 1 ( x ) = ω ( m + 1 ) ( ε ) x x 0 m + 1 ( m + 1 ) ! , ε ( 0 , 1 ) , A T Φ n ( x ) is the best approximation to ω ( x ) from Y. Therefore,
ω ( x ) A T Φ n ( x ) 2 2 ω ( x ) ω 1 ( x ) 2 2 = 0 1 ω ( x ) ω 1 ( x ) 2 d x = 0 1 ω m + 1 ( ε ) x x 0 m + 1 ( m + 1 ) ! 2 d x M 2 [ ( m + 1 ) ! ] 2 0 1 x x 0 2 m + 2 d x 2 M 2 S 2 m + 3 [ ( m + 1 ) ! ] 2 ( 2 m + 3 )
The upper bound on ω ( x ) A T Φ n ( x ) 2 can be obtained by taking the square root. □

4.2. Numerical Example

The accuracy of the shifted Chebyshev polynomial algorithm will be illustrated by a numerical example. The following equation is a dimensionless equation, and the coefficients in the equation can be any value in the interval range without practical significance. The specific equation is given in below (See Appendix A for detailed calculation procedure):
10,000 2 ω ( x , t ) t 2 + 0.01 4 ω ( x , t ) x 4 + 0 . 01 C D t α 4 ω ( x , t ) x 4 0.1 2 ω ( x , t ) x 2 ω ( x , t ) x 2 0 . 01 C D t α 2 ω ( x , t ) x 2 ω ( x , t ) x 2 + 0.001 ω ( x , t ) t = f ( x , t )
The boundary conditions and initial conditions are written as follows:
ω ( 0 , t ) = ω ( 1 , t ) = ω ( 0 , t ) x = ω ( 1 , t ) x = 0 . ω ( x , 0 ) = ω ( x , 0 ) t = 0 . Where α = 0.5 , x [ 0 , 1 ] , t [ 0 , 1 ] .
f ( x , t ) = 20,000 x 2 ( 1 x ) 2 + 0.24 t 2 + 0.24 Γ ( 3 ) Γ ( 3 α ) t 2 α 0.1 2 12 x + 12 x 2 t 6 2 x 6 x 2 + 4 x 3 2 0.01 2 12 x + 12 x 2 2 x 6 x 2 + 4 x 3 2 Γ ( 7 ) Γ ( 7 α ) t 6 α + 0.002 t x 2 2 x 3 + x 4
The analytical solution of Equation (49) can be found as:
ω ( x , t ) = x 2 ( 1 x ) 2 t 2
When n = 4 , Equation (49) is solved by the shifted Chebyshev polynomial algorithm. The numerical solution is found by using a MATLAB program and matching point method. ω n ( x , t ) is the numerical solution, and e ( x , t ) is the absolute error:
e ( x , t ) = | ω ( x , t ) ω n ( x , t ) |
In Figure 2a–c are the analytical solution, numerical solution, and absolute error of Equation (49), respectively. By observing Figure 2c, it can be found that the change trend of the absolute error is: viewed from the x-axis, it is symmetrically distributed. From the perspective of the t-axis, the absolute error fluctuates less and is relatively stable in the range of t [ 0 , 0.9 ] and increases sharply at t [ 0.9 , 1 ] . However, from a macro point of view, the absolute error is small. It can be seen in Figure 2c that the maximum magnitude of error is 10 8 , and the minimum magnitude is 10 9 . Therefore, the shifted Chebyshev polynomial algorithm has a relatively high accuracy and can effectively solve such fractional order differential equations. This is then used to solve the problem for the governing equations of viscoelastic micro-beams.
Table 1, Table 2, Table 3 and Table 4 list the absolute errors between the numerical solution and the analytical solution when the derivative α = 0.2 , α = 0.5 , α = 0.8 , and α = 1 , respectively. When α = 0.2 , the series of absolute errors are 10 7 and 10 8 ; when α = 0.5 , the series of absolute errors are 10 9 and 10 10 , and when α = 0.8 , the series of absolute errors are 10 8 and 10 9 . We can see that, when calculating the fractional derivative equation, the fitting effect of the numerical solution is better. By observing Table 4, when α = 1 (first order derivative), the absolute error is relatively large, and the fitting effect is not ideal.
By comparing the absolute errors of solutions of fractional differential equations and integer order differential equations, for equations of this type, the numerical solutions obtained by fractional order differential equations have a better fitting effect than integer order differential equations.

5. Numerical Simulation

This section presents numerical simulations of the nonlinear viscoelastic micro-beam governing equation. All of the following studies will be based on the parameter values provided in Table 5. The MATLAB program for solving the fractional order differential equation of the same form as Equation (49) has been constructed in Section 4. Now, it is necessary to substitute the material parameters listed in Table 5 into the MATLAB program to find the numerical solution of the deflection first, i.e., the solution of Equation (49). After that, they are solved separately according to the relationship between deflection and stress, strain, and potential energy. The program only needs to run for about 20 s, which is fast and easy to solve. In addition, the solved deflection solution is fully consistent with the actual situation of the micro-beam, where the deflection occurs at the boundary conditions, and initial conditions are zero, and the most obvious deflection occurs at the middle of the micro-beam.

5.1. Effect of Viscous Damping Coefficient on Deflection of the Micro-Beam

Figure 3 shows the effect of different viscous damping coefficients on the deflection of the micro-beam. Because the effect of the viscous damping coefficient on the deflection is not very obvious intuitively, we apply a simple harmonic load ( f = 100 ( x 2 ) c o s 1.7 t ) with a frequency of Ω = 1.7 to reach the resonant frequency of the micro-beam (through multiple experiments, it is determined that the resonance frequency is between 1.69–1.71), so that the amplitude reaches the maximum value. In this way, the observation effect will be more intuitive and clear. The deflection change of the micro-beam is symmetrical from the middle position of the micro-beam. In other words, under the action of simple harmonic load, the deflection is the most obvious in the middle of the micro-beam, and the closer to the fixed points at both ends, the smaller the deflection. By observing Figure 3, it can be found that: when the viscous damping coefficient a = 0.1 , the maximum deflection is 5.99 × 10 4 cm ; when the viscous damping coefficient a = 0.5 , the maximum deflection is 5.61 × 10 4 cm ; when the viscous damping coefficient a = 0.1 , the maximum deflection is 5.39 × 10 4 cm . It can be found that the smaller the viscous damping coefficient, the greater the deflection of the micro-beam. This is consistent with the properties of viscoelastic damping. The larger the damping coefficient, the stronger the ability to resist the deformation of the material, resulting in smaller deflection. This is also consistent with the conclusion obtained in [49], that is, as the viscous damping coefficient increases, the amplitude of the longitudinal impedance decreases.

5.2. Effect of Length Scale Parameters on Deflection of the Micro-Beam

In this section, different micro-beams with length scale parameters are analyzed as research objects, and the applied load is f = 100 ( x 2 ) c o s t . Figure 4 shows the deflection variation of the micro-beam with five length scale parameters in the middle position. By observing Figure 4, it can be found that, as time goes on, the deflection increases. The micro-beam with the smallest length scale parameter ( l / h = 96 ) has the largest deflection, and the deflection reaches 5.8 × 10 5 cm at t = 1 s . When the micro-beam length scale parameter is 192, the deflection is the smallest. In addition, the deflection is 2.6 × 10 5 cm at t = 1 s . Under the same load, the deflection of the micro-beam decreases with the increase of the length scale parameters at the same time. It can be predicted that the stiffness of the structure increases as the length scale parameter increases. This is compared with the deflection of micro-beams under different length scale parameters shown in [50], and the conclusion is found to be consistent.

5.3. Effect of Different Simple Harmonic Loads on the Deflection of the Micro-Beam

The effect of different simple harmonic loads on the deflection of the micro-beam is investigated with a viscous damping coefficient of a = 0 . Since the use of simple harmonic forces is more frequent in practical engineering, we therefore mainly study the deflection variation of micro-beam when B ( x ) in f ( x , t ) = B ( x ) c o s Ω t is constant, primary, secondary, and tertiary functions, respectively, and Ω = 0.75 , 1 , 1.25 .
The four plots in Figure 5 show the deflection on the micro-beam at the moment of time t = 0.5 s for the four cases where B ( x ) in the simple harmonic load f ( x , t ) is constant (Figure 5a), primary function (Figure 5b), secondary function (Figure 5c), and tertiary function (Figure 5d), respectively. From the four graphs, it can be found that the deflection of the micro-beam increases as angular frequency increases, when B ( x ) takes the same value. We can guess that Ω = 0.75–1.25 is a stage that gradually tends to the resonance frequency. The rule at this stage is that the deflection increases with frequency. This is consistent with the reality. Moreover, this is also in line with the deflection curve of the Euler–Bernoulli beam [51,52], which shows that the analysis of the micro-beam deflection in this article is correct. Under the load in Figure 5d, at t = 0.5 s, the deflection produced by Ω = 1 and Ω = 1.25 is very close. Under the load in Figure 5c, the deflection produced when Ω = 0.75 and Ω = 1 is quite different. There is a maximum difference of 3 × 10 6 cm in the middle of the beam. This provides a new idea for the future research on simple harmonic loads, which is beneficial to expanding the idea of solving practical engineering problems.

5.4. Stress and Strain of the Micro-Beam

From the stress–strain relationship equation:
σ = E ( z ) ε + η α ε t α
Stress and strain are positively related. Figure 6a,b show the variation of strain and stress of the micro-beam at different locations at different times under the simple harmonic load f ( x , t ) = 100 ( x 2 ) c o s t . Through the three-dimensional graph, it can be clearly seen that the trend of stress and strain is consistent, showing a positive relationship. Under the simple harmonic load, the stress is symmetrical at x = 2.5 × 10 4 m in the middle of the micro-beam. At t = 0 s , the stress and strain of the micro-beam are zero at any position, which conforms to the initial condition of the governing equation of the viscoelastic micro-beam. As the loading time increases, the force resisting the deformation of the micro-beam also increases. The stress and strain values of the micro-beam are the smallest at the middle position, and the stress and strain gradually increase when approaching the two ends. The stress reaches the maximum at x = 0.625 × 10 4 m and x = 4.375 × 10 4 m , which are 29.01 Pa and 28.78 Pa , respectively. This shows that the micro-beam has the weakest ability to resist deformation in the middle position, and the closer to the end point, the stronger the ability to resist deformation. This is consistent with the theory of Euler–Bernoulli beams with both ends fixed. It is shown that the shifted Chebyshev polynomial algorithm is capable of solving the stress and strain problems in viscoelastic materials.

5.5. Potential Energy Change of the Micro-Beam

The previous sections compare the physical quantities of viscous damping coefficient, length scale parameters, micro-beam deflection, and strain with the conclusions in other papers. It can be verified that the algorithm proposed in this paper is effective for the numerical simulation of viscoelastic micro-beams. Therefore, in this section, we will numerically simulate the potential energy of micro-beams to provide more reference data for the study of viscoelastic materials.
The variation of potential energy for the viscoelastic micro-beam subjected to the simple harmonic load f = 100 ( x 2 ) c o s t will be studied in this section. The potential energy here is the energy generated by the elastic deformation of the micro-beam. The relationship between potential energy and deflection can be found from Equation (11). According to the algorithm proposed in this paper, the deflection occurred by the micro-beam is found. Then, the potential energy generated by the micro-beam during the motion is found. The three-dimensional relationship diagram can be obtained as shown in Figure 7.
Figure 7 shows the potential energy values generated by the micro-beam at different times and at different positions under the simple harmonic load f = 100 ( x 2 ) c o s t . By looking at the three-dimensional diagram on the left and the planar diagram on the right, we can see that the generated potential energy is symmetrical with respect to the micro-beam x = 2.5 × 10 4 m . The generated potential energy decreases and then increases from the middle of the micro-beam towards the two end points. The maximum potential energy can reach 25.68 J . The potential energy in the middle position x = 2.5 × 10 4 m is 7.13 J . The position with the smallest potential energy is 0.21 J at x = 1 × 10 4 m and x = 4 × 10 4 m . The potential energy starts to increase sharply about 6.25 × 10 5 m from the ends of the micro-beam. The potential energy generated increases with the time of application of the simple harmonic load.
See Table A1 for explanations of professional terms.

6. Conclusions

In this paper, a numerical algorithm for solving nonlinear-fractional order differential equations is proposed. The features of the algorithm are mainly reflected in the following three aspects. First, the algorithm can directly find the numerical solution of the unknown function in the time domain, and only needs to substitute the operator matrix form of the fractional order differential equation into the Matlab program, which is easier to operate. Secondly, the numerical solution of the fractional order differential equation is solved with high accuracy using this algorithm. This can be verified in the numerical example, where the error between the analytical and numerical solutions is of the order of 10 9 . In addition, the numerical simulations of micro-beams in viscoelastic materials is investigated in this paper, and the following conclusions are obtained. This has important practical implications for the study of viscoelastic materials:
  • The constitutive equations of the nonlinear-fractional-order viscoelastic micro-beam are first established and dimensionless. The numerical examples are used to demonstrate the effectiveness and accuracy of the algorithm;
  • Through numerical analysis, it is found that viscous damping has the effect of resisting deformation. The larger the viscous damping coefficient, the stronger the ability to resist deformation;
  • The algorithm is used to compare the deflection of micro-beams under different length scale parameters. It can be found that the deflection of the microbeam is inversely proportional to the length scale parameter;
  • The deflection changes of the viscoelastic micro-beam were obtained using this algorithm. The effect of different simple harmonic loads on the deflection variation of the micro-beam was studied. It can be found that, when the frequency of the simple harmonic load approaches the first resonance region, the deflection will increase as the frequency increases. It is also consistent with the actual situation;
  • Using this algorithm, we calculated the stress and strain of the micro-beam. The calculated stress and strain conform to the properties of viscoelastic materials, and the two are also proportional. The reliability of the algorithm is verified;
  • Using this algorithm, the change in potential energy of a viscoelastic micro-beam under the simple harmonic load was calculated. It is found that the value of the potential energy of the micro-beam is symmetrically distributed in the middle of the micro-beam, with the potential energy decreasing and then increasing. The potential energy increases sharply at the two end points of the micro-beam.
This paper provides a research direction for the research of fractional viscoelastic constitutive model and dynamic analysis of viscoelastic materials. In the future, further research can be carried out on this basis. The research directions are as follows:
  • The parameters in the paper are all selected from the references, and experiments can be carried out in the future to use the experimental data for kinetic analysis;
  • In the future, the algorithm in this paper can be used to compare the fitting effects of different fractional models to the same viscoelastic material micro-beam;
  • In the future, more complex viscoelastic materials can be numerically simulated using variable fractional orders.

Author Contributions

Writing—original draft, A.Y., Q.Z., J.Q., Y.C. (Yuhuan Cui) and Y.C. (Yiming Chen). All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the National Natural Science Foundation of China (52074126) and the Hebei Provincial Natural Science Foundation of China (E2022209110).

Data Availability Statement

No data was used for the research described in the article.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

This part explains the rationality of Equation (51) as a solution of Equation (49).
From Formula (51), we know that:
ω ( x , t ) t = 2 x 2 ( 1 x ) 2 t
2 ω ( x , t ) t 2 = 2 x 2 ( 1 x ) 2
ω ( x , t ) x = 2 ( x 3 x 2 + 2 x 3 ) t 2
2 ω ( x , t ) x 2 = 2 ( 1 6 x + 6 x 2 ) t 2
3 ω ( x , t ) x 3 = 12 ( 2 x 1 ) t 2
C D t α 4 ω ( x , t ) x 4 = C D t α 24 t 2 = 24 Γ ( 3 ) Γ ( 3 α ) t 2 α
C D t α 2 ω ( x , t ) x 2 ω ( x , t ) x 2 = C D t α 2 12 x + 12 x 2 2 x 6 x 2 + 4 x 3 2 t 6 = 2 12 x + 12 x 2 2 x 6 x 2 + 4 x 3 2 Γ ( 7 ) Γ ( 7 α ) t 6 α
Substituting Equations (A1)–(A7) into the left side of the equal sign in Equation (49), Equation (49) becomes:
20,000 x 2 ( 1 x ) 2 + 0.24 t 2 + 0.24 Γ ( 3 ) Γ ( 3 α ) t 2 α 0.1 2 12 x + 12 x 2 t 6 2 x 6 x 2 + 4 x 3 2 0.01 2 12 x + 12 x 2 2 x 6 x 2 + 4 x 3 2 Γ ( 7 ) Γ ( 7 α ) t 6 α + 0.002 t x 2 2 x 3 + x 4 = f ( x , t )
It can be obtained that Equation (A8) is equal to Equation (50).
In addition, for boundary conditions:
ω ( 0 , t ) = x 2 ( 1 x ) 2 t 2 | x = 0 = 0
ω ( 1 , t ) = x 2 ( 1 x ) 2 t 2 | x = 1 = 0
ω ( 0 , t ) x = ω ( x , t ) x | x = 0 = 2 ( x 3 x 2 + 2 x 3 ) t 2 | x = 0 = 0
ω ( 1 , t ) x = ω ( x , t ) x | x = 1 = 2 ( x 3 x 2 + 2 x 3 ) t 2 | x = 1 = 0
The initial conditions are as follows:
ω ( x , 0 ) = x 2 ( 1 x ) 2 t 2 | t = 0 = 0
ω ( x , 0 ) t = 2 x 2 ( 1 x ) 2 t | t = 0 = 0

Appendix B

Table A1. Nomenclature.
Table A1. Nomenclature.
SymbolExplanation
c D x α Caputo fractional order derivative operator
c D x β fractional integration
Across-sectional area
ρ density
Emodulus of elasticity
Tkinetic energy
U p potential energy
ttime
xposition
Voccupied area volume
σ stress
ε strain
mdeviation part of the couple stress tensor
γ symmetric curvature tensor
zvertical coordinate
θ rotation vector
μ shear modulus of elasticity
Llength scale parameter
IMoment of inertia
aViscous damping coefficient
η d viscoelastic coefficient
hheight
Φ n ( u ) , Φ n ( x ) , Φ n ( t ) Family of shifted Chebyshev polynomials
R i ( s ) Chebyshev polynomials
R ¯ i ( u ) , R ¯ i ( x ) , R ¯ j ( t ) Shifted Chebyshev polynomials
C n , U , H , H α , W Coefficient matrix
Q x ( m ) Integer order operator matrix
Q t ( α ) fractional operator matrix
ω ( x , t ) analytical solution
ω n ( x , t ) numerical solution
e ( x , t ) absolute error
f ( x , t ) load
Ω frequency

References

  1. Antil, H.; Brown, T.S.; Khatri, R.; Onwunta, A.; Verma, D.; Warma, M. Optimal Control, Numerics, and Applications of Fractional PDEs. arXiv 2021, arXiv:2106.13289. [Google Scholar]
  2. Sheikh, N.A.; Chuan Ching, D.L.; Khan, I.; Ahmad, A.; Ammad, S. Concrete Based Jeffrey Nanofluid Containing Zinc Oxide Nanostructures: Application in Cement Industry. Symmetry 2020, 12, 1037. [Google Scholar] [CrossRef]
  3. Alcoutlabi, M.; Martinez-Vega, J.J. Application of fractional calculus to viscoelastic behavior modelling and to the physical ageing phenomenon in glassy amorphous polymers. Polymer 1998, 39, 6269–6277. [Google Scholar] [CrossRef]
  4. Sun, H.G.; Yong, Z.; Baleanu, D.; Wen, C.; Chen, Y.Q. A new collection of real world applications of fractional calculus in science and engineering. Commun. Nonlinear Sci. Numer. Simul. 2018, 64, 213–231. [Google Scholar] [CrossRef]
  5. Song, Y.; Shateyi, S. Inverse Multiquadric Function to Price Financial Options under the Fractional Black—Scholes Model. Fractal Fract. 2022, 6, 599. [Google Scholar] [CrossRef]
  6. Rexma Sherine, V.; Chellamani, P.; Ismail, R.; Avinash, N.; Britto Antony Xavier, G. Estimating the Spread of Generalized Compartmental Model of Monkeypox Virus Using a Fuzzy Fractional Laplace Transform Method. Symmetry 2022, 14, 2545. [Google Scholar] [CrossRef]
  7. Enelund, M.; Lesieutre, G.A. Time domain modeling of damping using anelastic displacement fields and fractional calculus. Int. J. Solids Struct. 1999, 36, 4447–4472. [Google Scholar] [CrossRef]
  8. Tarasov, V.E. Generalized Memory: Fractional Calculus Approach. Fractal Fract. 2018, 2, 23. [Google Scholar] [CrossRef] [Green Version]
  9. Zguaid, K.; Alaoui, F.; Boutoulout, A. Regional observability for linear time fractional systems. Math. Comput. Simul. 2021, 185, 77–87. [Google Scholar] [CrossRef]
  10. Cb, A.; Th, A.; Ypl, B.; Ga, C.; Msod, E.; Hj, F.; Aaa, G. Analytical solutions of fractional wave equation with memory effect using the fractional derivative with exponential kernel. Results Phys. 2021, 25, 104196. [Google Scholar]
  11. Shen, L.J. Fractional derivative models for viscoelastic materials at finite deformations. Int. J. Solids Struct. 2019, 190, 226–237. [Google Scholar] [CrossRef]
  12. Atanackovic, T.M.; Pilipovic, S. Zener Model with General Fractional Calculus: Thermodynamical Restrictions. Fractal Fract. 2022, 6, 617. [Google Scholar] [CrossRef]
  13. Nikan, O.; Avazzadeh, Z. Numerical simulation of fractional evolution model arising in viscoelastic mechanics. Appl. Numer. Math. 2021, 169, 303–320. [Google Scholar] [CrossRef]
  14. Hendy, M.H.; Amin, M.M.; Ezzat, M.A. Two-dimensional problem for thermoviscoelastic materials with fractional order heat transfer. J. Therm. Stress 2019, 42, 1298–1315. [Google Scholar] [CrossRef]
  15. Bonfanti, A.; Kaplan, J.L.; Charras, G.; Kabla, A. Fractional viscoelastic models for power-law materials. Soft Matter 2020. [Google Scholar] [CrossRef]
  16. Lewandowski, R.; Wielentejczyk, P. Nonlinear vibration of viscoelastic beams described using fractional order derivatives. J. Sound Vib. 2017, 399, 228–243. [Google Scholar] [CrossRef]
  17. Duan, J.S.; Hu, D.C.; Chen, Y.Q. Simultaneous Characterization of Relaxation, Creep, Dissipation, and Hysteresis by Fractional-Order Constitutive Models. Fractal Fract. 2021, 5, 36. [Google Scholar] [CrossRef]
  18. Usuki, T.; Suzuki, T. Dispersion curves for a viscoelastic Timoshenko beam with fractional derivatives. J. Sound Vib. 2012, 331, 605–621. [Google Scholar] [CrossRef]
  19. Catania, G.; Fasana, A.; Sorrentino, S. Finite element analysis of vibrating non-homogeneous beams with fractional derivative viscoelastic models. IFAC Proc. Vol. 2006, 39, 280–285. [Google Scholar] [CrossRef]
  20. Xu, J.; Chen, Y.; Tai, Y.; Xu, X.; Chen, N. Vibration analysis of complex fractional viscoelastic beam structures by the wave method. Int. J. Mech. Sci. 2019, 167, 105204. [Google Scholar] [CrossRef]
  21. Podlubny, I. Fractional Differential Equations. In Proceedings of the Mathematics in Science and Engineering, San Diego, CA, USA, 15 January 1999. [Google Scholar]
  22. Loghman, E.; Bakhtiari-Nejad, F.; Ali, K.E.; Abbaszadeh, M.; Amabili, M. Nonlinear vibration of fractional viscoelastic micro-beams. Int. J. Non-Linear Mech. 2021, 137, 103811. [Google Scholar] [CrossRef]
  23. Cao, J.; Wang, Z.; Wang, Z. A Uniform Accuracy High-Order Finite Difference and FEM for Optimal Problem Governed by Time-Fractional Diffusion Equation. Fractal Fract. 2022, 6, 475. [Google Scholar] [CrossRef]
  24. Loghman, E.; Kamali, A.; Bakhtiari-Nejad, F.; Abbaszadeh, M. Nonlinear free and forced vibrations of fractional modeled viscoelastic FGM micro-beam. Appl. Math. Model. 2021, 92, 297–314. [Google Scholar] [CrossRef]
  25. Attia, M.A.; Melaibari, A.; Shanab, R.A.; Eltaher, M.A. Dynamic Analysis of Sigmoid Bidirectional FG Microbeams under Moving Load and Thermal Load: Analytical Laplace Solution. Mathematics 2022, 10, 4797. [Google Scholar] [CrossRef]
  26. Eltayeb, H.; Kılıçman, A.; Bachar, I. On the Application of Multi-Dimensional Laplace Decomposition Method for Solving Singular Fractional Pseudo-Hyperbolic Equations. Fractal Fract. 2022, 6, 690. [Google Scholar] [CrossRef]
  27. Birol; Bayram, M. Analytical approximate solution of time-fractional Fornberg—Whitham equation by the fractional variational iteration method. Alex. Eng. J. 2014, 53, 911–915. [Google Scholar]
  28. Duan, J.S.; Chaolu, T.; Rach, R.; Lu, L. The Adomian decomposition method with convergence acceleration techniques for nonlinear fractional differential equations. Comput. Math. Appl. 2013, 66, 728–736. [Google Scholar] [CrossRef]
  29. Abbas, M.; Bibi, A.; Alzaidi, A.S.M.; Nazir, T.; Majeed, A.; Akram, G. Numerical Solutions of Third-Order Time-Fractional Differential Equations Using Cubic B-Spline Functions. Fractal Fract. 2022, 6, 528. [Google Scholar] [CrossRef]
  30. Chen, Y.M.; Wei, Y.Q.; Liu, D.Y.; Yu, H. Numerical solution for a class of nonlinear variable order fractional differential equations with Legendre wavelets. Appl. Math. Lett. 2015, 46, 83–88. [Google Scholar] [CrossRef]
  31. Jin, S.; Xie, J.; Qu, J.; Chen, Y. A Numerical Method for Simulating Viscoelastic Plates Based on Fractional Order Model. Fractal Fract. 2022, 6, 150. [Google Scholar] [CrossRef]
  32. Chen, Y.; Han, X.; Liu, L. Numerical Solution for a Class of Linear System of Fractional Differential Equations by the HaarWavelet Method and the Convergence Analysis. Comput. Model. Eng. Sci. 2014, 97, 391–405. [Google Scholar]
  33. Roodgar Saffari, P.; Sher, W.; Thongchom, C. Size Dependent Buckling Analysis of a FG-CNTRC Microplate of Variable Thickness under Non-Uniform Biaxial Compression. Buildings 2022, 12, 2238. [Google Scholar] [CrossRef]
  34. Shen, G.; Xia, Y.; Hu, P.; Zheng, G. Construction of peridynamic beam and shell models on the basis of the micro-beam bond obtained via interpolation method. Eur. J. Mech. A Solids 2021, 86, 104174. [Google Scholar] [CrossRef]
  35. Al-Basyouni, K.S.; Tounsi, A.; Mahmoud, S.R. Size dependent bending and vibration analysis of functionally graded micro beams based on modified couple stress theory and neutral surface position. Compos. Struct. 2015, 125, 621–630. [Google Scholar] [CrossRef]
  36. Ysla, B.; Tx, A. Free vibration of the one-dimensional piezoelectric quasicrystal microbeams based on modified couple stress theory. Appl. Math. Model. 2021, 96, 733–750. [Google Scholar]
  37. Zhang, Z.J.; Zhang, Q.C.; Li, F.C.; Yang, J.W.; Liu, J.W.; Liu, Z.Y.; Jin, F. Modal characteristics of micro-perforated sandwich beams with square honeycomb-corrugation hybrid cores: A mixed experimental-numerical study. THin-Walled Struct. 2019, 137, 185–196. [Google Scholar] [CrossRef]
  38. Shishkina, E.; Sitnik, S. Basics of fractional calculus and fractional order differential equations—ScienceDirect. Transmutations, Singular and Fractional Differential Equations with Applications to Mathematical Physics; Academic Press: Cambridge, MA, USA, 2020; pp. 53–84. [Google Scholar]
  39. Rossikhin, Y.A.; Shitikova, M.V. On fallacies in the decision between the Caputo and Riemann–Liouville fractional derivatives for the analysis of the dynamic response of a nonlinear viscoelastic oscillator. Mech. Res. Commun. 2012, 45, 22–27. [Google Scholar] [CrossRef]
  40. Scherer, R.; Kalla, S.L.; Tang, Y.; Huang, J. The Grünwald–Letnikov method for fractional differential equations. Comput. Math. Appl. 2011, 62, 902–917. [Google Scholar] [CrossRef] [Green Version]
  41. Tp, A.; Atd, B. Numerical solution method for multi-term variable order fractional differential equations by shifted Chebyshev polynomials of the third kind. Alex. Eng. J. 2021, 61, 5145–5153. [Google Scholar]
  42. Farokhi, H.; Ghayesh, M.H.; Amabili, M. Nonlinear dynamics of a geometrically imperfect microbeam based on the modified couple stress theory. Int. J. Eng. Sci. 2013, 68, 11–23. [Google Scholar] [CrossRef]
  43. LeonardMeirovitch. Fundamentals of Vibrations; Waveland Pr Inc.: Long Grove, IL, USA, 2010. [Google Scholar]
  44. Zhang, Y.C.; Jun-Zheng, W.U.; Zhang, C.Y.; Zhang, N.H. Free Vibration of a Viscoelastic-Elastic Laminated Microcantilever Beam. Chin. Q. Mech. 2017. [Google Scholar] [CrossRef]
  45. Karunakar, P.; Chakraverty, S. Shifted Chebyshev polynomials based solution of partial differential equations. Sn Appl. Sci. 2019, 1, 285. [Google Scholar] [CrossRef]
  46. Akrami, M.H.; Atabakzadeh, M.H.; Erjaee, G.H. The operational matrix of fractional integration for shifted Legendre polynomials. Iran. J. Sci. Technol. Trans. Sci. 2013, 37, 439–444. [Google Scholar]
  47. Atabakzadeh, M.H.; Akrami, M.H.; Erjaee, G.H. Chebyshev operational matrix method for solving multi-order fractional ordinary differential equations. Appl. Math. Model. 2013, 37, 8903–8911. [Google Scholar] [CrossRef]
  48. Zhao, F.; Huang, Q.; Xie, J.; Li, Y.; Ma, L.; Wang, J. Chebyshev polynomials approach for numerically solving system of two-dimensional fractional PDEs and convergence analysis. Appl. Math. Comput. 2017, 313, 321–330. [Google Scholar] [CrossRef]
  49. Meng, K.; Cui, C.; Liang, Z.; Li, H.; Pei, H. An Analytical Solution for Longitudinal Impedance of a Large-Diameter Floating Pile in Soil with Radial Heterogeneity and Viscous-Type Damping. Appl. Sci. 2020, 10, 4906. [Google Scholar] [CrossRef]
  50. Sheng, G.G.; Wang, X. Nonlinear forced vibration of functionally graded Timoshenko microbeams with thermal effect and parametric excitation. Int. J. Mech. Sci. 2019, 155, 405–416. [Google Scholar] [CrossRef]
  51. Wang, L.; Chen, Y.; Cheng, G.; Barrière, T. Numerical analysis of fractional partial differential equations applied to polymeric visco-elastic Euler–Bernoulli beam under quasi-static loads. Chaos Solitons Fractals 2020, 140, 110255. [Google Scholar] [CrossRef]
  52. Samayoa, D.; Kryvko, A.; Velázquez, G.; Mollinedo, H. Fractal Continuum Calculus of Functions on Euler–Bernoulli Beam. Fractal Fract. 2022, 6, 552. [Google Scholar] [CrossRef]
Figure 1. Bending deformation of viscoelastic micro-beams under simple harmonic load f ( x , t ) .
Figure 1. Bending deformation of viscoelastic micro-beams under simple harmonic load f ( x , t ) .
Fractalfract 07 00204 g001
Figure 2. Numerical example results.
Figure 2. Numerical example results.
Fractalfract 07 00204 g002
Figure 3. Effect of different viscous damping coefficients on deflection of the micro-beam.
Figure 3. Effect of different viscous damping coefficients on deflection of the micro-beam.
Fractalfract 07 00204 g003
Figure 4. Deflection of micro-beams of different length scale parameters.
Figure 4. Deflection of micro-beams of different length scale parameters.
Fractalfract 07 00204 g004
Figure 5. Deflection of the micro-beam under different simple harmonic loads.
Figure 5. Deflection of the micro-beam under different simple harmonic loads.
Fractalfract 07 00204 g005
Figure 6. Stress and strain in the viscoelastic micro-beam under f ( x , t ) = 100 ( x 2 ) c o s t .
Figure 6. Stress and strain in the viscoelastic micro-beam under f ( x , t ) = 100 ( x 2 ) c o s t .
Fractalfract 07 00204 g006
Figure 7. Potential energy change of the viscoelastic micro-beam.
Figure 7. Potential energy change of the viscoelastic micro-beam.
Fractalfract 07 00204 g007
Table 1. Absolute error of each point when α = 0.2 .
Table 1. Absolute error of each point when α = 0.2 .
e ( x , t ) t = 0.1 t = 0.3 t = 0.5 t = 0.7
x = 0.1 5.21 × 10 8 8.94 × 10 8 3.25 × 10 8 5.41 × 10 7
x = 0.5 6.41 × 10 8 8.37 × 10 8 1.78 × 10 8 1.46 × 10 7
x = 0.9 1.10 × 10 7 8.82 × 10 8 1.01 × 10 7 1.46 × 10 7
Table 2. Absolute error of each point when α = 0.5 .
Table 2. Absolute error of each point when α = 0.5 .
e ( x , t ) t = 0.1 t = 0.3 t = 0.5 t = 0.7
x = 0.1 1.77 × 10 9 7.44 × 10 9 1.76 × 10 9 5.53 × 10 9
x = 0.5 3.32 × 10 10 8.26 × 10 9 1.15 × 10 9 6.82 × 10 9
x = 0.9 1.77 × 10 9 7.54 × 10 9 1.68 × 10 9 5.53 × 10 9
Table 3. Absolute error of each point when α = 0.8 .
Table 3. Absolute error of each point when α = 0.8 .
e ( x , t ) t = 0.1 t = 0.3 t = 0.5 t = 0.7
x = 0.1 7.01 × 10 9 1.34 × 10 9 1.52 × 10 8 6.17 × 10 9
x = 0.5 4.85 × 10 8 1.27 × 10 9 1.65 × 10 8 6.85 × 10 9
x = 0.9 7.04 × 10 8 1.34 × 10 9 1.54 × 10 9 5.17 × 10 9
Table 4. Absolute error of each point when α = 1 .
Table 4. Absolute error of each point when α = 1 .
e ( x , t ) t = 0.1 t = 0.3 t = 0.5 t = 0.7
x = 0.1 1.80 × 10 3 2.80 × 10 3 1.33 × 10 2 3.57 × 10 2
x = 0.5 1.90 × 10 3 3.00 × 10 3 1.46 × 10 2 4.81 × 10 2
x = 0.9 3.20 × 10 3 4.70 × 10 3 1.05 × 10 2 6.62 × 10 2
Table 5. Values of the parameters.
Table 5. Values of the parameters.
E (GPa)h (m) ρ (kg/m 3 ) μ (GPa) η d (GPa)l (m)
21 4.166 × 10 6 7850 0.956 ( μ / E ) α 5 × 10 4
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Yang, A.; Zhang, Q.; Qu, J.; Cui, Y.; Chen, Y. Solving and Numerical Simulations of Fractional-Order Governing Equation for Micro-Beams. Fractal Fract. 2023, 7, 204. https://0-doi-org.brum.beds.ac.uk/10.3390/fractalfract7020204

AMA Style

Yang A, Zhang Q, Qu J, Cui Y, Chen Y. Solving and Numerical Simulations of Fractional-Order Governing Equation for Micro-Beams. Fractal and Fractional. 2023; 7(2):204. https://0-doi-org.brum.beds.ac.uk/10.3390/fractalfract7020204

Chicago/Turabian Style

Yang, Aimin, Qunwei Zhang, Jingguo Qu, Yuhuan Cui, and Yiming Chen. 2023. "Solving and Numerical Simulations of Fractional-Order Governing Equation for Micro-Beams" Fractal and Fractional 7, no. 2: 204. https://0-doi-org.brum.beds.ac.uk/10.3390/fractalfract7020204

Article Metrics

Back to TopTop