分数阶微积分理论下粘弹性矩形薄膜自由振动

Free Vibration of the Viscoelastic Rectangular Films under the Fractional Calculus Theory

摘要: 薄膜振动问题一直都备受学者关注，经典的薄膜理论是建立在整数阶微积分上来刻画薄膜振动的波动方程，其在建模过程中忽略了材料记忆性的特征。本文针对粘弹性薄膜材料特征具有时间记忆性，提出了一种新的薄膜振动模型，解决了建立在整数阶微积分理论下波动方程无法准确刻画材料的时间记忆性难题。分数阶微积分模型所反映出来的性质与其整个发展史密切相关，本文基于分数阶微积分理论。将Westerlund提出的分数阶模型运用到粘弹性薄膜自由振动中，得到薄膜自由振动的分数阶波动方程。结合矩形和圆形薄膜初边界条件，建立了粘弹性薄膜自由振动的场系方程。运用分离变量法，拉普拉斯变换对矩形薄膜波动方程进行求解。结果表明：在矩形薄膜的自由振动中，体现影响时间因子的分数阶阶数 对振型的影响非常明显。

Abstract: The problem of thin film vibration has always attracted scholars. The classical thin film theory is based on the integer order calculus to characterize the wave equation of thin film vibration, which ignores the characteristics of material memory in the modeling process. In this paper, a new film vibration model is proposed for the temporal memory of viscoelastic film materials, which solves the problem of time memory that the wave equations based on integer calculus theory can not accurately characterize materials. The properties reflected by the fractional-order calculus model are closely related to its entire development history, and this paper is based on the theory of fractional-order calculus. The fractional order model proposed by Westerlund is applied to the free vibration of viscoelastic film to obtain the fractional order wave equation of free vibration of thin film. Combining the initial conditions and boundary conditions of rectangular and circular films, a field system equation for the free vibration of viscoelastic films is established. Using the separation variable method, the Laplace transform solves the rectangular thin film wave equation. The results show that in the free vibration of the rectangular film, the fractional order   of the influencing time factor has a very obvious effect on the mode shape.

1. 引言

膜结构是20世纪后期逐渐出现的一种新结构，其材料是重要组成部分。就目前来说对薄膜自由振动的研究不多，经典的薄膜理论是建立在整数阶微积分上来刻画薄膜振动的波动方程。整数阶导数模型所刻画出来的是某时刻的变化或某种性质。本文研究的矩形薄膜是粘弹性材料，粘弹性材料的力学性质依赖时间记忆性，整数阶微积分模型不能准确刻画材料的性质，由此便出现了一些分数阶模型来补充整数阶模型的不足之处，分数阶Kelvin-Voigt、分数阶Maxwell、分数阶Zener和分数阶Thonson模型 [1] 等。分数阶微积分所反应出来的性质与该现象的整个发展历史有关 [2]。分数阶导数的定义有两种，其一是Riemann和Liouville [3] 提出的，它是一个给定函数与幂次函数的卷积的导数。另外一个是Caputo [4] 提出来的，它是一个给定函数与幂次函数的局部导数的卷积。其中后者克服了前者的强奇异性，Caputo提出的分数阶导数更适合解决实际问题。相关学者针对分数阶做了一些研究，如Sharma和Kumar [5] 用分数阶勒让德函数对对流直翅片的分数模型进行了数值研究。Lai等 [6] 基于分数阶模型(FOM)对电动汽车锂离子电池荷电状态(SOC)和功率状态(SOP)进行理论分析。Carcione [7] 将时间上的分数阶导数转化为空间上的分数阶导数，推导出时间上为二阶偏导数和空间上为Riesz分数阶导数的常Q波动方程。推导出了时间上为二阶偏导数和空间上为Riesz分数阶导数的常Q波动方程。常Q波动方程可以避免对时间分数阶偏导数的计算，减少了变量的存储，同时对于空间上的Riesz分数阶导数，可以采用快速傅里叶变换进行计算，相比于时间上的分数阶波动方程，其计算量大为简化。国内一些学者在分数阶微积分方面有了一些成果。刘林超等 [8] 建立了分数阶导数模型的粘弹性桩振动方程，振动方程的数值解是用数值积分的方法求得的。通过分析表明：波速受到分数阶阶数的显著影响。骆文和 [9] 研究了分数阶导数模型的粘弹性土层扭转振动控制方程式，通过求解进行一系列的分析，进而发现：桩顶复刚度受土体的本构模型参数、桩土的形变模量比与桩的长径比受到分数阶阶数的显著影响。Zhu [10] 等人推导了一阶速度压力–应变形式的近似常Q (NCQ)波动方程。它是由时间上的二阶偏导数和空间上解耦的分数阶Laplace算子表示的，其中解耦的分数阶Laplace算子分别对应相速度的频散和振幅的衰减，文章中己经将NCQ波动方程运用到逆时偏移中，并取得了较好的成像效果。至此，已经知道用分数阶微分模型来刻画粘弹性材料的力学性质有着重大的现实意义。利用分数阶导数具有的时间记忆性，来刻画粘弹性材料力学特性是一种很好的模型。具体可以参考Mainardi的书籍 [11]，书中详细研究了运用分数阶微积分理论去讨论波在粘弹性介质中的传播，进一步揭示了时间分数阶波动方程在动力学所表现的行为。

