基于面积坐标的四边形杂交有限元法

高何金雨 张世全

摘 要    ： 本文针对二维线弹性问题提出了一种基于面积坐标的新型杂交应力四边形有限元法AGQ- LQ  6. 该方法基于广义Hellinger-Reissner变分原理，位移逼近采用含内部位移的四节点广义协调元，应力逼近则采用九参数线性应力模式. 数值算例表明，本文构造的有限元既能保持面积坐标广义协调元对网格畸变不敏感及粗网格精度较高的优点，又能有效克服泊松闭锁现象.

关键词 ：线弹性问题； 四边形面积坐标方法； 杂交应力有限元； 泊松闭锁现象

中图分类号 ： O241.82 文献标识码 ：A DOI ：  10.19907/j.0490-6756.2023.041002

A hybrid stress quadrilateral finite element method based on area coordinates

GAO He-Jin-Yu，  ZHANG Shi-Quan

（School of Mathematics， Sichuan University， Chengdu 610064， China）

In this paper， we propose a new hybrid stress quadrilateral finite element method AGQ- LQ  6 for the two-dimensional linear elasticity problems by using the area coordinates. In this method， a 4-node generalized conforming element with internal displacements for the displacement approximation and a 9-parameter linear stress mode for the stress approximation are adopted based on the generalized Hellinger-Reissner variational principle. Numerical experiments show that the method is insensitive to mesh distortion. Meanwhile， it is also shown that this method can keep high precision even under coarse mesh and avoid the so-called Poisson-locking phenomenon effectively.

Linear elasticity problem; Quadrilateral area coordinates; Hybrid stress finite element; Poisson-locking phenomenon

（2010 MSC 65M60）

1 引 言

上世纪六十年代以来，如何针对平面弹性问题构造具有较高精度的低阶四边形有限元成为工程计算领域的一个研究热点. 其中，基于等参坐标的四边形元法一直有广泛的应用，其中最为典型的代表是4节点等参双线性元. 然而，在一些标准的考核算例中，这种最低阶的四边形位移元的数值精度较低，对网格畸变敏感，并且会发生泊松闭锁现象. 因此，在不增大计算规模的前提下，如何提高低阶四边形元的性能一直广受关注.

1973年，Wilson等  ［1］在等参双线性元基础上引入非协调内部位移，构造了计算精度更高的Wilson非协调元. 之后，Wilson等  ［2］又对Wilson元进行了改进，使其在泊松闭锁测试算例中也能算出一致收敛的结果. 此外，基于广义变分原理的杂交应力有限元框架也是改进四边形位移元的性能的有效途径， 其中的H-R变分原理  ［3，4］因为将位移和应力分别作为独立变量的特点而成为研究者们推导混合元/杂交元的一种普遍方法. 之后，在各类不同类型的微分方程有限元求解中，不同的改进方法，如质量集中杂交应力有限元法  ［5］，组合杂交有限元法  ［6］，弱Galerkin有限元法  ［7］，增强混合有限元法  ［8］等都有合理的运用，并且这些方法都能在不增加计算规模的情况下达到提高数值精度的效果.

1999年，Long等  ［9-11］提出了采用四边形面积坐标代替传统的等参坐标系来构造广义协调四边形单元的方法，称为“四边形面积坐标（QAC）方法”. 与之前的两类改进方法不同，QAC方法以坐标变换为切入点，替换掉常用的等参坐标系，转而使用一种新型面积坐标，其与笛卡尔坐标之间的变换关系是线性的，可相互显式表出.我们知道，对于等参双线性元而言，一方面，在等参变换中笛卡尔坐标可以由等参坐标的双线性形式表示出来，而等参坐标却无法用笛卡尔坐标的有限项表示出来; 另一方面，其单元刚度矩阵积分式中含有雅可比逆的行列式，只能通过数值积分近似计算，而QAC方法却可以克服上述缺点.

本文基于杂交应力有限元框架，采用Zhou和Xie在文献［12，13］中提出的应力模式对四边形面积坐标下的广义协调元进行改进，进而提出了一种基于面积坐标的杂交应力元AGQ- LQ  6. 在该方法中，单元的位移逼近我们采用含内部位移的四节点广义协调元AGQ6，应力逼近则采用九参数线性应力模式. 最后，我们以标准的数值算例验证了该方法的数值性能.

