膜结构是近年国内百富策略白菜网较多的一种新型空间结构形式, 该类结构因跨度大、结构轻巧、形式多变而逐步得到推广百富策略白菜网;但同时也表现出较为明显的轻柔特点, 即膜结构在风荷载作用下会产生较大的变形和振动, 通常可达到结构跨度的1/50左右, 大幅的变形和振动反过来会影响结构表面的风压分布情况, 从而形成较为明显的流固耦合效应, 这也是膜结构抗风设计中的关键和难点问题[1,2].

为研究膜结构与气流之间的耦合振动问题, 通常可借用的方法主要有3类:一是结构物理风洞, 进行缩尺比下的模型风洞试验;二是借用高性能计算平台, 实现膜结构风振的数值风洞模拟;三是在实际膜结构工程上进行膜面的风振响应监测研究.卫东等[3]进行了较大足尺张拉膜结构的荷载试验, 重点研究了张拉膜结构振动与空气耦合引起的附加质量优先及其阻尼比值的变化规律;张其林等[4]进行了世博轴膜结构的气弹模型风洞试验, 对其进行了不同风场下的风振响应特性分析;王彬等[5]进行了膜结构风致响应的分区耦合算法研究, 结合ANSYS和CFD软件计算平台, 初步实现了鞍形膜结构的风振响应数值模拟, 但未给出膜结构风振过程中响应的量化指标;黄国辉等[6]利用非线性随机振动理论进行了张力膜结构时域的风振响应研究, 获取了张力膜结构的风激动力响应特性以及风振响应与结构参数的关系.

本文通过结合现有成熟商业软件, 初步建立计算平台, 实现了结构与流场的相互耦合计算, 从而成功实现张拉膜结构在风荷载作用下的气弹响应.

随着计算机硬件的飞速发展和数值计算软件的进一步完善, 采用数值模拟方法进行膜结构的风振响应分析越来越凸显其优势, 相对风洞试验只能局限于特定的结构形式, 数值模拟则可较为方便地实现各类结构形式以及不同结构参数下的模拟, 大大提高了研究的范围和效率.为实现膜结构风振耦合响应的数值模拟, 需重点解决好两个问题:首先是流场域的模拟和结构动力学的模拟 (CSD) , 目前该部分的数值模拟技术已相对较成熟;其次是流场域和结构域界面的匹配及物理信息的交换, 尤其在膜结构产生较大位移的前提下, 如何保证流场域与结构域边界条件的协调一致是耦合数值模拟的难点.

膜结构所使用的膜材料本身不具有刚度和形状, 只有对膜面施加张力才能获得结构承载所需的刚度和形状, 因此, 结构的求解为先施加结构的预张力, 然后迭代求解获得满足该预应力状态和边界条件下的膜面形状.膜结构非线性变形的振动控制方程为

对于工程领域的结构绕流场模拟, 通常属于大雷诺数下的湍流流动, 由于目前直接进行湍流场的模拟还不现实, 故有必要通过选取合适的湍流模型来考虑气流的湍流效应.常用的湍流模型有基于时间平均的雷诺平均模型 (RANS) 和基于空间平均大涡模型 (LES) .大涡模拟直接模拟湍流的大尺度涡旋, 同时对小尺度涡旋通过模型加以考虑, 该方法能够获得比传统RANS更好的模拟结果.大涡模拟直接求解过滤后的流体控制方程为[7]

根据Smagorinsky的标准SGS模型, 亚格子应力具有如下形式:

式中:

ΔxΔyΔz表示网格沿x、y、z轴的尺寸;Cs为常数, 取为0.1.

求解流固耦合问题的困难在于两相界面处于不同运动描述方式下的相互协调问题, 通常对于结构运动描述采用Lagrange方法, 网格点随结构一起运动;而针对流体的运动描述采用Euler方法, 网格点不随流体运动.为使得结构运动与流体运动协调一致, 最有效的方法是采用ALE (Arbitrary Lagrang-ian-Eulerian) 方法, 该方法通过引入流场域的网格速度项来考虑由于边界变形而导致流场域网格的更新.由于膜结构在风荷载作用下产生较为明显的变形, 故导致在新的边界条件下流场域的网格必须随着边界同步更新.如何快速高质量地获取新边界条件下的流场域网格划分成为流固耦合分析的关键.

动态边界条件下流场网格的更新通常采用弹簧法, 该方法将流场网格的每条边都等效为具有一定弹性系数的弹簧, 当边界移动后, 弹簧系统能够反馈边界的移动并传递到流场各个位置从而获得新边界下的网格划分.已有较多商用软件处理流固耦合问题时也常采用该方法, 该方法简便易用, 但选择约束不当在处理边界的大变形时易出现网格过小的夹角或翻转, 导致流场求解的不收敛, 限制了该方法百富策略白菜网范围, 特别是大位移大变形条件下的网格重生成难题一直限制了软件的百富策略白菜网范围.因此, 本文对该方法进行一定的优化, 将弹簧刚度的表达形式优化为如下形式[8]:

