当前, 建筑物的风振计算可分为两类:即建筑物的自身振动对脉动风力变化 (产生附加气动力) 有显著影响和没有影响。后者可以通过先求得静止建筑物的风力荷载, 再进行结构的动力响应分析的方法进行[1]。对于索膜结构等柔性结构而言, 由于结构在风力作用下的变形相对结构本身尺寸相当大时, 附加气动力的作用不可忽略, 因此必须采用考虑流固耦合的计算与分析。

流固耦合计算又分为两种:一种称之为紧耦合 (或强耦合) 方法, 是对流体与固体的联立方程式进行一起求解。例如, 流体与固体通过有限元方法进行离散化后再组合形成一个矩阵进行求解。这个方法的主要缺点是不能使用已有的得到广泛百富策略白菜网的流体和结构求解器, 而且关于积分的推导也是非常冗长乏味的。因此, 在实际百富策略白菜网中, 复杂的流场跟弹性变形较大的结构之间的相互作用的模拟就变得非常费时费力。另一种称之为松耦合 (或弱耦合) 方法, 即在每一个时间步内, 通过流体计算得到风压、风力并将其输入到结构动力计算中, 求得结构的位移、速度及加速度等动力响应。然后, 基于这些响应在流体计算中对建筑周围的计算网格进行移动从而改变建筑表面的边界条件, 重新计算结构表面的风荷载。松耦合方法最大的优势是可以利用现有的通用流体和结构软件, 并且可以分别对每一个软件单独地制定合适的求解方法, 而流体以及固体的求解器不需要改写, 可以通过联合现有的软件的方法从而得到耦合问题的一个精确的求解。此外, 这种方法还可以针对一个CFD程序选择多个CSD程序。因此可以根据工程上的不同需求, 进行流体和结构模型的相互交换, 使这种方法更具有一般性。但是, 松耦合方法也存在计算效率和收敛性等问题的限制[2]。近年来, 随着功能强大、具有二次开发功能的商用CFD (Computational Fluid Dynamics) 和FEM (Finite Element Mode) 软件以及功能更强大的网格剖分工具的出现, 基于松耦合方法的流固耦合数值模拟技术的实现逐步成为可能[3,4]。

本文根据三角形重心坐标系原理提出了基于松耦合计算的流固耦合网格协调插值新方法, 提高了流固耦合的计算效率;同时, 为解决由于柔性结构加速度过大而导致的求解不稳定问题, 通过在1个时间步上进行流体及结构间反复交替求解来得到稳定的时程解。最后, 通过与风洞试验的对比验证了本文方法的可靠性。

影响流固耦合计算精度的因素主要有以下几个方面:流体计算的精度、结构计算的精度和流固耦合界面的插值精度。而流固耦合界面插值精度又包括了节点力的插值和节点位移的插值两大问题。前者通常可以采用“邻近搜索”等方法实现, 而后者往往由于不同网格体系信息量的不对称导致插值计算的困难[5]。

目前常用的搜索算法主要有以下几种:邻近搜索, 该方法的困难在于耦合域可能分布于不同程序的好几个位置上, 并且不同的程序的耦合面之间还可能存在着初始距离。强力搜索, 即在所有的单元上执行一次循环, 显而易见, 这种方法的效率是比较低的。此外, 还有八叉树搜索、超前搜索、桶式搜索等。所有这些算法都基于不同的假设, 并且使用了不同的数据结构, 并且都具有相应的复杂性, 因此给流固耦合的工程百富策略白菜网带来的效率上的困难。

到目前为止, 用于流固耦合网格更新的方法很多, 如弹簧近似法[6,7,8]等。这些算法虽然可以保证网格单元的拓扑关系不发生改变, 但在处理大变形问题时容易发生单元交错, 从而导致运算失败。虽然网格实时重划[9]计算可以解决上述问题, 但该方法耗费机时, 且在网格重画过程中将破坏原有网格的拓扑结构, 从而带来插值误差。

在索膜结构流固耦合计算中通常使用三角形单元。虽然这些三角形都是空间的, 但三角形却是一个天生的2D物体, 原因是三角形的三个点在同一平面上。图1中三角形代表FEM模型中的一个结构单元, A、B、C分别为结构单元的节点, 点P代表落在三角形ABC内部的CFD网格的节点。假设以点A作为起点, 那么点B相当于在AB方向移动一段距离得到, 而点C相当于在AC方向移动一段距离得到。则对于平面内任意一点, 可由如下方程来表示:

P=A+u× (C-A) +v× (B-A) (1)

整理方程 (1) 得到

P-A=u× (C-A) +v× (B-A) (2)

令v0=C-A, v1=B-A, v2=P-A则,

v2=u×v0+v×v1 (3)

等式 (3) 两边分别点乘v0和v1得到两个等式:

