薄膜结构因重量轻、几何形状复杂、自振频率较一般结构低等特点, 导致风荷载成为此类结构设计的控制荷载之一。国内外许多文献[1,2]例证了此类结构完全或者部分倒塌的主要原因之一在于:对荷载大小和结构反应的错误估计, 尤其是风荷载的作用。由于风荷载具有动力特性, 导致结构极易发生大幅振动, 响应复杂, 因此, 在薄膜结构荷载分析时, 通常采用的静力风荷载计算无法完全正确地反映动力影响, 有必要对薄膜结构进行风致响应动力分析。

到目前为止, 虽然非线性动力时程分析理论已经比较成熟, 但在薄膜结构动力响应分析方面的工作还很少。在薄膜结构风致效应分析时, 需要综合考虑大变形引起的几何非线性和膜材褶皱引起的物理非线性的影响。非线性动力时程分析中, 常采用直接积分法, 如Wilson-θ和Newmark等, 本文采用稳定性和精确性较好的Newmark法。在求解动力方程时, 考虑到薄膜结构平面外刚度小, 在荷载作用下产生大位移的几何非线性问题, 本文用收敛性较好的全Newton-Raphson迭代法进行迭代计算。另外, 当薄膜结构位移较大时, 易发生膜材的褶皱, 处理褶皱的方法主要有张力场理论、可变泊松比的方法、修正变形梯度的方法和分叉理论等。冯虹等[3]指出利用第一和第二主应力来判别膜材的受力状态, 并根据其受力状态来修正单元的本构关系。谭峰[4]采用主应力-主应变判别准则和修正本构矩阵的方法对薄膜结构的褶皱进行了分析。本文采用适合薄膜结构褶皱分析的修正本构矩阵法解决物理非线性问题。在上述理论基础上, 通过典型马鞍形薄膜结构的案例分析, 总结薄膜结构振动特性以及其在不同平均风速和不同边界条件下的动力反应特性。

结构动力平衡方程为:

MU⋅⋅+CU˙+KU=P(t) (1)ΜU⋅⋅+CU˙+ΚU=Ρ(t)(1)

因式 (1) 对于任一时刻均成立, 所以结构在t+Δt时刻的平衡方程可表示为:

MU⋅⋅t+Δt+CUΜU⋅⋅t+Δt+CU⋅t+Δt+Kt+ΔtUt+Δt=Pt+Δt (2)⋅t+Δt+Κt+ΔtUt+Δt=Ρt+Δt(2)

设t时刻的位移Ut、速度U˙tU˙t和加速度U⋅⋅tU⋅⋅t已知, 据此求t+Δt时刻的位移Ut+Δt、速度U˙t+ΔtU˙t+Δt和加速度U⋅⋅t+ΔtU⋅⋅t+Δt。

由于薄膜结构在荷载作用下产生的位移较大, 会表现出几何非线性, 需利用迭代法求解每一时间步的动力平衡方程而得到Ut+Δt。常采用的迭代方法有全Newton-Raphson和修正Newton-Raphson迭代法等。全NR迭代法在一个荷载增量步内, 每次迭代求解前均更新刚度矩阵和非平衡载荷, 计算精度较高, 收敛性较好。但修正NR迭代法在一个荷载增量步内的各次迭代中均采用初始刚度矩阵, 仅更新非平衡载荷。比较两种方法, 本文采用全NR迭代法。

得到当前时间步的位移增量后, 需判别膜材是否发生了褶皱, 本文采用主应力-主应变判别准则[4]来判别膜单元的受力状态, 即:

(1) 当σ2>0时, 膜材处于纯拉状态;

(2) 当σ2≤0且e1>0时, 膜材处于单向褶皱状态;

(3) 当e1≤0时, 膜材处于双向褶皱状态。

其中, σ2和e1分别代表最小主应力和最大主应变。若出现褶皱的膜单元, 本文根据修正本构矩阵法[4]修正刚度矩阵后进行下一时间步的迭代分析。该方法首先确定主应力和主应力方向应变的关系, 然后利用单向褶皱时σ2=0和双向褶皱时σ1=0, σ2=0得到新的主应力和主应力方向应变的本构关系。具体方法见参考文献[4], 这里不再赘述。

振动形态分析是风荷载动力分析的基础, 因此下面给出一个马鞍形薄膜结构在简单动荷载作用下的基本形态, 总结薄膜结构振动形式的基本特点。

一个水平投影边长为6 m的方形平面马鞍形薄膜结构, 高点与低点的高差为2 m, 四边约束, 有限元模型如图1 (a) 所示。膜材采用PVC聚酯纤维膜, 密度为1.2 kg/m2, 经向弹性模量E1=900 MPa, 纬向弹性模量E2=600 MPa, 剪切模量G=20 MPa, 两泊松比分别为v12=0.4, v21=0.6, 两方向预应力σx=σy=2.0 kN/m。钢索采用5×13规格的钢索, 直径32 mm, 密度为每延米2 kg, 弹性模量E=1.9×105 MPa, 钢索预拉力为20 kN。图1 (b) 是该结构的找形结果图。