后文的安排如下. 在第2节中我们给出四边形面积坐标的定义及面积元的构造. 在第3节中，我们基于杂交应力有限元框架给出改进的面积元的离散格式. 我们在第4节中给出数值算例. 第5节总结了本文所得结果.

2 四边形面积坐标

2.1 四边形的特征参数

四边形单元面积坐标的构造基于三角形单元面积坐标构造的思想. 我们首先介绍四边形面积坐标需要用到的形状特征参数.如图1所示， A  是四边形单元的面积，1，2，3，4 为单元的四个顶点，以图中的顺序标号，  A′ 和  A″ 分别是三角形Δ124和Δ123的面积.参数的定义式为

g  1=  A′ A ， g  2=  A″ A ， g  3=1- g  1， g  4=1- g  2

（1）

由式（1） 可以看到，参数  g  3 和  g  4 能分别用  g  1， g  2 表示，且参数  g  1， g  2 在一个凸四边形中的取值范围是 0< g  1<1， 0< g  1<1 .

注1   我们称  g  1 和  g  2 为形状特征参数. 当  g  1， g  2 取一些特殊值时，图1中的单元会转变为一些有特定形状的单元：

（1） 当  g  1= g  2= 1 2   时，四边形单元就会退化为一个平行四边形；

（2） 当  （g  1 -g  2） （1-g  1 -g  2）=0 时，四边形单元就会退化为梯形；

（3） 当  g  1 g  2 （1-g  1）（ 1-g  2）=0 时，四边形单元就会退化为三角形.

2.2 四边形面积坐标的定义

如图2所示，令 P  为四边形单元内的任意一点，则 P 的面积坐标 （ L  1， L  2， L  3， L  4） 定义为  ［10］

L  i=  A  i A ，i=1，2，3，4  （2）

其中 A 是该四边形单元的面积，  A  i （ i =1，2，3，4）是由点 P 分别和单元其中两个相邻顶点构成的三角形的面积，单元四个顶点的坐标表示为：

节点1：  （ g  2，1- g  2，0，0）;

节点2：  （ 0，1-g  1， g  1，0）;

节点3：  （0，0，1- g  2， g  2） ;

节点4：  （1- g  1，0，0， g  1） .

四边形面积坐标也可以用笛卡尔坐标表示. 若图2中单元的四个顶点的笛卡尔坐标依次为 （ x  1， y  1） ， （ x  2， y  2） ， （ x  3， y  3） ， （ x  4， y  4） ，那么4个三角形的面积  A  i （ i =1，2，3，4）可以通过行列式计算得到

A  i= 1 2   1 x y 1  x  j  y  j 1  x  k  x  k   ，

i，j，k=1，2，3 i，j，k=2，3，4 i，j，k=3，4，1 i，j，k=4，1，2   .

再将  A  i 代入式（2）可得

L  i= 1 2A   a  i+ b  ix+ c  iy ，i=1，2，3，4   （3）

其中 A 表示该四角形单元的面积，   a  i， b  i ，c  i 以及 A 的表达式为

a  i= x  j y  k -x  k y  j ，

b  i= y  j -y  k ，

c  i= x  k -x  j ，

A=△124+△324=△132+△134 ，

i，j，k=1，2，3 i，j，k=2，3，4 i，j，k=3，4，1 i，j，k=4，1，2  .

进一步，四边形面积坐标还可以用等参坐标 （ξ，η） 表示. 如图3所示，笛卡尔坐标（ x ，  y ）和等参坐标 （ξ，η） 之间存在下述等参映射变换  F  K ：  K ︿ =  -1，1   2→K ：

x  y    =F   K（ξ，η）= 1 4 ∑  4  i=1 （1+ ξ   iξ）（1+ η   iη）   x  i  y  i  ，

其中 （ x  i， y  i） 表示在四边形四个顶点的笛卡尔坐标，

ξ  1  ξ  2  η  1  η  2    ξ  3  ξ  4  η  3  η  4   =   -1 1 -1 -1   1 -1 1 1     （4）

