Mindlin 矩形微板的热弹性阻尼解析解1)

马航空 周晨阳 李世荣

(扬州大学建筑科学与工程学院，江苏扬州 225127)

引言

微/纳谐振器作为微/纳电机系统(MEMS/NEMS)中的基本器件被广泛用作频率控制器件、逻辑开关、传感器、驱动器、陀螺仪、能量收集器等，在生物医学、自动控制、航空电子设备、国防科技等众多领域有着十分广泛的应用.为了设计和制造高性能的微/纳电机系统，需要谐振器在工作时耗能越小越好，即具有更高的品质因子.然而，在微/纳电机系统中不可避免地存在各种耗能机制.除了空气阻尼、支撑阻尼、挤压模阻尼等外部能量耗散外还存在内部能量耗散，例如热弹性阻尼(thermoelastic damping,TED)，或内摩擦[1].实践证明，可以通过完美的设计和制造工艺最大限度地降低或消除外阻尼引起的耗量能散.但是，TED 是系统所固有的能耗机制，是在结构周期振动时材料内部的热弹性耦合变形引起的内部阻尼，它不能通过改善外部条件而消除.因此，TED 有时会成为微电机系统的主要能耗形式，或影响品质因子的重要因素.要在理论上分析微/纳结构在振动中的热弹性阻尼，就需要基于热--弹耦合动力学理论建立结构振动的数学模型，通过联立求解热--弹耦合的运动方程和热传导方程获得系统的动力响应，从而求得表征热弹性阻尼的品质因子.TED 与微结构的材料性质、几何尺寸、支承条件以及环境温度等诸多因素密切相关.因此，定量地分析和预测TED 对高品质微/纳谐振器的研究和设计具有重要意义.

大多数谐振器力学模型可以简化为弹性梁或弹性板.最早关于谐振器热弹性阻尼的理论研究当属Zener 的工作[1].早在1938 年，Zener 首先基于Euler-Bernoulli (E-B)梁理论和单相耦合的准一维热传导方程求得了矩形截面微梁谐振器的TED近似解析解[1]，被称为著名的Zener 公式.2000 年，在Zener工作的基础上Lifshitz 和Roukes(L-R)[2]首次给出了E-B 微梁的准一维热传导方程的解析解，进而求解振动方程的特征值问题，获得了更为精确的矩形截面微/纳梁谐振器逆品质因子和频移的解析解，被称为著名的L-R 公式.Zener 和L-R 的杰出工作对后来关于微结构的TED 的进一步研究具有重要的参考和借鉴作用，他们的成果在后来研究者的论文中被高频次引用.研究者采用不同的结构变形理论、热传导理论和分析方法开展了关于微/纳梁板谐振器的热弹性阻尼研究.关于逆品质因子的计算主要采用能量法[1]和复频率法[2].

基于经典板理论，文献[3-15]报道了均匀各向同性材料弹性薄板谐振器的TED[3-15].Nayfeh 和Youns[3]采用摄动法求解了引入刚度阻尼的振动方程，获得了受静电载荷和残余面内应力共同作用的矩形微板TED 的近似解析解.Sun 等[3-5]、Ali 和Mohammadi[6]采用复频率法分别研究了微圆/环板的TED，给出了L-R 形式的解析解.Salajeghe 等[7]在运动方程的建立中考虑了变形的几何非线性，定量分析了轴对称振动微圆板谐振器的振幅对热弹性阻尼的影响规律.Li 等[8]采用准一维热传导方程，通过计算一个周期内弹性板损失的能量与结构内存储的最大弹性势能之比求得了矩形和圆形薄板谐振器L-R 形式的逆品质因子解析解.在此基础上Fang 等[9-10]分别采用二维和三维热传导方程研究了轴对称自由振动的圆板和矩形板谐振器的TED,分析了热传导方程的简化对热弹性阻尼值的预测精度的影响.考虑热传导过程中温度热流运动的延滞效应(或波动效应),文献[11-14] 利用广义热传导理论(或非傅里叶热传导理论)研究了微/纳板谐振器的热弹性耦合振动响应.Sharma 等[11]基于Lord-Shulman 广义热传导理论和经典板理论研究了非轴对称自由振动圆板谐振器的热弹性阻尼，分析了延滞时间参数对热弹性阻尼的影响.在此基础上，Sharma 和Grover[12]还进一步研究了由于空隙率变化而引起的延滞效应对板式谐振器的影响.Guo 等[13]采用双向延滞广义热传导模型研究了圆板谐振器的热弹性阻尼，采用复频率法获得逆品质因子的L-R 形式解析解.数值结果发现，在广义热传导理论下逆品质因子随着厚度的减小出现极小值.Grover 等[14]基于Kelvin-Voigt 材料模型下的广义热黏弹性理论研究了均匀微圆板谐振器的阻尼特性，定量地分析了机械松弛时间和热松弛时间对逆品质因子的影响.

