使用扩展的无网格有限元法对结构性能的数值仿真

摘要

本发明公开了在压缩和/或近不可压缩区域数值仿真工程产品的结构性能的系统、方法和软件产品。无网格扩展有限元法(ME-FEM)用于该数值仿真。ME-FEM要求采用包含多个有限元的FEM模型表示工程产品。在ME-FEM中使用的有限元通常是低阶有限元。FME模型中的每个有限元是由位于该元素的域中的至少一个无网格扩展(ME)节点扩展的。每个ME节点对于其所属的元素具有独立于其他角节点的自由度。将基于位移的第一阶凸无网格近似应用到该ME节点。该凸无网格近似在每个元素的边界具有克罗内克-δ特性以允许应用本质边界条件。该ME-FEM元的梯度矩阵满足积分约束条件。ME-FEM插值是元素级的无网格插值，其在不可压缩界点是离散无发散的。

说明书

技术领域

本发明一般涉及工程产品或部件的计算机辅助工程分析(CAE)，特别涉及 使用扩展的无网格有限元法(ME-FEM)对结构性能进行数值仿真的方法。

背景技术

传统意义上讲，FEM或有限元分析(finite element analysis，FEA)是工 程师或者科学家最常用来对与诸如三维非线性结构设计和分析等复杂系统相 关的工程问题进行建模和求解的一种计算机辅助工程工具。有限元分析的名称 源于对所关注的目标物体的几何特征进行描述的方式。随着现代数字计算机的 出现，有限元分析已在有限元分析软件中实现。基本上，有限元分析软件提供 有关于几何描述的基于网格的模型、以及该模型中每一点处的相关材料性能。 在此模型中，所分析的系统的几何特征被表示为各种大小的实体、壳体和梁， 这些被称为元。各元的顶点被称为节点。该模型是由有限数目的元组成的，这 些元都被分配有一材料名称以便于与材料性能相关联。因此，该模型表示了被 分析的目标物体所占的物理空间以及它的周围环境。有限元分析软件随后涉及 一列出了每种材料性能(例如，应力-应变本构方程、杨氏模量、泊松比和热 传导率)的表格。此外，指定了目标物体的边界条件(即负荷、物理约束等)。 遵循此种方式，建立了目标物体及其环境的模型。

然而，在某些条件下，FEM其自身也存在缺陷。例如，在闭联集 (continuum)的数值相容条件和物理相容条件并不完全相同时，FEM并非有 利的。特别地，在拉格朗日计算中，FEM可能经历网格畸变，这将完全终止 计算或导致精确度急剧变差。另外，FEM在求解与高梯度或是独特的局部特 征相关的问题时，需要非常精细的网格，而这样的计算是非常昂贵的。在被称 为任意拉格朗日-欧拉(Arbitrary Lagrangian Eulerian，ALE))表示的方法中， 其目的是独立于材料移动网格以最小化网格畸变。然而，对于较大应变或是高 速碰撞仿真来说，该畸变通常还是会造成一些数值误差。此外，在数值仿真中， 网格可携带固有偏差(inherent bias)，而该固有偏差将妨害计算过程。另一实 施例是应变局部化问题，其将严重影响网格对齐灵敏度。因此，通过一组无网 格约束的节点来离散闭联集可有效减少计算量。

因此，近年来也采用无网格方法消除FEM的这些缺陷。然而无网格法自 身也存在问题。例如，与FEM相比，无网格法的计算成本通常很高。而且， 在无网格法中，应用本质边界条件(essential boundary condition)不但困难， 而且计算成本很高。

无论是FEM还是无网格法，均能用于各种环境负载下的每个可能结构性 能的数值预测，例如称为体积闭锁的数值表现，其用于在可压缩和/或近不可 压缩区域(如橡胶、经历大塑性变形的材料)仿真结构性能。图1A示出了体 积闭锁效应的简化二维示例。元素(element)102的四个侧边中有两个侧边是 固定的。因此，节点112是空间固定的，其不能在数值仿真中移动。然而这并 不表示物理现实(physical reality)。

可使用下述一组数学推导过程描述体积闭锁：

把占据中多边形有界域Ω的平面均匀的各向同性性线性材料体看做具 有边界将该问题约束到较小变形且将无穷小应变张量ε定义成位移u 的函数：

$ϵ\left(u\right)=\frac{1}{2}(▿u+{▿u}^T)---\left(1\right)$

对于预定质量力(body force)f，Ω中的控制平衡方程式为

$▿·\sigma +f=0---\left(2\right)$

在此，σ是对称柯西应力张量。4阶弹性张量表示为C，给出Ω中本构方程 (constitutive equation)如下：

σ＝Cε＝2με+λ(trε)I    (3)

在此，I是4阶的恒等式张量。符号μ和λ为拉梅(Lamé)常数，其通过 下列方程与杨氏模量E和泊松比v相关：

$\mu =\frac{E}{2(1+v)},$$\lambda =\frac{\mathrm{vE}}{(1+v)(1-2v)}---\left(4\right)$

简而言之，假定该位移满足Γ_D上齐次狄利克雷(Dirichlet)边界条件u＝0， 且牵引力t在纽曼边界Γ_N上以Γ＝Γ_D∪Γ_N和Γ_D∩Γ_N＝0分布。该问题的变形是找 出位移u∈V，这样

$A(u,v)=l\left(v\right)\forall v\in V---\left(5\right)$

在此，空间由索伯列夫(Sobolev)空间H¹(Ω)中在边界上为0 的函数组成，并且可以定义如下：

V(Ω)＝{v：v∈H¹，在Γ_D上v＝0}    (6)

方程5中的双线性型A(·，·)和线性函数l(.)定义如下：

$A:\mathrm{VxV}\to A(u,v)=2\mu {\int }_{\Omega }ϵ\left(u\right):ϵ\left(v\right)\mathrm{d\Omega }+\lambda {\int }_{\Omega }(▿·u)(▿·v)\mathrm{d\Omega }---\left(7\right)$

上述方程中的双线性型A(·，·)在V上是对称的、V-椭圆形并且连续的。使 用科恩不等式和Lax-Milgram定理可获得该问题的唯一解。

接着可采用方程(5)的变形公式在有限维子空间上采用公式表示 标准伽辽金(Galerkin)法以找出u^h∈V^h，使得

$A(u^h,v^h)=l\left(v^h\right)\forall v^h\in V^h---\left(9\right)$

方程(7)右手侧的第二项类似于如斯托克斯问题中的经典补偿法(penalty  method)中的补偿因子。当λ→∞(或者v→0.5)，方程(7)中的能量保持有限， 则需执行如下约束条件

对于u∈V，$▿·u=0---\left(10\right)$

与有限维空间中类似，我们可获得

对于u^h∈V^h，$▿·u^h=0---\left(11\right)$

对位移场使用低阶逼近获得的方程(9)的解并不满足方程(11)中的约 束条件，且当v→0.5时，不能均匀收敛。这是近不可压缩问题中的体积闭锁。 换句话说，当近似空间V^h不是富余到可使得该近似能够验证无发散条件 时，发生体积闭锁。图1B中示出了使用EFG(现有技术的无网格法) 的3×3离散的8个非物理模式。

可采用方程(9)的左手侧(LHS)的特征值分析来进一步揭示基于有限 元和无网格法的低阶位移中的闭锁现象。对一个具有2×2高斯积分点的矩形 四边形双线性(Q1)有限元进行了特征值分析。图1C和1D获得并示出了第 一和第二非物理闭锁模式。当v逼近0.5时，该特征值变得极大，这是不能物 理预测的。因此，可使用4个矩形Q1元素从特征值分析获得8个非物理闭锁 模式。与4个Q1有限元模式中的结果类似，通过图1B中所示的无元伽辽金 法使用具有等于1.5的归一化节点支撑的一阶近似可生成相同数量的非物理闭 锁模式。高斯正交定律的节点支撑和阶数的增加并不能完全消除EFG法中的 体积闭锁。非物理模式的数量保持不变。

基于前述缺陷、不足和问题，需要提出一种扩展的计算机执行的方法， 用于在压缩和/或近不可压缩区域数值仿真结构性能。

发明内容

本发明公开了一种系统、方法和软件产品，用于在压缩和/或近不可压缩 区域数值仿真工程产品的结构性能。