本文基于分数阶微积分理论，来研究粘弹性矩形薄膜自由振动，并用图呈现的方式来观察影响时间因子的分数阶阶数 $\alpha$ 对薄膜振动振幅的影响。

2. 基础知识

本文所研究的粘弹性材料矩形薄膜振动模型是建立在Caputo意义的分数阶微积分理论基础之上。确定模型在求解分数阶方程过程中，会遇到一些特殊函数，因此首先给出一些相关函数的基本定义。

Mittag-Leffler函数相关定义及其性质：

定义1：单参数的Mittag-Leffler函数 [12]：

$E_{\alpha }\left(z\right)={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}\frac{z^k}{\Gamma \left(\alpha k+1\right)}}$ (1)

显然， $E_1\left(z\right)={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}\frac{z^k}{\Gamma \left(k+1\right)}={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}\frac{z^k}{k!}={\text{e}}^z}}$。

定义2：双参数的Mittag-Leffler函数 [12]：

$E_{\alpha ,\beta }\left(z\right)={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}\frac{z^k}{\Gamma \left(\alpha k+\beta \right)}}$ (2)

以上定义中的 $\alpha >0$， $\beta \in ℂ$， $z\in ℂ$。其中， $\Gamma$ 是Gamma函数。特别的， $E_{\alpha }\left(0\right)=1$， $E_{\alpha ,\beta }\left(0\right)=1$。

性质1：Mittag-Leffler函数的Laplace变换 [13] 公式形式如下：

(3)

定义3：Caputo微分定义 [4]：

${}_a{}^{C}D{}^{\alpha }f\left(t\right)=\frac{1}{\Gamma \left(n-\alpha \right)}{\displaystyle {\int }_0^t{\left(t-\tau \right)}^{n-\alpha -1}{f}^{\left(n\right)}\left(\tau \right)\text{d}\tau },\left(n-1<\alpha \le n,t>a\right)$ (4)

性质2：Caputo分数阶微分算子的拉普拉斯变换 [14] 为

$L\left\{{}_a{}^{C}D{}_t^{\alpha }f\left(t\right);s\right\}=s^{\alpha }F\left(s\right)-{\displaystyle \underset{k=0}{\overset{n-1}{\sum }}{s}^{\alpha -k-1}{f}^{\left(k\right)}}\left(0\right),\left(n-1<\alpha \le n\right)$ (5)

3. 问题模型

