膜结构是一种非传统的结构形式, 其造型优美, 极富现代气息;结构轻巧, 有极强的空间跨越能力;易于建造搬迁, 有明显的经济效益。因此膜结构建筑成为近几十年蓬勃发展的一种新型的大跨度空间形式。

膜结构的设计也迥异于传统的刚性结构。由于膜材的柔性特征, 其本身没有抗压刚度和抗弯刚度, 需赋予一定的预拉力才能形成确定的空间曲面形状和抵抗外荷载的能力, 成为真正的结构。这种在一定预应力作用下, 找出的一个既符合建筑美观, 又满足边界条件和力学平衡的过程就是找形分析过程。

1970年, E·Haug和G·H·Powell提出了一种基于Newton-Raphson非线性迭代的索膜结构找形方法。它是针对索膜结构具有强烈的几何非线性的特点, 在小应变、大位移的情况下, 首先将膜结构离散成由节点和三角形单元构成的空间结构, 设定一初始应力分布, 采用拉格朗日法建立非线性方程组, 结合边界条件迭代求解。由于是非线性求解, 无法避免收敛问题, 往往可以通过分段提升的方法来保证收敛性。现有的一些通用有限元软件, 如ANSYS, SAP2000, ADIANA等都有用这种方法来进行膜结构的找形。

非线性有限元法找形的基本思路是先将索膜结构离散成若干单元, 然后针对索膜结构的小应变、大位移特点, 百富策略白菜网几何非线性理论, 建立以节点位移为基本未知量的非线性有限元方程组, 最后用迭代计算方法并结合边界条件求解。由于初始平衡状态是纯力学平衡问题, 与所采用的材料无关, 故在计算过程中采用小弹性模量法, 以便使结构自由变形, 达到平衡。

膜单元采用三节点九自由度的三角形平面单元, 索单元采用二节点六自由度的只拉不压空间直杆单元。忽略材料非线性, 考虑几何非线性, 有:

{Δσ}=D{Δε}{Δσ}=D{Δε}

{Δε}=([BL+BNL]){ΔU}{Δε}=([BL+BΝL]){ΔU}

式中{Δσ}{Δσ}为单元应力增量矢量;{Δε}{Δε}为单元应变增量;D为材料本构矩阵, BL为线性应变位移转换矩阵;BNL为非线性应变位移转换矩阵;{ΔU}{ΔU}为单元节点位移增量矢量。根据虚功原理, 推导出非线性位移法找形的UL方程:

([KL+KNL]){ΔU}={R}−{F}([ΚL+ΚΝL]){ΔU}={R}-{F}

式中:{R}{R}为t+Δt时刻的荷载等效节点力矢量, 此项在找形过程中可忽略不计;{F}{F}为t时刻的单元应力节点力矢量;{KL}{ΚL}为线性应变增量刚度矩阵;{KNL}{ΚΝL}为非线性几何刚度矩阵或初应力矩阵。

ANSYS是一种基于非线性有限元思想的通用有限元软件, 可用于索膜结构的找形分形。其基本分析原理是:先用小弹性模量技术, 将目标节点提升到指定高度, 用支座移动法进行初步找形, 目标点固定, 其它点连动, 得到结构的近似平衡形状。在此几何位形基础上更新节点坐标, 释放预应力, 重新设定索膜结构的真实材料参数和预应力, 进行自平衡迭代求解。循环若干次, 释放掉不平衡力, 直至应力分布均匀度达到要求[3]、[4]。

由于索膜结构中索单元和膜单元均只能承受拉力, 故在单元选取时分别采用Link10和Shell41单元来模拟索单元和膜单元。在ANSYS中不能直接输入初始应力, 采用施加初始应变和降温的方法分别给索、膜施加初始预应力。

索单元的初始张力的施加:

ε=T/ (EcA)

式中:T为索的预张力;Ec取一个很小的弹性模量, 一般取真实模量的1/1000～1/5000;A为索的截面面积。

膜面预应力的施加:

先将膜面的参考温度设为0°, 则施加压力的温度荷载可用下式计算:

Δt=σEmαmh (1)Δt=σEmαmh(1)

式中:σ为膜面预张力, Ec取一个很小的弹性模量, 一般取膜材真实模量的1/1000～1/1500, αm为热膨胀系数 (可人为设定, 一般设为α=1) ;h为膜材厚度。

