A Numerical Method for Simulating Viscoelastic Plates Based on Fractional Order Model

Abstract

In this study, an efficacious method for solving viscoelastic dynamic plates in the time domain is proposed for the first time. The differential operator matrices of different orders of Bernstein polynomials algorithm are adopted to approximate the ternary displacement function. The approximate results are simulated by code. In addition, it is proved that the proposed method is feasible and effective through error analysis and mathematical examples. Finally, the effects of external load, side length of plate, thickness of plate and boundary condition on the dynamic response of square plate are studied. The numerical results illustrate that displacement and stress of the plate change with the change of various parameters. It is further verified that the Bernstein polynomials algorithm can be used as a powerful tool for numerical solution and dynamic analysis of viscoelastic plates.

1. Introduction

Plate and plate structure are widely used in many realms of mechanical, building and aerospace [1]. In addition, the plate is also a key component in aerospace engineering, which bears strong and sudden power including vibration. They are usually combined with viscoelastic materials to reduce the applied vibration and have the effect of damping [2]. In order to apply plate structure more widely in real life, many scholars are committed to the research of plate vibration. Early scholars analyzed the linear vibration of plates. Ziaee [3] used the Ritz method to study linear vibration of nanoplates. The effects of different parameters on polysilicon microplate were discussed. Cadou et al. [4] considered the linear vibration of the plate and testified availability of the method used for calculating eigenvalues of linear problems. Gradually, the research direction of some scholars began to change to nonlinear plate. Based on the semi-analytical method, Babahammou et al. [5] researched linear and nonlinear vibration of plates under the condition of full line or partial line support. Cho [6] proposed an algorithm for analyzing nonlinear vibration problems. The parameters of nonlinear free vibration characteristics of composite plates were discussed. Quan et al. [7] established the governing equation of sandwich plate vibration by using shear deformation theory and analyzed the nonlinear vibration of plate by Galerkin method and fourth-order Runge Kutta method. Although there have been a lot of research results, there are still problems in viscoelastic plates.

Viscoelastic materials have both elasticity and viscosity and are widely used for passive vibration isolation in engineering structures and applications owing to their light weight and high intensity. Therefore, it is essential to model viscoelastic materials properly [8,9]. In recent years, scholars have proposed various models to describe viscoelastic material. In order to simulate the elastic and viscous properties of materials at the same time, the fractional viscoelastic model is replaced by the classical viscoelastic model. The behavior of the system cannot only be appropriately described by fractional order model with less parameters but also fitted by fractional order operator [10]. Therefore, fractional order has a wide range of applications, especially in control systems. Using the multi-switch synchronization method, Pan et al. [11] considered the sliding-mode combinatorial synchronization of fractional-order chaotic systems under double random disturbances. Zhang et al. [12] introduced fractional order into sliding mode control of the system. The nonlinear term is estimated by using radial basis function neural network. Zhang et al. [13] provided a set of criteria for fractional order systems stability and verified the efficiency of controllers with numerical examples. With the development of viscoelastic material structure, some scholars gradually use fractional order to model viscoelastic plate. Rouzegar et al. [14] derived the governing equations of viscoelastic plate with Voigt viscoelastic model. The variation of amplitude and frequency of fractional viscoelastic plates under external excitation was given. Permoon et al. [15] discussed the natural frequency and characteristics of viscoelastic plate. Three constitutive models were compared by fitting a curve. Praharaj et al. [16] employed fractional damping derivative model to simulate plate structure. How the different orders of fractional order affect the vibration response of the plate was considered by combining the two methods to solve the differential equations. Ai et al. [17] established the finite element equations of stiffened plate by using fractional merchant model. The influences of altitude and overall arrangement on time-varying behavior of plates were analyzed.

The research on the numerical solution of viscoelastic plates is equivalent to the solution of fractional governing equation. In other words, numerical analysis of the viscoelastic plates not only needs to establish the material behavior equation but also an effective numerical method to approximate fractional governing dynamic equations. The common methods for solving the mentioned equations include Laplace transform [18], Fourier transform [19], Galerkin method [20], meshless method [21], multi-scale method [22] and variational iteration method [23]. Due to the large amount of calculation and difficulty to obtain the inverse transform, these methods are very difficult for obtaining solutions of this kind of equation directly in the time domain. However, plate differential governing equations take critical part in the wide application of engineering science. So, in recent years, polynomial approximation method has been widely used in solving fractional differential equations. Wang et al. [24] used shifted Legendre polynomials to approximate variable fractional differential equations. The dynamic response of viscoelastic pipe conveying fluid was analyzed. Hashim et al. [25] proposed shifted Chebyshev polynomials of the second kind to solve approximate solutions of time-delay variable fractional differential equations. Cao et al. [26] studied a significant method based on fractional rheological model to solve viscoelastic column problems. The fractional differential equation was solved by shifted Chebyshev wavelet function. Compared with the above polynomials, Bernstein polynomials have the advantages of simple structure and perfect properties. Therefore, it is widely used in solving differential equations and practical engineering applications [27]. More and more scholars use Bernstein polynomials to solve all kinds of differential equations. By using Bernstein polynomials, Khan et al. [28] obtained the calculation results of fractional constitutive equation. In this method, the coupled system is transformed into algebraic equations by operator matrix. Heydari et al. [29] proposed Bernstein polynomials to approximate advection diffusion reaction equation with fractal fractional derivative. Chen et al. [30] utilized Bernstein polynomials to solve a series of variable fractional order differential equations. Different Bernstein polynomials matrices were derived and used to transform initial equations into discrete nonlinear equations. However, the above studies are limited to solving fractional differential equations using Bernstein polynomials. Few studies directly solve such equations in time domain and analyze the dynamic behavior combined with three-dimensional plates.

