功能梯度材料明德林矩形微板的热弹性阻尼1)

2022-07-10 13:13李世荣

力学学报 2022年6期

李世荣

(南通理工学院土木工程系,江苏南通 226002)

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

引言

谐振器作为微/纳电机系统中的重要器件被广泛用作传感器、驱动器、陀螺仪、能量收集器、频率控制器、逻辑开关等.而大多数谐振器的力学模型可以简化为弹性微梁或微板结构.高性能的微/纳电机系统需要谐振器在工作时具有低耗能或高品质因子的优点.然而,谐振器在工作中不可避免地存在各种耗能机制.例如热弹性阻尼或内摩擦[1]就是一种主要的内耗能机制.热弹性阻尼是结构在一个振动周期内材料内部的热弹性耦合变形引起的内部能耗,它不能通过外部条件的改善而消除.因此,热弹性阻尼有时将会成为微电机系统的主要能耗形式.因此,精确地分析和预测热弹性阻尼对高品质微/纳谐振器的研究和设计具有重要意义.

Zener[1]首先采用能量方法和单向耦合的热传导方程给出了细长矩形截面微梁谐振器的逆品质因子解析解.Bishop 和 Kinra[2-3]首先把Zener[1]关于均匀材料微梁谐振器热弹性阻尼分析的能量法推广到了具有非完善界面复合材料层合微纳梁/板谐振器的热弹性阻尼研究,并通过计算一个振动周期内由不可逆传热产生的总热量与弹性总弹性势能之比给出了微层合结构的逆品质因子的解析解.与上述能量法不同,Lifshitz 和Roukes[4]通过联立求解单向耦合的结构振动方程和热传导方程获得了细长微梁热弹耦合自由振动的复频率,首次采用复频率法获得了微梁谐振器热弹性阻尼的解析解.此后,能量法[1-3]和复频率法[4]则被广泛地应用于微/纳谐振器热弹性阻尼的理论分析和求解中.研究者们采用不同的结构变形理论、不同热传导理论和不同的求解方法开展了不同结构形式的微/纳谐振器热弹性阻尼的理论研究.其中,许多工作是关于微板谐振器的热弹性耦合振动及其热弹性阻尼的理论研究[5-26].

在已有关于微板谐振器热弹性阻尼的理论研究中,结构振动方程几乎都是基于经典板理论,即基尔霍夫(Kirchhoff) 薄板理论建立的[5-12,4-23].考虑材料性质为均匀、各向同性的情况,采用复频率法,文献[5-9]分别求得了矩形和圆(环)形微板谐振器热弹性阻尼的解析解,分析了静电载荷和残余应力、边界条件、几何尺寸以及非线性几何变形对热弹性阻尼的影响规律.Li 等[10]、Fang 等[11-12]分别基于一维、二维和三维热传导理论建立了微板单向耦合的热传导数学模型,并采用能量法求得了矩形和圆形薄板谐振器的热弹性阻尼.最近,基于考虑横向剪切变形的明德林(Mindlin)板理论,马航空等[13]首次给出了四边简支中厚度矩形微板热弹性阻尼的解析解.研究发现随着微板厚度的增大,明德林微板的热弹性阻尼最大值单调减小,而基尔霍夫微板的热弹性阻尼峰值却保持不变.

考虑热传导过程中温度与热流运动间的延滞效应(或波动效应),部分作者[14-16]还利用广义热传导理论研究了微板谐振器的热弹性耦合振动响应.例如,文献[14-15]基于Lord-Shulman 热传导理论求得了微板谐振器热弹性阻尼解析解,分析了延滞时间参数和孔隙率变化对热弹性阻尼的影响.Guo 等[16]采用双向延滞广义热传导模型研究了圆板谐振器的热弹性阻尼.Grover[17]采用Kelvin-Voigt 材料模型下的广义热黏弹性理论研究了均匀微圆板谐振器的热弹性阻尼.最新的研究成果还包括:Ghugh 和Partap[18]在微板的热弹性耦合振动模型建立中同时考虑了热松弛和微拉伸 (micro stretch) 效应,研究结果表明微拉伸增强了微板热弹性阻尼的临界值;Wang 和 Li[19]综合考虑微板的尺度效应、热流矢量的延滞效应、非局部效应以及记忆效应,分别采用复频率法和能量法求得了基尔霍夫板理论热弹性阻尼的解析解.通过大量的数值结果详细地分析了微结构的尺度参数、热松弛时间参数以及边界条件对热弹性阻尼的影响规律.

