充气膜结构[1]是由柔性膜材构成的一种充气预应力结构体系, 包括气承式和气囊式两种, 这种结构体系与张拉式膜结构不同之处在于引入预应力方式不同:张拉式膜结构通过张拉边界引入预应力;充气膜结构边界不动, 通过在膜面上施加气体压力引入预应力, 使体系具有一定形状, 产生能够承受外荷载的结构刚度。充气膜结构兴起于20世纪70年代, 主要用于大型体育场馆和展览馆等公共民用建筑, 随着新型高性能膜材的出现和制作、施工工艺的改善, 充气膜结构在21世纪又重新得到膜结构学者的重视, 百富策略白菜网范围也扩展到了工业和航空航天领域, 被百富策略白菜网于新能源领域的沼气池 (气承式) 和空间可展薄膜天线 (气囊式) 。与充气膜结构的百富策略白菜网相比, 相关研究明显不足, 目前主要集中在找形分析[2]、结构分析和剪裁分析[3]。

现阶段充气膜结构设计过程主要包括初始几何态、找形平衡态、弹性平衡态、预应力平衡态、荷载作用态和剪裁展开态等6个形态, 这6个形态分别对应找形分析、弹性化分析、预张力分析、荷载分析和裁剪展开分析等5个分析过程。目前通用软件[3?5] (如ANSYS) 将弹性平衡态作为无应力状态, 在弹性平衡态的基础之上重新施加预张力, 陈务军[6]等 (2008) 发现将真实材料属性赋予单元以后会发生形态“漂移”和预应力松弛现象, 即弹性化以后得到的预应力态与找形状态不一致。对于一般小型项目这一差异可以被忽略, 对于大型工程项目或者精度要求很高的项目就必须解决这一问题。

为解决形态“漂移”现象, 需要求解真实无应力状态, 而不是简单的将弹性平衡态当作无应力态。本文在弹性平衡态和预应力态之间增加一个真实无应力状态, 称为零应力态。从弹性平衡态到零应力态, 预应力释放, 结构刚度退化, 刚度矩阵奇异。充气膜结构零应力态还没有一种有效方法来求解, 陈务军[7]提出了一个索网结构近似零应力态求解方法。本文将充气膜结构膜面等效为膜线, 由弹性平衡态预张力和材料参数, 计算膜线无应力长度和伸长量;然后根据杆系结构平衡矩阵理论[8?10], 利用协调矩阵M-P广义逆求解节点位移, 逆向叠加求出新状态位形, 即零应力态位形。这一求解过程满足物理本构方程和几何方程, 为满足平衡方程需要采用非线性迭代算法, 从而保证预应力态、找形平衡态和弹性平衡态形态一致。本文将该方法称为非线性协调矩阵广义逆法。

本文将充气膜结构分析全过程概括为9个状态和7个分析过程。以气承式膜结构为例说明充气膜结构分析过程, 如图1。9个状态:初始几何态、找形平衡态、弹性平衡态、零应力态、预应力态、荷载作用态、构形零应力态、随机成形态和构形预应力态, 其中初始几何态、找形平衡态、弹性平衡态、零应力态、预应力态、荷载作用态为数值模拟分析态, 构形零应力态、随机成形态和构形预应力态是实际物理态, 弹性平衡态和预应力态一般为设计目标态, 通过零应力态可以将数值模拟分析态和实际物理态建立内在联系, 实现数理模型和物理模型统一。7个分析过程依次为:找形分析、弹性化分析、松弛分析、预张力分析、荷载分析、裁剪分析和充气成形分析。

找形分析是获得给定初始几何和预张力的找形平衡态, 不考虑材料本构方程, 主要的基本找形方法有力密度法、动力松弛法、小杨氏模量法。弹性化分析是将获得的找形平衡态物理化, 将膜材真实材料属性参数 (刚度EA) 赋予对应单元, 此过程其余状态参数不变, 实现找形分析的物理模型化, 得到弹性平衡态。根据物理化的零应力态, 先进行预张力分析, 获得预应力平衡态, 最后进行荷载分析。裁剪分析是由预应力态通过裁切、展开和应变补偿以后得到膜面展开形态。充气成形分析是膜片缝合完成以后充气过程中的状态分析。

由上述分析可知, 弹性平衡态已知, 零应力态和预应力态未知, 尚无法进行荷载分析, 因此, 必须先求出零应力态。从弹性平衡态到零应力态是逆问题, 称为应力释放分析或者松弛分析。表1给出充气膜结构全过程分析计算假设, 表示满足条件, 表示不满足条件, ?表示尚不确定条件。由表1可知, 松弛分析是研究难点。

