薄膜结构以其独特的造型、绚丽的色彩和充满张力的外观成为最引人关注的空间结构之一。与传统的结构形式不同, 薄膜结构具有自重轻、覆盖跨度大、形态丰富多彩、造价相对低廉、施工速度快、结构抗震性能好、有适度透光性等特点, 可以在大跨度、大空间的公共建筑中实现人们所追求的那种自然、和谐、开敞、明亮的内部效果。

构成薄膜的柔性膜材只能承受面内拉力, 因此由面外荷载作用下产生的弯矩和剪力将通过膜材的变形转换成面内拉力。当膜材受到面外荷载作用时, 变形会产生面内压力, 由此使膜材产生局部褶皱。膜材表面的局部褶皱在整个结构上沿某个方向扩展, 导致结构本身产生局部失稳而影响结构整体的稳定性, 产生最不利的变形状态, 并影响结构的整体设计性能。特别是对于大跨度轻质的薄膜结构而言, 褶皱的出现使薄膜表面的光滑度降低, 降低了在外荷载作用下的受力性能;同时, 褶皱的出现使薄膜刚度发生变化, 改变荷载传递的路径, 尤其对薄膜结构本身的动力特性以及在风荷载作用下的气动特性产生很大影响[3]。

目前研究褶皱的方法主要有张力场理论[4?6]、可变泊松比的方法[7?10]、修正本构矩阵方法[11]、修正变形梯度的方法[12?13]和屈曲理论[14?17]等。其中, 只有基于屈曲理论的壳单元方法, 考虑了薄膜的弯曲刚度和褶皱屈曲模态, 能够得到膜材褶皱的具体信息。Suzuki等[14]采用Kármán方程和Rayleigh Ritz法模拟了不考虑弯曲刚度的圆形张拉薄膜结构, 并与考虑弯曲刚度而不含初应力的分析结果进行比较, Suzuki的工作仅限于对结构屈曲模态的分析, 而没有考虑薄膜结构的非线性特征。Rimrott等[15]考虑膜材抗弯刚度的影响, 分析了挂毯在其两个角点约束受重力作用时的褶皱情况, 并与实验值进行比较;结果表明考虑薄膜的弯曲刚度和自身重力的影响时, 薄膜会产生一系列如余弦波形的褶皱变形。Miyamura[16]将屈曲理论与有限元方法相结合分析了圆形薄膜扭转问题, 该方法精确地得到了褶皱的形状, 且与实验结果相吻合, 但是Miyamura的分析中没有考虑褶皱数量和幅度。Su等[17]利用特征值屈曲理论对一个张紧的方形Kapton薄膜的褶皱进行了分析, 并指出膜材本身的弯曲刚度是产生局部褶皱的原因。Wong等[18]基于屈曲理论, 采用有限元方法, 分析了平面膜材在面内剪切荷载以及面内拉压荷载作用下平面膜结构褶皱的形成和发展, 其结果与实验数据吻合的较好;分析中, 作者假定膜材存在小量弯曲刚度, 并考虑了初始缺陷的以及面外形变的影响。

综上所述, 由于薄膜材料存在局部抗弯刚度, 其在褶皱模拟中不可忽略。褶皱的幅值、波长和方向等具体信息对分析薄膜结构的受力性能具有重要影响。所以, 本文将采用薄壳单元分析张拉膜材和气枕式充气膜结构在几种荷载作用下的褶皱形态。同时, 由于褶皱问题的强几何非线性特征, 边界条件往往也很复杂, 不得不做较多的简化和假设。因此, 理论分析的结果具有较大的局限性, 难以得出整个变形过程的精确解, 在实际百富策略白菜网中受到较多的限制。随着有限元理论的日益完善和计算机技术的飞速发展, 采用有限元数值计算方法可以得到满足一定精度要求的数值结果。

膜材在褶皱产生过程中, 几何变形较大, 且褶皱的产生与其形成过程有关。基于该特点, 采用更新Lagrangian方法的薄壳单元对膜材在荷载作用下的变形进行描述。

单元类型选择的依据主要是对薄膜褶皱的准确描述。在褶皱计算中, 薄壳单元满足两个重要的假设:1) 沿厚度方向应力为零;2) 中性面的材料纤维在薄壳变形过程中保持直线形状, 但不一定垂直与变形后的板壳中性面。

