基于等几何分析的薄壳静力学分析

刘慧善，朱灯林

(河海大学机电工程学院，江苏 常州 213022)

HUGHES等[1]提出的等几何分析方法(isogeometric analysis，IGA)以CAD模型的NURBS/B-Spline[2-3]函数作为基函数来构造力学分析中的位移场，对复杂模型进行精确建模，可以采用稀疏单元划分获得较高的分析精度，提高了运算效率，同时便于CAD/CAE集成，特别适用于复杂结构和板壳结构的分析和优化设计。

板壳结构广泛应用于各种工程实际中，如机械结构、建筑结构等。对板壳结构的力学性能分析可采用三维实体单元式壳单元进行计算。HUGHES[4]提出了利用高阶多项式对板壳单元插补的方法，进行板壳结构的数学建模。但是高阶多项式插补存在着效率低、建模能力差等问题。针对复杂板壳结构的等几何分析[5-6]，DORNISCH[7]推导了等几何分析过程中坐标系的建立方法以及坐标系的相互转换。KANG[8]提出了裁剪单元方法，对复杂曲面模型进行了分析。KIENDL[9]提出了等几何Kirchhoff-love壳单元，并用弯曲片方法处理多片组合壳体问题。詹双喜[10]根据Kirchhoff-love壳理论推导了薄壳自由振动分析等几何方法的列式。BENSOND[11]结合Kirchhoff-love和Reissner-Mindlin两种壳单元特点发展了混合单元。

本文对Reissner-Mindlin(RM)板壳结构的等几何分析建模和位移的插补方法进行了研究，分别采用了两种位移插补方法：一种是根据HUGHES[4]提出的高阶多项式插补方法推导的RM壳单元IGA位移插补方法，另一种是根据KANG[8]提出的位移插补方法推导的RM壳单元IGA位移插补方法。运用这两种方法对相关实例进行了分析计算，比较了两种插补方法的精确度和计算效率。研究结果可为板壳结构的等几何分析研究和应用提供参考。

1 基于IGA的壳体几何模型构造

1.1 壳体的定义

壳体结构通常采用其参考面和厚度方向的参数对其进行几何定义，如图1所示，其数学模型如下：

(1)

图1 壳体定义

1.2 RM壳体坐标系

其中

式中x,ξ,x,η表示对参数坐标系求偏导。

图2 RM壳局部坐标系(lamina[1]坐标系)

由此可得从全局坐标系到壳局部坐标系的转换矩阵：

式中：i,j分别代表矩阵的行和列。

在传统壳结构分析[14-16]中，壳内任一点的位移可由参考面沿总体坐标系x,y,z方向的3个位移分量ua,va,wa及对应的θa1,θa2绕节点坐标系转动所确定。由于传统壳单元节点位于参考面上，因此各节点的局部坐标系按文献[1]中的方法确定即可。但是在IGA中，由于控制节点通常不在参考平面内，这就给节点处的坐标系的建立以及后续位移场的插补带来了困难。为了在每个控制节点处建立相应的局部坐标系，本文采用如下的插补方法。

(2)

(3)

式中：np为参考面上控制节点数；Ω为参数空间积分区域。

图3 IGA控制节点坐标系(fiber[1]坐标系)

则有：

(4)

根据式(4)可得：

(5)

2 壳体位移插补方法

(6)

(7)

(8)

第二种是根据KANG在文献[8]中提出的位移插值方法推导得到的RM壳单元IGA位移插补方法。在插值过程中直接采用角度进行插值，然后再将转角转换为各个方向的位移。具体表达式如下：

(9)

3 刚度矩阵的求解

采用第一种位移插补方法时，应变位移函数矩阵定义如下:

B=[B1B2…Ben]

(10)

应变可以写成如下形式：

(11)

对式(6)求导可得：

应变位移函数矩阵的形式如下：

(12)

其中

(13)

(14)

其中：

采用第二种位移插补方法时，应变位移函数矩阵定义与式(10)～式(12)一样。对式(9)求导可得：

应变位移函数矩阵中，对式(13)的求解与前面所述相同，式(14)的具体计算过程如下：

4 积分简化

由于板壳结构在厚度方向上的尺寸与其他方向相差较大，在分析过程中可能会出现剪切锁死或者零势能的情况。为避免这些情况的发生，通常采用简化积分的方法，即在有限元刚度积分的过程中对部分自由度采用高阶高斯积分，其余部分采用低阶高斯积分[1]。

对于壳单元而言，在对应变函数矩阵进行积分时，通常对横向剪切应力和膜应力积分部分进行简化[1]，简化方法如下：

式中：×表示不需要精确的积分;•表示需要精确积分。

5 实例分析

基于上述两种插补方法，首先在MATLAB中分别编程对壳体结构进行静力学分析，并将分析结果与ANSYS分析结果进行比较以证明分析结果的正确性。然后对两种位移插补方法的分析结果进行比较，得出两种插补方法各自的优缺点。

5.1 悬臂薄板