首先考察结构的中心点处有竖向初始速度时, 膜面的振动情况。假设结构中心点 (25点) 竖向初始速度0.5 m/s, 不考虑阻尼, 25点、17点和9点 (图1中用星号标出的点) 的竖向位移如图2所示。

由图2可见, 在初始条件的作用下, 25点从开始便向上产生位移, 0.015 s时该点位移达到最大值29 mm, 而此时17点、9点的位移分别为0.05 mm和0.005 mm, 即在25点产生最大位移的时刻, 17点、9点基本上还处于静止状态。在0.033 s时, 25点回到了平衡位置并开始向下运动时, 17点具有1.9 mm的向上位移, 并且继续向上移动, 9点位移仍然很小 (0.009 mm) 。当17点在0.056 s达到最大位移值7 mm时, 25点又从负的最大位移值回到了零位移的位置, 此时9点的位移只有0.24 mm。

图3给出了膜面在不同时刻的变形形状, 由这些图可以看出结构的运动是从结构中心点向周围逐渐传递的。在此后的过程中, 三点的峰值交错分布。以上结果表明, 由于膜材为柔性材料, 其平面外的抗剪能力很低, 垂直膜面作用力的传递速度较慢, 使整个膜面并非一起上下振动, 而呈波浪状, 结构的波动效应明显。

然后分析结构在冲击荷载作用下的响应, 假设结构的中心点 (25点) 处有100 N的突加集中力, 然后力持续作用。先进行静力分析, 后对结构在此冲击荷载作用下的动力响应进行分析, 并将二者的结果进行比较。

静力分析得到结构的最大竖向位移为25点, 最大值为21.8 mm, 而单元最大应力为2.05 kN/ m。

动力分析时阻尼为C=α1M+α2K, 其中阻尼系数分别为α1=0.1/ (1/f1+1/f2) , α2=0.1 (f1+f2) 。由薄膜结构自振特性分析程序[5]得到结构的第一频率f1=4.25, 第二频率f2=6.54, 从而可计算出α1=0.26, α2=0.009。取时间步长为0.01 s, 25点竖向位移时程曲线如图4所示, 其最终平衡位移为21.7 mm, 与静力计算的结果21.8 mm基本一致。25点在整个过程中位移最大值为27 mm, 位移放大系数为27/21.7=1.24。

图5为应力最大值随时间变化的曲线, 最终平衡时最大应力为2.06 kN/m, 与静力分析的结果2.05 kN/m基本一致。在整个过程中最大应力的最大值为2.09 kN/m, 应力放大系数为2.09/2.05=1.02。以上结果表明, 结构在冲击荷载作用下的平衡值与静力分析结果一致。结构的位移放大系数与应力放大系数不一致, 结构的非线性效应明显。

为了研究薄膜结构在风荷载作用下的动力响应特性, 讨论风速和边界条件对马鞍形薄膜结构动力响应特性的影响, 下面利用编制的薄膜结构动力分析程序对典型的马鞍形薄膜结构进行风荷载动力分析。

选取平面边长为10 m的马鞍形薄膜结构为计算模型, 膜材仍采用PVC聚酯纤维膜, 膜中索采用5×13规格的钢索, 结构最高点与最低点的高差是4 m。由于几何非线性的影响, 薄膜结构的刚度在变形过程中是变化的, 频率也随之变化[5], 因此, 首先需要确定初始形态。本文以形态分析得到的薄膜结构为初始状态进行风荷载非线性动力分析, 形态分析结果如图6所示。

风荷载采用下面表达式计算获得:

F=Cp12ρv20 (3)F=Cp12ρv02(3)

其中:F为作用于结构表面的风压力;Cp为风压系数, 1/2ρv2002为速压。本文采用Davenport谱生成了脉动风速时程 (图7) , 将脉动风速与平均风速相加得到风速时程曲线, 然后计算速压。对于风压系数, 本文利用FLUENT计算得到 (图8) , CFD模拟具体参数主要参考文献[6]。最后结合两者求得结构上各点的风荷载值。

在平均风速为5 m/s的风荷载作用下, 结构的动力响应如图9所示。此时, 结构中心点 (结点61) 的竖向位移随时间的变化曲线与风速的变化具有明显的一致性[图9 (a) , 图中虚线为静力分析结果], 结构明显地表现出受迫振动。动力分析最大位移值为66.4 mm, 是静力分析结果 (22.2 mm) 的3.0倍。由图9 (b) 可以看出, 结构中心位置单元111的第一主应力和第二主应力变化趋势完全相反。此时动力分析得到的第一主应力最大值 (2.174 kN/m) 只是静力分析结果 (2.098 kN/m) 的1.04倍。应力放大系数比位移放大系数小, 说明:变形较小时, 结构呈现出弱化性。