考虑材料性质沿单轴方向阶梯变化，Bishop和Kinra[15-16]最早采用能量法开展了关于复合材料层合微/纳谐振器的热弹性阻尼研究.在表面边绝热条件和内部界面处的非完善协调条件下，求得了准一维热传导方程的解析解，并通过计算一个振动周期内由不可逆传热产生的总热量与弹性总弹性势能之比给出了谐振器的逆品质因子[15].并定量地分析了对称铺设的三层矩形微板的弹性热动阻尼特性[16].Sun 等[17]基于准一维热传导方程，采用复频率法研究了对称铺设的三层微圆板在轴对称自由振动下的热弹性阻尼.采用与文献[16]中相同的能量方法，忽略温度梯度在面内的变化，Zuo 等[18]研究了双层微板的热弹性阻尼，给出了周边夹紧矩形板和圆板的逆品质因子解析解.考虑温度梯度在面内的变化，Liu 等[19-20]采用格林函数法分别求解了二维和三维热传导方程，利用能量法获得了双层圆板和非对称铺设三层矩形板在给定边界条件下的逆品质因子的级数形式解析解.这里需要说明的是，文献[18-20]通过引入物理中面的概念消去了几何中面的面内位移，简化了几何方程.但是，当材料性质变化关于几何中面非对称时板内将会产生热薄膜力，而热薄膜力对热弹性阻尼具有贡献[29].显然，在上述文献中忽略了热薄膜力对热弹性阻尼的影响.

近年来，亦有少量的文献研究了材料性质沿厚度连续变化的功能梯度材料(functionally graded material,FGM)微/纳梁板结构的热弹性阻尼.Azizi 等[21]研究了由硅和压电材料为组分材料复合而成的两端夹紧的矩形截面FGM微梁横向自由振动时的TED.采用Galerkin 法求得了FGM 微梁的复频率进而获得反映热弹性阻尼的逆品质因子.数据结果表明可以通过压电材料组分的合理调节来降低系统的TED.基于一阶剪切变形理论和单向耦合的准一维热传导方程，Emami 和Alibeigloo[22]研究了材料性质沿横向幂函数变化的FGM 微梁的TED.将热传导方程的系数和温度场同时展开成关于横向坐标的泰勒级数形式，获得了变系数热传导方程的级数形式的解析解，进而获得了两端简支FGM 微梁的TED.在此基础上，文献[22]的作者中又基于一阶剪切变形板理论和修正的应变梯度理论，研究了四边简支功能梯度中厚度矩形板的热弹性阻尼[23].但是，通过研究发现其中给出的关于均匀材料方板的热弹性阻尼的计算结果是错误的.由此无法肯定其中关于功能梯度材料微板热弹性阻尼数值结果的可靠性.考虑材料性质在厚度方向按指数函数分布特殊情形，Zhong 等[24]采用解析方法获得了功能梯度E-B 微梁的热弹性阻尼解析解.但是，作者在应变的计算中却直接忽略了由于材料性质在横向的非均匀性而导致的拉--弯耦合变形.许新等[25-26]基于单向耦合的准一维热传导方程研究了材料性质在横向任意连续变化的FGM E-B 微梁的热弹性阻尼.发展了求解横向非均匀微梁的复杂变系数热传导方程的分层均匀化方法，从而获得FGM 微梁热弹性阻尼的半解析解,其中精确地考虑了热轴力对TED 的影响.最近又将分层均匀化方法成功地应用于功能梯度微板的热弹性阻尼的求解[27-28].研究了材料梯度指数、环境温度、边厚比、振动模态等对热弹性阻尼的影响规律，首次定量地分析了热薄膜力对热弹性阻尼的影响程度.

文献调研表明，截至目前关于均匀微板热弹性阻尼的解析解是基于Kirchhoff经典弹性板理论推导出的，其中忽略了横向剪切变形对振动响应的影响.本文基于考虑一阶剪切变形效应的Mindlin 中厚板理论和单向耦合热传导理论建立四边简支均匀材料微板的热--弹耦合自由振动的控制微分方程.消去两个独立的转角变量将振动方程转化为只包含挠度函数和热弯矩的四阶偏微分方程.忽略温度梯度在面内的变化，在上下表面绝热边界条件下求得了用挠度函数表示的热传导方程的解析解.从而将包含热弯曲内力的结构振动方程转化为只包含挠度振幅的偏微分方程.利用特征值问题的数学相似性,在四边简支条件下获得了用均匀无阻尼Kirchhoff微板的固有频率表示的Mindlin 矩形微板的复频率，从而由复频率法给出了反映热弹性阻尼水平的逆品质因子解析解.通过数值结果定量地分析了剪切变形对热弹性阻尼值的影响程度.

