粘弹性材料下纳米弦的振动

Vibration of Nano Strings in Viscoelastic Materials

摘要: 本文给出了一个具有Caputo时间分数阶导数的纳米弦振动的波动方程在不同边界条件下的解析解。根据Mittag-Leffler函数和Sturm-Liouville问题的一组完整的本征函数，采用分离变量法、拉普拉斯变换法，得到在初始条件、边界条件以及外力条件特殊情况下的一些结果，表明整数阶方程的相应解是时间分数阶方程的特殊情况，且该方程可用于描述复杂介质或粘弹性介质中的记忆特性。

Abstract: In this paper, an analytical solution of the wave equation of nano string vibration with Caputo time-fractional derivative under different boundary conditions is given. According to the Mittag-Leffler function and a complete set of eigenfunctions of Sturm-Liouville problem, some results under the special conditions of initial conditions, boundary conditions and external force conditions are obtained by using the separated variable method and Laplace transform method. It shows that the corresponding solution of the integer-order equation is a special case of the time-fractional order equation, and the equation can be used to describe the memory characteristics in complex media or viscoelastic media.

1. 引言

近些年来随着材料科学的不断发展，涌现出许多具有粘弹性性质的材料。针对粘弹性材料而言，粘弹性材料的高阻尼性与高耗能性，利用整数阶模型无法很好地描述粘弹性材料。对大部分情况而言，整数阶模型以及分数阶模型均可用于描述材料的力学特性，相较于经典的整数阶模型，分数阶模型有着很强的应用性与优越性。这是由于分数阶模型可以利用少量的参数变量阐述粘弹性材料的一些力学特征，从而可以更好地体现出真实化的粘弹性介质的力学特点。整体而言，整数阶微积分理论是研究的基础理论，分数阶微积分理论是基于整数阶微积分理论推导出来的。在不同研究背景下，分数阶模型可以退化为整数阶模型，也可以推广得到更加广义的模型，因而分数阶模型是一个十分重要的模型。

振动现象无处不在，会对现实社会的生产生活带来一定的影响。振动现象的影响存在两面性，我们可以对振动现象的有利影响进行利用。研究一些机械系统的振动可以使我们了解它们的特性，规避它们的缺点，提升它们的工作效率，因此研究一些材料的振动是非常重要的。

随着越来越多的粘弹性材料被应用到不同的研究领域中，人们对粘弹性材料的力学特性也越来越关注。粘弹性材料既含有弹性性质又含有粘性性质，由于粘弹性材料的这种特性，无法利用单一的弹性或粘性关系来很好地描述其力学性质。粘弹性介质在波的振动传播过程中，由于介质具有粘性一定程度上会消耗动能，达到缓冲减震的效果，利用传统的整数阶微分波动方程 [1] 来进行波场模拟，不能很好地解释这些现象。因此，越来越多的研究学者采用分数阶的微分波动方程进行模拟。随着分数阶导数理论不断完善，更简洁更高效的求解方法也不断增加，分数阶导数理论被应用在很多工程领域中，因此越来越重视利用分数阶导数形式描述粘弹性材料的波动方程。这不仅是对经典力学的理论补充，而且推动了粘弹性表面力学的发展。

分数阶微积分理论及其应用已经成为科研工作者关注的热点。Gemant [2] 提出利用分数阶理论进一步阐述分数阶微分算子在部分粘弹性流体材料中的必要性与重要性。Stiassnie [3] 利用分数阶微积分的理论知识推导粘弹性材料的本构模型，研究结果表明相应模型的性能与实验结果拟合程度较好。Bagley和Torvik [4] 基于一种五参数双分数阶模型去描述粘弹性材料的物理性特征，并在此基础上进一步推导得出粘弹性阻尼结构运动方程。Carcione等人 [5] [6] 基于广义Zener模型的线性粘弹性理论，利用频谱松弛的机制来描述本构关系，通过引入记忆变量表示松弛过程，建立了线性粘弹性介质的波动方程，并对粘弹性介质波场特征进行了分析。孙成禹等 [7] [8] 基于常Q理论推导构建在粘弹性介质中的分数阶黏声波与粘弹性波动方程，利用优化交错网络有限差分法，对不同模型下进行正演模拟和偏移模拟。闰启方 [9] 等人对分数阶Kelvin模型描述的粘弹性材料的松弛特性、蠕变特性和动态力学行为进行了研究。李想 [10] 等人详细介绍了几种处理粘弹性波传播问题的分析方法，重点讲解Laplacc变换法以及Laplacc变换在粘弹性波中的应用。袁宇彤 [11] 等人研究了粘弹性地基上薄板的波动和振动问题，讨论了基于分数导数理论的粘弹性地基模型上薄板弯曲波的传播特性以及固有频率对地基的依赖特性。

本文研究了不同边界条件的Caputo时间分数阶波动方程，考虑的时间分数阶波动方程是一般整数阶波动方程的推广，分析了不同边界和初始条件以及外力不同形式的例子。

2. 预备知识

2.1. Riemann-Liouville分数阶导数

Riemann-Liouville分数阶导数定义在科学研究中应用十分广泛，是一种通过先积分再微分的方式构建的分数阶导数。即以Riemann-Liouville分数阶积分定义为基础，由于导数和积分互为逆运算，推导形成Riemann-Liouville分数阶导数。Riemann-Liouville分数阶积分定义分为左侧分数阶积分和右侧分数阶积分，其中具体定义形式分别为：

