苏海东，韩陆超，颉志强

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

1 研究背景

梁板壳结构分析通常基于一些被实验验证的基本变形假设(已归纳为具有很高模拟精度的近似理论)[1]：对于经典的薄梁板壳，有Euler-Bernoulli梁理论和Kirchhoff-Love板壳理论，要求原垂直于中面的横截面在变形过程中始终保持为平面且仍垂直于中面，忽略横向剪应变；对于需要考虑横向剪应变的厚梁板壳，采用Timoshenko梁理论和Reissner-Mindlin板壳理论，只要求横截面保持为平面。当然，前者也可当作后者在厚度变薄时的特殊情况。此外，沿厚度方向产生的直接应变可忽略。

从薄梁板壳假设出发，前人推导了不同于实体分析的控制方程[1]，最终都归结为关于挠度的4阶微分方程，由此带来在有限元数值分析中要求插值函数满足C1连续性(即导数连续)的难题，至今都未能得到很好的解决[2-3]。

进一步对于曲梁和曲壳而言，由于存在初始曲率，其控制方程的推导本身就比较复杂(特别是曲壳[4])，因而一般不采用由控制方程建立积分弱形式的标准方式，而是用直梁或平板单元近似地模拟曲梁或曲壳，容易产生几何误差，进而带来力学分析上的误差。

笔者基于独立覆盖流形法和分区级数解思想，提出梁板壳结构数值计算的新方法[5-7]，针对一般的中面参数方程的几何描述，将中面的精确几何变化在应变计算中反映为中面局部坐标系的方向余弦关于整体坐标的导数运算，首次实现了精确几何的曲梁和曲壳分析，但目前仅限于厚梁板壳假设[6-7]，还未在经典的薄梁板壳假设中实现，这是该方法的基础理论部分尚未完成的工作。另外，对于各种中面的曲线或曲面描述，其相关的几何公式都需要进行一遍复杂繁琐的人工推导。

本文针对上述问题开展研究。首先简要介绍独立覆盖流形法及分区级数解的思想，并讨论其近似函数的C1连续性，然后以平面薄曲梁为例详细给出计算公式，并简述薄曲壳的计算过程，最后通过典型算例验证方法的有效性。

2 独立覆盖流形法及分区级数解简介

1991年，数值流形方法[8-9](以下简称流形法)的发明者石根华将现代数学的流形思想引入数值计算：如图1(a)所示，一系列相互重叠的覆盖用于分割求解域(要求覆盖整体上包含求解域)，使整体分布复杂的物理场在每个覆盖区域内简单化，以利于局部近似函数的逼近。

从2011年开始，笔者在流形法中首次引入“独立覆盖”[10]，即覆盖的独立区域(非重叠部分)，如图1(b)所示，每个覆盖中的“蓝色”区域就是独立覆盖，其近似函数通常采用多项式等完备级数。每个覆盖中的“黄色”区域是该覆盖与周边覆盖相重叠的部分，在此重叠区域内，采用单位分解函数φi(在流形法中又称为权函数)将各覆盖中的级数Vi联系起来[11]，以保证整体近似函数V的连续性,即

(1)

式中n为重叠的覆盖数。

笔者进而提出独立覆盖流形法[12]，强调独立覆盖是计算分析的主要对象：如图1(c)所示，独立覆盖占据一个覆盖的主要面积，重点研究独立覆盖内的级数，覆盖之间的重叠区域形成面积尽可能小(可只占覆盖面积的1%)的窄“条形”，仅起到连接各覆盖级数的作用，并以有限单元的线性形函数作为权函数，比如两个覆盖之间采用一维有限元的形函数[2]，即:

(2)

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

图1 流形覆盖示意图Fig.1 Manifold covers

笔者在文献[12]中讨论了该方法的收敛性，指出：随着各分区级数的阶次升高，整体近似函数逐步逼近真实解；不仅场函数本身是收敛的，其导数也是收敛的；收敛性不受分区形状的影响，也不受各覆盖在条形重叠区域是否错位连接的影响。因此发展出了基于任意网格(任意形状、任意连接和任意加密)的数值计算[13]，以及无需人工参与的自动计算[14]。