由式（3）可知，面积坐标  L  i（x，y） 是单元的局部坐标，通过等参变换改写为以 （ξ，η） 为变量的表达式  L  i（ξ，η） .我们将 x（ξ，η） 和 y（ξ，η） 代入式（3）并整理就可以得到

L  1= 1 4 （1-ξ）  g  2（1-η）+ 1- g  1 （1+η）  ，

L  2= 1 4 （1-η）［ （1-g  2） 1-ξ +

（1- g  1）（1+ξ）］ ，

L  3= 1 4  1+ξ   g  1（1-η）+ 1- g  2  1+η    ，

L  4= 1 4 （1+η）  g  1（1-ξ）+ g  2（1+ξ）

（5）

文献［11］中介绍了四边形单元，以及任意点 P 的面积坐标和笛卡尔坐标之间的一阶导数的变换，面积坐标沿边  L  i=0 （ i =1，2，3，4）的任意幂函数的线积分公式，以及面积坐标的任意幂函数的面积积分公式.

注2   与等参坐标相比，面积坐标下的幂函数能给出精确积. 从计算的角度看，相比于等参坐标单元刚度矩阵的数值积分精度更好，结果更准确.

3 基于广义协调元的杂交应力法

本节基于杂交应力有限元框架对面积坐标四边形元进行一定的改进.

3.1 记号与基本概念

考虑如下的线弹性问题：

-div  σ=f ， σ=Dε（u）  in Ω，

σ·n |    Γ   1= T -  ，  u |    Γ   0=0 ， on   Γ =  Ω =   Γ   1 - ∪   Γ   0 -    （6）

其中Ω  ∈R  2 是有界开集; u 代表位移; σ 是应力张量; ε（u）= 1 2 （SymbolQC@

+ SymbolQC@

T）u 表示应变; D 代表弹性模量矩阵; f 是规定的体力，即分布载荷，  T -  为给定的表面牵引力;    Γ   1 -  为Ω的边界   Ω中表面牵引力的部分. 在基于杂交变分原理的广义协调有限元推导中，位移场表示为 v= v  c+ v  I ， 其中  v  c 表示单元的协调位移，  v  I 表示非协调的内部位移， τ 表示假定的分片独立应力，位移和应力逼近的有限维空间   U  h= U  h  c U  h  I 以及  V  h 分别满足

U  h  c∈ U  c =  v∈  ∏  K∈ T  h  H  1 K    2；v |    Γ   0=0  ，

U  h  I   K  span{bubble functions}，

V   h∈V=∏  K∈ T   h H（ div ;K），

其中  T  h=∪{K} 表示Ω上四边形剖分的全体;   H  s（K），s≥0 是通常的整数阶Sobolev空间;  H（ div ;K） 的定义为

H（ div ;K）= τ=   τ  11， τ  22， τ  12   T∈ L  2  K  3   .

值得注意的是，文献［14，15］基于组合杂交变分原理引入了“能量协调性”的概念，对不相容内部位移   v  I 进行了研究，通过能量相容条件

∮    Kτ·n· v   I d s=0，  τ，  v   I   （7）

并基于修正的Hellinger-Reissner变分原理对面积坐标下的四边形有限元进行了改进.

3.2 九参数杂交应力格式

当我们采用含内部位移的四节点位移构造杂交应力有限元时，线弹性问题式（6）对应的Hellinger-Reissner变分原理为  ［12，14］：

Π   σ   h， u   h =  inf   v∈ U   h   sup   τ∈ V   h  Π （τ，v）   （8）

其中能量泛函Π （τ，v） 具体表示如下：

Π （τ，v）=∑  K ［- 1 2 ∫   Kτ· D  -1τ dΩ +

∫   Kτ·ε（v） dΩ -∮   Γ  1∩  K T - · v   c d s-

∮    Kτ·n· v   I d s-∫   Kf·v dΩ ］.

注意到

∫   Kτ·ε（v） dΩ =∫   Kτ·ε（ v   c） dΩ +  ∫   Kτ·ε（ v   I） dΩ ， v∈ U   h  I，

