广义有限单元法研究进展

吴继华，王志非 （同济大学土木工程学院，上海 200092）

1 概 述

传统的有限单元法（FEM）理论成熟，并且有许多商业软件可用，是目前在工程上应用最广泛的数值计算方法。但传统FEM对单元划分依赖度高，在处理大变形或裂纹扩展等问题时很难划分出或自动生成合理的单元以满足计算精度要求。无网格法克服了FEM对网格的依赖性，在求解涉及网格畸变、网格移动等问题时有明显的优势。在无网格法的研究中，单位分解法的出现使得在近似函数中可以加入能够反映待求解问题特性的函数，从而也就催生了广义有限单元法（GFEM）。GFEM采用分片任意高阶多项式甚至是任意函数的线性组合作为逼近空间，因此GFEM可以通过提高广义结点函数的阶数和引入特定的局部逼近函数，在不增加结点个数的前提下，有效地提高总体近似解的精度。近年来GFEM的理论基础和应用都得到了迅猛发展。

2 广义有限单元法的研究历史及现状

GFEM的概念在国外最早由Babuska I等人提出[1]，但只有单位分解法[2]出现后才开始了系统的GFEM 研究，T·Strouboulisd等将FEM的插值函数作为单位分解法的单元分解函数，可见GFEM是FEM和无网格法中的单位分解法的综合[3]。Babuska I、C·A·Duarte、Barros FB等人发表了一系列论文[4~7]，详细阐述了GFEM的基本原理、收敛性和误差估计等问题，并将GFEM应用在三维结构力学、动态裂纹发展、结构非线性等领域，使GFEM的理论基础和应用都得到迅速发展。我国学者梁国平等吸收了流形方法的思想，将每个结点只有一个位移未知量的拉格朗日型插值空间推广为每个结点有任意多个广义位移的函数展开式，也提出了广义结点有限元法[8]，并从数学上论证了它的可行性。栾茂田[9][10]等人对上述方法进行了深入研究，给出了平面问题一阶和二阶GFEM的具体形式。程尧舜、辛娅云[11]对栾茂田等人的方法进行了改进，将坐标原点移到考察结点上，给出了插值函数的新形式，使得广义自由度的意义得到明确，边界条件也变得容易处理。程尧舜、李凤岭[12]详细阐述了广义自由度置零法、罚函数法等处理GFEM位移边界条件的方法，并提出了一种新的处理GFEM位移边界条件的方法，提高了计算精度。

GFEM相比较传统的FEM而言，具有其自身的优越性，但是GFEM每个结点有多个自由度，增加了计算量。程尧舜、李刚[13]针对这一缺陷提出了一种新型的GFEM，称之为半弥散单元法（SDEM）。SDEM利用考察结点的邻近结点的函数值来表示结点的广义自由度，使得每个结点的自由度数和FEM相同，在系统自由度数不增加的前提下可以大大提高计算精度。S·Rajendran[14]结合FEM中形函数和无网格法中的单位分解法也提出了一种新的四边形单元，这种方法和SDEM有一定的相似之处，也能使得自由度数降低，但是和SDEM相比计算量更大，四边形单元在单元划分时适应性也比三角形单元低。黄书彪[15]通过引进更高阶的插值函数对SDEM改进，提高了算法精度，并且提出了新的处理位移边界条件的方法，通过改变已知位移边界结点的邻近点，使得在SDEM中可以使用FEM相同的方法来处理位移边界条件。B·R·Zhang、S·Rajendran[16][17]将新提出的四边形单元应用于非线性问题和动力问题，得到了较好的结果。黄亮[18]将SDEM应用于求解二维稳态热传导问题中，数值算例的计算结果表明，在相同网格划分下SDEM的精度远高于FEM，体现了SDEM在计算中的优势。

3 半弥散单元法的基本原理

设单元内的真实场函数是一个二次完全多项式：

式中：a、b、c、d、e、f均为常数。那么结点i 处的场函数值和场函数的一阶偏导数值如下：

下面分别用GFEM和无网格法的思想构造局部近似场函数。首先考查GFEM中插值函数的特性。如图1所示，三结点三角形单元的3个结点编码按逆时针依次为i、j、m，为了减小条件数，将坐标原点移到考察结点上，则未知场的二阶单元插值函数表示为如下形式：

式中：x、y 为坐标；Fh为表示单元局部场函数；di1、di2、di3为广义自由度；Li为传统有限单元法中三结点三角形单元的插值函数（见图1）。

由式(4)可知，GFEM中每个结点的自由度为FEM的3倍，如果广义自由度不做处理，计算量将大大增加。考察单元插值函数与真实场函数的联系，将结点i的坐标代入式(4)可得到结点i处场函数值Fh=di1，因此di1可确定为：

图1 三角形单元

将式(5)和式(6)代入式(4)，并利用FEM中三结点三角形单元的行函数Li的性质可得：