$I_{a^{+}}^{\alpha }f\left(x\right)=\frac{1}{\Gamma \left(\alpha \right)}{\displaystyle {\int }_a^x{\left(x-\tau \right)}^{\alpha -1}f\left(\tau \right)\text{d}\tau }$， $\left(x>a,\alpha >0\right)$ (1)

$I_{b^{-}}^{\alpha }f\left(x\right)=\frac{1}{\Gamma \left(\alpha \right)}{\displaystyle {\int }_x^b{\left(x-\tau \right)}^{\alpha -1}f\left(\tau \right)\text{d}\tau }$， $(b>x,\alpha >0)$ (2)

为了更加简洁，也可以将Riemann-Liouville分数阶积分定义为

$J_a^{\alpha }f\left(x\right)=\frac{1}{\Gamma \left(\alpha \right)}{\displaystyle {\int }_0^x{\left(x-\tau \right)}^{\alpha -1}f\left(\tau \right)\text{d}\tau }$， $\left(x>a,\alpha >0\right)$ (3)

其中 $a<x<b$， $J_a^{\alpha }$ 是 $\alpha$ 阶的Riemann-Liouvill分数阶积分算子。

利用导数和积分的可逆性，可得Riemann-Liouville分数阶导数定义 [12] 为

$f^{\left(\alpha \right)}\left(x\right)={\left(\frac{\text{d}}{\text{d}x}\right)}^n\frac{1}{\Gamma \left(n-\alpha \right)}{\displaystyle {\int }_a^x{\left(x-\tau \right)}^{n-\alpha -1}f\left(\tau \right)\text{d}\tau }$ (4)

2.2. Caputo分数阶导数

通过Riemann-Liouville分数阶导数定义可以知道，Riemann-Liouville分数阶导数存在有超奇异性。所以，Riemann-Liouville分数阶导数并不利于在工程与物理建模中进行大量使用。而为了解决这个问题，Caputo也提出了弱奇异的分数阶导数的定义，即Caputo分数阶导数。通过Caputo分数阶导数的弱奇异性，可以处理分数阶导数中分数阶的初值问题。

Caputo分数阶导数，定义 [13] 为：

${}_0{}^{C}D{}_t^{\alpha }f\left(t\right)=\{\begin{array}{l}\frac{1}{\Gamma \left(m-\alpha \right)}{\displaystyle {\int }_0^t\frac{f^m\left(\tau \right)}{{\left(t-\tau \right)}^{\alpha +1-m}}\text{d}\tau ,m-1<\alpha \le m,}\\ \frac{{\text{d}}^{m}f\left(t\right)}{\text{d}{t}^m},\alpha =m.\end{array}$ (5)

在有界域 $x\in \left[0,l\right]$ 以及勒贝格可积函数的空间中得到了这种分数阶微分方程的解，

$L\left(0,\infty \right)=\left\{f:{\Vert f\Vert }_1={\displaystyle {\int }_0^{\infty }\left|f\left(t\right)\right|\text{d}t}<\infty \right\}$ (6)

它的拉普拉斯变换 [14] 为：

$\mathcal{L}\left[D_t^{\alpha }f\left(t\right)\right]=s^{\alpha }\mathcal{L}\left[f\left(t\right)\right]-{\displaystyle \underset{k=0}{\overset{m-1}{\sum }}{f}^{\left(k\right)}\left(0_{+}\right)s^{\alpha -1-k}}$， $m-1<\alpha \le m$ (7)

由于多孔介质力学、飞流牛顿流体力学、黏弹性力学研究 [15] 的需要，分数阶导数的研究再度引起了广泛的重视，相关的研究可参考 [16] [17] [18] [19] [20]。由于反常扩散所具有的历史依赖与全域相关的特征，恰好可以由分数阶导数来表示。因此，它在时间分数弛豫和振荡过程，以及时间分数扩散和波动过程 [18] 有着广泛的应用。

2.3. Mittag-Leffler函数

单参数的Mittag-Leffler函数 [21]： $E_{\alpha }\left(z\right)={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}\frac{z^k}{\Gamma \left(\alpha k+1\right)}}$， $\alpha >0$

双参数的Mittag-Leffler函数 [21]： $E_{\alpha ,\beta }\left(z\right)={\displaystyle \underset{k=0}{\overset{\infty }{\sum }}\frac{z^k}{\Gamma \left(\alpha k+\beta \right)}}$， $\alpha >0$， $\beta >0$

注意到 $E_{\alpha ,1}\left(z\right)=E_{\alpha }\left(z\right)$，同时Mittag-Leffler函数是指数函数，双曲函数以及三角函数的推广， $E_{1,1}\left(z\right)={\text{e}}^z$， $E_{2,1}\left(z^2\right)=\cosh \left(z\right)$， $E_{2,1}\left(-z^2\right)=\cos \left(z\right)$， $E_{2,2}\left(-z^2\right)=\frac{\sin \left(z\right)}{z}$，对于Mittag-Leffler函数 [22]，以下公式成立：

${\displaystyle {\int }_0^x{t}^{\alpha -1}}{E}_{\alpha ,\alpha }\left(-a{t}^{\alpha }\right){\left(x-t\right)}^{\beta -1}{E}_{\alpha ,\beta }\left(-b{\left(x-t\right)}^{\alpha }\right)\text{d}t=\frac{E_{\alpha ,\beta }\left(-b{x}^{\alpha }\right)-E_{\alpha ,\beta }\left(-a{x}^{\alpha }\right)}{a-b}{x}^{\beta -1}$ (8)

