膜结构重量轻, 刚度小, 在风荷载作用下会产生较大位移、速度及加速度, 并有可能显著改变结构周围空气的绕流情况。因此, 分析其风致动力效应时, 应当考虑流固耦合的作用, 这也一直是众多学者努力的方向。

目前, 研究膜结构与风环境流固耦合作用FSI (Fluid Structure Interaction) 的方法可分为三类。武岳等[1]采用频域方法, 生成符合给定功率谱的随机风速时程, 根据膜结构体形系数将其转化为外力施加于膜结构, 并在时域内求解结构动力方程, 求解过程中通过风与结构的相对速度不断修正风荷载, 以此来考虑膜结构与风环境的“耦合”作用。该方法思路清晰, 理论成熟, 考虑了结构运动所引起的“卸载”和“超载”作用, 即引入了气动阻尼的作用, 但尚未能考虑由于结构运动引起的风场变化, 因此是一种半耦合模型。杨庆山等[2]提出了一种简化的气弹模型, 在结构运动方程中, 通过引入耦合参数 (附加质量、气动阻尼、气承刚度等) 修正结构运动方程, 来考虑膜结构与风环境的耦合作用, 耦合参数一般由经验解析公式或风洞实验确定。该方法全面地考虑了结构与风的耦合作用, 但由于膜结构体型复杂, 形态各异, 很多膜结构的气动耦合参数还不能通过理论计算得到, 或者其精度尚待于实验验证。Gluck[3]等基于CFD (Computational Fluid Dynamics) 和CSD (Computational Structure Dynamics) , 建立一种数值风洞模型, 通过中间数据交换平台传递CFD与CSD耦合数据, 间接地实现了膜结构与风环境的数值耦合作用。该方法可全面考虑结构与风场的耦合作用, 且不受模型复杂性约束, 但受CFD中湍流模拟等数值手段的影响, 其计算精度尚需试验鉴定。随着CFD理论及计算机硬件的发展, 该方法无疑具有良好的研究前景, 成为一种前沿的研究领域。

本文采用数值方法, 对典型马鞍面封闭式张拉薄膜结构的风致耦合作用进行了数值模拟, 得到了其响应时程, 对其风致响应特性作了分析, 并讨论耦合作用对其振动特性的影响。

一般而言, 薄膜结构风致耦合作用的数值模拟方法分为直接耦合和迭代耦合。直接耦合 (Direct Interaction) 也称为强耦合, 是指在每一时间步内对流体控制方程和结构控制方程联立求解, 所有的流体变量和结构变量, 通过求解单一的耦合方程同时得到更新。该方法概念清晰, 计算精度较高, 但直接耦合所导出的方程组非常庞大, 计算迭代需要巨大的数据存储和计算耗时, 且需要修正常规的流体和结构求解器, 编程负担重。直接耦合限于小规模耦合问题求解, 且在构造控制方程过程中往往对问题进行简化, 因此其百富策略白菜网有一定的局限性。迭代耦合 (Iterative Interaction) 又称为分区耦合 (Partitioned Interaction) , 依次求解流体和结构两个物理场, 其中一个物理场的结果作为外荷载或边界条件施加于另一个物理场。该方法准确性较难控制, 实现过程相对复杂, 但弱耦合分区求解占用内存较少, 且便于利用现有两个领域的研究成果, 可以根据不同领域的求解问题采用不同的求解方法, 适用范围较广, 因此, 在流固耦合数值模拟中为众多研究人员所采用[3,4,5], 本文采用迭代耦合方法作为数值研究的手段。

一个迭代耦合下 FSI 求解系统由三个主要部分组成:流体求解系统, 结构求解系统及一个用于以上两个求解系统的信息交换平台。采用迭代算法直至整个耦合方程的达到收敛解为止。