根据一方面，将一种扩展的无网格有限元法(ME-FEM)用于该数值仿真。 ME-FEM要求采用包括多个有限元(如，平面、实体等)的FEM模型来表示 工程产品。在ME-FEM中使用的有限元通常是低阶有限元(也就是，仅具有 角节点的有限元)。在FEM模型中的每个有限元都由至少一个附加节点扩展， 该附加节点被称为位于元素域中的无网格节点或是无网格扩展(ME)节点(也 就是，位于所述元素的外边界以内)。

根据另一方面，每个ME节点对于其所属的元素具有独立于其他角节点的 自由度。将被称作基于位移的第一阶无网格凸近似(convex meshfree  approximation)的特定方法应用到该ME节点。该无网格凸近似在元素的边界 具有克罗内克-δ(Kronecker-delta)特性，并且是用专门的被称为广义无网格 法的无网格处理来构建。ME-FEM元素的梯度矩阵满足积分约束条件 (integration constraint)。该ME-FEM插值是元素智能(element-wise)无网格 插值，其在不可压缩界点处离散无发散。无发散特性是通过四边形元素的无网 格扩展提供的。对于三角形和四面体形元素，无发散特性是通过无网格扩展和 应变平滑过程共同提供的。

根据另一方面，基于位移的第一阶凸无网格近似使用幂函数、反正切函数、 双曲正切函数或Renyi函数生成的基础函数。

可选择地，在另一方面，有限元可分成第一组和第二组。仅扩展属于第一 组的有限元(也就是，增加一个或多个ME节点)。换句话说，ME-FEM可包 括常规有限元(也就是，第一组)和无网格扩展有限元(也就是，第一组)。

根据又一方面，该ME-FEM元配置成避免在仿真近不可压缩区域中的结 构性能中的压力振荡。为了实现该目标，将无发散ME-FEM插值与区域加权 的积分相结合以在应变和压力上提供平滑效应。该区域加权的积分还特别在三 角形和四面体形元素中提供无发散特性。

根据再一方面，使用改进的胡海昌-鹫沣(Hu-Washizu)变形原理创建具 有平滑应变场的离散方程。

使用ME-FEM元，使得压缩和/或近不可压缩区域中的结构性能的数值仿 真能够以合理的精确度执行，从而有助于用户(也就是，工程师、科学家等) 改进工程产品(也就是，汽车、飞机等)的设计。

通过结合附图详阅接下来对实施例的详细描述，本发明的其它目标、特征 和优点将是显而易见的。

附图说明

本发明的这些以及其它特征、方面和优点通过以下描述、权利要求和附图 能得到更好的理解：

图1A是用以证实体积闭锁的有限元的示意图；

图1B是3x3离散的各个非物理模式的示意图；

图1C和图1D是Q1双线性有限元的两个非物理模式的示意图；

图2A是将实例域划分或三角划分成数个有限元的示意图；

图2B是根据本发明一个实施例的、在5-节点无网格扩展有限元(四边形 元素)中使用的示例型等参映射法则的示意图；

图2C是根据本发明一个实施例的、在4-节点无网格扩展有限元(三角形 元素)中使用的示例型等参映射法则的示意图；

图2D是根据本发明一个实施例的、在5-节点无网格扩展有限元(四面体 形元素)中使用的示例型等参映射法则的示意图；

图3A和3B是根据本发明一个实施例的、示例型无网格扩展有限元的形 状函数的示意图；

图4是根据本发明实施例的、内部节点和其相邻元素之间的关系的示意 图；

图5A和5B示出了根据本发明实施例的、使用原始的无网格法的非物理 模式和使用ME-FEM的现实模式；

图6A和6B是在本发明的实施例中使用的示例型区域加权的积分方法的 示意图；

图7是在包含区域加权的积分的四个ME-FEM元中的8个无闭锁模式的 示意图；

图8是根据本发明一个实施例的、使用ME-FEM在压缩和/或近不可压缩 区域中数值仿真结构性能的典型方法的流程图；

图9是典型的计算机系统的主要部件(salient component)原理框图，可 在其中执行本发明的实施例。

具体实施方式

本发明涉及使用ME-FEM数值仿真结构性能的方法和系统，在ME-FEM 中使用的一个技术被称为广义无网格(GMF)近似方法或步骤。

广义无网格(GMF)近似

引入第一阶GMF凸近似来近似ME-FEM的位移场。执行无网格凸近似 在每个有限元的边界应用弱克罗内克-δ特性，从而使得本质边界条件很容易处 理。GMF近似的基本思路是在谢泼德函数(Shepard function)和其线性相容 性的满意度中引入扩展基础函数。扩展基础函数的选择确定了GMF近似是否 具有凸性(convexity property)。节点组的凸包(convex hull) conv(Λ)定义如下：

GMF法用于构建给定平滑函数u的凸近似，该平滑函数u为 且具有满足多项式再生特性(polynomial reproduction  property.)的生长函数Ψ_i：conv(Λ)→R。典型的多项式再生特性表示成 在多维中的第一阶GMF近似可写成如下表达式：

${\Psi }_i(x,{\lambda }_r)=\frac{{\psi }_i}{\psi }=\frac{{\phi }_a(x;X_i){\Gamma }_i(x_i,{\lambda }_r)B(X_i,{\lambda }_r,\Theta )}{{\Sigma }_{j=1}^n{\phi }_a(x;X_j){\Gamma }_j(X_j,{\lambda }_r)B(X_j,{\lambda }_r,\Theta )}---\left(13\right)$

该表达式服从线性约束条件

$R_r(x,{\lambda }_r)=\underset{i=1}{\overset{n}{\Sigma }}{\Psi }_i{X}_i=0---\left(14\right)$

在此：

ψ_i＝φ_a(x；X_i)Γ_i(X_i，λ_r)B(X_i，λ_r，Θ)    (15)

$\psi ={\Sigma }_{i=1}^n{\psi }_i={\Sigma }_{i=1}^n{\phi }_a(x;X_i){\Gamma }_i(X_i,{\lambda }_r)B(X_i,{\lambda }_r,\Theta )---\left(16\right)$

X_i＝x-x_i    (17)

φ_a(x；X_i)：节点i的加权函数，具有支撑大小supp(φ_a(x；X_i))＝a_i，

Γ_i(X_i，λ_r)，GMF近似的基础函数，

x：内部点的坐标(也就是，固定点)，

x_i：节点i的坐标，

n：固定点x的支撑大小a(x)中的节点数量，

λ_r(x)(r＝1，2，Λ，m)：需要决定的约束参数，

m：约束的数量(在一维中m＝1，在二维中m＝2，在三维中m＝3)，

B(X_i，λ_r，Θ)：混合函数和Θ是内部函数或变量。

在GMF近似中，单位分解的特性是自动满足方程(13)中的归一化。引 入混合函数B(X_i，λ_r，Θ)以修正基础函数，用于生成具有特定目的(例如满足非 凸近似的弱克罗内克-δ(Kronecker-delta)特性)的特定GMF近似。在构建第 一阶凸近似的过程中，混合函数B(X_i，λ_r，Θ)是均一的(unity)。此外，通过找到 满足Eq.(14)的λ完成GMF近似，该方程Eq.(14)是确保GMF近似以具有线 性相容性的约束方程。为了在方程(13)中在任何固定点x确定λ，非线性基 础函数需要用到寻根(root-finding)算法。

GMF近似的凸性由在方程(13)中的正基础函数的选择确定。表格1概 述了通常用于生成凸或凹GMF近似的各种基础函数。基础函数或其他单调增 或减函数的组合可看做是其他种类的GMF函数。

表格1中的GMF(exp)近似已知为具有香农熵的无网格扩展(ME)近似， 而该GMF(MLS)近似是如移动最小二乘法(MLS)或再生核(RK)近似法。此 外，GMF(Renyi)近似是具有Renyi熵的ME近似。此外，GMF近似可延展到 高阶，其可用于改进无网格法的精确度。使用反正切函数构造的凸GMF近似 可在ME-FEM法中的位移近似中使用，并且在表格1中由GMF(atan)表示。 将具有矩形支撑的立方样条窗(三次样条的窗口)选作方程(15)中的加权函 数。该使用矩形支撑是因为其可以适合四边形元素等参映射域。

表格1-在GMF映射中的基础函数示例

GMF近似的扩展有限元插值

令M_h＝∪_eQ_e为图2A中示出的域Ω202到凸四边形元素204(只示出少量) 的规则划分(regular division)。每个元素Q_e 204包含4个角节点N_i，1≤i≤4 206a-d和至少一个无网格扩展(ME)节点208K_i，1≤i≤nb。该数量nb是元素 Q_e204中扩展点的数量。该元素Q_e 204也可通过以如下矢量形式表示的等参 映射变形成基准正方形214中的图像：