其中 $\alpha \ne b$，同时，Mittag-Leffler函数的拉普拉斯变换对于求解分数阶微分方程 [23] 非常重要，

(9)

Srivastava和Tomovski提出了 $\left({\epsilon }_{a+;\alpha ,\beta }^{w;\gamma ,k}\phi \right)\left(x\right)$ 积分算子 [24]，定义如下：

$\left({\epsilon }_{a+;\alpha ,\beta }^{w;\gamma ,k}\phi \right)\left(x\right)={\displaystyle {\int }_a^x{\left(x-t\right)}^{\beta -1}}{E}_{\alpha ,\beta }^{\gamma ,k}\left(\omega {\left(x-t\right)}^{\alpha }\right)\phi \left(t\right)\text{d}t$ (10)

其中， $E_{\alpha ,\beta }^{\gamma ,k}\left(z\right)$ 是广义四参数Mittag-Leffler函数 [25]，它由以下形式给出：

$E_{\alpha ,\beta }^{\gamma ,k}\left(z\right)={\displaystyle \underset{n=0}{\overset{\infty }{\sum }}\frac{{\left(\gamma \right)}_{kn}}{\Gamma \left(\alpha n+\beta \right)}}\cdot \frac{z^n}{n!}$， $z,\beta ,\gamma \in C$ (11)

注意到 $E_{\alpha ,\beta }^{\gamma }\left(z\right)$ 是(11)式的特殊情况， $E_{\alpha ,\beta }^{\gamma ,k}\left(z\right)=E_{\alpha ,\beta }^{\gamma }\left(z\right)$ 以及 $E_{\alpha ,\beta }^{1,1}\left(z\right)=E_{\alpha ,\beta }\left(z\right)$。

3. 弦的一般方程

给定一根两端固定且拉紧的均匀的柔软的弦，其长度为L，密度为 $\rho$。假定：弦是均匀的，柔软的，且是作微小横振动的。选取如图1所示的坐标系，弦的平衡位置为轴，两端分别固定在和处， $u\left(x,t\right)$ 表示弦上横坐标为x的点在t时刻沿垂直x轴方向的位移。

在弦上任取一小弦弧 $\stackrel{⌢}{M_1{M}_2}$，设其弧长是 $\Delta s={\displaystyle {\int }_{x_1}^{x_2}\sqrt{1+{\left(\frac{\partial u}{\partial x}\right)}^2}}\text{d}x$，由假定，弦是作微小横振动的， $u_x^2$ 与1相比可忽略不计，于是 $\Delta s=x_2-x_1$，这样我们可以认为这段弦在振动过程中并未伸长，因此由胡克定律知，弦在每一点所受的张力在运动过程中保持不变，我们分别把在点M₁，M₂处的张力记作T₁，T₂，它们的方向分别沿着弦在点M₁，M₂处的切线方向。由假定，弦只作横向振动，因此张力在x轴方向的分量的代数和为零，即有 $T_1\cos {\alpha }_1-T_2\cos {\alpha }_2=0$，由于微小振动： ${\alpha }_1\to 0$， ${\alpha }_2\to 0$，因此有 $T_1=T_2$，也就是说，张力也不随地点而异，综上，张力是常数记为T₀。张力在方向的代数和为 $T_0\sin {\alpha }_1-T_0\sin {\alpha }_2=T_0\left(\sin {\alpha }_1-\sin {\alpha }_2\right)$，由于微小振动， $\sin {\alpha }_1\approx \tan {\alpha }_1={\frac{\partial u}{\partial x}|}_{x=x_1}$， $\sin {\alpha }_2\approx \tan {\alpha }_2={\frac{\partial u}{\partial x}|}_{x=x_2}$，

Figure 1. Vibration diagram of string

图1. 弦的振动示意图

于是有 $T_0\left({\frac{\partial u}{\partial x}|}_{x=x_1}-{\frac{\partial u}{\partial x}|}_{x=x_2}\right)=T_0{\frac{{\partial }^{2}u}{\partial x^2}|}_{x=\xi }\left(x_2-x_1\right)$， $\left(x_1<\xi <x_2\right)$。另一方面，由于所取弦段很小，其上每点的加速度相差也不会太大，因此可用其中一点 $\eta$ 处的加速度 ${\frac{{\partial }^{2}u}{\partial t^2}|}_{x=\eta }$ 代替，于是该小段弦的质量与加速度的乘积为： $\rho \left(x_2-x_1\right){\frac{{\partial }^{2}u}{\partial t^2}|}_{x=\eta }$， $\left(x_1<\eta <x_2\right)$，应用牛顿第二定律，有 $\rho {\frac{{\partial }^{2}u}{\partial t^2}|}_{x=\eta }\left(x_2-x_1\right)=T_0{\frac{{\partial }^{2}u}{\partial x^2}|}_{x=\xi }\left(x_2-x_1\right)$，于是 $\rho \frac{{\partial }^{2}u}{\partial t^2}=T_0\frac{{\partial }^{2}u}{\partial x^2}$，这就得到弦的自由横振动方程：

$\frac{{\partial }^{2}u}{\partial t^2}=a^2\frac{{\partial }^{2}u}{\partial x^2}$， $a=\sqrt{\frac{T_0}{\rho }}$ (12)

若还有外力 $F\left(x,t\right)$ 作用到弦上，则得到弦的强迫横振动方程：

