虽然ETFE薄膜已经发明30余年, 但主要用于医药卫生行业.而在建筑领域虽有近20年百富策略白菜网历史, 只是1995年后被成功百富策略白菜网于气囊膜, 才获得快速发展和广泛关注.地标性建筑有:英国为迎接新世纪而建的“伊甸园 (Eden of Garden) ”、德国为2006年世界杯而建的“安联运动场 (Allianz Arena) ”、以及我国为2008奥运会而兴建的“水立方”和“鸟巢”等[1], 之外在欧洲还有众多ETFE膜建筑.

除却ETFE气囊膜制造技术原因外, 国内外仅见关于ETFE薄膜基本特性和百富策略白菜网的综述文献[2].但无基于ETFE薄膜的气囊膜设计原理、数值分析方法及结构特性方面的研究.文[3]利用Ansys软件分析了“水立方”中典型气囊膜结构特性.

本文首先介绍气囊膜设计分析方法, 随后介绍基于约束条件力密度法的气囊膜非线性数值分析算法以及EasyVol数值分析功能.重点研究了气囊膜单元在风吸、风压载荷作用下, 不同约束条件对数值分析结果的影响及分析过程的数值特征变化.

薄膜自身无定形, 随预应力作用而产生稳定平衡形态.张拉膜由机械张拉产生预应力 (mechanically stressed) , 充气膜则由充气产生预应力 (pneumatic stressed) .预应力起到膜成形、维持形态、抵抗外载荷和作用, 是膜抗力体系重要因素, 而不是单纯载荷或作用.文[4]规定预应力为永久荷载.文[5,6]规定膜预应力作为永久荷载进行膜设计, 但未对荷载组合系数、抗力分项系数作出合理规定.对较刚性结构体系, 如预应力网壳, 预应力等同为普通永久荷载影响较小.对形态完全由预应力决定的轻结构体系, 预应力仅等同为永久荷载进行设计有待商榷和深入研究, 特别是封闭建筑空间的充气膜 (气囊膜或气承膜) .ASCE17-96采用分项系数的极限状态设计方法, 规定了气承膜内气压与永久荷载不同的组合系数, 便是充分认识到气压既是抗力因素又是永久荷载, 而不同于一般永久荷载.ASCE17-96规定基本荷载组合包括[7]:

(a) 1.4D-1.0P0;

(b) 0.9D-1.6P0;

(c) 1.2D+ (0.2Lr或1.6Sl或1.6Rl) -1.0P0;

(d) 1.2D+1.3W+ (0.1Lr或0.5Sl或0.5Rl) -1.0P0;

(e) 1.2D-1.3W+ (0.1Lr或0.5Sl或0.5Rl) -1.6P0;

(f) 0.9D-1.3W-1.1Pm;

(g) 1.2D+1.0E-1.0P0;

(h) 0.9D-1.0E-1.6P0.

其中, D为永久荷载, P0为正常工作气压 (一般0.2kPa～0.5kPa, 常0.3kPa左右) , Pm最大工作气压, Sl雪荷载, Lr活荷载, Rl雨荷载, W风荷载, E地震作用.Rl雨荷载 (刚性屋面体系一般不考虑) 常被忽略, 但对膜设计却应予以重视.

气囊膜可按下列设计表达式进行设计:

γ0S(0P+P0)≤R(0P+P0) (1)γ0S(0Ρ+Ρ0)≤R(0Ρ+Ρ0)(1)

式中, γ0结构重要性系数, S(0P+P0)S(0Ρ+Ρ0)为将P0作为荷载的荷载效应组合设计值 (可参考ASCE17−96[7]),R(0P+P0)17-96[7]),R(0Ρ+Ρ0)为与气压成函数关系的膜抗力设计值.

式 (1) 中0P为大气压, 气囊内气压值Pi=0P+P0, 气囊外为Pe=0P, 内外气压差ΔP=P0, 正是由于这个气压差给膜“成形, 并产生抗力”.标准大气压值为101.325kN/m2, 大气压0P随所在位置地球纬度、高度、温度等有所变化, 常取100kPa[1,8].