For these reasons, it is necessary to combine the fractional order model with a new calculation method to settle the above problem. Therefore, the Bernstein polynomials algorithm is proposed to numerically simulate differential equations of the plate and analyze the effects of parameters on the numerical solutions of the displacement and stress. This algorithm has good applicability to calculate the fractional governing equation of three-dimensional plates in the time domain. It is also fit for dynamic analysis of viscoelastic plates. In addition, this technology will supply a new approach for the numerical study of viscoelastic plates.

The paper is structured as follows: The preliminary knowledge of Caputo derivative is recommended in Section 2. Section 3 presents the governing equation of three-dimensional viscoelastic plate by using the fractional constitutive model. Section 4 derives the matrices of Bernstein polynomials. Demand problem is expressed by various orders differential operator matrices. In Section 5, the availability of the proposed algorithm is verified by two methods. Section 6 discusses displacements of plate under various conditions. Finally, conclusions of this paper are obtained in Section 7.

2. Preliminaries

Next, several properties and mathematical preliminaries of fractional calculus are given, which will be applied in the following sections.

Definition 1.

The fractional Caputo derivative of order α is defined as [31]

$$\left(D_t^{\alpha }f\right)\left(t\right)=\left\{\begin{array}{c}\frac{1}{\mathsf{\Gamma }\left(r-\alpha \right)}\left[{\int }_0^t\frac{f^{\left(r\right)}\left(\tau \right)}{{\left(t-\tau \right)}^{\alpha +1-r}}d\tau \right],-1\le \alpha <r,\hfill \\ \frac{d^r}{d{t}^r}f\left(t\right),\alpha =r.\hfill \end{array}\right.$$

(1)

where $D_t^{\alpha }$ is Caputo fractional differential operator, $0<\alpha \le 1$ is fractional order, f is integrable and continuous on $(0,+\infty )$, $\mathsf{\Gamma }\left(·\right)$ is Gamma function and has $\mathsf{\Gamma }\left(z\right)={\int }_0^{\infty }{e}^{-t}{t}^{z-1}dt$.

Based on the Caputo derivative, we have

$$D_t^{\alpha }{t}^r=\left\{\begin{array}{c}0,r=0\\ \frac{\mathsf{\Gamma }\left(r+1\right)}{\mathsf{\Gamma }\left(r+1-\alpha \right)}{t}^{r-\alpha },r=1,2,3,\dots \end{array}\right.$$

(2)

Three properties of Caputo fractional differential are as follows, for $\lambda ,\mu \in R,0<\alpha \le 1,$C is a constant, and $f\left(t\right)\in C^1\left(R\right)$ [26,32]

$$\begin{array}{c}\left(1\right)D_t^{\alpha }C=0\hfill \\ \left(2\right)D_t^{\alpha }\left(Cf\left(t\right)\right)=C{D}_t^{\alpha }f\left(t\right)\hfill \\ \left(3\right)D_t^{\alpha }\left[\lambda f\left(t\right)+\mu g\left(t\right)\right]=\lambda D_t^{\alpha }f\left(t\right)+\mu D_t^{\alpha }g\left(t\right)\hfill \end{array}$$

(3)

3. Governing Equation of Fractional Viscoelastic Plate

In this part, a fractional viscoelastic plate shown in Figure 1 with sides of a and b, thickness of h is studied. This rectangular plate is fixedly supported on four sides. The motion equation of the plate is as follows [33]

$$\rho h\frac{{\partial }^{2}u}{\partial t^2}+\eta D_t^{\alpha }u+D{\nabla }^{4}u-h\left(\frac{{\partial }^2\psi }{\partial y^2}\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^2\psi }{\partial x^2}\frac{{\partial }^{2}u}{\partial y^2}-2\frac{{\partial }^2\psi }{\partial x\partial y}\frac{{\partial }^{2}u}{\partial x\partial y}\right)=F\left(x,y,t\right)$$

(4)

where $\rho$ is mass density, $\eta$ is constant damping for the fractional derivative element model, ${\nabla }^4=\frac{{\partial }^4}{\partial x^4}+\frac{{\partial }^4}{\partial y^4}+2\frac{{\partial }^4}{\partial x^2\partial y^2}$ is biharmonic operator, $D=\frac{E{h}^3}{12(1-v^2)}$ is flexural stiffness, E is Young’s modulus, v is the Poisson’s ratio, $\psi$ is the Airy stress function and $F(x,y,t)$ is transverse force excitation.

From the above formulas, the governing equation of fractional viscoelastic plate is expressed as [34,35]

$$\rho h\frac{{\partial }^{2}u}{\partial t^2}+\eta D_t^{\alpha }u+\left(\frac{E{h}^3}{12(1-v^2)}-h\right)\left(\frac{{\partial }^{4}u}{\partial x^4}+\frac{{\partial }^{4}u}{\partial y^4}\right)+2\left(\frac{E{h}^3}{12(1-v^2)}+h\right)\frac{{\partial }^{4}u}{\partial x^2\partial y^2}=F\left(x,y,t\right)$$

(5)

The boundary conditions of plate with four fixed edges are

$$\begin{array}{c}u\left(x,y,t\right)\left|{}_{x=0,a}\right.=\frac{\partial u\left(x,y,t\right)}{\partial x}\left|{}_{x=0,a}\right.=0\hfill \\ u\left(x,y,t\right)\left|{}_{y=0,b}\right.=\frac{\partial u\left(x,y,t\right)}{\partial y}\left|{}_{y=0,b}\right.=0\hfill \end{array}$$

(6)

4. Numerical Algorithm for Bernstein Polynomials

In this section, Bernstein polynomials algorithm is introduced to approximate unknown functions. The governing equations are transformed from differential operator matrices into algebraic equations.

4.1. Bernstein Polynomials

Definition 2.

The definition of Bernstein polynomials of degree m in $[0,1]$ is [36]

$${\mathcal{B}}_{i,m}\left(x\right)=\left(\begin{array}{c}m\\ i\end{array}\right)x^i{\left(1-x\right)}^{m-i},0\le i\le m.$$

(7)

The following formula is written as

$${\mathcal{B}}_{i,m}\left(x\right)=\left(1-x\right){\mathcal{B}}_{i,m-1}\left(x\right)+x{\mathcal{B}}_{i-1,m-1}\left(x\right),i=0,1,\dots ,m.$$

(8)

Expand the binomial ${\left(1-x\right)}^{m-i}$ to obtain the following formula

$$\begin{array}{c}{\mathcal{B}}_{i,m}\left(x\right)=\left(\begin{array}{c}m\\ i\end{array}\right)x^i{\left(1-x\right)}^{m-i}={\displaystyle \sum _{k=0}^{m-i}}{\left(-1\right)}^k\left(\begin{array}{c}m\\ i\end{array}\right)\left(\begin{array}{c}m-i\\ k\end{array}\right)x^{i+k},i=0,1,\dots ,m.\hfill \end{array}$$

(9)

For expanding the scope of x, Bernstein polynomials in $[0,R]$ is written as

$$B_{i,m}\left(x\right)=\left(\begin{array}{c}m\\ i\end{array}\right)\frac{x^i{\left(R-x\right)}^{m-i}}{R^m}={\displaystyle \sum _{k=0}^{m-i}}{\left(-1\right)}^k\left(\begin{array}{c}m\\ i\end{array}\right)\left(\begin{array}{c}m-i\\ k\end{array}\right)\frac{x^{i+k}}{R^{i+k}}$$

(10)

where R is any positive integer.

A matrix $\phi \left(x\right)$ consisted by Bernstein polynomials is

$$\phi \left(x\right)={\left[B_{0,m}\left(x\right),B_{1,m}\left(x\right),\dots ,B_{m,m}\left(x\right)\right]}^{\mathrm{T}}=HK\left(x\right)$$

(11)

where

$$H={\left[h_{i,j}\right]}_{i,j=0}^m,h_{i,j}=\left\{\begin{array}{c}\frac{{\left(-1\right)}^{j-i}\left(\begin{array}{c}m\\ i\end{array}\right)\left(\begin{array}{c}m-i\\ j-i\end{array}\right)}{R^j},j\ge i\\ 0,j<i\end{array}\right.$$

(12)

$$K\left(x\right)={\left[1,x,\dots ,x^m\right]}^{\mathrm{T}}$$

(13)

4.2. Function Approximation

For the displacement function $u(x,y,t)$ of Equation (5), it is expressed as [37]

$$u(x,y,t)=u(x,y)u\left(t\right)$$

(14)

where $u\left(x,y\right)\in L^2\left(\left[0,R\right]\times \left[0,S\right]\right)$ and $u\left(t\right)\in L^2\left(\left[0,T\right]\right)$.

For any continuous function $u\left(x\right)\in L^2\left(\left[0,R\right]\right)$ with one variable, it can be expressed by Bernstein polynomials as

$$u\left(x\right)={\displaystyle \sum _{i=0}^{\infty }}{n}_i{B}_{i,m}\left(x\right)\approx {\displaystyle \sum _{i=0}^m}{n}_i{B}_{i,m}\left(x\right)=N^{\mathrm{T}}\phi \left(x\right)$$

(15)

where $N^{\mathrm{T}}=\left[n_0,n_1,\dots ,n_m\right]$.

The inner product is calculated as

$$〈u\left(x\right),{\phi }^{\mathrm{T}}\left(x\right)〉=N^{\mathrm{T}}〈\phi \left(x\right),{\phi }^{\mathrm{T}}\left(x\right)〉=N^{\mathrm{T}}{\left[o{}_{i,j}\right]}_{i,j=0}^m=N^{\mathrm{T}}O$$

(16)

where $o_{i,j}={\int }_0^R{B}_{i,m}\left(x\right)B_{j,m}\left(x\right)dx,N^{\mathrm{T}}=〈u\left(x\right),{\phi }^{\mathrm{T}}\left(x\right)〉O^{-1}$.

Similarly, the function of two variables $u\left(x,y\right)\in L^2\left(\left[0,R\right]\times \left[0,S\right]\right)$ can be approximated as

$$u\left(x,y\right)={\displaystyle \sum _{i=0}^{\infty }}{\displaystyle \sum _{j=0}^{\infty }}{w}_{i,j}{B}_{i,m}\left(x\right)B_{j,m}\left(y\right)\approx {\displaystyle \sum _{i=0}^m}{\displaystyle \sum _{j=0}^m}{w}_{i,j}{B}_{i,m}\left(x\right)B_{j,m}\left(y\right)={\phi }^T\left(x\right)W\phi \left(y\right)$$

(17)

where $W={\left[w_{i,j}\right]}_{i,j=0}^m$.

The function of t is approximated as

$$u\left(t\right)={\displaystyle \sum _{k=0}^{\infty }}{n}_k{B}_{k,m}\left(t\right)\approx {\displaystyle \sum _{k=0}^m}{n}_k{B}_{k,m}\left(t\right)=N^{\mathrm{T}}\phi \left(t\right)$$

(18)

where $\phi \left(t\right)={\left[B_{0,m}\left(t\right),B_{1,m}\left(t\right),\dots ,B_{m,m}\left(t\right)\right]}^{\mathrm{T}}$.

Taking Equations (17) and (18) into Equation (14), the displacement function is written as

$$u\left(x,y,t\right)=u\left(x,y\right)u\left(t\right)={\phi }^{\mathrm{T}}\left(x\right)W\phi \left(y\right)N^{\mathrm{T}}\phi \left(t\right)$$

(19)

4.3. Differential Operator Matrix of Bernstein Polynomials

4.3.1. Integer Differential Operator Matrix

Definition 3.

$G_x^1$ satisfaction ${\phi }^{\prime }\left(x\right)=G_x^1\phi \left(x\right)$ is the 1th operator matrix of Bernstein polynomials. The derivation procedure is as follows

$$\begin{array}{c}\hfill {\phi }^{\prime }\left(x\right)={\left(HK\left(x\right)\right)}^{\prime }=H{K}^{\prime }\left(x\right)=HPK\left(x\right)=HP{H}^{-1}\phi \left(x\right)=G_x^1\phi \left(x\right)\end{array}$$

(20)

where $P={\left[p_{i,j}\right]}_{i,j=0}^m,p_{i,j}=\left\{\begin{array}{c}i,i=j+1\\ 0,i\ne j+1\end{array}\right.,$ and $G_x^1=HP{H}^{-1}$.

Definition 4.

$G_x^2$ satisfaction ${\phi }^{″}\left(x\right)=G_x^2\phi \left(x\right)$ is the 2th operator matrix of Bernstein polynomials. The derivation procedure is as follows

$$\begin{array}{cc}\hfill {\phi }^{″}\left(x\right)& ={\left(HK\left(x\right)\right)}^{″}=H{\left(K^{\prime }\left(x\right)\right)}^{\prime }={\left(HP{H}^{-1}\phi \left(x\right)\right)}^{\prime }\hfill \\ \hfill & =HP{H}^{-1}{\left(\phi \left(x\right)\right)}^{\prime }=\left(HP{H}^{-1}\right)\left(HP{H}^{-1}\right)\phi \left(x\right)=G_x^2\phi \left(x\right)\hfill \end{array}$$

(21)

where $G_x^2={\left(HP{H}^{-1}\right)}^2$.

According to Equations (20) and (21), the matrices of integer order differential operator are obtained

$$\begin{array}{c}\hfill {\phi }^m\left(x\right)={\left(\phi \left(x\right)\right)}^m={\left(HP{H}^{-1}\right)}^m\phi \left(x\right)=G_x^m\phi \left(x\right)\end{array}$$

(22)

$$\begin{array}{c}\hfill {\phi }^m\left(y\right)={\left(\phi \left(y\right)\right)}^m={\left(HP{H}^{-1}\right)}^m\phi \left(y\right)=G_y^m\phi \left(y\right)\end{array}$$

(23)

$$\begin{array}{c}\hfill {\phi }^m\left(t\right)={\left(\phi \left(t\right)\right)}^m={\left(HP{H}^{-1}\right)}^m\phi \left(t\right)=G_t^m\phi \left(t\right)\end{array}$$

(24)

The items in Equation (5) are written as

$$\begin{array}{cc}\hfill \frac{{\partial }^{2}u\left(x,y,t\right)}{\partial t^2}& \approx \frac{{\partial }^2\left({\phi }^{\mathrm{T}}\left(x\right)W\phi \left(y\right)N^{\mathrm{T}}\phi \left(t\right)\right)}{\partial t^2}={\phi }^{\mathrm{T}}\left(x\right)W\phi \left(y\right)N^{\mathrm{T}}\frac{{\partial }^2\phi \left(t\right)}{\partial t^2}\hfill \\ \hfill & ={\phi }^{\mathrm{T}}\left(x\right)W\phi \left(y\right)N^{\mathrm{T}}{\left(HP{H}^{-1}\right)}^2\phi \left(t\right)\hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill \frac{{\partial }^{4}u\left(x,y,t\right)}{\partial x^4}& \approx \frac{{\partial }^4\left({\phi }^{\mathrm{T}}\left(x\right)W\phi \left(y\right)N^{\mathrm{T}}\phi \left(t\right)\right)}{\partial x^4}=\frac{{\partial }^4{\phi }^{\mathrm{T}}\left(x\right)}{\partial x^4}W\phi \left(y\right)N^{\mathrm{T}}\phi \left(t\right)\hfill \\ \hfill & ={\phi }^{\mathrm{T}}\left(x\right){\left({\left(HP{H}^{-1}\right)}^{\mathrm{T}}\right)}^{4}W\phi \left(y\right)N^{\mathrm{T}}\phi \left(t\right)\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill \frac{{\partial }^{4}u\left(x,y,t\right)}{\partial y^4}& \approx \frac{{\partial }^4\left({\phi }^{\mathrm{T}}\left(x\right)W\phi \left(y\right)N^{\mathrm{T}}\phi \left(t\right)\right)}{\partial y^4}={\phi }^{\mathrm{T}}\left(x\right)W\frac{{\partial }^4\phi \left(y\right)}{\partial y^4}{N}^{\mathrm{T}}\phi \left(t\right)\hfill \\ \hfill & ={\phi }^{\mathrm{T}}\left(x\right)W{\left(HP{H}^{-1}\right)}^4\phi \left(y\right)N^{\mathrm{T}}\phi \left(t\right)\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill \frac{{\partial }^{4}u\left(x,y,t\right)}{\partial x^2\partial y^2}& \approx \frac{{\partial }^4\left({\phi }^{\mathrm{T}}\left(x\right)W\phi \left(y\right)N^{\mathrm{T}}\phi \left(t\right)\right)}{\partial x^2\partial y^2}=\frac{{\partial }^4\left({\phi }^{\mathrm{T}}\left(x\right)W\phi \left(y\right)\right)}{\partial x^2\partial y^2}{N}^{\mathrm{T}}\phi \left(t\right)\hfill \\ \hfill & ={\phi }^{\mathrm{T}}\left(x\right){\left({\left(HP{H}^{-1}\right)}^{\mathrm{T}}\right)}^{2}W{\left(HP{H}^{-1}\right)}^2\phi \left(y\right)N^{\mathrm{T}}\phi \left(t\right)\hfill \end{array}$$

(28)

4.3.2. Fractional Differential Operator Matrix

Definition 5.

$G_t^{\alpha }$ satisfaction $D_t^{\alpha }\phi \left(t\right)=D_t^{\alpha }HK\left(t\right)=H{D}_t^{\alpha }K\left(t\right)=HQ{H}^{-1}\phi \left(t\right)=G_t^{\alpha }\phi \left(t\right)$ is the αth operator matrix of Bernstein polynomials.

$$\mathit{where} Q={\left[q_{i,j}\right]}_{i,j=0}^m,q_{i,j}=\left\{\begin{array}{c}\frac{\mathsf{\Gamma }\left(i+1\right)}{\mathsf{\Gamma }\left(i+1-\alpha \right)}{t}^{-\alpha },i=j,i\ge 1\hfill \\ 0,else\hfill \end{array}\right.,\mathit{and}{G}_t^{\alpha }=HQ{H}^{-1}.$$

The partial differential term in Equation (5) is formulated as

$$\begin{array}{cc}\hfill D_t^{\alpha }u\left(x,y,t\right)& =D_t^{\alpha }\left({\phi }^{\mathrm{T}}\left(x\right)W\phi \left(y\right)N^{\mathrm{T}}\phi \left(t\right)\right)={\phi }^{\mathrm{T}}\left(x\right)W\phi \left(y\right)N^{\mathrm{T}}{D}_t^{\alpha }\phi \left(t\right)\hfill \\ \hfill & ={\phi }^{\mathrm{T}}\left(x\right)W\phi \left(y\right)N^{\mathrm{T}}\left(HQ{H}^{-1}\right)\phi \left(t\right)\hfill \end{array}$$

(29)

4.4. Discretizationdiscretization Governing Equation

The governing equation can be expressed by various orders differential operator matrices as

$$\begin{array}{c}\rho h{\phi }^{\mathrm{T}}\left(x\right)W\phi \left(y\right)N^{\mathrm{T}}{\left(HP{H}^{-1}\right)}^2\phi \left(t\right)+\eta {\phi }^{\mathrm{T}}\left(x\right)W\phi \left(y\right)N^{\mathrm{T}}\left(HQ{H}^{-1}\right)\phi \left(t\right)\hfill \\ +\left(\frac{E{h}^3}{12\left(1-v^2\right)}-h\right)\left({\phi }^{\mathrm{T}}\left(x\right){\left({\left(HP{H}^{-1}\right)}^{\mathrm{T}}\right)}^{4}W\phi \left(y\right)N^{\mathrm{T}}\phi \left(t\right)\right)\hfill \\ +\left(\frac{E{h}^3}{12\left(1-v^2\right)}-h\right)\left({\phi }^{\mathrm{T}}\left(x\right)W{\left(HP{H}^{-1}\right)}^4\phi \left(y\right)N^{\mathrm{T}}\phi \left(t\right)\right)+2\left(\frac{E{h}^3}{12\left(1-v^2\right)}+h\right)\hfill \\ \times \left({\phi }^{\mathrm{T}}\left(x\right){\left({\left(HP{H}^{-1}\right)}^{\mathrm{T}}\right)}^{2}W\phi \left(y\right){\left(HP{H}^{-1}\right)}^2{N}^{\mathrm{T}}\phi \left(t\right)\right)=F\left(x,y,t\right)\hfill \end{array}$$

(30)

The boundary conditions are converted into

$$\begin{array}{c}{\phi }^{\mathrm{T}}\left(x\right)W\phi \left(y\right)N^{\mathrm{T}}\phi \left(t\right)\left|{}_{x=0,a}=0\right.\hfill \\ {\phi }^{\mathrm{T}}\left(x\right)W\phi \left(y\right)N^{\mathrm{T}}\phi \left(t\right)\left|{}_{y=0,b}=0\right.\hfill \\ {\phi }^{\mathrm{T}}\left(x\right)\left(HP{H}^{-1}\right)W\phi \left(y\right)N^{\mathrm{T}}\phi \left(t\right)\left|{}_{x=0,a}=0\right.\hfill \\ {\phi }^{\mathrm{T}}\left(x\right)W\left(HQ{H}^{-1}\right)\phi \left(y\right)N^{\mathrm{T}}\phi \left(t\right)\left|{}_{x=0,b}=0\right.\hfill \end{array}$$

(31)

Configure variable $\left(x,y,t\right)$ as discrete variable $\left(x_i,y_j,t_k\right)$. Matrices W and N areis gained by MATLAB software. Thus, the initial equation is solved.

5. Error Analysis and Mathematical Example

For the sake of proving accuracy and effectiveness of the mentioned Bernstein polynomials algorithm, the following error analysis is carried out.

5.1. Error Bound

Theorem 1.

If $u\left(t\right)\in C^{m+1}\left[0,T\right]$ and $u_0\left(t\right)\in Y$ is the best approximation of $u\left(t\right)$. $\{B_{0,m}\left(x\right)$, $B_{1,m}\left(x\right),B_{2,m}\left(x\right),\dots ,B_{m,m}\left(t\right)\}\subset L^2\left(\left[0,T\right]\right)$,

$Y=\mathrm{span}\left\{B_{0,m}\left(x\right),B_{1,m}\left(x\right),B_{2,m}\left(x\right),\dots ,B_{m,m}\left(x\right)\right\}$. Then the expression of the error is [38,39]

$${∥u\left(t\right)-u_0\left(t\right)∥}_2={∥{\epsilon }_{u\left(t\right)}∥}_2<\frac{V}{\left(m+1\right)!}\frac{T^{\frac{2m+3}{2}}}{\sqrt{\left(2m+3\right)}}$$

(32)

where $V={max}_{\left(t\right)\in \left[0,T\right]}\left|\frac{{\partial }^{m+1}u\left(t\right)}{\partial t^{m+1}}\right|$.

Proof. 

The Taylor expansion of $u_1\left(t\right)$ is

$$u_1\left(t\right)=u\left(0\right)+\left(t\frac{\partial }{\partial t}\right)u\left(0\right)+\frac{1}{2!}{\left(t\frac{\partial }{\partial t}\right)}^{2}u\left(0\right)+\dots +\frac{1}{m!}{\left(t\frac{\partial }{\partial t}\right)}^{m}u\left(0\right)$$

(33)

where $u_1\left(t\right)\in Y$

Because $u_0\left(t\right)$ is the best approximation of $u\left(t\right)$, there is

$$\begin{array}{cc}\hfill {∥{\epsilon }_{u\left(t\right)}∥}_2& ={∥u\left(t\right)-u_0\left(t\right)∥}_2\le {∥u\left(t\right)-u_1\left(t\right)∥}_2\hfill \\ \hfill & ={\left({{\int }_0^T\left(\frac{{\left(t\frac{\partial }{\partial t}\right)}^{m+1}u\left(\xi \right)}{\left(m+1\right)!}\right)}^{2}dt\right)}^{\frac{1}{2}}\hfill \\ \hfill & ={\left({{\int }_0^T\left(\frac{{\displaystyle \sum _{k=0}^{m+1}}\left(\begin{array}{c}m+1\\ k\end{array}\right)\frac{{\partial }^{m+1}u\left(\xi \right)}{\partial t^{m+1}}{t}^{m+1}}{\left(m+1\right)!}\right)}^{2}dt\right)}^{\frac{1}{2}}\hfill \\ \hfill & \le \frac{V}{\left(m+1\right)!}{\left({\int }_0^T{t}^{2m+2}dt\right)}^{\frac{1}{2}}\hfill \\ \hfill & \le \frac{V}{\left(m+1\right)!}\frac{T^{\frac{2m+3}{2}}}{\sqrt{2m+3}}\hfill \end{array}$$

(34)

where $\xi \in \left[0,t\right]$ and $k=0,1,2,\dots ,m+1$.

Therefore, Theorem 1 is proved. Similarly, when $u\left(x,y\right)\in C^{m+1}\left[0,R\right]\times \left[0,S\right]$, it can be proved in the same way based on Bernstein polynomials. The results show that the proposed algorithm is precise and efficacious for approximating unknown functions of three variables. □

5.2. Mathematical Example

The accuracy of the algorithm is verified by a mathematical example. It is represented by the following equation. The parameters in the mathematical example can be any values and have no realistic significance. The specific equation is as follows

$$0.1\frac{{\partial }^{2}u}{\partial t^2}+0.6D_t^{\alpha }u+0.01\left(\frac{{\partial }^{4}u}{\partial x^4}+\frac{{\partial }^{4}u}{\partial y^4}\right)+0.06\frac{{\partial }^{4}u}{\partial x^2\partial y^2}=F\left(x,y,t\right)$$

(35)

The boundary conditions are

$$\begin{array}{c}u\left(x,y,t\right)\left|{}_{x=0,1}\right.=\frac{\partial u\left(x,y,t\right)}{\partial x}\left|{}_{x=0,1}\right.=0\hfill \\ u\left(x,y,t\right)\left|{}_{y=0,2}\right.=\frac{\partial u\left(x,y,t\right)}{\partial y}\left|{}_{y=0,2}\right.=0\hfill \end{array}$$

(36)

where $\alpha =0.75$, $x\in [0,1]$, $y\in [0,2]$ and $t\in [0,1]$.

The exact solution is

$$u\left(x,y,t\right)=x^2{\left(1-x\right)}^2{y}^2{\left(2-y\right)}^2{t}^2$$

(37)

Substituting exact solution into Equation (35), mathematical equation is derived as

$$\begin{array}{c}F\left(x,y,t\right)=0.1\left(x^2{\left(1-x\right)}^2{y}^2{\left(2-y\right)}^2\right)+0.6\frac{\mathsf{\Gamma }\left(3\right)}{\mathsf{\Gamma }\left(3-\alpha \right)}{t}^{2-\alpha }\left(x^2{\left(1-x\right)}^2{y}^2{\left(2-y\right)}^2\right)\hfill \\ +0.01\left(24y^2{\left(2-y\right)}^2{t}^2+24x^2{\left(1-x\right)}^2{t}^2\right)+0.06\left(12x^2-12x+2\right)\left(12y^2-24y+8\right)t^2\hfill \end{array}$$

(38)

The Bernstein polynomials algorithm with the number of items $m=4$ is used to solve the proposed mathematical example. The numerical solution is $u_n\left(x,y,t\right)$. The absolute error $e_n$ is

$$e_n\left(x,y,t\right)=\left|u_n\left(x,y,t\right)-u\left(x,y,t\right)\right|$$

(39)

Figure 2a,b is the analytical and numerical solutions, respectively, at $t=1$. It can be seen that numerical solution and analytical solution are remarkably unanimous. Figure 2c shows absolute error and its minimum order of magnitude can reach ${10}^{-5}$. It can be proved that this method is high accuracy and its availability for simulating governing equations of fractional viscoelastic plates.

6. Numerical Analysis

The equation of the fractional quadrilateral clamped square plate is numerically analyzed. The plate displacement under external load is being investigated. The effect of different side lengths of the square plate on plate displacement is analyzed. Dynamic behavior of fractional plate with three boundary conditions is studied. In addition, the effect of thickness of the plate on the stress is also studied. In all subsequent studies, the time parameter is always considered as $t=1\mathrm{s}$.

Table 1. Geometric properties of viscoelastic plate materials [8].

Physical Quantity	Symbol	Value	Dimension
Fractional order	$\alpha$	0.75	1
Length	a	2	$\mathrm{m}$
Width	b	2	$\mathrm{m}$
Thickness	h	0.02	$\mathrm{m}$
Density of the plate	$\rho$	7850	$\mathrm{kg}·{\mathrm{m}}^3$
Poisson’s ratio	v	0.3	1
Young’s modulus	E	$2.1\times {10}^5$	$\mathrm{MPa}$
Damping coefficient	$\eta$	$5\times {10}^{-3}$	1

shows geometric characteristics of viscoelastic plate materials in dynamic analysis.

6.1. Influence of Different Simple Harmonic Loads on Plate Displacement

When taking the parameters in Table 1 for research, different simple harmonic loads are exerted to clamped square plate. Its form is $F=\vartheta sin\left(0.01t\right)$. Figure 3 is the dynamic response of a square plate. It can be found that the displacement of plate is symmetric at $x=y=1\mathrm{m}$ and reaches the maximum at the center point. With an increase in simple harmonic load coefficient $\vartheta$, the displacement of square plate also increases. Furthermore, when the load condition is $F=\vartheta cos\left(\frac{\pi }{4}t\right)$, the results of

Table 2. Displacement numerical solutions under simple harmonic loads.

$\mathit{u}\left(\mathit{x},\mathit{y},\mathit{t}\right)$	$\mathit{x}=0$	$\mathit{x}=0.5$	$\mathit{x}=1$	$\mathit{x}=1.5$	$\mathit{x}=2$
$F=0.002cos\left(\frac{\pi }{4}t\right)$	0	$0.01463$	$0.02522$	$0.01463$	0
$F=0.004cos\left(\frac{\pi }{4}t\right)$	0	$0.03052$	$0.05262$	$0.03052$	0
$F=0.006cos\left(\frac{\pi }{4}t\right)$	0	$0.04653$	$0.08024$	$0.04653$	0

are obtained through fixed width $y=1\mathrm{m}$. Displacement also increases with load and is symmetrical at the midpoint of the length.

6.2. Influence of Side 0 of the Plate on Plate Displacement

In this part, the influences of different side lengths of square plates on plate displacement are studied. Figure 4 is the change of displacement with side length of the plate under load $F=0.08sin\left(0.01t\right)$.

It can be seen that the displacement of the square plate increases with side length. This is consistent with the results of Reference [40], but this paper can use the proposed algorithm to solve the problem directly in the time domain. In Reference [40], a weak formal equation was constructed based on Hamilton’s principle and a four node rectangular plate element was used to discretize the region. When different external loads were used, the vibration response of the laminated plate was studied by refined plate theory finite element approach. The results indicate that Bernstein polynomials algorithm is an efficacious tool for solving fractional differential equations of three-dimensional plates. Displacement solution obtained by this algorithm has high accuracy.

6.3. Influence of Boundary Conditions on Plate Displacement

Figure 5 is about the influence of three boundary conditions on displacement of viscoelastic plates. CCCC, SSCC and SSSS denote completely simply supported plate, simply supported clamped plate and completely clamped plate in proper order. The expressions for the three of them are

$$\begin{array}{c}u\left(x,y,t\right)\left|{}_{x=0,a}\right.=\frac{\partial u\left(x,y,t\right)}{\partial x}\left|{}_{x=0,a}\right.=0\hfill \\ u\left(x,y,t\right)\left|{}_{y=0,b}\right.=\frac{\partial u\left(x,y,t\right)}{\partial y}\left|{}_{y=0,b}\right.=0\hfill \end{array}$$

(40)

$$\begin{array}{c}u\left(x,y,t\right)\left|{}_{x=0,a}\right.=\frac{{\partial }^{2}u\left(x,y,t\right)}{\partial x^2}\left|{}_{x=0,a}\right.=0\hfill \\ u\left(x,y,t\right)\left|{}_{y=0,b}\right.=\frac{\partial u\left(x,y,t\right)}{\partial y}\left|{}_{y=0,b}\right.=0\hfill \end{array}$$

(41)

$$\begin{array}{c}u\left(x,y,t\right)\left|{}_{x=0,a}\right.=\frac{{\partial }^{2}u\left(x,y,t\right)}{\partial x^2}\left|{}_{x=0,a}\right.=0\hfill \\ u\left(x,y,t\right)\left|{}_{y=0,b}\right.=\frac{{\partial }^{2}u\left(x,y,t\right)}{\partial y^2}\left|{}_{y=0,b}\right.=0\hfill \end{array}$$

(42)

As summarized from the figures, the square plate with the minimum constraints has the largest displacement, and the displacement of plate decreases with increase in constraints. This is consistent with the findings in Reference [40]. In this paper, the three-dimensional diagram is used to more intuitively show the change of displacement with constraints. This means that the displacement of the plate can be reduced by increasing the constraint of boundary conditions. Therefore, the proposed algorithm provides a theoretical basis for study of vibration analysis of viscoelastic plates.

6.4. Influence of Plate Thickness on Stress

The variation of stress with plate thickness will be analyzed. The expression for stress is

$$\sigma \left(x,y,t\right)=\eta D_t^{\alpha }\frac{{\partial }^{2}u\left(x,y,t\right)}{\partial x^2}$$

(43)

The stress distribution of the viscoelastic plate under external load $F=0.08sin\left(0.01t\right)$ is shown in Figure 6. The displacement of the plate will increase by the thickness of the plate. The conclusion is consistent with Reference [2]. By using the generalized multi-axial Maxwell model and Hamilton’s principle, they obtained the viscoelastic constitutive equation and proposed an effective equal geometric analysis formula for solving the nonlinear vibration problem of the viscoelastic plate. This verifies correctness of the numerical results. Therefore, the proposed algorithm is suitable for solving and studying stress problems of viscoelastic plates.

7. Conclusions

A new method for calculating fractional differential equations is presented in this paper. The differential operator matrices of Bernstein polynomials are used to approximate the displacement function directly in the time domain. The governing equations are transformed into algebraic equations for the solution. Error analysis and numerical results prove the correctness and effectiveness of the mentioned algorithm. In addition, the dynamic response of the three-dimensional viscoelastic plate is analyzed.

1. The displacement of viscoelastic plate is gained directly in the time domain by Bernstein polynomials algorithm. The unknown function is approximated by the operator matrix form of the mentioned polynomials.

2. The displacement of plate increases with the increase in the simple harmonic load coefficient and side length of square plate. The smaller the boundary condition constraints, the larger the viscoelastic plate displacement.

3. The effect of the plate thickness on stress is discussed. As the thickness of the plate increases, the stress value also increases.

4. The proposed algorithm can be applied to more complex differential governing equations, such as variable fractional viscoelastic plates and shells.