膜结构兴起于20世纪初, 是现代柔性空间结构的典型代表[1], 具有造型优美、通透性好和施工便利等优点, 工程百富策略白菜网日益广泛[2].但实际工程中, 由于膜结构形状复杂, 设计不合理、膜面初始缺陷等因素常会导致膜面局部排水不畅、产生并加剧袋状效应, 最终导致膜结构积水破坏.近年来全球灾害天气增多, 袋状效应导致的膜结构破坏事故越来越多[3], 给社会生产和人民生命安全造成重大损失, 因此有必要对膜结构的袋状效应进行深入研究.

国内外对膜结构的研究主要集中在找形分析、受力分析和裁剪分析等方面[4].近年来建筑膜结构研究热点主要有最小曲面[5]、考虑流固耦合的风振分析[6]和长期荷载作用下膜结构性能[7]等.工程实例表明, 强降雨是造成膜结构破坏的重要因素, 有关膜结构在强降雨作用下的袋状效应研究尚处于探索阶段.2008年, 罗永赤[3]对膜面形态进行判断, 给出了防止膜结构积水的合理控制条件;2012年, 吴明儿[8]对ETFE膜片进行了积水极限承载力研究.2016年, Brown等[9]根据水力学理论创建了用于模拟任意形状三维膜面上积水分布的计算工具, 结合传统有限元软件, 进行了膜结构在积水作用下的动态分析.

本文以充气膜为例, 采用向量式有限元法对强降雨作用下具有初始缺陷的膜结构产生袋状效应的过程进行研究, 探讨了基于向量式有限元的膜单元内力公式、膜结构初始缺陷施加方式、强降雨模拟方法以及膜结构破坏准则等问题, 建立了膜结构积水袋状效应的分析流程, 并通过自编MATLAB程序, 对一具有初始缺陷的正方形充气膜结构在强降雨作用下从袋状效应形成直至结构破坏的全过程进行模拟, 验证了本文方法的正确性和有效性.

向量式有限元 (Vector Form Intrinsic Finite Element, VFIFE) [10,11,12]是一种基于空间和时间点值描述、广义向量力学理论的新型有限元分析方法.该方法将构件离散成有质量的质点和质点间无质量的单元, 通过引入虚拟的逆向运动, 扣除单元的刚体位移得到单元纯变形位移, 再由材料本构关系求解单元节点内力, 得到质点不平衡力, 根据牛顿第二定律, 运用中央差分法求解质点位移, 从而得到整个体系的行为.向量式有限元法质点位移公式为:

式中, xin+1表示第n+1时间步初时刻质点位移, xin、 为第n时间步初时刻质点位移和速度, mi为质点i的质量, Fin为质点受到的不平衡力, Fin=fnext_i+fnint_i, fnext_i为质点外力, fnint_i为单元传给质点的内力, h为时间步长, ζ为阻尼因子, 

三角形常应变薄膜单元 (Constant Strain Triangle, CST) 是最常用的薄膜单元[13], 单元面内应力应变为常量分布.向量式有限元CST薄膜单元基本公式的推导流程主要包括单元节点纯变形位移计算、单元节点内力计算、质点合力集成和质点位移计算[14].

三角形CST薄膜单元的变形和内力满足虚功方程:

式中 是变形坐标系下节点i的节点位移和节点内力向量.单元节点位移向量可记为 , 在变形系下有 , 可简写为

引入传统有限元形函数 , 可得单元上任意位置点的位移向量 为:

由膜单元的变形位移分布, 可得应变分布向量:

式中 , Bm是应变-位移关系矩阵, 形同传统有限元, Bm*是由Bm去掉 对应的行和列后得到的简化矩阵.

假设材料为线弹性, 本构矩阵为D, 可得单元应力分布为:

根据应力应变关系, 可计算膜单元的变形虚功:

式中, 是单元的初始应力矩阵.

根据虚功方程, 可求得节点内力向量:

式中, , 其余节点内力分量 可由静力平衡条件求得.

在向量式有限元计算中, 根据主应力-主应变准则判断膜材是否出现褶皱, 褶皱处进行本构矩阵修正, 实现褶皱的检测和模拟[14].膜面的松弛和褶皱判定采用主应力-主应变准则[15].设膜单元主应力σ1≥σ2, 主应变ε1≥ε2:当σ2>0时, 膜单元为纯拉状态, 膜面不产生褶皱;当ε1≤0时, 膜单元为松弛状态, 膜面产生双向褶皱;除上述两条以外的其他情况, 膜单元处于单向拉伸状态, 膜面产生单向褶皱.在向量式有限元中, 可根据张力场理论采用修正本构关系矩阵的方法进行褶皱的处理[14].具体处理方法详见文献[14].