(v2) ·v0= (u×v0+v×v1) ·v0 (4)

(v2) ·v1= (u×v0+v×v1) ·v1 (5)

求解上述方程 (4) 、 (5) 得到:

u=(v1⋅v1)(v2⋅v0)−(v1⋅v0)(v2⋅v1)(v0⋅v0)(v1⋅v1)−(v0⋅v1)(v1⋅v0)u=(v1⋅v1)(v2⋅v0)-(v1⋅v0)(v2⋅v1)(v0⋅v0)(v1⋅v1)-(v0⋅v1)(v1⋅v0) (6)

v=(v0⋅v0)(v2⋅v1)−(v0⋅v1)(v2⋅v0)(v0⋅v0)(v1⋅v1)−(v0⋅v1)(v1⋅v0)v=(v0⋅v0)(v2⋅v1)-(v0⋅v1)(v2⋅v0)(v0⋅v0)(v1⋅v1)-(v0⋅v1)(v1⋅v0) (7)

如果系数u或v为负值, 那么相当于朝相反方向移动, 即BA或CA。因此, 如果:u≥0, v≥0, 且u+v≤1则可以判断点P位于三角形ABC内部。上述公式 (6) 、 (7) 中仅包含少量的四则运算, 因此在计算效率非常高。

基于上述坐标系统, 原先位于一个结构三角形单元内部的流体网格会随着结构的变形依然贴附于原三角形内。并且随着三角形在平面内的变形, 三角形内部的点仅发生相对位置的变形, 而变形前后的比例关系保持不变, 如图2所示。

在进行流固耦合计算时, 流体网格一般总要比结构网格密, 即一个结构单元 (以三角形单元为例) 内总会包含多个流体网格节点。因此, 以一个三角形结构单元建立的三角形重心坐标系去定位其中所包含的各个流体网格节点, 既能保证各节点在变形前后的相对位置保持不变, 同时又能提高计算效率。

工程界已经可以通过使用相对成熟的有限元软件对结构在风荷载下的应力、变形等效应进行精确计算。同时, 也能够通过使用计算流体动力学软件对结构所受的风荷载进行准确估计。然而, 对于进行流固耦合计算的方法和工程百富策略白菜网却少之又少。传统的松耦合在一个时间步内仅需进行一个循环的流体计算及结构动力响应计算, 这种方法对于像空气与建筑物这样相对质量比大的情况有效。但是, 当由建筑物变形产生的附加空气质量不能忽略时, 以上方法会由于在结构上生成过大加速度从而导致求解不稳定。文献[2]指出可通过在1个时间步上进行流体及结构间反复交替、不断修正求解, 解决以往柔性结构流固耦合计算中的收敛性问题, 但未对具体方法进行说明。本文以此为思想, 提出了一种解决流固耦合计算收敛性的方法, 具体过程如下:

图3所示为传统松耦合算法。图中Fi为i时刻结构所受荷载, Fi+1为下一时刻i+1时刻需要加载的量。按照传统算法, 即使荷载步被细分, 如图中所示分为n个荷载步, 每次增加dF/n的量, 最终和一次加载的效果是一样的。图4为改进方法, 即荷载步被细分, 但Fi+1也随时进行更新, 从第二个荷载子步开始, 每次增量dF/n变成dF′/n、dF″/n、dF?/n…以此类推。因此, 最终的加载效果如图4中曲线所示。事实上, 这种算法更符合流固耦合的实际发生机理, 同时也使数值计算的运行性能得以改善, 使计算过程更加容易收敛。

Huang[10]提出了一种新的信号处理方法—经验模态分解方法 (Empirical Mode Decomposition, EMD) 。Huang[11]又将该方法进行了一些改进。该方法从本质上讲是对一个信号进行平稳化处理, 其结果是将信号中不同尺度的波动或趋势逐级分解开来, 产生一系列具有不同特征尺度的数据序列, 每一个序列称为一个本征模函数 (Intrinsic Mode Function, IMF) 分量。最低频的IMF分量通常情况下代表原始信号的趋势或均值。作为一种百富策略白菜网, EMD分解方法可以有效地提取一个数据序列的趋势或去掉该数据序列的均值。

EMD分解方法的基本思想是:加入一个原始数据序列X (t) 的极大值或极小值数目比上跨零点 (或下跨零点) 的数目多2个 (或2个以上) , 则该数据序列就需要进行平稳化处理。最后, 原始的数据序列即可由一系列IMF分量以及一个均值或趋势项表示:

X(t)=∑j=1nCj(t)+rn(t)X(t)=∑j=1nCj(t)+rn(t) (8)

由于每一个IMF分量是代表一组特征尺度的数据序列, 该过程实际上将原始数据序列分解为各种不同特征波动的叠加。