为研究粘弹性材料下矩形薄膜的波动方程，建立如图1所示的模型，考虑边长为a和b且边固定的矩形薄膜的自由振动。

传统的整数阶二维波动方程：

$\frac{{\partial }^{2}u}{\partial t^2}=c^2\left(\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}\right)$ (6)

Figure 1. Free vibration of the rectangular membrane

图1. 矩形薄膜的自由振动

本文中，考虑的是粘弹性材料，所以采用分数阶模型来刻画矩形薄膜的波动方程。由固体力学可知，理想弹性体的应力-应变关系满足胡克定律： $\sigma \left(t\right)~\epsilon \left(t\right)$，牛顿流体满足牛顿定律： $\sigma \left(t\right)~\frac{\text{d}\epsilon \left(t\right)}{\text{d}t}$。如果把胡克定律改写成这样的形式： $\sigma \left(t\right)~\frac{{\text{d}}^0\epsilon \left(t\right)}{{\text{d}}^{0}t}$，由此很容易得到处于弹性和牛顿流体之间的材料的本构关系应该是： $\sigma \left(t\right)~\frac{{\text{d}}^{\beta }\epsilon \left(t\right)}{{\text{d}}^{\beta }t}\left(0\le \beta \le 1\right)$，这是Scott Blair [15] 和Gerasimor [16] 提出的革命性观点。这不仅仅是刻画了弹性体和牛顿流体，并且包含了粘弹性材料。经典力学的基本力学本构关系：刚体的牛顿第二定律： $F=m\frac{{\text{d}}^{2}x}{\text{d}{t}^2}$，该本构关系涉及到了加速度。在描述多项介质的力学行为时，这个经典模型存在一定的缺陷。从微分方程建模的观点看，Westerlund [17] 建议采用以下简单统一的分数阶模型来描述： $F=\rho \frac{{\text{d}}^{\alpha }x}{\text{d}{t}^{\alpha }}$， $\left(0\le \alpha \le 2\right)$。根据Westerlund分数阶模型，粘弹性材料薄膜自由振动的加速度可以这样表示： $\frac{{\text{d}}^{\alpha }u}{\text{d}{t}^{\alpha }}\left(1\le \alpha \le 2\right)$。所以可将整数阶波动方程(6)转化为粘弹性材料下的分数阶波动方程，形式如下：

$\frac{{\partial }^{\alpha }u}{\partial t^{\alpha }}=c^{２}\left(\frac{{\partial }^{２}u}{\partial x^{２}}+\frac{{\partial }^{２}u}{\partial y^{２}}\right),1\le \alpha \le 2$ (7)

其中c表示波的传播速度。

结合初边界条件来考虑，容易得到定解问题：