目前百富策略白菜网较多的是基于层面和纤维理论的Hughes-Liu单元 (HL单元) [19]和基于局部坐标变换的Belytschko-Wong-Chiang单元 (BWC单元) [20]。HL单元是由八节点的六面体单元退化而成四节点的四边形薄壳单元, 其每个节点对应实体元的2个节点, 且相邻2个节点由长度不变的纤维顺序连接起来, 因此HL单元的每个结点具有5个自由度 (缺少沿厚度方向的旋转自由度) 。在厚度不变的壳元层面内, 根据不同数量的高斯积分点来决定壳元输出应力情况。BWC单元为单点积分的双线性四节点四边形壳单元, 每个节点具有全部6个自由度。BWC单元考虑了翘曲的影响, 采用映射矩阵对剪切变形的计算进行修正, 在几何方程中考虑平动自由度曲率项, 消除单元翘曲刚度带来的影响。BWC单元的运动学和动力学关系采用速率形式描述, 具有更好的数值稳定性, 因此本文采用BWC单元进行褶皱模拟。

薄膜褶皱模拟方法主要有静态隐式有限元法和动态显式有限元法两种。由于隐式算法在求解时存在收敛性问题, 同时计算速度变慢, 使得其在求解复杂张拉膜材和充气式薄膜问题时十分困难;而显式算法采用中心差分法进行动态问题的时域积分, 其不存在收敛性问题, 但对时间步长的选取有一定限制, 否则会产生数值解的非稳定。在褶皱问题分析中, 对于板壳元, 时间步ΔT满足:

式 (1) 中:ν为材料泊松比;α为时间步因子, 根据不同材料泊松比进行选取;为材料声速;Lmin为单元的最小特征长度。根据单元形式的不同, Lmin选取有几种方式:

其一, ;对三角形和四边形单元, β分别取1和0;As为单元面积;L1?4为单元的边长。另外, 当采用四边形单元时, 可采用;D1?2分别为单元对角线长度。如果采用三角形单元亦可采, min (L1, L2, L3, L4, +β1020) ]针对不同的材料泊松比选取不同的β值。

本文采用上述几种长度尺度选取标准, 基于不同的单元类型和数量, 比较得到最小时间步长作为ΔT保证计算稳定。

膜材的面外刚度很小, 其对初始缺陷非常敏感, 荷载一旦作用于膜材时褶皱便出现, 随着荷载的增加, 大量褶皱出现在膜材上。因此本文选取适当几阶模态进行线性叠加作为膜材的初始缺陷, 用以准确的计算终态褶皱形式。初始缺陷以膜材面外变形的形式施加于膜材上, 满足:

式中:Δz为面外缺陷值;wi为第i阶模态;φi为折减系数;N为选取的模态数量。

基于显式有限元理论进行静力问题分析需采用准静态 (quasi-static) 方法, 其实现方式有2种:时间缩放 (time-scaling) 和质量缩放 (mass-scaling) 。时间缩放通过减小时间步数, 降低加载速率而增大计算时间;质量缩放是指通过增加单元的质量而增大计算时间步从而减少计算时间的方法。通过增加质量来减少计算时间, 虽然会改变模型质量和单元质心位置, 但能够显著提高计算效率, 减少计算时间。在准静态分析中, 为了确保计算稳定性, 要求变形产生的动能值较小而势能值较大。因此, 本文采用计算中膜材的动能与总应变能之比小于5%, 确保惯性力引起的动力效应足够小[21], 满足工程计算要求。

由于局部作用的荷载或变形引起膜材的褶皱会在膜面内大量产生和传递;同时显式算法较小的时间步要求有利于褶皱产生与发展过程, 因此本文采用显式算法进行褶皱问题模拟。

某矩形膜材长为0.380m, 高为0.128m, 厚度为0.025×10-3m的膜材, 在有限元分析时取E=3530MPa;ν=0.31;ρ=1500kg/m3。膜材由7904个四边形单元构成, 如图1所示。为了保证计算稳定性, 将底边所有节点的自由度约束, 左右两边只保持y方向自由, 后沿y正方向施加0.5×10-3m的初始位移模拟预应力 (此处选择较大的预应力保证较少的局部模态出现) 。在预应力状态下, 将左右两边的约束解除, 将顶边除x方向以外的自由度完全约束并沿x正方向施加3×10-3m的位移荷载。根据Wong[22]的实验结果, 膜材的局部变形发生在特征值为0.2Hz附近, 所以本文选择特征值为0.2Hz附近的4阶特征模态进行叠加作为初始缺陷, 此处折减系数iφ取厚膜材厚度的12.5%, 即3.125×10-6。荷载分析前, 需沿y轴正方向施加小量的位移荷载模拟预应力, 保证计算稳定性 (本文取0.05×10-3m) , 再沿x正方向施加3×10-3m的位移, 同时对膜材左右两边界的转动自由度施加约束。

