薄膜结构在风荷载的作用下通常会产生较大的变形和振动, 这种大幅的变形和振动反过来也会影响到其表面风压分布, 产生所谓的“流固耦合”效应。目前, 国内外在这一领域的研究均处于起步阶段。武岳和孙晓颖[1,2]建立了适用于索膜结构流固耦合风振分析的数值风洞方法。百富策略白菜网该方法, 实现了对二维悬臂板和单向柔性屋盖的流固耦合效应模拟。本文在文献[2,3]的理论框架下, 综合运用计算流体力学软件FLUENT6.0和自行编制的膜结构动力分析程序MDLFX, 在Compaq Visual Fortran 6.5环境下搭建了薄膜结构流固耦合的CFD数值模拟平台[4], 并就数值模拟中的关键技术问题进行了系统研究, 编制了相应的计算程序。基于该软件平台, 对单向柔性屋盖和鞍形膜结构屋盖的流固耦合进行了数值模拟, 验证了方法的有效性。

薄膜结构在风荷载作用下的耦合振动问题, 在理论上可归结为不可压缩粘性流体与几何非线性弹性体之间的非定常耦联振动问题。对这一问题的求解包括流体域、结构域和网格域三个计算模块: (1) 流体域主要是模拟近地面大气边界层风场, 属于钝体空气动力学范畴; (2) 结构域主要是模拟张拉膜结构的风致动力响应, 属于几何非线性弹性体的大位移、小应变受迫振动问题; (3) 网格域主要是以ALE技术为基础的动态网格计算问题以及流体-结构网格之间的数据传递问题。基于上述求解策略, 本文构造了薄膜结构流固耦合的数值模拟平台如图1所示。

近地面风可假设为低速、不可压缩及粘性的牛顿流体, 其基本控制方程为表达动量守恒的N-S方程和质量守恒的连续方程, 其一般形式可以写为

ddt∫Vddt∫VρΦ dV+∫S[ρ (u-ug) Φ-ΓΦgradΦ]·n dS

=∫VSΦdV (1)

对于连续方程, 取Φ=1;对于动量方程, 取Φ={u, v, w};方程中的扩散系数ΓΦ和源项SΦ需根据Φ来选取。

对于式 (1) 需要引入湍流模型使流场均值化, 对难以分辨的小尺度涡在均值化过程中加以模式化, 使其在湍流模型中体现出来。目前均值化方法有两种: (1) 系统平均, 对应湍流模型为雷诺应力方程模型 (RANS) ; (2) 空间滤波, 对应湍流模型为大涡模拟 (LES) 。

本文采用URANS模型来描述瞬态流场, 以实现非定常流固耦合的数值模拟。该想法是通过Fluent6.0提供的UDF编程实现的。

结构的运动方程为

[M]{x¨(t)}+[C]{x?(t)}+[K]{x(t)}={P(t)} (2)[Μ]{x¨(t)}+[C]{x?(t)}+[Κ]{x(t)}={Ρ(t)}(2)

式中 [M]为集中质量矩阵, [C]采用Rayleigh阻尼, [K]为结构的切线刚度矩阵, {P (t) }表示由流体域施加给结构的荷载。

膜结构在风振响应过程中的几何形态是不断变化的, 必须及时为CFD计算提供这一信息, 以保证CFD网格能随着结构的变形而适应至新的位置, 这就需要引入动网格技术。本文综合运用代数法和迭代法的优势, 编制了动网格变形程序, 有效地实现了薄膜结构流固耦合运算中动网格的更新。基本思想[7] 如图2所示, 在流体域网格划分时, 将计算域划分成多个子块, 首先利用弹簧法确定出各子块之间的角点坐标, 即假设各子块的角点之间通过弹簧相连, 当处于耦合边界上的角点移动时, 通过弹簧网络系统的平衡可以确定其他角点的新坐标;一旦所有角点的位置确定后, 就可将角点信息传递给各个子块;然后对于每个子块, 采用基于弧长的无限插值法TFI (Transfinite Interpolation Method) 生成和改变这一子块的面网格及内部的体网格。