1 问题的控制方程

1.1 运动方程

考虑长度a、宽度为b、厚度为h的等厚度矩形微/纳板.选取直角坐标系，板在空间的区域定义为：0 ＜x＜a，0 ＜y＜b，-h/2 ＜z＜h/2 (见图1).基于Mindlin 板理论[29-33]，板内任意一点的位移分量可表示为

其中，t为时间；u，v和w分别为板内任意一点沿x，y和z轴方向的位移分量；φx和φy为中面法线分别绕y轴和x轴的转角.

图1 微板的几何尺寸和坐标示意图Fig.1 Schematic illustration of a micro plate and the coordinate system

将式(1)代入弹性力学的几何方程可得板的应变场

根据胡克定律得到应力场

其中，θ(x,y,z,t)=T(x,y,z,t)-T0为变温场，T为瞬态温度，T0为初始温度.E，v和α 分别为弹性模量、泊松比和热膨胀系数.这里的温度场是由热--弹耦合振动引起的.

将弹性力学空间问题应力形式的运动方程在厚度进行静力等效，并利用上下表面的应力边界条件可得用等效内力和内力矩表示的Mindlin 板[29-33]自由振动微分方程

其中Ii为等效惯性参数，具体表示为

ρ 为板的质量密度.式(4)～式(6)中的等效内力和内力矩可分别表示为

将式(8)代入式(4)和式(5)，分别关于坐标x和y求导后相加，利用式(6)可得

将式(8d)和式(8e)代入式(6)可得

最后，将式(11)代入式(9),可得只用横向位移w0表示的运动方程

上式中热弯矩由热--弹耦合产生的变温场确定.

1.2 热传导方程

根据单向耦合的热弹性动力学理论，可给出微/纳板的准一维热传导方程[3-8]

其中，κ，C分别是热传导系数和比热，e是体积应变.Mindlin 板的体积应变可表示为

将式(14)代入式(13)，可得单向耦合的准一维热传导方程

2 热弹性阻尼的求解

2.1 自由振动响应

假设位移场合温度场的动态响应同步，则系统的热--弹耦合自由振动响应可假设为调和形式

可得无量纲控制方程

其中

其中结构振动方程通过无量纲的热弯矩振幅mT与热传导方程耦合.

2.2 热传导方程的求解

方程(20)的通解为

微谐振器一般在恒温真空环境内工作,可认为没有外界的热流输入.而由于热--力耦合振动在谐振器内部产生的变温场幅度很小,因而热流强度也很低.于是可假设在微板的上下表面没有热流通过.这样，在微板的上下表面通常给出绝热边界条件[2-28]

利用边界条件(24)可确定系数A=qK/[pcos(p/2)]和B=0.于是得到特解

2.3 热弹性阻尼

将式(25)代入热弯矩的计算公式可得

将式(26)代入式(19),并利用式(11)的无量纲形式

可得只用无量纲挠度表示的结构振动微分方程

微分方程(29)可记为标准形式

其中

另外，由弹性薄板理论可得不考虑热弹性阻尼的Kirchhoff板的振动微分方程[29]

对于四边简支矩形板，Kirchhoff薄板理论下方程(32)的边界条件可表示为

可以证明(具体证明见附录)，对于四边简支Mindlin 矩形板，振动微分方程(30)的边界条件亦可表示为

由微分方程边值问题(30)，(34)与问题(32)，(33)之间的相似性，可得它们的特征值之间的关系

于是，将式(31)代入式(35)，可得特征方程

其中R=μ00-(μ2+csf).由式(27)知f(Ω)是关于频率Ω 的超越函数，因此方程(36)具有很强非线性.为了简化求解,可近似地令f(Ω)=f(Ω0)[2-14]，其中Ω0为等温(不考虑热弹性阻尼)Mindlin 微板的固有频率.于是，由式(36)可容易求得用无阻尼Kirchhoff板的频率表示的具有热弹性阻尼的Mindlin 微板的复频率

其中，四边简支等温Kirchhoff矩形微板的固有频率的解析解为

在式(37)中令f=0 可得

式(39)给出了等温Mindlin 微板的固有频率与相应的Kirchhoff板的固有频率之间的解析关系.将文献[29]中的式(12.2.16)进行无量纲运算后所得结果与式(39)完全相同.