采用分区级数解的一大优势是，可事先构造出满足物理场整体或局部特性的级数形式。笔者提出的梁板壳数值计算新方法，对于多项式级数中包含厚度方向坐标的各项，仅在模拟板壳的薄膜位移(或梁的轴线向位移)时保留至厚度方向坐标的1阶项，就可将厚梁板壳的基本变形假设固化在近似函数中[5-7]。同时在文献[5]的直梁分析中，尝试将所有含厚度方向坐标的级数项去掉，实现了细长梁(以下也称为薄梁)的Euler-Bernoulli梁理论，对于其所要求的C1连续性，当时的解释是“基于导数的收敛性，近似函数的导数趋向于连续的真实解”，本文补充说明如下。

图2 2个覆盖及其条形连接Fig.2 Two covers and a stripbetween them

在条形重叠区域内有

V=φ1V1+φ2V2。

(3)

不难验证C0连续性。条形区域关于x方向的一阶导数为

在完全收敛的情况下，在条形区域内由于真实解V*的连续性，有V1=V2→V*。即使在未完全收敛的情况下，考虑到条形区域很窄，则在达到一定阶次时，条形处的V1和V2非常接近，即此处V1≈V2，因此以上情况都有

(5)

3 薄曲梁的分区级数表达及应变矩阵计算

如图3所示，在整体直角坐标系x-y中，由曲线参数方程描述中面坐标(x0,y0)，建立沿中面曲线变化的局部正交坐标系xi-yi[6]，其中xi沿梁的轴线方向(中面曲线的切向)，yi沿厚度方向(中面曲线的法向)，xi轴在整体坐标系下的方向余弦为cosαxi、cosβxi，yi表示梁上的点到中面的距离，yi轴的方向余弦为cosαyi、cosβyi。

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

根据梁保持平截面的假设，局部坐标系下的两个方向位移为[2]

(6)

其中，对于满足Euler-Bernoulli理论的薄曲梁，有[4]

(7)

式中：一维多项式f1(xi)和f3(xi)分别表示梁中面的轴线方向和厚度方向的位移;f2(xi)=θ(xi)表示截面转角；R为中面点处的曲率半径；h为弧长微分系数。

因此只需考虑f1(xi)和f3(xi)，如图4所示，对于多项式级数，在轴线和厚度方向去掉所有关于yi的项，只让关于xi的一维级数(图4中的虚线左侧)参与计算(式(8)),就能通过式(6)和式(7)实现薄梁的变形假设，而无需推导薄梁控制方程。

(8)

式(6)—式(8)写成矩阵形式,则有

(9)

分别记cosαxi、cosβxi为cosα、cosβ，同时考虑两个局部坐标轴的垂直关系，将局部坐标下的位移转换到整体坐标下的位移，即

(10)

并记其中的函数公式(对于第q项)为

二维应变子矩阵(对于第q项)为

(12)

其中，由链式法则，各项对整体坐标的偏导数为

(13)

形成3×3数组D。其中的数组运算(式(14)中j=1～3)为

(14)

然后D数组的每一项再与C数组的每一项按式(14)的方式进行组合。最后依次放入式(12)的对应位置，形成应变子矩阵Bq，并代入刚度矩阵公式计算单元刚度矩阵[6]。

对于连接覆盖1和覆盖2的条形区域，求关于整体坐标x或y的偏导数时，Bq矩阵每一项系数c还要与权系数φi(i=1,2)进行组合，形成应变子矩阵Biq，则有

(15)

其他关于积分方法等内容参见文献[6]。另外，取R→∞和h=1可退化到计算细长直梁。

4 曲梁的几何计算公式

(1)给定中面曲线的参数方程为

(16)

式中t为参数。取局部坐标xi=t，yi=r，r为曲梁上的各点到中面的距离。

(2)关于t的1—3阶导数为：

(3)弧长的微分系数h表达式为

(18)

(4)曲率半径R为

(19)

(5)方向余弦为:

(20)

(6)局部坐标关于整体坐标的导数。曲梁上的任一点为