根据上面所述, 百富策略白菜网有限元软件ANSYS进行索膜结构找形的流程为:

(1) 选择单元类型, 一般索采用Link10, 膜采用Shell41, 设定索和膜的虚拟弹性模量和实常数值;

(2) 建立初始的平面几何模型, 进行单元网格划分, 通过初应变法及降温法为索膜单元施加初始预应力;

(3) 为各固定点施加位移约束, 将控制点一次性提高到目标高度, 打开大变形和应力刚化开关, 用几何非线性法求解, 得到初步找形结果;

(4) 更新节点坐标, 固定支座控制点, 恢复真实材料常数, 重新设定索膜的预应力状态, 关掉大变形开关, 求解;

(5) 循环若干次, 观察应力分布云图, 直至应力分布均匀度达到要求。

为了验证ANSYS找形方法的准确性, 对两个基本曲面-悬链面和马鞍面进行了分析。

根据上述原理与步骤, 用ANSYS软件对悬链面进行找形。

取内圆半径10 m, 外圆半径50 m, 厚度为1 mm的圆环膜面张拉成高度h=22.943悬链面。取1/4结构进行分析, 内外圆周固定, 内圆周提升到指定高度。膜材的弹性模量取2.55E8, 找形时取虚拟弹模为真实值的1/1000。

网格划分与找形后形状如图1a、1b所示。

从ANSYS输出的等效Von Mises应力云图 (图2) 可见, 最小应力1.88 kN, 最大应力为2.03 kN, 大部分区域为2.0 kN。同样取其中一条半径上的非约束节点进行分析, 将ANSYS输出的结果与解析解比较, 列于表1。

从表1中可见, 找形误差均在7%以内, 基本满足精度要求。分析原因, 由式 (1) 可知, 在ANSYS非线性有限元法找形过程中, 由于采用降温法来给膜面施加预应力, 导致了小弹性模量的存在, 约束了膜面的完全自由变形, 所以不可能找出完全的等应力曲面, 这就造成了与解析解间的误差比较大。由结果可见, ANSYS找形结果网格分布比较均匀。

取标准马鞍面的方程为z=2+ (x2-y2) /25, 边长为10 m, 高低点高差为4 m。找形中膜面初始预张力为2 kN, 膜材厚为1 mm。用ANSYS进行找形分析, 四条边固定, 膜材的弹性模量取2.55E8, 找形时取虚拟弹模为真实值的1/1000。

网格划分及找形结果如图3a、3b所示。经过几次不平衡力释放后, 最小应力为1.99 kN, 最大应力为2.0 kN, 等效Von Mises应力云图如4所示。在此例中, 因膜面比较平缓, 曲率变化不大, 故虽有弹性模量的存在, 应力的分布还是基本均匀的。

同样, 取两高点1、11连成的对角线上的所有节点进行分析, 与解析解的比较如表2所示。从表中可见, 两者的误差是很小的。

由以上两个算例可以得出结论, ANSYS找形方法对于曲率变化较小的情况 (标准抛物面) 有较高的计算精度, 而对于曲率变化较大的情况 (悬链面) , 计算精度较低。且由于在ANSYS找形中, 首次计算采用是大位移小弹模计算, 二次以后的平衡迭代采用真实弹模小位移计算, 弹性模量的存在在一定程度上约束了单元的自由变形, 所以找形后的网格较为均匀, 但应力情况出现了不均匀分布, 实际情况应力分布不可能是完全均匀的。可见, 上述用ANSYS软件找形的思路与步骤是正确的, 可用于进一步的工程算例分析。

本文介绍了索膜结构找形分析的有限元法, 针对索膜结构极强的几何非线性特点, 介绍了非线性有限元的基本思想。另外, 介绍了利用通用有限元软件ANSYS进行找形的方法。用悬链面及双曲抛物面这两个算例进行了验证, 证明ANSYS找形也是实用可行的。并发现ANSYS找形方法对于曲率变化较小的情况 (标准抛物面) 有较高的计算精度, 而对于曲率变化较大的情况 (悬链面) , 计算精度较低。分析原因, 在ANSYS找形过程中, 由于采用降温法来给膜面施加预应力, 引入了小弹性模量, 约束了膜面的完全自由变形, 所以不可能找出完全的等应力曲面, 这就造成了与解析解间的误差比较大。