x＝F_e(ξ)                 (18)

在此x＝[x，y]^T且ξ＝[ξ，η]^T.

映射函数F_e220映射基准正方形214到常规非规则元素Q_e204。为了 解释和描述上的简明，图2B示出了5节点无网格扩展(ME)有限元204，其 中第5节点K₁208位于该5节点无网格扩展(ME)有限元204的中心(也就 是，对该例子来说，nb等于一(1))。使用不同数量的节点和不同的无网格凸 近似的各种无网格扩展方法都可用来实现基本类似的目的。

由于元素的形状函数可通过无网格凸近似构建，形状函数显示下列凸特 性：

Ψ_i(ξ_j)＝δ_ij  1≤i，j≤4

${\Psi }_5\left(\xi \right)>0\forall \xi \in \overline{Q_e}---\left(19\right)$

${\Psi }_5\left(\xi \right)=0\forall \xi \in \partial \overline{Q_e}$

假定每个节点分配给在基准正方形214中具有相同矩形支撑 supp(φ_a(ξ；ξ_i))＝a_i＝a的加权函数，方程(15)中的该张量积加权函数φ_a(x；X_i)如 下定义。

${\phi }_a(\xi ;{\xi }_i):\overline{Q_e}\to \left[\mathrm{0,1}\right]{\phi }_a(\xi ;{\xi }_i)=\phi \left(\frac{\xi -{\xi }_i}{a}\right)\phi \left(\frac{\eta -{\eta }_i}{a}\right)---\left(20\right)$

基于上述假定，在外形函数和其导数中，存在某些对称特性。使用5节点 ME-FEM形状函数时，方程(18)中的等参映射函数可如下定义：

对于全部的$\xi \in \overline{Q_e},F_e:\overline{Q_e}\to Q_e,$$x=F_e\left(\xi \right)=(F_e^1,F_e^2)=(\underset{i=1}{\overset{5}{\Sigma }}{x}_i{\Psi }_i(\xi ,\eta ),\underset{i=1}{\overset{5}{\Sigma }}{y}_i{\Psi }_i(\xi ,\eta ))---\left(21\right)$

采用上述符号，位移场中的下述近似空间可定义如下：

对于全部的$Q_e\in M_h\},$

在此，表示该空间包含位于基准正方形214中的 一组基础函数。

ME有限元中的第5或ME节点提供函数，该函数与传统有限元方法中的 不相容应变元素或增强型应变元素中的bubble函数相似。然而，与第五节点 相关的基础函数并不与和角节点函数相关的基础函数正交，这是因为其依赖于 无网格形状函数。这样，该第5节点不能由ME有限元中的静凝聚法消除。与 传统或者经典bubble函数的系数未知不同，第5节点具有物理位移。在计算 ME有限元的质量力矢量时，必须考虑第5节点。根据方程(21)，在第5节 点上的等参映射为：

$(x_5,y_5)=F_e\left(\mathrm{0,0}\right)=(\underset{i=1}{\overset{5}{\Sigma }}{x}_i{\Psi }_i\left(\mathrm{0,0}\right),\underset{i=1}{\overset{5}{\Sigma }}{y}_i{\Psi }_i(0,0))---\left(23\right)$

因此，可确定第5节点的全局同等坐标(global co-coordinate)为：

$(x_5,y_5)=(\frac{\underset{i=1}{\overset{4}{\Sigma }}{\Psi }_i\left(\mathrm{0,0}\right)x_i}{(1-{\Psi }_5\left(\mathrm{0,0}\right))},\frac{\underset{i=1}{\overset{4}{\Sigma }}{\Psi }_i\left(\mathrm{0,0}\right)y_i}{(1-{\Psi }_5\left(\mathrm{0,0}\right))})=(\frac{\underset{i=1}{\overset{4}{\Sigma }}{\Psi }_i\left(\mathrm{0,0}\right)x_i}{\underset{i=1}{\overset{4}{\Sigma }}{\Psi }_i\left(\mathrm{0,0}\right)},\frac{\underset{i=1}{\overset{4}{\Sigma }}{\Psi }_i\left(\mathrm{0,0}\right)y_i}{\underset{i=1}{\overset{4}{\Sigma }}{\Psi }_i\left(\mathrm{0,0}\right)})=(\frac{\underset{i=1}{\overset{4}{\Sigma }}{x}_i}{4},\frac{\underset{i=1}{\overset{4}{\Sigma }}{y}_i}{4})---\left(24\right)$

其为四边形的质心(也就是，5-节点ME有限元204)。

对于4节点无网格扩展有限元(三角形元素224，该三角形元素224具有 三个示出为实心点的角节点以及一个示出为圆环的ME节点)。图2C中示出 了示范性的等参映射机制，可见，对于四边形元素204，圆环形支撑230可用 于三角形元素224以替代长方形支撑。

此外，本发明可应用到图2D中示出的三维元素(也就是，四面体元素 224)，其中示出示范性的等参映射机制。

图3A示出了角节点206a-d中的一个的示范性形状函数，而图3B示出了 无网格扩展(ME)节点208(也就是，中节点或第5节点)的示范性形状函 数。这些形状函数是5节点ME-FEM元素204的GMF(atan)近似，该5节点 ME-FEM元204在边界显示了克罗内克-δ特性。ME-FEM元素的形状函数沿 着显示凸近似特性的元素边降低到标准双线性有限元形状函数。应注意，5节 点ME-FEM元204中的形状函数(下列方程(79)中列出)的导数不能沿着 元素边简并成有限元的导数，这是因为ME或第5节点208的导数的非零值。

雅可比矩阵J_ξ通过下式定义成映射函数F_e的梯度

$J_{\mathrm{ij}}=\frac{\partial F_e^i\left(\xi \right)}{{\partial \xi }_j},i,j=\mathrm{1,2}---\left(25\right)$

或是以矩阵形式

$J_{\xi }=G_F{G}_F^T---\left(26\right)$

其中，矩阵G_F和G_x定义如下：

$G_F={▿}_{(\xi ,\eta )}{\left[{\Psi }_i(\xi ,\eta )\right]}_e=\left(\begin{array}{ccc}\frac{{\partial \Psi }_1}{\partial \xi }& K& \frac{{\partial \Psi }_5}{\partial \xi }\\ \frac{{\partial \Psi }_1}{\partial \eta }& K& \frac{{\partial \Psi }_5}{\partial \eta }\end{array}\right)---\left(27\right)$

$G_x=\left(\begin{array}{ccc}{x}_1& \Lambda & x_5\\ y_1& \Lambda & y_5\end{array}\right)---\left(28\right)$

${▿}_{(\xi ,\eta )}=\left(\begin{array}{c}\frac{\partial }{\partial \xi }\\ \frac{\partial }{\partial \eta }\end{array}\right)---\left(29\right)$

由于无网格形状函数的有理数特性，难以获取映射函数F_e的可逆性(或 雅可比的非奇异性)的分析证据。并且该分析证据需要采用反函数定理加以验 证。在工程应用中，更常见的是在用于积分ME元素204的高斯积分点验证该 可逆性。

5节点ME-FEM四边形元素204的等参映射数值保持全局元素区域A_e， 在此换句话说，虽然ME-FEM元素形状 函数和其导数是有理函数且节点支撑大小可调，但是使用2*2高斯积分算法的 det(J_ξ)的和保持全局元Q_e204的区域。

证据

从方程(26)，可使用下面给出的传统的2×2高斯积分获得全局元素Q_e的 区域：

$A_e=\underset{J=1}{\overset{4}{\Sigma }}\det {\left(J_{\xi }\right)}_J=\underset{J=1}{\overset{4}{\Sigma }}\det \left(G_F\left({\xi }_{\mathrm{gJ}}\right)G_x^T\right)$

$=\underset{J=1}{\overset{4}{\Sigma }}(\underset{I=1}{\overset{5}{\Sigma }}\frac{{\partial \Psi }_I\left({\xi }_{\mathrm{gJ}}\right)}{\partial \xi }{x}_I·\underset{I=1}{\overset{5}{\Sigma }}\frac{{\partial \Psi }_I\left({\xi }_{\mathrm{gJ}}\right)}{\partial \eta }{y}_I-\underset{I=1}{\overset{5}{\Sigma }}\frac{{\partial \Psi }_I\left({\xi }_{\mathrm{gJ}}\right)}{\partial \xi }{y}_I·\underset{I=1}{\overset{5}{\Sigma }}\frac{{\partial \Psi }_I\left({\xi }_{\mathrm{gJ}}\right)}{\partial \eta }{x}_I)---\left(30\right)$