在荷载作用的“瞬时状态”, 假设气囊膜内气体质量不变, 即不与外界发生气体和热交换, 温度不变.此时, 气体满足气体热力学玻伊耳定律 (Boyle's Law) , 有

PiVT=constΡiVΤ=const1PiV=const2 (2)

式中, 气囊内气压Pi=0P+P0, V为气囊内气体体积, T为气体温度.

载荷分析时, 由于气囊发生大位移, 气囊体积发生改变, 相应内气压产生变化, 从而导致气囊抗力水平改变.因此, 在结构非线性分析过程的每一个循环步内, 还应满足

P1iV1=P2iV2=const2 (3)

并进行双重循环迭代计算, 从而求出稳定收敛解.否则, 求出的内压、膜应力存在误差, 以至于结论不正确[3].

当气囊膜使用过程随季节温度变化, 大气压Pe=0P会出现较小变化 (可忽略不计) , 但由于温度改变导致内气压变化, 进而显著影响气囊膜结构特性.

P1iV1T1=P2iV2T2=(P1i+ΔP)(V1+ΔV)T1+ΔT (4)Ρi1V1Τ1=Ρi2V2Τ2=(Ρi1+ΔΡ)(V1+ΔV)Τ1+ΔΤ(4)

式中, T1、T2 (°K) 为温度.

如年温差ΔT=±20°K, 取T1=288.15°K (15℃) , 0P≈100kPa, P0=0.3kN/m2, 并假设体积变化忽略 (实际气囊体积有较小被动变化) , 则气压变化为ΔP=±6.962kPa, 远超过一般气囊膜内气压最大值和最小值, 因此, 需要控制设备主动对气囊内气压调节, 泵出或泵入气体以维持相对稳定的正常工作气压.

基于连续介质力学的薄膜线弹性理论, 文献[9]给出了详实的气囊膜形态 (任意曲面、球面、柱面、截锥面、组合曲面等) 在对称雪载或非对称雪载下变形、内力解析公式, 这是气囊膜分析理论的经典文献.

气囊膜在载荷作用下会产生大位移, 出现强几何非线性, 线性分析理论不满足工程设计要求.张拉膜分析理论文献多, 但充气膜分析理论较少.气压力始终法向作用于膜面, 随气囊形态改变, 法向量也改变, 这是非保守力.根据非线性连续介质力学理论, 文[10]考虑气囊内空气为不可压缩介质, 建立了统一的大位移有限元列式和增量迭代非线性分析方法.文[11]采用动能阻尼 (kinetic damping) 动态松弛法由增量列式增量迭代求解气囊膜.文[12]首次建立了张拉膜找形与载荷非线性分析的一般理论“力密度法”, 然后提出了充气膜分析理论“受约束条件的力密度法”, 并经过长期研究开发成EasyVol[13].

基于非线性找形分析理论基础, Linkwitz和Schek (1971) 首先提出了力密度法 (Force Density Method) , 随后经Schek、Linkwitz、Argyris、Grundig等 (1974) 逐渐发展完善, 至今仍然是欧洲特别是德国最流行的索网和张拉膜结构找形分析方法.膜结构分析力密度法的第一步是采用膜线单元模拟膜面, 然后对膜线网格体系进行分析, 最后等代求膜面应力.

设等代的膜线网中自由节点pi(xi‚yi‚zi)pi(xi‚yi‚zi), 边界约束节点(xfi‚yfi‚zfi)(xfi‚yfi‚zfi), 分别构成n维向量Xn、Yn、Zn和nf维向量Xf、Yf、Zf.且设膜线长lj和膜线轴力sj构成m维向量l、s.膜线网中两相邻节点间坐标差向量u、v、w可以用杆元—节点矩阵[C]=[CnCf][C]=[CnCf]表示为[12,14]

u=CnXn+CfXf, v=CnYn+CfYf, w=CnZn+CfZf (5)

式中, 杆元—节点矩阵[C]=[CnCf][C]=[CnCf]定义为

C(j,b)=⎧⎩⎨⎪⎪+1 i(j)=b−1 k(j)=b0 其他C(j,b)={+1i(j)=b-1k(j)=b0其他

杆元j有节点i、k, 则[C]为+1、-1构成的m×(n+nf)m×(n+nf)稀疏矩阵.

设U、V、W、L分别为向量u、v、w、l构成的对角阵, 膜线网体系节点平衡方程可写为

CtnUL-1s=px, CtnVL-1s=py, CtnWL-1s=pz (6)

记q=L-1s, qj=sj/lj表示膜线j单元的力密度, 因此, 平衡方程 (6) 式可改写为

CtnUq=px, CtnVq=py, CtnWq=pz (7)

记Q=diag(q)(q), 则Uq=Qu、Vq=Qv、Wq=Qw, 代替式 (7) 对应量, 然后将式 (5) 代入式 (7) , 并整理可得

CtnQCnx+CtnQCfxf=pxCtnQCny+CtnQCfyf=

pyCtnQCnz+CtnQCfzf=pz (8)

对膜线网体系, 根据初始条件 (初始位置坐标、力密度、节点荷载) , 可由 (8) 式求出体系平衡状态位形.无外荷载时, 根据力密度得到对应平衡形态, 就是“找形分析”.大型膜结构体系, (8) 式中系数矩阵常为大型稀疏矩阵, 因此, 采用直接解法效率低.Easy软件采用稀疏矩阵压缩存储, 并由共轭梯度迭代法求解[14].

对气囊膜 (充气膜) 仅根据初始力密度和初始基本几何条件还不能求出唯一稳定平衡形态, 还需要设定合理气压、封闭体积、最小曲面表面积等.

假设充气膜封闭体积为常数, 即V(x‚y‚z)=const=aV(x‚y‚z)=const=a.又因气囊膜位形与力密度存在对应关系, 位形坐标可表示为力密度的函数, 即x=x(q)x=x(q)、y=y(q)y=y(q)、z=z(q)z=z(q).因此, 气囊膜封闭体积为常数可表示为V∗(q)=V(x(q)‚y(q)‚z(q))=const=a.V*(q)=V(x(q)‚y(q)‚z(q))=const=a.

当采用一般力密度第一次找形分析后, 得到一个形态以及对应封闭体积, 但此体积不满足给定体积约束, 需求新位形.设初始力密度q(0)q(0), 对应体积V(q(0))≠aV(q(0))≠a, 则求q(1)q(1)使其对应稳定形态体积V∗(q(1))=aV*(q(1))=a.假使q(1)=q(0)+Δqq(1)=q(0)+Δq, 则计算出Δq便可满足V∗(q(1))=a.V*(q(1))=a.

对V∗(q(1))=V∗(q(0)+Δq)=aV*(q(1))=V*(q(0)+Δq)=a在q(0)q(0)位形按照Taylor级数展开, 并取线性项, 可得到

V∗(q(0))+∂V∗(q(0))∂qΔq=a (9)V*(q(0))+∂V*(q(0))∂qΔq=a(9)

记Gt=∂V∗(q(0))∂q,r=−V∗(q(0))Gt=∂V*(q(0))∂q,r=-V*(q(0)), 则 (9) 式可简写为

GtΔq=r+a (10)

式 (10) 为欠定方程, 可求最小范数最小二乘解, 即使ΔqtΔq+atPa→Min, P为加权对角占优阵.由此可解得

Δq=G(GtG+P−1)−1r (11)Δq=G(GtG+Ρ-1)-1r(11)

由式 (11) 可求出力密度增量, 进而得到新的形态.判断新形态是否满足约束条件, 否则继续迭代分析, 直到满足收敛准则.

对给定约束条件, 气压为常数, 即P0(x(q)‚y(q)‚z(q))=const=p0Ρ0(x(q)‚y(q)‚z(q))=const=p0.可按照相似的方法, 建立基于气压为目标函数的约束条件, 并计算相应平衡形态.气压始终法向作用于膜面, 膜面形态改变, 对应膜面法向量也变化, 气压法向应随之修正.

充气膜在找形分析时约束条件可以为体积常数或气压常数, 而载荷分析时则常为封闭体积内无气体 (定质量) 与热交换 (定温度) , 即PV=Const, 单纯体积、气压约束只能理解为特殊载荷工况.

EasyVol模块是Easy软件中专门用于充气膜结构分析的模块, 它基于EasyForm、EasySan模块和基本集成环境, 自身还包括两个模块:找形 (VolFin) 和载荷分析 (VolStan) .EasyVol可进行气囊膜和气承膜找形、载荷分析, 以及与梁、杆单元组合空间网格体系的结构分析.本节和上节算法功能均由EasyVol实现.

ETFE薄膜总体上有一个塑性较小段, 强度较低, 反复载荷下仍呈弹塑性, 因此ETFE适宜内力较小、动力效应较小的结构形式.单元化、网格化是ETFE气囊膜设计主要思想, 而菱形气囊膜单元具有较好形态稳定与承载能力, 以及灵活的几何设计, 如“安联”体育场用菱形气囊单元[1,14].下面详细研究菱形气囊膜单元初应力态、风压、风吸载荷时, 不同约束条件的数值分析特征与结构特性.

图1为菱形气囊膜单元, 对角长12m×6m, 初应力1.0kN/m, 膜线网格0.2×0.2, 约束条件:V=25.000m3.找形得到平衡形态, 气囊矢高0.7313m, 上下对称, 气囊厚1.4626m, 体积24.979067m3, 内气压0.318135kPa, PV=2505.853434kN·m, 体积与约束达到收敛标准.当约束条件为P0=0.318135kPa, 平衡态体积24.979067m3、PV=2502.670802kN·m;当约束条件为PV=2505.853434kN·m, 平衡态体积24.947342m3.可以发现:采用任意协调约束条件, 可得到一致的找形结果, 差异缘于数值分析本身.但需特别指出:PV=const必须是包含大气压, 而不仅内气压, 即P=0P+P0.

取ETFE-200NJ薄膜, 弹性刚度25kN/m[14].无外荷载下仅内气压、自应力态分析, 所得形态与找形结果一致, 双向应力均1kN/m, 且三种约束条件对计算结果无影响.

基于初应力态分析结果, 在上表面加风压0.55kPa.表1为风压荷载作用下数值分析过程采用不同约束条件时迭代过程数值特征.V=const和PV=(0P+P0)V=constΡV=(0Ρ+Ρ0)V=const约束所得结果差异较小, 前者位形略高, 上、下表面矢高差分别为7.3mm、2.3mm, 气囊膜最高点纵横剖面如图2.上表面下降挠度为255.7mm、263.1mm, 下表面下降挠度271.8mm、274.0mm, y向挠跨比4.26%～4.57%.两者应力基本相同, 上表面x向最大应力0.7kN/m, y向0.4kN/m;下表面x向最大应力1.3kN/m, y向1.9kN/m.内气压均增大, 由0.318135kPa分别增加为0.663414kPa、0.660385kPa.当约束条件P=(0P+P0)=constΡ=(0Ρ+Ρ0)=const, 气囊被压凹 (如图2) 为悬吊状, 形态失稳.虽然此时气囊仍具有承载能力, 但已失结构设计意义.上下表面矢高-709.5mm、-753.6mm, 挠度分别为1440.8mm、22.2mm, 气囊仅43.9mm.上表面最大应力0.6kN/m, 下表面x、y向最大应力分别为1.0kN/m、1.1kN/m.

基于初应力态分析结果, 在上表面加风吸0.55kPa.表2为风吸荷载作用下数值分析过程采用不同约束条件时迭代过程数值特征.V=const和PV=(0P+P0)V=constΡV=(0Ρ+Ρ0)V=const约束所得结果相同, 气囊膜最高点纵横剖面如图3.上表面被吸上升265.1mm, 下表面上升242.7mm, 气囊厚度增22.4mm, 内气压降低较大-0.21014kPa, 仅0.107995kPa.气压降低将显著减小气囊形态稳定性, 容易晃荡.上表面x向最大应力1.3kN/m, y向1.8kN/m;下表面x向最大应力0.8kN/m, y向0.5kN/m.当约束条件P=(0P+P0)=constΡ=(0Ρ+Ρ0)=const, 气囊被吸最高, 上表面上升388.0mm, y轴挠跨比6.47%, 下表面下降22.2mm, 体积增6.747m3 (27.01%) , 气囊增厚402.2mm (27.5%) , 原厚1462.6mm.上表面x向最大应力1.5kN/m, y向2.3kN/m;下表面x向最大应力1.0kN/m, y向1.1kN/m.约束条件V=const、P= (0P+P0) =const表明气囊为开口 (open) , 由增加、减少内部气体达到计算平衡, 这可以理解为使用动态过程.而PV= (0P+P0) V=const则表示气囊封闭, 由气压、体积、荷载作用时形成平衡.气囊膜荷载非线性分析时必须考虑内部气体作为抗力对形态平衡迭代过程分析的约束条件, V=const、PV=(0P+P0)V=constΡV=(0Ρ+Ρ0)V=const两个约束条件分析结果相近或差异极小.P=(0P+P0)=constΡ=(0Ρ+Ρ0)=const分析过程内气压值不变, 这是仅将内气压作为荷载进行分析时而不自觉采用的约束条件, 由此可见分析结果的巨大差异.

根据ASCE17-96规定, 气囊膜设计分析应充分考虑内气压是特殊抗力因素和长期载荷因素, 在设计时内气压与一般永久荷载采用不同的荷载组合系数.在计算分析过程应基于内气压为抗力因素, 进行非线性分析过程必须考虑内气压变化符合气体玻伊耳定律及对平衡形态和数值计算收敛的影响, 任何仅将内气压单纯作为荷载 (即P=(0P+P0)=constΡ=(0Ρ+Ρ0)=const) 而无其他约束条件进行迭代非线性分析都不合理, 所得分析结果不正确.

致谢:感谢Dr.Dieter Strobel、Dr.Peter Singer对分析研究过程指导和帮助, 以及德国Hi-Technet公司提供用于科研的Easy软件.