我们借助分部积分公式可得

∫   Kτ·ε（ v   I） dΩ -∮    Kτ·n· v   I d s=

-∫   K div τ· v   I dΩ .

泛函式（8）中的  Π  可以改写为

Π （τ，v）=∑  K ［- 1 2 ∫   Kτ· D  -1τ dΩ +

∫   Kτ·ε（ v   c） dΩ -∫   K div τ· v   I dΩ -

∮   Γ   1∩ K T - · v   c d s-∫   Kf·v dΩ ］.

采用含内部位移的四节点广义协调元AGQ6  ［9］的位移场，以及九参数线性应力模式，我们提出一种面积坐标下的杂交应力四边形有限元AGQ- LQ  6.   K∈ T  h ，用  U  h  c 表示面积坐标下相应的节点的位移空间，

U  h  c={ v  c； v  c |  K∈ span { L  1， L  2， L  3， L  4，

（ L  3- L  1），（ L  4- L  2）}  2} ，

即，   v  c∈ U  h  c ，具有如下形式

v  c |  K=   N  1 0  0  N  1   N  2 0  0  N  2   N  3 0  0  N  3   N  4 0  0  N  4   q  c= N  c q  c ，

其中节点位移

q  c=   u  1 v  1 u  2 v  2 u  3 v  3 u  4 v  4   T ，

形函数  N  i 为

N  o  i=-  g  k 2 + L  i+ L  j+ ξ  i η  i g  kP  （9）

i，j，k=1，2，3 i，j，k=2，3，4 i，j，k=3，4，1 i，j，k=4，1，2  .

其中

P= 3  L  3- L  1   L  4- L  2 -  L  3- L  1   g  2- g  3  1+ g  1 g  3+ g  2 g  4

-   g  1- g  2   L  4- L  2 - 1 2   g  2 g  4- g  1 g  3  1+ g  1 g  3+ g  2 g  4 .

用  U  h  I 表示面积坐标下相应的内部非协调位移的空间，

U  h  I={ v  I； v  I |  K∈ span    L  1 L  3， L  4 L  2   2，

K∈ T  h} ，

即，   v  I∈ U  h  I 有如下形式

v  I |  K=   N  λ1 0  0  Nλ  1   N  λ2 0  0  N  λ2   q  I= N  I q  I .

其中内部位移  q  I∈   R  4 ，形函数  N  λi 为

N  λ1= L  1 L  3 ，   N  λ2= L  2 L  4  （10）

那么整个位移空间可以记为  U  h  W= U  h  c  U  h  I .   v∈ U  h  W 具有如下形式

v |  K= （v  I+ v  c） |  K=［  N  cN  I］  q  c  q  I  .

注3   当单元形状退化为平行四边形时，有

g  1= g  2= 1 2 ， L  3- L  1= ξ 2 ， L  4- L  2= η 2  .

则形函数式（9）和式（10）会变为

N  o  i= 1 4  1+ ξ  iξ  1+ η  iη ， i=1，2，3，4，

N  λ1= 1 16 （1- ξ  2） ，

N  λ2= 1 16 （1- η  2） .

此时节点位移的形函数和内部位移的形函数与Q6元相同，在形式上只相差了常数倍.

当使用面积坐标时，我们用  V  h  1 表示相应的分片线性应力子空间，

V  h  1={τ；τ |  K∈ span {1，（ L  3- L  1），

（ L  4- L  2）}  3} ，

其中， τ∈ V  h  1 具有形式 τ |  K=  Φ   β   β  （τ） ，  β  （τ）∈   R  9 ，

Φ   β=  1     L  3- L  1 1    L  4- L  2  L  3- L  1 1    L  4- L  2  L  3- L  1     L  4- L  2   .

注4   面积坐标中只有两个是相互独立的坐标分量，所以在分片应力中， （ L  3- L  1）和（ L  4- L  2） 可以看作两个独立的分量.

AGQ- LQ  6基于下述修正的Hellinger-Reissner 变分原理：

Π    A-L  σ   h， u   h =  inf   v∈ U   h  w   sup   τ∈ V   h  1   Π    A-L（τ，v）  （11）