第一阶相容条件给出：

$G_F\left({\xi }_{\mathrm{gi}}\right)G_{\xi }^T=I,i=\mathrm{1,2,3,4}---\left(31\right)$

$G_{\xi }=\left(\begin{array}{ccccc}-1& 1& -1& -1& 0\\ -1& -1& 1& 1& 0\end{array}\right)---\left(32\right)$

使用方程(79)中的对称特性和第一阶相容条件一起可获得A_e的下列显 式表达：

$A_e=\frac{(-x_1+x_2+x_3-x_4)·(-y_1-y_2+y_3+y_4)-(-y_1+y_2+y_3-y_4)·({-x}_1-x_2+x_3+x_4)}{2}---\left(33\right)$

A_e是全局元素Q_e204的精确区域(exact area)且独立于ME或第5节点 208的同等坐标或定位和节点支撑大小a。此外，ME-FEM公式(ME-FEM  formulation)在狄利克雷边界值问题的伽辽金近似中保持线性精确，也就是：

${\int }_{\Omega }{B}_I\mathrm{d\Omega }=0\mathrm{or}{\int }_{\Omega }{b}_I\mathrm{d\Omega }=0\forall$内节点I∈Ω∪Γ    (34)

其已知为无网格伽辽金法中的积分约束(IC)，在方程(34)中，B_I是 下式给出的梯度矩阵：

$B_I\left(\xi \right)=\left(\begin{array}{cc}\frac{{\partial \Psi }_I\left(\xi \right)}{\partial x}& 0\\ 0& \frac{{\partial \Psi }_I\left(\xi \right)}{\partial y}\\ \frac{{\partial \Psi }_I\left(\xi \right)}{\partial y}& \frac{{\partial \Psi }_I\left(\xi \right)}{\partial x}\end{array}\right)---\left(35\right)$

且，$b_I\left(\xi \right)=\left(\begin{array}{c}\frac{{\partial \Psi }_I\left(\xi \right)}{\partial x}\\ \frac{{\partial \Psi }_I\left(\xi \right)}{\partial y}\end{array}\right)$是节点I的形状函数导数的矢量。

证据

首先考虑内节点是全局元素Q_e204的中央节点K₁208的情况，在每个元 素中，

${\int }_{Q_e}{b}_5\mathrm{d\Omega }={\int }_{-1}^1{\int }_{-1}^1{b}_5\det \left(J_{\xi }\right)\mathrm{d\xi d\eta }$

$={\int }_{-1}^1{\int }_{-1}^1{J}_{\xi }^{-1}·{\bar{b}}_5\det \left(J_{\xi }\right)\mathrm{d\xi d\eta }$

$={\int }_{-1}^1{\int }_{-1}^1{J}_{\xi }^{*}·{\bar{b}}_5\mathrm{d\xi d\eta }$

$=\left(\begin{array}{c}\underset{j=1}{\overset{4}{\Sigma }}(\left(\underset{i=1}{\overset{5}{\Sigma }}{y}_i\frac{{\partial \Psi }_i\left({\xi }_{\mathrm{gj}}\right)}{\partial \eta }\right)\frac{{\partial \Psi }_5\left({\xi }_{\mathrm{gj}}\right)}{\partial \xi }-\left(\underset{i=1}{\overset{5}{\Sigma }}{y}_i\frac{{\partial \Psi }_i\left({\xi }_{\mathrm{gj}}\right)}{\partial \xi }\right)\frac{{\partial \Psi }_5\left({\xi }_{\mathrm{gj}}\right)}{\partial \eta })\\ \underset{j=1}{\overset{4}{\Sigma }}(\left(\underset{i=1}{\overset{5}{\Sigma }}{x}_i\frac{{\partial \Psi }_i\left({\xi }_{\mathrm{gj}}\right)}{\partial \eta }\right)\frac{{\partial \Psi }_5\left({\xi }_{\mathrm{gj}}\right)}{\partial \xi }-\left(\underset{i=1}{\overset{5}{\Sigma }}{x}_i\frac{{\partial \Psi }_i\left({\xi }_{\mathrm{gj}}\right)}{\partial \xi }\right)\frac{{\partial \Psi }_5\left({\xi }_{\mathrm{gj}}\right)}{\partial \eta })\end{array}\right)$

$=\left(\begin{array}{c}\underset{i=1}{\overset{5}{\Sigma }}\left(y_i\underset{j=1}{\overset{4}{\Sigma }}(\frac{{\partial \Psi }_i\left({\xi }_{\mathrm{gj}}\right)}{\partial \eta }\frac{{\partial \Psi }_5\left({\xi }_{\mathrm{gj}}\right)}{\partial \xi }-\frac{{\partial \Psi }_i\left({\xi }_{\mathrm{gj}}\right)}{\partial \xi }-\frac{{\partial \Psi }_5\left({\xi }_{\mathrm{gj}}\right)}{\partial \eta })\right)\\ \underset{i=1}{\overset{5}{\Sigma }}\left(x_i\underset{j=1}{\overset{4}{\Sigma }}(\frac{{\partial \Psi }_i\left({\xi }_{\mathrm{gj}}\right)}{\partial \eta }\frac{{\partial \Psi }_5\left({\xi }_{\mathrm{gj}}\right)}{\partial \xi }-\frac{{\partial \Psi }_i\left({\xi }_{\mathrm{gj}}\right)}{\partial \xi }\frac{{\partial \Psi }_5\left({\xi }_{\mathrm{gj}}\right)}{\partial \eta })\right)\end{array}\right)$

$=\left(\begin{array}{c}\underset{i=1}{\overset{5}{\Sigma }}{y}_i{g}_{i5}\\ \underset{i=1}{\overset{5}{\Sigma }}{x}_i{g}_{i5}\end{array}\right)---\left(36\right)$

在此，

${\bar{b}}_I={▿}_{(\xi ,\eta )}{\Psi }_I=\left(\begin{array}{c}\frac{{\partial \Psi }_I(\xi ,\eta )}{\partial \xi }\\ \frac{{\partial \Psi }_I(\xi ,\eta )}{\partial \eta }\end{array}\right)---\left(37\right)$

$J_{\xi }^{*}\left(\xi \right)=\det \left(J_{\xi }\right)J_{\xi }^{-1}\left(\xi \right)---\left(38\right)$

$g_{\mathrm{iI}}=\underset{j=1}{\overset{4}{\Sigma }}(\frac{{\partial \Psi }_i\left({\xi }_{\mathrm{gj}}\right)}{\partial \eta }\frac{{\partial \Psi }_I\left({\xi }_{\mathrm{gj}}\right)}{\partial \xi }-\frac{{\partial \Psi }_i\left({\xi }_{\mathrm{gj}}\right)}{\partial \xi }\frac{{\partial \Psi }_I\left({\xi }_{\mathrm{gj}}\right)}{\partial \eta })---\left(39\right)$

使用方程(79)中的对称特性，计算

${\int }_{{\Omega }_e}{b}_5\mathrm{d\Omega }=0---\left(40\right)$

接着，我们考虑数个相邻元素相连的内节点在此上标“e1“是指元 素数量，而下标“3”是指图4中示出的特定元素的节点数。与以上推导类似， 我们可获得：

${\int }_{{\Omega }_e}{b}_3^{e1}\mathrm{d\Omega }={\int }_{-1}^1{\int }_{-1}^1{b}_3^{e1}\det \left(J_{\xi }^{e1}\right)\mathrm{d\xi d\eta }$