用隐函数求导法分别对x和y求偏导，并求解方程组，令g1=x′t-r(cosβ)′t，g2=y′t+r(cosα)′t，g=g1cosα+g2cosβ，则：

(22)

综上所述，在给定中面参数方程后，人工推导仅限于式(17)中相对简单的各阶求导，其他运算都由程序自动完成。因此，本文方法解决了以往在几何公式推导上的难题，对于任意给定的中面参数方程都具有通用性。

另外，本文参考文献[13]，采用边界条的方式严格施加边界条件。如式(23)所示的在边界xi=xc处的边界覆盖级数(一般取为2阶多项式即可，此处l为边界条的宽度)，若施加简支条件，则仅需令系数a0=b0=0，其余各项在xi=xc处自动为0，从而严格满足边界条件。若在xi处施加固支条件，考虑到式(7)中的f3(xi)的导数对ui的影响，还需令b1=0。

(23)

5 薄曲壳计算简述

如图5所示的曲壳中面，xi和yi为正交曲线坐标，zi表示壳体上的点(整体坐标为(x,y,z))到中面的距离[7]。给定中面参数方程为

(24)

图5 整体直角坐标系下的曲壳及其局部坐标Fig.5 A curved shell and its local coordinates underthe global rectangular coordinate system

按薄板壳的Kirchhoff-Love理论，定义局部坐标系下位移覆盖函数为

(26)

其中，对于正交的局部坐标系，两个方向的转角为[4]

式中：h1、h2分别为xi、yi两个局部坐标的拉梅系数；R1、R2分别为两个局部坐标方向的曲率半径(对于正交的局部坐标系，本文为两个主曲率半径)[4]，其表达式为

(28)

其中：

(29)

将局部坐标系下的位移转化到整体坐标系下的位移，并记其中的函数公式(对于第q项)为

(30)

式中：fq=fq(xi,yi)为关于(xi,yi)的第q项多项式；cosαxi、cosβxi、cosγxi分别是xi轴关于整体坐标x、y、z的方向余弦，其他以此类推。

薄曲壳计算的关键是三维应变子矩阵Bq[7]。整体的计算思路与上两节的薄曲梁相似，只是过程更复杂一些：同样要拆解成数组运算并进行组合，最后放入Bq矩阵的相应位置；在几何运算中，需要事先推导中面方程关于局部坐标的1—3阶导数公式(曲率半径的微分需要3阶导数[15])，之后的运算都由程序自动完成。除了式(28)和式(29)以外，其他几何量的表达式见文献[7]。

本文参考文献[13]，将曲梁和曲壳的计算程序改写成标准的二维独立覆盖及条形计算方式，而不采用以往基于强制约束形成独立覆盖及条形的方式[5-7]，这样可以大幅减少计算量，并能实现任意网格划分，以及通过边界条严格施加边界条件。

6 算 例

采用薄梁板壳假设，将文献[6]和文献[7]的曲梁和曲壳主要算例重新计算一遍。

6.1 常曲率的曲梁算例——圆形曲梁

如图6所示，梁的矩形截面厚度为0.01 m。左、右两端采用边界条方式施加固端约束。法向集中力F=1 kN作用于梁的中点。弹性模量E=105kN/m2，泊松比μ=0。荷载作用点处的位移及与细密网格的有限元结果对比见表1。

图6 圆形曲梁及其变形Fig.6 Circular curved beam and its deformation

表1 荷载作用点处的位移比较(圆形曲梁)Table 1 Comparison of the displacements of loadingaction point of the circular-curved beam

6.2 变曲率的曲梁算例——椭圆形曲梁

设椭圆的长轴a=2 m，短轴b=1 m，梁的矩形截面厚度为0.01 m。弹性模量E=108kN/m2，泊松比μ=0。如图7所示，只用1个独立覆盖，左端全约束，右端分别承受水平力和竖向力F=1 kN。荷载作用点处的位移及其与细密网格的有限元结果对比见表2。

图7 椭圆形曲梁及其变形Fig.7 Ellipse-curved beam and its deformation

