索膜结构是新型大跨空间结构形式之一, 具有质量轻、刚度较小等特点[1], 这些特点决定了其风敏感性.由此风荷载成为沿海地区这类结构的设计及防灾减灾分析的控制荷载[2,3].近年来, 索膜结构在国内外大型场馆建设中得到广泛百富策略白菜网, 人们对这类结构的风振响应及其抗风设计方法开展了一系列研究.

目前, 国内外对大跨索膜结构风振响应的研究方法主要有风洞实验、频域分析、时域分析和基于计算力学的气动弹性力学模型分析等[4,5,6,7].风洞实验成本较高并且受各种试验设备的限制, 在许多情况下无法实验模拟.频域法的随机振动理论仅仅适用于线性结构的随机振动.流固耦合的气动弹性力学模型又称数值风洞技术, 目前还不是十分完善[8,9].因此, 索膜结构抗风分析的较好方法是时域分析方法.

时域分析方法可以考虑频域分析中难以处理的非线性效应以及风与结构的相互作用, 得到结构风振响应的全过程, 这对结构的主被动设计是非常有意义的.本研究结合随机模拟时程法, 编写了索膜结构非线性风振响应分析程序, 计算了实际索膜结构的风振响应, 通过分析参数变化对结构风振响应的影响, 揭示了索膜结构风振响应的主要特征.

根据伯努利方程, 并且考虑风与结构的耦合作用, 结构上的风荷载P (t) 可用下式计算:

式中:ρa为空气质量密度;A为荷载作用面积;Cp为风压分布系数;v (t) 为风速;u?u? (t) 为结构振动速度.在实际计算中, 一般结构振动速度u?u? (t) 远小于风速v (t) , 因此可以忽略高阶小量12CpρaAu?2(t).12CpρaAu?2(t).

索膜离散结构节点处的动力平衡方程是一个复杂的非线性动力方程, 本研究采用Newmark逐步积分法和Newton-Raphson迭代理论, 推导出结构在t+Δt时刻增量形式的平衡方程:

式中:M为质量矩阵;D为阻尼矩阵;K为刚度矩阵;u为位移向量;R表示任一单元的节点力列阵;下标为t+Δt时表示t+Δt时刻对应的量.

假设方程 (2) 的解近似为ui−1t+Δtt+Δti-1.在ui−1t+Δtt+Δti-1点进行Taylor级数展开并忽略高次项, 则得到对应的近似线性公式:

Pt+Δt−Ri−1t+Δt−Mu¨i−1t+Δt−Dt+Δtu?i−1t+Δt−Ki−1t+Δtui−1t+Δt=Ρt+Δt-Rt+Δti-1-Μu¨t+Δti-1-Dt+Δtu?t+Δti-1-Κt+Δti-1ut+Δti-1=

Mu¨′t+Δt+Dt+Δtu?′t+Δt+Ki−1t+ΔtΜu¨′t+Δt+Dt+Δtu?′t+Δt+Κt+Δti-1ut+Δt-ui−1t+Δtt+Δti-1 (3)

式中:u¨u¨′t+Δt、u?u?′t+Δt分别为加速度和速度对位移的一阶导数;上标i-1表示迭代步数.

设t～t+Δt时间段内任一时刻的加速度u¨u¨介于u¨u¨t和u¨u¨t+Δt之间, 即

u¨=u¨t+δ(u¨t+Δt−u¨t), 0≤δ≤1 (4)u¨=u¨t+δ(u¨t+Δt-u¨t),0≤δ≤1(4)

而u¨u¨t+Δt在t时刻的一阶泰勒级数展开式为

u?t+Δt=u?t+u¨Δt (5)u?t+Δt=u?t+u¨Δt(5)

将式 (4) 代入式 (5) , 得到

u?t+Δt=u¨t+(1−δ)u¨tΔt+δu¨t+ΔtΔt (6)u?t+Δt=u¨t+(1-δ)u¨tΔt+δu¨t+ΔtΔt(6)

同理得到

ut+Δt=ut+u?tΔt+12u¨Δt2 (7)ut+Δt=ut+u?tΔt+12u¨Δt2(7)

选取不同的控制参数α, 由式 (4) 可得

u¨=u¨t+2α(u¨t+Δt−u¨t), 0≤α≤1 (8)u¨=u¨t+2α(u¨t+Δt-u¨t),0≤α≤1(8)

将式 (8) 代入式 (7) 得到

u¨t+Δt=a0(ut+Δt−ut)−a2u?t−a3u¨t (9)u¨t+Δt=a0(ut+Δt-ut)-a2u?t-a3u¨t(9)

式中:a0=1αΔt2a0=1αΔt2;a2=1αΔta2=1αΔt;a3=12α−1.a3=12α-1.

将式 (9) 代入式 (6) 可得