其中

Π    A-L（τ，v）=∑  K ［- 1 2 ∫   Kτ· D  -1τ dΩ +

∫   Kτ·ε（v） dΩ -∮   Γ   1∩ K T - · v   c d s-

∫   Kf· v   c dΩ ］.

因此，与式 （11） 相应的有限元格式为：

求   σ  h， u  h= u  c+ u  I ∈U  h  c  U  h  I  ，使得

a  σ  h，τ -b τ， u  h =0， τ∈ V  h  1 ，

b  σ  h， v  I+ v  c =f  v  c +g  v  c ，

v= v  I+ v  c∈ U  h  W  （12）

其中

α（σ，τ）=∑∫   Kσ· D  -1τ dΩ ，

b（τ，v）=∑∫   Kτ·ε（v） dΩ ，

f（v）=∑∫   Kf·v dΩ

g（v）=∑∮   Γ   1∩ K T - · v   c d s.

下面我们给出AGQ- LQ  6单元刚度矩阵的推导，单元应变表示为

ε（v）=ε（ v  c）+ε（ v  I）= B  c q  c+ B  c q  I=

［  B  cB  I］   q  c  q  I   ，

B  c =    B  c1    B  c2      B  c3    B  c4   ，

B  ci =     N  0  i  x  0 0    N  0  i  y     N  0  i  y     N  0  i  x   ， i=1，2，3，4，

B  I =    B  I1    B  I2   ，

B  Ii =     N  λi  x  0 0    N  λi  y     N  λi  y     N  λi  x    i=1，2 .

将式（12）中的积分项写为

a（σ，τ）   （K）=∫   Kσ· D  -1τ dΩ =

β   σ    T∫  1  -1  ∫  1  -1   Φ   β   T D  -1  Φ    β J  d ξ d η β   σ

β   σ    TH β   σ，

b（τ，v）   （K）=∫   Kτ·ε（v） dΩ =

β   σ    T∫  1  -1  ∫  1  -1   Φ    β   T［  B   cB   I］ J  d ξ d η   q   c  q   I

（ β   σ）   T［  G   cG   I］ β   σ，

f（ v   c）   （K）+ g（ v   c）   （K）=∫   Kf· v   c dΩ +

∮   Γ   1∩ K T - · v   c d s ［ Q   c0］   q   c  q   I  .

从式（12）中我们得到单元 K 上应力与位移的关系为：

β  σ= H  -1［  G  cG  I］   q  c  q  I   .

再代入到变分形式（12）的第二个等式中可得

G  T  c H  -1 G  c  G  T  c H  -1 G  I  G  T  I H  -1 G  c  G  T  I H  -1 G  I     q  c  q  I  =    Q  c  T 0    （13）

对内部位移参数采取与Wilson 元类似的静态冷凝方法，在式（13）中消去

q  I=-   G  T  I H  -1 G  I   -1 G  T  I H  -1 G  c q  c ，

则可将式（12） 转化为仅含有节点参数的常规单元刚度方程  K  8×8 q  c=  Q  c  T，  其中8×8的单元刚度矩阵具体表示为

K  8×8= G  T  c H  -1 G  c- G  T  c H  -1 G  I   G  T  I H  -1 G  I   -1

G  T  I H  -1 G  c .

单元刚度矩阵的显式表达式可以通过面积积分公式获得. 因为数值积分法对计算机编码更方便，所以我们也可以通过等参坐标求出单元刚度矩阵. 与 4 节点等参元的积分形式相比，在面积坐标下形函数中没有雅可比的逆，通过面积坐标积分公式便可以计算得到 ［K］ 的精确值.

3.3 收敛性分析

我们首先介绍网格的正则细分. 对于每一个四边形 K∈ T  h ，若  h  K， h′  K 和  θ  K  i 分别为单元 K 的直径，最小边长度和与顶点相关的角，则  T  h 是一个正则细分需满足以下条件  ［16］：

存在常数  γ′ 和 γ 使得下列不等式对所有 K∈ T  h 都一致成立：

h  K≤  γ′h′  K ，

max   1≤i≤4   cos  θ  K  i ≤γ≤1 .