图9 (c) 是风速最大时结构的变形图, 结构的变形与风压系数的分布较为一致, 表现为B点附近的膜面变形较大。此时迎风面处应力分布比较均匀, 背风面处的应力分布较为复杂, 应力较大和应力较小的单元均分布在背风面处[图9 (d) ], 与风压系数的分布具有一致性。而且结构在风荷载作用下没有出现褶皱单元。

当平均风速为30 m/s时, 结构的响应如图10。由于出现了较多的双向褶皱单元, 位于结构中心处节点61的位移及单元111的应力随时间的变化复杂, 变化幅度较大, 出现较多的峰值[图10 (a) 和 (b) ]。动力分析得到的结构中心点最大竖向位移 (1.056 m) 是静力分析结果 (0.338 m) 的2.7倍。动力分析中单元111第一主应力最大值达到了82.781 kN/m, 是静力分析结果 (10.421 kN/m) 的7.9倍。

在47.8 s时结构的变形如图10 (c) 所示, 膜面上的竖向位移有正有负, 形状极为不规则, 此时结构中共有41个褶皱的单元[图10 (d) ], 其中1个为单向褶皱, 40个为双向褶皱。由图10 (e) 可以看出, 由于出现了较多的双向褶皱单元, 此时结构中的应力分布也极不均匀, 第一主应力值最小的只有-4.417 kN/m, 最大的达到了125.896 kN/m, 应力大小差别较大。并且结构中单元的第一主应力最大为125.896 kN/m, 已经超过了聚酯纤维膜材的C级抗拉强度, 说明此时膜材已经破坏。

由褶皱单元数量随时间变化的曲线[图10 (f) ]可以看出, 结构中的褶皱单元数量较多, 多数时间有30个褶皱单元, 最多时达到了110个, 超过了单元总数的一半。

由上面对柔性边界结构的分析可以得知, 当风荷载较小时结构位移和应力随时间的变化与风速的变化是一致的, 结构明显地表现为受迫振动, 结构的变形及应力分布与风压系数表现出一致性;当荷载较大时, 结构中出现较多的褶皱单元, 特别是较多的双向褶皱单元, 对结构的受力性能影响较大, 表现为单元的褶皱使结构弱化, 结构的位移和应力随时间变化比较复杂, 出现较多的峰值, 结构的变形不规则, 应力分布极不均匀等。

薄膜结构边界条件主要分为三种:刚性边界, 周围各边均固定在刚性结构上;柔性边界, 周围各边均为柔性索, 只有几个点与刚性结构相连;考虑支承柔性边界, 通过钢杆和地拉索将结构固定在基础上, 此时应将钢杆和地拉索组成的支承结构与薄膜结构一起分析。

对上述三种不同边界的薄膜结构进行风荷载非线性动力分析。柔性边界薄膜结构的计算模型见图6;图11 (a) 为刚性边界薄膜结构的计算模型, 它是以柔性边界结构的找形结果图为初始形状, 将其四边固定进行形态分析所得;图11 (b) 为考虑支承柔性边界薄膜结构的计算模型, 四个角点处各布置有一个竖直钢杆和两根地拉索, 结构最高点距地面10 m, 最低点距地面6 m。其中, 膜材采用PVC聚酯纤维膜, 膜中索采用5×13规格的钢索;八根地拉索规格为5×19, 直径35 mm, 预应力为51.47kN。两个高点处, 钢立杆的截面面积为6 566 mm2, 预压力为104.52 kN。两低点处, 钢立杆的截面面积为3 054 mm2, 预压力为73.79 kN。

分析结果表明:三种边界的薄膜结构受力性能基本相似。但结构在不同的边界条件下, 当风荷载较大时出现的褶皱单元数量有区别, 由于刚性边界结构的约束较多, 不易出现褶皱的单元, 在较大荷载作用下结构中褶皱的单元较少, 且多为单向褶皱[如图12 (a) ], 而两种柔性边界结构则较易出现褶皱的单元, 特别是较多的双向褶皱单元[如图12 (b) ]。

本文采用动力分析的Newmark法, 考虑几何非线性的全Newton-Raphson迭代法和考虑物理非线性的修正本构矩阵法, 实现了薄膜结构动力全过程分析。通过算例分析表明薄膜结构振动形态具有波动特性和非线性的特性。另外, 分析了不同风速和不同边界的情况下薄膜结构的非线性动力反应。结果表明:当风荷载较小, 不发生褶皱, 结构反应与荷载分布相关程度较高;当风荷载较大, 褶皱单元较多, 结构反应与荷载分布相关程度较低。而且不同边界条件的马鞍形薄膜结构出现的褶皱单元数量有区别。

通过本文的分析还发现, 当风荷载特别大时, 结构的变形会很大, 薄膜会产生褶皱, 工程师设计分析时应考虑褶皱影响。而且结构的风压系数也应该有较大的变化, 有必要重新计算变形后结构的风压系数, 即此时需要考虑薄膜结构与风荷载的耦合作用。要想解决这一跨专业的复杂问题, 则需要将上面的动力分析与流体力学相结合, 这是下一步薄膜结构分析需要解决的问题。