宋彦琦,石博康,李向上

(中国矿业大学(北京)力学与建筑工程学院,北京 100083)

随着科技的不断发展,工程环境愈加复杂严苛,人们对材料的性能要求也随之提升.由于传统的材料已经不能满足科研和工程要求,研究者将目光转向了复合材料.材料属性的不连续性导致传统复合材料有明显的缺陷,易于失效.而功能梯度复合材料(functionally graded material,FGM)在微观结构上连续变化,能够消除不同组分材料的界面效应,具有高强度、高韧性,可减少应力集中等力学性能[1],被广泛应用于机械、化工、电子、核工程和航空航天等领域.在某些严苛的工程环境中,功能梯度板(functionally graded plate,FGP)常会产生非线性变形,因此研究FGP板的非线性变形特征是有意义的.

国内外学者对功能梯度复合材料板的非线性弯曲进行了大量研究.贾金正等[2]采用物理中面的概念,研究了横纵荷载作用下梯度指数和长高比对功能梯度梁变形情况的影响.王雪等[3]采用逐步加载技术,研究了热-机械荷载联合作用下功能梯度梁的荷载-挠度曲线.Khabbaz等[4]基于高阶剪切变形理论研究了功能梯度板的非线性弯曲问题.上述研究均引入了Green应变张量,Green应变张量虽然不会在大变形时产生虚假应变,但没有与之对应的转动张量,故该理论是有缺陷的[5].陈至达[5]提出了S-R和分解定理,克服了Green应变张量的缺点,能够确定唯一的应变.

另外,部分学者基于有限元法对功能梯度板的非线性变形进行了研究.Hassan等[6]研究了SiC-Al功能梯度板的非线性弯曲问题,并计算了塑性区的大小.Nam等[7]基于高阶剪切变形理论,利用多边形有限元法求解了Al-Al2O3板的非线性弯曲问题.利用有限单元法求解大变形问题时,网格会产生严重的变形,严重影响了计算结果的精度.而无网格法可以部分或者彻底地消除网格,保证结果的精度[8-9].1994年,Belytschko等[10]提出无网格伽辽金(element-free Galerkin,EFG)法后,无网格法得到了迅速的发展,并被应用到众多领域.叶志强等[11]和刘世宇等[12]利用基于全方向导数的无网格法,对二维翼型绕流问题进行了求解.Li等[13]利用EFG法对金属成型过程中的接触冲击大变形问题进行了分析.张弛等[14]利用无网格法对变厚度功能梯度材料圆板的自由振动问题进行了分析.陈海昌等[15]基于一阶剪切变形理论,应用无网格法对不锈钢-陶瓷功能梯度板弯曲问题进行了计算.陈祖芳等[16]、罗丹[17]和宋彦琦等[18]推导了2D-S-R-EFG离散方程,求解了二维梁结构的非线性变形问题,宋彦琦等[19]又提出了三维情况下的S-R无网格法.

综上所述,目前研究普遍应用的Green应变张量没有相应的转动张量,而有限元法在求解大变形问题时严重影响结果的精度.另外,S-R和分解定理在功能梯度板非线性变形的问题中的应用较少.如果将S-R和分解定理与无网格法结合起来,就可以建立更加完善的大变形问题求解方法.因此,本工作基于S-R和分解定理,推导出3D-S-R-EFG离散方程,并利用MATLAB编写了S-R无网格法程序,对功能梯度材料板的非线性变形弯曲进行了研究.

1 功能梯度材料模型

功能梯度板由陶瓷和金属混合制成,长度为a,宽度为b,厚度为h,板的上表面受均布荷载q,示意图如图1所示.板的材料参数沿板的厚度方向变化,符合幂律型分布,

图1 功能梯度板示意图Fig.1 Schematic diagram of functionally graded plate