$\{\begin{array}{l}\frac{{\partial }^{\alpha }u}{\partial t^{\alpha }}=c^{２}\left(\frac{{\partial }^{２}u}{\partial x^{２}}＋\frac{{\partial }^{２}u}{\partial y^{２}}\right),\left(0\le x\le a,0\le y\le b,0\le t\right)\\ {u|}_{x=0}={u|}_{x=a}={u|}_{y=0}={u|}_{y=b}=0\\ {u|}_{t=0}=f\left(x,y\right), {\frac{\partial u}{\partial t}|}_{t=0}=g\left(x,y\right)\end{array}$ (8)

4. 定解问题的求解

求解方程(8)，采用分离变量法，令：

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

将式(9)代入分数阶波动方程(8)中的控制方程，并令其比值为 $-{\omega }^2$，得到，

$T^{\left(\alpha \right)}\left(t\right)V\left(x,y\right)=c^2\left(V_{xx}\left(x,y\right)+V_{yy}\left(x,y\right)\right)T\left(t\right)=-{\omega }^2$ (10)

其中 $\omega$ 为待求常数。

再将方程(10)进一步处理得到：

$\frac{T^{\left(\alpha \right)}\left(t\right)}{T\left(t\right)}=c^2\left(\frac{V_{xx}\left(x,y\right)+V_{yy}\left(x,y\right)}{V\left(x,y\right)}\right)=-{\omega }^2$ (11)

得到如下两个方程：

$T^{\left(\alpha \right)}\left(t\right)+{\omega }^{2}T\left(t\right)=0$ (12)

$V_{xx}\left(x,y\right)+V_{yy}\left(x,y\right)+c^2{\omega }^{2}V\left(x,y\right)=0$ (13)

由边界条件 ${u|}_{x=0}={u|}_{x=a}={u|}_{y=0}={u|}_{y=b}=0$，则有如下边界条件：

$\{\begin{array}{l}V\left(0,y\right)=V\left(a,y\right)=0\\ V\left(x,0\right)=V\left(x,a\right)=0\end{array}$ (14)

再对方程(13)进行分离变量，令：

$V\left(x,y\right)=X\left(x\right)Y\left(y\right)$ (15)

代入方程(13) 分离变量，令其比值为 $-{\mu }^2$ 得到：

$\frac{X^{″}\left(x\right)}{X\left(x\right)}+\frac{Y^{″}\left(y\right)}{Y\left(y\right)}+c^2{\omega }^2=-{\mu }^2$ (16)

根据边界条件(14)，则式(16)满足边界条件：

$\{\begin{array}{l}X\left(0\right)=X\left(a\right)=0\\ Y\left(0\right)=Y\left(b\right)=0\end{array}$ (17)

此时，得到两个固有值问题：

(I) $\{\begin{array}{l}{X}^{″}\left(x\right)+{\mu }^{2}X\left(x\right)=0\\ X\left(0\right)=0,X\left(a\right)=0\end{array}$ 和(II) $\{\begin{array}{l}{Y}^{″}\left(y\right)+\left(c^2{\omega }^2-{\mu }^2\right)Y\left(y\right)=0\\ Y\left(0\right)=0,Y\left(b\right)=0\end{array}$

方程(I)的通解为：

$X\left(x\right)=A\cos \mu x+B\sin \mu x$ (18)

并且满足边界条件 $X\left(0\right)=X\left(a\right)=0$。

因为 $X\left(0\right)=0$，所以 $A\cos \mu 0=0$，而 $\cos \mu 0=1$，所以 $A=0$，此时可以将微分方程解的形式简化为：

$X\left(x\right)=B\sin \mu x$ (19)

又因为 $X\left(a\right)=0$，所以 $B\sin \mu a=0$，考虑到 $B\ne 0$，所以只能是 $\sin \mu a=0$。所以 $\mu =\frac{n\pi }{a}$，( ${\mu }^2=\frac{n^2{\pi }^2}{a^2}$, $n=1,2,3,\cdots$ )。

所以方程(I)的解为：

$X_n\left(x\right)=\sin \frac{n\pi }{a}x$ (20)

方程(II)的通解为：

$Y\left(y\right)=A\cos \sqrt{c^2{\omega }^2-{\mu }^2}y+B\sin \sqrt{c^2{\omega }^2-{\mu }^2}y$ (21)

满足边界条件 $Y\left(0\right)=Y\left(b\right)=0$。

同方程(I)，因为 $Y\left(0\right)=0$，所以可以得到 $A=0$，即微分方程(II)的解形式简化为：

$Y\left(y\right)=B\sin \sqrt{c^2{\omega }^2-{\mu }^2}y$ (22)

又因为 $Y\left(b\right)=0$，所以 $B\sin \sqrt{c^2{\omega }^2-{\mu }^2}b=0$，考虑到 $B\ne 0$，所以只能是 $\sin \sqrt{c^2{\omega }^2-{\mu }^2}b=0$。进而有 $\sqrt{c^2{\omega }^2-\mu }=\frac{m\pi }{b}$，( $c^2{\omega }^2-{\mu }^2=\frac{m^2{\pi }^2}{b^2}$, ${\omega }^2=\frac{1}{c^2}\left(\frac{m^2{\pi }^2}{b^2}+\frac{n^2{\pi }^2}{a^2}\right)$, $m=1,2,3,\cdots$ )。

所以方程(II)的解为：

$Y_n\left(x\right)=\sin \frac{m\pi }{b}x$ (23)

所以微分方程(I)和(II)构成的固有值问题对应的固有值为：

${\omega }^2=\frac{{\pi }^2}{c^2}\left(\frac{n^2}{a^2}+\frac{m^2}{b^2}\right)$ (24)

固有函数为：

$V_{mn}\left(x,y\right)=\sin \frac{n\pi }{a}x\sin \frac{m\pi }{b}y$ (25)

下面求解关于时间的分数阶微分方程，对方程(15)进行拉普拉斯变换作用：

$L\left[T^{\left(\alpha \right)}\left(t\right)+{\omega }^{2}T\left(t\right)\right]=L\left[T^{\left(\alpha \right)}\left(t\right)\right]+{\omega }^{2}L\left[T\left(t\right)\right]=0$ (26)

再根据Caputo分数阶导数的拉普拉斯变换式(8)有：

$s^{\alpha }L\left[T\left(t\right)\right]-s^{\alpha -1}T\left(0\right)-s^{\alpha -2}{T}^{\prime }\left(0\right)+{\omega }^{2}L\left[T\left(t\right)\right]=0$ (27)

进而有：

$L\left[T\left(t\right)\right]=\frac{s^{\alpha -1}T\left(0\right)+s^{\alpha -2}{T}^{\prime }\left(0\right)}{s^{\alpha }+{\omega }^2}$ (28)

根据初始条件， ${u|}_{t=0}=f\left(x,y\right),{\frac{\partial u}{\partial t}|}_{t=0}=g\left(x,y\right)$ 容易得到 $T\left(0\right)$ 和 $T^{\prime }\left(0\right)$ 为常数。

对方程(28)进行拉普拉斯逆作用，根据Mittag-Leffler函数的拉普拉斯变换，有：

$T\left(t\right)=L^{-1}\left[\frac{s^{\alpha -1}T\left(0\right)+s^{\alpha -2}{T}^{\prime }\left(0\right)}{s^{\alpha }+{\omega }^2}\right]=E_{\alpha ,1}\left(-{\omega }^2{t}^{\alpha }\right)T\left(0\right)+t{E}_{\alpha ,2}\left(-{\omega }^2{t}^{\alpha }\right)T^{\prime }\left(0\right)$ (29)

所以该微分方程的通解为：

$T_{mn}\left(t\right)=E_{\alpha ,1}\left(-{\omega }_{mn}^2{t}^{\alpha }\right)T_{mn}\left(0\right)+t{E}_{\alpha ,2}\left(-{\omega }_{mn}^2{t}^{\alpha }\right){T^{\prime }}_{mn}\left(0\right)$ (30)

其中， $T_{mn}\left(0\right)$ 和 ${T^{\prime }}_{mn}\left(0\right)$ 是任意常数，具体表达式：

$T_{mn}\left(0\right)=\frac{1}{{\Vert X_n\left(x\right)\Vert }^2{\Vert Y_n\left(y\right)\Vert }^2}{\displaystyle {\int }_0^a{\displaystyle {\int }_0^{b}f\left(x,y\right)X_n\left(x\right)Y_n\left(y\right)\text{d}x\text{d}y}}$ (31)

${T^{\prime }}_{mn}\left(0\right)=\frac{1}{{\Vert X_n\left(x\right)\Vert }^2{\Vert Y_n\left(y\right)\Vert }^2}{\displaystyle {\int }_0^a{\displaystyle {\int }_0^{b}g\left(x,y\right)X_n\left(x\right)Y_n\left(y\right)\text{d}x\text{d}y}}$ (32)

其中，

${\Vert X_n\left(x\right)\Vert }^2={\displaystyle {\int }_0^a{\left(\sin \frac{n\pi }{a}x\right)}^2\text{d}x}=\frac{a}{2}$ (33)

${\Vert Y_n\left(y\right)\Vert }^2={\displaystyle {\int }_0^b{\left(\sin \frac{m\pi }{b}x\right)}^2\text{d}x}=\frac{b}{2}$ (34)

所以，

$T_{mn}\left(0\right)=\frac{4}{ab}{\displaystyle {\int }_0^a{\displaystyle {\int }_0^{b}f\left(x,y\right)X_n\left(x\right)Y_n\left(y\right)\text{d}x\text{d}y}}$ (35)

${T^{\prime }}_{mn}\left(0\right)=\frac{4}{ab}{\displaystyle {\int }_0^a{\displaystyle {\int }_0^{b}g\left(x,y\right)X_n\left(x\right)Y_n\left(y\right)\text{d}x\text{d}y}}$ (36)

故原定解问题的解有如下形式：

$u\left(x,y,t\right)={\displaystyle \underset{m=1}{\overset{\infty }{\sum }}{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\sin \frac{n\pi x}{a}}}\sin \frac{m\pi y}{b}\left(T_{mn}\left(0\right)E_{\alpha ,1}\left(-{\omega }_{mn}^2{t}^{\alpha }\right)+{T^{\prime }}_{mn}\left(0\right)t{E}_{\alpha ,2}\left(-{\omega }_{mn}^2{t}^{\alpha }\right)\right)$ (37)

其中， $T_{mn}\left(0\right)$ 和 ${T^{\prime }}_{mn}\left(0\right)$ 的表达式是(38)和(39)。

验证理论解的正确性，先进行数值上的分析：假设初始条件 $f\left(x,y\right)=xy\left(x-a\right)\left(y-b\right)$， $g\left(x,y\right)=0$，波速 $c=1$。则分数阶波动方程位移解的表达式是：

$u\left(x,y,t\right)={\displaystyle \underset{m,n=1}{\overset{\infty }{\sum }}\sin \frac{n\pi x}{a}}\sin \frac{m\pi y}{b}\left(T_{mn}\left(0\right)E_{\alpha ,1}\left(-{\omega }_{mn}^2{t}^{\alpha }\right)\right)$ (38)

系数：

$\begin{array}{c}{T}_{mn}\left(0\right)=\frac{4}{ab}{\displaystyle {\int }_0^a{\displaystyle {\int }_0^{b}f\left(x,y\right)X_n\left(x\right)Y_n\left(y\right)\text{d}x\text{d}y}}\\ =\frac{4}{ab}{\displaystyle {\int }_0^a{\displaystyle {\int }_0^{b}xy\left(x-a\right)\left(y-b\right)}}\sin \frac{n\pi }{a}x\sin \frac{m\pi }{b}y\text{d}x\text{d}y\\ =\frac{16a^2{b}^2}{n^3{m}^3{\pi }^6}\left(1-\cos n\pi \right)\left(1-\cos m\pi \right)\end{array}$ (39)

所以位移解具体形式是：

$u\left(x,y,t\right)={\displaystyle \underset{m,n=1}{\overset{\infty }{\sum }}\frac{16a^2{b}^2}{n^3{m}^3{\pi }^6}\left(1-\cos n\pi \right)\left(1-\cos m\pi \right)\sin \frac{n\pi x}{a}}\sin \frac{m\pi y}{b}{E}_{\alpha ,1}\left(-{\omega }_{mn}^2{t}^{\alpha }\right)$ (40)

当 $\alpha =2$ 时，该定解问题成了整数阶方程，对应的解为：

$u\left(x,y,t\right)={\displaystyle \underset{m,n=1}{\overset{\infty }{\sum }}\frac{16a^2{b}^2}{n^3{m}^3{\pi }^6}\left(1-\cos n\pi \right)\left(1-\cos m\pi \right)}\sin \frac{n\pi x}{a}\sin \frac{m\pi y}{b}\cos \frac{\sqrt{a^2{m}^2+b^2{n}^2}}{ab}\pi t$ (41)

所以只要说明当 $\alpha =2$ 时， $E_{\alpha ,1}\left(-{\omega }_{mn}^2{t}^{\alpha }\right)$ 和 $\cos \frac{\sqrt{a^2{m}^2+b^2{n}^2}}{ab}\pi t$ 相等即可证明所求解的分数阶位移解是正确的。化简 $E_{\alpha ,1}\left(-{\omega }_{mn}^2{t}^{\alpha }\right)$，其中 ${\omega }_{mn}^2={\pi }^2\left(\frac{n^2}{a^2}+\frac{m^2}{b^2}\right)$ 则有：

$\begin{array}{c}{E}_{2,1}\left(-{\omega }_{mn}^2{t}^2\right)={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}\frac{{\left(-{\omega }_{mn}^2{t}^2\right)}^k}{\Gamma \left(2k+1\right)}}={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}\frac{{\left(-{\pi }^2\left(\frac{n^2}{a^2}+\frac{m^2}{b^2}\right)\right)}^k}{\left(2k\right)!}{t}^2{}^k}\\ ={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}{\left(-1\right)}^k\frac{{\left(\frac{n^2}{a^2}+\frac{m^2}{b^2}\right)}^k}{\left(2k\right)!}{\pi }^2{}^k{t}^2{}^k}\\ ={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}\frac{{\left(-1\right)}^k}{\left(2k\right)!}}{\left(\frac{n^2}{a^2}+\frac{m^2}{b^2}\right)}^k{\pi }^2{}^k{t}^2{}^k\\ ={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}\frac{{\left(-1\right)}^k}{\left(2k\right)!}}{\left(\frac{n^2{b}^2+m^2{a}^2}{a^2{b}^2}\right)}^k{\pi }^2{}^k{t}^2{}^k\\ ={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}\frac{{\left(-1\right)}^k}{\left(2k\right)!}}{\left(\frac{\sqrt{n^2{b}^2+m^2{a}^2}}{ab}\right)}^{2k}{\pi }^2{}^k{t}^2{}^k\\ =\cos \frac{\sqrt{a^2{m}^2+b^2{n}^2}}{ab}\pi t\end{array}$ (42)