$\frac{{\partial }^{2}u}{\partial t^2}=a^2\frac{{\partial }^{2}u}{\partial x^2}+f\left(x,t\right)$， $a=\sqrt{\frac{T_0}{\rho }}$， $f\left(x,t\right)=\frac{F\left(x,t\right)}{\rho }$。 (13)

4. 分数阶波动方程及其解

4.1. 时间分数阶波动方程

经典力学包括理论力学，弹性力学和流体力学。它们的基本力学关系如下所示：

弹性体的胡克定律： $F=kx$ (k表示弹性系数，x表示位移)

牛顿粘性流体的牛顿定律： $F=v\frac{\partial u}{\partial x}$ (v表示粘性系数，u表示速度，x表示空间变量)

刚体的牛顿第二定律： $F=m\frac{{\text{d}}^{2}x}{\text{d}{t}^2}$ (m表示质量，x表示位移，t表示时间)

我们注意到以上三个本构关系涉及了位移的0、1、2阶时间导数(位移、速度、加速度)。在描述多相介质(例如，软物质，粘弹性介质、复杂流体等流体、气体和固体的混合物)的力学行为时，这些经典模型并不准确。常用方法是将这些模型耦合成一个复杂模型，从分数阶微分方程建模的观点看，采用以下简单统一的分数阶模型来描述：

$F=\rho \frac{{\text{d}}^{\alpha }x}{\text{d}{t}^{\alpha }}$ $\left(0\le \alpha \le 2\right)$ (14)

我们考虑粘弹性弦的微小横振动，根据分数阶模型，粘弹性弦的加速度表示为 $\frac{{\text{d}}^{\alpha }u\left(x,t\right)}{\text{d}{t}^{\alpha }}$ $\left(1\le \alpha \le 2\right)$。和推导弹性弦的振动方程类似，唯一的区别在于弦的材料不同，因此两类不同材料的弦在振动时的加速度有所不同。应用牛顿第二定律有 $\rho {\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial t^{\alpha }}|}_{x=\eta }\left(x_2-x_1\right)=T_0{\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}|}_{x=\xi }\left(x_2-x_1\right)$，于是 $\rho \frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial t^{\alpha }}=T_0\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}$，因此粘弹性弦的自由横振动方程：

$\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial t^{\alpha }}=a^2\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}$， $a=\sqrt{\frac{T_0}{\rho }}$ (15)

同样，在外力的作用下，粘弹性弦的强迫振动方程为：

$\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial t^{\alpha }}=a^2\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}+f\left(x,t\right)$， $a=\sqrt{\frac{T_0}{\rho }}$， $f\left(x,t\right)=\frac{F\left(x,t\right)}{\rho }$ (16)

其中， $\frac{{\partial }^{\alpha }u\left(x,t\right)}{\partial t^{\alpha }}$ 代表Caputo型时间分数阶导数，其中， $1<\alpha \le 2$，也可用 ${}_0{}^{C}D{}_t^{\alpha }u\left(x,t\right)$ 表示，以此来描述粘弹性纳米弦振动时的时间记忆特性。

现考虑一根粘弹性弦，其长度为l，密度为 $\rho$，其振动方程为：

${}_0{}^{C}D{}_t^{\alpha }u\left(x,t\right)=a^2\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}$， $t>0$， $0\le x\le l$， $1<\alpha \le 2$ (17)

边界条件：

${u\left(x,t\right)|}_{x=0}=0$， ${u\left(x,t\right)|}_{x=l}=t^3$ (18)

初始条件：

${u\left(x,t\right)|}_{t=0}=0$， ${\frac{\partial u\left(x,t\right)}{\partial t}|}_{t=0}=x\left(l-x\right)$ (19)

4.2. 解过程

将 $u\left(x,t\right)$ 用以下形式表示：

$u\left(x,t\right)=U\left(x,t\right)+v\left(x,t\right)$ (20)

选择函数 $v\left(x,t\right)$ 满足边界条件(18)：

${v\left(x,t\right)|}_{x=0}=0$， ${v\left(x,t\right)|}_{x=l}=t^3$ (21)

容易得到：

$v\left(x,t\right)=\frac{x}{l}{t}^3$ (22)

根据式(18)和式(21)，可得到 $U\left(x,t\right)$ 的边界条件：

${U\left(x,t\right)|}_{x=0}=0$， ${U\left(x,t\right)|}_{x=l}=0$ (23)

由初始条件(19)和关系式(20)得：

${U\left(x,t\right)|}_{t=0}=0-{v\left(x,t\right)|}_{t=0}$， ${\frac{\partial U\left(x,t\right)}{\partial t}|}_{t=0}=x\left(l-x\right)-{\frac{\partial v\left(x,t\right)}{\partial t}|}_{t=0}$ (24)

使用替换：

$U\left(x,t\right)=U_1\left(x,t\right)+U_2\left(x,t\right)$ (25)

以及(17)和(20)，有

${}_0{}^{C}D{}_t^{\alpha }\left[U_1\left(x,t\right)+U_2\left(x,t\right)\right]=a^2\frac{{\partial }^2\left[U_1\left(x,t\right)+U_2\left(x,t\right)\right]}{\partial x^2}+\tilde{f}\left(x,t\right)$ (26)

其中 $\tilde{f}\left(x,t\right)=a^2\frac{{\partial }^{2}v\left(x,t\right)}{\partial x^2}-{}_0{}^{C}D{}_t^{\alpha }v\left(x,t\right)$。