在剪切位移作用下褶皱区域分布如图1和图2所示。中间区域的褶皱与底边成45°角, 而在约束角部位置, 即左下角和右上角, 褶皱比较集中, 形成“扇形”集中区。图3 (a) 、图3 (b) 为沿y方向64×10-3m A-A截面处, 1.6×10-3m和3×10-3m剪切位移作用下的第一主应力和第二主应力值与Wong[22]的实验结果比较, 计算表明单元的第一主应力从0值开始迅速增大到稳定值 (22MPa) 附近平稳一段距离后再衰减为0;同一位置单元的第二主应力始终保持在0值附近振荡并保持较小值, 其范围在-0.25MPa?0.3MPa之间。与实验值相比, 1.6×10-3m剪切位移时截面的第一主应力值略小于实验值, 其最大误差约为10%;第二主应力值的均值与实验值非常接近。在3×10-3m剪切位移时截面的第一主应力值略与实验值十分接近, 其最大误差约为5%;第二主应力的均值与实验值接近而略大于0。

图4 (a) 、图4 (b) 为A-A截面处, 1.6×10-3m和3×10-3m剪切位移作用下的褶皱幅值和波长与实验结果比较, 计算表明单元的褶皱幅值在左右两端约束开始增大并迅速达到位移幅值最大后缓慢再增大。与实验值相比, 1.6×10-3m剪切位移时单元褶皱最大幅值与实验值一致, 如图4 (a) 所示, 在x轴50×10-3m位置褶皱最小幅值与实验结果相差50%而在其它位置的误差控制的较好, 在0?200mm位置的褶皱波长与实验结果相比错后约半个波长;而在200mm?400mm位置膜材单元的褶皱波长与实验结果吻合较好。3×10-3m剪切位移作用下单元褶皱情况与实验值相比如图4 (b) 所示, 其结论与1.6×10-3m的情况类似:在x轴50×10-3m位置褶皱最小幅值与实验结果相差50%而在其它位置的误差控制的较好;在0?50mm位置的褶皱波长与实验结果错后约半个波长, 在膜材其它部分单元的褶皱波长与实验结果吻合较好。同时在两种位移荷载工况下的褶皱数量均与实验值一致。

上述分析表明:采用壳单元模拟薄膜结构褶皱问题的方法是可行的, 通过数值分析得到褶皱分布情况、数量、方向以及应力状态均与实验结果吻合较好。

充气式薄膜结构作为薄膜结构的一种, 其依靠膜内外气压差成形来承受外荷载, 可分为气承式和气枕式两种形式。气承式是指直接向膜材所覆盖的使用空间内注入一定压力的空气使之成形;而气枕式是指向封闭气枕内充入一定压力的气体, 以形成具有一定刚度和形状的构件, 再由多个这种构件相互连接形成空间结构。本文对国家游泳中心[23]采用的气枕式充气式膜结构单元在荷载作用下的变形和褶皱情况进行研究。

国家游泳中心的气枕式充气膜结构, 其基本单元的投影平面为不规则五边形或六边形。对于气枕式充气膜结构而言, 需对封闭膜材内充入一定量的高压气体才能获得结构承载所必需的刚度和形状。因此, 计算过程首先是形态分析, 然后在通过形态分析得到结构的几何形状和相应的应力分布基础上进行荷载分析, 验证各种荷载组合作用下结构位移与膜面的应力是否在许可范围内, 确定薄膜单元是否处于松弛、褶皱的状态, 从而确保膜材既不会因变形过大而导致膜面松弛也不会因应力过大而导致破坏, 保证气枕的正常使用功能

由于受边界条件的影响, 充气式膜结构的膜内应力一般是不均匀的, 因此膜内的应力分布状态不能预先确定。同时对于建筑结构, 在其设计阶段并不能够明确给出具体空间形状, 而只能给出平面形状及某些控制点的高度 (如气枕的矢跨比一般为1/12?1/15) , 故此时其形与态均为未知, 不可同时求解。本文采用先找形后找态的分析方法确定[24]。先采用Ansys进行静力找形, 完成隐式分析得到模型初始应力态后, 将单元类型进行转换, 通过定义新的单元类型和材料, 创建显式分析所需的模型, 读入隐式求解结果在Ls-DYNA为气枕“充气”, 得到基于常温时气枕在指定工作压力下的应力状态, 完成形态分析。

国家游泳中心某气枕单元最大跨度为8.4m, 材料常数E=900MPa, ν=0.47, 厚度为0.001m采用Shell181纯膜单元模拟, 整个气枕由1454个三角形单元组成。在平面位置两层膜面的上下内表面分别施加沿外法方向不断增加的均布压力, 将气枕“吹起”。计算表明内压增至950Pa时, 得到的气枕矢高 (最高点高度) 为0.720m, 满足工程中常用的L/12?L/15 (L为气枕的最大跨度, 即8.4m) 矢高要求 (如图5所示) 。

