ANSYS-膜结构找形分析技术

一．前言

索膜结构是应用非常广泛的一种结构形式，由于索膜在无应力情况下没有刚度，不具有承载力和一定的形状，所以必须施加适当的预应力来使其产生足够的刚度并确定性状。

其设计中主要涉及到三个关键环节：找形、载荷分析、裁剪分析。

找形也叫形态分析，指的是给定预应力分布以及控制点(即约束点，通常为实际的支座点)坐标，通过适当的方法确定该预应力分布下索膜结构的平衡形态。

载荷分析是用来分析预应力索膜结构在外载荷作用下的应力、位移，确定其承载能力，以验证结构是否具有足够的刚度（此刚度为预应力刚度）以及在外载荷作用下是否会出现皱褶。

裁剪分析是将通过找形荷载荷分析确定的膜面（通常为不可展曲面）预应力释放，并根据几何拓扑理论进行适当的剖分，然后将剖开的膜片展开为平面，作为施工下料的依据，以保证施工安装后的膜面形状与分析得到的形状相吻合。

裁剪分析涉及到几何拓扑理论，单纯用有限元分析工具解决有一定困难，需要综合有限元工具以及其它一些特殊的技术来解决，需要有限元以外的专门研究与开发，所以本文不准备讨论ANSYS在这方面的应用。

载荷分析是一个预应力索膜结构的典型有限元分析，对于ANSYS没有任何困难，只要在给预应力作为初应力施加于相应单元，同时打开大变形效应，施加其它载荷，ANSYS就会考虑预应力的预应力刚度进行计算得到相应的结果。所以没有任何困难，本文也将不去讨论。

找形是载荷分析和裁剪分析的基础，是索膜设计的出发点，也是一个难点，需要找到给定预应力分布下的平衡形态，因为预先并不知道该形态，在初设形态下预应力一般不能平衡，需要通过适当的方法进行迭代计算来确定能够使预应力分布平衡的位移形态。本文将探讨这种计算方法，并给出ANSYS解决方案以及相应的验证算例。

二．ANSYS的找形方法

1．单元类型

   采用SHELL181和LINK180。其原因如下分析。

   ANSYS提供了膜单元（SHELL41）以及其它的壳单元（SHELL181、SHELL63等），膜单元考虑了膜的性质，不抗弯、不抗压。但SHELL41单元不能直接指定初应力，而在找形分析中初应力需要用来模拟预应力，是必不可少的，虽然可以通过温度应力来达到预应力的目的但在给定预应力的前提下很难确定什么样的温度分布可以产生相应的预应力，所以也不实用。

   SHELL181单元是一个壳单元，可以用于薄壳到厚壳结构，其单元公式中包含了抗弯刚度和抗压刚度。但我们可以通过采用极薄的壳来克服这些刚度。

首先实际的膜结构是没有弯曲刚度的，其受力形态为张拉力，通常膜结构的刚度都可以用面刚度给出即E×T的值，其中E为弹性模量，T为膜的厚度。所以我们可以设定SHELL181单元的厚度为很小的值，同时按比例调整材料的弹性模量，只要保证E×T的值不变就可以了。而在厚度很小的情况下壳单元的弯曲刚度会很小以致可以忽略不计，这样就达到了没有弯曲刚度的要求，根据经验，让单元厚度小于结构典型尺寸的千分之一就可以满足要求。

一般认为索膜是不抗压的，但不抗压的结果是索膜在面内压力作用下会起皱，是一种弯曲现象，所以从本质上来讲，不抗压的表面现象是由于没有弯曲刚度引起的（类似于屈曲），因此只要如上将壳的厚度变得很小，则壳实际上就没有抗压强度了，因为平衡态是一个空间曲面，一旦产生，结构就会产生弯曲（褶皱）。而且用这种单元还有一个好处，它保留的极小的弯曲刚度还会有利于计算稳定，为载荷分析中可能出现的皱褶情况下地计算收敛提供帮助，有助于模拟出皱褶形状（这正好反映了这样一种情况：实际的索膜结构弯曲刚度也不是完全为0）。

LINK180单元是一个杆索单元，其单元公式中就没有包括抗弯刚度，所以不必进行什么特别的处理就适合于模拟索，这里就不用多说了。

这两种单元都可以指定初应力，为找形分析提供了必备条件。

