基于独立覆盖流形法的板壳与实体单元刚性连接研究

苏海东,韦玉霞,韩陆超,颉志强

(长江科学院 材料与结构研究所，武汉 430010)

1 研究背景

梁板壳是极为常见的工程结构，其特点是厚度远小于其他方向的尺度。作为梁板壳数值分析的主要手段，有限元法应用非常广泛，但一直存在以下问题[1-3]：采用薄梁板壳理论，满足C1连续性(即导数连续)是构造单元插值函数的难点；厚梁板壳理论虽不要求C1连续，但退化到求解薄梁板壳时会出现自锁现象，采用减缩积分或多场方法等处理手段又会带来稳定性或复杂性问题；曲梁和曲壳由于控制方程较复杂(特别是曲壳)，一般难以通过控制方程来推导数值计算公式，通常用直梁或平板拼接成曲梁或曲壳，易产生几何误差进而带来力学分析误差；梁板壳单元与实体单元的连接也需特殊处理。

笔者应用独立覆盖流形法[4]，提出了梁板壳数值分析的分区级数解[5-9](以下简称“新方法”)：采用实体计算模式，只需设定各分区的多项式级数中含厚度方向坐标的级数项(或高阶项)不参与计算，就能真实模拟薄或厚梁板壳的基本变形理论，无需推导梁板壳控制方程及其数值计算公式，避免了自锁现象，近似函数在收敛性上满足C1连续。对于任意给定中面参数方程和厚度分布的曲梁和曲壳，首次将中面几何变化反映为随动的曲线坐标及其方向余弦关于整体坐标的导数运算，实现了精确几何描述下的力学分析，避免了由直梁或平板模拟曲梁或曲壳造成的几何误差。因此，新方法不仅避开了梁板壳有限元法几乎所有的基础性问题(除了与实体单元的连接尚未研究以外)，而且展现了精确几何描述的曲梁和曲壳分析的独特优势。

实际工程常常遇到三维实体和板壳组成的结构。在有限元法中，由于薄板壳与实体不同的控制方程以及板壳单元特有的转角自由度等原因，板壳的处理方式与实体有很大差别，因而两类单元的连接需要特殊处理，常用的方法有[1,3]：约束方程法(可采用罚函数法或直接引入法)和过渡单元法。在网格划分上还要求两者在连接处的单元结点必须匹配。以上问题都会带来一定的不便。

新方法在实体计算模式的统一框架下，板壳与实体之间的连接是自然的。本文在现有的板壳分析和实体分析基础上，研究板壳单元与实体单元的刚性连接，展现两者网格可以不匹配所带来的方便性，并初步展示曲壳与实体相交曲线的精确几何。另外，还要解决新方法目前在泊松比不为0时的误差问题。

2 梁板壳数值分析的分区级数解简介

2.1 独立覆盖流形法与分区级数解

在石根华先生发明的数值流形方法[10]基础上，笔者提出了构造高阶近似函数的新途径——独立覆盖流形法[4,11]：借助现代数学的流形覆盖将求解域划分为一块块分区，在每个分区内用多项式等完备级数直接逼近局部真实场函数，相邻覆盖通过窄条形(三维为窄带状)的重叠区域保证整体近似函数的连续性。

图1 3个独立覆盖及条形重叠区域Fig.1 Three independentcovers and theiroverlapping strips

如图1所示以3个覆盖情况为例，每个带颜色的大块分区称为一个独立覆盖(即覆盖的独立区域)，独立覆盖是计算分析的主要对象，重点研究独立覆盖内的级数。而覆盖间的窄条形重叠区域(可只占覆盖面积的1%～5%)仅起到连接各覆盖级数的作用。通常情况下为2个相邻覆盖之间的连接(2个以上的多个覆盖连接处由于面积很小，其近似函数的构造不必严格要求)，设2个覆盖级数分别为V1和V2，在条形重叠区域内采用单位分解法[12]形成近似函数V，表示为

(1)

式中2个覆盖之间的单位分解函数φ1和φ2为一维有限元形函数，分别表示为

(2)

式中ξ为条形厚度方向的局部坐标(以条形中点为原点)。这样，由各覆盖分区的级数联合形成的整体近似函数，通过伽辽金法逼近真实解，称为基于流形思想的“分区级数解”(这是独立覆盖流形法的内核)，可作为微分方程的数值解答。