杂交应力元AGQ- LQ  6的收敛需要满足如下更严格的细分条件：

在网格参数 h=  max   K  h  K  趋于0的情况下，单元 K 的两条对角线中点连线的距离  d  K 是  O（ h  2） 的.

当AGQ- LQ  6单元形状退化为平行四边形时，有下式成立：

g  1= g  2= 1 2 ， L  3- L  1= ξ 2 ， L  4- L  2= η 2 .

位移形函数式（9）和式（10）以及分片线性应力矩阵 Φ   β 写为

N  o  i= 1 4  1+ ξ  iξ  1+ η  iη ， i=1，2，3，4，

N  λ1= 1 16 （1- ξ  2） ，

N  λ2= 1 16  1- η  2 ，

Φ   β=  1 ξ 2  η 2     1 ξ 2  η 2     1 ξ 2  η 2    .

此时，AGQ- LQ  6元和文献［17］中的 LQ  6元的位移逼近空间以及应力逼近空间相同. 换句话说，根据 LQ  6元的收敛性分析，我们就能推导出AGQ- LQ  6元对于问题（6）有限元解的存在唯一性证明，并得到与 LQ  6元一样的误差估计，即在  L  2 范数意义下位移逼近能达到2阶精度，应变逼近能达到1阶精度，应力逼近能达到1阶精度.

当AGQ- LQ  6形状不是矩形的时候，单元的位移逼近我们采用的是广义协调元AGQ6的位移场. 针对广义协调元的收敛，Shi  ［17］通过FEM检验给出过证明. 因为非平行四边形情况下的杂交元AGQ- LQ  6较为复杂，本文没有给出严格的理论证明，但通过一系列工程算例来验证了单元的收敛性，并展示的数值算例结果中.结果显示，AGQ- LQ  6元的位移逼近可以达到2阶精度，应变逼近可以达到1阶精度，应力逼近可以达到1阶精度，符合预期的理论结果.

注5   在后面的数值算例中，我们均采用二分法的细分方法，即连接四边形对边中点，这种网格细分方法满足了杂交元AGQ- LQ  6的收敛条件.

4 数值算例

我们选取一些代表性数值考核算例，通过比较新方法与几个相关的四节点四边形元  ［9，12，18］的数值结果来验证了改进效果.

例4.1  （网格畸变测试） 在这个测试中，我们将悬臂梁剖分为两个四边形单元，如图4所示. 单元扭曲程度由扰动参数 e 来决定， 尖端点 A 的最大位移数据见表1，其中的杨氏模量 E =1500，泊松比 ν =0.25.

表1中的数据显示，四边形面积元 AGQ6不受单元扭曲程度影响，改进的面积元AGQ- LQ  6也能做到不受参数 e 的影响.

例4.2  （Cook平面应力板问题） 板的左端固定，右端施加均匀分布的单位载荷，杨氏模量 E  = 1，泊松比 ν = 1/3，见图5. 我们的目的是通过斜网格剖分来检验四边形单元的性能. 鉴于此模型没有解析解，我们提供的参考值为相对精细网格下的计算结果. 在 C 点的最大位移， A 点的最大主应力和 B 点的最小主应力的数值结果见表2.

表2中的数据显示，在各种网格情况下，改进的面积元AGQ- LQ  6的位移和应力结果比AGQ6元和 LQ  6元更接近参考值.

例4.3  （非多项式解析解的平面应变测试） 泊松locking-free测试的算例基于悬臂梁的模型来说明. 我们考虑更为复杂的初值条件和解析解为非多项式的情况. 考虑区域Ω=［0，10］×［-1，1］，杨氏模量 E  =1500，泊松比 初始值ν =0.3，如图6. 令体力为

f=E π  2   cos  πx  sin  πy  - sin （πx） cos （πy）   .

则问题的精确解为

u= u=（1+ν） cos （πx） sin （πy）  -2 （1-ν）  2xy v=- 1+ν  sin  πx  cos  πy   +  1-ν   2 x  2+ν 1+ν   y  2-1   ，

σ=E  -π sin  πx  sin  πy -2y 0 0  sin  πx  sin  πy   .

