膜结构是一种轻柔空间结构, 该类结构对地震具有良好的适应性, 而对风作用非常敏感, 风荷载是该类结构的控制荷载, 如图1所示。膜材在风荷载作用下会有较大的位移、速度和加速度反应, 使得结构发生剧烈的形状改变, 而此变化又会反过来影响风场的分布, 从而引起膜结构和风场之间的耦合作用[1]。因此要对膜结构的风荷载以及风振响应做出准确计算, 就需要考虑流固耦合的影响。对于膜结构风振中流固耦合问题的研究方法主要以风洞试验和数值模拟为主。

采用数值模拟方法研究膜结构风振中的流固耦合方法主要基于有限体积法、边界元法等。有限体积法需要对整个流场进行单元划分, 而边界元法仅仅在流场边界上离散计算区域。膜结构作为一种曲面结构, 其与流体的耦合恰恰发生在流场的边界, 边界解可以直接用于膜结构变形计算, 这恰好利用了边界元的优势。Sygulski[2]利用边界元法并考虑膜结构流固耦合效应对膜结构进行稳态分析, 不过距应用边界元法对膜结构流固耦合全过程进行数值模拟有很大距离。Xu等[3]将边界元法应用在膜结构与风耦合计算中, 在计算中假定流场为势流场, 为边界元法对膜结构的流固耦合研究打下了良好基础。Li等[4]采用弱耦合方法计算了马鞍面膜结构在势流场中的流固耦合, 由于流场为势流, 其计算结果与风洞试验结果有较大差异, 但在传统边界元方法中引入了非协调边界元技术, 可以很好地解决边界元数值求解过程中流体区域的尖角和尖边问题。

以上基于边界元法的工作限于势流理论的计算方法, 即假定空气的流动为非粘性。虽然空气的粘度较小, 但是如果完全忽略空气粘度的影响, 仍然会对计算结果产生影响, 甚至得到错误的结果。

在最近十余年涡方法已经被证明适用于不可压缩、粘性流的高精度计算[5,6,7];由于土木建筑中空气经常被假定为不可压缩、具有粘度的流体, 由此可知采用该方法计算膜结构流固耦合是可行的。同时该方法为无网格化, 可以避免对复杂的计算域进行网格划分。

基于以上原因, 本文将非协调边界元法引入涡方法即改进涡方法对膜结构的流固耦合进行计算, 通过与气弹风洞试验结果的比较分析, 验证了方法的正确性。

基于改进涡方法对膜结构的流固耦合进行计算时, 首先需要将整体计算域分解为两个部分:数字边界层区域, 即靠近壁面极薄的区域, 对于此区域的流场应用涡量-流函数为变量的Prandtl方程进行模拟;内部区域, 即整体计算域中除去数字边界层的区域, 控制方程为以涡量-流函数为变量的Navier-Stokes方程。

对三维Navier-Stokes方程进行变换, 使用涡量-流函数为变量代替初始变量得到三维不可压缩、粘性流动的涡量动力学方程为:

式中:x= (x, y, z) 为无量纲坐标向量, u (x, t) = (u, v, w) 为无量纲速度向量, t为无量纲时间, υ为雷诺数的倒数。

下面分别阐述内部计算域和数字边界层计算域。

对方程 (1) 可以采用粘性分解方法将每个时间步分解为两个部分:采用欧拉方程计算涡量场的变化, 再运用粘性扩散方程对涡量场的变化进行修正。

涡方法的本质是将整体连续的涡量场离散化为若干个涡元, 由此可以得到涡元的无网格化欧拉方程的Lagrangian形式解, 并采用时间积分方法可以得到欧拉方程的解。扩散方程的Green函数解同时为具有0均值和方差为2t/Re的三维Gaussian分布函数。因此, 每个时间步的粘性扩散过程可以表示为每个涡量单元的随机走步过程。

