用曲面有限单元建立的膜结构分析理论

用曲面有限单元建立的膜结构分析理论

     叶小兵         吴健生  

  [提  要]：  本文首次采用曲面膜单元分析膜结构，在非线性有限元理论

    基础上详细推导了全部公式。对新出现的问题和已往平面膜单元分析理论

    中各主要技术环节给出了新的解决办法和改进措施。 

    关键词：  膜结构，  曲面有限元，  几何非线性有限元理论

一  引  言

    膜结构是近一段时间内获得迅猛发展的一种新型空间结构，它由于只承受张力，使得材料受力性能得到充分发挥，而且还有以下几点有利因素：(1)重量轻，允许跨越大空间和设计成各种外形；(2)施工周期短，维护方便，造价低廉；(3)可满足多方面的使用要求，易与自然环境融为一体，被国际建筑界誉为二十一世纪的建筑。

    膜结构的设计可分为三个步骤：(1)找出一个初始平衡形状；(2)各种荷载组合下的力学分析以保证安全；(3)裁剪制作。发达国家从六十年代起开始提出多种计算方法，到目前为止以有限元法为最先进，最普遍被采用的方法。而单元类型皆为三角形平面常应变单元，该方法是从刚性板壳大变形理论移植过来的。

    从本文可以看到膜结构作为只能抗拉的软壳体是不适宜采用这种平面单元的，因为对于刚性壳体来说，这种平板单元可以看成平面应力单元和平板弯曲单元的组合，其单元刚阵可以由这两种单元刚阵合并而成。而膜结构作为软壳体是不能抗弯的，只能靠薄膜曲面的曲率变化，从而引起膜表面中内力重分布来抵抗垂直于曲面的外荷载。如果还是采用这种只有平面内应力的板单元，则应变的线性部分将不反映平面外z方向位移的影响，这导致单元不包含z方向节点反力，就每个单元来说静力是不平衡的。所幸的是应变的非线性部分考虑了z向位移的影响，使得各单元合并起来的总的平衡方程通过不

断迭代能近似达到平衡，缺点是需要过多的平面内位移来满足平衡的要求，而实际情况是只需要一定的平面外和平面内的位移及曲率变化就可以了。

    考虑到这些，本文首次采用曲面膜单元，应变的线性部分引入了z向位移及单元的曲率和扭率，非线性部分仍然保留z向位移的影响项。这样无论是每个单元还是各单元合并后的平衡方程都能很容易满足，迭代次数大为减少，而变形结果也更符合真实情况。

而且由于单元内各点应力都不相同，据此判断皱折是否出现会更为精确。最后求出的每个单元的曲率和扭率对于判断初始找形的正误和优劣以及裁剪下料都能提供很多非常有用的信息。

二  三维连续介质大变形下非线性几何方程的推导

及针对膜结构的简化

 1. 为了描述膜结构大的变形，建立如图2-1所示坐标系

在三维连续介质中，取整体笛卡尔直角坐标系，任意点的位置可以用位置矢量R来表示。坐标轴方向单位矢量

                  R

图2-1

    在膜曲面上取一组正交曲线坐标系,于是R又可表示为,

沿坐标曲线切线方向的基本矢量

  2.  在正交曲线坐标系中大变形下的Green应变张量表达式

    当膜变形后，膜中任意点的位移矢量，可以用单位矢量活动标架

    上的分量来表示。其中 为三维正交曲线坐标系的拉梅参数。

                                             (1)

变形后的位置矢量,变形后的基本矢量其中    线元变形后成为新线元,在变形前的坐标系中比较二者，得到Green应变张量表达式：

                                 (3) 在给出的具体表达式之前，针对后面将采用的有限元数值方法的特点，对(3)式作一些简化处理：

  1)由于有限单元相对来说较小较扁平，所以可以假设膜曲面上的正交曲线坐标