$=\left(\begin{array}{c}\underset{j=1}{\overset{4}{\Sigma }}(\left(\underset{i=1}{\overset{5}{\Sigma }}{y}_i^{e1}\frac{{\partial \Psi }_i\left({\xi }_{\mathrm{gj}}\right)}{\partial \eta }\right)\frac{{\partial \Psi }_3\left({\xi }_{\mathrm{gj}}\right)}{\partial \xi }-\left(\underset{i=1}{\overset{5}{\Sigma }}{y}_i^{e1}\frac{{\partial \Psi }_i\left({\xi }_{\mathrm{gj}}\right)}{\partial \xi }\right)\frac{{\partial \Psi }_3\left({\xi }_{\mathrm{gj}}\right)}{\partial \eta })\\ \underset{j=1}{\overset{4}{\Sigma }}(\left(\underset{i=1}{\overset{5}{\Sigma }}{x}_i^{e1}\frac{{\partial \Psi }_i\left({\xi }_{\mathrm{gj}}\right)}{\partial \eta }\right)\frac{{\partial \Psi }_3\left({\xi }_{\mathrm{gj}}\right)}{\partial \xi }-\left(\underset{i=1}{\overset{5}{\Sigma }}{x}_i^{e1}\frac{{\partial \Psi }_i\left({\xi }_{\mathrm{gj}}\right)}{\partial \xi }\right)\frac{{\partial \Psi }_3\left({\xi }_{\mathrm{gj}}\right)}{\partial \eta })\end{array}\right)---\left(41\right)$

$=\left(\begin{array}{c}{y}_2^{e1}{g}_{23}^{e1}+y_4^{e1}{g}_{43}^{e1}+y_5^{e1}{g}_{53}^{e1}\\ x_2^{e1}{g}_{23}^{e1}+x_4^{e1}{g}_{43}^{e1}+x_5^{e1}{g}_{53}^{e1}\end{array}\right)$

方程(41)中的结果独立于节点的同等坐标和其在元素e1中的对角线 节点(diagonal node)。对于邻近元素e1和e2，

${\int }_{Q_{e1}}{b}_3^{e1}d\Omega +{\int }_{Q_{e2}}{b}_4^{e2}\mathrm{d\Omega }=\left(\begin{array}{c}{y}_2^{e1}{g}_{23}^{e1}+y_4^{e1}{g}_{43}^{e1}+y_5^{e1}{g}_{53}^{e1}\\ x_2^{e1}{g}_{23}^{e1}+x_4^{e1}{g}_{43}^{e1}+x_5^{e1}{g}_{53}^{e1}\end{array}\right)+\left(\begin{array}{c}{y}_1^{e2}{g}_{14}^{e2}+y_3^{e2}{g}_{34}^{e2}+y_5^{e2}{g}_{54}^{e2}\\ x_1^{e2}{g}_{14}^{e2}+x_3^{e2}{g}_{34}^{e2}+x_5^{e2}{g}_{54}^{e2}\end{array}\right)---\left(42\right)$

使用关系式和可将方程(42)写作

${\int }_{Q_{e1}+Q_{e2}}{b}_I^{e1}+b_I^{e2}\mathrm{d\Omega }=\left(\begin{array}{c}{y}_4^{e1}{g}_{43}^{e1}+y_5^{e1}{g}_{53}^{e1}\\ x_4^{e1}{g}_{43}^{e1}+x_5^{e1}{g}_{53}^{e1}\end{array}\right)+\left(\begin{array}{c}{y}_3^{e2}{g}_{34}^{e2}+y_5^{e2}{g}_{54}^{e2}\\ x_3^{e2}{g}_{34}^{e2}+x_5^{e2}{g}_{54}^{e2}\end{array}\right)---\left(43\right)$

应注意，来自图4。接着，对于多个元素nl，可容易地验证：

${\int }_{\Omega }{b}_I\mathrm{d\Omega }=\underset{n=1}{\overset{\mathrm{nl}}{\Sigma }}{\int }_{-1}^1{\int }_{-1}^1{b}_I^n\det \left(J_{\xi }^n\right)\mathrm{d\xi d\eta }=0---\left(44\right)$

无闭锁无网格扩展有限元公式

无发散无网格扩展有限元插值

从方程(9)中的近不可压缩问题可知，由基于标准位移的伽辽金法获得 的精确解u和唯一解u^h之间的误差可以能量误差的形式估计和表示如下：

${\left|\right|u-u^h\left|\right|}_E^2\le C_0\underset{v^h\in V^h}{\inf }{\left|\right|u-v^h\left|\right|}_E^2$

$\le C_0{\left|\right|u-{\Pi }_{h}u\left|\right|}_E^2$

$=C_{0}A(u-{\Pi }_{h}u,u-{\Pi }_{h}u)$(45)

$=C_0(2\mu {|u-{\Pi }_{h}u|}_{1,\Omega }^2+\lambda {\left|\right|▿·(u-{\Pi }_{h}u|)\left|\right|}_{0,\Omega }^2)$

方程(45)中的能量形式定义为

${\left|\right|u\left|\right|}_E={\left(A(u,u)\right)}^{\frac{1}{2}}---\left(46\right)$

H_h：H¹(Ω)→V^h(Ω)是与方程(21)中的等参映射F_e220相关的插值算子， 且定义如下：

在此，是映射算子，其保持多项式的次数≤1；也就是 且是上多项式的次数≤1的空间。||·||_m和|·|_m分别是 阶数m的索伯列夫范数或半范数。方程(45)中的常量C₀并不依赖于元素大 小h和被考虑的函数。应注意，采用众所周知的Céa引理在方程(45)中的第 一不等式来估计该误差。我们进一步假定解u∈H²(Ω)是规律的，且使用第一阶 近似的近似特性来获得全局近似误差

${\left|\right|u-u^h\left|\right|}_E^2\le C_0(2\mu {|u-{\Pi }_{h}u|}_{1,\Omega }^2+\lambda {\left|\right|▿·(u-{\Pi }_{h}u)\left|\right|}_{0,\Omega }^2)$

$\le C_{01}\mu h^2{\left|u\right|}_{2,\Omega }^2+C_{10}{\mathrm{\lambda h}}^2{|▿·u|}_{1,\Omega }^2---\left(48\right)$

$\le C_{11}{h}^2({\left|u\right|}_{2,\Omega }^2+\lambda {|▿·u|}_{1,\Omega }^2)$

从方程(48)可清楚地知晓，当网格在可压缩的情况下是精细的时，解 u^h将收敛于精确解u。然而，当λ→∞时，u^h的收敛可能不能通过低阶近似获 得。对于将要在不可压缩界点中受到约束的能量误差形式而言，需要使用高阶 近似严格执行无发散条件或是通过混合/简化积分/改进变量形式策略弱化执行 无发散条件并由上界-下界条件(inf-sup condition)予以保证。

与弱化执行无发散条件不同，可严格执行无发散条件，但是在高斯积分点 410(图4中示出为三角形)逐点执行该无发散条件，也就是

在高斯点，$▿·u^h\left({\xi }_{\mathrm{gi}}\right){|}_{Q_e}\to 0,\lambda \to \infty \forall Q_e\in M_h,i=\mathrm{1,2,3,4}---\left(49\right)$

发散的近似是在高斯点410估计的，其可表示为：

$▿·u^h\left({\xi }_{\mathrm{gi}}\right)=\mathrm{tr}(▿u^h\left({\xi }_{\mathrm{gi}}\right))=\mathrm{tr}(J_{\xi }^{-1}·{▿}_{(\xi ,\eta )}{u}^h\left({\xi }_{\mathrm{gi}}\right))=\mathrm{tr}({\left(G_F{G}_x^T\right)}^{-1}·{▿}_{(\xi ,\eta )}{u}^h\left({\xi }_{\mathrm{gi}}\right))---\left(50\right)$

在此，ξ_gi＝(ξ_gi，η_gi)，i＝1，2，3，4是高斯点410。方程(50)中的项可 另写成：

${▿}_{(\xi ,\eta )}{u}^h\left({\xi }_{\mathrm{gi}}\right)=\left(\begin{array}{ccc}\frac{{\partial \Psi }_1\left({\xi }_{\mathrm{gi}}\right)}{\partial \xi }& L& \frac{{\partial \Psi }_5\left({\xi }_{\mathrm{gi}}\right)}{\partial \xi }\\ \frac{{\partial \Psi }_1\left({\xi }_{\mathrm{gi}}\right)}{\partial \eta }& L& \frac{{\partial \Psi }_5\left({\xi }_{\mathrm{gi}}\right)}{\partial \eta }\end{array}\right)\left(\begin{array}{cc}{u}_1& v_1\\ M& M\\ u_5& v_5\end{array}\right)=G_F\left({\xi }_{\mathrm{gi}}\right)u_e,i=\mathrm{1,2,3,4}---\left(51\right)$

在此，u_e是元素Q_e204的节点位移的集合，并给出如下：

$u_e^T=\left(\begin{array}{ccc}{u}_1& \Lambda & u_5\\ v_1& \Lambda & v_5\end{array}\right)---\left(52\right)$