笔者在文献[4]中指出：独立覆盖流形法的收敛性是在各覆盖采用的完备级数逼近真实解(及其导数)的基础上建立的；此收敛性不受分区形状的影响，也不受各覆盖在条形重叠区域是否错位连接的影响。因此发展了基于任意网格(包括“任意形状”、“任意连接”和“任意加密”的3个“任意”的特性)的数值计算[13]，以及误差控制下的自适应分析和自动计算[14]。

分区级数的多样性是独立覆盖流形法的一个显著优势，在实体计算模式下，选择多项式级数的合适形式，以反映梁板壳的整体变形特征[5-9]，从而避开梁板壳有限元法的诸多问题。

2.2 梁板壳分析的分区级数解

图2 整体直角坐标系下的曲梁局部坐标Fig.2 A curved beam and its local coordinates underthe global rectangular coordinate system

(3)

(4)

如图2所示，引入转换矩阵，将曲梁各点的局部坐标位移转换到整体坐标位移u1和v1，而

(5)

新方法分析曲梁的关键点在于：应变计算在对整体坐标求导时，不仅考虑局部坐标位移为变量，而且将转换矩阵中的方向余弦也作为变量，如式(6)所示。

(6)

(7)

式中：D为弹性矩阵；Ωe为单元积分区域。

式(6)涉及一些几何运算：在给出中面参数方程的条件下，只需推导中面方程关于参数坐标的导数公式(Timoshenko梁理论需2阶导数，Euler-Bernoulli梁理论涉及R求导因而需3阶导数)，代入一系列几何公式后就可得到式(6)中的各项(包括由曲线坐标的多项式定义的局部坐标位移以及方向余弦)关于整体坐标的导数，从而实现曲梁的精确几何。

以上曲梁分析方法在考虑方向余弦为常量且曲率半径无穷大时，退化为直梁。三维曲壳(含平板)分析类似[7-8]，可模拟考虑横向剪切的Reissner-Mindlin板壳理论和不考虑横向剪切的Kirchhoff-Love薄板壳理论。

3 板壳与实体的刚性连接

新方法在实体计算模式的统一框架下，板壳与实体之间的连接是自然的，只不过板壳采用局部坐标系，实体采用整体坐标系。

仍以平面问题来说明，如图3(a)所示的平面曲梁和实体单元的连接。设连接处的曲梁单元为独立覆盖1，实体单元为独立覆盖2。曲梁单元采用局部坐标系，其局部坐标位移如式(3)所示。实体单元采用整体坐标系，其位移为

(8)

式中：ank、bnk分别为待求的多项式系数；p为设定的多项式最高阶次。实际计算中对每个单项式也要进行归一化处理。与之相对应的应变为

(9)

式中：B2为应变矩阵；a2为式(8)中待求的多项式系数组成的向量；φ2对于独立覆盖而言为1。

图3 平面曲梁与实体单元的连接Fig.3 Connection between the curved beamelement and the solid element in the plane

本文提出的梁板壳与实体的刚性连接方式，是以梁板壳为主体形成覆盖重叠区域。将图3(a)中的连接部位局部放大如图3(b)所示，从梁与实体的分界线向上及向下形成覆盖重叠区域，上半部分位于梁上，下半部分将梁延伸至实体内，从而得到图3(b)中斜划线所表示的覆盖重叠区域。其中，重叠区域的大小(即梁延伸至实体内的长度)还需满足只占覆盖面积1%～5%的要求。

根据式(1)，覆盖重叠区域的位移函数为

(10)

覆盖重叠区域内的刚度子矩阵为

(11)

式中：Ω为覆盖重叠的积分区域；i=1，2；j=1，2。弹性矩阵D在上半部分采用梁的材料参数，下半部分采用实体的材料参数。

式(11)采用曲梁的局部坐标，以曲梁方式进行积分[6]：分解为沿中轴线方向的一维Gauss积分(积分点数与级数阶次有关)，以及在上述的积分点处沿梁的厚度方向的一维Gauss积分(常用2个积分点)。在每个积分点处计算被积函数，再乘以相应的积分权系数[3]。其中，由曲梁的局部坐标通过中面方程及方向余弦转换成整体坐标,即

(12)

代入式(8)中的实体多项式级数。式(2)中的ξ也可用曲梁局部坐标表示，但要注意用链式法则求其关于整体坐标的导数[7]。

为便于操作，与文献[13]采用同样的方式，不必事先形成覆盖重叠区域，而是如图3(a)所示，首先不考虑覆盖重叠，而以梁和实体的分界线向上和向下分别计算梁和实体独立覆盖的刚度矩阵。然后由程序自动插入图3(b)的覆盖重叠区域，在每个积分点处计算式(11)的刚度子矩阵，同时扣除按各自独立覆盖求出的刚度系数，以防止积分区域的重复。