气承式充气膜的找形可参考传统有限元的找形方法, 即求结构最终平衡态的过程[16].采用向量式有限元法进行充气膜结构找形的流程如下:

(1) 输入节点和单元信息, 建立充气膜平面投影的几何模型;

(2) 设定时间步长、阻尼比和材料属性, 施加边界约束;

(3) 在初始时间步令膜单元的初始应力矩阵, 其中σ0为初始预应力.施加气压荷载p0, 求解质点位移, 得到膜结构形态和应力分布;

(4) 以上一步的膜结构形态和应力分布为计算构型, 施加气压荷载p0, 求解质点位移, 得到膜结构形态和应力分布;

(5) 重复步骤 (4) 直至膜结构应力分布水平满足目标要求.

由于膜结构设计不合理或受局部荷载作用, 膜面会出现初始凹陷, 这类初始缺陷使得积水或积雪荷载进一步堆积, 从而形成袋状效应.本文通过在膜结构顶部施加集中力 (近顶部区域施加均布荷载) 的方式, 模拟结构顶部微小凹陷.

实际工程中一般采用每小时降水量用于结构计算[17].表1为文献[17]给出的降雨等级分类表.其中, 强降雨指降雨强度I≥64mm/h的大暴雨.袋状效应模拟中, 先将每小时降雨深度换算成每一时间步长的降雨深度, 再求汇集到凹陷内的积水深度, 可得相应水压力荷载.

膜结构设计不合理、膜面存在初始缺陷常常会导致局部排水不畅, 形成局部凹陷.降雨过程中, 水往凹陷处堆积, 局部凹陷范围和深度逐渐增大, 进一步加剧了袋状效应.在强降雨的持续作用下, 汇集到凹陷内的雨水越来越多, 膜面变形和应力增大速度越来越快, 短时间内可因膜面撕裂而破坏.本文以具有初始缺陷的充气膜为计算模型, 对袋状效应的发展过程进行数值模拟, 基本假定如下:

(1) 膜材处于线弹性阶段;

(2) 假设膜材密度恒定且不计膜材自重;

(3) 降雨过程中, 雨强保持不变且雨滴均匀下落;

(4) 凹陷内降落的雨水全部沿膜面滑落至底部汇集, 滑落时间忽略不计;

(5) 凹陷外降落的雨水全部滑落至地面, 不参与计算.

积水荷载模拟的主要工作为确定每一时间步膜面积水半径和深度.图1和图2分别为第n时间步和第n+1时间步积水荷载沿y=0的剖面示意图.

图1曲线DA1B1C1B′1A′1D′为第n时间步初时刻的膜面形状, A1 (A′1) 为凹陷的边缘点, C1为凹陷的最低点, B1 (B′1) 为积水平面与膜面的交点.R1、r1、h1分别为凹陷半径、积水半径、积水深度.图2曲线DA2B2C2B′2A′2D′为第n+1时间步初时刻的膜面形状, 各符号含义同上.充气膜第n时间步末时刻的膜面形状, 为第n+1时间步初时刻的膜面形状.假设微小的膜面变形不引起积水深度变化, 则第n+1时间步的积水深度h2为h1+Δh, 其中Δh为由降雨增加的积水深度.Δt时间内, 凹陷内降落的雨水全部汇集至凹陷底部, 近似形成以r1为底面半径的圆柱体, 则Δh为该圆柱体的高度.

设降雨强度为I, 新增积水深度为:

则第n+1时间步的积水深度:

设积水平面与膜面相交处单元倾角为θ, 则第n+1时间步的积水半径:

对积水平面以下、积水半径内的膜单元施加积水荷载pw (静水压强) :

式中, ρ为水的密度, g为重力加速度, hm为膜单元形心与积水平面的高差.

膜材的破坏是膜结构分析中的一个重要且复杂的问题, 包括几何非线性、材料非线性和接触界面非线性等非线性特征[18].本文?用Mises应力达到失效应力值作为膜材断裂失效的判断准则[19], 即当Mises应力达到容许值时认为膜材已破坏并终止运算.

吴明儿等[8]验证了积水荷载作用下膜面可达到稳定状态.本文所设定的稳定状态判定标准为:膜面变形稳定即凹陷最低点的位置基本保持不变, 同时积水至膜面最高点.

