张力膜结构是一种大跨柔性结构, 风荷载是其主要控制荷载。在风荷载作用下它会产生较大的变形, 而结构的变形又改变了作用于其上的风场, 从而形成了结构场与风场之间的耦合。计算表明[1], 若不考虑这种耦合效应, 分析结果会产生较大的误差。

由于涉及两种不同物理场的计算, 故而在计算中都是尽量将其中一种场加以简化。如流体工程师会将结构场简化为刚体或具有弹性支座的刚体, 反之结构工程师则趋向于将流体场简化为静荷载或用随机振动理论 (谐波合成法) 、统计理论 (AR线性自回归法) 等形成作用于结构上各点在时间及空间上相关的风荷载时程。对于与风场的耦合, 虽然可以用荷载项、附加质量、气动阻尼、气承刚度等因素加以考虑, 但结果并不是很好, 并且难以形成适用于各种结构形式的计算方法。而风洞试验方法由于成本高昂无法进行大规模研究, 并且风洞实验仅局限于某些特定的结构。随着计算机硬件及数值计算方法的飞速发展, 使得在计算机上综合运用计算流体力学 (CFD) 和计算结构力学 (CSD) , 即数值风洞方法来模拟结构及其周围绕流的运动过程成为一种可以选择的有效方法。

结构场的绕流属于大雷诺数、低速流动, 所以必须对湍流进行恰当的模拟。目前对湍流的模拟模型有大涡模拟 (LES) 、各种k-ε模型及k-ω模型。LES模拟的效果较好, 但其对网格划分的要求高, 计算量很大, 在进行三维分析时目前计算机的容量难以达到这个要求。标准k-ε模型对具有驻点及分离点的流动的模拟效果不好, 故在本文的计算中采用RNG k-ε湍流模型 (以下简称RNG模型) 。该模型对扩散方程的修正项改进了湍流动能的扩散, 所以对具有驻点及分离点流动的模拟结果较好, 比较适合于结构场的绕流问题。此外, RNG模型是两方程模型, 计算量与标准k-ε模型相当, 可以实现在普通微型机中模拟三维结构场与流场间的耦合问题。文献[2]对悬索桥桥面二维截面耦合风场的计算表明, 在各种湍流模型中, RNG模型及LES模型的计算结果与实验值最为接近, 而前者的计算量远小于后者。

对于不可压缩粘性流体, RNG模型控制方程如下:

连续性方程:

Δ⋅u=0 (1)Δ⋅u=0(1)

纳维-斯托克斯方程:

ρu?+ρu⋅Δu=−Δp+2(μ+ρCμk2ε)Δ⋅s (2)ρu?+ρu⋅Δu=-Δp+2(μ+ρCμk2ε)Δ⋅s(2)

湍流动能方程:

ρk?+ρu⋅Δk=(μ+ρCμk2σkε)(Δ⋅Δ)k+ρG−ρε (3)ρk?+ρu⋅Δk=(μ+ρCμk2σkε)(Δ⋅Δ)k+ρG-ρε(3)

消散率方程:

ρε?+ρu⋅Δε=(μ+ρCμk2σεε)(Δ⋅Δ)ε+C1εkρG−C2ρε2k (4)ρε?+ρu⋅Δε=(μ+ρCμk2σεε)(Δ⋅Δ)ε+C1εkρG-C2ρε2k(4)

式中:u, p, ρ, μ分别为平均速度, 平均静压力, 流体密度及动力粘性系数。

应变率s及动能积ρG分别定义为:

s=12(Δu+(Δu)T)ρG=2ρCμk2εs∶Δus=12(Δu+(Δu)Τ)ρG=2ρCμk2εs∶Δu

在RNG模型中, 方程 (4) 中的系数C1定义为:

C1=1.42−η(1−η/η0)/(1+βη3)C1=1.42-η(1-η/η0)/(1+βη3)

其中η=Sk/ε‚S=2s∶s−−−−−√η=Sk/ε‚S=2s∶s;Cμ、C2、σk、σε、β、η0为常数, 分别采用经验取值:Cμ=0.085, C2=1.68, σk=0.72, σε=0.72, β=0.012, η0=4.38[2]。

采用Taylor-Galerkin方法对控制方程中的非稳定项展开, 即可得到显式的有限元列式[2]。这种离散格式隐含了流线迎风的耗散作用, 具有较高的精度和稳定性, 可以有效地用于大雷诺数流动的求解。方程的求解采用逐次代换法, 这种方法在雷诺数较大的情况下求解纳维-斯托克斯方程时, 较牛顿法具有更大的收敛半径。

结构场控制方程:

Mu¨(t)+Cu?(t)+K(u(t))u(t)=F(t) (5)Μu¨(t)+Cu?(t)+Κ(u(t))u(t)=F(t)(5)