EMD分解的主要目的之一是进行Hilbert变换, 进而得到Hilbert谱。在对每一个IMF分量Cj (t) 作Hilbert变换之后, 得到一个变换平面内的数据序列C～j(t)C～j(t):

C～j(t)=1πP∫Cj(t′)t−tdtC～j(t)=1πΡ∫Cj(t′)t-tdt (9)

其中P为Cauchy主值, 由Cj (t) 和C～j(t)C～j(t)可以构成一个复序列Zj (t) :

Zj(t)=Cj(t)+iC～j(t)=aj(t)ejθj(t)Ζj(t)=Cj(t)+iC～j(t)=aj(t)ejθj(t) (10)

其中:

aj(t)=[C2j(t)+C～2j(t)]12θj(t)=arctanC～j(t)Cj(t) (11)aj(t)=[Cj2(t)+C～j2(t)]12θj(t)=arctanC～j(t)Cj(t)(11)

得到的瞬时频率为:

ωj(t)=dθj(t)dtωj(t)=dθj(t)dt (12)

因此, 原始数据序列可以表示为

X (t) =Re∑j=1naj(t)ej∫wj(t)dt∑j=1naj(t)ej∫wj(t)dt (13)

如果是在环境激励下而不是冲击荷载下采集到随机相应时程数据, 在对EMD法所得到的结果运用Hilbert变换前, 必须运用随机减量法 (Random Decrement Technique, RDT) 对各个IMF分量进行自由衰减响应的提取。该方法通过时间平均, 从随机振动信号中提取自由衰减响应。随机减量技术在使用上简单易行, 物理意义明确, 对激励只有定性要求, 且可利用随机扰动来激励, 同时只需对响应信号进行处理, 因此该方法在结构模态参数识别等许多领域中都得到了成功的百富策略白菜网[12,13]。

本文采用流固耦合松耦合计算方法, 联合百富策略白菜网经验模态分解法、随机减量法和希尔伯特变换法, 对文献[14]中大跨度柔性屋盖结构在四周封闭和突然开口两种情况下的气动阻尼和气承刚度进行研究, 同时探讨数值计算的准确性。

本次试验所采用的是能考虑结构和来流之间相互耦合作用的气动弹性模型。模型屋面为长宽均等于600 mm的正方形, 底裙高200 mm, 高跨比H/L=1/3, 如图5所示。采用1 mm厚的铝合金板作为屋面的材料, 弹性模量 E=7.1×1010 N/m2, 密度ρ=2 420 kg/m3, 泊松比γ=0.31。同时, 为了研究大跨度平屋面在四周封闭和迎风面开孔情况下的风振响应特征, 在模型迎风面正中开一100 mm×100 mm方孔。对于四周封闭模型, 方孔由一块同等大小的木板封闭。

风洞试验在南京航空航天大学603研究所的NH-2低速风洞中进行的。模拟B类地貌大气边界层, 地貌粗糙度系数α=0.16。试验共采用6只加速度传感器, 传感器在屋面板上的布置如图6所示。信号的采集与分析由靖江东华测试技术开发有限公司研制的DH5935/DH5936动态信号测试系统完成。封闭模型的采样频率为500 Hz, 开孔模型的采样频率1 000 Hz, 采样持续时间一般都超过90 s。试验风速为:8 m/s、9 m/s、10 m/s、11 m/s、12 m/s。

在数值模拟中按1∶1尺寸建立与风洞试验模型一致的计算模型。整个计算流域范围确定为宽8 m, 高3 m, 长度20 m的空间区域, 通过设置对称边界条件模拟整个大气环境并有效减小数值风洞的网格数量。屋面设置为动网格边界, 墙面和地面设置为固定网格边界。模型位于距离入口约6 m, 距离数值风洞出口约14 m。计算网格如图7所示。

非稳态计算选用大涡模拟。采用随机过程模拟方法得到脉动风速, 通过将风速作无散度化处理保证连续方程的成立, 进而作为入口风速。本文采用文献[15]中的谐波叠加法, 在数值风洞入口边界的每一单元中心位置处进行风速时程模拟, 程序中引入FFT算法以提高计算效率。

用上述模拟得到的时间序列还不能直接代入到实际的风压进行计算, 这是因为这些脉动风速的时间序列并不一定能保证连续性方程的成立, 因此, 必须将边界条件中的风速作无散度化处理, 即

u(S)n+1i=u(R)n+1i−1Δt1ρ∂ΔPn∂xiui(S)n+1=ui(R)n+1-1Δt1ρ∂ΔΡn∂xi (14)

1Δt1ρ∂2ΔPn∂x2i=∂u(R)n+1i∂xi1Δt1ρ∂2ΔΡn∂xi2=∂ui(R)n+1∂xi (15)