在流固耦合作用问题中, 由于耦合界面是可变形的, 一般用任意拉格朗日-欧拉ALE (Arbitrary Lagrangian-Eulerian) 坐标系描述流体流动, 流体流动的未知量既包括通常的流体变量 (压力、速度) , 也包括结构运动引起的流体网格运动速度。不可压缩流体ALE坐标下的流体控制方程为

连续方程:

vi, i=0 (1)

动量方程:

ρ(∂vi∂t+(vj−wj)vi,j)=τij,j+fBiρ(∂vi∂t+(vj-wj)vi,j)=τij,j+fiB (2)

本构方程:

τij=-p δij+2μ eij (3)

式中 ρ为流体密度, vi和wi分别为xi方向流体流动速度和网格运动速度, τij为应力张量, fBiiB为体力, p为压力, δij为Kronecker符号, μ为广义流体粘度, eij=(vi,j+vj,i)/2eij=(vi,j+vj,i)/2为速度应变张量。

利用上述流体控制方程求解流体流动的算法称为直接模拟DNS (Direct Numerical Simulation) , 其求解类似于近地层风场的高雷诺数湍流流动需要十分强大的计算资源, 目前阶段较难实行。因此有必要引入相关的湍流模型以节约计算资源。除DNS以外, 现行的湍流数值模拟方法主要有两种:基于时间平均的雷诺平均模型RANS (Reynolds Averaged NS) 和基于空间平均的大涡模型LES (Large Eddy Simulation) 。LES对大尺度漩涡进行直接求解, 通过亚格子尺度模型SGS (Sub-Grid Scale) 考虑小尺度涡对大涡的影响。大涡模拟以求解瞬时流场为出发点, 可体现压力场和速度场的脉动状况, 具有较高的精确性。因此得到更多学者的研究与百富策略白菜网, 正在不断的走向成熟[6]。LES的主要研究内容之一是亚格子模型的改进, 目前百富策略白菜网较为广泛的亚格子模型是基于涡粘假设的Smagorinsky亚格子模型以及动态亚格子模型 (Dynamic Subgrid Model) [7]。薄膜结构体型多样且复杂, 其周围流场时变性质显著, 相应的结构风致响应随机特性也十分强烈。鉴于上述时变性与随机性, 本文采用大涡模型来求解瞬态湍流流动。

引入Smagorinsky-lily LES亚格子模型:

μ=μ0+μt=μ0+2√ρ Λ2sDμ=μ0+μt=μ0+2ρΛs2D (4)

式中 μ0为流体动力粘度, μt为湍动粘度, Λs=min(κ d,kDV1/3)Λs=min(κd,kDV1/3)为亚格子混合长度, κ为Von Karman常数0.41, d为距壁面的距离, kD为无量纲模型常数, V为网格体积, D=eijeij−−−−√D=eijeij为有效变形张量。

薄膜结构在风荷载作用下的动力响应带有强烈的几何非线性特性, 其特点可归纳为小应变大变形。因此结构域求解考虑几何非线性的结构运动方程:

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

式中 M为质量矩阵, C为阻尼矩阵, K为刚度矩阵, u⋅⋅(t)‚u?(t)u⋅⋅(t)‚u?(t)和u(t)u(t)分别为结构运动的加速度、速度和位移, F(t)F(t)为流场对结构的作用向量。

迭代耦合算法中, 流体域与结构域之间数据的传递, 通过流固耦合边界平台实现。而流固耦合界面应满足的两个基本相容条件, 亦即迭代收敛判据为

1) 位移相容条件:dˉˉf=dˉˉsdˉf=dˉs

2) 应力平衡条件:n⋅τˉf=n⋅τˉsn⋅τˉf=n⋅τˉs

式中dˉˉfdˉf和dˉˉsdˉs分别为耦合界面处流体网格与结构网格的位移, τˉfτˉf和τˉsτˉs分别为耦合界面上流体与结构应力, n为外法线方向矢量。