根据Helmholtz定律可以将流场中的速度矢量和速度梯度分解为由涡量场引起的部分与由势流场引起的部分, 由涡量场引起的速度可以由Biot-Savart定律求得。内部区域任意点的势流速度矢量可以采用边界元方法进行计算, 将边界离散为NE个单元则有:

利用法向速度边界条件, 可以通过边界元方程求解出所有边界点的速度势。根据边界元方程可知由边界点的速度势可以求解出内部点的速度及速度梯度。由于膜结构的边界为尖角, 所以在采用边界元方法求解势流时, 其尖角处的法向速度会产生突变。本文采用非协调单元法处理流场区域的尖角, 它将所有节点均取在光滑边界处, 从而避免了法向速度的不连续现象。

对于二维问题如图2所示, 假设初始划分单元时尖角处的单元节点为1, 单元12边界点的局部坐标分别为ξ1=-1, ξ2=1, 对应的差值函数为N1 (ξ) =0.5 (1-ξ) , N2 (ξ) =0.5 (1+ξ) , 由此得到该单元的法向速度为q=N1 (ξ) q1+N2 (ξ) q2, q1, q2分别为节点1、2的法向速度。而加入非协调边界元之后, 将单元边界点1'取在单元12内部距离, 距离尖角L/4处, 从而得到边界单元的新局部坐标为ξ1=-0.5, ξ2=1。利用更新坐标之后的差值函数的性质可以得到新的差值函数表达式为由此得到加入非协调边界元之后的单元法向速度q'=N'1 (ξ) q1+N'2 (ξ) q2。而三维单元同样可以由类似方法得到单元的插值函数。

经过以上方程的求解可以得到边界点的速度势, 由边界元方程的定义可知, 区域内点y的速度和速度梯度可以由边界点的速度势通过边界元方程得到:

式中:q, j=φ, jni, φ, j= (1/4πr3) (3r, ir, j+δij) , q, jl=φ, jlni, φ, jl= (3/4πr4) (-5r, ir, jr, l+δijr, l+δilr, j+δjlr, i) 。由式 (4) 与式 (5) 形式可知如采用传统边界元方法求解内部点的速度和速度梯度, 当内部点靠近边界时会产生强奇异积分与超强奇异积分。关于超强奇异积分国内外学者也做了相关研究[8,9,10], 本文采用文献[9-10]介绍的正则化算法对其奇异性进行消除, 其基本思想是将强奇异面积分和超强奇异面积分转化为沿单元围道上一系列线积分, 而高斯数值积分能够有效地计算这些线积分, 从而解析地计算出强奇异积分和超强奇异积分。

本节所有的变量都将采用局部坐标系表示。z向为平面的法向并且指向流体内部, z=0表示固壁表面, xy平面为固壁表面切平面, 整体流场及数字边界层的计算示意图如图3所示。

由于内部区域流场涡量的对流和扩散及应用法向非穿透条件之后, 会在壁面引起切向速度 (u+, v+) , 而静止壁面速度为 (u-, v-) = (0, 0) , 因此会在极小的薄层内产生速度阶跃, 即为涡量的产生。涡层强度为:

式中和它的离散形式即为面涡, u、v分别为x、y方向的速度。方程 (6) 定义了满足壁面非滑移边界条件的力学机理。同样可以被离散为一系列长方形的涡列, 边长分别为htxi和htyi, 在点 (xi, yi, 0) 处具有涡量强度值如果涡列强度大于涡列最大强度γmax (MSS) , 则将其分解为Nsi个涡列, 每个涡列具有的涡量为

一旦涡量在壁面产生, 它在数字边界层内的发展和变化可以由涡量-流函数为变量的Prandtl方程控制。同样对Prandtl方程进行时间积分, 可以得到涡列在数字边界层内的流动方程;而一旦有涡列进入固壁区域, 则舍去, 并在以固壁表面为对称轴的数字边界层内部位置产生一新的与原涡列性质相同的涡列。