2．找形方法步骤

  2.1．建立初设形态的模型

  初设形态由支座点控制，以支座点为控制关键点，建立初设形态几何模型（包括面、线）并划分单元，注意具有不同预应力的单元最好指定不同的材料，以方便指定初应力。单元应该采用三角形形状，因为索膜结构空间曲面可能扭曲，用四边形时计算过程中可能由于单元扭曲出现问题。

  2.2．施加约束（通常在支座点）。基于索膜的特点，仅约束平动自由度。

  2.3．指定初应力。需注意的是，ANSYS初应力基于单元坐标系。另外，如果结构中每种材料具有同样的初应力，则可以直接用命令ISTRESS完成初应力施加，如果每种材料具有不同的初应力，则需要通过初应力文件来施加，可以通过命令流自动生成初应力文件。具体参见后面的算例命令流。

  2.4．打开大变形效应，进行计算。通常初设形态下预应力不会平衡，会产生位移，预应力会由于位移而释放。也就是说，计算结束后应力状态将和预应力不一致。

  2.5．开始迭代计算，如下：

     2.5.1．开始一个新分析，按照前一个形态的计算位移结果更新体型，采用UPGEOM命令。

     2.5.2．设定初应力（同前一个形态，由于开始新分析时，初应力会归零，所以需要重新设置）

     2.5.3. 判断是否收敛。由于平衡的形态下预应力不会引起任何位移，所以可以用结构内的最大位移来判断是否收敛，即若当前形态下产生的最大位移小于初设形态下最大位移的某一个百分比或小于某一个适当的绝对值，则认为收敛，一般来说，这个百分比可取为0.1%－1.0%，另外绝对值可取为1e-4－1e-6，这个绝对值是为了防止初始形态即为平衡形态的情况，具体取值取决于精度要求。只要两个条件满足了一个，即可认为收敛

     2.5.4. 如果收敛，则最后更新一次体型，完成找形，新的体型则为平衡态体型。如果没有收敛，回到2.5步继续迭代

三．找形分析实例

1．悬链面薄膜结构找形

 悬链面为一个上下端固定的等应力面，是由悬链线绕Z轴旋转得到的面，其方程为：

                                                        (1)

本例中取a=5m, b=30m, 可计算出h=12.3894m。以这三个数据可以确定悬链面的上下边缘，通过找形找出平衡态。

膜材为面内各向同性材料，其刚度为：E×T＝2.36×105 N/m, 泊松比为0.4

悬链面的平衡态形状与预应力大小无关（当然预应力大小会对结构的受力特性产生影响，这不是本文的讨论范围），本例取sx×T=sy×T=2×104 N/m.

计算中壳单元的厚度采用0.0001m，则材料的弹性模量应该为E=2.36×109 N/m2，预应力大小应该为sx =sy =2×108 N/m2。由于模型对称，采用了1/4对称模型进行分析。

分析采用的命令流为文件membrane.txt，其有限元模型如图1所示。这个例子在迭代计算33次后收敛，图2为收敛后的体型，即平衡态体型，图3和图4为平衡态的预应力分布，可以看到已经非常接近要求的预应力值，其误差可以忽略不计。

                                                                                图1: 1/4有限元模型

                                                                                    图2: 平衡态体型

                                                                                  图3: 平衡态X预应力分布

                                                                                     图4: 平衡态Y预应力分布

表1为边界线上节点坐标位置与理论值的比较，可以看到二者是非常接近的。

表1：ANSYS计算位置与理论位置的比较

节点	1	56	57	58	59	60	61	62	63	64
X	5.000	5.395	6.354	7.614	9.008	10.463	11.946	13.440	14.940	16.440
计算Y	12.389	10.505	8.870	7.550	6.484	5.601	4.852	4.204	3.633	3.124
理论Y	12.389	10.414	8.788	7.477	6.420	5.545	4.804	4.162	3.597	3.092
节点	65	66	67	68	69	70	71	72	38
X	17.947	19.453	20.958	22.465	23.971	25.477	26.984	28.492	30.000
计算Y	2.661	2.238	1.850	1.489	1.154	0.839	0.544	0.264	0
理论Y	2.634	2.216	1.831	1.474	1.142	0.831	0.538	0.262	0

2．马鞍形索网结构找形

图5所示为一个索网结构的初设形态，其轴线方向长86.8m，圆弧形拱高度为28.49m，拱宽为97.0m，四周的索端均固定，在轴向索和环向索中具有一定预应力时，索网会呈现出马鞍形的平衡态，本例将找出这个平衡形状。