式中:P为板内任一点的材料参数(弹性模量E、泊松比μ、剪切模量G等);Pc和Pm分别是陶瓷和金属的材料参数;Vc为陶瓷的体积分数;n为陶瓷的体积分数指数.

图2给出了体积分数指数n变化时,陶瓷的体积分数随厚度变化的情况.

图2 体积分数随板的厚度变化情况Fig.2 Variation of the volume fraction through the thickness of the plate

2 S-R和分解定理

采用自然拖带系法[5]描述物体的运动,选取两个参考系:固定参考系{Xi}(i=1,2,3)和嵌含在变形体中的自然拖带坐标系{xi}(i=1,2,3).自然拖带坐标系如图3所示.质点变形前后的位置矢量分别为r和R,位移矢量为u,

图3 自然拖带坐标系Fig.3 The co-moving coordinate

定义变形前后的基标矢量:

对式(2)求导,可得

式中:

其中为一阶逆变分量uj对xi的协变导数为初始拖带系中的第二类Christoffel符号.

由上分析可以得到变形前后基矢的转换关系为

式中:F是局部坐标系基矢的变换函数;为Kronecker符号.

S-R和分解定理如下:给定一个物理可能的位移函数,此函数在形变体点集域内是单值连续的,处处有一阶导数,则此运动变换可以唯一分解为正交与对称两个子变换的直和.正交变换体现点集的转动,而对称变换体现点集的形变[5],即

式中:S表示应变张量;R表示转动张量.

3 基于更新拖带坐标系法的增量变分方程

3.1 增量变分方程

利用虚功率原理并结合增量方法,可得

式中:t+ΔtV j‖i表示在t+Δt时刻的速度t+ΔtV i关于t+Δt时刻拖带系协变基矢t+Δtgi的协变导数;t+ΔtR为外力虚功,

考虑准静力学问题,根据速度梯度和应力的线性近似有

将式(12)、(13)代入式(10),并结合Δuj|k=ΔttV j‖k,可得到

3.2 更新拖带坐标系法

图4所示为更新拖带坐标系法.由图4可知:在初始时刻t0,质点位于点P,其初始拖带坐标系的基矢为在t时刻,质点移动到点P′处,其初始拖带系基矢变为定义为t时刻的初始拖带系的基矢,即在t～t+Δt时刻的增量步上,对初始拖带系进行了一次更新,则t时刻瞬时拖带系tgi的tV i‖j可表示为

图4 更新拖带坐标系法Fig.4 The updated co-moving coordinate method

将式(15)代入式(14),可得

式中:

引入速率型物性方程

式中:

利用式(17),可得

忽略体力矩的影响,并令初始系为直线直角坐标系,可得[20]

则式(23)可简化为

式(25)即可作为全局弱势的无网格法的增量变分方程.

4 基于全局弱势的无网格Galerkin法

根据无网格Galerkin法,利用移动最小二乘法(moving least squares,MLS)对t时刻初始拖带系下节点的速度tV进行插值,则有

式中:tφk为t时刻节点k处的MLS形函数;tVk是t时刻k节点的速度参数.采用罚函数法施加边界条件,对式(25)进行修正,可得

对于tV i‖j,有

将式(26)代入式(27),可得

故可将式(17)写为

式中:

对于[Bs]和[Br],有

注意到矩阵[Bs]和[Br]中含有,令初始拖带系和直线直角坐标系同胚,则有这大大降低了计算的复杂度,也体现了更新拖带坐标系法的优势.

将式(30)代入式(27),并利用Δt[V]=[Δu],可得

式中:

5 数值算例

5.1 各向同性板的非线性弯曲

表1为本工作中无网格法数值计算的参数设定.为了检验本工作的方法和程序的正确性,对受均布荷载的四边简支方形薄板进行计算.令FGP两种材料的弹性模量均为5.38×1010Pa,泊松比均为0.3.根据式(1),该板为各向同性板.板的长度和宽度均为25.4 cm,厚度为2.54 cm,节点分布为13×13×3.该算例被广泛应用于验证薄板非线性弯曲计算方法的准确性.