充气膜结构膜面模拟方式一般分为2种:一种是将膜面划分成等应变三角形单元;另一种是将膜面划分成四边形网格[3,11], 用网格线 (膜线) 通过等截面刚度原则来代替膜面。本文采用第2种方式, 即用膜线来代替膜面。图2为辐射式膜面网格, 网格由节点和膜线组成, 粗实线为膜线, 填充区域为膜线所代表的膜面, 粗虚线为膜线代表的膜面宽度 (平均宽度) 。设膜面经向纬向弹性模量分别为Emw/ (k N·m-1) 、Emf/ (k N·m-1) , 膜线沿经纬向正交布置成网格, 则膜线经纬向线刚度/k N为:

将所有膜面单元转化成膜线, 则充气膜结构就转化成为预应力网格结构。与传统网格结构不同的是, 该网格结构膜线只承受拉力, 不能承受压力。作用于膜面上的气体压力根据从属面积原则分配到相应网格节点上, 节点荷载等于单位面积压力乘以从属面积, 方向为从属膜面法向。

充气膜结构弹性平衡态位形Xt、膜线长tL、轴向刚度EmA、张力tT和体系拓扑均为已知。假设膜线小应变、线弹性, 则膜线应力[3]:

其中:为膜线应力;Em为弹性模量;为弹性应变。

同时, 膜线应力又可写成:

其中:tt为膜线弹性平衡态张力;A为膜线等效面积。

膜线弹性应变可表示为:

式中:0l为膜线无应力长度;tl为膜线弹性平衡态长度。

通过式 (2) ~式 (4) 可得出无应力长度:

单元长度变化量为 (35) l (28) lt-l0, 将式 (5) 代入得到:

弹性平衡态一定, 膜线无应力长度0l不变, 而单元长度变化量 (35) l需要通过当前状态计算, 不断变化。

3.3零应力态求解

平衡矩阵理论以小变形假定为前提, 主要百富策略白菜网于杆系预应力网格结构自应力模态和机构位移模态求解, 从而判定体系的类型。弹性平衡态膜面预张力远大于自重产生的张力, 所以, 从弹性平衡态求解零应力态不考虑膜材自重。

充气膜结构中任一自由节点i与节点j、节点k相连, 膜线张力t、当前长度l和节点荷载fi如图3, 弹性平衡态节点平衡方程[3,10]为:

充气膜结构预应力引入机理是通过内部气压使膜面产生张力, 从而形成结构刚度。零应力态求解是要将引入的预应力去掉, 这样就必须把内部气压撤出, 使整个体系回归到无应力状态。在撤出内部气压的过程中, 所有自由节点都将产生符合变形协调条件的节点位移。按式 (7) 对所有自由节点运算, 并写成矩阵形式:

其中:N为充气膜结构体系自由节点数;b为膜线单元数。

将式 (8) 简写为:

式中:[A]为平衡矩阵;{t}为膜线内力矢量, 由膜面张力转化而来;{f}为节点荷载矢量, 由膜面气体压力转化而来。

同样, 在小变形假设条件下, 可建立充气膜结构协调方程:

式中:[B]称为协调矩阵;{u}为节点位移矢量;{∆l}为膜线轴向变形矢量。

根据虚功原理有:

将式 (11) 代入式 (10) , 则充气膜结构协调方程写成矩阵形式:

弹性平衡状态已知, 所以协调矩阵[B]已知, 当膜线长度改变量向量{∆l}求出以后, 就可利用式 (12) 求解节点位移。网格结构的协调矩阵一般为长方阵, 不可逆, 需要利用协调矩阵的MoorePenrose逆求节点位移的最小二乘解, 式 (12) 最小二乘解可写为:

式中, ([A]T) (10) 为[A]T的M-P广义逆。

则充气膜结构新状态坐标为:

零应力态对应的膜线张力应为零, 但由于数值计算并不绝对解析零, 而是一个极小值, 因此, 可以用下式表示零应力态的判定准则或计算收敛准则:

在撤出内部气压以后, 如果通过式 (14) 求得的状态不能完全满足式 (15) , 就需要在当前位形重新计算协调矩阵[B]和向量{∆l}, 然后用式 (13) 求解新节点位移向量, 采用式 (14) 生成新的状态位形, 并重新检验式 (15) 是否满足。重复上述过程, 直到满足迭代终止条件式 (15) 为止。计算流程图如图4。