与其投影平面上相对应的正交坐标 并无甚区别。故单元内每一点有

 不妨再设。

  2)两个位置矢量的二次导数的乘积如忽略不计，近似取为零。

  3)只在线性项中考虑曲率和扭率对的影响，而在非线性项中忽略不计曲率和扭率对的影响。

    由于膜单元只能受拉不能抗压抗弯，所以只需考虑,将(2)式和微分几何中的Weingarteu公式及Gauss求导公式代入(3)式，并考虑以上的简化处理后得到：

                                                   (4) 

其中， ，为曲面在该点处的曲率和扭率，为垂直于膜曲面方向的位移。

       三  每个膜单元的曲率和扭率的求法

    每个膜单元的表面作为空间曲面，它的方程可以表示成径向形式：

                                    (5)

由微分几何学知识可知，曲面的曲率和扭率如(2)所示，其中H1,H2为拉梅常数。

          (6)

        图3-1

   考虑到第2节中的简化假设(1),取,则曲面方程可表示为 z=z(x,y)

                      (7)

                            (8)

四  有限元位移模式和坐标变换式

    采用6节点三角形曲面等参元对解决这个问题是非常必要的也是最简便的，因为3个角节点和3个边中节点的坐标可以唯一确定这个扁壳单元的曲率和扭率，而等参元可以方便的将曲边三角形单元转化为直边三角形单元来处理。

    以3个角节点形成的平面为局部坐标系的xy平面，垂直于这个平面的方向为局部坐标系的z方向。投影于xy平面上的曲边6节点三角形单元的位移插值函数取面积坐标

为参数，6个节点的位移插值函数分别是：

                                            (9)

    在有限元单元刚度矩阵形成过程中需考虑3次坐标系变换：

1)  第2节(4)式中的是沿着曲线坐标系方向的，需转换成单元局部坐系oxyz中的。转换关系式为：

                                      (10)

2)  局部坐标系坐标由等参元的参数坐标即面积坐标来代替。

3)  局部坐标系转换到整体坐标系中

五  薄膜结构有限元分析公式

    修正的拉格郎日描述(U,L法)是以t时刻为参考构形，来求解时刻的各变量，根据能量原理，可得线性化后的虚功方程：

     (11)

式(1)中，分别为定义在t时刻构形上的格林应变的线性部分和非线性部分。

    时刻外力虚功；       时刻构形所占据的体积

 代入第4节中的位移插值函数后，可将(11)式化为：                                 (2)

    其中        

其中分别为线性应变矩阵和非线性应变矩阵，分别为柯西应力矩阵和柯西应力矢量。

    根据(2)式可求得节点位移增量U,由于式(2)是线性化处理后所得，故将时刻的近似位移代入平衡方程后必存在一个不平衡力：

                             (3)

消去  的过程即为平衡代的过程，为加速收敛过程，还使用了松驰因子:0<1,迭代公式可表示为:

          (4)

其中 为切线刚度矩阵，本文采用完全的Newton-Raphson法迭代。六  关于皱折的处理

    应用非线性有限元理论，按照第5节中推导的公式(2)-(4)可以求得每一增量步中的变形状态和应力状态。考虑到薄膜材料只有抗拉刚度这一特点，在每一增量步后必须要校核应力状态，以防出现压应力即wrinkling现象。本文用单元中精度最高的3个高斯点的主应力来判断该单元是否出现皱折。

(1) 若3个点的，则无皱折。

(2) 若有一个点的,则该单元失效。

(3) 若只有一个点的,则该单元产生皱折并以该点的方向为单向        

    受拉方向，改变应力应变关系矩阵 C

(4) 若有2个或3个点出现单向受拉状态，如果受拉方向都相同，则以该方向为整个单   

    元的单向受拉方向。如果受拉方向不同但相差不大，则取平均值为整个单元的单向

    受拉方向。如果受拉方向相差过大，则按单元失效处理。

按照处理后的应力状态和应力应变关系矩阵C来计算第5节中的切线刚度矩阵和不平衡力矢量。

七  结束语

    本文采用从力学上讲更合理也更精确的曲面膜单元取代以往广泛应用的平面膜单元理论，详细推导了理论公式，提出了数值计算的具体方案。从实际算例的结果看，两种理论确有差别，因篇幅所限，在以后的系列文章中详细论述算例分析和模型实验。

参考文献