针对材料性质沿单轴方向的阶梯式变化,部分作者研究了由不同均匀材料分层组成的层合微板的热弹性阻尼特性.Sun 等[20]采用复频率法研究了对称铺设的三层微圆板横向轴对称振动下的热弹性阻尼.采用文献[3]中的方法,忽略温度梯度在面内的变化,Zuo 等[21-22]分别研究了双层和三层微板的热弹性阻尼,发现随频率的变化热弹性阻尼会有多个峰值出现.Liu 等[23-24]采用格林函数法分别求解了二维和三维热传导方程,利用能量法获得了双层微圆板和非对称铺设三层矩形微板热弹性阻尼的级数形式解析解.然而,为了消去结构振动方程中的拉−弯耦合项,文献[20-24]都不约而同地采用了物理中面法将几何中面的面内位移用挠度来表示,从而大大简化了在数学上的求解难度.但是,在材料性质变化关于几何中面不对称的情况下,热−弹耦合振动不但会产生热弯矩,还会产生热薄膜力.而他们在简化中完全忽略了热薄膜力对热弹性阻尼的贡献.

考虑材料性质沿着厚度连续变化,最近亦有少数作者研究了功能梯度材料微板的热弹性阻尼[25-27].基于一阶剪切变形理论和应变梯度理论,文献[25]研究了材料性质沿横向幂函数变化的功能梯度微板的热弹性阻尼.应用泰勒级数展开方法,获得了变系数热传导方程的级数形式的解析解,进而获得了四边简支功能梯度中厚度矩形板的热弹性阻尼.基于基尔霍夫经典板理论和准一维热传导理论,Li 等采用解析方法研究了功能梯度圆形微板[26]和矩形微板[27]的热−弹耦合自由振动响应,并成功地将分层均匀化方法应用于复杂变系数热传导方程的求解.获得了材料性质沿厚度按幂函数变化的功能梯度微板的热弹性阻尼解析解.首次定量地分析了热薄膜力对热弹性阻尼的影响程度.

本文拟在前期工作[13,27]的基础上进一步研究横向剪切变形对功能梯度中厚度微板热弹性阻尼的影响规律.基于明德林板理论和单向耦合的准一维热传导理论建立材料性质沿厚度连续变化的功能梯度矩形微板热−弹耦合自由振动数学模型.采用分层均匀化方法求解复杂变系数热传导微分方程,利用结构振动特征值问题在数学上的相似性求得用基尔霍夫均匀微板的无阻尼固有频率表示的四边简支功能梯度材料微板的复频率,进而采用复频率法提取表征热弹性阻尼的逆品质因子.最后,以材料性质按照幂函数变化的两种陶瓷−金属 (Al-SiC 和 Ni-Si3N4) 功能梯度微板为例,给出热弹性阻尼的数值结果,定量分析剪切变形、梯度指数以及几何尺寸对热弹性阻尼的影响规律.

1 问题的控制方程

1.1 运动方程

图1 为微板的几何尺寸和坐标示意图,基于明德林板理论[13,25,28],可给出下列功能梯度中厚度板的位移场

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

其中,t为时间;u0,v0和w0分别为几何中面上一点分别在x,y和z方向的位移分量;φx和 φy为中面法线分别绕y轴和x轴的转角.

根据弹性力学的几何方程可得微板的应变场

对于线弹性材料由胡克定律可得到应力分量

其中,θ (x,y,z,t)=T(x,y,z,t)−T0为变温场,是由热−弹耦合振动产生的;T为瞬态温度,T0为初始温度;E,ν 和 α 分别为弹性模量、泊松比和热膨胀系数,它们都是坐标z的函数.由于功能梯度板的泊松比通常相比其他材料参数变化很小,为了后面方便计算,将其看作常数(取组分材料的平均值).

