梯度有限元法
开篇：润墨网以专业的文秘视角，为您筛选了一篇梯度有限元法范文，如需获取更多写作素材，在线客服老师一对一协助。欢迎您的阅读与分享！
【摘 要】功能梯度材料具有连续梯度变化的微观结构，其力学性能也因此呈梯度变化。为了能够使用有限元法对其进行力学分析，梯度单元因此出现，能够适应单元内材料属性的剧烈变化。研究了一种可行的有限单元，编制了四节点四边形梯度单元的有限元程序，计算了相关算例，并同弹性理论的精确解进行比较。
【关键词】功能梯度材料；有限元法；梯度单元
0 引言
功能梯度材料（Functionally Graded Materials， FGM）是20世纪80年代末出现的一种新型复合材料，具有连续梯度变化的微观结构，其属性也因此呈梯度变化，可以拥有强度和韧性两种特性[1]。根据FGM制备方法的不同，可使其表现出各向同性或者各项异性的材料特性。
有限元法是当今进行工程分析的最广泛通用的方法，也是功能梯度材料结构分析方法的发展方向之一。现有的商用有限元软件的一个缺陷就是不能够适当处理对单元内材料属性的剧烈变化。虽然能够使用细化的单元来保证精确度，但是为了能够在使用稀疏的单元下保证精确度，梯度有限元法因此出现，能实现单元内的属性变化。
1 梯度有限元法算法实现
在有限元法的等参单元中，使用形函数来表示单元内部的坐标等[2]。以等参四边形四节点平面单元为例，单元内任一点的位移可以表示为：
u=■Niui
其中ui为节点的位移矢量，Ni为形函数，在此
Ni=（1+ξξi）（1+ηηi）/4，i=1，2，3，4
其中（ξ，η）是在区间[-1，1]内的固有坐标系，（ξi，ηi）表示i节点的局部坐标。根据应变和位移的关系，可以得到应变的表达式：
ε=Bu
其中B是关于形函数Ni偏导数的几何矩阵。
应力的应变的关系可以通过应力-应变矩阵D（x）表示为：
σ=D（x）ε
对于一般的均匀材料来说，每个单元的D（x）矩阵都是由一些材料属性常数表示的。然而，对于FGM来说，其材料属性在空间上呈梯度变化，所以，每个单元的D（x）矩阵都应该是与坐标相关的函数。Santare MH[3]和Paulino GH[4-5]研究了梯度单元本构阵D（x）的表示方法。Santare MH是直接计算高斯点的弹性模量值来计算本构阵。而Paulino GH通过等参单元的基本概念，使用相同的形函数Ni来表示弹性模量和泊松比：
E=■NiEi，？淄=■Ni？淄i
在实际计算过程中，可以假设泊松比？淄在整个区域上为常数。这种计算方法的精度要高于Santare MH所采用的直接取高斯点的属性值的计算方法。
有限元法中单元刚度矩阵Ke定义为节点位移矢量u和节点载荷f e之间的关系矩阵：
f e=Keu
同时根据虚功原理有：
■
对于平面单元来说：
■
一般情况下，Ni以及■，■等皆为？孜，η的函数，因此B、D也为？孜，η的函数。同时利用雅可比矩阵J可以将面积元素ds也应以d？孜，dη表示为：
ds=dxdy=Jd？孜dη
此时，单元的刚度矩阵可以写成在自然坐标系内进行积分的形式：
■
有限元法中都采用数值积分求此式的近似值。等参单元计算中，为了减少计算点的数目和便于编制程序，采用高斯数值积分方法。
本文梯度有限元程序使用Fortran开发，MSC/Patran作为前后处理的程序
2 算例一
算例一模型的几何尺寸及约束如图1（a）所示，为无量纲模型。几何尺寸为9×9，分成9×9个单元。材料弹性模量沿x方向梯度变化，分布函数为E（x）=E0eβx，其中β=0.2。整个区域上泊松比为常量0.3。考虑如图1（b）和（c）的两种受力情况：均布载荷和弯矩载荷。
（a）几何尺寸
（b）均布载荷 （c）弯矩载荷
图1 算例模型
Kim JH和Erdogan F[6]推导了对于当y方向为无穷大时的理论解，作为梯度单元正确性的参考。
图2 均布载荷时沿x轴的σy分布情况
图2显示了均布载荷时，沿x轴的σy的计算结果。图3是弯矩载荷时的计算结果。实线为根据弹性理论解获得的结果曲线，‘*’为梯度单元求解的结果。可以发现，无论是均布载荷或是弯矩载荷，梯度有限元法的结果与弹性理论解非常吻合，具有相当好的计算精度。
图3 弯矩载荷时沿x轴的σy分布情况
3 算例二
Venkataraman S和Sankar BV[7]使用弹性理论求解了梯度圆环板在轴对称载荷和非轴对称载荷情况下的静力问题。
对于一个内外半径分别为r1和r2的圆环板模型，内外边界上具有轴对称的轴向均布载荷p1和p2。假设圆环板是各项同性材料，弹性模量在径向梯度变化，分布函数为E=E0rβ，在整个区域上泊松比为常数。
可以求得应力分布的理论解为：
■
其中A1、A2、λ1、λ2都是求解得到的常数。
在此，
r1=5，r2=15
p1=0，p2=1
E0=1，ν=0.25
β=1
求得A1=5.4458，A2=0.3632，λ1=-1.5，λ2=0.5。
同时，使用梯度单元法对其进行分析，取1/4模型进行建模，径向为10个单元，圆周上也为10个单元。
图4所示为弹性理论和梯度单元求解得到的径向应力分布。可以发现，梯度单元法求得的结果与理论解还是相当吻合的，但在圆环板的内边界上，σr有非常大的误差，可以认为在内边界上的结果是错误的。引起此处这么大错误的原因可能在于此处边界条件的处理上。
图4 圆环板径向应力分布
4 结论
一些简单外形的模型可以采用均匀化的单元、赋予不同单元不同的属性，利用现有的有限元程序进行分析。在弹性模量变化较剧烈的地方，可以采用细分网格的方式提高精确度。但是对于今后实际中采用功能梯度材料制造的零件来说，其外形不会有那么简单。此时，只需要指定材料属性变化函数就能进行有限元分析的梯度单元法的优势就体现出来了，将比均匀化单元的办法更加方便、结果也更加精确。
【参考文献】
[1]Miyamoto Y， Kaysser W A， Rabin B H， et al. Functionally Graded Materials： Design， Processing and Applications[M]. Kluwer Academic Publishers， 1999.
[2]Zienkiewicz OC， Taylor RL. The finite element method， Vol.1：The basis[M]. The fifth edition. Butterworth-Heinemann， 2000.
[3]Santare MH， Lambros J. Use of Graded Finite Elements to Model the Behavior of Nonhomogeneous Materials[J]. Journal of Applied Mechanics， 2000，67：819-822.
[4]Kim JH， Paulino GH. Isoparametric graded finite elements for nonhomogeneous isotropic and orthotropic materials[J]. Journal of Applied Mechanics， 2002，69：502-514.
[5]Buttlar WG， Paulino GH， Song SH. Applied of Graded Finite Elements for Asphalt Pavements[J]. Journal of Engineering Mechanics， 2006， 132： 240-249.
[6]Erdogan F， Wu BH. The Surface Crack Problem for a Plate with Functionally Graded Properties[J]. ASME Journal of Applied Mechanics， 1997， 64：449-456.
[7]Venkataraman S， Sankar BV. Elasticity analysis and optimization of a functionally graded plate with hole[C]//44th AIAA/ASME/ASCE/AHS Structures， Structural Dynamics， and Materials Conference. Norfolk， Virginia， 2003.