表2 荷载作用点处的位移比较(椭圆形曲梁)Table 2 Comparison of the displacements of loadingaction point of an ellipse-curved beam

6.3 球面壳算例

计算如图8(a)所示的第1象限内的球面壳，半径为1 m，θ=π/4～π/2，φ=0～π/2。壳体厚度为0.01 m，左端(φ=0)和右端(φ=π/2)为固支约束，球面壳上表面承受法向均布压力P=10 kN/m2。弹性模量E=106kN/m2，泊松比μ=0。图8(b)显示了用于计算球面壳参考解的划分为5 000个平板单元的有限元网格。

图8 球坐标及第1象限内的球面壳有限元网格Fig.8 Spherical coordinates and finite element meshesof a spherical shell in the first quadrant

本文在φ-θ的投影平面上分别采用3种网格：网格a——1个独立覆盖网格；网格b——均匀划分为5×3(φ向×θ向)网格；网格c——均匀划分为9×5网格。计算得到的各方向最大位移见表3，出现在(θ=π/2，φ=π/4)点处，均与细密网格的有限元解非常接近。以其中的网格a为例，本文方法的自由度为417，而文献[7]采用厚壳假设时为669。

表3 球壳最大位移计算结果比较Table 3 Comparison of the maximum displacements ofthe spherical shell

将球壳两端改为简支约束，采用网格a计算，则9阶情况下的最大位移u=v=-0.654 cm，w=-0.321 cm，而细密网格的有限元解为u=v=-0.658 cm，w=-0.320 cm。

程序也可用于计算平板。采用网格c，平板范围为x=(0～π/2)m，y=(π/4～π/2)m，考虑四周简支约束，其他条件不变。按薄板理论，采用4阶多项式计算得到中点位移为0.460 2 m，与有限元参考解0.460 5 m非常接近。

与文献[6]和文献[7]相比，上述各算例在两种变形假设下的计算结果几乎一致，本文基于薄梁板壳假设可节省自由度约30%(对于平板弯曲算例，还可去掉所有关于薄膜位移的自由度)，但单元计算要复杂一些。

7 结 语

本文在精确几何曲梁和曲壳的前期研究基础上，完成了基于Euler-Bernoulli梁理论和Kirchhoff-Love板壳理论的薄曲梁和曲壳分析，不仅补上了梁板壳数值计算新方法在基础理论上的最后一块拼图，而且解决了几何公式推导复杂的问题(人工只需推导中面方程关于参数坐标的最多3阶导数)。以下再次总结新方法的特点和优势：

(1)采用实体计算模式，只需设定各分区的多项式级数中含厚度方向坐标的级数项不参与计算(或仅在薄膜位移中保留厚度方向坐标的1阶项)，就能模拟薄或厚梁板壳的2种基本假设，而无需梁板壳控制方程及推导相应的数值计算公式。

(2)在收敛性方面具备C1连续，且用厚梁板壳假设计算薄梁板壳时无闭锁现象，也没有通常实体计算中厚度远小于其他尺度导致的数值病态。

(3)对于任意给定中面参数方程和厚度分布的曲梁和曲壳(目前对于薄曲壳假设情况，中面参数方程暂定为正交曲线坐标描述)，两种基本假设情况均实现了精确几何描述下的力学分析，避免了由直梁或平板单元近似模拟曲梁或曲壳造成的几何误差。

因此，基于独立覆盖“分区级数解”的梁板壳数值计算新方法不仅避免了现有方法的诸多难题，而且展现了在精确几何曲梁和曲壳分析方面的独特优势。近期将进行梁板壳与一般实体单元连接以及梁板壳相互连接的算例研究，然后再结合独立覆盖流形法在任意网格划分上的优势开展梁板壳自适应分析以及自动计算研究，并体现其高阶级数收敛快、自由度少的优势。

本文研究的另一个重要意义在于：对于有限元法而言，构造C1连续的近似函数一直存在比较大的困难，而独立覆盖流形法具有导数自然收敛的特性，无需特殊操作就可实现C1连续，这为4阶微分方程的数值求解打下了很好的基础。

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