张靖华，赵幸幸，李世荣

(1.兰州理工大学理学院，甘肃 兰州730050；(2.扬州大学建筑科学与工程学院，江苏 扬州 225127)

Hamilton体系下功能梯度梁的热冲击动力屈曲分析*

张靖华1，赵幸幸1，李世荣2

(1.兰州理工大学理学院，甘肃 兰州730050；(2.扬州大学建筑科学与工程学院，江苏 扬州 225127)

在Hamilton体系下，基于Euler梁理论研究了功能梯度材料梁受热冲击载荷作用时的动力屈曲问题；将非均匀功能梯度复合材料的物性参数假设为厚度坐标的幂函数形式，采用Laplace变换法和幂级数法解析求得热冲击下功能梯度梁内的动态温度场：首先将功能梯度梁的屈曲问题归结为辛空间中系统的零本征值问题，梁的屈曲载荷与屈曲模态分别对应于Hamilton体系下的辛本征值和本征解问题，由分叉条件求得屈曲模态和屈曲热轴力，根据屈曲热轴力求解临界屈曲升温载荷。给出了热冲击载荷作用下一类非均匀梯度材料梁屈曲特性的辛方法研究过程，讨论了材料的梯度特性、结构几何参数和热冲击载荷参数对临界温度的影响。

功能梯度材料；Euler梁；热冲击；辛方法；动力屈曲

功能梯度材料(functionally graded materials, FGM)是由两种或以上不同材料制成的非均匀复合材料，具有以梯度形式连续变化的材料组份，使得其物理性能也沿梯度方向连续变化，可使构件中的应力集中降到最小[1-2]。由热障FGM所制成的结构常工作在高温梯度环境中，受静态或动态热载荷作用，所引起的热屈曲和热过屈曲吸引了很多力学工作者的研究兴趣[1-11]。

FGM结构的静态或动态热屈曲或热过屈曲问题，由于控制方程为变系数的微分方程而很难解析求解，目前一般都采用数值法近似求解，例如传统的有限元法、差分法、摄动法、打靶法等都用于该问题的求解。S.M.N.Mehrian等[3]基于L-S耦合热弹理论，研究了径向热冲击载荷作用下FGM环板的动态稳定特性。B.Mirzavand等[4]结合物性参数的温度依赖性研究了压电FGM圆柱壳的热冲击动态屈曲问题，用Budiansky稳定判别准则求解获得动态屈曲温度；同时B.Mirzavand等[5]又在文献[4]研究的基础上，基于高阶剪切变形壳理论，求解获得相同结构的动态后屈曲平衡路径和屈曲温度。以上将常规的有限元法或有限差分法用于FGM结构的分析，由于材料性质宏观上的非均匀性，求解时需要划分大量的单元或网格，计算工作量巨大。L.S.Ma[6]基于剪切变形板理论，采用打靶法研究了FGM圆板在机械和热载荷共同作用下的非线性弯曲及过屈曲；S.R.Li等[7-8]采用打靶法分别研究了可伸长FGM Timoshenko梁和FGM缺陷板在热、机载荷同时作用下的屈曲与后屈曲，指出对于边界固支的FGM梁和完善板，即使作用非均匀升温载荷，其变形仍为分叉屈曲，而缺陷FGM板在热载荷作用下无分叉屈曲。但采用打靶法时，对于外形规则结构的静态屈曲容易实现，对于动力屈曲问题的偏微分方程却无法求解。N.L.Shegokar等[9]采用摄动法研究了压电FGM梁在热、电、机多种载荷同时作用下的过屈曲；M.Shariyat[10]考虑预应力和初始缺陷的影响研究了热-力冲击功能梯度圆柱壳的动态屈曲特性，并由改进的Budiansky稳定判别准则求得屈曲临界载荷；K.J.Sohn[11]研究了功能梯度材料板受热和空气动力载荷作用时的静态稳定和颤振特性。然而这些研究都是通过近似求解高阶微分或者是偏微分方程而实现，对于静态稳定性问题尚可，但是对于动力稳定性问题，由于需考虑载荷随时间的变化特性等重要因素，若仍然采用经典弹性力学方法求解，难度将很大。并且求解中还常采用“模态”展开法将偏微分方程退化为常微分方程，该法很难判定是否还存在其他形式屈曲模态。

相比之下，在Hamilton体系[12]下基于辛几何方法[13]，利用分离变量、辛本征展开等方法研究结构的稳定特性，避免了经典弹性力学方法需要求解偏微分方程组的瓶颈，容易得到问题的解答。并且将临界载荷与屈曲模态的求解归结为辛本征值与本征解的求解，无需选择屈曲判别准则。刘淼[14]针对FGM结构，提出了的辛空间有限元-时间子域法的列式方法，该方法中引入了动量作为基本变量,通过对时间域离散，可将动力学问题的求解退化为求解线性方程组，通过分析FGM矩形板的动力响应，表明了辛空间有限元-时间子域法的计算精度高、实用性强。辛方法在结构稳定性方面的应用，褚洪杰等[15]在Hamilton体系下，采用辛本征解展开法，研究了均匀梁的热屈曲和热过屈曲；J.B.Sun[16]采用辛方法研究了热和压缩载荷作用下FGM圆柱壳的静态屈曲。目前，采用辛方法求解FGM结构热冲击稳定性的研究成果极为少见。本文中，建立热冲击下功能梯度Euler梁动态热屈曲的辛方法求解过程，借助于Hamilton体系的完备性解析求解完备的屈曲模态空间，并探索临界载荷、屈曲模态与辛体系中的辛本征值、本征解的相互对应关系，并用数值模拟方法讨论其变化规律和影响因素。

1 基本方程

考虑一长为l，厚为h，宽为b的矩形截面FGM梁。梁的上表面为纯金属，下表面是纯陶瓷，中间由陶瓷和金属依照变化的体积含量连续过渡复合而成。选取变形前的轴线作为x轴，原点在左端面形心处。梁两端固定，无初始变形及速度，在下表面受均布热冲击载荷作用，上表面与外部环境进行热交换，研究其热冲击屈曲特性。

假设FGM组分材料的体积分数沿厚度方向以幂函数形式连续变化，以k表示陶瓷的体积分数指数(0≤k<∞)，其等效物性参数(弹性模量E、热膨胀系数α、密度ρ、热传导系数K以及热容C等)基于Voigt等应变假设，表示为成分体积分数的函数[2,8]。通常材料的泊松比ν变化很小，为了简化计算，假定ν为常数[2,8,17]。

1.1 正则方程

基于Euler梁理论，考虑FGM为线性热弹性，则FGM梁表示能量的Lagrange函数为：

(1)

式中：Ek和Ep分别为系统的动能和势能；w(x,t)为梁轴线上x点的挠度，t为时间；κ为曲率；I为FGM梁单位长度的质量，D为FGM梁的刚度系数；NT和MT分别为有温度T引起的热轴力和热弯矩。I、D、NT、MT定义为：

(2a)

(2b)

(2c)

式中：T(z,t)为温度的变化，简称升温。

(3)

根据Hamilton原理，利用Hamilton体系下的对偶正则方程，由上式可得：

(4)

同时，考虑两端固定的FGM梁，挠度和转角为0的边界条件在Hamilton体系下写为：

(5)

式(4)～(5)分别表示Hamilton体系下辛空间中功能梯度梁热冲击动力稳定问题的对偶正则方程以及边界条件，方程中含有未知的随温度变化的热轴力，我们首先求得非均匀功能梯度梁中的瞬态温度场方可求解该方程组。

1.2 热传导方程及其求解

考虑功能梯度梁的下表面从t=0时刻开始加升温载荷T(-h/2,t)=ΔT·f(t)，ΔT为载荷幅值，f表示热载荷随时间变化的规律。功能梯度梁的上表面与外界进行热交换，换热系数用hr表示，该一维热传导问题的热传导方程及初边值条件由傅立叶热传导定律描述如下：

(6a)

(6b)

具体计算时考虑指数函数形式的升温载荷，即f(t)=1-e-at，其中a为热载荷参数。式(6a)为变系数偏微分方程，直接求解较困难。采用Laplace变换法，令ψ(z,s)=L[T(z,t)]，将式(6a)和(6b)对时间t进行Laplace变换可得：

(7a)

(7b)

(7c)

式(7)为关于ψ(z,s)的变系数常微分方程和初边值条件，其中s为复变量。假定方程(7a)有级数解为如下形式：

(8)

2 正则方程的求解

(9)

式中：X=x/l,Q=q/l,θ=NTL2/D。方程(9)为线性齐次方程，其通解为：

(10)

式中：C1、C2、C3、C4为待定常数。根据量纲一边界条件Q|X=0=Q|X=1=0, (∂Q/∂X)|X=0=(∂Q/∂X)|X=1，可由式(10)得：

(11)

功能梯度梁屈曲的条件是式(11)具有非零解，此时它的系数行列式必为0，将该行列式展开并化简，即可得本问题的分叉条件：

(12)

求解特征方程(12)的根可得无量纲特征值θ1,θ2,…,θn,…，将特征值θn依次代入θ=NTl2/D中，并结合式(2c)求解可获得屈曲载荷幅值ΔT。

对于线性齐次代数方程组(11)，采用归一化方法，令积分常数C4=1，并求解可得其余3个常数:

(13)

然后将以上常数代入式(10)中，可以得到FGM梁的第n阶屈曲模态：

(14)

3 数值结果及讨论

3.1 屈曲模态

采用Newton-Raphson法求解功能梯度材料梁的屈曲分叉条件(式(11))，可得一组量纲一特征值:θ1=39.47，θ2=80.76，θ3=157.91，θ4=238.72，θ5=355.31，θ6=475.60，…。将所得无量纲特征值θn依次代入式(12)可得各阶屈曲模态。FGM梁的前四阶热冲击屈曲模态如图1～4所示。可以看出，不同的热冲击载荷作用下功能梯度梁的失稳构形也不同。

图1 第一阶临界屈曲模态Fig.1 First order critical buckling mode

图2 第二阶屈曲模态Fig.2 Second order critical buckling mode

图3 第三阶屈曲模态Fig.3 Third order critical buckling mode

图4 第四阶屈曲模态Fig.4 Fourth order critical buckling mode

3.2 屈曲升温载荷

为了说明本文理论推导、数值计算的正确性和辛方法求解FGM梁热屈曲的有效性，令体积分数指数k=0将FGM梁退化为均匀陶瓷梁，采用辛方法计算了静态热载荷作用下两端固定Euler梁的临界屈曲载荷，并将数值结果与文献[7]中Euler梁和考虑剪切效应的Timoshenko梁的打靶法计算结果同时列于表1中，其中选取了文献[7]所采用的量纲一热载荷参数λcr=L2αmTm。由表可见，本文中采用辛方法所得结果与文献[7]采用打靶法计算所得结果极其接近，说明辛方法计算结果正确可靠。表格中两种梁理论的结果之间的差别反映了剪切变形对临界热载荷的影响。显然由于考虑剪切变形时梁的柔性增加，临界载荷较低，但从表中也发现剪切变形的影响随着FGM梁的长细比的增加而逐步减弱，所以对于细长梁可将剪切变形忽略而看作Euler梁。

表1 陶瓷梁的静态热屈曲量纲一临界温度

以下求解升温载荷时，选取组分材料分别为陶瓷SiC和金属Ni，两种材料的物理性能参数请见文献[17]。若无特别指出，热冲击载荷参数给定为a=10，作用时间为Δt=5 s，FGM梁的几何尺寸为h=1 cm,l=40 cm；换热系数为hr=50。表2首先列出了给定不同体积分数指数k，即不同体积含量的FGM梁热冲击时的前三阶屈曲升温ΔT。

表2 功能梯度梁的各阶屈曲升温

图5 FGM梁的1、2阶屈曲升温随k的变化Fig.5 Variations of first and second order buckling temperature rise with k

图5进一步绘出了FGM梁受热冲击时第一阶和第二阶屈曲升温ΔT随体积分数指数k的变化关系曲线。由表2和图5可见，随着模态阶数的增高，屈曲载荷值增大，且模态阶数的变化对屈曲载荷值有显著的影响。对于本文中给定的指数函数型的热冲击载荷作用下，FGM梁的屈曲升温介于均质陶瓷和均质金属梁的屈曲升温之间，且随着组份材料的体积分数指数k的增大而减小，即随着k的增大，结构的强度降低，承受热冲击的能力减小。同时可见，当k<2时减小的幅度较大，而当k>2时曲线较平缓，对结构的强度及屈曲升温的影响也变小。

表3列出了给定不同换热系数时体积含量不同的FGM梁受热冲击时的临界屈曲升温(ΔT)l。由表可见，不同换热系数时FGM梁的临界载荷值基本不变，即换热系数对临界屈曲升温基本无影响。

表3 不同换热系数(hr)时FGM梁的临界屈曲升温

表4列出了不同长细比的FGM梁受热冲击时的临界屈曲升温，可以看出，随着长细比λ的逐渐增大，受热冲击时的临界屈曲升温(ΔT)l明显减小，这是由于其弯曲刚度减小的缘故。

表4 不同长细比(λ)下FGM梁的临界屈曲升温(ΔT)l

图6 不同a时FGM梁的临界屈曲升温(ΔT)lFig.6 Variations of critical temperature for specified values of a

图6为给定不同载荷参数a时FGM梁热冲击屈曲的临界升温(ΔT)l随体积分数指数k的变化关系曲线。由图可见，临界屈曲升温随着a的逐渐增大而缓慢减小。当参数a<5时，a对临界温度有影响；但参数a>5时由于临界屈曲升温基本不随a的增大而变化，所以基本无影响；不同参数a下，FGM梁受热冲击时的临界屈曲升温随组份参数k的变化趋势都相同。

表5列出了载荷作用时间Δt不同时体积含量不同的FGM梁受热冲击时的屈曲临界屈曲升温(ΔT)l。由表可见，功能梯度梁的临界屈曲升温随载荷作用时间的增加而逐渐减小。当Δt<5 s时，随作用时间的延长有较大变化，当Δt>5 s时变化很小，最终随时间的增大而趋于一定值。其原因在于，热冲击载荷作用时间越长，梁内的温度分布越均匀，对屈曲临界屈曲升温的影响也逐渐消失。

表5 热冲击载荷作用时间(Δt)不同时FGM梁的临界屈曲升温(ΔT)l

4 结 论

通过对热冲击载荷作用下FGM梁的动力屈曲特性的研究，得出以下结论：

(1)利用辛方法在Hamilton体系下求解FGM结构的热冲击屈曲是可行的，绕开了经典弹性力学解偏微分方程的瓶颈，容易得到问题的解答；

(2)热冲击下FGM梁随着体积分数指数的增大，强度降低，其临界屈曲升温介于均质陶瓷梁和均质金属梁的相应结果之间，随着梯度材料的体积分数指数的增大而逐渐减小；

(3)热冲击载荷作用时间对FGM梁的临界屈曲升温有较大影响，但介质换热系数和载荷变化参数 对FGM梁的临界屈曲升温影响不大。

[1] 仲政,吴林志,陈伟球.功能梯度材料与结构的若干力学问题研究进展[J].力学进展,2010,40(5):528-541. Zhong Zheng, Wu Linzhi, Chen Weiqiu. Pregress in the study on mechanics problems of functionally graded materials and structures [J]. Advances in Mechanics, 2010,40(5):528-541.

[2] 张靖华,潘双超,李世荣.热冲击下功能梯度圆板的动力屈曲[J].应用力学学报,2015,32(6):901-907. Zhang Jinghua, Pan Shuangchao, Li Shirong. Dynamic buckling of functionally graded circular plate under thermal shock [J]. Chinese Journal of Applied Mechanics, 2015,32(6):901-907.

[3] Mehrian S M N, Naei M H. Two dimensional analysis of functionally graded partial annular disk under radial thermal shock using hybrid Fourier-Laplace transform [J]. Applied Mechanics and Materials, 2013,436:92-99.

[4] Mirzavand B, Eslami M R, Shakeri M. Dynamic thermal postbuckling analysis of piezoelectric functionally graded cylindrical shells [J]. Journal of Thermal Stresses, 2010,33(7):646-660.

[5] Mirzavand B, Eslami M R, Reddy J N. Dynamic thermal postbuckling analysis of shear deformable piezoelectric FGM cylindrical shells [J]. Journal of Thermal Stresses, 2013,36(3):189-206.

[6] Ma L S, Wang T J. Relationships between the solutions of axisymmetric bending and buckling of functionally graded circular plates based on the third-order plate theory and the classical solutions for isotropic circular plates [J]. International Journal of Solids and Structures, 2004,41(1):85-101.

[7] Li S R, Zhang J H, Zhao Y G. Thermal post-buckling of functionally graded material Timoshenko beams[J]. Applied Mathematics and Mechanics, 2006,27(6):803-811.

[8] Li S R, Zhang J H, Zhao Y G. Nonlinear thermo-mechanical post-buckling of circular FGM plate with geometric imperfection [J]. Thin Walled Structures, 2007,45(5):528-536.

[9] Shegokar N L, Lal A.Thermo-electromechanically induced stochastic post buckling response of piezoelectric functionally graded beam [J]. International Journal of Mechanics and Materials in Design, 2014,10(3):329-349.

[10] Shariyat M. Dynamic thermal buckling of suddenly heated temperature-dependent FGM cylindrical shells under combined axial compression and external pressure [J]. International Journal of Solids and Structures, 2008,45(9):2598-2612.

[11] Sohn K J, Kim J H. Structural stability of functionally graded panels subjected to aero-thermal loads [J]. Composite Structure, 2007,82(3):317-325.

[12] 徐新生,段政,马源,等.辛方法和弹性圆柱壳在内外压和轴向冲击下的动态屈曲[J].爆炸与冲击,2007,27(6):509-514. Xu Xinsheng, Duan Zheng, Ma Yuan, et al. A symplectic method and dynamic buckling of elastic cylindrical shells under both axial impact and internal or external pressure[J]. Explosion and Shock Waves, 2007,27(6):509-514.

[13] 谈梅兰,吴光,王鑫伟.矩形薄板面内非线性分布载荷下的辛弹性力学解[J].工程力学,2008,25(10):50-53. Tan Meilan, Wu Guang, Wang Xinwei. Symplectic elasticity solutions for thin rectangular plates subjected to nonlinear distributed in plane loadings[J]. Engineering Mechanics, 2008,25(10):50-53.

[14] 刘淼.功能梯度材料结构的非传统Hamilton变分原理及其有限元法[D].上海:同济大学,2008.

[15] 褚洪杰,徐新生,林志华,等.弹性梁非线性热屈曲行为与辛本征解展开方法[J].大连理工大学学报,2011,51(1):1-6. Chu Hongjie, Xu Xinsheng, Ling Zhihua, et al. Nonlinear thermal buckling of elastic beams and expanding method of symplectic eigensolutions[J]. Journal of Dalian University of Technology, 2011,51(1):1-6.

[16] Sun J B, Xu X S, Lim C W. Buckling of functionally graded cylindrical shells under combined thermal and compressive loads [J]. Journal of Thermal Stresses, 2014,37(3):340-362.

[17] Zhang J H, Li G Z, Li S R. DQM based thermal stresses analysis of a FG cylindrical shell under thermal shock[J]. Journal of Thermal Stresses, 2015,38(9),959-982.

(责任编辑 王小飞)

Dynamic buckling analysis of functionally graded beam under thermal shock in Hamilton system

Zhang Jinghua1, Zhao Xingxing1, Li Shirong2

(1.SchoolofSciences,LanzhouUniversityofTechnology,Lanzhou730050,Gansu,China;2.SchoolofCivilScienceandEngineering,YangzhouUniversity,Yangzhou225127,Jiangsu,China)

Based on the Euler beam theory, the dynamic buckling of the functionally graded beam subjected to thermal shock was investigated in the Hamilton system. The material properties of the functionally graded beam were assumed to be graded in the thickness direction according to a simple power law distribution in terms of the volume fractions of the constituents. The transient temperature fields were solved analytically using the Laplace transform and power series method. It was shown that the dynamic buckling problem can be reduced to a zero-eigenvalue problem in the symplectic space, the buckling loading and the buckling mode of the FGM beam correspond to the generalized eigenvalue and eigen solution. The buckling mode and the buckling thermal axial forces can be obtained through bifurcation condition, and the buckling temperature rise of the FGM beam can be obtained by inverse solution. In this research, the solution process for dynamic buckling of the FGM beam subjected to thermal shock using the symplectic method were given, and the effects of the material constitution, geometric parameters and the parameters of thermal shock load on the critical temperature were discussed.

functionally graded materials; Euler beam; thermal shock; symplectic method; dynamic buckling

10.11883/1001-1455(2017)03-0431-08

2015-11-23；

2016-06-20

国家自然科学基金项目(11662008,11262010)

张靖华(1979- ),女,博士，副教授，zjhhrb@163.com。

O347.2 国标学科代码： 13015

A