板的运动方程可用等效内力和内力矩表示为[13,27]

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

ρ为板的质量密度.式(4)～式(8)中的等效内力和内力矩分别定义为

其中ks=5/6 为剪切修正系数.

将式(2)、式(3)代入式(10),可得用位移分量表示的等效内力和弯矩

其中NT和MT分别为热薄膜力和热弯矩,其定义为

方程(11)中的刚度系数由下式给出

将式(11)代入式(4)～式(8)可得用位移表示的运动方程.这是包含5 个基本未知函数的相互耦合的偏微分方程组,其中还包含由温度场来确定的热薄膜力和热弯矩.因此,联立求解这些微分方程在数学上是很困难的.本文通过消去面内位移和转角,最终将5 个方程转化为只用挠度表示的运动方程.首先将式(11a)～ (11c) 分别代入式(4)和式(5),并分别对两式关于坐标x和y求偏导数后相加,利用式(13)可得

类似地将式(11d)～ (11f)代入式(6)和式(7),利用式(8)可得

利用式(11g)和(11 h)可得方程(8)的位移形式

最后,将式(14)代入式(16),忽略面内应变 ε0的惯性项,并利用式(17)可得只用横向位移w表示的运动方程

上式中热弯矩和热薄膜力由热−弹耦合产生的变温场确定.

1.2 热传导方程

根据单向耦合的热弹性动力学理论,可给出功能梯度微板的热传导方程[27]

其中,k,C分别是热传导系数和比热,它们都是坐标z的已知函数;e是体积应变.板的体积应变可表示为

将上式代入式(20),并考虑到板内温度改变主要是由横向振动引起的,因此,为了方便求解可忽略温度场在面内梯度变化,可得单相耦合的准一维热传导方程[13,27]

2 热弹性阻尼的求解

2.1 自由振动响应

假设位移场和温度场的动态响应同步,则系统的热−弹耦合自由振动响应可表示为下列调和形式

这时振动方程(28)中没有热薄膜力,热传导方程(29)中的系数是常数,于是可容易得到其解析通解[13].

2.2 热传导方程的求解

利用文献[26-27,29]中的分层均匀化方法将功能梯度材料微板沿板厚平均划分为N层子区域,当分层数足够大时每层的材料性质可近似看作均匀的.这样,就可以将变系数微分方程(25)转化为分别定义在各分层上的一系列常系数的微分方程组

其中

考虑上下表面为绝热的功能梯度微板,温度边界条件以及层间界面处连续性条件可分别表示为

分析式(31)和齐次边界及连续性条件式(33),系数Aj和Bj可表示为

2.3 热弹性阻尼

热弹性阻尼要从振动方程(24)在给定边界条件下的复频率中提取.而要求解结构自由振动响应,就需要将方程(24)中的热薄膜力和热弯矩用挠度表示.为此,先将式(23)代入式(14),忽略关于中面应变 ε0的惯性力可得

再将式(35)代入式(27)积分可得

将式(23)代入式(17)可得

最后,将式(40b)、式(40c)和式(42)代入方程(24)可得以挠度的振幅表示的功能梯度明德林微板的结构振动方程

微分方程(43)可以写成标准形式

等温的均匀基尔霍夫微板的自由振动方程可以表示为[27]

考虑边界条件为周边简支,则方程(47)的边界条件可表示为

可以证明(见附录A),方程(45)在四边简支约束下的边界条件也可表示为

从而,由微分方程边值问题式(45)、式(49) 与式(47)、式(48)在数学上的相似性,可得它们特征值之间的关系

将式(46)代入式(50),可得关于固有频率代数方程

最后,采用复频率法可得用逆品质因子表示的微板的热弹性阻尼

其中 R e(ω) 和 I m(ω) 分别是复频率的实部和虚部.

3 数值结果与讨论

