薄板结构在使用过程中不可避免地要受到热载荷和力载荷的影响,因而有必要研究薄板在各种载荷下的屈曲行为.目前,Nguyen等研究了板壳结构的热屈曲问题[1-7];吴晓等研究了双弹性模量圆板的热屈曲问题[8];彭凡等研究了功能梯度薄板的热屈曲问题[9-12];Mathew等用有限元方法计算了薄板的热屈曲问题[13];Sun和Hsu指出板的厚度对热屈曲行为有重要的影响[14];Chang基于高阶剪切变形理论,采用有限元方法,研究了薄板的屈曲和热屈曲问题[15].但以上研究均未考虑热载荷和力载荷同时作用时,它们相互之间对结构失稳的影响,以致有些分析不准确.
薄板在航空航天、船舶等工程领域应用十分广泛,研究价值大.本文在前人研究的基础上,先推导了热力耦合屈曲临界载荷的表达式,然后验证这些公式的正确性,最后通过单一薄板和复杂薄板(如蜂窝)的算例分析了热载荷与力载荷对结构失稳的影响.
将温度载荷转化为与其相当的节点力,采用平面应力单元对薄板进行计算,应力-应变关系的表达式可写为
σ=D(ε-εT),
(1)
其中
(2)
E为板的弹性模量,μ为Poisson比.单元节点力的关系式可写为
(3)
B为应变-节点位移矩阵,h为板厚.将应力-应变关系式代入节点力关系式,可得
(4)
上式等号左端的第一项反映了力载荷作用下的节点力向量,第二项反映了温度载荷作用下的节点力向量.
基于有限元法,热屈曲前结构的状态方程为
Kd=P,
(5)
其中K为热屈曲前薄板结构的总体刚度矩阵,d为总体位移向量,总体节点力向量的关系式可写为P=PN+PT,其中PN和PT分别表示力载荷产生的节点力向量和温度载荷产生的节点力向量.由上式可求得屈曲前节点位移d,然后由节点位移可求得节点力.
板屈曲的能量变分形式的控制方程为
δΠ=δ(U+V)=0,
(6)
其中U为结构体的弹性势能,V为载荷产生的势能.整体结构的弹性势能可表示为
其中m为单元总数.将应力-应变关系代入上式,可得
(7)
对上式进行简化后,可得到弹性势能的表达式为
(8)
其中Ce为常数,Ke为单元刚度矩阵,Se为热载荷引起的热应力.
假定在一个初始应力状态下发生一个无限小的屈曲位移,薄板内力就基本保持不变.假定薄板内力基本保持不变,在热载荷作用下,薄板内力的表达式可写为其中λT为载荷参数,NxT,NyT,NxyT为由温度载荷引起的热屈曲前的板内力.假定薄板内力基本保持不变,在力载荷作用下,薄板内力的表达式可写为其中λN为载荷参数,NxN,NyN,NxyN为力载荷引起的屈曲前的薄板内力.在温度载荷或者力载荷的作用下,板的内力做的功为
(9)
将形状函数代入上式,并且进行变分运算,可得到
(10)
上式可写为
其中是由力载荷或者热载荷引起的几何刚度矩阵,它可写为
(11)
其中
将单元节点位移换为整体节点位移,并将单元刚度矩阵组装为整体刚度矩阵,就可得到总势能的变分方程为
δΠ=δ(U+V)=δdT(K-λKg)d=0.
(12)
由δdT的任意性可得(K-λKg)d=0.根据所求得的内力状态,就可写出结构屈曲的总几何刚度矩阵Kg,结构的几何刚度矩阵仅与结构的形状函数和屈曲前的薄板内力状态有关.由式(12)可得屈曲控制方程为
(K-λTKg,T-λNKg,N)d=0.
(13)
根据上面的表达式可求得式(13)的特征值,其最小的特征值λT,λN就是待求的临界载荷因子.临界载荷的表达式为
Pcr=λTP0,T+λNP0,N.
(14)
上面是热力耦合载荷下屈曲临界载荷的求解过程.在式(13)中,由于出现两个未知量λT,λN,所以求解时须给定其中一个未知量,来求解另一个未知量.关于式(13)的求解可以分为以下两种情况:
① 给定一个λN值求解λT;
② 给定一个λT值求解λN.
相应的控制方程可写为
(K′-λK′g)d=0,
(15)
其中
由上述的有限元公式,在MATLAB环境下编写有限元程序,对薄板的热力耦合屈曲进行求解,有限元程序的编写流程如图1所示.
图1 有限元编写流程
Fig. 1 The FEM programming process
取各向同性均匀材料的四边简支约束的矩形板,假设它的几何尺寸为a=50 mm,b=30 mm,h=1 mm,材料弹性模量取为E=40 GPa,Poisson比ν=0.3,热膨胀系数α=2E-6,以此来做屈曲分析.对矩形板进行载荷施加,分两种情况:第一种是仅受力载荷的作用,即释放矩形板上边界y方向的自由度,然后施加力载荷p=1 MPa;第二种是仅受温度载荷的作用,即在板四边简支约束的情况下,然后施加均匀分布的温度场.在这两种条件的作用下,得到表1和表2的结果.由表1和表2的结果可以看出,求解的精度与有限元网格的划分有关,网格划分越密,精度就越高,但在超过30×30时,精度提高不显著,所以在接下来计算薄板时采用30×30的网格进行计算.
表1 仅力载荷下的薄板一阶特征值(单位: MPa)
Table 1 First-order eigenvalues of the thin plate under mechanical load only(unit: MPa)
ANSYS result mesh5×510×1015×1520×2025×2530×30246.09result262250.7249248.5248.27248.15error6.1%1.6%1.3%1.1%0.9%0.8%
表2 仅均匀温度载荷下的薄板一阶特征值(单位: ℃)
Table 2 First-order eigenvalues of the thin plate under uniform temperature load(unit: ℃)
ANSYS result mesh5×510×1015×1520×2025×2530×30513.01result605.92553.33548.70547.30546.69546.36error15.2%7.2%6.3%6.2%6.1%6.1%
将表1和表2的数据与ANSYS的解进行比较,可见本研究的解与ANSYS的解符合良好,证明了本文计算程序的可靠性.
首先,以简单模型结构来做屈曲分析,以单一薄板为例.假设单一薄板结构的材料为金属材料(E=40 GPa,ν=0.3,α=2E-6),薄板尺寸为a=50 mm,b=30 mm,h=1 mm,以此为例来分析力载荷分量与热载荷分量对结构失稳的影响.在单一薄板上施加沿高度方向线性分布的热载荷,热载荷的施加需要给定温度场的分布形式,给定温度沿薄板y方向线性增长,由0增加到临界值Tcr,如图2所示.
图2 热力耦合载荷
Fig. 2 The thermo-mechanical coupling load
下面是力载荷与热载荷共同作用下的计算结果.在仅受到热载荷的作用时,单一薄板的一阶临界载荷为802 ℃;在仅受到力载荷作用时,单一薄板的一阶临界载荷为247.5 MPa.在力热载荷共同作用时,分两种情况:第一种是在给定力载荷由0~247.5 MPa的变化下, 求临界热载荷, 如图3所示; 第二种是在给定温度载荷由0~802 ℃的变化下, 求临界力载荷, 如图4所示.
在给定热载荷或力载荷的作用下,随着单一薄板长宽比和宽厚比的变化,图3和4给出了临界载荷的变化曲线.随着给定热载荷或力载荷的增加,临界载荷降低且几乎呈线性变化.在热载荷与力载荷共同作用下,热载荷与力载荷的增加,都会使得结构的抗失稳能力降低.同时,由图3和4可看出,随着长宽比的增大和宽厚比的减小,曲线逐渐向内侧移动,这说明了增大长宽比和减小宽厚比,会降低结构的抗失稳能力,其中宽厚比对结构抗失稳能力的影响要大于长宽比.
图3 给定力载荷下的临界热载荷
Fig. 3 Critical thermal loads under given mechanical loads
图4 给定热载荷下的临界力载荷
Fig. 4 Critical mechanical loads under given thermal loads
(a) p=0 MPa, T=1℃
(b) p=100 MPa, T=1℃(c) p=250 MPa, T=1℃
图5 热力耦合下的一阶屈曲模态图
Fig. 5 First-order buckling modal shapes of the thin plate under thermo-mechanical coupling loads
注 为了解释图中的颜色,读者可以参考本文的电子网页版本.
下面将从求解得到的特征向量中提取每个节点的扰度, 并将其作为颜色插值, 绘制单一薄板一阶屈曲的三个模态图(如图5所示), 来做不同载荷条件下的比较.这三个模态图是二维图像.
由于温度载荷为线性分布,靠近薄板上边界温度高,下边界温度低,所以波形较为靠近上边界.由图5可看出随着力载荷的增加,屈曲波形逐渐向下边界移动,表明力载荷分量对单一薄板失稳的影响越来越显著.
然后,以复杂的模型结构来做屈曲分析,以蜂窝为例.蜂窝可看作由六块离散的薄板组合成的周期性、重复性的结构,是一种复杂的薄板结构.蜂窝结构常常作为热防护以及隔热材料应用在航天器保护罩、飞机机头保护罩、飞机机翼壁板中,其使用环境常常处在热力耦合场中.假设蜂窝的单个薄板尺寸、 材料属性、 载荷施加方式与算例1的图2相同, 蜂窝结构如图6所示.
图6 蜂窝结构
Fig. 6 The honeycomb structure
图7 热力耦合载荷作用下给定力载荷的临界温度载荷
Fig. 7 Critical temperature loads under given mechanical loads
下面以六边形蜂窝结构为例,做一个加载模型分析.在仅受到力载荷作用时,一阶临界载荷为112.56 MPa;在仅受到热载荷作用时,一阶临界载荷为933.22 ℃.在力热载荷共同作用时,分两种情况: 第一种是在给定力载荷由0~112.56 MPa的变化下,求临界热载荷,如图7所示;第二种是在给定温度载荷由0~933.22 ℃的变化下,求临界力载荷,如图8所示.
图8 热力耦合载荷作用下给定温度载荷的临界力载荷
Fig. 8 Critical force loads under given thermo loads
算例1和算例2中,力载荷分量与热载荷分量近似呈线性相关,屈曲临界载荷值的表达式可以是T/Tcr+σ/σcr=1,其中T/Tcr体现温度载荷的作用,σ/σcr体现力载荷的作用,两者对临界载荷的贡献近似呈线性关系.
薄板结构在航空航天和船舶领域中应用广泛,并且在航空领域薄板一般作为主承力结构使用,在使用过程中不可避免地会受到热载荷与力载荷的影响,所以在工程结构设计中须考虑这些因素对结构抗失稳能力的影响.本文基于有限元方法,采用Rayleigh-Ritz法,研究了单一薄板和复杂薄板的热力耦合屈曲问题.在假设力载荷与热载荷同时加载的条件下,分析了力载荷分量与热载荷分量对失稳的影响.研究结果表明:屈曲临界载荷值的表达式是T/Tcr+σ/σcr=1,其中T/Tcr体现温度载荷的作用,σ/σcr体现力载荷的作用,两者对薄板临界载荷的贡献近似呈线性关系.当结构处于热力耦合环境中,热载荷分量与力载荷分量会降低其抗失稳能力,可按表达式T/Tcr+σ/σcr=1来改进结构的设计.上述研究可为飞机、船舶、航天器等提供结构设计的依据,为提高它们的稳定性提供技术支持.
[1] NGUYEN N D, NGUYEN T K, NGUYEN T N, et al. New Ritz-solution shape functions for analysis of thermo-mechanical buckling and vibration of laminated composite beams[J]. Composite Structures, 2018, 184(15): 452-460.
[2] SABZIKAR M, BOROUJERD Y. Thermal buckling of piezo-FGM shallow spherical shells[J]. Mechanica, 2013, 48(4): 887-899.
[3] 洪杰. 火焰筒结构局部热屈曲分析方法[J]. 北京航空航天大学学报, 2010, 36(2): 248-252.(HONG Jie. Local thermal buckling analysis method of combustor liner[J]. Journal of Beijing University of Aeronautics and Astronautics, 2010, 36(2): 248-252.(in Chinese))
[4] 袁武, 王曦, 宋宏伟, 等. 轻质金属点阵夹层板热屈曲临界温度分析[J]. 固体力学学报, 2014, 35(1): 1-7.(YUAN Wu, WANG Xi, SONG Hongwei, et al. Thermal buckling and its critical temperature analysis of sandwich panels with metal-truss core[J]. Chinese Journal of Solid Mechanics, 2014, 35(1): 1-7.(in Chinese))
[5] 李忱, 田雪坤, 王海任, 等. 薄球壳在均布外压与温度耦合作用下的热屈曲研究[J]. 应用数学和力学, 2015, 36(9): 924-935.(LI Chen, TIAN Xuekun, WAGN Hairen, et al. Study on thermal buckling of thin spherical shell under the coupling of reuniform external pressure and temperature[J]. Applied Mathematics and Mechanics, 2015, 36(9): 924-935.(in Chinese))
[6] 夏巍, 赵东伟, 冯宇鹏. 基于Mindlin横剪变形理论的功能梯度板热屈曲分析[J]. 应用力学学报, 2016, 33(1): 13-18.(XIA Wei, ZHAO Dongwei, FENG Yupeng. Thermal buckling analysis of functionally graded plates based on Mindlin’s transverse shear deformation theory[J]. Chines Journal of Applied Mechanics, 2016, 33(1): 13-18.(in Chinese))
[7] 朱永安, 王璠, 刘人怀. 考虑横向剪切的对称圆柱正交异性层合扁球壳的热屈曲[J]. 应用数学和力学, 2008, 29(3): 263-271.(ZHU Yongan, WANG Fan, LIU Renhuai. Thermal buckling of axisymmetrically laminated cylindrically orthotropic shallow spherical shells including transverse shear[J]. Applied Mathematics and Mechanics, 2008, 29(3): 263-271.(in Chinese))
[8] 吴晓, 赵均海, 黄志刚. 双模量材料圆板热弯曲及热屈曲的研究[J]. 应用力学学报, 2015, 32(4): 549-555.(WU Xiao, ZHAO Junhai, HUANG Zhigang. Study on thermal bending and thermal buckling of circular plates with double modulus materials[J]. Chinese Journal of Applied Mechanics, 2015, 32(4): 549-555.(in Chinese))
[9] 彭凡, 顾勇军. 热环境中功能梯度圆柱薄壳分岔屈曲的边界约束效应[J]. 固体力学学报, 2011, 32(5): 475-482.(PENG Fan, GU Yongjun. Effect of boundary constraints on bifurcation buckling of functionally graded material circular cylindrical shells in thermal environment[J]. Acta Mechanica Solida Sinica, 2011, 32(5): 475-482.(in Chinese))
[10] KOCATURK T, AKBAS S D. Post-buckling analysis of Timoshenko beams made of functionally graded material under thermal loading[J]. Structural Engineering and Mechanics, 2012, 41(6): 775-789.
[11] LEVYAKOV S V. Elastica solution for thermal bending of a functionally graded beam[J]. Acta Mechanica, 2013, 224(8): 1731-1740.
[12] JABERZADEH E, AZHARI M, BOROOMAND B. Thermal buckling of functionally graded skew and trapezoidal plates with different boundary conditions using the element-free Galerkin method[J]. European Journal of Mechanics A: Solids, 2013, 42: 18-26.
[13] MATHEW T C, SINGH G, RAO G V. Thermal buckling of cross-ply composite laminates[J]. Computers & Structures, 1992, 42(17): 281-287.
[14] SUN L X, HSU T R. Thermal buckling of laminated composite plates with transverse shear deformation[J]. Computers & Structures, 1990, 36(5): 883-889.
[15] CHANG J S. FEM analysis of buckling and thermal buckling of antisymmetric angle-ply laminates according to transverse shear and normal deformable high order displacement theory[J]. Computers & Structures, 1990, 37(6): 925-946.
李若愚(1994—),男,助理工程师(通讯作者. E-mail: 412364086@qq.com);
王天宏(1968—),男,教授,博士.