关系式(26)中的函数可以按以下方式分离：

${}_0{}^{C}D{}_t^{\alpha }{U}_1\left(x,t\right)=a^2\frac{{\partial }^2{U}_1\left(x,t\right)}{\partial x^2}$ (27)

${U_1\left(x,t\right)|}_{x=0}=0$， ${U_1\left(x,t\right)|}_{x=l}=0$ (28)

${U_1\left(x,t\right)|}_{t=0}=0-{v\left(x,t\right)|}_{t=0}$， ${\frac{\partial U_1\left(x,t\right)}{\partial t}|}_{t=0}=x\left(l-x\right)-{\frac{\partial v\left(x,t\right)}{\partial t}|}_{t=0}$ (29)

以及

${}_0{}^{C}D{}_t^{\alpha }{U}_2\left(x,t\right)=a^2\frac{{\partial }^2{U}_2\left(x,t\right)}{\partial x^2}+\tilde{f}\left(x,t\right)$ (30)

${U_2\left(x,t\right)|}_{x=0}=0$， ${U_2\left(x,t\right)|}_{x=l}=0$ (31)

${U_2\left(x,t\right)|}_{t=0}=0$， ${\frac{\partial U_2\left(x,t\right)}{\partial t}|}_{t=0}=0$ (32)

对于 $U_1\left(x,t\right)$，对等式(27)使用变量分离法，将 $U_1\left(x,t\right)$ 表示为两个函数的乘积 $U_1\left(x,t\right)=X\left(x\right)T\left(t\right)$，可以获得以下微分方程：

${}_0{}^{C}D{}_t^{\alpha }T\left(t\right)+\lambda T\left(t\right)=0$ (33)

$a^2\frac{{\text{d}}^{2}X\left(x\right)}{\text{d}{x}^2}+\lambda X\left(x\right)=0$ (34)

其中 $\lambda$ 是分离常数，函数 $X\left(x\right)$ 满足以下边界条件：

${X\left(x\right)|}_{x=0}=0$， ${X\left(x\right)|}_{x=l}=0$ (35)

带有边界条件(35)的方程(34)表示Sturm-Liouville特征值问题，具有特征值 ${\lambda }_n$ 和特征函数 $X_n\left(x\right)$，且在 $L^2\left[0,L\right]$ 空间中满足下列条件：

${\displaystyle {\int }_0^l{X}_n^2\left(x\right)\text{d}x}={\Vert X_n\left(x\right)\Vert }^2$ (36)

方程(34)的通解为：

$X\left(x\right)=C\sin \frac{\sqrt{\lambda }}{a}x+D\cos \frac{\sqrt{\lambda }}{a}x$ (37)

且满足 $X\left(0\right)=0,X\left(l\right)=0$。

当 $X\left(0\right)=0$ 时，由于 $\cos 0=1$，要求 $D=0$，即 $X\left(x\right)=C\sin \frac{\sqrt{\lambda }}{a}x$。

当 $X\left(l\right)=0$ 时，即 $X\left(l\right)=C\sin \frac{\sqrt{\lambda }}{a}l=0$，要求 $C\ne 0$，只能 $\sin \frac{\sqrt{\lambda }}{a}l=0$，即 $\lambda$ 只能取特定的值，称为特征值，即 $\frac{\sqrt{\lambda }}{a}l=n\pi$，故可得特征值

${\lambda }_n=\frac{n^2{\pi }^2{a}^2}{l^2}$， $n=1,2,\cdots$ (38)

对每个 ${\lambda }_n$，其相应的解 $X_n\left(x\right)$ 称为特征函数

$X_n\left(x\right)=\sin \frac{n\pi x}{l}$， $n=1,2,\cdots$ (39)

对(33)式两边进行拉普拉斯变换得：

$s^{\alpha }\mathcal{L}\left[T_n\left(t\right)\right]-{\displaystyle \underset{k=0}{\overset{m-1}{\sum }}{T}_n^{\left(k\right)}\left(0_{+}\right)s^{\alpha -1-k}}+{\lambda }_n\mathcal{L}\left[T_n\left(t\right)\right]=0$ (40)

整理得：

$\mathcal{L}\left[T_n\left(t\right)\right]={\displaystyle \underset{k=0}{\overset{m-1}{\sum }}{T}_n^{\left(k\right)}\left(0_{+}\right)\frac{s^{\alpha -1-k}}{s^{\alpha }+{\lambda }_n}}$ (41)

再由Mittag-Leffler函数的拉普拉斯逆变换可以得到：

$T_n\left(t\right)={\displaystyle \underset{k=0}{\overset{m-1}{\sum }}{T}_n^{\left(k\right)}\left(0_{+}\right)E_{\alpha ,k+1}\left(-{\lambda }_n{t}^{\alpha }\right)}$ (42)

其中 $T_n^{\left(k\right)}\left(0\right)$ 为广义傅里叶系数，又因 $1<\alpha \le 2$，上式中 $m=2$， $k=0,1$，有

$T_n^{\left(0\right)}\left(0_{+}\right)=\frac{1}{{\Vert X_n\left(x\right)\Vert }^{\text{2}}}{\displaystyle {\int }_0^l\left[0-{v\left(x,t\right)|}_{t=0}\right]X_n\left(x\right)\text{d}x}=0$ (43)