当计算域的涡量场和速度场确定之后, 压力场就可以通过对Poisson方程和Neumann边界条件联合求解得到:

式中:B为Bernoulli函数应用该式可以得到Bernoulli函数的分布, 进而得到压力场分布。

由于改进涡方法计算得到的压力为无量纲化的, 为得到对应现实模型的压力值还需要将其进行变换, 方程如下:

式中:p*为无量纲压力, p为压力, ρ为流场密度, V为来流速度。根据式 (8) 即可以将无量纲压力转化为作用在现实膜结构上的风场压力。

膜结构流固耦合过程中由于流场即风场引起的膜面压力可以由第一节介绍的改进涡方法求解得到, 而膜结构的整体刚度矩阵可以由ansys对应的模拟模型提取。通过单元节点压力与由ansys提取出的整体刚度矩阵联合求解, 可以得到该时间步结束时的膜结构表面单元的位移。由膜结构表面的位移及壁面无穿透条件又可以通过非协调边界元方法求解得到下一个时间开始时计算域内的势流场, 而采用改进涡方法又可以计算得到该时间步膜结构表面压力, 由此迭代求解得到膜结构的流固耦合全过程。

采用改进涡方法模拟膜结构流固耦合的流程为:

(1) 建立膜结构的ansys模拟模型, 提取该模型的整体刚度矩阵;

(2) 对整体计算域流场进行网格划分;

(3) 运行改进涡方法程序, 计算得到该时间步结束时的膜结构表面压力;

(4) 将膜结构表面压力与ansys提取刚度矩阵联合求解, 得到单元节点位移;

(5) 更新膜结构表面的节点坐标, 将其带入步骤2;

(6) 重复计算2~5步, 直至得到膜结构流固耦合的全过程。

本文选用方形平屋盖膜结构验证本文方法的正确性。模型气弹风洞试验是在湖南大学风工程试验研究中心的HD-2大气边界层风洞的高速试验段完成的。该风洞气动轮廓全长53 m、宽18 m, 其中高速试验段的尺寸为宽3 m、高2.5 m、长17 m。采用激光位移计进行位移测量, 选用朗斯公司LC0408T压电加速度传感器进行加速度测量, 两者均采用东华DH-5920信号采集及数据处理系统进行测量记录和数据处理。

本次试验对象为平屋盖膜结构, 采用两种材性的模型形式, 即索网覆膜模型与直接张拉膜模型。其中索网覆膜模型结构刚度全部由索网提供, 覆面材料仅起到挡风作用。索网材料选用轻质软弹簧, 覆面材料选用轻质不透风薄膜;而直接张拉膜模型按照膜面初始预张力的缩尺比在模型膜片四边均匀施加预张力使模型最终成形, 膜材选用轻薄高弹薄膜 (以上材料的材性试验见3.3节) 。确定模型尺寸为正方形边长840mm, 屋檐高度200 mm。根据模型尺寸和风洞试验段横截面尺寸 (3 m×2.5 m) 计算可得本文计算风向角的堵塞度为2.24%, 满足风洞试验要求。两种模型底座的底部和四周均采用18 mm厚高强度粘合木板, 满足试验对模型围护结构刚度的要求。两种风洞试验模型及索网覆膜模型的节点构造如图4所示。

Cheng等[11]和Letchford等[12]分别通过一组半球壳体模型的风洞试验, 指出当半球壳的雷诺数在2.3×105-2×106之间时, 模型表面的风压分布不再随雷诺数变化而变化, 即模型进入雷诺数自模区, 对于进入自模区的相似准则可以不予模拟。本试验沿用这一结论, 确定三种风速工况为5 m/s、10 m/s、15 m/s, 对应风洞试验的雷诺数分别为2.8×105、5.6×105, 8.4×105、根据文献[11-12]可知本文试验在三种风速下均已进入雷诺数自模区, 故放弃对雷诺数的模拟。弗劳德数反映了重力场对风振的影响, 但是膜结构自重非常轻, 且具有一定的预张力, 弗劳德数影响不大, 仅当膜面处于松弛状态时应当考虑。因此本试验在相似比设计时, 尽量选取较大的预张力来削弱弗劳德数不相似造成的误差影响。