数值计算时分别选取Ni-Si3N4和Al-SiC 两种由金属和陶瓷组分相复合而成的功能梯度微板.在表1 中列出了上述四种均匀组分材料在室温T0=300 K 下的物性参数.假设材料体积分数在厚度方向按幂函数连续变化,则等效材料性质参数可由混合律给出[26-29]

表1 功能梯度微板组分材料的物性参数值 (T0=300 K)[24,26]Table 1 Values of the parameters of the material properties of the constituents of FGM micro plate (T0=300 K) [24,26]

其中,n≥0 表示材料性质梯度变化指数;Pc和Pm分别表示纯陶瓷和纯金属材料的物性参数.在计算中选择金属板为参考均匀板.将式(54)分别代入式(9)、式(13) 可计算出刚度系数和惯性参数.将式(54)代入式(30)可得各分层内的常系数热传导方程.

关于分层均匀化方法的收敛性已在前期研究中[26-27,29]给以了检验和讨论.在后面的数值计算中取分层数N=500,这已达到非常高的精度.

首先,在n=0 的特殊情况下,表2 中给出了纯陶瓷(Si3N4)明德林正方形微板在不同振动模态和边厚比(a/h)下的热弹性阻尼值.并与相应基尔霍夫微板的结果进行了比较.由此可见,随着板厚和模态阶数的增加,两种板理论预测的热弹性阻尼之间的相对误差变大.在a/h=10 时,最大误差达到了10%.进一步在表3 中给出了具有不同梯度变化指数和不同边厚比的功能梯度(Ni-Si3N4)正方形微板(h=1 μm)在基频振动时的热弹性阻尼值.从两种板理论的预测结果比较可见,随着金属组分的增加,两种理论预测的值相对误差随之增大.上述数值结果定量地表明,随着板厚的增大剪切变形对热弹性阻尼的影响变得显著.

表3 两种板理论下功能梯度正方形微板的热弹性阻尼 (Q−1×104) 比较 (h=1 μm,一阶模态)Table 3 Comparison of TED (Q−1×104) in an FGM (Ni-Si3N4) square plate based on the two plate theories (h=1 μm,1st mode)

Err=×100%

为了更进一步考查剪切变形对热弹性阻尼的影响,在图2 中给出了不同边厚比对应的功能梯度正方形微板(Ni-Si3N4)的热弹性阻尼随着板厚连续变化的特性曲线.图中结果再次表明剪切变形对中厚板和厚板的影响显著.但是,在a=30h时两种板理论的预测值已经非常接近.这说明对于薄板,基尔霍夫板理论已给出十分满意的解答.

图2 具有不同边厚比 a/h 的功能梯度(Ni-Si3N4)微板的热弹性阻尼与板厚 h 的关系曲线(一阶模态)Fig.2 Curves of TED versus the plate thickness h of FGM (Ni-Si3N4) micro plate with different side-to-thickness ratio a/h (1st mode)

图3 中给出了边长为a=100 μm 的功能梯度(Al-SiC)正方形微板的热弹性阻尼随着材料性质梯度变化指数和板厚尺寸的变化规律.在厚度h<6 μm时,热弹性阻尼随着指数的增大(金属组分增加)先增大后减小,有峰值出现;当h>6 μm 后热弹性阻尼则为单调递增的.这种热弹性阻尼随指数增大在一定范围内发生的非单调性是由于材料性质非均匀性在这部分区间的急剧变化而引起的复杂热−弹耦合振动引起的.从图中可见,极大值发生在材料性质梯度变化最剧烈的区间内.这一结果也和所选取的两种组分材料的性质有关.对于Ni-Si3N4组分材料,由文献[27]中的结果可知在板很薄时热弹性阻尼随梯度指数增加会出现极小值.由此可见,材料性质的非均匀性以及组份材料的物性参数不同会导致这种复杂的热−弹耦合振动效应.由此可见非均匀材料谐振器热−弹耦合振动响应的复杂性.

图3 具有不同厚度的功能梯度 (Al-SiC) 正方形微板的热弹性阻 Q −1 尼与梯度指数 n 间的变化曲线(a=b=100 μm)Fig.3 Curves of TED Q −1 versus the gradient index n of an FGM (Al-SiC) micro square plate with different values of the thickness (a=b=100 μm)