式中:M为质量矩阵, C为阻尼矩阵, 采用瑞利阻尼, 阻尼比取为0.05, K (u (t) ) 为刚度矩阵, F (t) 为作用于结构上的风荷载。

张力膜结构具有强几何非线性, 对结构场控制方程的求解采用Newmark积分法及全牛顿叠代法。

位移协调条件:

uF=uS (6)uF=uS(6)

平衡条件:

TF=TS (7)ΤF=ΤS(7)

式中, uF、uS分别为流体场及结构场的边界位移;TF、TS分别为流体场及结构场的边界力。

对两种场的求解采用序列交错的方法:

(1) 对流体场及结构场分别建模;

(2) 对流体场施加初始边界条件, 进行一轮稳态叠代求解后得初始流场分布。对结构场施加预张力;

(3) 从第 (2) 步的终点时刻开始, 根据结构场的位移相容条件先对流场施加边界条件, 叠代后求得流场速度及压力分布;

(4) 将流场的压力施加于结构表面, 叠代求解结构场;

(5) 对 (3) 、 (4) 步进行时间循环。

为了较好地进行瞬态求解, 时间步长必须足够小, 对于膜这类柔性结构, 可以取结构最低模态周期的1/20～1/30。

对一双坡张力膜屋面气弹模型进行了风洞试验, 模型几何缩尺比为1/50, 尺寸如图1a所示。屋面膜结构采用绸布作为模型材料, 并用超薄型粘性薄膜贴于表面以防止漏风。材料密度为ρ=643.8 kg/m3, 泊松比为ν=0.4, 弹性模量为E=2.5×107 N/m2, 施加预张力后用502胶水固定于四周, 薄膜在两个方向的张力均为N=8 N/m, 模型四周封闭, 采用刚性壁模拟, 第1阶频率为12.25 Hz。场地模拟B类地貌大气边界层, 测点位置如图1b所示。试验中共采用了0°～90°间的六种风向角, 分为四周封闭、四周敞开和突然开孔等三种情形, 但本算例仅对0°风向角的四周封闭模型作了计算分析。试验风速有5、5.5、6、7、8 m/s等五种, 信号采样频率为1 000 Hz, 采样持续时间90 s。

图2所示为0°风向角时五种风速下各测点的加速度均方根曲线。

按照相似准则, 推导出气弹模型所对应原型的张拉刚度为Et=2.5×107 N/m, 密度ρ=643.8 kg/m3, 泊松比为ν=0.4, 薄膜张力N=20 kN/m。模型与原型的风速比及时间缩尺比为1/50−−√1/50, 加速度相似系数为1。在原型周围取150 m×250 m×50 m空气柱进行计算, 入口风速输入由线性自回归过滤器随机生成的风速时程[4,5], 风速谱采用Kaimal谱以考虑风谱沿高度方向的变化, 侧壁、底面及膜面采用无滑移条件, 顶面为对称条件 (垂直壁面风速为0) , 出口断面采用零压力条件, 入口处湍流强度取10%。

图3及图4显示了原型风速56.35 m/s时测点2及测点5的加速度时程和各自的功率谱密度曲线。

为与不考虑流固耦合时的计算结果比较, 本算例还作了如下计算:将屋面膜结构模拟为刚性壁, 仅进行流体场的计算, 在流场入口处输出以上计算中生成的风速时程并纪录屋面各结点处三个方向的风速时程, 利用伯努利方程将其转换为结点上的压力, 再对结构场进行时程分析。图5分别为迎风面测点2和背风面测点5的计算结果及与试验值的比较。

从图5可以看出, 除个别值外, 计算结果与试验值相比吻合较好, 随风速增加计算误差有加大的趋势, 且迎风面的计算结果好于背风面, 这可能是由于RNG模型对于钝体绕流的处理还不够完善所致。图5的比较还显示, 考虑与不考虑流体—固体耦合的计算结果是不一样的, 特别是在高风速区, 这种差别则更加明显。

(1) 从分析结果可以看出, RNG模型及流体场与结构场的序列交错求解方法物理意义明确, 可以较好地模拟张力膜结构的耦合风振问题, 这给大跨度张力膜结构的耦合风振响应计算提供了一种非试验性的方法, 与直接进行风洞试验相比, 虽然计算精度上不如后者, 但它成本低, 耗时少, 对于普通膜结构的风振分析已能达到足够的精度;

(2) 对于张力膜这类柔性结构, 不考虑流体—固体耦合的分析结果误差较大, 但计算时间少, 可用于普通结构的初步设计;

(3) 限于目前的硬件条件, 直接数值风洞模拟方法非常耗费机时, 如以上简单的数值算例, 其每一种工况的30 s时程计算时间约为10 h, 可见对于大型复杂结构计算时间则更长, 所以该方法不便于进行诸如参数分析等大规模计算, 从这个角度来看, 对简化流体—固体耦合力学计算模型的研究仍是相当重要和必不可少的。