A Mixed Finite Volume Element Method for Time-Fractional Damping Beam Vibration Problem

Abstract

In this paper, the transverse vibration of a fractional viscoelastic beam is studied based on the fractional calculus, and the corresponding scheme of a viscoelastic beam is established by using the mixed finite volume element method. The stability and convergence of the algorithm are analyzed. Numerical examples demonstrate the effectiveness of the algorithm. Finally, the values of different parameter sets are tested, and the test results show that both the damping coefficient and the fractional derivative have significant effects on the model. The results of this paper can be used for the damping modeling of viscoelastic structures.

1. Introduction

According to the definition of fractional derivative, the fractional partial differential equation has more advantages than the integer equation in the study of some memory processes, genetic properties and heterogeneous materials. Therefore, such equations are widely applied to the fractal and dispersion in porous media [1,2], non-Newtonian fluid [3], the anomalous diffusion [4], image and signal processing [5], electric conduction [6], oil seepage and the piping of the boundary layer effect [7,8], etc.

Due to the particularity and complexity of viscoelastic materials, the traditional integer-order model cannot describe the viscoelastic properties well, so the fractional-order operator is introduced to construct the constitutive model of viscoelastic materials. Gement [9] first proposed the fractional derivative constitutive model of viscoelastic materials in 1936. In recent years, Demir et al. [10,11] studied the influence of the damping term modeled by a fractional derivative on the dynamic analysis of beams with viscoelastic properties under the action of harmonic external forces. Reza et al. [12] studied the forced vibration of a fractional-order viscoelastic beam and discretized the equations into a set of linear ordinary differential equations by the Galerkin method. The nonlocal fractional-order viscoelastic model of a nanobeam resting on a viscoelastic foundation was studied by Cajic et al. [13], where the solution of the fractional-order differential equation with two fractional parameters and retardation times was given. Liu et al. [14] proposed a simple and universal residual calculation method for the stochastic response behaviors of axially moving viscoelastic beams under random noise excitation and fractional constitutive relation. Yu et al. [15] analyzed the application of the fractional derivative in a damping vibration analysis of a viscoelastic single-mass system. Faraji et al. [16] analyzed the size-dependent geometrically nonlinear free vibrations of fractional viscoelastic simply supported and clamped-free nanobeams. The Galerkin scheme was used to simplify the fractional integral–partial differential governing equation into a time-dependent fractional ordinary differential equation, which was then solved by the predictive correction method. Yang et al. [17] investigated the stability of an axially moving beam constituted by fractional-order material under parametric resonances, where the governing equations of the beam transverse vibration were derived and then the multi-scale method was used to analyze the equation. Cao et al. [18] obtained a simple analytical expression for a free vibration analysis of non-uniform and non-homogenous beams under different boundary conditions by using the asymptotic perturbation approach. Liang et al. [19] utilized the Adomian decomposition method to solve a linear differential equation with an arbitrary fractional derivative order which can describe a fractionally damped beam structure. Sansit et al. [20] used the fractional finite element model to study the nonlocal response of Euler–Bernoulli beams under different loads and boundary conditions and provided analytical expressions and finite element solutions for the nonlocal continuum model of the Euler–Bernoulli beams. Stempin et al. [21] established a spatial fractional Timoshenko beam model with a functionally graded material effect and gave the experimental verification.

In this article, we consider the time-fractional damping beam vibration problem