u?t+Δt=a1(ut+Δt−ut)−a4u?t−a5u¨t (10)u?t+Δt=a1(ut+Δt-ut)-a4u?t-a5u¨t(10)

式中:a1=δαΔta1=δαΔt;a4=δα−1a4=δα-1;a5=Δtδ2α−1δ2α-1.

将式 (9) 和式 (10) 对位移求导, 则有

{u¨′t+Δt=a0u?′t+Δt=a1 (11){u¨′t+Δt=a0u?′t+Δt=a1(11)

将式 (9) 至 (11) 代入式 (3) , 整理后得到

Kˉˉˉˉi−1t+ΔtΔui=Pˉˉˉi−1t+Δt (12)Κˉt+Δti-1Δui=Ρˉt+Δti-1(12)

式中:

Kˉˉˉˉi−1t+Δt=a0M+a2Dt+Δt+Ki−1t+Δt (13)Pˉˉˉi−1t+Δt=Pt+Δt−MΚˉt+Δti-1=a0Μ+a2Dt+Δt+Κt+Δti-1(13)Ρˉt+Δti-1=Ρt+Δt-Μ

a0ui−1t+Δt−ut)−a2u?t−ut+Δti-1-ut)-a2u?t-

a3u¨u¨t-Dt+Δta1ui−1t+Δtt+Δti-1-ut

Δui=ui−1t+Δtt+Δti-1-uit+Δtt+Δti (15)

当δ≥0.5, 且α≥0.25 (0.5+δ) 2时, Newmark算法无条件稳定.本研究取δ=0.5, α=0.25, 在每个时段内反复迭代计算, 直至收敛;从而可由前一时刻的位移、速度和加速度求出下一时刻的位移、速度和加速度.如此循环计算, 最终得到结构的风振响应全程 .

一般来说, 脉动风的卓越周期在1 min左右;高层建筑和高耸结构高阶振型的周期与此相差很远, 这些结构的风振响应是与一阶振型有关的.但是, 索膜结构的前几阶频率密集, 结构的风振响应不仅与其一阶振型有关, 还与其他低阶振型有关.因此, 本研究先进行该结构的自振特性计算, 再进行其风振响应分析.

本研究采用通用有限元程序对图1所示索膜结构进行自振特性计算.该索膜结构膜的厚度为1 mm, 泊松比μ=0.3, 质量密度ρ=1.25 kg/m2, 张拉刚度Et=250 kN/m, 预拉力S=30 N/cm, 脊索弹性模量E=2.1×105 MPa, 预拉力T1=120 kN, 长边索的预拉力T2=100 kN, 短边索的预拉力T3=30 kN.

通过计算分析, 得出结构前50阶自振频率和振型.膜的张拉刚度、膜预张力及结构跨度等结构参数改变对前4阶自振频率的影响见表1.

由表1可知, 索膜结构的自振频率较低, 且前几阶振型的频率分布密集;结构的自振频率随着其张拉刚度的增加略有增大, 随着膜预张力的增大而明显增大, 随着其跨度的增加而明显减小.因此, 索膜结构自振频率的主要影响因素是膜预张力和结构本身的跨度.

本研究采用自行编制的线性回归滤波器法 (AR法) 计算程序进行风速模拟. 模拟时水平风速谱选择Davenport提出的沿高度不变的水平风速谱:

Sv(f)=4Kvˉ210m2f(1+m2)4/3 (16)Sv(f)=4Κvˉ102m2f(1+m2)4/3(16)

竖向风速谱采用Lumley等[10]提出的经验公式:

Sw(n)=6fu2∗n(1+4f)2 (17)Sw(n)=6fu*2n(1+4f)2(17)

式中:vˉvˉ10为离地面10 m高处的平均风速;K为表面阻力系数; f为频率;m=1 200 fvˉ10m=1200fvˉ10;u*为摩擦速度 (剪切速度) ;n=f z/vˉ10,zn=fz/vˉ10,z为高度.

本研究计算的索膜结构的单元编号、节点编号及风压系数如图2～4所示, 风压系数是根据鞍形及抛物面索网结构选取的, 平均风速取15 m/s.结构在x和y两个垂直方向上的跨度分别为Lx=45 m, Ly=15 m;矢高为7.5 m.计算中积分时间步长取0.01 s, 结构响应总计算时间为120 s.

本研究仅列出迎风面和背风面部分节点及单元的计算结果.图5中给出索膜结构迎风面节点44和背风面节点47的位移响应时程, 从中可以看出膜面节点的振动均为受迫振动.迎风面风对结构的作用为上吸力, 风压系数较大, 节点44的位移向上, 振幅较大, 最大振幅达到0.5 m;而节点47处于背风面, 受到的风压力较小, 位移是向下的, 振幅也较小.