为了显示可通过使用ME-FEM近似在某些点上逐点地实现无发散条件， 使用特征值分析。使用5节点元素204Q_e＝[-1，1]×[-1，1]，该全局同等坐标与基 准同等坐标相符。换句话说，方程(26)中的雅可比矩阵变成且发散的表达式简并成简单形式。图5A示出了在矩形双线性4节点有限 元中的一个体积闭锁模式510，在此，全局同等坐标中的位移场给出如下：

[u，v]＝[c₁xy，0]    (53)

在方程(53)中，c₁是描述闭锁模式的变形的系数，其为节点位移的函数。 在评价全部4个高斯点410时，发散并不是无效的，并且在不可压缩界点 中，插值并不符合方程(49)的无发散条件。另一方面，可使用方程(50)表 达在高斯点410使用ME-FEM近似评价的处于相同变形模式的发散如下：

$▿·u^h\left({\xi }_{\mathrm{gi}}\right)=\mathrm{tr}(J_{\xi }^{-1}·{▿}_{(\xi ,\eta )}{u}^h\left({\xi }_{\mathrm{gi}}\right))$

对于i＝1，2，3，4，$=\mathrm{tr}(I·{▿}_{(\xi ,\eta )}{u}^h\left({\xi }_{\mathrm{gi}}\right))---\left(54\right)$

$=\underset{k=1}{\overset{5}{\Sigma }}{\Psi }_{k,\xi }({\xi }_{\mathrm{gi}},{\eta }_{\mathrm{gi}})u_k+\underset{k=1}{\overset{5}{\Sigma }}{\Psi }_{k,\eta }({\xi }_{\mathrm{gi}},{\eta }_{\mathrm{gi}})v_k$

根据方程(49)，ME-FEM元素Q_e204中的无发散插值为

$▿·u^h\left({\xi }_{\mathrm{gi}}\right)=0,i=\mathrm{1,2,3,4}---\left(55\right)$

或以矩阵形式

$\left(\begin{array}{cc}{\Psi }_{5,\xi }({\xi }_{g1},{\eta }_{g1})& {\Psi }_{5,\eta }({\xi }_{g1},{\eta }_{g1})\\ {\Psi }_{5,\xi }({\xi }_{g2},{\eta }_{g2})& {\Psi }_{5,\eta }({\xi }_{g2},{\eta }_{g2})\\ {\Psi }_{5,\xi }({\xi }_{g3},{\eta }_{g3})& {\Psi }_{5,\eta }({\xi }_{g3},{\eta }_{g3})\\ {\Psi }_{5,\xi }({\xi }_{g4},{\eta }_{g4})& {\Psi }_{5,\eta }({\xi }_{g4},{\eta }_{g4})\end{array}\right)\left(\begin{array}{c}{u}_5\\ v_5\end{array}\right)=-\left(\begin{array}{c}\underset{k=1}{\overset{4}{\Sigma }}({\Psi }_{k,\xi }({\xi }_{g1},{\eta }_{g1})u_k+{\Psi }_{k,\eta }({\xi }_{g1},{\eta }_{g1})v_k)\\ \underset{k=1}{\overset{4}{\Sigma }}({\Psi }_{k,\xi }({\xi }_{g2},{\eta }_{g2})u_k+{\Psi }_{k,\eta }({\xi }_{g2},{\eta }_{g2})v_k)\\ \underset{k=1}{\overset{4}{\Sigma }}({\Psi }_{k,\xi }({\xi }_{g3},{\eta }_{g3})u_k+{\Psi }_{k,\eta }({\xi }_{g3},{\eta }_{g3})v_k)\\ \underset{k=1}{\overset{4}{\Sigma }}({\Psi }_{k,\xi }({\xi }_{g4},{\eta }_{g4})u_k+{\Psi }_{k,\eta }({\xi }_{g4},{\eta }_{g4})v_k)\end{array}\right)$

上述线性系统由沿着在4个高斯点410和8个预定节点位移估计的已知形 状函数导数的2个未知值(u₅和v₅)和4个约束方程组成。然而，上述4个约 束方程并不是完全独立的。使用方程(79)，我们可将4个约束方程减少为两 个约束方程并通过下式获得2个未知值：

$\left(\begin{array}{c}{u}_5\\ v_5\end{array}\right)=\left(\begin{array}{c}\frac{c_1(-{\Psi }_{1,\xi }({\xi }_{g1},{\eta }_{g1})-{\Psi }_{3,\xi }({\xi }_{g1},{\eta }_{g1})+{\Psi }_{2,\xi }({\xi }_{g1},{\eta }_{g1})+{\Psi }_{4,\xi }({\xi }_{g1},{\eta }_{g1}))}{{\Psi }_{5,\eta }({\xi }_{g1},{\eta }_{g1})}]\\ 0\end{array}\right)---\left(57\right)$

上述结果证明ME-FEM元素中ME节点的扩展可精确产生在标准双线性 四边形元素中丢失的两个非闭锁模式。结果，全部的ME-FEM元素性能免于 体积闭锁。此外，该5节点ME-FEM元素的位移场在高斯点410是无发散且 逐点(divergence-free point-wise)。因为使用凸近似，方程(57)中的形状函 数导数Ψ_5，η(ξ_g1，η_g1)是非0的(也就是除了在ME节点自身以外，在元素中都是 非0的)，可唯一确定ME节点的位移。图5B介绍了使用ME-FEM法的特征 值分析结果。该模式形状520是免于体积闭锁的。靠近不可压缩界点的边界特 征值指示ME-FEM元是无体积闭锁的。来自特征值分析的结果可确保 ME-FEM元不包含0能量模式。因为方程(45)中的体积能量是被限制在不 可压缩界点内的，ME-FEM解的能量形式被认为是靠近O(h)的理论渐近速度 的。因此，使用下列式子计算压力：

$p^h=\lambda ▿·u^h---\left(58\right)$

其在不可压缩界点保持有限，且该L₂-范值误差估计保持

${\left|\right|▿·(u-u^h)\left|\right|}_{0,\Omega }\le c_{2}h{|▿·u|}_{1,\Omega }---\left(59\right)$

在某些情况下，虽然ME-FEM元素不能包含任何伪零能量模式，使用传 统高斯积分法则的ME-FEM公式在近不可压缩分析中经历压力振荡。使用传 统高斯积分法则的ME-FEM公式中的压力振荡的原因不同于中断压力/位移混 合有限元法中的压力振荡，该压力/位移混合有限元法中的压力振荡在集合的 压力方程中显示秩亏(rank-deficiency)。

ME-FEM法中的压力振荡是由于使用低阶近似的基于加权余差的伽辽金 法的固有特性。此外，ME-FEM法中的位移导数是间断越过元素边界的。采 用应变平滑过程来消除该振荡。可将某些后处理压力平滑技术用来产生平滑压 力分布。然而，还需要进行附加的计算，并且在非线性分析中无需保证鲁棒性 能。

为了抑制可能的压力振荡，开发了结合无发散ME-FEM插值的区域加权 的积分以为应变和压力以及无发散性能提供平滑效应。

区域加权的积分

将基于应变平滑概念的区域加权的积分引入到ME-FEM元素204以在应 变上提供平滑效应以抑制近不可压缩分析中的压力振荡。另外，区域加权的积 分在三角形和四面体形元素中提供无发散的ME-FEM插值。该应变平滑定义 如下：

在此，表示平滑的梯度算子。A_m601是平滑区域Ω_m602的区域。每个 元素604包含4个平滑区域602且每个平滑区域602具有来自图6A中示出的 两个高斯点610(示出为三角形)的应变贡献(strain contribution)。因此，第 一平滑区域302中的平滑应变表示为：

${\widetilde{ϵ}}^h=\frac{1}{2A_m}(ϵ({\xi }_{g1},{\eta }_{g1})\Delta \left(J_1\right)+ϵ({\xi }_{g2},{\eta }_{g2})\Delta \left(J_2\right))---\left(61\right)$

$A_m=\frac{\det \left(J_1\right)+\det \left(J_2\right)}{2}=\frac{\mathrm{\Delta t}{\left(G_F{G}_x^T\right)}_1+\Delta {\left(G_F{G}_x^T\right)}_2}{2}---\left(62\right)$

这可获得如下给出的一元模型中的平滑分散

这验证了不可压缩界点中的方程(49)的无发散条件。

对于图6B中示出的相邻元素614a-b，可从两个元素614a-b中的4个雅可 比行列式的和中计算出该平滑区A_m611如下：