此外, 如将膜结构表面作无滑移壁面处理, 还要求耦合界面上结构运动速度与流体运动速度相容, 即d−⋅s=vˉd-⋅s=vˉ。此式与位移相容条件统称为运动相容条件。

本文将一个时间层内的迭代求解称为内迭代, 一个内迭代循环分为三个阶段。 (1) 在某一时刻 (初始时刻或上一时间层) , 根据式 (1) ～式 (4) 求解第k次迭代下的流体域解向量Xkffk。 (2) 根据应力平衡条件为A n⋅τˉkf=n⋅τˉksAn⋅τˉfk=n⋅τˉsk, 将本次迭代解得的耦合界面处流体压力作为外荷载施加于结构, 由式 (5) 求解结构位移、速度等变量。 (3) 由位移相容性条件dˉˉkf=dˉˉksdˉfk=dˉsk更新耦合界面处流体节点位置, 并利用动网格求解技术更新其他流体节点坐标。内迭代收敛由位移相容条件以及应力相容条件同时控制:

rd=∥dˉˉks−dˉˉk−1s∥max{∥dˉˉks∥,ε0}≤εdrd=∥dˉsk-dˉsk-1∥max{∥dˉsk∥,ε0}≤εd (6)

rτ=∥τˉkf−τˉk−1f∥max{∥τˉks∥,ε0}≤ετrτ=∥τˉfk-τˉfk-1∥max{∥τˉsk∥,ε0}≤ετ (7)

式中 εd和ετ分别为位移和应力收敛容差, dˉˉksdˉsk和dˉˉk−1sdˉsk-1分别为本次与上次内迭代结构位移, τˉkfτˉfk和τˉk−1fτˉfk-1分别为本次与上次内迭代耦合面上的流体应力, ε0为预设常数 (10-8) 。

如本次 (第k次) 迭代同时满足式 (6) 和式 (7) , 则该时间层内迭代收敛, 进入下一时间层, 依次推进, 直至完成全部的薄膜结构风致耦合时程分析。

本文选取典型马鞍面封闭式张拉薄膜结构作为流固耦合分析模型, 如图1所示。模型平面投影为正四边形, 边长0.42 m, 四角点高低交错, 高点高度为0.2 m, 膜矢高f为0.075 m。膜材物理参数为:密度1000 kg/m3, 弹性模量5×105 N/m2, 厚度0.1 mm。一般而言, 薄膜结构只有在预应力作用下才具有初始刚度, 但为使薄膜与风的耦合作用最为明显, 本文选取施加预应力为零的工况。因膜厚度很小, 其受力可作为平面应力问题分析, 采用三维平面应力实体单元对结构进行有限元离散, 释放单元扭转自由度以模拟膜材“零”抗弯刚度的物理特性, 有限元模型如图2所示。

流体域方面, 采用有限元方法离散求解流固耦合微分方程组 (1) ～ (4) , 六面体结构网格总量约为18万, 单元插值函数为FCBI-C[8], 每个有限单元内设有27个积分点。采用原始变量SIMPLE算法求解压力速度耦合方程, 对流项采用二阶迎风格式离散以减小数值粘性的影响, 时间积分为二阶Runge-Kutta法, 时间步长0.01 s, 待计算累计一定时长, 初始条件的“瞬时效应”消散后开始统计, 统计计算时长为4.00 s。流体求解区域及细部网格分别如图3和图4所示。

采用均匀流入流面, 风速12 m/s, 风向如图所示。结构四周墙壁及流场底面 (地面) 为无滑移壁面, 流场域外侧壁及流场顶面采用滑移壁面条件, 薄膜屋顶处设置为流固耦合交界面, 用以流体域与结构域耦合迭代数据的传递。耦合分析以有限元商用软件ADINA为计算平台。

图5给出了数值与风洞实验中薄膜屋盖的平均风压系数, 可知数值解略小于风洞实验结果[9], 这可能由于风洞中带有一定的湍流强度, 与本文采用零湍流强度存在差异。总体上数值解能够较为准确地反映薄膜屋盖的风压分布特性。

