基于超级混合单元的复合材料层合板应力分析

刘艳红，张 莹

（中国民航大学航空工程学院，天津 300300）

传统混合有限元法（简称传统混合法）[1-5]的主要优点是：一方面，可同时求得应力和位移结果，而且应力结果的精度较高，一些算例表明其应力结果可接近精确解[1-2]；另一方面，不同单元同一结点的应力结果是连续的。但一般情况下，传统混合元系数矩阵的主对角线上有0元素。因此，不经特殊处理的传统混合元的数值结果不稳定。传统混合法中各种非常规的稳定元素技术方案所涉及的数学理论较复杂，不利于工程师掌握和应用[4]。另外，基于传统混合法的模型主要涉及具体的二维板壳类问题。有关三维工程问题的文献[5]表明，相应的稳定元素技术的理论过程更为复杂，因而很少有针对具体三维工程问题进行数值分析的文献。

无需任何稳定混合元素技术，基于两个变分原理（即最小势能原理和H-R变分原理）或广义混合变分原理可得到本质上稳定的广义混合元[6-7]，即系数矩阵主对角线上没有0元素。更进一步，文献[8]根据修正的H-R变分原理直接建立了广义的非协调辛元。

由于所采用的变分原理不同，广义混合元[6-8]的形式当然也不同。多个实例的数值分析表明：广义混合法[6-8]和被广泛应用的位移法或杂交应力法一样，用于分析线弹性静力学问题没有原则上的困难。不必讳言，广义混合法和传统混合法一样，在有限元控制方程中增加了应力变量。但正因如此，不需要从单元级依据材料的本构关系求解应力结果，可避免应力数值结果不准确的问题。非协调广义混合元的位移结果和应力结果收敛速度快。位移和应力变量的数值结果精确度高是文献[8]中非协调辛元最突出的特点。

超级单元的形式多种多样，构造方法种类繁多，如模态综合子结构法就是一种很常见的超级单元方法。在桁架或钢架结构分析中，人们常用“组合单元”这个名称来表述超级单元。一般情况下，超级单元法或组合单元法可统称为子结构法。钟万勰[9]提出的区段方法也可理解为子结构法。

研究以非协调辛元为基础，针对复合材料层合板的结构特点，给出了构建复合材料层合板超级混合单元的基本过程和步骤，并用实例证明其可靠性。

1 理论公式

1.1 修正的H-R变分原理与广义部分混合元列式

为了清楚地表示复合材料层合板超级部分混合元的基本理论，首先简要地给出非协调广义部分混合元的基本步骤。

式（2）中：σo=[σxzσyzσzz]T为平面外应力；Φ11、Φ21和Φ22是与材料参数相关的子矩阵；Δ

1、Δ

2和Δ

3是微分算子矩阵，即

平面内应力可由位移变量和平面外应力变量表示为

根据非协调位移元[13-14]理论，可将位移插值表达式分为协调部分和非协调部分，即

其中：diag[N]=[NeNeNe]，Ne为对应8结点六面体元等参形函数；位移参数向量qe=[uevewe]T；diag[Nr]=[NreNreNre]为对应单元内非协调位移参数向量re的形函数。

平面外应力σo变量可与位移u变量一样在相同的空间域内被离散，设其表达式为

平面外应力的结点参数向量pe=[σexzσeyzσezz]T。

根据文献[8]修正的H-R变分原理（1）离散泛函的变分结果，可得如下非协调广义部分混合元

1.2 广义部分混合元的辛对称性质

事实上，式（6）可改写为如下形式

为了方便，不妨设

根据辛代数的理论，设J为单位辛阵，即

式中：I24×24为一般意义下的单位阵，其下标表示维数。

不难证明以下关系

从式（7）可看出，其系数矩阵本身是不对称的，而在辛代数中，其系数矩阵辛对称。

1.3 超级混合元

假设层合板结构如图1（a）所示，其第i层中结点的排列如图1（b）所示，则可方便地分离式（7）中上下表面的结点，因而得到如下表达式

图1 层合板结构和第i层的结点排列Fig.1 Laminated plate structure and joint array on Level i

对式（11）交换第2行（列）和第3行（列）有

上式可简写为

式中