索的截面面积为67.4×10-6 m2，弹性模量1.7×1011 N/m2。本例的命令流见文件cable.txt，需要说明的是，由于通常情况下，索膜结构的平衡态只与预应力的相对分布有关而与具体值无关，所以本例中利用了这个特点，在迭代计算中将预应力一致增大了1000倍来加快收敛速度（当然最大倍数与具体问题有关，如果过大也可能引起计算不收敛，可以根据具体问题初设体型的位移确定），在收敛结束后再用实际的预应力进行一次计算即可。另外，采用了读入自动生成的初应力文件的方式给不同的索施加不同的预应力。

首先设定轴向索和环向索中均具有5000N的预张力（即74.2×106 N/m2的预应力）。图6为找形得到的平衡态，图7为这种平衡态的预应力分布，可以看到与设定预应力非常接近。

图8为索中预张力为2500N（37.1×106 N/m2的预应力）情况下找形得到的平衡态，其马鞍形顶部轴向线上最低点的高度11.278m，而图7中该点的高度为11.260m，考虑到迭代误差，二者可以认为是相等的，这说明平衡态与预应力大小无关。

但预应力相对分布对结构成形有很大影响，图9为轴向索与环向索预张力4:1（分别为5000N和1250N）时的平衡态，图10为轴向索与环向索预张力100:1（分别为5000N和50N）时的平衡态，可以看出，索膜结构形状的调整可以通过改变预应力的相对分布类获得。  

                                                                                     图5: 索网初设形态

                                                                                   图6: 等预张力平衡形态

                                                                                     图7: 平衡态预应力分布

                                                                                      图8: 预应力减小一半的平衡形态

                                                                                       图9: 张力比4:1

                                                                                           图10: 张力比100:1

3．伞形索膜组合结构找形

一柔性边界的六边形伞形索膜组合结构，初设形状如图11所示，其六边形平面的外接圆半径为5m，中心是一个半径0.5m的钢环，升起高度为3.5m，模型的边界布置6根边索，内部布置6根脊索，中央钢环和六边形六个角点完全固定。

膜材为面内各向同性材料，其刚度为：E×T＝2.36×105 N/m, 泊松比为0.4，膜材中的预应力为：sx×T=sy×T=2×103 N/m. 

计算中壳单元的厚度采用0.0001m，则膜材的弹性模量应该为E=2.36×109 N/m2，预应力大小应该为sx =sy =2×107 N/m2。

边索和脊索的截面面积均为67.4×10-6m2，弹性模量1.7×1011 N/m2。脊索中具有5000N的预张力（即74.2×106 N/m2的预应力），边索中具有10000N的预张力（即148.4×106 N/m2的预应力）

分析所用的命令流文件为mem-cable.txt，类似于例2，在采用了初应力文件施加不同的预应力，在迭代过程中将预应力等比例扩大了10倍来加快收敛速度。

图12为上述设定预应力下的平衡态，图13、图14为平衡态下的膜面两向预应力分布，图15为索中预应力分布，可以看到平衡态下的预应力满足设定要求。

                                                                                    图12: 平衡形态

                                                                             图13: 平衡态膜面x向应力分布

                                                                               图14: 平衡态膜面Y向应力分析

                                                                                图15: 平衡态索中应力分布

结构形状的改变可以通过改变膜和索的预应力比值来实现，从而得到形态各异的形态。在保持索的预张力不变的情况下，增加膜的预张力将增大结构的曲率，例如将膜的预应力增加至原来的1.5倍，即sx =sy =3×107 N/m2，得到的成形图如图16所示。反之，在保持索的预张力不变的情况下，降低膜的预张力将减小结构的曲率，使得曲面变得较为平坦，例如将膜的预应力减小到原来的一半，即sx =sy =1×107 N/m2，得到的成形图如图17所示。

                                                                          图16: 膜面预应力增加至原来的1.5倍

                                                                            图17: 膜面预应力减小至原来的一半

四．结论

本文针对建筑结构中常见的索膜结构的关键问题－找形，提出了ANSYS进行找形分析的解决方案，并用实例验证了这种方法的有效性。以此为基础，可以进一步使用ANSYS的强大的结构分析、多物理场分析技术非常方便有效地针对找形得到的索膜结构平衡形态，考虑该结构与其它结构的相互作用以及复杂载荷物理环境，从系统级、耦合的层次进行结构分析以校核其受力性能、安全性能，为设计提供依据。

有效地实现索膜找形分析也将进一步拓展ANSYS在建筑领域的应用前景。