所以当分数阶的阶数 $\alpha =2$ 时，分数阶波动方程和整数阶波动方程的位移解相等，这就证明所求的位移解是正确的。

再用matlab软件模拟出本征函数运动的振型图。首先模拟不加时间的振型图，即本征函数 $V_{mn}\left(x,y\right)=\sin \frac{n\pi }{a}x\sin \frac{m\pi }{b}y$ 对应的图。给定长 $a=2$， $b=1$， $c=1$ 时，固有值： ${\omega }^2=\frac{{\pi }^2}{c^2}\left(\frac{n^2}{a^2}+\frac{m^2}{b^2}\right)$。

固有频率(见表1)对应的振型图如图2~5所示：

Table 1. The natural frequency of vibration of a rectangular membrane

表1. 矩形薄膜振动的固有频率

结合图与表1来观察，发现变量m与n的改变，都对本征模的振型有着显著的影响，并且本征模的振型随着固有频率的增大而变得复杂。图2表明在 $y=1$， $x=0.5$ 时振幅达到峰值。图3表明在 $y=\pm 1.5$， $x=0.5$ 时振幅达到峰值。图4表明在 $y=0.5\pm 0.25$， $x=1$ 时振幅达到峰值。图5表明在 $y=0.5$， $x=1$ 时振幅达到峰值，且在 $y=0.17$， $x=1$ 以及 $y=0.83$， $x=1$ 附近达到峰值。在振幅达到峰值以后，会向周围逐渐降低。