$T_n^{\left(1\right)}\left(0_{+}\right)=\frac{1}{{\Vert X_n\left(x\right)\Vert }^{\text{2}}}{\displaystyle {\int }_0^l\left[x\left(l-x\right)-{\frac{\partial v\left(x,t\right)}{\partial t}|}_{t=0}\right]X_n\left(x\right)\text{d}x}=\frac{4l^2\left(1-{\left(-1\right)}^n\right)}{n^3{\pi }^3}$ (44)

得到：

$T_n\left(t\right)=\frac{4l^2\left(1-{\left(-1\right)}^n\right)}{n^3{\pi }^3}{E}_{\alpha ,2}\left(-{\lambda }_n{t}^{\alpha }\right)$ (45)

则方程(27)的解为：

$U_1\left(x,t\right)={\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\left[\frac{4l^2\left(1-{\left(-1\right)}^n\right)}{n^3{\pi }^3}\right]E_{\alpha ,2}\left(-\frac{n{}^2\pi {}^2{a}^2}{l^2}{t}^{\alpha }\right)}\sin \frac{n\pi x}{l}$。 (46)

这个和表示具有非零初始条件的纳米弦的自由振动问题的解，它是以特征函数集 $X_n\left(x\right)$ 为基础的函数 $U_1\left(x,t\right)$ 的广义傅里叶展开，且函数 $U_1\left(x,t\right)$ 满足与本征函数 $X_n\left(x\right)$ 相同的边界条件。

对于 $U_2\left(x,t\right)$，方程(30)的解同样可以通过使用特征函数 $X_n\left(x\right)$ 来表示，且函数 $U_2\left(x,t\right)$ 满足与本征函数 $X_n\left(x\right)$ 相同的边界条件：

$U_2\left(x,t\right)={\displaystyle \underset{n=1}{\overset{\infty }{\sum }}{u}_n\left(x\right)X_n\left(t\right)}$ (47)

以下面的形式展开函数 $\tilde{f}\left(x,t\right)$

$\tilde{f}\left(x,t\right)={\displaystyle \underset{n=1}{\overset{\infty }{\sum }}{\tilde{f}}_n}\left(t\right)X_n\left(x\right)$ (48)

其中 ${\tilde{f}}_n\left(t\right)=\frac{1}{{\Vert X_n\left(x\right)\Vert }^2}{\displaystyle {\int }_0^l\tilde{f}\left(x,t\right)X_n\left(x\right)\text{d}x}=\frac{2{\left(-1\right)}^{n+1}}{n\pi }{}_0{}^{C}D{}_t^{\alpha }{t}^3$。

通过使用关系式(30)，(34)以及(48)得：

${\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\left[{}_0{}^{C}D{}_t^{\alpha }{u}_n\left(t\right)+{\lambda }_n{u}_n\left(t\right)-{\tilde{f}}_n\left(t\right)\right]X_n\left(x\right)}=0$ (49)

上式满足：

${}_0{}^{C}D{}_t^{\alpha }{u}_n\left(t\right)+{\lambda }_n{u}_n\left(t\right)-{\tilde{f}}_n\left(t\right)=0$ (50)

对(50)式使用拉普拉斯变换：

$s^{\alpha }\mathcal{L}\left[u_n\left(t\right)\right]-{\displaystyle \underset{k=0}{\overset{m-1}{\sum }}{u}_n^{\left(k\right)}\left(0_{+}\right)s^{\alpha -1-k}}+{\lambda }_n\mathcal{L}\left[u_n\left(t\right)\right]-\mathcal{L}\left[{\tilde{f}}_n\left(t\right)\right]=0$ (51)

由于初始条件(32)可以得出， $u_n\left(0_{+}\right)=0$ 以及 $\frac{\partial u_n\left(0_{+}\right)}{\partial t}=0$，则式(51)变为以下形式：

$s^{\alpha }\mathcal{L}\left[u_n\left(t\right)\right]+{\lambda }_n\mathcal{L}\left[u_n\left(t\right)\right]-\mathcal{L}\left[{\tilde{f}}_n\left(t\right)\right]=0$ (52)

于是

$\mathcal{L}\left[u_n\left(t\right)\right]=\frac{1}{s^{\alpha }+{\lambda }_n}\mathcal{L}\left[{\tilde{f}}_n\left(t\right)\right]=\mathcal{L}\left[t^{\alpha -1}{E}_{\alpha ,\alpha }\left(-{\lambda }_n{t}^{\alpha }\right)\right]\mathcal{L}\left[{\tilde{f}}_n\left(t\right)\right]$ (53)

从上式可以看出 $u_n\left(t\right)$ 是两个函数的卷积，所以：

$u_n\left(t\right)={\displaystyle {\int }_0^t{\left(t-\tau \right)}^{\alpha -1}{E}_{\alpha ,\alpha }\left(-{\lambda }_n{\left(t-\tau \right)}^{\alpha }\right){\tilde{f}}_n\left(\tau \right)\text{d}\tau }$ (54)

根据式(47)以及Mittag-Leffler函数的性质，故可得方程(30)的解是：

$U_2\left(x,t\right)={\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\frac{2{\left(-1\right)}^{n+1}}{n\pi }\left[{\displaystyle {\int }_0^t{\left(t-\tau \right)}^{\alpha -1}{E}_{\alpha ,\alpha }\left(-\frac{n{}^2\pi {}^2{a}^2}{l^2}{\left(t-\tau \right)}^{\alpha }\right)}\left({}_0{}^{C}D{}_{\tau }^{\alpha }{\tau }^3\right)\text{d}\tau \right]\sin \frac{n\pi x}{l}}$ (55)