$$\left\{\begin{array}{cc}\left(a\right){\mu }_0^c{D}_t^{\alpha }u+u_{tt}+a^2{u}_{xxxx}=g(x,t),\hfill & (x,t)\in \Omega \times (0,T],\hfill \\ \left(b\right)u(x,0)=\phi \left(x\right),u_t(x,0)=\psi \left(x\right),\hfill & x\in \Omega ,\hfill \\ \left(c\right)u(0,t)=u(L,t)=0,u_{xx}(0,t)=u_{xx}(L,t)=0,\hfill & t\in [0,T].\hfill \end{array}\right.$$

(1)

In the model, $u(x,t)$ represents the transverse vibration displacement of the beam, $\Omega =(0,L)$, $0<T<\infty$, $a=\sqrt{EI/\left(\rho A\right)}$, where $EI$, $\rho$ and A indicate the bending stiffness, density and cross-sectional area of the beam, respectively, and $\mu$(>0) is the damping coefficient. In this paper, we assume that the parameters $\rho$, A, $\mu$, E and I are constants. $g(x,t)$ represents the force exerted on the beam, and $\phi \left(x\right)$ and $\psi \left(x\right)$ represent the displacement and velocity of the beam at the initial time, respectively, and $\phi \left(x\right)\in C^2(\Omega )$, $\psi \left(x\right)\in C^1(\Omega )$, $g(x,t)\in L^1(0,T;L^2(\Omega ))$. ${}_0^c{D}_t^{\alpha }u$ is the Caputo fractional derivative, which is defined as [22]

$${}_0^c{D}_t^{\alpha }u(x,t)=\frac{1}{\Gamma (n-\alpha )}{\int }_0^t{(t-s)}^{n-\alpha -1}\frac{{\partial }^{n}u(x,s)}{\partial s^n}dsn-1<\alpha <n.$$

(2)

It is difficult to find the analytical solution of fractional partial differential equations, so the numerical method for solving fractional partial differential equations is widely considered. A great deal of work has been performed on the numerical solutions of fractional partial differential equations, and different methods have been discussed. Guo et al. [23] and Gao et al. [24] used the finite difference method to study fractional partial differential equations. Li et al. [25] used the compact finite difference method to solve the 2D time-fractional convection diffusion equation of groundwater pollution problems. Jin et al. [26] gave an error estimate for the fractional parabolic equation semi-discrete finite element method. Liu et al. [27] studied the H1-Galerkin mixed finite element method for the time-fractional reaction–diffusion equation. Su et al. [28] studied higher-order compact finite-volume schemes for two-dimensional multinomial time-fractional diffusion equations. Youssri [29,30] proposed the orthogonal ultraspherical operation matrix algorithm for the fractal and fractional Riccati equation with generalized Caputo derivatives and two Fibonacci operation matrix pseudo-spectral schemes for the nonlinear fractional Klein–Gordon equation. Sabir et al. [31] studied the fractional mathematical model of breast cancer immune–chemotherapy based on neural networks and designed a stochastic framework to solve the fractional differential model.

In 1995, the mixed finite volume element (MFVE) method was proposed by Russell [32]. This scheme is widely used in practical problems because it can solve two unknowns at the same time and keep the local conservation of a physical quantity. Up to now, there are few papers that discuss the numerical methods of time-fractional damping beam vibration problems. In this paper, we apply the MFVE method to the vibration problem (1).

The arrangement of the article is as follows. In the second part, the vibration equation of the damped beam (1) is transformed into second-order equations by introducing intermediate variables. Then, the spatial derivative term is discretized by the MFVE method, and the time-fractional derivative is approximated by the $L1$ interpolation formula to construct the MFVE scheme of (1). The third part gives some lemmas required for proof. The stability and convergence analysis for the MFVE scheme are analyzed in the fourth and fifth parts, respectively. In the sixth part, the accuracy of the scheme is verified by two numerical examples, and the parameter sets are tested to verify the influence of the parameters on the model.

2. Fully Discrete MFVE Scheme

We establish the approximate format of MFVE by introducing two intermediate variables

$$\begin{array}{c}\hfill v(x,t)=u_t(x,t),w(x,t)=-a{u}_{xx}(x,t),\end{array}$$

$v(x,t)$ and $w(x,t)$ have actual physical significance and represent the velocity and bending moment of the beam during transverse vibration, respectively, then (1) can be written as the following equation

$$\left\{\begin{array}{cc}\left(a\right){\mu }_0^c{D}_t^{\alpha -1}v+v_t-a{w}_{xx}=g(x,t)\hfill & (x,t)\in \Omega \times (0,T],\hfill \\ \left(b\right)w_t+a{v}_{xx}=0\hfill & (x,t)\in \Omega \times (0,T],\hfill \\ \left(c\right)v(x,0)=\psi \left(x\right),w(x,0)={\phi }^{\prime \prime }\left(x\right),\hfill & x\in \Omega ,\hfill \\ \left(d\right)v(0,t)=v(L,t)=0,w(0,t)=w(L,t)=0,\hfill & t\in [0,T].\hfill \end{array}\right.$$

(3)

Multiply Equation (3) (a) and (3) (b) by $\varphi \in H_0^1(\Omega )$, then integrate over $\Omega$ to obtain the weak form equivalent to (3): find $(v,w):[0,T]\to H_0^1(\Omega )\times H_0^1(\Omega )$, such that

$$\left\{\begin{array}{cc}\left(a\right)\mu {(}_0^c{D}_t^{\alpha -1}v,\varphi )+(v_t,\varphi )+(a{w}_x,{\varphi }_x)=(g,\varphi ),\hfill & \forall \varphi \in H_0^1(\Omega ),\hfill \\ \left(b\right)(w_t,\varphi )-(a{v}_x,{\varphi }_x)=0\hfill & \forall \varphi \in H_0^1(\Omega ),\hfill \\ \left(c\right)v(x,0)=\psi \left(x\right),w(x,0)={\phi }^{\prime \prime }\left(x\right),\hfill & x\in \Omega ,\hfill \\ \left(d\right)v(0,t)=v(L,t)=0,w(0,t)=w(L,t)=0,\hfill & t\in [0,T],\hfill \end{array}\right.$$

(4)

where $H_0^1(\Omega )=\left\{f\right|f\in H^1(\Omega ),f{|}_{\partial \Omega }=0\}.$

Next, we introduce the semi-discrete MFVE scheme of (1). Let $0=x_0<x_1<x_2<\cdots <x_N=L$ be the primal partition of $\overline{\Omega }$, the matching dual partition is $0=x_0<x_{\frac{1}{2}}<x_{\frac{3}{2}}<\cdots <x_{N-\frac{1}{2}}<x_N=L$, where $x_{i+\frac{1}{2}}=\frac{x_i+x_{i+1}}{2}$, $(i=0,1,2,\cdots ,N-1)$.

Give the primal partition ${\Re }_h=\{A_i=[x_i,x_{i+1}];i=0,1,2\cdots N-1\}$ of the region $\Omega$, the diameter of unit $A_i$ is $h_i=x_{i+1}-x_i$, let $h=\underset{0\le i\le N-1}{max}{h}_i$. Suppose ${\Re }_h$ is a quasi-uniform partition, that is, there exists some positive constant $\kappa$ such that $h_i\ge \kappa h$, ($i=1,2,\cdots ,N-1$). The dual subdivision is defined as ${\Re }_h^{*}=\{A_i^{*}=[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}];i=1,2,\cdots ,N-1\}$, where $A_0^{*}=[x_0,x_{\frac{1}{2}}],A_N^{*}=[x_{N-\frac{1}{2}},x_N]$, $A_i^{*}$ forms the dual interval of node i. As for the boundary nodes, its dual interval is revised accordingly.

Then, we define the finite element space

$$\begin{array}{cc}\hfill & U_h=\{u_h\in C\left(\overline{\Omega }\right),u_h{|}_A\in P_1,\forall A\in {\Re }_h\},\hfill \\ \hfill & V_h=\{v_h\in L^2(\Omega ),v_h{|}_{A^{*}}\in P_0,\forall A^{*}\in {\Re }_h^{*}\},\hfill \\ \hfill & U_{0h}=\{u_h\in U_h,u_h\left(0\right)=u_h\left(L\right)=0\},\hfill \\ \hfill & V_{0h}=\{v_h\in V_h,v_h\left(0\right)=v_h\left(L\right)=0\}.\hfill \end{array}$$

where $U_h$ represent the linear finite element space corresponding to the primal subdivision ${\Re }_h$, $V_h$ is the constant function space of corresponding dual subdivision ${\Re }_h^{*}$.

We define an interpolation operator ${\Pi }_h^{*}:U_h\to V_h$ by

$$\begin{array}{c}\hfill {\Pi }_h^{*}{w}_h=\sum _{i=1}^{N-1}{w}_h\left(x_i\right){\aleph }_{A_i^{*}},\forall w_h\in U_h.\end{array}$$

${\aleph }_{A_i^{*}}$ represents the eigenfunction on $A_i^{*}$, i.e.,

$${\aleph }_{A_i^{*}}=\left\{\begin{array}{c}1x\in A_i^{*},\hfill \\ 0x\notin A_i^{*}.\hfill \end{array}\right.$$

By integrating (3) over $A_i^{*}$, we obtain

$$\left\{\begin{array}{c}\mu {\int }_{A_i^{*}}{}_0^c{D}_t^{\alpha -1}v(x,t)dx+{\int }_{A_i^{*}}{v}_t(x,t)dx-a[w_x(x_{i+\frac{1}{2}},t)-w_x(x_{i-\frac{1}{2}},t)]={\int }_{A_i^{*}}g(x,t)dx,\hfill \\ {\int }_{A_i^{*}}{w}_t(x,t)dx+a[v_x(x_{i+\frac{1}{2}},t)-v_x(x_{i-\frac{1}{2}},t)]=0.\hfill \end{array}\right.$$

(5)

Add all the elements together and notice that: for all $w\in L^2(\Omega )$ and ${\varphi }_h\in U_h$

$$\begin{array}{c}{\displaystyle (w,{\Pi }_h^{*}{\varphi }_h)=\sum _{i=1}^{N-1}{\varphi }_h\left(x_i\right){\int }_{A_i^{*}}wdx.}\hfill \end{array}$$

Define

$$\begin{array}{c}B(w,{\Pi }_h^{*}{\varphi }_h)=-{\displaystyle \sum _{i=1}^{N-1}}{\varphi }_h\left(x_i\right)a[w_x(x_{i+\frac{1}{2}},t)-w_x(x_{i-\frac{1}{2}},t)].\hfill \end{array}$$

Then, we have

$$\left\{\begin{array}{c}\mu {(}_0^c{D}_t^{\alpha -1}v,{\Pi }_h^{*}{\varphi }_h)+(v_t,{\Pi }_h^{*}{\varphi }_h)+B(w,{\Pi }_h^{*}{\varphi }_h)=(g,{\Pi }_h^{*}{\varphi }_h),\forall {\varphi }_h\in U_{0h},\hfill \\ (w_t,{\Pi }_h^{*}{\varphi }_h)-B(v,{\Pi }_h^{*}{\varphi }_h)=0,\forall {\varphi }_h\in U_{0h}.\hfill \end{array}\right.$$

(6)

Then, the corresponding semi-discrete MFVE scheme of problem (1) is: find $(v_h,w_h)\in U_{0h}\times U_{0h}$, such that

$$\left\{\begin{array}{c}\left(a\right)\mu {(}_0^c{D}_t^{\alpha -1}{v}_h,{\Pi }_h^{*}{\varphi }_h)+(v_{ht},{\Pi }_h^{*}{\varphi }_h)+B(w_h,{\Pi }_h^{*}{\varphi }_h)=(g,{\Pi }_h^{*}{\varphi }_h),\forall {\varphi }_h\in U_{0h},\hfill \\ \left(b\right)(w_{ht},{\Pi }_h^{*}{\varphi }_h)-B(v_h,{\Pi }_h^{*}{\varphi }_h)=0,\forall {\varphi }_h\in U_{0h}.\hfill \end{array}\right.$$

(7)

Now, let $0=t_0<t_1<\cdots <t_M=T$ be the subdivision of time interval $[0,T]$ with step length $\tau =T/M$, $t_n=n\tau ,n=0,1,\cdots ,M$. For a smooth function $\varphi$ on $[0,T]$, we denote ${\varphi }^n=\varphi \left(t_n\right),{\partial }_t{\varphi }^n=({\varphi }^n-{\varphi }^{n-1})/\tau$. We can use $L1$-formula [33] to approximate the time-fractional derivative at $t=t_n$ as follows

$$\begin{array}{cc}\hfill {}_0^c{D}_t^{\alpha }{u}^n={}_0^c{D}_t^{\alpha -1}{v}^n& =\frac{{\tau }^{1-\alpha }}{\Gamma (3-\alpha )}[v^n-\sum _{k=1}^{n-1}(a_{n-k-1}^{(\alpha -1)}-a_{n-k}^{(\alpha -1)})v^k-a_{n-1}^{(\alpha -1)}{v}^0]+R.\hfill \end{array}$$

where $a_l^{\left(\alpha \right)}={(l+1)}^{(1-\alpha )}-l^{(1-\alpha )},\left|R\right|\le \frac{1}{2\Gamma (2-\alpha )}[\frac{1}{4}+\frac{\alpha -1}{(2-\alpha )(3-\alpha )}\underset{t_0\le t\le t_n}{max}|v^{\prime \prime }\left(t\right)\left|{\tau }^{(3-\alpha )}\right].$

Denote $\delta =\Gamma (3-\alpha ){\tau }^{\alpha -1}$, $D_{\tau }^{\alpha -1}{v}^n=\frac{1}{\delta }[v^n-{\displaystyle \sum _{k=1}^{n-1}}(a_{n-k-1}^{(\alpha -1)}-a_{n-k}^{(\alpha -1)})v^k-a_{n-1}^{(\alpha -1)}{v}^0]$, we have

$$\begin{array}{c}\hfill {}_0^c{D}_t^{\alpha -1}v\left(t_n\right)=D_{\tau }^{\alpha -1}{v}^n+R.\end{array}$$

Then, we obtain the fully discrete MFVE scheme to find $\{v_h^n,w_h^n\}\in U_{0h}\times U_{0h}$, $n=1,2,\cdots ,M$, such that

$$\left\{\begin{array}{c}\left(a\right)\mu (D_{\tau }^{\alpha -1}{v}_h^n,{\Pi }_h^{*}{\varphi }_h)+({\partial }_t{v}_h^n,{\Pi }_h^{*}{\varphi }_h)+B(w_h^n,{\Pi }_h^{*}{\varphi }_h)=(g^n,{\Pi }_h^{*}{\varphi }_h),\forall {\varphi }_h\in U_{0h},\hfill \\ \left(b\right)({\partial }_t{w}_h^n,{\Pi }_h^{*}{\varphi }_h)-B(v_h^n,{\Pi }_h^{*}{\varphi }_h)=0,\forall {\varphi }_h\in U_{0h}.\hfill \end{array}\right.$$

(8)

3. Some Lemmas

In this part, we will give some necessary lemmas. Let $u_h\in U_{0h}$, we define the following norms

$$\begin{array}{cc}\hfill & \parallel u_h{\parallel }_{0,h}^2=\sum _{i=1}^{N-1}h{u}_i^2,\hfill \\ \hfill & \parallel u_h{\parallel }_{1,h}^2=\parallel u_h{\parallel }_{0,h}^2+{\left|u_h\right|}_{1,h}^2,\hfill \\ \hfill & |u_h{|}_{1,h}^2=\sum _{i=1}^{N}h{(\frac{u_i-u_{i-1}}{h})}^2.\hfill \end{array}$$

Lemma 1

([34]). On $U_h$, the norms $\mid ·{\mid }_{1,h}$ and $\mid ·{\mid }_1$;${\Vert ·\Vert }_{0,h}$ and $\Vert ·\Vert$;${\Vert ·\Vert }_{1,h}$ and ${\Vert ·\Vert }_1$ are equivalent, respectively. That is,

$$\begin{array}{cc}\hfill & C_1\parallel v_h{\parallel }_{0,h}\le \parallel v_h\parallel \le C_2{\parallel v_h\parallel }_{0,h},\forall v_h\in U_{0h},\hfill \\ \hfill & C_3\parallel v_h{\parallel }_{1,h}\le \parallel v_h{\parallel }_1\le C_4{\parallel v_h\parallel }_{1,h},\forall v_h\in U_{0h}.\hfill \end{array}$$

where $C_1,\cdots ,C_4$ are positive constants independent of $U_h$.

Lemma 2

([35]). For the bilinear form $B(·,{\Pi }_h^{*}·)$, it holds that

$$\begin{array}{cc}\hfill & B(v_h,{\Pi }_h^{*}{w}_h)=B(w_h,{\Pi }_h^{*}{v}_h),\forall v_h,w_h\in U_h,\hfill \\ \hfill & B(v_h,{\Pi }_h^{*}{v}_h)\ge \frac{1}{2}{\left|v_h\right|}_1^2,\forall v_h\in U_{0h},\hfill \end{array}$$

Lemma 3

([35]). $(·,{\Pi }_h^{*}·)$ satisfies

$$\begin{array}{cc}\hfill & (v_h,{\Pi }_h^{*}{w}_h)=(w_h,{\Pi }_h^{*}{v}_h),\forall v_h,w_h\in U_h,\hfill \\ \hfill & (v_h,{\Pi }_h^{*}{v}_h)\ge \frac{1}{4}{\Vert v_h\Vert }^2,\forall v_h\in U_h,\hfill \\ \hfill & (w,{\Pi }_h^{*}{v}_h)\le 3\Vert w\Vert ·\Vert v_h\Vert ,\forall w\in H^1(\Omega ),\forall v_h\in U_h.\hfill \end{array}$$

4. Stability Analysis for Fully Discrete MFVE Scheme

Theorem 1

([36]). For the scheme $\left(8\right)$, the following stable inequality holds, for sufficiently small τ

$$\begin{array}{c}\begin{array}{c}\hfill \Vert v_h^n{\Vert }^2+\Vert u_h^n{\Vert }^2\le C[\frac{t_n^{2-\alpha }}{\Gamma (3-\alpha )}{\Vert \psi \Vert }^2+\Vert {\phi }^{\prime \prime }{\Vert }^2+\Gamma (2-\alpha )t_n^{\alpha -1}\tau \sum _{k=1}^n\Vert g^k{\Vert }^2],\end{array}\hfill \end{array}$$

where $C>0$ is a constant free of two mesh parameters τ and h.

Proof. 

Take ${\varphi }_h=v_h$ in (8)(a), ${\varphi }_h=w_h$ in (8)(b), we have

$$\begin{array}{c}\hfill ({\partial }_t{w}_h^n,{\Pi }_h^{*}{w}_h^n)+({\partial }_t{v}_h^n,{\Pi }_h^{*}{v}_h^n)+\mu (D_{\tau }^{\alpha -1}{v}_h^n,{\Pi }_h^{*}{v}_h^n)=(g^n,{\Pi }_h^{*}{v}_h^n).\end{array}$$

(9)

Note the fact that

$$\begin{array}{cc}\hfill & ({\partial }_t{w}_h^n,{\Pi }_h^{*}{w}_h^n)\ge \frac{1}{2\tau }[(w_h^n,{\Pi }_h^{*}{w}_h^n)-(w_h^{n-1},{\Pi }_h^{*}{w}_h^{n-1})],\hfill \end{array}$$

(10)

$$\begin{array}{cc}\hfill & ({\partial }_t{v}_h^n,{\Pi }_h^{*}{v}_h^n)\ge \frac{1}{2\tau }[(v_h^n,{\Pi }_h^{*}{v}_h^n)-(v_h^{n-1},{\Pi }_h^{*}{v}_h^{n-1})],\hfill \\ \hfill & (D_{\tau }^{(\alpha -1)}{v}_h^n,{\Pi }_h^{*}{v}_h^n)\ge \frac{1}{\delta }[\frac{1}{4}\Vert v_h^n{\Vert }^2-\sum _{k=1}^{n-1}(a_{n-k-1}^{(\alpha -1)}-a_{n-k}^{(\alpha -1)})(v_h^k,{\Pi }_h^{*}{v}_h^n)-a_{n-1}^{(\alpha -1)}(v_h^0,{\Pi }_h^{*}{v}_h^n)]\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill & \ge \frac{1}{\delta }[\frac{1}{4}\Vert v_h^n{\Vert }^2-\sum _{k=1}^{n-1}(a_{n-k-1}^{(\alpha -1)}-a_{n-k}^{(\alpha -1)})(\frac{1}{4}\Vert v_h^k{\Vert }^2+9\Vert v_h^n{\Vert }^2)-a_{n-1}^{(\alpha -1)}(\frac{1}{4}\Vert v_h^0{\Vert }^2+9\Vert v_h^n{\Vert }^2].\hfill \end{array}$$

(12)

Using (10)–(12), we rewrite (9) as

$$\begin{array}{cc}\hfill & \frac{1}{2\tau }[(w_h^n,{\Pi }_h^{*}{w}_h^n)-(w_h^{n-1},{\Pi }_h^{*}{w}_h^{n-1})]+\frac{1}{2\tau }[(v_h^n,{\Pi }_h^{*}{v}_h^n)-(v_h^{n-1},{\Pi }_h^{*}{v}_h^{n-1})]+\frac{\mu }{4\delta }{\Vert v_h^n\Vert }^2\hfill \\ \hfill & \le \frac{\mu }{\delta }\sum _{k=1}^{n-1}(a_{n-k-1}^{(\alpha -1)}-a_{n-k}^{(\alpha -1)})(\frac{1}{4}\Vert v_h^k{\Vert }^2+9\Vert v_h^n{\Vert }^2)+\frac{\mu }{\delta }{a}_{n-1}^{(\alpha -1)}(\frac{1}{4}\Vert v_h^0{\Vert }^2+9\Vert v_h^n{\Vert }^2)+(g^n,{\Pi }_h^{*}{v}_h^n).\hfill \end{array}$$

which leads to

$$\begin{array}{cc}\hfill & \frac{1}{2\tau }(w_h^n,{\Pi }_h^{*}{w}_h^n)+\frac{1}{2\tau }(v_h^n,{\Pi }_h^{*}{v}_h^n)+\frac{\mu }{4\delta }\sum _{k=1}^n{a}_{n-k}^{(\alpha -1)}{\Vert v_h^k\Vert }^2\hfill \\ \hfill & \le \frac{1}{2\tau }(w_h^{n-1},{\Pi }_h^{*}{w}_h^{n-1})+\frac{1}{2\tau }(v_h^{n-1},{\Pi }_h^{*}{v}_h^{n-1})+\frac{\mu }{4\delta }\sum _{k=1}^{n-1}{a}_{n-k-1}^{(\alpha -1)}{\Vert v_h^k\Vert }^2\hfill \\ \hfill & +\frac{9\mu }{\delta }\Vert v_h^n{\Vert }^2+\frac{\mu }{4\delta }{a}_{n-1}^{(\alpha -1)}{\Vert v_h^0\Vert }^2+(g^n,{\Pi }_h^{*}{v}_h^n).\hfill \end{array}$$

(13)

Denote

$$\begin{array}{c}\hfill F^n=(w_h^n,{\Pi }_h^{*}{w}_h^n)+(v_h^n,{\Pi }_h^{*}{v}_h^n)+\frac{\tau \mu }{2\delta }\sum _{k=1}^n{a}_{n-k}^{(\alpha -1)}{\Vert v_h^k\Vert }^2.\end{array}$$

Multiplying $\left(13\right)$ by $2\tau$ and using Young inequality, we obtain

$$\begin{array}{cc}\hfill F^n& \le F^{n-1}+\frac{\tau \mu }{2\delta }{a}_{n-1}^{(\alpha -1)}\Vert v_h^0{\Vert }^2+\frac{18\tau \mu }{\delta }{\Vert v_h^n\Vert }^2+2\tau \mid (g^n,{\Pi }_h^{*}{v}_h^n)\mid \hfill \\ \hfill & \le F^0+\frac{\tau \mu }{2\delta }\sum _{k=1}^n{a}_{k-1}^{(\alpha -1)}\Vert v_h^0{\Vert }^2+\frac{18\tau \mu }{\delta }\sum _{k=1}^n{\Vert v_h^k\Vert }^2+2\tau \sum _{k=1}^n\mid (g^k,{\Pi }_h^{*}{v}_h^k)\mid \hfill \\ \hfill & \le \Vert v_h^0{\Vert }^2+\Vert w_h^0{\Vert }^2+\frac{\tau \mu }{2\delta }\sum _{k=1}^n{a}_{k-1}^{(\alpha -1)}\Vert v_h^0{\Vert }^2+\frac{18\tau \mu }{\delta }\sum _{k=1}^{n-1}\Vert v_h^k{\Vert }^2+\frac{18\tau \mu }{\delta }{\Vert v_h^n\Vert }^2\hfill \\ \hfill & +2\tau \sum _{k=1}^n(\frac{\delta }{\mu a_{n-k}^{(\alpha -1)}}\Vert g^k{\Vert }^2+\frac{\mu a_{n-k}^{(\alpha -1)}}{4\delta }\Vert v_h^k{\Vert }^2).\hfill \end{array}$$

Choosing $\tau$ to satisfy $\frac{18\tau \mu }{\delta }<\frac{1}{4}$, and using the discrete Gronwall’s lemma and Lemma 3, we have

$$\begin{array}{c}\hfill \Vert v_h^n{\Vert }^2+\Vert w_h^n{\Vert }^2\le C[\Vert v_h^0{\Vert }^2+\Vert w_h^0{\Vert }^2+\frac{\tau \mu }{2\delta }\sum _{k=1}^n{a}_{k-1}^{(\alpha -1)}\Vert v_h^0{\Vert }^2+\sum _{k=1}^n\frac{2\tau \delta }{\mu a_{n-k}^{(\alpha -1)}}\Vert g^k{\Vert }^2].\end{array}$$

(14)

Noting that $a_{n-k}^{(\alpha -1)}\ge (2-\alpha ){(n-k+1)}^{1-\alpha }\ge (2-\alpha )n^{1-\alpha }$ when $1\le k\le n$, we derive

$$\begin{array}{c}\hfill \frac{a_{n-k}^{(\alpha -1)}}{\delta }=\frac{a_{n-k}^{(\alpha -1)}}{\Gamma (3-\alpha ){\tau }^{\alpha -1}}\ge \frac{(2-\alpha )n^{1-\alpha }}{\Gamma (3-\alpha ){\tau }^{\alpha -1}}=\frac{t_n^{1-\alpha }}{\Gamma (2-\alpha )}.\end{array}$$

Thus, we find that

$$\begin{array}{c}\hfill \sum _{k=1}^n\frac{2\tau \delta }{\mu a_{n-k}^{(\alpha -1)}}\Vert g^k{\Vert }^2\le 2\Gamma (2-\alpha )t_n^{\alpha -1}\frac{\tau }{\mu }\sum _{k=1}^n{\Vert g^k\Vert }^2.\end{array}$$

(15)

From ${\displaystyle \sum _{k=1}^n}{a}_{k-1}^{(\alpha -1)}=n^{2-\alpha }$, the following estimate holds

$$\begin{array}{c}\hfill \frac{\tau \mu }{2\delta }\sum _{k=1}^n{a}_{k-1}^{(\alpha -1)}\Vert v_h^0{\Vert }^2=\frac{\mu t_n^{2-\alpha }}{2\Gamma (3-\alpha )}{\Vert v_h^0\Vert }^2.\end{array}$$

(16)

Collecting from $\left(15\right),\left(16\right)$ to $\left(14\right)$, it holds that

$$\begin{array}{c}\hfill \Vert v_h^n{\Vert }^2+\Vert w_h^n{\Vert }^2\le C[\frac{t_n^{2-\alpha }}{\Gamma (3-\alpha )}{\Vert \psi \Vert }^2+\Vert {\phi }^{\prime \prime }{\Vert }^2+\Gamma (2-\alpha )t_n^{\alpha -1}\tau \sum _{k=1}^n\Vert g^k{\Vert }^2].\end{array}$$

We complete the proof. □

5. Convergence Analysis for Fully Discrete MFVE Scheme

In this section, we will estimate the error for the fully discrete MFVE scheme. Firstly, we introduce the MFVE elliptic projection to analyze the error of the scheme: find $({\tilde{v}}_h^n,{\tilde{w}}_h^n):[0,T]\to U_{0h}\times U_{0h}$, satisfies:

$$\left\{\begin{array}{cc}B({\tilde{w}}_h^n-w^n,{\Pi }_h^{*}{\varphi }_h)=0,\hfill & \forall {\varphi }_h\in U_{0h},\hfill \\ B({\tilde{v}}_h^n-v^n,{\Pi }_h^{*}{\varphi }_h)=-({\tilde{w}}_h^n-w^n,{\Pi }_h^{*}{\varphi }_h),\hfill & \forall {\varphi }_h\in U_{0h}.\hfill \end{array}\right.$$

(17)

On the condition that ${\Re }_h$ is C-uniform subdivision, the MFVE elliptic projection is unique and satisfies [34,35]:

$$\begin{array}{c}\mid v^n-{\tilde{v}}_h^n{\mid }_1+\Vert w^n-{\tilde{w}}_h^n\Vert \le Ch(\Vert v^n{\Vert }_3+\Vert w^n{\Vert }_2),\hfill \\ \Vert v^n-{\tilde{v}}_h^n\Vert \le Ch(\Vert v^n{\Vert }_3+\Vert w^n{\Vert }_2).\hfill \end{array}$$

(18)

We split the errors

$$\begin{array}{cc}\hfill v_h^n-v^n& =(v_h^n-{\tilde{v}}_h^n)+({\tilde{v}}_h^n-v^n)={\xi }^n+{\eta }^n,\hfill \\ \hfill w_h^n-w^n& =(w_h^n-{\tilde{w}}_h^n)+({\tilde{w}}_h^n-w^n)={\rho }^n+{\theta }^n.\hfill \end{array}$$

Differentiating (18) on t, the estimate of ${\eta }_t^n$ and ${\theta }_t^n$ can be obtained by the same method as employed in Refs. [34,35].

$$\begin{array}{cc}\hfill \Vert {\eta }_t^n\Vert & =\Vert {\tilde{v}}_{ht}^n-v_t^n\Vert \le Ch(\Vert v_t^n{\Vert }_3+\Vert w_t^n{\Vert }_2),\hfill \\ \hfill \Vert {\theta }_t^n\Vert & =\Vert {\tilde{w}}_{ht}^n-w_t^n\Vert \le Ch(\Vert v_t^n{\Vert }_3+\Vert w_t^n{\Vert }_2).\hfill \end{array}$$

From (6), (8) and combining with the elliptic projection (17), we obtain the following error equations

$$\left\{\begin{array}{c}\left(a\right)({\partial }_t{\rho }^n,{\Pi }_h^{*}{\varphi }_h)+({\partial }_t{\theta }^n,{\Pi }_h^{*}{\varphi }_h)+B({\xi }^n,{\Pi }_h^{*}{\varphi }_h)-({\theta }^n,{\Pi }_h^{*}{\varphi }_h)\hfill \\ =(R_w^n,{\Pi }_h^{*}{\varphi }_h),\forall {\varphi }_h\in U_{0h},\hfill \\ \left(b\right)\mu (D_{\tau }^{\alpha -1}{\xi }^n,{\Pi }_h^{*}{\varphi }_h)+\mu (D_{\tau }^{\alpha -1}{\eta }^n,{\Pi }_h^{*}{\varphi }_h)+({\partial }_t{\xi }^n,{\Pi }_h^{*}{\varphi }_h)+({\partial }_t{\eta }^n,{\Pi }_h^{*}{\varphi }_h)-B({\rho }^n,{\Pi }_h^{*}{\varphi }_h)\hfill \\ =(R_v^n,{\Pi }_h^{*}{\varphi }_h)+\mu (R^n,{\Pi }_h^{*}{\varphi }_h),\forall {\varphi }_h\in U_{0h},\hfill \end{array}\right.$$

(19)

where

$$\begin{array}{c}\hfill R_w^n={\partial }_t{w}^n-w_t^n=\frac{1}{\tau }{\int }_{t_n-1}^{t_n}(t_{n-1}-s)w_{tt}\left(s\right)ds,R_v^n={\partial }_t{v}^n-v_t^n=\frac{1}{\tau }{\int }_{t_n-1}^{t_n}(t_{n-1}-s)v_{tt}\left(s\right)ds.\end{array}$$

Theorem 2.

On the condition that ${\Re }_h$ is quasi-uniform subdivision, if $(v^n,w^n)$ is a solution to problem (3) and $v,w$ satisfies the required regularity condition. Then, the solution $(v_h^n,w_h^n)\in U_{0h}\times U_{0h}$ of the fully discrete MFVE scheme (8) converges to $(v^n,w^n)$, and there exists a positive constant C which does not depend on the subdivision of ${\Re }_h$ meeting the following estimation

$$\begin{array}{cc}\hfill & \underset{0\le n\le T/N}{max}\Vert v_h^n-v^n\Vert +\underset{0\le n\le T/N}{max}\Vert w_h^n-w^n\Vert \le C(h+\tau ),\hfill \\ \hfill & \underset{0\le n\le T/N}{max}\mid v_h^n-v^n{\mid }_1\le C(h+\tau ).\hfill \end{array}$$

Proof. 

Choosing ${\varphi }_h={\rho }^n$ in (19) (a) and ${\varphi }_h={\xi }^n$ in (19) (b) to obtain

$$\begin{array}{cc}\hfill \mu & (D_{\tau }^{\alpha -1}{\xi }^n,{\Pi }_h^{*}{\xi }^n)+\mu (D_{\tau }^{\alpha -1}{\eta }^n,{\Pi }_h^{*}{\xi }^n)+({\partial }_t{\xi }^n,{\Pi }_h^{*}{\xi }^n)+({\partial }_t{\eta }^n,{\Pi }_h^{*}{\xi }^n)+({\partial }_t{\rho }^n,{\Pi }_h^{*}{\rho }^n)+({\partial }_t{\theta }^n,{\Pi }_h^{*}{\rho }^n)\hfill \\ \hfill & =({\theta }^n,{\Pi }_h^{*}{\rho }^n)+(R_v^n,{\Pi }_h^{*}{\xi }^n)+(R_w^n,{\Pi }_h^{*}{\rho }^n)+\mu (R^n,{\Pi }_h^{*}{\xi }^n).\hfill \end{array}$$

(20)

Note the fact that

$$\begin{array}{cc}\hfill & ({\partial }_t{\xi }^n,{\Pi }^{*}{\xi }^n)\ge \frac{1}{2\tau }[({\xi }^n,{\Pi }_h^{*}{\xi }^n)-({\xi }^{n-1},{\Pi }_h^{*}{\xi }^{n-1})],\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill & ({\partial }_t{\rho }^n,{\Pi }^{*}{\rho }^n)\ge \frac{1}{2\tau }[({\rho }^n,{\Pi }_h^{*}{\rho }^n)-({\rho }^{n-1},{\Pi }_h^{*}{\rho }^{n-1})].\hfill \end{array}$$

(22)

Substituting (21) and (22) into $\left(20\right)$ and using $L1$-formula, we can obtain

$$\begin{array}{cc}\hfill & \frac{1}{2\tau }[({\xi }^n,{\Pi }_h^{*}{\xi }^n)-({\xi }^{n-1},{\Pi }_h^{*}{\xi }^{n-1})]+\frac{1}{2\tau }[({\rho }^n,{\Pi }_h^{*}{\rho }^n)-({\rho }^{n-1},{\Pi }_h^{*}{\rho }^{n-1})]+\frac{\mu }{\delta }({\xi }^n,{\Pi }_h^{*}{\xi }^n)\hfill \\ \hfill & \le \frac{\mu }{\delta }\sum _{k=1}^{n-1}(a_{n-k-1}^{(\alpha -1)}-a_{n-k}^{(\alpha -1)})({\xi }^k,{\Pi }_h^{*}{\xi }^n)+\frac{\mu }{\delta }{a}_{n-1}^{(\alpha -1)}({\xi }^0,{\Pi }_h^{*}{\xi }^n)-[({\partial }_t{\eta }^n,{\Pi }_h^{*}{\xi }^n)+\mu (D_{\tau }^{\alpha -1}{\eta }^n,{\Pi }_h^{*}{\xi }^n)\hfill \\ \hfill & +({\partial }_t{\theta }^n,{\Pi }_h^{*}{\rho }^n)]+({\theta }^n,{\Pi }_h^{*}{\rho }^n)+(R_v^n,{\Pi }_h^{*}{\xi }^n)+(R_w^n,{\Pi }_h^{*}{\rho }^n)+\mu (R^n,{\Pi }_h^{*}{\xi }^n).\hfill \end{array}$$

(23)

Multiplying $\left(23\right)$ by $2\tau$, then using Lemma 3 and Young inequality, we obtain

$$\begin{array}{cc}\hfill & [({\xi }^n,{\Pi }_h^{*}{\xi }^n)-({\xi }^{n-1},{\Pi }_h^{*}{\xi }^{n-1})]+[({\rho }^n,{\Pi }_h^{*}{\rho }^n)-({\rho }^{n-1},{\Pi }_h^{*}{\rho }^{n-1})]+\frac{\tau \mu }{2\delta }{\Vert {\xi }^n\Vert }^2\hfill \\ \hfill & \le \frac{2\tau \mu }{\delta }\sum _{k=1}^{n-1}(a_{n-k-1}^{(\alpha -1)}-a_{n-k}^{(\alpha -1)})(\frac{1}{4}\Vert {\xi }^k{\Vert }^2+9\Vert {\xi }^n{\Vert }^2)+\frac{2\tau \mu }{\delta }{a}_{n-1}^{(\alpha -1)}(\frac{1}{4}\Vert {\xi }^0{\Vert }^2+9\Vert {\xi }^n{\Vert }^2)\hfill \\ \hfill & -2\tau [({\partial }_t{\eta }^n,{\Pi }_h^{*}{\xi }^n)+\mu (D_{\tau }^{\alpha -1}{\eta }^n,{\Pi }_h^{*}{\xi }^n)+({\partial }_t{\theta }^n,{\Pi }_h^{*}{\rho }^n)]+2\tau ({\theta }^n,{\Pi }_h^{*}{\rho }^n)\hfill \\ \hfill & +2\tau (R_v^n,{\Pi }_h^{*}{\xi }^n)+2\tau (R_w^n,{\Pi }_h^{*}{\rho }^n)+2\tau \mu (R^n,{\Pi }_h^{*}{\xi }^n).\hfill \end{array}$$

which leads to

$$\begin{array}{cc}\hfill & ({\xi }^n,{\Pi }_h^{*}{\xi }^n)+({\rho }^n,{\Pi }_h^{*}{\rho }^n)+\frac{\tau \mu }{2\delta }\sum _{k=1}^n{a}_{n-k}^{(\alpha -1)}{\Vert {\xi }^k\Vert }^2\hfill \\ \hfill & \le ({\xi }^{n-1},{\Pi }_h^{*}{\xi }^{n-1})+({\rho }^{n-1},{\Pi }_h^{*}{\rho }^{n-1})+\frac{\tau \mu }{2\delta }\sum _{k=1}^{n-1}{a}_{n-k-1}^{(\alpha -1)}\Vert {\xi }^k{\Vert }^2+\frac{\tau \mu }{2\delta }{a}_{n-1}^{(\alpha -1)}\Vert {\xi }^0{\Vert }^2+\frac{18\tau \mu }{\delta }{\Vert {\xi }^n\Vert }^2\hfill \\ \hfill & -2\tau [({\partial }_t{\eta }^n,{\Pi }_h^{*}{\xi }^n)+\mu (D_{\tau }^{\alpha -1}{\eta }^n,{\Pi }_h^{*}{\xi }^n)+({\partial }_t{\theta }^n,{\Pi }_h^{*}{\rho }^n)]+2\tau ({\theta }^n,{\Pi }_h^{*}{\rho }^n)\hfill \\ \hfill & +2\tau (R_v^n,{\Pi }_h^{*}{\xi }^n)+2\tau (R_w^n,{\Pi }_h^{*}{\rho }^n)+2\tau \mu (R^n,{\Pi }_h^{*}{\xi }^n).\hfill \end{array}$$

(24)

Choosing $G^n=({\xi }^n,{\Pi }_h^{*}{\xi }^n)+({\rho }^n,{\Pi }_h^{*}{\rho }^n)+\frac{\tau \mu }{2\delta }{\displaystyle \sum _{k=1}^n}{a}_{n-k}^{(\alpha -1)}{\Vert {\xi }^k\Vert }^2$, and substituting into $\left(24\right)$, we have

$$\begin{array}{cc}\hfill G^n& \le G^{n-1}+\frac{\tau \mu }{2\delta }{a}_{n-1}^{(\alpha -1)}\Vert {\xi }^0{\Vert }^2+\frac{18\tau \mu }{\delta }\Vert {\xi }^n{\Vert }^2-2\tau [({\partial }_t{\eta }^n,{\Pi }_h^{*}{\xi }^n)+\mu (D_{\tau }^{\alpha -1}{\eta }^n,{\Pi }_h^{*}{\xi }^n)\hfill \\ \hfill & +({\partial }_t{\theta }^n,{\Pi }_h^{*}{\rho }^n)]+2\tau ({\theta }^n,{\Pi }_h^{*}{\rho }^n)+2\tau (R_v^n,{\Pi }_h^{*}{\xi }^n)+2\tau (R_w^n,{\Pi }_h^{*}{\rho }^n)+2\tau \mu (R^n,{\Pi }_h^{*}{\xi }^n)\hfill \\ \hfill & \le G^0+\frac{\tau \mu }{2\delta }\sum _{k=1}^n{a}_{k-1}^{(\alpha -1)}\Vert {\xi }^0{\Vert }^2+\sum _{k=1}^n\frac{18\tau \mu }{\delta }{\Vert {\xi }^k\Vert }^2+\sum _{j=1}^7{M}_j.\hfill \end{array}$$

(25)

To analyze the right-hand terms of (25) in turn, we have

$$\begin{array}{cc}\hfill M_1=-2\tau \sum _{k=1}^n({\partial }_t{\eta }^k,{\xi }^k)& \le \sum _{k=1}^n(\frac{18\delta }{\tau \mu a_{n-k}^{(\alpha -1)}}\Vert {\int }_{t_{k-1}}^{t_k}{\eta }_t{ds\Vert }^2+\frac{\tau \mu a_{n-k}^{(\alpha -1)}}{8\delta }\Vert {\xi }^k{\Vert }^2)\hfill \\ \hfill & \le \sum _{k=1}^n(\frac{18\delta }{\mu a_{n-k}^{(\alpha -1)}}{\int }_{t_{k-1}}^{t_k}\Vert {\eta }_t{\Vert }^{2}ds+\frac{\tau \mu a_{n-k}^{(\alpha -1)}}{8\delta }\Vert {\xi }^k{\Vert }^2).\hfill \end{array}$$

$$\begin{array}{c}\hfill M_2=-2\tau \mu \sum _{k=1}^n(D_{\tau }^{\alpha -1}{\eta }^k,{\xi }^k)\le \sum _{k=1}^n\frac{18\tau \mu \delta }{a_{n-k}^{(\alpha -1)}}\Vert D_t^{(\alpha -1)}{\eta }^k{\Vert }^2+\sum _{k=1}^n\frac{\tau \mu a_{n-k}^{(\alpha -1)}}{8\delta }{\Vert {\xi }^k\Vert }^2.\end{array}$$

$$\begin{array}{cc}\hfill M_3=-2\tau \sum _{k=1}^n({\partial }_t{\theta }^k,{\rho }^k)& \le \sum _{k=1}^n(9\Vert {\int }_{t_{k-1}}^{t_k}{\theta }_t{ds\Vert }^2+\frac{1}{4}\Vert {\rho }^k{\Vert }^2)\hfill \\ \hfill & \le \sum _{k=1}^n(9\tau {\int }_{t_{k-1}}^{t_k}\Vert {\theta }_t{\Vert }^{2}ds+\frac{1}{4}\Vert {\rho }^k{\Vert }^2).\hfill \end{array}$$

$$\begin{array}{cc}\hfill M_4=2\tau \sum _{k=1}^n({\theta }^k,{\rho }^k)& \le \sum _{k=1}^n(36{\tau }^2\Vert {\theta }^k{\Vert }^2+\frac{1}{4}\Vert {\rho }^k{\Vert }^2).\hfill \end{array}$$

$$\begin{array}{cc}\hfill M_5=2\tau \sum _{k=1}^n(R_v^k,{\xi }^k)& \le \sum _{k=1}^n(\frac{72\tau \delta }{\mu a_{n-k}^{(\alpha -1)}}\Vert R_v^k{\Vert }^2+\frac{\tau \mu a_{n-k}^{(\alpha -1)}}{8\delta }\Vert {\xi }^k{\Vert }^2)\hfill \\ \hfill & \le \sum _{k=1}^n\frac{72{\tau }^2\delta }{\mu a_{n-k}^{(\alpha -1)}}{\int }_{t_{k-1}}^{t_k}\Vert v_{tt}{\Vert }^{2}ds+\sum _{k=1}^n\frac{\tau \mu a_{n-k}^{(\alpha -1)}}{8\delta }{\Vert {\xi }^k\Vert }^2.\hfill \end{array}$$

$$\begin{array}{cc}\hfill M_6=2\tau \sum _{k=1}^n(R_w^k,{\rho }^k)& \le 36{\tau }^2\sum _{k=1}^n\Vert R_w^k{\Vert }^2+\frac{1}{4}\sum _{k=1}^n{\Vert {\rho }^k\Vert }^2\hfill \\ \hfill & \le 36{\tau }^3\sum _{k=1}^n{\int }_{t_{k-1}}^{t_k}\Vert w_{tt}{\Vert }^{2}ds+\sum _{k=1}^n\frac{1}{4}{\Vert {\rho }^k\Vert }^2.\hfill \end{array}$$

$$\begin{array}{cc}\hfill M_7=2\tau \mu \sum _{k=1}^n(R^k,{\xi }^k)& \le \sum _{k=1}^n\frac{72\delta \tau \mu }{a_{n-k}^{(\alpha -1)}}\Vert R^k{\Vert }^2+\frac{\tau \mu }{8\delta }\sum _{k=1}^n{a}_{n-k}^{(\alpha -1)}{\Vert {\xi }^k\Vert }^2\hfill \\ \hfill & \le \sum _{k=1}^n\frac{72\tau \mu \delta }{a_{n-k}^{(\alpha -1)}}{\tau }^{6-2\alpha }(\underset{0\le s\le T}{max}\Vert v_{tt}{\left(s\right)\Vert )}^2+\sum _{k=1}^n\frac{\tau \mu a_{n-k}^{(\alpha -1)}}{8\delta }{\Vert {\xi }^k\Vert }^2.\hfill \end{array}$$

Substitute all the above estimates into (25), we obtain

$$\begin{array}{cc}\hfill G^n& \le G^0+\frac{\tau \mu }{2\delta }\sum _{k=1}^n{a}_{n-1}^{(\alpha -1)}\Vert {\xi }^0{\Vert }^2+\frac{18t_n^{\alpha -1}\Gamma (2-\alpha )}{\mu }{\int }_0^T\Vert {\eta }_t{\Vert }^{2}ds+18\mu \tau t_n^{\alpha -1}\Gamma (2-\alpha ){\Vert D_{\tau }^{\alpha -1}{\eta }^k\Vert }^2\hfill \\ \hfill & +36\tau {\int }_0^T\Vert {\theta }_t{\Vert }^{2}ds+36{\tau }^2\sum _{k=1}^n\Vert {\theta }^k{\Vert }^2+\frac{72{\tau }^2{t}_n^{\alpha -1}\Gamma (2-\alpha )}{\mu }{\int }_0^T\Vert v_{tt}{\Vert }^{2}ds+36{\tau }^3{\int }_0^T{\Vert w_{tt}\Vert }^{2}ds\hfill \\ \hfill & +72\mu {\tau }^{6-2\alpha }T{t}_n^{\alpha -1}(\underset{0\le s\le T}{max}\Vert v_{tt}{\left(s\right)\Vert )}^2+\sum _{k=1}^n\frac{\tau \mu a_{n-k}^{(\alpha -1)}}{2\delta }\Vert {\xi }^k{\Vert }^2+\sum _{k=1}^n\frac{18\tau \mu }{\delta }\Vert {\xi }^k{\Vert }^2+\frac{3}{4}\sum _{k=1}^n{\Vert {\rho }^k\Vert }^2.\hfill \end{array}$$

Noticing ${\rho }^0={\xi }^0=0$, and the following inequality [37]

$$\begin{array}{cc}\hfill |D_{\tau }^{\alpha -1}{v}^k|& \le |D_{\tau }^{\alpha -1}{v}^k-{}_0^c{D}_t^{\alpha -1}{v}^k|+|{}_0^c{D}_t^{\alpha -1}{v}^k|\hfill \\ \hfill & \le O\left({\tau }^{3-\alpha }\right)+\underset{0\le s\le T}{max}\left|{}_0^c{D}_t^{\alpha -1}v\left(s\right)\right|.\hfill \end{array}$$

Choosing $\tau$ to satisfy $\frac{18\tau \mu }{\delta }<\frac{1}{4}$, by the error estimation of elliptic projection and using the discrete Gronwall’s Lemma, we obtain

$$\begin{array}{cc}\hfill \Vert {\rho }^n{\Vert }^2+{\Vert {\xi }^n\Vert }^2& \le C[h^2{\int }_0^T(\Vert v_t{\Vert }_3+\Vert w_t{\Vert }_2{)}^{2}ds+{\tau }^2\sum _{k=1}^n{\Vert {\theta }^k\Vert }^2\hfill \\ \hfill & +h^2(O\left({\tau }^{3-\alpha }\right)+\underset{0\le s\le T}{max}\Vert {}_0^c{D}_t^{\alpha -1}v\left(s\right){\Vert }_3+\underset{0\le s\le T}{max}\Vert {}_0^c{D}_t^{\alpha -1}w\left(s\right){{\Vert }_2)}^2\hfill \\ \hfill & +{\tau }^2{\int }_0^T\Vert v_{tt}{\Vert }^{2}ds+{\tau }^3{\int }_0^T\Vert w_{tt}{\Vert }^{2}ds+{\tau }^{6-2\alpha }(\underset{0\le s\le T}{max}\Vert v_{tt}{\left(s\right)}^{\Vert )2}].\hfill \end{array}$$

(26)

Next, we estimate $\mid v_h^n-v^n{\mid }_1$, choosing ${\varphi }_h={\xi }^n$ in (19)(a)

$$\begin{array}{c}\hfill ({\partial }_t{\rho }^n,{\Pi }_h^{*}{\xi }^n)+({\partial }_t{\theta }^n,{\Pi }_h^{*}{\xi }^n)-({\theta }^n,{\Pi }_h^{*}{\xi }^n)+B({\xi }^n,{\Pi }_h^{*}{\xi }^n)=(R_w^n,{\Pi }_h^{*}{\xi }^n).\end{array}$$

Using Lemmas 2 and 3 and Young inequality

$$\begin{array}{cc}\hfill \frac{1}{2}{\left|{\xi }^n\right|}_1^2& \le C\parallel {\partial }_t{\rho }^n{\parallel }^2+C\parallel {\partial }_t{\theta }^n{\parallel }^2+C\parallel {\theta }^n{\parallel }^2+C\parallel R_w^n{\parallel }^2+\frac{1}{C}{\parallel {\xi }^n\parallel }^2.\hfill \end{array}$$

(27)

To estimate the right-hand side of (27), we have

$$\begin{array}{cc}\hfill & {\displaystyle \parallel {\partial }_t{\rho }^n{\parallel }^2=\parallel \frac{1}{\tau }{\int }_{t_{n-1}}^{t_n}{\rho }_t{\left(s\right)ds\parallel }^2\le \frac{1}{\tau }{\int }_{t_{n-1}}^{t_n}{\parallel {\rho }_t\left(s\right)\parallel }^{2}ds.}\hfill \\ \hfill & {\displaystyle \parallel {\partial }_t{\theta }^n{\parallel }^2=\parallel \frac{1}{\tau }{\int }_{t_{n-1}}^{t_n}{\theta }_t{\left(s\right)ds\parallel }^2\le \frac{1}{\tau }{\int }_{t_{n-1}}^{t_n}{\parallel {\theta }_t\left(s\right)\parallel }^{2}ds.}\hfill \\ \hfill & {\displaystyle \parallel R_w^n{\parallel }^2=\parallel \frac{w^n-w^{n-1}}{\tau }-w_t^n{\parallel }^2=\parallel \frac{1}{\tau }{\int }_{t_{n-1}}^{t_n}(t_{n-1}-s)w_{tt}\left(s\right){ds\parallel }^2\le \tau {\int }_{t_{n-1}}^{t_n}{\parallel w_{tt}\left(s\right)\parallel }^{2}ds.}\hfill \end{array}$$

The above inequality leads to

$$\begin{array}{cc}\hfill |{\xi }^n{|}_1^2& {\displaystyle \le C(\frac{1}{\tau }{\int }_{t_{n-1}}^{t_n}\parallel {\rho }_t{\left(s\right)\parallel }^{2}ds+\frac{1}{\tau }{\int }_{t_{n-1}}^{t_n}\parallel {\theta }_t{\left(s\right)\parallel }^{2}d\tau +\tau {\int }_{t_{n-1}}^{t_n}\parallel w_{tt}{\left(s\right)\parallel }^{2}ds)+\frac{1}{C}{\parallel {\xi }^n\parallel }^2.}\hfill \end{array}$$

Combining (26) and the error estimation of MFVE elliptic projection, we obtain

$$\begin{array}{c}\hfill |{\xi }^n{|}_1\le C(h+\tau ).\end{array}$$

(28)

Finally, by applying (27), (28) and triangle inequality, the theorem is concluded. □

6. Numerical Examples

In this section, we will present two numerical examples to verify our MFVE method. The numerical results show the efficiency and accuracy order of the proposed scheme. The time-fractional damping beam vibration equation is considered as follows

$$\left\{\begin{array}{cc}\left(a\right){\mu }_0^c{D}_t^{\alpha }u+u_{tt}+a^2{u}_{xxxx}=g(x,t),\hfill & (x,t)\in (0,L)\times (0,T],\hfill \\ \left(b\right)u(x,0)=\phi \left(x\right),u_t(x,0)=\psi \left(x\right),\hfill & x\in (0,L),\hfill \\ \left(c\right)u(0,t)=u(L,t)=0,u_{xx}(0,t)=u_{xx}(L,t)=0,\hfill & t\in [0,T].\hfill \end{array}\right.$$

(29)

Here, we consider steel material with uniform cross section and uniform mass, $a=\sqrt{EI/\left(\rho A\right)}$, where the material density $\rho =7850$ kg/m${}^2$, the elastic modulus $E=1.96\times {10}^{11}$ Pa, the cross-section area $A=1.5\times {10}^{-2}$ m${}^2$ and the cross-section moment of inertia $I=1.25\times {10}^{-5}$ m${}^{-2}$.

Example 1.

We choose $T=1,L=1$, $\mu =1$ in (29), the external force applied on the beam is $g(x,t)=(12t^2+a^2{\pi }^4{t}^4+\frac{\Gamma \left(5\right)}{\Gamma (5-\alpha )}{t}^{4-\alpha })sin\left(\pi x\right)$, then we can obtain the exact solution $u(x,t)=t^{4}sin\left(\pi x\right)$. $v_h,w_h$ are solved by MFVE scheme, $u_h$ can be obtained by $v_h$ using the backward Euler method. We fix the spatial step size $h=1/1000$, select the time-step length $\tau =1/10,1/20,1/40,1/80$ and give the error results of $u,v,w$ in

Table 1. $L^2$-norm errors and temporal convergence order of MFVE method.

$\mathit{\alpha }$	$\mathit{\tau }$	$\parallel \mathit{u}-{\mathit{u}}_{\mathit{h}}{\parallel }_{0,\mathit{h}}$	Order	$\parallel \mathit{v}-{\mathit{v}}_{\mathit{h}}{\parallel }_{0,\mathit{h}}$	Order	$\parallel \mathit{w}-{\mathit{w}}_{\mathit{h}}{\parallel }_{0,\mathit{h}}$	Order
$\alpha =1.3$	1/10	3.0223$\times {10}^{-1}$	−	4.0686$\times {10}^{-1}$	−	6.3215	−
	1/20	1.4641$\times {10}^{-1}$	1.046	2.0832$\times {10}^{-1}$	0.965	3.1727	0.995
	1/40	7.1948$\times {10}^{-2}$	1.025	1.0506$\times {10}^{-1}$	0.988	1.5877	0.999
	1/80	3.5650$\times {10}^{-2}$	1.013	5.2749$\times {10}^{-2}$	0.994	7.9404$\times {10}^{-1}$	1.000
$\alpha =1.7$	1/10	3.0233$\times {10}^{-1}$	−	4.0756$\times {10}^{-1}$	−	6.3245	−
	1/20	1.4641$\times {10}^{-1}$	1.046	2.0820$\times {10}^{-1}$	0.979	3.1726	0.995
	1/40	7.1945$\times {10}^{-2}$	1.025	1.0505$\times {10}^{-1}$	0.991	1.5877	0.999
	1/80	3.5648$\times {10}^{-2}$	1.013	5.2763$\times {10}^{-2}$	0.996	7.9393$\times {10}^{-1}$	1.000

and

Table 2. $H^1$-harf-norm errors and temporal convergence order of MFVE method.

$\mathit{\alpha }$	$\mathit{\tau }$	$|\mathit{u}-{\mathit{u}}_{\mathit{h}}{|}_{1,\mathit{h}}$	Order	$|\mathit{v}-{\mathit{v}}_{\mathit{h}}{|}_{1,\mathit{h}}$	Order	$|\mathit{w}-{\mathit{w}}_{\mathit{h}}{|}_{1,\mathit{h}}$	Order
$\alpha =1.3$	1/10	9.4947$\times {10}^{-1}$	−	1.2782	−	1.9859×10	−
	1/20	4.5997$\times {10}^{-1}$	1.046	6.5446$\times {10}^{-1}$	0.966	9.9673	0.995
	1/40	2.2603$\times {10}^{-1}$	1.025	3.3004$\times {10}^{-1}$	0.988	4.9881	0.999
	1/80	1.1200$\times {10}^{-1}$	1.013	1.6572$\times {10}^{-1}$	0.994	2.4946	1.000
$\alpha =1.7$	1/10	9.4979$\times {10}^{-1}$	−	1.2804	−	1.9869×10	−
	1/20	4.5996$\times {10}^{-1}$	1.046	6.5407$\times {10}^{-1}$	0.969	9.9672	0.995
	1/40	2.2602$\times {10}^{-1}$	1.025	3.3001$\times {10}^{-1}$	0.987	4.9878	0.999
	1/80	1.1200$\times {10}^{-1}$	1.013	1.6576$\times {10}^{-1}$	0.993	2.4942	1.000

. The results show that the order of time convergence is approximately 1, which is consistent with the theoretical results in Theorem 2.

Next, we chose the spatial step $h=\frac{1}{N}$, time step $\tau =\frac{1}{M}$. When $M=N^2$, the errors and spatial convergence orders are shown in

Table 3. $L^2$-norm errors and spatial convergence order of MFVE method.

$\mathit{\alpha }$	h	$\parallel \mathit{u}-{\mathit{u}}_{\mathit{h}}{\parallel }_{0,\mathit{h}}$	Order	$\parallel \mathit{v}-{\mathit{v}}_{\mathit{h}}{\parallel }_{0,\mathit{h}}$	Order	$\parallel \mathit{w}-{\mathit{w}}_{\mathit{h}}{\parallel }_{0,\mathit{h}}$	Order
$\alpha =1.3$	1/8	4.4899$\times {10}^{-1}$	−	6.6632$\times {10}^{-2}$	−	1.2098	−
	1/16	1.1115$\times {10}^{-2}$	2.014	1.6676$\times {10}^{-2}$	1.999	3.0108$\times {10}^{-1}$	2.007
	1/32	2.7719$\times {10}^{-3}$	2.004	4.1700$\times {10}^{-3}$	2.000	7.5184$\times {10}^{-2}$	2.002
	1/64	6.9254$\times {10}^{-4}$	2.001	1.0426$\times {10}^{-3}$	2.000	1.8790$\times {10}^{-2}$	2.000
$\alpha =1.7$	1/8	4.4911$\times {10}^{-2}$	−	6.6675$\times {10}^{-2}$	−	1.2104	−
	1/16	1.1117$\times {10}^{-3}$	2.014	1.6683$\times {10}^{-2}$	1.999	3.0118$\times {10}^{-1}$	2.007
	1/32	2.7722$\times {10}^{-3}$	2.004	4.1717$\times {10}^{-3}$	2.000	7.5200$\times {10}^{-2}$	2.002
	1/64	6.9260$\times {10}^{-4}$	2.001	1.0430$\times {10}^{-3}$	2.000	1.8793$\times {10}^{-2}$	2.001

and

Table 4. $H^1$-harf-norm errors and spatial convergence order of MFVE method.

$\mathit{\alpha }$	h	$|\mathit{u}-{\mathit{u}}_{\mathit{h}}{|}_{1,\mathit{h}}$	Order	$|\mathit{v}-{\mathit{v}}_{\mathit{h}}{|}_{1,\mathit{h}}$	Order	$|\mathit{w}-{\mathit{w}}_{\mathit{h}}{|}_{1,\mathit{h}}$	Order
$\alpha =1.3$	1/8	1.4015$\times {10}^{-1}$	−	2.0800$\times {10}^{-1}$	−	3.7764	−
	1/16	3.4863$\times {10}^{-2}$	2.007	5.2305$\times {10}^{-2}$	1.992	9.4436$\times {10}^{-1}$	2.000
	1/32	8.7046$\times {10}^{-3}$	2.002	1.3095$\times {10}^{-2}$	1.998	2.3610$\times {10}^{-1}$	2.000
	1/64	2.1755$\times {10}^{-3}$	2.001	3.2750$\times {10}^{-3}$	2.000	5.9026$\times {10}^{-2}$	2.000
$\alpha =1.7$	1/8	1.4019$\times {10}^{-1}$	−	2.0812$\times {10}^{-1}$	−	3.7781	−
	1/16	3.4869$\times {10}^{-2}$	2.007	5.2328$\times {10}^{-2}$	1.992	9.4465$\times {10}^{-1}$	2.000
	1/32	8.7057$\times {10}^{-3}$	2.002	1.3101$\times {10}^{-2}$	1.998	2.3615$\times {10}^{-1}$	2.000
	1/64	2.1756$\times {10}^{-3}$	2.001	3.2762$\times {10}^{-3}$	2.000	5.9035$\times {10}^{-2}$	2.000

, respectively. The above table shows that the displacement, bending moment and velocity of the vibration beam in the sense of $L^2$ norm and $H^1$-harf norm are more approximate than the theoretical estimates, which proves the effectiveness of the MFVE scheme.

Figure 1 and Figure 2 show the function images of the numerical solution and the exact solution at the last time point. It is proved that the numerical solution fits the exact solution effectively.

Figure 3 and Figure 4 show the contour plots of the numerical solution and the exact solution of displacement, respectively, and Figure 5 and Figure 6 show the contour plots of the numerical solution and the exact solution of bending moment, respectively. It can be seen from the figure that the numerical solution approximates the exact solution at different grid points, which proves the effectiveness of the MFVE method.

Example 2.

We choose $\phi \left(x\right)=sin\left(\pi x\right)$, $\psi \left(x\right)=0$, $T=5$, $L=1$, $\alpha =1.3$ and $x=0.5$ in (29). Suppose that the beam is in free vibration, that is, $g(x,t)=0$, different μ values were taken to verify the influence of material damping on beam vibration. The obtained results are shown in Figure 7, from which it can be concluded that the greater the damping coefficient of the material, the faster the vibration attenuation rate of the beam.

Next, we fixed $\mu =3$, changed the value of $\alpha$ and observed the vibration curve at the midpoint of the damping beam. It can be seen from Figure 8 that the attenuation rate of the beam vibration decreases with the increase in the order of the fractional derivative. In addition, it is generally shown that the peaks of these curves gradually increase and shift to the right as the order of the fractional derivatives increases.

7. Conclusions and Suggestions

In this paper, a mixed finite volume element method is proposed to solve the fractal-order damped beam vibration equation. By introducing two auxiliary variables with practical significance, the original fourth-order problem is transformed into a second-order equation system. The stability and convergence of the scheme are analyzed. The numerical examples demonstrate the effectiveness of the proposed method, and it can be seen that the larger the damping coefficient and the smaller the order of the fractional derivative, the faster the attenuation frequency of the beam vibration.

Although there are some other methods which can be used to solve such problems, the MFVE method shows its advantages: (i) The feature of the finite volume element scheme is retained, so the local conservation of physical quantities can be preserved. (ii) Two physical quantities with practical significance can be solved at the same time; thus, the computing cost is reduced. (iii) Compared with the finite element method, the space smoothness requirement is lower.

In the future work, we can apply this method to other types of beam vibration equations, such as beam vibration equations with structural damping, non-uniform beam vibration equations, etc. At the same time, other methods can be used to discretize the time derivatives to improve the accuracy of the method.