当考虑本征模与时间因子结合：

本征模 $V_{mn}\left(x,y\right)=\sin \frac{n\pi }{a}x\sin \frac{m\pi }{b}y$ 结合时间因子 $E_{\alpha ,1}\left(-{\omega }_{mn}^2{t}^{\alpha }\right)$ 以后，就表示各个平面驻波分量的振动。如下一系列图所展示，高表示膜振动的幅度。膜振动是在一个时间段内是一个连续变化的过程，但是没有办法将这一运动变化过程详尽展示在论文里，所以我们选择特定时刻 $t=1$ 来观察，当矩薄膜振动到指定时刻就会停止，截取此时的图形，就有如下所示的一系列图形。改变分数阶阶数 $\alpha$ 的值会观察到不同的振型图(图6~8)。

Figure 2. The mode shape when m is equal to 1 and n is equal to 1

图2. 当m取1，n取1时的振型

Figure 3. The mode shape when m is equal to 2 and n is equal to 1

图3. 当m取2，n取1时的振型

Figure 4. The mode shape when m is equal to 1 and n is equal to 2

图4. 当m取1，n取2时的振型图

Figure 5. The mode shape when m is equal to 1 and n is equal to 3

图5. 当m取1，n取3时的振型图

Figure 6. The mode shape when m is equal to 3, n is equal to 2, t is equal to 1, $\alpha$ is equal to 1.9