TFI对于每个子块, 采用基于弧长的无限插值法进行求解, 可分为以下几个步骤:① 将计算区域各表面的网格点参数化;② 计算所有边和角点的变形;③ 执行一维, 二维和三维无限插值法求解变形;④ 将变形和原始网格累加起来得到最终网格。

无限插值法 (TFI) 能有效地求解小位移问题, 因此适用于解决单个子块中的网格运动问题。而各子块之间的角点运动属于大位移问题, 很难由代数方法决定, 所以引入弹簧的概念, 假设各分区角点之间通过弹簧相连。弹簧刚度反比于该边的边长, 表达如下:

ki +1/2, j, k=1[ (xi +1, j, k-xi, j, k) 2+ (yi +1, j, k-yi, j, k) 2+

(zi +1, j, k-zi, j, k) 2]1/2 (3)

由于CFD和CSD计算对所需网格密度要求不同, 导致耦合界面上的网格不相匹配。通常CFD要求的网格密度比CSD密得多, 由此产生耦合界面上两套非匹配网格之间的数据传递问题如图3所示, 这在数学上是双向插值问题。本文基于薄板样条法, 编制了插值计算程序, 有效地实现流固耦合交界面上CFD和CSD网格之间的数据传递。

薄板样条插值是一个多变量插值问题, 即从待配准的两幅图像中 (图像维数为d) , 标定n对对应标记点pi和qi, 在一个合适的Hilbert空间H上寻找连续变换g:Hd→Rd, 满足: (1) 最小化一个给定的泛函J (g) :H→R; (2) 实现如下插值:qi=g (pi) 。

通过薄板样条矩阵[g], 可将CSD网格上的变形向量{uS}转换成CFD网格上的变形向量{uF}, CFD求解器在新的状态下求解流场, 再将CFD网格上的流体荷载向量{FF}按式 (4) 转换成CSD网格上的荷载向量{FS}。

{uF}=[g]{uS}, {FS}=[g]T{FF} (4, 5)

单向柔性屋盖的计算模型如图4所示, 在屋盖上等距离布置9个测点, 取屋盖高度H=10 m, 跨度L=40 m, 入口风速V=20 m/s, 模拟B类地貌, 流动雷诺数Re=1.38×107, 无量纲时间步长Δt=0.005, 屋盖结构采用索单元模拟, 单位长度质量g=5 kg/m, 张拉刚度E t=5.5×106 N/m, 预张力T=20 kN/m, 根据计算, 结构前两阶自振频率分别为0.79 Hz和1.58 Hz。

表1给出了平屋盖刚性模型和弹性模型各测点风压系数的均值与方差, 可以看出:① 屋盖表面的风荷载以吸力为主, 靠近屋盖前缘的部位风吸力较大, 随着距离的加大, 风吸力逐渐减弱。② 同刚性模型相比, 弹性模型的平均风压在屋盖前缘明显降低, 后缘明显增大, 但脉动风压在整体上呈减小趋势。如图5所示, 本文的数值计算结果和文献[2,3]的数值模拟结果比较, 可以发现两者在变化趋势上吻合较好, 而在具体数值上还存在一定的差异。

图6给出了平屋盖采用本文数值模拟方法 (考虑流固耦合作用) 和随机振动时域分析方法 (不考虑流固耦合作用) 求得的各点最大位移;可以看出, 两者的屋盖最大位移差值约为15%左右, 而以往采用的随机振动时域分析方法计算结果偏大, 可能会造成设计偏于保守。

以菱形平面鞍形膜结构为例, 计算模型如图7所示, 四周封闭, 取跨度L=28 m, 矢跨比f/L=1/12, 预张力T=2.5 kN/m, 薄膜厚度为1 mm, 单位面积的质量g=1.25 kg/m2, 张拉刚度为E t=2550 N/cm, G t=800 N/cm, 泊松比为0.6。来流风向为45°, 入口风速V=20 m/s, 模拟B类地貌, 时间步长Δ t=0.005 s。根据计算, 结构前两阶自振频率分别为2.26 Hz和3.03 Hz。