由积分算子式(27)以及Mittag-Leffler函数的性质，(55)式有以下形式：

$U_2\left(x,t\right)={\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\frac{2{\left(-1\right)}^{n+1}}{n\pi }\left({\epsilon }_{0+;\alpha ,\alpha }^{-{\lambda }_n;1,1}{}_0{}^{C}D{}_{\tau }^{\alpha }{\tau }^3\right)\left(\tau \right)\sin \frac{n\pi x}{l}}$ (56)

这个和表示具有非零初始条件以及在外力下的纳米弦的强迫振动问题的解。最后，在 $1<\alpha \le 2$ 的情况下，通过关系式(22)，(46)以及(56)得到方程(17)在非齐次边界条件下的解析解：

$\begin{array}{c}u\left(x,t\right)=v\left(x,t\right)+U_1\left(x,t\right)+U_2\left(x,t\right)\\ =\frac{x}{l}{t}^2+{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\frac{4l^2\left(1-{\left(-1\right)}^n\right)}{n^3{\pi }^3}{E}_{\alpha ,2}\left(-\frac{n{}^2\pi {}^2{a}^2}{l^2}{t}^{\alpha }\right)}\sin \frac{n\pi x}{l}\\ +{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\frac{2{\left(-1\right)}^{n+1}}{n\pi }\left({\epsilon }_{0+;\alpha ,\alpha }^{-{\lambda }_n;1,1}{}_0{}^{C}D{}_{\tau }^{\alpha }{\tau }^3\right)\left(\tau \right)\sin \frac{n\pi x}{l}}\end{array}$ (57)

5. 算例分析

例一，将方程(17)的非齐次边界条件 ${u\left(x,t\right)|}_{x=l}=t^3$ 换作齐次边界条件 ${u\left(x,t\right)|}_{x=l}=0$，其余条件不变，它的解为：

$u\left(x,t\right)={\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\frac{4l^2\left(1-{\left(-1\right)}^n\right)}{n^3{\pi }^3}{E}_{\alpha ,2}\left(-\frac{n{}^2\pi {}^2{a}^2}{l^2}{t}^{\alpha }\right)}\sin \frac{n\pi x}{l}$ (58)

注意到 $E_{2,1}\left(-z^2\right)=\cos \left(z\right)$，(58)式还有如下形式：

$u\left(x,t\right)={\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\frac{4l^2\left(1-{\left(-1\right)}^n\right)}{n^3{\pi }^3}\cos \left(\frac{n\pi a}{l}t\right)}\sin \frac{n\pi x}{l}$ (59)

对于 $a=2$，解(58)，(59)由图2给出。

Figure 2. For $x\in \left[0,1\right]$. $a=2$, the graphical representation of the solution

图2. 对于 $x\in \left[0,1\right]$， $a=2$，解的图形表示

从图3可以看出对于 $x\in \left[0,1\right]$， $t=0.15$， $a=2$， $\alpha =1.85$ (点线)， $\alpha =1.9$ (点虚线)， $\alpha =1.95$ (虚线林)， $\alpha =2$ (实线)，随着时间分数阶阶数 $\alpha$ 的增大， $u\left(x,t\right)$ 的图像越接近。

例二，将方程(17)的非齐次边界条件 ${u\left(x,t\right)|}_{x=l}=t^3$ 换作齐次边界条件 ${u\left(x,t\right)|}_{x=l}=0$，同时给方程(17)再加上外力项 $f\left(x,t\right)=t^{\beta -1}{E}_{\alpha ,\beta }\left(-b{t}^{\alpha }\right)$， $1<\beta \le 2$，b为常数，其余条件不变，得出的解为(图4)：

Figure 3. For $x\in \left[0,1\right]$, $t=0.15$, $a=2$, the graphical representation of the solution

图3. 对于 $x\in \left[0,1\right]$， $t=0.15$， $a=2$，解的图形表示

$\begin{array}{c}u\left(x,t\right)={\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\frac{4l^2\left(1-{\left(-1\right)}^n\right)}{n^3{\pi }^3}{E}_{\alpha ,2}\left(-\frac{n{}^2\pi {}^2{a}^2}{l{}^2}{t}^{\alpha }\right)}\sin \frac{n\pi x}{l}\\ +2t^{\beta -1}\sin l{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\frac{1-{\left(-1\right)}^n}{n\pi }}\frac{E_{\alpha ,\beta }\left(-b{t}^{\alpha }\right)-E_{\alpha ,\beta }\left(-\frac{n{}^2\pi {}^2{a}^2}{l^2}{t}^{\alpha }\right)}{\frac{n{}^2\pi {}^2{a}^2}{l^2}-b}\sin \frac{n\pi x}{l}\end{array}$ (60)

Figure 4. For $x\in \left[0,2\right]$, $a=2$, $\beta =1.8$, $a=b=2$, the graphical representation of the solution

图4. 对于 $x\in \left[0,2\right]$， $a=2$， $\beta =1.8$， $a=b=2$，解的图形表示

Figure 5. For $x\in \left[0,2\right]$, $a=2$, $\beta =1.8$, $a=b=2$, the graphical representation of the solution

图5. 对于 $x\in \left[0,2\right]$， $a=2$， $\beta =1.8$， $a=b=2$，解的图形表示