先由式(39)求得Ω0，然后代入式(27)可得到f=f(Ω0),最后将其代入式(37)即可得到热--弹耦合自由振动Mindlin 微板的复频率.根据复频率法[2-8]即可得到用逆品质因子表示的Mindlin 微板的热弹性阻尼解析解

其中Re(Ω)和Im(Ω)分别是复频率的实部和虚部.

可以证明，对于简支多边形微板式(33)和式(34)仍然成立[29].因此，式(37)和式(40)也可以计算直边多边形中厚度微板的复频率和热弹性阻尼.

如果在式(36)中令μ2=cs=0，f=f()，则可得到Kirchhoff微板的复频率[3-8]

利用式(27)和式(41)可得文献[8]中 Kirchhoff微板热弹性阻尼L-R[2]形式的解析解

3 数值结果与讨论

数值计算时微板的材料分别选取碳化硅(SiC)、氮化硅(Si3N4)、金属铝(Al)和金属镍(Ni)[23,29],具体材料参数见表1.环境温度为T0=300 K.

表1 微板的材料性质(T0=300 K)[23,25]Table 1 Material properties of the micro plate(T0=300 K)[23,25]

首先，针对碳化硅(SiC)方板，在表2 中给出了在不同振动模态下，基于Mindlin 板理论得到的热弹性阻尼(40)与文献[8]中采用经典理论所得到的解析解(42)的比较.由表中的数据可见，随着微板的边/厚比的减小(或边长的减小)以及模态阶数的增大Mindlin 微板的逆品质因子与Kirchhoff微板的逆品质因子之间的相对误差Error=100%(-逐渐增大.在a/h=10 和基频振动下，两种理论预测结果之间的相对误差仅为2.66%,而当微板以高阶模态(r,s)=(3,3)振动时，相对误差却达到了8.94%.对应同样的几何参数、材料性质参数和振动模态，在表3 中给出了文献[23]中的数值结果.与表2 对比后可见，文献[23]中的结果与本文结果相差甚远.而且在表3 中利用文献[8]中的Kirchhoff微板的解析公式(42)计算得到的数值结果与本文利用同样公式计算结果亦有很大差别.特别是在边厚比a/h=10、振动模态为(3,3)的条件下，本文所得Mindlin 板和Kirchhoff板的TED 分别为8.04×10-5和8.83×10-5，相对误差为8.94%；而在文献[23]中的相应结果则分别为6.4×10-5和1.30×10-4，两种理论的预测值之间的相对误差竟然达到了50.7%.而对于长厚比a/h=5 的正方形微板,在振动模态为(3,3)的条件下，本文预测的两种理论下的热弹性阻尼之间的相对误差也仅为28.4%.

表2 两种板理论下正方形形陶瓷(SiC)微板的热弹性阻尼(Q-1×104)比较(a=b,h=1 μm)Table 2 Comparison of TED(Q-1×104)of the square micro plate of ceramic(SiC)based on both of the Kirchhoffand Mindlin plate theories(a=b,h=1 μm)

表3 具有不同边/厚比的正方形微板在前六阶模态下的热弹性阻尼比较TED(Q-1×103)(h=1 μm)[23]Table 3 Comparison of TED(Q-1×103)of square microplate with different side-to-thickness ratios for the first six modes(h=1 μm)[23]

为了进一步验证本文解析解的可靠性，选择a=b=100 μm 的陶瓷(SiC)微板，在图2 中分别给出了由式(34)和式(40)预测的TED 和采用有限元法得到数值结果的比较.有限元法是基于全耦合热弹性动力学理论、采用空间六面体单元获得的.由图中结果可见，在板厚较小时二者吻合很好.随着厚度的增大，有限元解答与本文解析解稍有偏离.这是因为在板理论中，将侧面上的应力边界条件沿着厚度静力等效为几何中面周边上的合力矩边界条件，根据圣维南原理随着厚度的增大板的自然边界条件与空间弹性力学的应力边界条件的差别将会增大.实际上，位移边界条件的差别也在增大.

图2 分别由本文方法和有限元法预测的正方形陶瓷(SiC)微板的热弹性阻尼值的比较Fig.2 Comparison between the values of TED of square ceramic micro plate by FEM and present approach,respectively

为了定量地分析剪切变形对热弹性阻尼预测值的影响，给定不同边/厚比a/h，在图2 中绘出了正方形金属(Ni)板谐振器在基频振动时的热弹性阻尼随板厚的连续变化曲线.并在表4 中列出了图2 中热弹性阻尼的最大值Q-1max和对应的临界厚度hcr.数值结果表明，Kirchhoff板理论的热弹性阻尼预测值始终大于Mindlin 板理论的预测值.两种理论预测的热弹性阻尼之间的差值在临界厚度附近十分显著.另外，随着微板的边/厚比的增大，Mindlin 微板的热弹性阻尼最大值单调增加，而Kirchhoff微板的热弹性阻尼最大值却保持不变.从表4 中数据可见，Mindlin 微板的临界厚度大于Kirchhoff微板的临界厚度.为了更加清晰地反映两种理论预测值的差值的变化，在图3 中绘出了差值-随厚度的连续变化曲线.图中结果再次清楚地表明在临界厚度附近差值变化最大.但是，在不考虑热弹性阻尼的情况下由式(39)可知Mindlin 板的无量纲固有频率只与边厚比有关.在给定边厚比后无量纲频率为常数.因此，两种理论预测的等温板的无量纲故有频率差值为常数.

表4 具有不同边/厚比正方形微板(Ni)的最大热弹性阻尼和相应的临界厚度(一阶模态)Table 4 The maximum TED and related critical thickness of square micro plate(Ni)for different side-to-thickness ratios(the first mode)

图3 两种板理论下正方形金属(Ni)微板的热弹性阻尼随厚度的变化曲线(一阶模态)Fig.3 Curves of the TED versus the plate thickness of a square micro metal plate(Ni)based on the two plate theories(in the first mode)

图4 中绘出了表1 中所列4 种材料的正方形微板在一阶模态振动时的热弹性阻尼随厚度的连续变化曲线.其中给出了两种板理论预测结果的对比.从图中可见，金属微板的热弹性阻尼明显大于陶瓷微板的热弹性阻尼.并且两种理论下金属微板的最大热弹性阻尼之间的差值比陶瓷微板显著.图5 给出了a/h=10 时，具有不同长宽比的矩形陶瓷(SiC)微板在基频振动下的TED 与厚度的关系曲线.结果表明，随着长宽比的增大(相当于弯曲刚度增大)两种板理论预测的热弹性阻尼的最大值之间的差值显著增大，而经典理论下的最大值却保持不变.图6 则为a/h=10 的中等厚度板金属(Al)微板在前四阶振动模态下的热弹性阻尼随厚度的变化曲线.其中的变化规律与图5 相似.随着振动模态阶数的增大，板的固有频率增大.这相当于弯曲刚度在增加，从而导致剪切变形对热弹性阻尼的影响加大.

图4 给定不同边厚比时，两种板理论预测的正方形金属(Ni)微板的热弹性阻尼差值随厚度的变化曲线(一阶模态)Fig.4 Difference between the values of TED of a square metal(Ni)micro plate evaluated by the two plate theories varying with the thickness for some specified values of a/h(in the first mode)

图5 一阶模态下不同材料正方形微板的热弹性阻尼与板厚之间的关系曲线Fig.5 Relationship curves of TED versus the plate thickness of the square micro plates of different materials in the first mode

图6 具有不同长宽比(a/b)的矩形陶瓷(SiC)微板的热弹性阻尼与板厚之间的关系曲线(一阶模态)Fig.6 TED versus the plate thickness in ceramic(SiC)micro rectangular plate with different length to width ratio(in the first mode)

图7 不同振动模态下金属(Al)微板的热弹性阻尼随厚度的关系曲线(a/h=10)Fig.7 TED versus the plate thickness in a square metal(Al)micro plate in different vibration modes(a/h=10)

4 结论

基于一阶剪切变形理论和单向耦合准一维热传导理论，采用解析方法研究了四边简支中厚度矩形微板谐振器的热弹性阻尼.利用微分方程边值问题之间的相似性，推导出了用无阻尼Kirchhoff微板的固有频率表示的Mindlin 微板热--弹耦合自由振动复频率的解析解，从而采用复频率法给出了表征Mindlin 微板热弹性阻尼的逆品质因子解析解.分别选定两种金属和两种陶瓷材料微板，通过参数变化的数值结果定量地分析了剪切变形对热弹性阻尼预测值的影响规律.数值结果表明，经典板理论预测的热弹性阻尼值大于一阶剪切变形板理论的预测值.通过研究还发现，两种理论的预测的热弹性阻尼值的差别在微板的临界厚度附近最为显著.

致谢根据审稿专家的建议,在修改稿中增加了本文解析解与有限元法所得数值结果的比较.扬州大学建筑科学与工程学院的张志超博士和他指导的硕士研究生曹静同学在有限元数值计算中给予作者重要的技术帮助.在此表示感谢.