结构振动时会带动周围部分空气随之运动, 从而使膜结构的实际质量增加, 增加的部分即为结构的附加质量。不同于大多数土木工程结构, 由于膜结构的质量较轻, 附加质量与膜结构的质量可能在同一个数量级, 甚至更大。王磊[13]提出了可以计算二维薄膜附加质量的简化附加质量模型, 并通过试验对比验证了此简化模型的有效性。对各阶模态, 每一振型区域单位面积的薄膜附加质量ma为:

式中:α为附加质量系数, 圆形膜结构取0.65;ρ为空气密度;l为薄膜特征尺寸, 为该振型区域的内切圆直径。一般膜材密度为1 000 kg/m3左右, 厚度为1 mm, 则膜材质量为1 kg/m2。取空气密度ρ=1.29 kg/m3, 若原型结构尺寸l为30 m, 根据式 (9) 则第一阶模态的附加质量为0.65ρl=25 kg/m2, 附加质量则为膜材质量的25倍。显然对于轻柔的膜结构而言, 附加质量远超过薄膜自身质量, 故本试验对膜结构的质量相似稍作放松。

对于膜结构, 结构的刚度由两部分组成:一是由几何变形引起的几何刚度, 与膜面曲率和预张力有关;一是由材料应变引起的弹性刚度, 与材料的本构特性有关。通常情况下, 几何刚度要比弹性刚度大得多, 因此对模型几何刚度的模拟更为重要, 即本试验对于模型本身需要满足几何相似和初始预张力相似条件。

综上, 膜结构风洞试验的相似比选取如表1所示。

采用量纲分析法, 可知风压缩尺比λp为1/4、时间缩尺比λt为1/10, 进一步可得到可见本试验还满足斯托罗哈数St和欧拉数Eu相似。

根据3.2节分析可知, 膜材应该选择质量轻、同时其预张力范围应满足弹性段模型缩尺比要求的材料。索网覆膜模型中索网用直径3 mm的软弹簧代替, 覆膜选用厚度0.1 mm的TPU薄膜。直接张拉膜模型中薄膜选用厚度0.3 mm的乳胶膜。本文对弹簧和乳胶膜材两种材料各制作三组试件, 试件尺寸如图5所示。在东南大学材料实验室使用电子万能试验机CMT5105进行常温标准拉伸试验, 得到材料的力-位移曲线如图6所示。

由图6可见, 两种材料的三组试件力-位移曲线吻合, 说明所选材料性能稳定。根据表1可知索网覆膜模型对弹簧施加的预张力为20 N, 而由图6可见当拉力小于23 N时弹簧处于线弹性范围。对图6中三组试件计算求均值, 可得其刚度k=0.15 N/mm, 弹性模量为E=4.67×106Pa;同样根据表1可知直接张拉膜模型对膜材施加的单位预张力为100 N/m, 对应于图5中乳胶膜拉伸试件尺寸的拉力为7 N, 而由图6中乳胶膜的材性试验可知, 当拉力小于10 N时, 乳胶膜处于线弹性范围。试验测得乳胶膜的泊松比μ=0.38。

限制于实验条件, 对每个模型布置3个激光位移计, 测点布置及风向角如图7所示, 激光位移计布置在试验模型内部, 以避免对外部流场的影响, 如图8所示。风洞试验流场选用均匀流场, 试验选用三种风速工况:5 m/s、10 m/s、15 m/s。模拟均匀流场时, 试验段不放置任何湍流发生器, 风洞试验如图9所示。

根据量纲分析法可知:

式中:fm、fp分别为模型和原型的特征频率, tm、tp分别为模型和原型的时间, Dm、Dp分别为模型和原型的几何尺寸, Um、Up分别为风洞试验的风速和实际风速。本试验模型的几何缩尺比为1/20, 试验风速比为1/2, 实际风场的截止频率为1.5 Hz (根据相同模型在B类地貌下的Kaimal谱曲线的变化趋势 (图10) 可以看出, 此频率处的能量已经很小) , 代入上式可得fm=15 Hz。根据模型的前若干阶自振频率 (表2、表3) 可知, 模型振动频段集中在100 Hz以内。根据奈奎斯特定理:对于任意频带受限的信号, 用不低于两倍信号的最高频率对其采样, 不会丢失信息[14], 因此试验中的采样频率应不小于200 Hz。采样时长一般以保证对应于原型值不小于10 min为宜, 代入式 (11) 可得tm=60 s。故本试验信号采样频率定为500 Hz, 采样时长定为100 s。

根据表1风洞试验的三种风速工况, 确定三种工况下的特征风速分别为5 m/s、10 m/s、15 m/s, 入口处的无量纲风速恒定为1;特征长度根据膜结构边长选为840 mm;三种工况下计算模型的雷诺数根据特征风速、特征长度分别确定为2.8×105、5.6×105, 8.4×105。建立膜结构流固耦合的整体计算域 (图11 (a) 所示) , 整体计算域共划分738个矩形单元, 740个节点, 其中膜结构被均匀的划分为100个矩形单元 (图11 (b) 所示) ;应用非协调边界元之后, 节点增加到932个, 膜结构的尖角尖边处节点均取在了单元内部 (图11 (c) 所示) 。

由第2节介绍可知, 采用改进涡方法对膜结构进行流固耦合全过程模拟, 需要确定膜结构的整体刚度矩阵, 故本文在ansys中建立与两种风洞试验模型相对应的有限元模型。

Ansys中使用shell41单元模拟直接张拉膜模型, 并采用降温法对shell41单元施加预张力。由表1确定的膜结构单位长度预应力及张拉成形的乳胶膜材厚度0.18 mm可知, 风洞试验模型的膜材沿厚度方向的单位预应力为在ansys有限元模型中调整施加在shell41单元上的温度条件, 直至shell41单元中沿厚度方向的预应力同样达到0.556×106Pa。

由于改进涡方法计算中以整体膜结构作为计算对象, 而索网覆膜模型的刚度全部由弹簧提供, 膜面仅起挡风作用, 无法直接提取该类膜结构的整体刚度矩阵, 故需要建立与索网覆膜模型相对应的整体膜结构ansys有限元模型。本文根据自振频率近似相等原则建立对应模型:首先建立与风洞试验对应的ansys有限元模型, 轻质软弹簧和膜面分别采用link10单元和shell41单元模拟 (如图12 (a) 所示) , 由于膜面不受力, 所以只需要通过在link10单元中施加与表1中弹簧预张力相等的初始应变即可以确定ansys索网覆膜有限元模型, 进而进行模态分析。然后将索网覆膜计算模型转化为整体膜计算模型:采用shell41单元对整体膜面进行模拟 (如图12 (b) 所示) , 通过降温法对shell41单元施加预应力。调整预应力以调整模型的自振频率, 直至其达到与原索网覆膜模型相近。

试验中采用初位移法对风洞试验模型进行初始激励以测得模型各阶模态, 即对测点1 (图7) 施加向下初位移。Ansys有限元模型自振频率与试验模型自振频率如表2、表3所示。

由3.2节分析可知:对于轻柔的膜结构而言, 由空气引起的附加质量可能不容忽视。而本文在建立ansys有限元模型时并没有考虑空气附加质量的影响, 所以表2、表3中膜结构试验值普遍小于计算值。但是由于本文有限元模型的预张力已经与风洞试验模型的预张力一致, 所以ansys有限元模型中提取的刚度矩阵即为对应风洞试验模型的刚度矩阵。