薄膜结构在风荷载作用下竖向位移响应的平均值与均方根值等值线分别如图6和图7所示。

由平均位移等值线图6可知, 结构四周封闭, 建筑薄膜表面整体受负压风吸力作用, 结构整体上扬, 且由于四边固定, 膜面中部位移最大, 其值约为膜边长度的1/24。结构位移变化梯度比较均匀, 且沿两对角线连线的云图分布基本呈对称分布, 但沿两低点位移变化梯度较两低点梯度大。

需要指出的是, 由图5 (a) 可知屋盖上风荷载具有明显的非对称特性, 结构位移却沿对角线对称, 如图6所示。一方面, 受边界条件的影响, 结构周边刚度大而中点刚度最小, 较小的风压即可能引起中部较大的位移;另一方面, 膜面上沿AB运动的“行波”也增大了中部的位移 (如下文所述) 。

由位移均方根等值线图7可知, 高点对角线方向的位移均方根值明显高于低点对角线方向, 且III区附近位移均方根值明显高于其他部位, 这说明结构在该区域附近的振动幅度较大。分析可知, 均匀气流在结构迎风边处 (AB连线) 发生分离并产生分离漩涡, 这些分离涡极不稳定性, 随着时间的推进向下游脱落。高点A附近的气流碰撞与分离更为剧烈, 使结构I区的振动幅度较II区大;另一方面, 由于未对薄膜施加预应力, 向下游运动的漩涡卷带薄膜运动, 在膜面形成“行波”, 而结构在对角线方向的刚度更弱, 膜面“行波”总是沿对角线方向由I区向III区运动, 并在III区附近消散, “行波”的鞭梢效应使薄膜位移在III区附近被放大, 如图7所示。

为了研究薄膜结构上各部分的振动规律, 在膜面上选取平行于顺风向膜边的1～3点以及高点对角线方向1、4、5点研究膜面各部位的振动特性, 样本点位置如图1所示。

图8给出了1～3点脉动位移的幅值谱曲线。比较可知, 1、2两点的卓越振动频率基本相同, 约等于2 Hz, 但1点的振动幅值要明显高于2点, 正如前文有关图7的分析所述。3点的卓越振动频率较高, 约为8 Hz, 但其振动幅值较1、2点小, 且该点振动频率与1点的高阶振动频率相对应。

图9给出了高点对角线上三样本点脉动位移的幅值谱曲线。比较可知, 1、4、5三点的卓越振动频率依次增高, 5点的振动频率约为9 Hz, 4点次之, 约为6 Hz。5点振动幅值明显高于其他两点, 这正如对图7的分析, “行波”在III区附近的鞭梢效应显著增大了该区域的振动幅值与频率。

本文采用迭代耦合的方法对封闭式马鞍形薄膜结构的风致耦合效应进行了数值模拟研究, 数值结果表明, 结构的平均位移与结构振动形式、边界条件及风荷载分布有关, 而脉动位移与风场的涡形、演化规律以及结构形式的关系明显, 特别当薄膜上由于气流作用出现“行波”效用时, 空气与结构的相互作用亦愈加明显, 这也体现了数值耦合分析的优越性。在低预应力条件下, 薄膜上各部分的振动特性有明显区别, 气流分离、漩涡脱落强烈以及“行波”尾端鞭梢效应明显的区域, 脉动位移幅值及其振动频率较大。

影响薄膜结构风致耦合效应的因素还有:来流的平均风速、风向角、脉动风的特性、结构自身形状及特征尺寸和薄膜的预应力等。本文作为数值耦合方法的初探性研究, 尚未考虑过于复杂的影响因素, 但文中的研究证明了数值方法对于耦合作用效应研究的可行性, 并对薄膜的风致耦合效应做了一些有益的探讨, 为后续研究复杂流场及结构形式的影响做好铺垫。