以上表述的梁与实体连接方式的一大特点是梁可以任意插入到实体中，梁与实体的网格没有必要匹配。这是基于独立覆盖流形法的覆盖任意连接的特性，可以证明，当相邻覆盖的级数V1和V2均收敛于真实解时，覆盖之间的重叠区域也收敛于真实解。

空间板壳与三维实体的连接情况类似，以板壳为主体形成覆盖重叠区域，将板壳深入到实体内部，只不过式(11)的积分还要多一个维度，即首先沿着板壳中面与实体的交界线采用一维Gauss积分(参考文献[7])。整个积分在板壳的投影面上进行，映射到准确的空间位置，可以反映交界线为曲线情况下的精确几何(见下面的算例2)。

4 板壳三维弹性矩阵的修正

在前期研究中，鉴于板壳是三维实体的特殊情况，因此采用了三维实体的弹性矩阵。在泊松比为0的板壳算例中，与采用细密网格的有限元参考解非常吻合[7-8]。但当泊松比不为0时，比如文献[7]中的球面壳算例，泊松比为0.2时最大位移约有2%的误差。本文所涉及的板壳与实体连接，泊松比效应也会有一定的影响，必须解决此问题。

研究发现，这是由三维实体弹性矩阵引起的。在梁板壳基本变形理论中，厚度方向的应力可以忽略，类似于中面上的平面应力问题。因此，有限元板壳单元的弹性矩阵与实体单元不同，无论是薄膜变形还是弯曲变形，均采用了平面应力的弹性矩阵。

参考平面应力问题，在中面上定义的局部坐标系下，本文将新方法的弹性矩阵修正为式(13)形式。

(13)

式中：E为弹性模量；μ为泊松比。

定义转换矩阵[16]

T=

(14)

式中：l1、m1、n1分别为壳体中面的第一个正交曲线坐标矢量e1的方向余弦；l2、m2、n2为第二个正交曲线坐标矢量e2的方向余弦；l3、m3、n3为壳体中面法向e3的方向余弦。将式(13)转换到整体坐标系下[16]，则有

D=TTD′T。

(15)

然后再代入刚度矩阵式(7)和式(11)。

对于曲壳中面参数方程的一般性定义，比如中面方程为z=f(x,y)的一般函数形式，可能形成非正交的曲线坐标e1和e2，需要按e′2=e3×e1转换为正交坐标后再计算式(14)。

经过以上修正后，泊松比不为0的板壳算例也与细密网格的有限元参考解符合得很好(见下面的算例2)。

5 算例验证

5.1 算例1——悬臂梁

文献[3]中的典型算例，如图4所示的变截面悬臂梁，截面高度为0.7 cm的一段采用实体单元，截面高度为0.1 cm的一段采用板单元。弹性模量为3×105MPa，泊松比为0。文献[3]采用有限元法，实体划分为10个16或20结点的六面体单元，板划分为7个8结点的矩形板单元，其中在板和实体连接处采用了1个过渡板单元，与同为0.1 cm高度的实体单元相连接(网格结点必须匹配)。与全部采用板单元的模型进行了对比，结果很吻合，说明了有限元法采用过渡单元的有效性。

本文仅划分1个独立覆盖实体单元和1个独立覆盖板单元，均采用3阶多项式。按图4的布置，板单元位于实体侧面的中间位置，可见两者在连接处的网格不匹配。由于文献[3]未给出具体数值，本文采用商用有限元软件MARC划分96个板单元的结果进行对比，见图5，两者符合得很好。

图4 变截面悬臂梁[3]Fig.4 Cantilever beam with variable sections[3]

图5 悬臂梁挠度结果对比Fig.5 Comparison of thedeflections of cantilever beam

初步分析，本算例中的弯曲变形，不管对于实体还是板而言，采用3阶多项式都足以表达真实解，因此即使两者网格不匹配(板在实体的侧面连接)，在覆盖重叠区域也会逼近真实解，从而得到高精度的计算结果。

5.2 算例2——球面壳与实体的连接

文献[7]中的球面壳如图6所示，参数方程为x=Rsinθcosφ,y=Rsinθsinφ,z=Rcosθ。θ=π/4～π/2,φ=0～π/2，半径R为1 m，厚度为0.01 m。φ=0和φ=π/2两端固支，壳体上表面承受均布压力10 kPa。弹性模量为106kPa，泊松比为 0.2。