由泊松方程中得到的压力梯度修正项ΔPn代入到上述方程得到修正后的速度u(S)n+1ii(S)n+1, 为了避免流量的不平衡影响计算的收敛, 在实际计算中还使用了流量修正的方法, 即在计算的每个时间步上对入口处的总流量进行修正, 让流动的进出口流量保持一致。修正速度的大小为:

Δv=−(∑i=1nuiAi−S′)/AΔv=-(∑i=1nuiAi-S′)/A (16)

式中:ui是入口处每个单元上的脉动速度, Ai是该单元所代表的面积, A是风速入口处的总面积, S′是该时间步出口处的流量。入口处每一单元上的脉动速度经过无散度化后还需叠加上修正速度 Δv, 这样最后得到的脉动风速时程才可以代入到风荷载数值模拟计算的边界条件上去。

在入口边界的每个网格中心点处输入一条模拟得到的风速时程。同时将上述风速数据经过湍流度和流量修正, 得到真正实用的非稳态计算入口边界。图8所示为4 m高处一网格中心点位置经过湍流度和流量修正前后的风速时程和脉动风速谱的比较。从图中可以看出, 经过修正, 尤其是经过湍流度修正高处的脉动风能量出现明显下降。

有限元模型较为简单, 无需建立墙体等其他结构, 而仅需建立屋面模型。有限元模型前5阶模态与风洞试验模型基本吻合, 如表1所示。

如前文所述, 经验模态法可以提取时程信号前若干阶本征模函数C1 (t) ～Cn (t) , 其中第一个本征模函数C1 (t) 是从时程中分解出来的振幅最大、频率最高的波动, 依次下去的各内在模函数, 振幅逐渐变小、频率逐渐变低。

本文提取了测点4在封闭和开孔两种工况下下在8 m/s, 9 m/s, 10 m/s, 11 m/s, 12 m/s五种不同风速下的加速度时程响应数据, 联合运用HHT和RDT法, 可得到加速度时程的自由衰减曲线。对于封闭模型而言, 第一、二阶本征模函数还没能把主要信号提取出来, 而第三阶则成功提取了响应的典型频率。对于开孔模型, 第五阶本征模函数提取了响应的典型频率, 而从第六阶开始, 波动的振幅很小、频率极低, 可能是由于数据采样频率不够高和波动衔接等原因造成的噪声信号。限于篇幅, 在此不列出所有自由衰减时程曲线, 而由典型本征模函数运用RDT法得到的加速度时程自由衰减曲线如图10所示。

对图10所示的加速度时程自由衰减信号进行Hilbert变换, 可得到图11所示的加速度时程振幅的对数值lna的原始曲线。

在对lna进行线性拟合后, 可计算出两种计算工况在不同风速下屋盖结构的自然频率f和阻尼比ξ, 结果见表2。

从以上图表可以看出:① 数值计算结果的自然频率要高于风洞试验, 原因在于数值风洞模型是理想密封的 (即使对于迎风面开孔工况, 四周墙体是完全密闭的) , 而试验模型总存在着一定的气体泄漏, 结果导致结构气承刚度的降低。② 数值计算和风洞试验结果虽然在时程曲线上的吻合度不是非常准确, 但lna曲线的斜率即阻尼比的结果还是吻合得很好, 这说明采用本文提出的流固耦合计算方法可以较准确地计算出结构风振阻尼。③ 开孔状态各个风速下的阻尼要比封闭状态时的大, 这是由于开孔工况屋面振幅更大, 从而使得风振的气动阻尼也更大。因此对于发生大变形的结构振动而言, 气动阻尼的影响不可忽略。④ 由于气动阻尼比ξa与结构阻尼比ξs之和为总的阻尼比ξ, 所以, 结构阻尼比ξs和总的阻尼比ξ已知, 就可以得到气动阻尼比ξa, ξa=ξ-ξs。

(1) 根据三角形重心坐标系原理提出了基于松耦合计算的流固耦合网格协调插值新方法, 提高了流固耦合的计算效率。

(2) 为解决由于柔性结构加速度过大而导致的求解不稳定问题, 通过在1个时间步上进行流体及结构间反复交替求解来得到稳定的时程解。

(3) HHT变换现已被证明是一种很好的方法而运用于结构气动阻尼的确定。本文通过HHT变换, 很好地将振动的干扰模态分离出来, 并结合RDT法, 得到了理想的自由衰减曲线。

(4) 计算结果显示四周封闭建筑的阻尼比远小于迎风面开孔建筑的阻尼比, 且四周封闭状态下气动阻尼随风速变化小, 其对总阻尼的影响较小;而迎风面开孔状态下气动阻尼较大且随风速的增大而增大。因此, 对于柔性且阻尼较小的结构, 特别对于开孔结构, 气动阻尼的影响不可忽略。