式中:Kij为网格的等效弹簧刚度;d为距离改进系数, 主要保证边界附近的网格不发生大的变形;Vmin为单元体积改进系数, 主要保证网格划分不出现负体积;sinθ为角度改进系数, 主要保证网格不出现过小的夹角;lij为网格的边长;α, β, γ为加强因子, 本文取为2.

本文建立了膜结构气弹响应数值模拟的计算平台, 其中结构计算模块和流体计算模块分别调用商业软件ANSYS、CFX完成, 流场域与结构域边界和动态网格的更新通过自编的平台完成.该平台能够实现大位移边界下的网格动态更新, 能够匹配结构与流场的边界一致, 实现结构与流场的弱耦合求解.

为验证该平台数值模拟的可行性, 本文对工程中常用的鞍形膜结构进行了风振响应分析.膜结构边长为10m, 弹性模量900MPa, 泊松比0.47, 密度800kg/m3, 厚度1mm, 预应力2~4kN/m, 鞍形膜空间布置如图1所示, 有限元结构模型如图2所示.

在膜结构周围建立一个空间的三维数值风洞区域, 尺寸为600m×200m×80m;入口处脉动风速的生成采用谐波合成法获得, 风速谱采用Kaimal谱, 入口处湍流度强度取为15%;膜面及地面采用无滑移壁面边界, 侧面及顶面设置为对称边界, 出口设置为无压力梯度边界, 如图3所示;计算时间步长取为0.01s, 计算时长150s, 膜结构风振过程的风压系数分布变化云图如图4所示, 由图可见, 气流沿着迎风面膜面尖角形成旋涡脱落的发展轨迹.

图5为风速20m/s、风向角为0°时, 膜面中心位置处的位移响应曲线及对应的频谱图, 由图可见, 膜面风振位移由多个模态组成, 说明膜结构的振动是由多模态耦合叠加而成.大跨空间结构不同区域的风振系数差别较大, 由于张拉膜结构的整体结构形式较为简单, 故本文将整体结构只划分为一块进行整体风振系数的定义, 该系数的定义通过将节点位移的响应作为权值进行统计平均, 以方便进行不同模拟工况下的结果对比, 针对整体膜结构的位移定义一致的风振系数为

式中:σi为膜面节点位移响应的标准差;g为位移响应峰值因子, 取为2.5;ui-为膜面节点位移响应的平均值;β为膜结构的整体风振系数.该表达式对上下同乘节点位移的平均响应进行统计分析, 是为了将节点平均位移作为权值因子引入统计中, 使得平均位移大的区域对整块风振系数的贡献更大, 同时也避免某些区域平均位移很小而出现统计风振系数异常大的情况[9].

张拉膜结构风振系数随风速变化的相关曲线如图6所示, 由图可见, 风振系数随风速的增加略有增大, 基本呈线性增长关系, 说明随着风速的增大, 膜面位移随之增大, 结构的几何非线性表现得更为突出.张拉膜结构风振系数随风向角变化的相关曲线如图7所示, 由图可见, 风振系数为0°、90°时较大, 而在45°时相对较小, 说明张拉膜在斜风作用下风振响应相对较小, 其原因在于斜风作用下气流更易沿着膜面运动而产生相对较小的漩涡脱落, 从而降低了膜面的风振响应.

张拉膜结构风振系数随预张力变化的相关曲线如图8所示, 由图可见, 风振系数基本随预张力的增加呈线性增长关系, 预张力的增大使得膜面的刚度增加, 但张拉膜结构整体的风振响应减小, 原因在于风振系数的定义是极值响应与均值响应的比值, 说明在低张力作用下, 膜结构几何非线性更明显, 从而产生相对更大的平均响应, 降低了整体结构的风振系数.张拉膜结构风振系数随矢跨比变化的相关曲线如图9所示, 由图可见, 风振系数随矢跨比的增加而略有增大, 随着矢跨比的增加, 张拉膜外形更趋于平缓, 膜面整体的风振响应减小, 说明低矢跨比下张拉膜结构在风荷载作用下表现出了更强的几何非线性, 平均位移也相对更为明显, 从而降低了结构的风振系数.

本文通过结合现有成熟商业软件, 初步建立计算平台, 实现了结构与流场的相互耦合计算, 从而成功实现张拉膜结构在风荷载作用下的气弹响应.证明了利用数值模拟技术求解膜结构的流固耦合问题是完全可行的, 可补充风洞试验的一些缺陷, 大大扩展了膜结构风致振动的研究手段和效率.

典型张拉膜结构的风致振动响应数值模拟结果表明, 膜结构在气流作用下的振动是多个模态的耦合风致振动;同时由于结构位移响应较为明显, 振动过程也表现出了较为明显的几何非线性而不是简单的线弹性风致振动.风速、风向角、预张力和矢跨比等因素都会对张拉膜结构的风致振动产生一定的影响, 其变化规律需通过具体数值模拟分析, 张拉膜结构整体风振系数通常在1.6~1.8.

百富策略白菜网CFD技术求解膜结构的流固耦合问题无疑是一项有前途但又面临巨大挑战的工作, 无论是在高雷诺数下湍流的模拟, 还是大位移运动边界下的网格动态更新, 或是数值求解的高效高精度收敛, 都有待进一步研究.