US－FE－LSPIM四边形单元性能进一步研究★

贾 程

(盐城工学院土木工程学院，江苏 盐城 224051)

0 引言

有限单元法已经成为工程领域中一个重要的分析工具。但是，传统的等参元在网格畸变时精度下降给人们带来很大的麻烦。在计算大变形问题时，由于精度下降还需要重新划分网格。为提高有限单元法的精度，寻找更先进的方法，学者们做了大量的研究工作。提出了QM6单元、增强假设应变法(EAS)、四边形面积坐标法和无网格法(EFG)[1]等。

最近，RAJENDRAN 等[2]提出的 FE－LSPIM 四边形单元属于杂交单元。该单元结合有限单元法和无网格法的优点，具有高精度和Kronecker－delta性质。但是，该单元施加单元整体边界需满足的位移不方便，需要用罚函数或拉格朗日乘子法。文献[3]改进了它在单元整体边界位移施加中的不足，提出US－FE－LSPIM四边形单元。本文进一步研究US－FE－LSPIM四边形单元在网格畸变、剪切自锁和体积自锁时的性能。

1 US-FE-LSPIM四边形单元

1.1 FE－LSPIM四边形单元的形函数

单元内位移u写成:

其中，N′=[N1N2N3N4];Ni为四边形等参元形函数，i=1，2，3，4。

ui(x，y)为节点位移函数，i=1，2，3，4，在该节点处等于其位移值，即:ui(xi，yi)=ui，i=1，2，3，4，位移函数 ui(x，y)由 i点的支持域内的节点值运用最小二乘点插值法(LSPIM)得到:

其中，Φi=[Φi1Φi2Φi3… Φi

N]，i=1，2，3，4。

Φi为LSPIM法的关于节点i的形函数矩阵;ui为节点i支持域内节点的位移参数向量;N为节点i支持域内的所有节点数。对于该单元内的其他节点，N可能不同。一个单元所有节点支持域的合集构成了一个单元的支持域Ω。

节点支持域和单元支持域的定义如图1所示。

设单元支持域Ω内的节点数为M，则1≤N≤M，方程(3)就可以写成单元支持域的形式:

其中，Φ为相应于单元支持域的LSPIM单元形函数矩阵;u为单元支持域内所有节点的x方向的位移向量。

图1 支持域的定义

将式(4)代入式(1)得:

从式(5)中，得到该单元的形函数矩阵:

Φ中的元素Φi可以由i点的支持域内的节点值运用最小二乘点插值法(LSPIM)得到:

其中，¯Q，¯q定义见文献[2];1为所有元素均为1的列向量。则利用式(6)进而可以求出单元的形函数矩阵ψ。

1.2 US－FE－LSPIM 四边形单元

线弹性体在平衡状态下的虚功方程为:

其中，b为体力;t为面力;δu为虚位移场;δε为相应的虚应变场;σ为真实的应力场。

为了重构一个完整的二次位移场，形函数需满足以下方程:

FE－LSPIM四边形单元形函数ψ能够满足式(9)所有的方程[2]，而常规的四边形等参元形函数只能满足低阶的条件，即p，q取0或1。

由于方程(8)中的虚位移可以是任意的，只要它满足本质边界条件和域内连续性条件。常规的四边形等参元形函数插值的位移场δ¯un由于能够满足单元间C0阶连续和单元内C1阶连续的要求，故是个合适的选择。

在传统的等参元中，σ由等参元形函数插值的应力¯σ代替。但是由于不能满足式(9)中的高阶条件，导致等参元在网格畸变时精度下降。为了提高在单元二次位移域下的计算精度，选择FE－LSPIM插值的应力场作为试函数，即

将以上的虚位移和应力函数代入式(8)，整理可得:

2 数值算例

2.1 划分为五个畸变网格的悬臂梁

如图2所示，悬臂梁被划分成五个畸变的网格，弹性模量E=1 500，泊松比υ=0.25，厚度t=1。考虑两种荷载工况:a.弯矩 M作用下的纯弯曲;b.横向剪力作用下的线性弯曲。A点的竖直位移列于表1。

图2 划分为五个畸变网格的悬臂梁

表1 悬臂梁A点的竖直位移

从表1可以看出US－FE－LSPIM四边形单元在畸变的网格下具有较高的精度。

2.2 剪切自锁问题

为研究剪切自锁问题，考虑如图3所示悬臂梁。弹性模量E=107，泊松比 υ =0.25，厚度 t=0.1。弯矩 M=h，长宽比 L/h 分别取3，30，300和3 000，使用2×2高斯积分。梁端的竖直位移的计算结果列于表2。

图3 悬臂梁剪切自锁问题

表2 悬臂梁端部位移

从表2中可以看出，在弯矩作用下，随着长宽比的增大，四边形等参元Q4的精度急剧下降，表现出剪切锁定现象。QM6单元和US－FE－LSPIM单元在长宽比3～300时，精度变化不大;在长宽比为3 000时，QM6单元精度下降较大，而 US－FE－LSPIM单元精度略有下降，显示了较强的抗剪切自锁性能。

2.3 体积自锁问题

如图4所示的悬臂梁，受端部弯矩作用，考虑平面应变问题。弹性模量 E=1 500，泊松比 υ 分别取 0.25，0.49，0.499，0.499 9。e分别取0和3。A点的竖直位移计算结果列于表3。

图4 悬臂梁体积自锁问题

表3 悬臂梁体积自锁问题A点的竖直位移

从表3中可以看出，Q4四边形单元在e等于0和3时都表现出明显的体积自锁现象。和QM6单元一样，US－FE－LSPIM单元无体积自锁现象。但在e=3时，US－FE－LSPIM单元的精度高于QM6单元。

3 结语

本文进一步研究了US－FE－LSPIM单元的静力性能。通过数值算例表明，该单元在畸变网格下具有较高的精度，优于Q4，QM6单元，并且该单元显示了较强的抗剪切自锁和体积自锁时的能力。

[1] BELYTSCHKO，LU Y Y，Gu L.Element free Galerkin methods[J].International Journal for Numerical Methods in Engineering，1994，37(2):229－256.

[2] RAJENDRAN S，ZHANG B R.A“FE－Meshfree”QUAD4 element based on partition of unity[J].Computer Methods in Applied Mechanics and Engineering，2007，197(1):128－147.

[3] 贾 程，陈国荣，陈卉卉.US－FE－LSPIM 四边形单元及其在几何非线性问题中的应用[J].计算力学学报，2011，28(5):785－791.

[4] Chen XM，Cen S，Long YQ，et al.Membrane elements insensitive to distortion using the quadrilateral area coordinate method[J].Computers and Structures，2004，82(1):35－54.

[5] 铁摩辛柯，古地尔.弹性理论[M].第3版.徐芝纶，译.北京:高等教育出版社，1990.