对式（13）求和，可得第i层的控制方程为

假设考虑如图1（a）所示的4层结构，根据层合板中各层间上下表面的平面外应力和位移的关系（即层间上下表面的平面外应力和位移相等的关系），对式（14）求和，可得总体结构的控制方程为

根据文献[6-8]中的求解方法，自然也可对控制方程（15）同时引入平面外应力和3个位移变量的边界条件。求得位移和平面外应力结果后，平面内应力的求解方法与文献[6-8]相同。

2 数值算例

考虑如图 2（a）所示的 2 层复合材料板[0°/90°]，已知板的长和宽a=b=1.0，厚度h=2h1=0.1。材料参数[15]为：E11=10E22=10E33，G12=G13=0.6E33，G23=0.5E33，μ12=μ13=μ23=0.25。板上表面（z=0.1处）受垂直向下的均布载荷q=0.1作用。四边简支（SSSS）：在x=0和x=a处，σx=w=v=0；在y=0和y=b处，σy=w=u=0。

图2 复合材料板、坐标系与层数划分Fig.2 Composite laminates,coordinate system and layer division

厚度方向的层数划分如图2（b）所示。上下2层中各4薄层，每层厚度s1=0.01125h；接近中性面有上下2个薄层，每层厚度s2=0.005h。所以，总体上板被划分为10层。考虑到结构的对称性，有限元的数值分析采用1/4模型。网格划分为12×12×10，分别采用了广义部分混合元及超级混合元方法计算应力及位移，两种算法所消耗的时间如表1所示。

表1 两种方法对比分析Tab.1 Comparative analysis between two methods

因为每层4个薄层的单元及材料相同，另一薄层的单元及材料也相同，因此，上下2层的4个薄层分别只需分析1层，并考虑接近中性面的2个薄层，所以超级混合单元法仅对576个单元进行数值积分。就此实例而言，从表1中可看出，超级混合单元法大约可节省40%的时间。可见采用超级混合元分析层合结构静力学问题的成效明显。

图3～图8为SME（代表超级混合元）中位移变量和部分应力变量与精确解（EXACT）[15]的对比结果，其中最大误差发生于σxz（a/2，b/8，z）沿厚度的分布这一变量上，如图5所示，不超过0.8%。由于材料不同，图7表明平面内应力σxx在2层材料的界面之间不连续，由此可见所研究的复合材料层合板的超级混合元同样能反映这一客观特性。图3～图8中，x3指z方向，u1指位移 u，u3指位移 w，数字 1、2、3 分别指 x、y、z方向。

图3 位移u（a/8，b/2，z）沿厚度方向的分布Fig.3 Distribution of displacement u（a/8，b/2，z）along thickness direction

图4 位移w（a/2，b/2，z）沿厚度方向的分布Fig.4 Distribution of displacement w（a/2，b/2，z）along thickness direction

图5 应力σxz（a/2，b/8，z）沿厚度方向的分布Fig.5Distribution of stress σxz（a/2，b/8，z）along thickness direction

3 结语

图6 应力σzz（a/2，b/8，z）沿厚度方向的分布Fig.6Distribution of stress σzz（a/2，b/8，z）along thickness direction

图7 应力σxx（a/2，b/2，z）沿厚度方向的分布Fig.7Distribution of stress σxx（a/2，b/2，z）along thickness direction

图8 应力σxy（a/8，b/8，z）沿厚度方向的分布Fig.8Distribution of stress σxy（a/8，b/8，z）along thickness direction

考虑到复合材料层合板的特点，提出了基于非协调辛元的超级混合单元。超级混合单元的组装仅涉及矩阵的加法，因而运算快捷。算例表明，用超级混合元分析层合结构静力学问题的效果明显。但必需说明：对于不同的问题，节省时间不仅与结构本身的形式相关，还与所划分的相同层数多少有关。对于复杂结构，可能没有简单结构的效果明显。但对于大型的复杂结构而言，计算模型所用时间较长，只要网格划分得当，速度提高的效果也可能非常有意义。

如果要获得压电材料类复杂层合结构静力问题的高精度结果，那么广义部分混合法是可选的方法之一，而所提出的超级混合单元法同样可用于这类问题的数值分析。