本文求得的零应力态为解析零应力态, 将该求解零应力态方法定义为非线性协调矩阵广义逆法。该方法完全满足平衡方程、物理方程和几何协调方程, 避免了直接从弹性平衡态到预应力态的位形“漂移”和预应力松弛现象, 保证了结构安全和形态精确。

另外, 充气膜结构膜线在预应力状态时都处于拉伸状态, 从弹性平衡态到零应力态为逆运算, 所以式 (13) 中的向量{∆l}应全部取负值。

为验证方法的正确性, 选用单一膜线单元为例, 如图5, 单元与x轴重合, 节点1固定、节点2自由, 节点2上作用等效节点荷载f21=1.0k N。在弹性平衡态, 节点坐标分别为1 (0, 0, 0) 和2 (1, 0, 0) , 节点坐标向量{Xt}={0, 0, 0, 1, 0, 0}T, 膜线长度l12=1.0m, 预张力为t12=1.0k N, 膜线轴向刚度EmA=350k N。

撤出自由节点2上的荷载便可得到该结构零应力态。由式 (11) 得无应力长度l0=0.9972m, 单元长度改变量∆l=0.0028m, 协调矩阵B=AT=[1, 0, 0], 协调矩阵M-P逆B+=[1, 0, 0]T, 根据式 (14) 求得零应力态节点坐标向量{X0}={0, 0, 0, 0.9972, 0, 0}T, 节点2产生x向位移-0.0028m, 零应力态索长为0.9972m, 与无应力长度相同。此时单元张力为零, 此状态即为撤除节点荷载以后的理想零应力态从而证明了本文求解零应力态方法的正确性。

半球形气承式膜结构弹性平衡态如图6, 整体坐标系圆心与球心重合, 周边节点全部固定, 半径R=3m, 体积按公式V= (2/3) πR3计算, 找形分析保持体积为56.55m3不变, 气体压力为0.69k Pa。

弹性化时采用正交异性膜材, 经向弹性模量为500k N/m, 纬向弹性模量为300k N/m, 弹性平衡态网格大小不同, 所以各膜线刚度不同, 张力也不同。

将0.69k Pa内压全部撤出, 各膜线从有应力长度收缩到无应力长度, 自由节点产生符合协调变形条件的位移, 膜线张力变为零, 得到零应力态。为保证解析零应力态尽可能接近理想零应力态, 取迭代终止条件满足:经过6次迭代计算以后, 最大膜线张力为2.14×10-13k N, 近似为0, 此状态可视为理想零应力态, 证明了算法的高效性。为比较系统从弹性平衡态到零应力态形态变化, 将节点 (编号1~编号9) 坐标及相应位移列出 (见表2) ;将膜线 (编号 (1) ~编号 (8) ) 两种平衡形态张力在表3中列出。

从位形变化来看, 从弹性平衡态到零应力态, 越靠近x-y平面, 节点位移越大, 最大位移绝对值为0.059m。零应力态单元长度与无应力长度在表4中列出, 保留到小数点后4位时两者无差异。

气枕形充气膜结构弹性平衡态如图7所示, 气枕平面为矩形, 长15m, 宽9m。整体坐标系原心与几何形心重合, 周边节点全部固定, 找形分析保持内压为0.3k Pa不变。弹性化时采用正交异性膜材, 经向弹性模量为500k N/m, 纬向弹性模量为300k N/m。

将0.3k Pa内压全部去除, 迭代终止条件满足:||{t}||≤0.001N和||[A]{t}||≤0.001N, 经过6次迭代求得解析零应力态最大膜线张力为5.3×10-9k N。将气枕上半部节点 (编号1~编号14) 坐标及位移在表5列出;将膜线 (编号 (1) ~编号 (8) ) 两种平衡形态张力在表6中列出。从位形变化来看, 最大位移绝对值为0.169m。

充气膜结构分为气承式和气囊式两种, 本文将充气膜结构分析全过程分为9个状态和7个分析过程, 通用软件对充气膜结构分析, 没有求解零应力态, 而是采用弹性平衡态作为近似零应力态, 这样会产生形态“漂移”和预应力松弛现象。由于柔性预应力膜结构形态复杂、分析理论难, 认识不够深入等原因, 目前还没有有效求解零应力态方法。在零应力态求解过程中, 因弹性平衡态已经是稳定张力平衡形态, 单元符合线弹性、小应变假定。本文用膜线单元模拟膜面, 将膜结构转化成为网格结构, 然后利用杆系结构平衡矩阵理论, 提出了一种求解零应力态的方法, 称为非线性线性协调矩阵广义逆法, 并编制了相应求解程序。算例分析表明, 本文算法正确、高效, 从而解决了一个充气膜结构分析的理论难题。