根据鞍形屋盖结构布置形式及风压分布特点, 将鞍形屋盖分为角区 (Ac, Bc, Cc, Dc) 、边区 (Ae, Be, Ce, De) 和中区 (Am, Bm, Cm, Dm) 等共12个区域如图8所示。图9给出了45°风向下刚性模型和弹性模型各分区的风压系数 (括号里的表示脉动风压系数, 括号外的表示平均风压系数) 。可以看出, 与刚性模型相比, 弹性模型屋盖表面的平均风压整体上呈增大趋势, 约增大了7%左右, 而脉动风压变化幅度不大, 可见对于柔性结构, 仅通过刚性模型来确定结构的风压分布, 可能偏于不安全。将本文的计算结果与文献[3]中风洞测压试验结果 (图10) 比较, 可以发现二者的风压变化趋势吻合较好。文献[3]中气弹模型试验结果与刚性模型试验相比, 整个屋面的平均风吸力约增大10%左右, 而脉动风压的变化幅度很小。由于试验还存在一定的误差, 因此, 本文的数值模拟结果与文献[3]中试验结果在具体数值上还存在一定的差异。

图11给出了采用两种方法求得的屋盖各点位移最大值分布图, 可以看出以往采用的随机振动时域分析方法 (不考虑流固耦合作用) 计算结果偏小, 可能会造成设计偏于不安全。这主要是因为45°风向下来流垂直于屋盖边缘, 屋盖曲面形状对绕流的影响较大, 屋盖的双向曲率对气流的旋涡脱落的影响均很大。考虑流固耦合作用时, 屋盖在风吸力作用下逐渐抬起, 气流分离作用增强, 风吸力增大, 而脉动风压的变化很小, 因此45°风向下鞍形屋盖的流固耦合效应较为明显。

根据文献[2,3]和本文的研究可以发现, 流固耦合的作用使结构的振动特性和绕流特性都发生了明显变化:① 结构在平均风荷载作用下产生的静态变形, 反过来改变了平均风压在结构表面的分布, 导致刚性模型与弹性模型间存在差异;② 脉动风使结构在平衡位置附近发生瞬态振动, 风与结构之间的相互作用使结构振动形态很大程度上受到旋涡脱落的影响, 并不是按照其自振频率振动。因而可以考虑分别从平均响应和动力响应入手研究结构的流固耦合效应。

在文献[3]的理论框架下, 综合运用计算流体力学软件FLUENT6.0和自行编制的膜结构动力分析程序MDLFX, 在Compaq Visual Fortran6.5环境下搭建了薄膜结构流固耦合的数值模拟平台, 并基于该软件平台, 对单向柔性屋盖和鞍形膜结构屋盖进行了流固耦合风致振动的数值模拟。主要结论如下。

(1) 基于弱耦合分区求解策略, 建立了三维薄膜结构流固耦合的数值模拟平台。数值模拟平台主要包括两部分:几何建模和流固耦合数值模拟 (包含了CFD计算模块、CSD计算模块以及数据交换平台Interface三部分) 。

(2) 综合运用代数法和迭代法, 编制了动网格变形程序, 有效实现了薄膜结构流固耦合运算中动网格的更新。在流体域网格划分时, 将计算域划分成多个子块, 首先利用弹簧法确定出各子块之间的角点坐标, 然后对于每个子块, 采用基于弧长的无限插值法TFI生成和改变这一子块的面网格及内部的体网格。

(3) 采用薄板样条法, 编制了插值计算程序, 有效实现了流固耦合交界面上CFD和CSD网格之间的数据传递。

(4) 流固耦合作用使结构的振动特性和绕流特性都发生了明显变化:① 结构在平均风荷载作用下产生的静态变形, 反过来改变了平均风压在结构表面的分布, 导致刚性模型与弹性模型间存在差异, 因此对于薄膜结构, 仅通过刚性模型来确定结构风压是不可靠的;② 风与结构之间的相互耦合作用使得结构振动形态很大程度上受到旋涡脱落的影响, 并不是按照其自振频率振动。