如图4所示，悬臂薄板尺寸为10m×10m×0.1m，其中一边固定，在相对的另一边中点处受F=1kN的集中载荷，材料弹性模量E=1×107MPa，泊松比υ=0.3。

图4 悬臂薄板模型

图5为悬臂板模型采用ANSYS分析得到的结果。图6是采用第一种插补方法通过MATLAB编程分析得到的悬臂薄板位移云图。图7为第二种插补方法通过MATLAB编程分析得到的悬臂薄板位移云图。图8所示为载荷加载处位移随着单元数变化的曲线。图9所示为分析时间与单元数目的关系曲线。由图5～9可知，对于平面薄板而言，两种插补方法得到的结果基本与ANSYS分析结果一致，第二种插补方法分析效率更高。

图5 悬臂薄板ANSYS位移云图

图6 第一种插补方法位移云图

图7 第二种插补方法位移云图

图8 两种插补方法位移收敛曲线图

图9 两种插补方法分析时间图

5.2 半圆环

如图10所示，半圆环半径为10m，两底端固定，圆环中心处受集中载荷F=1kN的作用，材料弹性模量E=1×107MPa，泊松比υ=0.3。

图10 半圆环模型

图11为半圆环模型采用ANSYS分析得到的结果。图12为第一种插补方法通过MATLAB编程分析得到的半圆环位移云图。图13为第二种插补方法通过MATLAB编程分析得到的半圆环位移云图。图14为载荷施加处位移随单元数目的变化曲线。图15为分析时间与单元数目的关系曲线。由图11～15可知，两种插补方法在最终均收敛，且收敛结果与ANSYS接近，同样，也是第二种插补方法求解效率更高。

图11 半圆环ANSYS位移云图

图12 第一种插补方法位移云图

图13 第二种插补方法位移云图

图14 两种插补方法位移收敛曲线图

图15 两种插补方法分析时间图

6 结束语

本文讨论了RM壳单元的等几何分析计算建模方法，分析了基于等几何分析的两种插补方法在收敛性、效率及在分析精度方面的差异。第一种插补方法具有较高的分析精度，但在分析过程中随着单元数目的增多，对控制节点坐标系统求解所消耗的时间也随之增多；第二种插补方法在分析过程中省去了求解控制点坐标系统的过程，因而插补模型二求解效率更高。

两种插补方法在采用单元数目较少的情况下，结果有一定差异但差异不大，而在单元划分足够细时，两种方法得到的结果基本一样，但第二种插补方法的计算效率更高。因此在分析板壳结构时，为了提升分析效率可直接采用第二种位移插补方法。

参考文献：

[1] HUGHES T J R . Isogeometric Analysis:Toward Integration of CAD and FEA[M].New York：Wiley,2009.

[2] PIEGL Les,TILLER Wayne.非均匀有理B样条[M].赵罡，穆国旺，王拉柱,译. 2版.北京：清华大学出版社,2010,1-314.

[3] SEVILLA R,MÉNDEZ S F,HUERTA A.3D NURBS- enhanced finite element method (NEFEM)[J].Int. J. Numer. Methods Engrg,2008,88:103-125.

[4] HUGHES T J R . The Finite Element Method:Linear Static and Dynamic Finite Element Analysis[M]. Englewood Cliffs,New Jersey:Prentice-Hall Inc, 1987:383-405.

[5] 吴紫俊,黄正东,左兵权,等. 等几何分析研究概述[J]. 机械工程学报,2015(5):114-129.

[6] 陈涛，莫蓉，万能. 等几何分析中Dirichlet边界条件的配点施加方法[J]. 机械工程学报. 2012(5): 157-164.

[7] DORNISCH W.An efficient and robust rotational formulation for isogeometric Reissner-Mindlin shell elements[J].Computer Methods in Applied Mechanics and Engineering, 2016,303:1-34.

[8] KANG Pilseong.Isogeometric shape optimization of trimmed shell structures[J].Struct Multidisc Optim ,2016,53:825-845.

[9] KIENDL J.Isogeometric shell analysis with Kirchhoff-love elements[J].Computer Method in Applied Mechanics and Engineering,2009,198:3902-3914.

[10] 詹双喜. 薄壳结构自由振动分析的等几何方法[C]//第16届北方七省市区力学学会学术会议.河南：郑州大学出版社，2016:237-244.

[11] BENSOND J.Blended isogeometric shells[J].Computer Method in Applied Mechanics and Engineering,2013,255:133-146.

[12] 申文斌.张量分析与弹性力学[M].北京：科学出版社， 2016.

[13] GREEN A E, ZERNA W.Theoretical Elasticity [M]. 2nd ed. New York:Dover Publications,INC.

[14] 王勖成.有限单元法[M].北京：清华大学出版社,2008,334-417.

[15] 曾攀.有限元分析及应用[M].北京：清华大学出版社,2004,144-269.

[16] MOAVENI Saeed.有限元分析——ANSYS理论与应用 [M].王崧,译.3版.北京：电子工业出版社,2013.