图6中给出部分膜单元的张力响应时程.由此可看出, 图2中的膜单元10位于迎风面支座附近, 风吸力很大, 膜单元的张力变化明显, 膜张力的峰值达到初始张力 (30 N/cm) 的2.5倍左右, 因此支座处属于该索膜结构的危险部位.当风荷载较大时, 索膜结构可能从这里最先破坏.背风面膜单元21的张力变化较小, 基本上在初始张力附近波动.

图7是部分索单元拉力的响应时程.从图7 (a) 可知, 长边索单元154的拉力变化较明显, 峰值可达到预拉力 (100 kN) 的2.5倍左右.由图7 (b) 看出, 脊索单元168的拉力变化较平缓, 但其拉力值均达到预拉力 (120 kN) 的2倍左右.因此, 索单元的拉力值较大, 在设计时应采用较大的安全系数, 避免钢索发生屈服.短边索拉力的变化规律类似于长边索拉力, 这里不再赘述.

通过计算分析, 发现索膜结构的风振响应与诸多因素有关.本研究考虑风向 (纵向风, 横向风及竖向风) 、风速、膜预张力等参数变化对结构风振响应的影响.

在实际情况中, 风可以从四面八方吹向结构.本研究第2节考虑的是纵向风, 以下分析横向风下的结构响应, 并与纵向风作用时的响应进行对比.选取横向风的风压系数如图8所示.

图9中给出不同水平风下结构风振响应的计算结果, 从位移响应图可以看出:结构在纵向风作用下, 位移均值在左半部分迎风面较大, 在右半部分背风面较小;在横向风作用下, 位移均值也是在迎风面较大, 在背风面较小.从膜张力响应图可以看出, 在纵向风作用下, 结构左半部分受吸力, 右半部分受压力, 结构的位移较大, 膜单元张力均值较大;在横向风作用下, 膜单元张力均值较小.位于支座附近的多数膜单元张力较大, 节点位移也较大, 而张力均值低于10N/cm的少数膜单元的节点振动方向与膜片初始曲率的方向相反, 出现了褶皱, 但不影响结构的整体性能.

从图9 (c) 中可以看出, 风向对本身刚度较大的索的拉力影响不大:在横风向作用下, 迎风面的长边索及脊索单元的拉力均值略有增大, 其他索单元的拉力均值基本不变化.

现有数据[2]表明, 实际风向与地面存在±10°的夹角, 竖向风对结构响应有一定的影响.笔者自编基于AR法的风速模拟程序, 模拟计算了结构的竖向风速时程, 进而计算了竖向风作用下结构的风振响应, 部分计算结果如图10所示.

从图10中可以看出, 结构在竖向风作用下的位移响应均值明显小于水平风下的响应, 膜单元张力响应均值相差较远, 不在一个数量级上.但是, 竖向风对结构风振响应仍有一定的影响:在竖向风作用下位移及膜张力响应均值为水平风作用下的13%左右;对长边索和脊索的拉力均值有较明显的影响.

风速对结构的响应有影响.本研究选取10、15、20、25 m/s 4种风速进行计算, 计算结果如图11所示.

由图11可以看出, 风速即风压对结构响应的影响较大.随着风速的增大, 结构的位移、膜单元张力呈非线性增长;索拉力均值大多呈非线性增长, 索单元168的索拉力均值呈非线性减小.多数情况下, 迎风面的节点及单元响应增长较快, 背风面的增长则相对较慢.

本研究选取20、25、30、35 N/cm 4种膜预张力参数, 讨论膜预张力变化对结构响应的影响.计算结果如图12所示.

由图12可以看出, 随着膜预张力的增大, 结构整体刚度变大, 结构的位移均值大多呈下降趋势, 而节点33的呈上升趋势, 迎风面节点的位移下降趋势较明显, 而背风面的下降不明显;此外, 索拉力均值变化不大.

本研究在索膜结构的自振特性的有限元分析基础上, 用自编程序对实际索膜结构进行了风振响应分析, 通过参数分析, 得出以下结论:

(1) 索膜结构的自振频率较低, 分布密集, 前几阶振型频率非常接近;随着索膜结构张拉刚度的增加, 自振频率略有增大;随着膜预张力的增大, 自振频率明显增大;随着结构跨度的增加, 结构自振频率明显减小.

(2) 风速对结构风振响应的影响较大.风速增大, 结构的位移、膜单元张力均值呈非线性增长, 索拉力均值大多呈非线性增长.多数情况下, 迎风面的响应增长较快, 背风面的增长较慢.

(3) 水平风对结构的位移影响很大, 在竖向风作用下位移及膜张力响应均值为水平风作用下的13%左右;对长边索和脊索的拉力均值有较明显的影响.

(4) 膜预张力与结构的位移均值呈相反关系, 且对迎风面节点的影响比背风面大, 但对索拉力均值的影响不大.

目前缺乏实际索膜结构的风洞实验数据, 脉动风与结构响应的耦合作用机理还需深入的研究.