采用改进涡方法计算时无量纲时间间隔Δt为0.05、最大涡列强度MSS为0.5、涡泡计算半径为0.1、数字边界层厚度b取槡总计算步均为500步。由于初始时刻时, 整体流场并没有充分发展, 所以计算结果取为100步至500步之间的位移均值与标准差。采用改进涡方法计算两种平屋盖膜结构模型得到的位移值与风洞试验位移值对比如表4、表5所示。

由表4、表5可知三种风速工况采用本文方法对膜结构流固耦合效应计算得到的位移均值结果大部分与风洞试验结果吻合, 基本在15%以内;而位移标准差的相对误差则较大, 基本在30%以内。原因如下: (1) 本文在计算膜结构位移时采用拟静力计算方法, 没有考虑动力学方程中的惯性力项与阻尼力项, 故一般而言反映平均风压效应的位移平均值误差较小, 而对于反映脉动风效应的位移标准差则误差增大。 (2) 改进涡方法中采用涡元的随机走步模拟涡量的扩算, 由于其随机走步的特性, 致使壁面的压力分布会产生随机涨落, 故无法精确模拟表征脉动风压特性的位移标准差。 (3) 均匀流场中, 来流未经阻碍而直接作用于膜结构表面, 湍流成分较少, 由此风吸力的平均效应显著大于脉动效应, 故位移标准差本身的基数较小。

因为采用了相同的计算参数, 两种模型的流固耦合全过程计算时间均约为25个小时。现行fluent软件无法单独对膜结构进行流固耦合计算, 需要与其他有限元软件联合求解膜结构的流固耦合过程。假定本文试验模型为刚性模型, 使用fluent软件进行其表面风压模拟时, 采用不同湍流模型对粘性流场进行模拟的计算时间不尽相同, 但达到稳定状态的时间一般最少在30小时以上。如果fluent软件再与其它有限元软件联合求解膜结构的流固耦合全过程, 则计算时间会更长, 甚至无法估计。由此可知本文将GMRES算法引入改进涡方法, 可以大幅度缩减计算时间。本文提出的算法不仅可以计算膜结构的流固耦合全过程, 而且需要的计算时间远小于fluent软件。

膜材在风荷载作用下会有较大的位移、速度和加速度反应, 使得结构发生形状改变, 而此变化又会反过来影响风场分布, 从而引起膜结构和风场之间的耦合作用。本文提出一种模拟膜结构流固耦合全过程的方法, 改进涡方法。该方法可以求解得到每个时间步的单元节点位移, 并由该时间步的单元节点位移求解下一时间步的流场分布, 由此模拟膜结构的流固耦合全过程。该方法具有以下特点:

(1) 在整体计算域中将Navier-Stokes控制方程转换为涡量-流函数方程, 流场中涡量的各种特性均可以由数学方法进行模拟。计算中将整体计算域分解为采用涡泡模拟的内部区域和用涡列模拟的数字边界层区域。计算过程中将非协调边界元方法引入涡方法整体计算过程中, 可以很好地解决计算域中的尖角尖边问题。在势流计算中采用正则化方法对近边界点的速度和速度梯度进行求解, 消除了计算过程中的奇异性。

(2) 通过改进涡方法在每个时间步求解得到的单元节点压力可以与由ansys提取出的刚度矩阵联合求解得到单元节点位移, 而通过改进涡方法又可以求解出下一个时间步由于节点位移而产生变化的整体计算域的流场分布, 以此达到模拟膜结构流固耦合全过程的目的。

(3) 将GMRES算法引入改进涡方法, 可以在保证计算精度的前提下, 大幅度缩减计算时间。

(4) 本文方法同风洞试验结果相比, 大部分位移平均值相对误差基本在15%以内;位移标准差相对误差基本在30%以内。