图6. 当m取3，n取2， $t=1$， $\alpha =1.9$ 时的驻波分量的振动

Figure 7. The mode shape when m is equal to 3, n is equal to 2, t is equal to 1, $\alpha$ is equal to 1.95

图7. 当m取3，n取2， $t=1$， $\alpha =1.95$ 时的驻波分量的振动

Figure 8. The mode shape when m is equal to 3, n is equal to 2, t is equal to 1, $\alpha$ is equal to 2

图8. 当m取3，n取2， $t=1$， $\alpha =2$ 时的驻波分量的振动

通过以上振型图发现当本征模结合时间因子以后，分数阶的阶数 $\alpha$ 的值，对驻波分量的振动图形有显著的影响，变化更为复杂。整体变化趋势是当时间因子中的变量分数阶的阶数 $\alpha$ 逐渐递增到2时，驻波分量的振幅也逐渐变大，这刚好与 $E_{\alpha ,1}\left(-{\omega }_{mn}^2{t}^{\alpha }\right)$ 关于 $\alpha$ 是递增函数相互吻合。且当 $\alpha =2$ 时，驻波分量的运动画面和文献 [18] 整数阶所呈现的图一致，因此验证了理论解的正确性。

5. 结论

基于分数阶微积分理论，本文主要讨论了粘弹性矩形薄膜的自由振动问题。运用Westerlund分数阶模型将经典的二维波动方程结合初边界条件，建立了粘弹性矩形薄膜自由振动的场系方程。用分离变量法，拉普拉斯变换对方程进行求解。先后用数值分析和软件绘图的方法验证理论解的正确性。同时用软件绘图来研究矩形薄膜振动振型图和体现时间因子影响的分数阶阶数对驻波分量的振幅影响。结果表明固有值对本征模的振型影响很大，振型图随着固有值的增大而变得复杂。同时当振幅达到峰值以后会随着自变量的变化而逐渐降低。当本征模与时间因子结合以后，分数阶阶数 $\alpha$ 对矩形薄膜振动的振幅影响十分显著，不能忽略。当时间因子中的变量分数阶阶数 $\alpha$ 增大时，驻波分量的振幅也逐渐升高。本论文对工程领域的研究提供了一定的理论指导，与此同时，也丰富了粘弹性材料力学性质的研究内容。

基金项目

国家自然科学基金(11862014)。