为确定该形状在工作内压下膜内应力状态, 将通过Ansys找形得到的模型导入LS-DYNA, 将Shell181纯膜单元转换为Shell163 BWC壳单元[28], 气枕的体积为:

式中:xi为i方向坐标;nix为i方向坐标与单元法线方向的夹角。与此时体积对应的气枕内压可以由气体状态方程确定, 其与气体内能、密度有关:

式中:P为气体压力;ρ为气体密度;e为气枕内气体的内能, 其与气体体积的关系为lne= (1-k) ln V;k=cP/cv, 其中cP为定压比热, cv是定容比热。给定气枕内部空气的温度以及此状态下的定压比热与定容比热值后, 定义在该条件下的质量流量为气枕进行充气。

气枕在300°K常温条件下的绝对工作压力下 (标准大气压力105Pa与相对工作压力之和) , 假定气枕内的压力分布是完全均匀的。由质量流量确定相对工作压力约为580Pa, 此时对应气枕体积为38.842m3, 此时上下膜面的最大主应力均为3.036MPa。

为清晰展示气枕在集中荷载作用下的褶皱情况, 在施加垂直于薄膜表面向下作用的集中荷载的同时, 沿坐标轴x方向施加剪切荷载, 其作用位置如图6所示。

在上层膜外表面中央处施加向下的集中荷载200N及沿x方向的剪切荷载2k N后的变形情况如图7 (a) 、图7 (b) 所示, 上表面沿竖直方向最大位移为0.147m;下表面沿竖直方向最大位移约为0.0365m。此时内压约为740Pa, 比找形状态增加了160Pa;体积变为38.796m3, 比找形状态减少了0.046m3, 此时气枕内压与体积之积 (3908309Pa·m3) 与初始状态气枕内压与体积之积 (3902101Pa·m3) 约相等 (相差小于0.16%) 。气枕上下表面应力分布如图8 (a) , 图8 (b) 所示, 其中上表面的最大应力出现在荷载作用处的两侧, 约为7.519MPa;由于内压的增大, 作用于下表面的气压力使主应力增大, 约为4.496MPa。气枕在受到与内压反向荷载作用时, 上表面向下运动, 相对于初始形态膜面的高度有所下降, 由于剪切荷载的作用, 部分单元受压产生了褶皱 (如图9所示) ;计算得到气枕的内部压力变化很小, 所以内压对下表面的作用产生位移亦很小, 下表面的应力分布变化亦很小, 同时荷载的分布形式没有改变, 所以没有褶皱产生。

柔性材料的褶皱问题在其受力分析中无法回避, 基于屈曲理论, 采用薄壳单元模拟膜材在荷载作用下褶皱情况的研究基础还比较薄弱。通过平面张拉膜材与气枕式充气薄膜结构算例发现:

(1) BWC单元考虑了翘曲作用, 同时对剪切变形进行了修正, 完全适用于薄膜材料的褶皱分析。

(2) 通过对平面张拉薄膜的计算, 本文发现:膜材中间区域的褶皱模拟情况最理想, 但是在接近约束和自由边的位置, 褶皱模拟的结果与实验值区别较大, 同时由于初始缺陷的不确定性造成数值模拟褶皱的波长与位置还有实验结果有差别, 这是由于膜材自由边和约束位置的初始状态难于确定, 即自由边可能存在卷曲而约束位置存在残余应力的等原因造成。

(3) 膜材在平面内剪切和拉压荷载作用下的褶皱模拟与实验结果的吻合较好, 但由气枕式充气膜结构的算例可知, 考虑面外荷载作用时, 褶皱的模拟还缺乏实验结果的支持。同时考虑整个体育场的所有气枕单元实际风荷载作用下的变形情况是非常困难的, 而且气枕单元的流线型外形和光滑的薄膜材料决定了其在风荷载作用下受到的剪切作用较小而垂直于膜材表面的荷载较大, 本文仅给出剪切荷载和垂向荷载联合作用的工况, 得到对应的褶皱情况, 还须对气枕在真实风场作用下的褶皱情况进行深入分析。

通过对气枕式充气薄膜结构的褶皱模拟, 发现:由于薄膜材料质轻且薄, 气枕内空气质量与膜材本身的质量处在同一数量级, 由此产生的气弹效应不可忽略, 所以考虑气枕内气体的刚度和阻尼效应对膜材本身的受力状态以及褶皱的影响是未来工作重点。