$A_m=\frac{\Delta \left(J_2^{e1}\right)+\Delta \left(J_3^{e1}\right)+\Delta \left(J_1^{e2}\right)+\Delta \left(J_4^{e2}\right)}{2}$

$=\frac{\Delta {\left(G_F{G}_x^T\right)}_2^{e1}+\Delta {\left(G_F{G}_x^T\right)}_3^{e1}+\Delta {\left(G_F{G}_x^T\right)}_1^{e2}+\Delta {\left(G_F{G}_x^T\right)}_4^{e2}}{2}$(64)

在此，上标e₁和e₂为共享相同元素边m和两个连接节点I_s1和I_s2616a-b 的邻近元素614a-b。与方程(63)的推导构成类似，可证明邻近元素614a-b 中的平滑发散满足不可压缩界点中的无发散条件。使用方程(64)，可给 出如下平滑区域Ω_m中的应变-位移矩阵的加权平均值：

${\tilde{b}}_{\mathrm{Ix}}=\frac{{\left(b_{\mathrm{Ix}}\right)}_2^{e1}\Delta \left(J_2^{e1}\right)+{\left(b_{\mathrm{Ix}}\right)}_3^{e1}\Delta \left(J_3^{e1}\right)+{\left(b_{\mathrm{Ix}}\right)}_1^{e2}\Delta \left(J_1^{e2}\right)+{\left(b_{\mathrm{Ix}}\right)}_4^{e2}\left(J_4^{e2}\right)}{2A_m}---\left(66\right)$

${\tilde{b}}_{\mathrm{Iy}}=\frac{{\left(b_{\mathrm{Iy}}\right)}_2^{e1}\Delta \left(J_2^{e1}\right)+{\left(b_{\mathrm{Iy}}\right)}_3^{e1}\Delta \left(J_3^{e1}\right)+{\left(b_{\mathrm{Iy}}\right)}_1^{e2}\Delta \left(J_1^{e2}\right)+{\left(b_{\mathrm{yx}}\right)}_4^{e2}\left(J_4^{e2}\right)}{2A_m}---\left(67\right)$

在此，和中的分量表示为：

$\left(\begin{array}{c}{\left(b_{\mathrm{Ix}}\right)}_i^{e1}\\ {\left(b_{\mathrm{Ix}}\right)}_i^{e1}\end{array}\right)=J_i^{{e1}^{-1}}\left(\begin{array}{c}\frac{{\partial \Psi }_I^{e1}({\xi }_{\mathrm{gi}},{\eta }_{\mathrm{gi}})}{\partial \xi }\\ \frac{{\partial \Psi }_I^{e1}({\xi }_{\mathrm{gi}},{\eta }_{\mathrm{gi}})}{\partial \eta }\end{array}\right),$i＝2和3    (68)

$\left(\begin{array}{c}{\left(b_{\mathrm{Ix}}\right)}_i^{e2}\\ {\left(b_{\mathrm{Ix}}\right)}_i^{e2}\end{array}\right)=J_i^{{e2}^{-1}}\left(\begin{array}{c}\frac{{\partial \Psi }_I^{e2}({\xi }_{\mathrm{gi}},{\eta }_{\mathrm{gi}})}{\partial \xi }\\ \frac{{\partial \Psi }_I^{e2}({\xi }_{\mathrm{gi}},{\eta }_{\mathrm{gi}})}{\partial \eta }\end{array}\right),$i＝1和4    (69)

非常容易验证图4中示出的多个元素nl，内节点I的平滑梯度矩阵满足该 积分约束，也就是

${\int }_{\Omega }{\tilde{b}}_I\mathrm{d\Omega }=\underset{n=1}{\overset{\mathrm{nl}}{\Sigma }}{\int }_{-1}^1{\int }_{-1}^1{\tilde{b}}_I^n\Delta \left(J_{\xi }^n\right)\mathrm{d\xi d\eta }=\underset{n=1}{\overset{\mathrm{nl}}{\Sigma }}{\int }_{-1}^1{\int }_{-1}^1{b}_{I^{\prime }}^n\Delta \left(J_{\xi }^n\right)\mathrm{d\xi d\eta }=0---\left(70\right)$

因此，该区域加权积分保持线性精确性从而使得在可压缩情况下ME-FEM 公式有效，例如通过斑片试验(在有限元分析中第一个基本要求)。

改进的胡海昌-鹫津公式

为了在伽辽金近似中引入应变平滑公式，可考虑并给出胡海昌-鹫沣变形 原理的弱解如下：

在此，u是位移场，是假定应变场，且是使用假定应力场从方程(3) 中的本构方程计算的应力场。该项δW_ext表示质量力f和牵引力t作用到自然边 界Γ_N上的外部虚工。

通过假定压力场与位移梯度和平滑应变域的差分正交，采用假定的 应变法来消除方程(71)的RHS上的第二项。使用方程(60)，我们可以发现

其精确满足了该正交条件。方程(72)中的指数nm表示Ω(Ω＝∪_mΩ_m)中 平滑区域Ω_m的总数。

在使用方程(72)从方程(71)消除应力分量后，仅基于应变和位移场获 得改进的胡海昌-鹫沣公式的弱解：

${\mathrm{\delta U}}_{\mathrm{HW},\mathrm{mod}}(u,\widetilde{ϵ})={\int }_{\Omega }\delta {\widetilde{ϵ}}^{T}C\widetilde{ϵ}\mathrm{d\Omega }-\delta W_{\mathrm{ext}}---\left(73\right)$

每个元素中的位移由ME-FEM插值近似并可表示为：

$u^h=\underset{I=1}{\overset{5}{\Sigma }}{\Psi }_I{u}_I---\left(74\right)$

可近似并给出应变如下：

${\widetilde{ϵ}}^h=\underset{I=1}{\overset{\mathrm{NP}}{\Sigma }}{\tilde{B}}_I{u}_I---\left(75\right)$

在此，是方程(65)中定义的平滑梯度矩阵且5≤NP≤8是在每个平滑区 域Ω_m中相关节点数量。通过在方程(73)中引入位移和应变近似，获得下列 离散控制方程：

$f^{\mathrm{ext}}={\int }_{\Omega }{\Psi }^T\mathrm{fd\Omega }+{\int }_{{\Gamma }_N}{\Psi }^T\mathrm{td\Gamma }---\left(78\right)$

因为ME-FEM元素插值在元素边界保持克罗内克-δ特性，可以规范方式 (例如有限元分析)应用本质边界条件。

当使用包括应变平滑的ME-FEM法再次分析单个元素的特征值分析时， 可获得图5B中示出的基本同一的无闭锁模式。

同样地，图7中呈现了多个元素的结果。与图1B中示出的EFG结果相 比，具有应变平滑的ME-FEM的8个特征值被限制在近不可压缩界点中，且 对应的变形模式为非闭锁。

形状函数导数的对称特性

在5节点ME-FEM元中的形状函数导数的对称特性表示如下：

Ψ_1，ξ(ξ_g1，η_g1)＝-Ψ_2，ξ(ξ_g2，η_g2)＝-Ψ_3，ξ(ξ_g3，η_g3)＝Ψ_4，ξ(ξ_g4，η_g4)

Ψ_1，η(ξ_g1，η_g1)＝Ψ_2，η(ξ_g2，η_g2)＝-Ψ_3，η(ξ_g3，η_g3)＝-Ψ_4，η(ξ_g4，η_g4)

Ψ_2，ξ(ξ_g1，η_g1)＝-Ψ_1，ξ(ξ_g2，η_g2)＝-Ψ_4，ξ(ξ_g3，η_g3)＝Ψ_3，ξ(ξ_g4，η_g4)

Ψ_4，η(ξ_g1，η_g1)＝Ψ_3，η(ξ_g2，η_g2)＝-Ψ_2，η(ξ_g3，η_g3)＝-Ψ_1，η(ξ_g4，η_g4)

Ψ_2，η(ξ_g1，η_g1)＝Ψ_1，η(ξ_g2，η_g2)＝-Ψ_4，η(ξ_g3，η_g3)＝-Ψ_3，η(ξ_g4，η_g4)

Ψ_4，ξ(ξ_g1，η_g1)＝-Ψ_3，ξ(ξ_g2，η_g2)＝-Ψ_2，ξ(ξ_g3，η_g3)＝Ψ_1，ξ(ξ_g4，η_g4)

Ψ_3，ξ(ξ_g1，η_g1)＝-Ψ_4，ξ(ξ_g2，η_g2)＝-Ψ_1，ξ(ξ_g3，η_g3)＝Ψ_2，ξ(ξ_g4，η_g4)

                                                                                       (79)