由式(7)可知，单元插值函数可以精确的表达二阶真实场函数。其中参数di1即为该结点的场函数值，di2,di3分别该结点场函数对x、y的偏导数值的一半，即：

其次，从无网格法的思想出发，构造结点i附近的局部近似场函数，取二阶表达式如下：

将结点i的坐标和场函数值代入上式可得：

在结点i附近选取5个邻近点，将邻近点坐标和场函数值代入式(9)可得一组线性方程组：

式中，xij、yij是第j个邻近点的位置坐标；Fij为第j个邻近点函数值，j为邻近点编号，j取1、2、3、4、5。参数中已由式(10)确定，另外5个参数可以由线性方程组(11)求得。因此参数已经可以完全由结点的场函数值确定。

考察结点局部近似场函数与真实场函数的联系，结点i处场函数值和一阶导数值为：

由结点i局部近似场函数真实场函数求得的函数值和场函数的导数值应该和由求得的值相同，因此：

比较(8)和(13)两式，可以建立两种局部近似场函数中广义自由度的联系：

根据上式，GFEM的广义自由度可以由无网格法中广义自由度来表示，而无网格法中广义自由度可以由结点的场函数值来表达。因此，GFEM中的结点广义自由度可以由结点邻近点的场函数值解出。和一般的GFEM不同，SDEM正是利用这一联系，在求解之前已经把结点广义的广义自由度用结点邻近点的场函数值表示出来，使得自由度大大减少。

4 结 论

GFEM就其本质来讲是在FEM的结点上增加新的自由度，再匹配以可以表现问题特性的插值函数，这就使得在遇到复杂问题的时候GFEM可以不增加结点而只需要增加自由度以达到足够的精度，避免了传统有限单元法的网格重构过程。相比FEM，广义结点有限元具有其自身的优越性，但是每个结点有多个自由度，这增加了数据的储存量，增加了计算量。SDEM的提出给广义有限单元法的研究带来了新的思想，利用考察结点的邻近结点的函数值来表示结点的广义自由度，使得每个结点的自由度数和FEM相同，在系统自由度数不增加的前提下可以大大提高计算精度。

[1]Babushka I,Osborn.Generalized finite element methods:their performan-ce and their relation to mixed methods[J].SIAM Journal of Numerical Analysis,1983(3).

[2]Babuska I,Melenk JM.The partition of unity method[J].Int.J.Numer.Meth,Engng,1997(40).

[3]Strouboulis T,Babuska I,Copps K.Thedesign and analysisof the Genera-lized Finite Element Method[J].Comput.Methods Appl.Mech.Engrg.2000(181).

[4]Duarte CA,Babuska I,Oden JT.Generalized finite element methods for three-dimensional structural mechanics problems[J].Computer&Struc-tures,2000(77).

[5]Duarte CA,Hamzeh ON,Liszka TJ,et al.A generalized finite element method for thesimulation of three-dimensional dynamic crack propagation[J].Comput.Methods Appl.Mech.Engrg,2001(15-17).

[6]Strouboulis T,Babuska I,Hidajat R.The generalized finite element meth-od for Helmholtzequation:Theory,computation,and open problems[J].Comput.Methods Appl.Mech.Engrg,2006(37-40).

[7]Barros FB,Proenca SPB,de Barcellos CS.Generalized finiteelement met-hod in structural nonlinear analysis:ap-adaptivestrategy[J].Computational Mechanics,2004(2).

[8]梁国平,何江衡.广义有限元方法——一类新的逼近空间[J].力学进展,1995(4).

[9]栾茂田,田荣,杨庆.广义结点有限元法[J].计算力学学报,2000(2).

[10]田荣,栾茂田,杨庆.高阶形式广义结点有限元法及其应用[J].大连理工大学学报,2000(4).

[11]辛娅云.改进的广义结点有限元法[D].上海:同济大学,2006.

[12]李凤岭.广义结点有限元法中位移边界条件实施研究[D].上海：同济大学,2009.

[13]李刚.半弥散单元法[D].上海：同济大学，2007.

[14]S·Rajendran,B·R·Zhang.A“FE-meshfree”QUAD4 element based on partition of unity[J]Comput.Methods Appl.Mech.Engrg,2007(197).

[15]黄书彪.改进的半弥散单元法[D].上海：同济大学，2008.

[16]B·R·Zhang,S·Rajendran.“FE-Meshfree”QUAD4 element for freevibrationanalysis[J]Comput.Methods Appl.Mech.Engrg.2008(197).

[17]S·Rajendran,B·R·Zhang.A partition of unity-based‘FE-meshfree’QUAD4 element for geometric non-linear analysis[J]Int.J.Numer.Meth.Engng,2010(82).

[18]黄亮.半弥散单元法在维稳热传问题中的应用[D].上海：同济大学，2013.