壳体底部(z=0 m)连接一块长、宽各为1.5 m、高度为0.5 m的实体基座，基座的底部全约束，如图7所示，采用MARC软件建立细密网格的三维有限元模型用于计算参考解，其中壳体单元2 500个，实体单元5 556个。

图6 球面壳Fig.6 A spherical shell图7 球面壳与实体基座的有限元模型Fig.7 The FEM model of thespherical shell with a solid base

如图8所示，在参数坐标φ-θ的投影平面上，划分5×3个独立覆盖。首先不考虑实体基座，在均布压力作用下，最大位移发生在壳体底面的中点处(φ=π/4)，整体坐标x、y、z这3个方向的位移分别为u=v=-6.458 mm(-6.459 mm)，w=-3.066 mm(-3.069 mm)(括号内为MARC参考解,采用了2 500个壳体有限单元,下同)。可见在5×3独立覆盖的网格密度下，计算精度很高，相对误差<0.1%。本算例的泊松比不为0，如果采用以往的三维弹性矩阵，则最大位移为u=v=-6.338 mm，w=-3.015 mm，约有2%的相对误差。可见本文对弹性矩阵修正后的效果。

图8 φ-θ投影平面上的5×3个独立覆盖Fig.8 5×3 independent covers on φ-θ projection plane

再考虑实体基座。首先，实体仅为1个独立覆盖，如图9(a)所示的x-y平面上的网格1，交线为精确几何的圆弧段(对应φ方向的5个独立覆盖)，可见壳体网格与实体网格的不匹配。壳体和实体均采用6阶多项式。在受到基座的约束后，最大位移发生在壳体顶面中点处。当实体的材料参数与壳体相同时，计算得到的最大位移为u=v=-1.720 mm(-1.732 mm)，w=-1.241 mm(-1.247 mm)，可见与有限元参考解吻合良好，相对误差<0.7%。

将实体基座的弹性模量降低10倍，为105kPa时，最大位移为u=v=-1.689 mm(-1.637 mm)，w=-1.222 mm(-1.192 mm)，与有限元参考解约有3%的相对误差。初步分析，基座刚度降低10倍，造成连接部位附近的变形较为复杂，1个实体独立覆盖难以很好地反映。因此需要对实体进行网格加密。计算结果如下：

如图9(b)所示的网格2，在x-y平面内将其均匀划分为4块，计算得到的最大位移为u=v=-1.666 mm，w=-1.209 mm，仍约有1.7%的相对误差。

如图9(c)所示的网格3，均匀划分为16块，u=v=-1.653 mm，w=-1.200 mm，相对误差降至1%以下。

进一步探究误差是否还会下降，在连接处局部加密网格，如图9(d)所示的网格4，则u=v=-1.636 mm，w=-1.190 mm，相对误差降到0.1%左右。

图9 实体基座在x-y平面上的网格划分及圆弧交线Fig.9 Mesh division of the solid base and the arcintersection on x-y plane

图10(a)给出球面壳在φ=π/4的中间断面处的位移对比，图10(b)为顶面x=0～0.5 m线上的位移对比，可见符合得很好。

图10 球面壳位移对比Fig.10 Comparison of the displacements ofspherical shell

综上所述，本文提出的板壳与实体的刚性连接方法是可行的，且板壳单元和实体单元的网格可以不匹配，这样非常有利于前处理工作，在网格划分达到一定密度的情况下能得到高精度的计算结果。

6 结 语

在笔者前期提出的梁板壳数值分析的分区级数解方法基础上，本文研究了板壳与实体单元的刚性连接方法：由于板壳也采用了实体计算模式，因此与实体之间通过覆盖重叠区域自然连接，只不过板壳采用局部坐标系，实体采用整体坐标系，最终都要转换到整体坐标系下进行统一的连接及分析；基于独立覆盖流形法的覆盖任意连接的特性，将板壳插入到实体中形成覆盖重叠区域，同时采用板壳积分方式进行积分；实体单元与板壳单元可以各自划分网格，在连接处不要求匹配，这样非常有利于前处理工作，在网格划分达到一定密度的情况下能得到高精度的计算结果，即在连接处能准确反映变形，将来可以采用误差分析方式进行自适应分析。本文顺带修正了新方法的三维弹性矩阵。

本文还通过球面壳与实体的水平表面相交的算例，初步展示了板壳与实体相交曲线的精确几何。下一步将研究空间曲壳-曲壳、曲壳-实体连接的一般情况，即2个空间曲面(实体表面也为曲面)之间交线的精确几何，进一步体现新方法在这一点上的优势。

致谢：感谢石根华先生的指导！