Ψ_3，η(ξ_g1，η_g1)＝Ψ_4，η(ξ_g2，η_g2)＝-Ψ_1，η(ξ_g3，η_g3)＝-Ψ_2，η(ξ_g4，η_g4)

Ψ_5，ξ(ξ_g1，η_g1)＝-Ψ_5，ξ(ξ_g2，η_g2)＝-Ψ_5，ξ(ξ_g3，η_g3)＝Ψ_5，ξ(ξ_g4，η_g4)

Ψ_5，η(ξ_g1，η_g1)＝Ψ_5，η(ξ_g2，η_g2)＝-Ψ_5，η(ξ_g3，η_g3)＝-Ψ_5，η(ξ_g4，η_g4)

Ψ_1，ξ(ξ_g1，η_g1)＝Ψ_1，η(ξ_g1，η_g1)

Ψ_2，ξ(ξ_g1，η_g1)＝Ψ_4，η(ξ_g1，η_g1)

Ψ_5，ξ(ξ_g1，η_g1)＝Ψ_5，η(ξ_g1，η_g1)

在此，ξ_gi＝(ξ_gi，η_gi)i＝1，4为高斯点610。

参照图8，示出了使用根据本发明的一个实施例的无网格扩展有限元法在 压缩和/或近不可压缩区域中数值仿真结构性能的典型方法800的流程图。该 方法800可在计算机系统执行的软件或是应用模块中实现。

方法800始于步骤802，在此接收工程产品的FEM模型(例如，汽车、 飞机或组件)。该FEM模型包含低阶FEM元素，如双线性四边形元素或三角 形元素。该低阶元素仅包含角节点而不包含任何边中部边(mid-edge)节点或 中间节点。

接着，在步骤804，通过增加至少一个无网格扩展(ME)节点到每个FEM 元素的子集来创建无网格扩展FEM(ME-FEM)模型。该ME节点配置成使 用通用无网格近似来消除前述段落中描述的体积闭锁。该子集可包括整个或一 部分FEM模型。换句话说，ME-FEM和FEM元素的混合模型也可以用于该 仿真。

接着该方法800在步骤806使用区域加权积分技术执行应变平滑，以避免 潜在的压力振荡问题(数值问题)。应用该技术是为了确保在不使用后处理平 均技术的基础上仿真结果更加贴近现实。

最后，在步骤808，在ME-FEM模型创建以后，获得工程产品的数值仿 真结构性能。这可通过使用具有安装在计算机系统中的应用模块的ME-FEM 模型实现。该应用模块配置成执行无网格扩展有限元分析。该无网格扩展有限 元分析使用基于位移的第一阶凸无网格近似和元素素级(element wise)的插 值。该插值在每个ME-FEM元素的边界具有克罗内克-δ特性以允许应用本质 边界条件。该基于位移的第一阶凸无网格近似配置成消除体积闭锁。该数值仿 真结构性能可用于辅助用户(如科学家、工程师等)改进工程产品。

根据一个方面，本发明涉及一个或多个计算机系统，其能够执行在此描述 的功能。图9中示出了计算机系统900的例子。计算机系统900包括一个或多 个处理器，如处理器904。该处理器904连接到计算机系统内部通信总线902。 根据该示例性计算机系统来对各种软件实施例进行实施。在阅读该说明后，如 何使用其它计算机系统和/或计算机体系来实施本发明，对于本领域技术人员 而言是显而易见的。

计算机系统900还包括主内存908，优选随机存取存储器(random access  memory，RAM)，还可能包括辅助存储器910。该辅助存储器910可以包括， 例如一个或多个硬盘驱动器912和/或一个或多个可移动存储驱动器914，典型 的有软盘驱动器、磁带驱动器、光盘驱动器等等。可移动存储驱动器914以众 所周知的方式从可移动存储单元918读取数据和/或向其写入数据。可移动存 储单元918，典型的有软盘、磁带、光盘等，通过可移动存储驱动器914读取 和写入数据。应当理解，可移动存储单元918包括一计算机可用存储介质，其 中存储有计算机软件和/或数据。

在另一实施例中，辅助存储器910可包括允许将电脑程序或其它指令加载 到计算机系统900中的其它类似装置。这些装置可能包括，例如，可移动存储 单元922和接口920。这样的例子可能包括程序盒和盒式接口(如视频游戏设 备中存在的)、可移动内存芯片(如可擦除可编写只读式存储器(EPROM)、 通用串行总线架构(USB)快闪存储器或PROM)及相关插座，以及其它允许 将软件和数据从可移动存储单元922传送到计算机系统900的可移动存储单元 922和接口920。一般来说，计算机系统900是由操作系统软件来控制和协调 的，它执行诸如流程调度、内存管理、网络和I/O服务的任务。

也可能存在连接到总线902的通信接口924。通信接口924允许在计算机 系统900和外部设备之间传输软件和数据。通信接口924的实例可能包括调制 解调器、网络接口(如以太网卡)、通信端口、符合个人计算机存储卡国际协 会(Personal Computer Memory Card International Association，PCMCIA)标准 的插槽和卡等。通过通信接口924传输软件和数据。计算机900通过根据一套 特定的规则(即协议)，通过数据网络来与其它计算设备进行通信。常见的协 议之一是经常在因特网上使用的TCP/IP(传输控制协议/互联网协议)。一般 来说，通信接口924将数据文件打包成较小的数据包并通过数据网络传输该数 据包，或者对接收的数据包重新组装来获得原始数据文件。此外，通信接口 924处理每个数据包的地址部分，以便使其到达正确的目标端口，或者以拦截 最终目标为计算机900的数据包。在本文中，术语“计算机程序介质”、“计算 机可读介质”、“计算机记录介质”和“计算机可用的介质”被用来泛指诸如可 移动存储驱动器914(例如，闪存驱动器)的介质，和/或安装在硬盘驱动器 912中的硬盘。这些计算机程序产品是用于为计算机系统900提供软件的装置。 本发明涉及这样的计算机程序产品。

计算机系统900可能还包括输入/输出(I/O)接口930，用于来计算机系 统900提供显示器、键盘、鼠标、打印机、扫描仪和绘图仪等等。

计算机程序(也称为计算机控制逻辑)作为应用模块906存储在主存储器 908和/或辅助存储器910中。计算机程序也可通过通信接口924来接收。这种 计算机程序，在执行时，使计算机系统900能够执行本文所讨论的本发明的功 能。特别是，计算机程序，在执行时，使处理器904能够执行本发明的功能。 因此，这些计算机程序象征计算机系统900的控制器。

在通过软件实施本发明的实施例中，软件可以存储在计算机程序产品中， 并可以使用可移动存储驱动器914、硬盘驱动器912或通信接口924加载到计 算机系统900中。应用模块906，在被处理器904执行时，使处理器904执行 本发明所述的功能。

为完成想要做的任务，主存储器908可以载入有一个或多个应用模块906， 在无论有否通过输入/输出接口930输入使用者指令的情况下，此模块可以被 一个或多个处理器904执行。在运行时，当至少有一个处理器904执行应用模 块906其中一个时，计算出结果并存储在辅助存储器910(即硬盘驱动器912) 中。分析结果根据用户的指令通过I/O接口930以文本或图形演示的方式报 告给客户。

虽然本发明是根据特定实施例进行描述的，但这些实施例仅仅是对本发明 进行说明，而非限制的。对具体公开的示例性实施例进行的改动或改变会对本 领域的普通技术人员起到暗示和提醒作用。例如，虽然仅示出并描述了5节点 四边形元素(也就是具有四(4)角节点加一(1)无网格扩展节点的元素)以 解释ME-FEM数学工作的。其他类型的ME-FEM元素也可用于实现本发明的 基本相同的目的，例如6节点四边形元素(也就是具有四(4)角节点+二(2) 无网格扩展节点的元素)或四节点三角形元素。另外，虽然仅示出了并描述了 二维元素(四边形元素)来示范本发明是如何工作的，但是本发明也可在三维 元素(例如图2D中示出的四面体元素)中实现。用于三角形和四面体形元素 的ME-FEM与四边形元素的ME-FEM类似。总之，不能用所述的和特指的典 型实施例来限定本发明的保护范围，并且，所有可轻易地对本领域普通技术人 员能起到暗示和提醒作用的所有改动仍应涵盖在本专利申请的精神和所附加 的权利要求的范围之内。