我们采用二分法进行网格细分. 表3～表5分别给出了位移的误差    u- u  h   0   u   0 ，应变的误差   ε（u）-ε  u  h    0   ε（u）   0  和应力的误差    σ- σ  h   0   σ   0  在不同网格下的数值结果. 值得注意的是， ECQ  4， LQ  6和AGQ6的位移和应力误差阶数结果从20×4的网格开始算.因在初始网格下的位移误差偏大，这对阶数的计算有所影响.

从表中数据，我们可以看到：在 ν →0.5时，改进的面积元AGQ- LQ  6的位移逼近精度能达到2阶，应变和应力逼近精度能达到1阶，且相比于 LQ  6和AGQ6的数值结果计算精度有一定的提高.特别地，粗网格下的数值结果比原单元数值性能明显提升.

5 结 论

本文针对二维线弹性问题提出了一种基于面积坐标的杂交应力四边形有限元AGQ- LQ  6. 数值算例显示，该方法既能保持面积坐标广义协调元对网格畸变的不敏感性和杂交应力元精度，又能有效克服泊松闭锁现象.

参考文献：

［1］   Wilson E L， Taylor R L， Doherty W P，  et al . Incompatible displacement modes ［M］. New York： Academic Press， 1973.

［2］  Wilson E L， Taylor R L， Beresford P J. A nonconforming element for stress analysis ［J］. Int J Numer Meth Eng， 1976， 10： 1211.

［3］  Hellinger E. Die allgemeinen ansatze der mechanik der kontinua， enz ［J］. Math Wis， 1914， 4： 602.

［4］  Reissner E. On a variational theorem in elasticity ［J］. J Math Phys， 1950， 29： 90.

［5］  陈豫眉， 周亚婧. 基于 Maxwell 模型的线性粘弹性波传播问题的质量集中杂交应力有限元法［J］. 四川大学学报： 自然科学版， 2021， 58： 061001.

［6］  付卫， 王皓， 张世全. 特征值问题的组合杂交有限元方法［J］. 四川大学学报： 自然科学版， 2017， 54： 708.

［7］  刘轩宇， 罗鲲， 王皓. 抛物型积分微分方程的新型全离散弱 Galerkin 有限元法［J］. 四川大学学报： 自然科学版， 2020， 57： 830.

［8］  张田田， 张百驹. 非线性弹性问题的增强混合有限元方法［J］. 四川大学学报： 自然科学版， 2021， 58： 011004.

［9］  Chen X M， Cen S， Long Y Q，  et al . Membrane elements insensitive to distortion using the uadrilateral area coordinate method ［J］. Comput Struct， 2004， 85： 35.

［10］  Long Y Q， Li J X， Long Z F，  et al . Area coordinates used in quadrilateral elements ［J］. Comput Mech， 1999， 19： 533.

［11］ Long Z F， Long Y Q， Li J X，  et al . Some basic formulae for area coordinates used in quadrilateral elements ［J］. Comput Mech， 1999， 19： 841.

［12］ Zhou T X， Xie X P. Optimization of stress modes by energy compatibility for 4-node hybrid quadrilaterals ［J］. Math Numer Sin， 2004， 59： 293.

［13］ Zhou T X. Stabilized hybrid finite element methods based on combination of saddle point principles of elasticity problem ［J］. Math Comp， 2003， 72： 1655.

［14］ Zhou T X. Finite element method based on combination of ‘saddle pointvariational formulations ［J］. Sci China E， 1997， 30： 285.

［15］ Xie X P. An accurate hybrid macro-element with linear displacements ［J］. Comm Numer Meth Eng， 2005， 21： 1.

［16］ Yu G Z， Xie X P， Carstensen C. Uniform convergence and a posteriori error estimation for assumed stress hybrid finite element methods ［J］. Comput Meth Appl Mech Eng， 2011， 200： 2421.

［17］ Shi Z C. The F-E-M-Test for convergence of nonconforming finite elements ［J］. Math Comput， 1987， 49： 391.

［18］ Sumihara K， Pian T H H. Rational approach for assumed stress finite elements ［J］. Numer Meth Eng， 1984， 20： 1685.