表1 无网格法参数设定Table 1 Essential parameters for S-R-EFG

图5给出了板中面中点的无量纲挠度与荷载的关系.从图中可以看出:在初始的小变形阶段,3种方法的计算结果吻合度高,不同理论对结果的影响较小,这在一定程度上说明了本方法的准确性;在非线性(大变形)阶段,本方法的计算结果与Zhu等[21]和Zhao等[22]的结果相比较小.因为Zhu等[21]采用了一阶剪切变形理论和基于滑动Kriging插值的无网格局部Petrov-Galerkin法,而Zhao等[22]采用了基于一阶剪切变形理论的无网格里兹法,但一阶剪切变形理论会引入一些假设,如假定沿厚度方向的正应变为0,降低了结果的准确性.此外,文献[21-22]均引入了von K´arm´an应变来体现非线性变形,而von K´arm´an应变从本质上来讲,仍然属于Green应变.Green应变张量没有推导出相应的转动张量,是存在理论缺陷的.本工作中的S-R无网格法推导了转动速率和速度的关系,并将其引入到了系统方程中,考虑了转动对变形的影响,故能得出更准确的结果.根据以上讨论,该算例能够说明本工作中S-R无网格法的准确性,并且也可以解释与其他结果的差异.

图5 各向同性板的无量纲挠度和荷载的关系Fig.5 Non-dimensional central deflection versus load parameter of isotropic plate

5.2 功能梯度板的非线性弯曲

表2为算例中FGP材料的力学参数[21].

表2 功能梯度板的力学参数Table 2 Mechanical parameters of FGP

图6给出了FGP在四边简支约束下无量纲挠度和无量纲荷载之间的关系.与文献[21-23]中的结果进行对比,可以看出,在采用13×13×3节点进行计算时,本工作的方法能够得到收敛的结果.与文献[21-22]相同,文献[23]中也引入了von K´arm´an应变来体现非线性,故本工作中的结算结果与其他结果的差异也是可以解释的.

图6 功能梯度板的无量纲挠度和荷载的关系Fig.6 Non-dimensional central deflection versus load parameter of FGP

下面研究体积分数指数对FGP非线性弯曲的影响.图7所示为FGP的挠度和体积分数指数n的关系.由图7可以看出,随着n的增加,FGP的刚度下降,板的无量纲挠度增加,这与Zhu等[21]、Zhao等[22]和Do等[23]得出的规律是相同的.图8所示为在不同体积分数指数下板中面中点的无量纲挠度和宽厚比的关系.从图中可以看出:当FGP的宽厚比从10增加到15时,板的中面挠度迅速减小;但是当板的宽厚比增加到20时,板的中面挠度只产生了很小的变化.因此可以预测,随着板的宽厚比增加,板的中面无量纲挠度会趋近于一个定值.

图7 不同体积分数指数下功能梯度板的无量纲挠度和荷载的关系Fig.7 Non-dimensional central deflection versus load parameter of FGP for the various volume fraction exponents

图8 不同宽厚比情况下功能梯度板的无量纲挠度和荷载的关系Fig.8 Non-dimensional central deflection versus load parameter of FGP for different lengthto-thickness ratios

6 结论

本工作以S-R和分解定理为基础,基于全局弱势的无网格伽辽金法得到了三维情况下的S-R-EFG离散方程,利用编制的3D-S-R-EFG程序求解了Ti-6Al-4V-Al2O3板的大挠度弯曲问题,证明了该方法在求解FGP非线性弯曲问题时的收敛性,并研究了体积分数指数n和宽厚比对FGP的无量纲挠度的影响.研究发现,板的无量纲挠度会随着n的增加而增加,而随着板的宽厚比的增加,板的无量纲挠度会趋近一个定值.这表明了S-R无网格法在求解FGP大挠度弯曲问题时的合理性.