引 言
对于通常的位移有限元法而言，为了提高收敛速度和(或)克服各类自锁现象，前人提出了许多不同的处理方法，其中最具有代表性的一类方法是假设应变、应力或非协调位移中的一个或多个独立的场量．一般来说，假设的未知场量在单元内，从而可以在单元级上消除这些未知量，最终形成单元刚度矩阵．以上可参考基于H-R变分原理的应力(应变)假设[1]、假设自然应变(ANS)理论[2]、增强假设应变(EAS)理论[3]和非协调位移法[4]等文献．
通常的位移法需要利用本构关系来计算应力．由于应变-位移是微商关系，微商只能转化为差商进行近似计算，这是应力结果精度不如位移结果精度的主要原因之一．所以为了提高应力结果的精度，需要额外的应力恢复技术．
经典混合单元[5-7]的有限元模型可同时求解位移和应力，两类变量的结果自然连续，并且收敛速度快、结果精度高．但是，基于经典混合单元控制方程是非正定的(主对角线上存在零元素)，因而若直接采用经典混合单元求解具体问题，需要借助非常规的稳定技术对单元列式进行改造，或采用特殊的求解技术，其数学理论比较复杂[6-7]．最近，为了改善这个问题，在不需要任何单元稳定技术方案的前提下，文献[8]提出了一种三维的8结点非协调广义混合单元，位移和应力变量均为有限元模型的线性方程组的未知量，可同时引入位移和应力变量的边界条件求解，理论简单实用．该单元最突出的特点是避免了传统混合法中系数矩阵主对角线上有零元素的问题，在不要求额外的稳定混合元方案的前提下，成功地解决了传统混合法中位移和应力结果震荡的问题．实践表明这种单元的位移和应力结果的收敛性好、精度非常高．但是，关于该单元非协调部分的积分比较复杂，需要较多的计算时间．
在假设自然应变方法中，应变场直接与某些点处位移相关的应变相关联，因此不会出现独立的自由度．文献[9]提出了一种应用于4结点板、壳单元以及8结点实体单元的增强假设应变理论，讨论了其与H-R变分原理的一般等价性，并通过数值实验证明了该方法的优异性能．
本文在不牺牲非协调广义混合单元优点的情况下，结合文献[9]的增强假设应变理论，建立一种更为高效的非协调广义混合单元，并加以实例验证．
1 基 本 理 论
1.1 最小势能原理和H-R变分原理
假设位移边界条件事先满足
且不考虑体积力，势能原理可表示为[10]
u)TC(
(1)
式中,表示求微分；u=[u1 u2 u3]T为位移向量；C为材料刚度矩阵；V为结构的体积；S为结构的表面积；
为作用在表面Sσ上已知的应力．
H-R变分原理可表达为[11]
(2)
式中，σ=[σ13 σ23 σ33 σ11 σ22 σ12]T是应力向量；S是材料柔度矩阵．
1.2 协调部分插值函数
考虑线弹性力学问题中的8结点协调单元，设位移场u和应力场σ可以用相同的形函数表示：
u=Nqe，
(3)
σ=Npe，
(4)
其中，N为形函数矩阵，具体表达式可参考文献[8]；qe为单元结点的位移向量；pe为单元结点的应力向量．
1.3 改进的非协调位移模型
对于8结点实体单元，采用增强假设应变理论[9]替代文献[8]中的非协调部分插值函数：
(5)
其中，re是单元内部应变参数向量；M是附加应变场的插值函数矩阵，
(6)
式中，J为Jacobi矩阵；J0为ξ=η=ζ=0时的Jacobi矩阵；
为6×6的矩阵；Mξ为形函数插值矩阵；Mξ和T0的具体形式可参考文献[9]．
2 改进的非协调广义混合单元模型
2.1 改进的非协调位移模型
根据文献[9]，对于非协调单元，应变ε=u可以表达为协调部分(N)qe与非协调部分Mre之和：
ε=u=(N)qe+Mre．
(7)
与文献[8]相同，将式(7)代入式(1)，可得
(8)
其中
N)TC(
N)TCMdV,
值得说明的是：依据非协调位移法的处理方法[12]，若不考虑体积力，式(8)中的
项可以被忽略．
由δΠP(re)=0,可得
(9)
将式(9)代入式(8)消去re，并对qe进行变分，可以得到Euler方程：
Rqqqe=fq,
(10)
其中
同样地，将式(7)和式(4)代入式(2)，可得
(11)
其中
将式(9)代入式(11)消去re，可得
(12)
其中
将由δΠHR(pe,qe)=0得到的Euler方程和式(10)相加联立，可得
(13)
对全部单元求和，得到非协调广义混合单元的整体刚度矩阵：
(14)
其中
R11=∑Rpp， R12=∑Rpq， R22=∑Rqq， f=∑fp．
2.2 位移和应力边界条件的引入
与文献[8]引入边界条件的方法相同，假定α为表面结点已知的位移和应力向量，其值由给定的作用在表面Sσ上的应力
和表面Su上的位移
决定．交换式(14)的行列并进行化简，可得改进的广义混合有限元法的最终求解方案为
(15)
其中，
和
分别为未知的结点应力和位移参数向量．
3 改进方法的优越性
文献[8]中位移模型中非协调部分的形式为[4,12]
u=Nqe+Nrre,
(16)
其中，Nr为形函数矩阵，具体形式可参考文献[13]．
在形成单元刚度矩阵的过程中需要对每个单元的Nr进行8次积分操作，然后再组装形成最终的总体刚度矩阵，这是繁琐且消耗时间的计算过程．在目前的方法中，非协调部分为式(5)的形式，仅需要对Mξ进行1次积分操作便可以得到单元刚度矩阵，简化了计算过程，并且从理论上讲可提高计算效率．
4 数 值 算 例
4.1 经典梁问题分析
考虑文献[14]中的悬臂梁，如图1、2所示．几何尺寸a=10.0，b=h=2.0，无量纲的材料参数E=1 500，μ=0.25．右端分别施加纯弯矩和剪切载荷．5种不同的网格如图3～7所示．
图1 纯弯矩载荷(载荷1)图2 剪切载荷(载荷2)
Fig. 1 The pure bending load (load case 1)Fig. 2 The shear force load (load case 2)
图3 网格a图4 网格b
Fig. 3 Mesh aFig. 4 Mesh b
图5 网格c图6 网格d
Fig. 5 Mesh cFig. 6 Mesh d
图7 网格e
Fig. 7 Mesh e
表1～4分别列出了两种载荷作用下，点A竖直位移u3和点B平面内应力σ11的值，其中，QS11-1，QS11-2及精确解的值均摘自文献[14]．
表1 点A的位移u3(载荷1)
Table 1 Displacements u3 at point A(load case 1)
meshQS11-1QS11-2presentexacta100.0100.0100.0100.0b89.092.296.9100.0c78.881.491.6100.0d77.077.085.7100.0e92.292.998.3100.0
表2 点B的应力σ11(载荷1)
Table 2 Stresses σ11 at point B(load case 1)
meshQS11-1QS11-2presentexacta-3 000.0-3 000.0-3 000.0-3 000.0b-2 619.5-2 710.8-2 853.7-3 000.0c-2 327.2-2 413.1-2 635.2-3 000.0d-2 348.0-2 348.0-2 654.0-3 000.0e-3 006.6-3 015.1-2 998.1-3 000.0
表3 点A的位移u3(载荷2)
Table 3 Displacements u3 at point A(load case 2)
meshQS11-1QS11-2presentexacta96.195.9101.0102.6b86.589.598.4102.6c78.781.594.5102.6d78.578.489.0102.6e94.194.9102.3102.6
表4 点B的应力σ11(载荷2)
Table 4 Stresses σ11 at point B(load case 2)
meshQS11-1QS11-2presentexacta-3 375.0-3 375.0-3 205.8-3 375.0b-2 925.6-3 024.8-2 979.1-3 262.5c-2 847.6-2 950.0-2 869.1-3 375.0d-2 735.0-2 735.7-2 821.0-3 150.0e-4 125.3-4 138.2-4 108.1-4 050.0
从表1～4的数据可以看出，本文位移和应力结果均比杂交应力单元QS11-1和QS11-2的结果精确．同时，改进的广义混合有限元方法对于模型网格的几何形状变化的敏感度相比于杂交应力单元更低，即使是不规则的网格尺寸，本文的方法也可得到精确度更高的结果．
4.2 收敛性及计算时间要求分析
考虑厚度均匀的平板，其几何尺寸a=b=1.0， h=0.10，无量纲的材料参数E=210， μ=0.3．边界条件：x1=0，a时，σ11=u2=u3=0；x2=0，b时，σ22=u1=u3=0．上表面施加垂直向上的载荷，σ0=1.0，下表面无载荷．4个模型的网格尺寸分别为4×4×6，8×8×6，12×12×6，16×16×6．
图8 u1(0,b/2,0)收敛速度图9 u3(a/2,b/2,h/2)收敛速度
Fig. 8 Comparison of u1(0,b/2,0)Fig. 9 Comparison of u3(a/2,b/2,h/2)
图10 σ13(0,b/2,h/2)收敛速度图11 σ33(a/2,b/2,h/2)收敛速度
Fig. 10 Comparison of σ13(0,b/2,h/2)Fig. 11 Comparison of σ33(a/2,b/2,h/2)
图12 σ11(a/2,b/2,h)收敛速度图13 σ12(0,0,0)收敛速度
Fig. 12 Comparison of σ11(a/2,b/2,h)Fig. 13 Comparison of σ12(0,0,0)
式(17)给出了不同网格尺寸的误差：
(17)
其中，Selement为该网格尺寸下的有限元解．分别以单元平面方向尺寸的对数和误差为横纵坐标作图，图8～13表示在4种有限元网格下各个变量的收敛情况，其中缩写符号NC3D8表示ABAQUS中8结点非协调位移单元，精确解Sexact由文献[15]中的方法得出．
图8～13表明，目前方法得到的结果误差及收敛情况均优于非协调的位移单元，且与文献[8]相差不大．
表5表示网格为12×12×6时本文方法与文献[8]方法计算时间的对比．当网格尺寸相同时，本文方法的计算时间小于文献[8]的方法，证明了本文方法可以提高计算效率．
表5 计算效率
Table 5 The computational efficient
methodref. [8]presenttime t/s2.881.98
4.3 复合材料层合板分析
考虑厚度均匀的复合材料层合板[0°/90°]，几何尺寸a=b=1.0，h=0.10，无量纲的材料参数E11=25E22=25E33，G12=G13=0.5E33，G23=0.2E33，μ12=μ13=μ23=0.25．边界条件为一对边固支一对边简支：x1=0，a时，u1=u2=u3=0；x2=0，b时，σ22=u1=u3=0．假设上表面施加垂直向上的载荷p=1.0sin(πx/a)sin(πy/b)，下表面的3个面外应力分量为0．4个模型的网格尺寸分别为6×6×12，8×8×12，10×10×12，12×12×12，精确解摘自文献[16]．
表6 复合材料层合板位移和应力结果
Table 6 Displacements and stresses of a composite laminate
variable6×6×128×8×1210×10×1212×12×12exact[16]εerror/%u1(a/4,b/2, h)-0.978-0.983-0.977-0.977-0.9760.10u3(a/2, b/2, h/2)0.6420.6460.6470.6480.6490.15σ13(a/8, b/2, h/2)1.6441.7751.6551.5821.5920.63σ33(a/2, b/2, h/2)0.6440.6450.6450.6450.640.78σ11(a/2, b/2, 0)-4.783-4.734-4.720-4.714-4.6531.31σ12(a/8, 0, 0)0.2110.2160.2160.2180.2211.36
表6为复合材料层合板不同网格下各个物理参数的数值结果，其中εerror是网格为12×12×12时的百分误差，由式(17)给出．表6数据表明，改进的非协调广义混合单元的结果贴近精确解，误差都在合理范围之内，并且当网格尺寸比较粗糙时，也可以得到可接受的结果．
5 结 论
文献[8]中提出的广义混合单元的非协调部分的积分操作比较复杂，消耗大量计算时间．本文结合增强假设应变理论，建立了新的8结点非协调广义混合单元，非协调部分只需要一次数值积分．因此，本文的非协调广义混合单元大大地节省了计算时间，在资源要求方面优于文献[8]．数值实例表明，改进的非协调广义混合单元结果稳定，精确度高且计算速度快．
基于与板壳理论相关的广义变分原理，同样可以构建高效的广义混合壳元．
[1] LEE S W, RHIU J J. A new efficient approach to the formulation of mixed finite element models for structural analysis[J]. International Journal for Numerical Methods in Engineering, 1986, 23(9): 1629-1641.
[2] DVORKIN E N, BATHE K J. A continuum mechanics based four-node shell element for general non-linear analysis[J]. Engineering Computations, 1984, 1(1): 77-88.
[3] SIMO J C, RIFAI M S. A class of mixed assumed strain methods and the method of incompatible modes[J]. International Journal for Numerical Methods in Engineering, 1990, 29(8): 1595-1638.
[4] TAYLOR R L, BERESFORD P J, WILSON E L. A non-conforming element for stress analysis[J]. International Journal for Numerical Methods in Engineering, 1976, 10(6): 1211-1219.
[5] BREZZI F, FORTIN M. Mixed and Hybrid Finite Element Methods[M]. New York: Springer-Verlag, 2011.
[6] QIU W, DEMKOWICZ L. Variable order mixed h-finite element method for linear elasticity with weakly imposed symmetry II: affine and curvilinear elements in 2D[J]. Mathematics, 2010. DOI: arXiv:1011.0485.
[7] GOPALAKRISHNAN J, GUZM
N J. Symmetric nonconforming mixed finite elements for linear elasticity[J]. SIAM Journal on Numerical Analysis, 2011, 49(4): 1504-1520.
[8] QING G, MAO J, LIU Y. Generalized mixed finite element method for 3D elasticity problems[J]. Acta Mechanica Sinica, 2018, 34(2): 371-380.
[9] ANDELFINGER U, RAMM E. EAS-elements for two-dimensional, three-dimensional, plate and shell structures and their equivalence to HR-elements[J]. International Journal for Numerical Methods in Engineering, 1993, 36(8): 1311-1337.
[10] 田宗漱, 卞学鐄. 多变量变分原理与多变量有限元方法[M]. 北京: 科学出版社, 2011.(TIAN Zongshu， PIAN T H H. Variational Principle With Multi-Variables and Finite Element With Multi-Variables[M]. Beijing: Science Press, 2011.(in Chinese))
[11] REISSNER E. On a variational theorem in elasticity[J]. Studies in Applied Mathematics, 1950, 29(1/4): 90-95.
[12] LIOU W J, SUN C T. A three-dimensional hybrid stress isoparametric element for the analysis of laminated composite plates[J]. Computers & Structures, 1987, 25(2): 241-249.
[13] 陈万吉. 一个高精度八结点六面体单元[J]. 力学学报, 1982(3): 262-271.(CHEN Wanji. A high-precision eight-node hexahedron element[J]. Chinese Journal of Theoretical and Applied Mechanics, 1982(3): 262-271.(in Chinese))
[14] CHEUNG Y K, WANJI C. Isoparametric hybrid hexahedral elements for three dimensional stress analysis[J]. International Journal for Numerical Methods in Engineering, 1988, 26(3): 677-693.
[15] PAGANO N J. Exact solutions for rectangular bidirectional composites and sandwich plates[J]. Journal of Composite Materials, 1970, 4(1): 20-34.
[16] VEL S S, BATRA R C. Analytical solution for rectangular thick laminated plates subjected to arbitrary boundary conditions[J]. AIAA Journal, 2012, 37(11): 1464-1473.