从图5可以看出对于 $x\in \left[0,2\right]$， $a=2$， $\beta =1.8$， $a=b=2$， $t=0.15$ (实线)， $t=0.2$ (虚线林)， $t=0.25$ (点虚线)， $t=0.3$ (点线)，随着时间t的增大， $u\left(x,t\right)$ 的图像在逐渐减小。

例三：对于方程(17)的解(57)：

$\begin{array}{c}u\left(x,t\right)=v\left(x,t\right)+U_1\left(x,t\right)+U_2\left(x,t\right)\\ =\frac{x}{l}{t}^2+{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\frac{4l^2\left(1-{\left(-1\right)}^n\right)}{n^3{\pi }^3}{E}_{\alpha ,2}\left(-\frac{n{}^2\pi {}^2{a}^2}{l^2}{t}^{\alpha }\right)}\sin \frac{n\pi x}{l}\\ +{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\frac{2{\left(-1\right)}^{n+1}}{n\pi }\left[{\displaystyle {\int }_0^t{\left(t-\tau \right)}^{\alpha -1}{E}_{\alpha ,\alpha }\left(-\frac{n{}^2\pi {}^2{a}^2}{l^2}{\left(t-\tau \right)}^{\alpha }\right)}\left({}_0{}^{C}D{}_t^{\alpha }{\tau }^3\right)\text{d}\tau \right]\sin \frac{n\pi x}{l}}\end{array}$

由Caputo分数阶导数的定义，有 ${}_0{}^{C}D{}_t^{\alpha }{t}^p=\frac{\Gamma \left(p+1\right)}{\Gamma \left(p-\alpha +1\right)}{t}^{p-\alpha }$， $n-1<\alpha <n$， $p>n-1$， $p\in R$，所以 ${}_0{}^{C}D{}_t^{\alpha }{\tau }^3=\frac{\Gamma \left(3+1\right)}{\Gamma \left(3-\alpha +1\right)}{\tau }^{3-\alpha }$，所以解的第三项有以下形式：

$\frac{\Gamma \left(3+1\right)}{\Gamma \left(3-\alpha +1\right)}{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\frac{2{\left(-1\right)}^{n+1}}{n\pi }}\left[{\displaystyle {\int }_0^t{\left(t-\tau \right)}^{\alpha -1}{E}_{\alpha ,\alpha }\left(-\frac{n{}^2\pi {}^2{a}^2}{l^2}{\left(t-\tau \right)}^{\alpha }\right){\tau }^{3-\alpha }\text{d}\tau }\right]\sin \frac{n\pi x}{l}$ (61)

当 $a=2$ 时，做一定的变换得到：

${\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\frac{6{\left(-1\right)}^{n+1}}{n\pi }}\left[t^3{E}_{2,3}\left(-\frac{n{}^2\pi {}^2{a}^2}{l^2}{t}^2\right)-{\displaystyle \underset{k=0}{\overset{\infty }{\sum }}\frac{{\left(-\frac{n{}^2\pi {}^2{a}^2}{l^2}\right)}^k{t}^{2k+3}}{\left(2k+3\right)\Gamma \left(2k+2\right)}}\right]\sin \frac{n\pi x}{l}$ (62)

带入解(57)，得到

$\begin{array}{c}u\left(x,t\right)=v\left(x,t\right)+U_1\left(x,t\right)+U_2\left(x,t\right)\\ =\frac{x}{l}{t}^2+{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\frac{4l^2\left(1-{\left(-1\right)}^n\right)}{n^3{\pi }^3}{E}_{\alpha ,2}\left(-\frac{n{}^2\pi {}^2{a}^2}{l^2}{t}^{\alpha }\right)}\sin \frac{n\pi x}{l}\\ +{\displaystyle \underset{n=1}{\overset{\infty }{\sum }}\frac{6{\left(-1\right)}^{n+1}}{n\pi }}\left[t^3{E}_{2,3}\left(-\frac{n{}^2\pi {}^2{a}^2}{l^2}{t}^2\right)-{\displaystyle \underset{k=0}{\overset{\infty }{\sum }}\frac{{\left(-\frac{n{}^2\pi {}^2{a}^2}{l^2}\right)}^k{t}^{2k+3}}{\left(2k+3\right)\Gamma \left(2k+2\right)}}\right]\sin \frac{n\pi x}{l}\end{array}$ (63)

由图6给出其图形表示。

(a) (b) (c) (d)

Figure 6. For $x\in \left[0,1\right]$, $a=2$, $a=1$, the graphical representation of the solution

图6. 对于 $x\in \left[0,1\right]$， $a=2$， $a=1$，解的图形表示

从图6可以看出对于 $x\in \left[0,1\right]$， $a=2$， $a=1$，图6(a)~(d)分别表示 $t\in \left[0,1\right]$， $t\in \left[0,0.1\right]$， $t\in \left[0,0.01\right]$， $t\in \left[0,0.001\right]$ 时解的图形，随着时间t的增大 $u\left(x,t\right)$ 震荡的幅度越大。

6. 结论

本文研究纳米弦振动在不同边界条件的时间分数阶波动方程。利用分离变量法和拉普拉斯变换法求解该波动方程，通过对不同边界和外力以及初始条件下的算例进行模拟计算，并对一些特殊的初值条件和边界条件的运行结果进行具体分析，得出了解决问题的仿真图形结果，随着分数阶阶数 $\alpha$ 和时间t的变化， $u\left(x,t\right)$ 有着不同的变化。

基金项目

国家自然基金项目(11862014)。

NOTES

^*通讯作者。