图4 中绘出了具有不同材料梯度指数的正方形微板(a=b=100 μm) 的热弹性阻尼随着板厚的连续变化的特性曲线.结果表明,热弹性阻尼随着厚度的增大先单调增加,达到极大值后单调减小.但是,对于两种不同组分材料的微板,热弹性阻尼极大值的变化规律并不相同.Ni-Si3N4微板的热弹性阻尼最大值随梯度指数增大单调增加,而Al-AiC 微板的热弹性阻尼的最大值的变化却并非单调.为了更加清晰地显示这种变化规律,图5 给出两种微板的最大热弹性阻尼与梯度指数n的连续变化曲线.图中可见,同等条件下Al-SiC 微板的−h曲线是局部非单调的.

图4 具有不同梯度指数的功能梯度正方形微板的热弹性阻尼与板厚之间的关系曲线(a=b=100 μm)Fig.4 Curves of TED versus the plate thickness of an FGM square micro plate with different values of the gradient index (a=b=100 μm)

图5 正方形功能梯度微板的最大热弹性阻尼值随材料梯度指数变化的特性曲线Fig.5 Characteristic curves of the maximum TED versus the material gradient index of FGM rectangular micro plate

图6 中给出了对应不同长宽比 (a/b) 的功能梯度 (Ai-SiC) 矩形微板在两种板理论下的热弹性阻尼随板厚的变化曲线.结果表明,长度固定后随着宽度的减小,基尔霍夫微板的热弹性阻尼的峰值保持不变,而明德林微板的峰值有所减小.然而,使得热弹性阻尼取得峰值的临界厚度hcr随着长宽比的增大而减小.

为了反映振动模态对热弹性阻尼的影响规律,图7 给出了梯度指数n=0.5 的功能梯度 (Al-SiC) 正方形微板分别以四前阶模态振动时的热弹性阻尼与厚度关系曲线.其中的变化特性与图6 类似.结果表明,微板抗弯刚度的增大(长度给定,宽度减小;振动模态阶数增大) 对热弹性阻尼的峰值影响甚微,但是对临界厚度的影响显著.

图6 不同长宽比的FGM (Al-SiC)矩形微板的TED 随板厚变化的特性曲线(n=0.5, a=100 μm)Fig.6 Curves of the TED versus the plate thickness of an FGM(Al-SiC) rectangular micro plate with different length-to-wideness ratio (n=0.5, a=100 μm)

图7 对应不同振动模态的FGM (Al-SiC)正方形微板的热弹性阻尼随板厚变化的特性曲线(n=0.5, a=b=100 μm)Fig.7 TED versus of the plate thickness of an FGM (Al-SiC) square micro plate corresponding to different vibration modes(n=0.5, a=b=100 μm)

4 结论

基于明德林一阶剪切变形板理论和准一维单向耦合的热传导理论建立了材料性质沿厚度按连续变化的功能梯度矩形微板热−弹性耦合自由振动控制微分方程.采用分层均匀化方法求得了复杂变系数热传导方程半解析解.在四边简支边界条件,利用结构振动微分方程边值问题的相似性,获得了用相同尺寸和边界条件的参考均匀无阻尼基尔霍夫微板的固有频率表示的功能梯度明德林微板的复频率,从而由复频率法求得了代表热弹性阻尼的逆品质因子半解析解.通过数值算例定量地分析了材料梯度变化指数、组份材料物性参数、几何尺寸以及振动模态对功能梯度微板热弹性阻尼的影响规律.通过与基尔霍夫薄板理论的预测结果进比较,定量地分析了剪切变形对热弹性阻尼的影响程度.结果表明,基尔霍夫板理论预测的热弹性阻尼值大于明德林板理论的预测值.而且随着板厚的增大二者的差别逐渐变得显著.两种理论预测的热弹性阻尼随梯度指数和板厚的变化趋势是相同的.

附录

考虑简支边界条件,在x=0,a处有