采用向量式有限元方法, 模拟强降雨过程中, 膜面因积水荷载不断下凹, 积水量随之逐渐增大的动态过程.首先进行充气膜的找形, 并在顶部施加集中荷载, 形成初始凹陷.然后施加积水荷载, 第一步假设微小凹陷处积满水, 其他时间步根据膜面形态, 按照4.2节所述方法确定积水荷载, 进而求出该时间步末时刻的膜面形态, 作为下一时间步的初始形态, 直至膜面破坏或稳定.膜结构袋状效应模拟分析流程见图3.

选用一正方形气承式充气膜结构, 模拟强降雨作用下的袋状效应.正方形边长L=10m, 起始平面状态和网格划分情况如图4所示, 边界条件为刚性固定.膜材选用ETFE材料, 膜厚t=0.2mm, 弹性模量Ex=Ey=703.67MPa, 泊松比为0.42.设定极限应力等于膜材的第一屈服强度σ=16.3MPa.在计算范围内, 膜材为线弹性各向同性材料[20].

充气膜成形时, 先在初始平面状态下施加0.3MPa初应力作为结构的初始内力.然后施加150Pa的内压作为外荷载, 运用向量式有限元法求得结构最终的平衡状态.如图5所示, 充气膜成形后矢高f=0.594m.

采用在顶部区域施加集中力的方式形成初始凹陷.为避免应力集中, 实际计算中将集中力分散为小范围内的均布力.对距充气膜中心点0.268m范围内的膜单元施加竖直向下、共计300N的集中力.经计算, 初始凹陷深度为0.028m, 初始凹陷半径为0.870m.初始凹陷产生前、后的示意图见图6、图7, 其中z向坐标值放大5倍.

本例降雨强度取I=100mm/h, 等级为大暴雨.在充气膜形成初始凹陷后, 第一时间步假设积满水, 计算充气膜的变形和应力.之后根据降雨强度和凹陷情况, 确定每一时间步积水的范围和深度, 并计算充气膜的变形和应力.当降雨时长约32min时, 膜单元最大等效应力达到设定的极限应力, 结构破坏.

图8为降雨不同时刻膜面等效应力图.从图中可以看出初始凹陷很小, 基本不影响膜面应力的均匀分布;降雨16min后, 膜面凹陷范围和应力略有变化;降雨后期凹陷范围和应力迅速增大.应力呈袋状分布, 中心区域应力逐渐增大直至达到极限应力.

图9为降雨不同时刻充气膜1-1剖面上各节点z向坐标图, 横坐标表示不同节点距离中心点的距离, 纵坐标表示节点的竖向坐标.从图上可以直观看出:随着时间的推移, 袋状变形的范围和深度不断增大, 且增大速度越来越快, 直至结构破坏.

图10为降雨不同时刻充气膜中心节点O的z向坐标图.可以看出中心节点在积水作用下向下移动的速度越来越快, 结构非线性特征明显.临近破坏时, 中心节点迅速下移, 也反映了袋状效应后期导致的破坏速度之快.

图11和图12是袋状凹陷不同时刻的深度和半径变化图.可以看出, 袋状效应发生后, 随着积水不断汇入, 膜面凹陷的深度和半径逐渐增大.图11显示, 凹陷深度增大速度随时间越来越快.图12中0~8min、12min~20min和24min~28min凹陷半径基本不变, 仅有凹陷深度的增加, 说明膜面凹陷范围增大过程中有一定的平稳期;8min~12min、20min~24min和28min~32min凹陷半径快速增大, 说明膜面凹陷范围增大过程中有一定的快速增长期;平稳期和快速增长期交替出现.

图13和图14是不同时刻袋内的积水深度和半径变化图.可以看出, 袋状效应发生后, 降落的雨水不断汇入膜面凹陷内, 积水深度和半径逐渐增大, 且增大速度不断加快.

本文研究了膜结构积水袋状效应的基本假定、数值模拟流程及计算方法, 采用向量式有限元法对具有初始缺陷的充气膜结构进行了强降雨作用下从袋状效应产生到结构破坏的全过程模拟, 得到如下结论:

(1) 采用向量式有限元法分析膜结构的袋状效应全过程是有效的;

(2) 具有初始凹陷的充气膜在强降雨作用下, 凹陷形如袋状, 凹陷的范围和深度随时间增加会迅速增大, 结构由于袋状效应最终破坏;

(3) 为避免强降雨作用下膜结构产生袋状效应而破坏, 膜结构设计和使用中应尽量避免局部凹陷.如已出现局部凹陷, 应及时采取人工排水等措施避免袋状效应的进一步加剧.