流固耦合问题是研究变形体在流场作用下的各种行为以及固体变形对流场影响的一门学科.流固耦合力学的重要特征是两相介质之间的相互作用, 即变形体在流体载荷作用下会产生变形或运动, 而变形或运动又反过来影响流场, 从而改变流体载荷的分布和大小.正是这种相互作用将在不同条件下产生形形色色的流固耦合现象.最近的几十年中, 人们逐渐地发现流体与结构的相互作用在许多土木工程问题中都起着非常重要的作用.例如悬索桥、高层建筑、输电塔以及大屋面的轻型索、膜结构等.其中索、膜结构等三维体系的流固耦合问题十分复杂, 这方面的研究大大落后于桥梁结构方面的类似研究.文献[1]提出了一种简化的气动弹性模型, 用一些系数来模拟流体与结构之间的气动力作用, 同时还引

入了附加质量、气动阻尼和气承刚度等新概念.这些代表流固耦合作用的参数或系数需由实验或半理论半实验的方法来确定.目前, 这种简化理论还只能初步解决极其简单的问题, 离实用尚远.现今广泛使用的是耗时费钱的风洞试验, 而风洞试验的一个主要缺陷就是很难快速地改变结构模型的几何图形, 以对不同的结构形式进行分析.因此非常需要一种有效的数值风洞, 使其能够直接运用到土木工程的设计和分析过程中.于是, 借鉴桥梁、飞机等工程领域的经验, 综合运用计算流体动力学和计算结构动力学的方法, 开发适用于大跨柔性结构的数值风洞技术, 以便对流体-结构的耦合效应进行更精确的数值模拟, 无疑是一种有前途的研究方法.这种方法的理论难度较大, 但十分吸引人, 已成为当前的前沿研究领域, 也已取得了一定进展.

从总体上来看, 流固耦合问题按其耦合机理可分为两大类:直接耦合和分区耦合.第一种方法是通过单元矩阵或荷载向量把耦合作用构造到控制方程中, 然后对控制方程直接求解, 在构造控制方程过程中常常不得不对问题进行某些简化, 其优点是概念比较清晰, 缺点是计算准确程度较难保证[2].第二种方法是分别对流体动力方程和结构动力方程依次求解, 通过把第一个物理场的结果作为外荷载施加于第二个物理场来实现两个场的耦合, 其优点是可以根据不同物理场的特点分别制定求解方案, 计算准确程度较高, 缺点是计算过程比较复杂[3,4].直接法比较适用于对耦合场的理论分析, 分区法比较适用于对耦合场的数值计算[5].

为实现对薄膜结构与风之间耦合作用的数值模拟, 本文作者拟采用分区耦合算法, 即将流体计算与结构计算分开处理, 通过一些参数的传递来实现两个分区的耦合.

对于流体, 要求解的质量守恒方程和动量守恒方程为

∂ρ∂t+∂∂xi(ρui)=0   (1)∂∂t(ρui)+∂∂xj(ρuiuj−τij)=−∂p∂xi+si   (2)$\begin{array}{l}\frac{\partial \rho }{\partial t}+\frac{\partial }{\partial x{}_i}(\rho u{}_i)=0(1)\\ \frac{\partial }{\partial t}(\rho u{}_i)+\frac{\partial }{\partial x{}_j}(\rho u{}_{i}u{}_j-\tau {}_{ij})=-\frac{\partial p}{\partial x{}_i}+s{}_i(2)\end{array}$

式中, t为时间;xi, xj为笛卡尔坐标 (i, j=1, 2, 3) ;ui, uj分别为流体在xi, xj方向的绝对速度分量; p为压力;ρ为密度;τij为应力张量分量;si为动量源项分量.

采用“求解压力耦合方程组的半隐式方法 (SIMPLE算法) [6]”.该方法所得到压力修正方程为

ai,jp′i,j=ai+1,jp′i+1,j+ai−1,jp′i−1,j+  ai,j+1p′i,j+1+ai,j−1p′i,j−1+b′i,j   (3)$\begin{array}{l}a{}_{i,j}{p}^{\prime }{}_{i,j}=a{}_{i+1,j}{p}^{\prime }{}_{i+1,j}+a{}_{i-1,j}{p}^{\prime }{}_{i-1,j}+\\ a{}_{i,j+1}{p}^{\prime }{}_{i,j+1}+a{}_{i,j-1}{p}^{\prime }{}_{i,j-1}+b^{\prime }{}_{i,j}(3)\end{array}$

式中, ai, j为给定系数;p′i, j 是压力修正值; b′i, j 是由于不正确的速度场所导致的“连续性”不平衡量.

借助隐式时间积分方法, 在每个时间步内进行迭代, 直到取得本时间步的收敛解, 然后转入下一时间步继续重复上述过程.

湍流流动是一种高度非线性的复杂流动, 但人们已经能够通过某些数值方法对湍流进行模拟, 取得与实际比较吻合的结果.关于湍流的数值计算, 这里主要采用Launder和Spalding于1972年提出的基于Reynolds平均法的标准k-ε模型.

标准的k-ε方程中含有两个方程, k方程和ε方程, 即

ρ∂k∂t+ρ∂(kui)∂xi=∂∂xj[(μ+μtσk)∂k∂xj]+  Gk+Gb−ρε−YM+Sk   (4)ρ∂ε∂t+ρ∂(εui)∂xi=∂∂xj[(μ+μtσε)∂ε∂xj]+  C1εεk(Gk+C3εGb)−C2ερε2k+Sε   (5)$\begin{array}{l}\rho \frac{\partial k}{\partial t}+\rho \frac{\partial (ku{}_i)}{\partial x{}_i}=\frac{\partial }{\partial x{}_j}\left[\left(\mu +\frac{\mu {}_t}{\sigma {}_k}\right)\frac{\partial k}{\partial x{}_j}\right]+\\ G{}_k+G{}_b-\rho \epsilon -Y{}_{{\rm M}}+S{}_k(4)\\ \rho \frac{\partial \epsilon }{\partial t}+\rho \frac{\partial (\epsilon u{}_i)}{\partial x{}_i}=\frac{\partial }{\partial x{}_j}\left[\left(\mu +\frac{\mu {}_t}{\sigma {}_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x{}_j}\right]+\\ C{}_{1\epsilon }\frac{\epsilon }{k}(G{}_k+C{}_{3\epsilon }G{}_b)-C{}_{2\epsilon }\rho \frac{\epsilon {}^2}{k}+S{}_{\epsilon }(5)\end{array}$

式中, k为湍动能;ε为耗散率;Gk为由于平均速度梯度引起的湍动能k的产生项;Gb为由于浮力引起的湍动能k的产生项;YM为可压湍流中脉动扩张的贡献;C1ε、C2ε和C3ε为经验常数;σk和σε分别为与湍动能k及耗散率ε对应的Prandtl数;Sk和Sε为用户定义的源项;μ是动力粘度;μt是湍动粘度.

另外, 对于所有的CFD问题都需要有边界条件.所谓的边界条件, 是指在求解区域的边界上所求解的变量或其导数随地点和时间的变化规律.只有给定了合理边界条件的问题, 才可能计算得出流场的解.因此, 边界条件是使CFD问题有定解的必要条件.针对本文的具体情况, 所采用的边界条件为:

(1) 流动进口边界条件.采用速度进口边界, 用来定义流动进口处的风速以及湍流参数湍动能k及耗散率ε.

(2) 流动出口边界条件.采用压力出口边界, 在流动分布的详细信息未知, 但边界的压力值已知的情况下, 使用压力出口边界条件.取参考压力为一个大气压0.1 MPa.假定垂直于出口方向的速度分量是零梯度的.

(3) 对称边界条件.指所求解的问题在物理上存在对称性, 应用对称边界条件可避免求解整个计算域, 从而使求解规模缩减至整个问题的一半.在对称边界上, 垂直边界的速度取为0, 而其他物理量的值在该边界内外是相等的, 即计算域外紧邻边界节点的值等于对应的计算域内紧邻边界节点的值.

(4) 壁面边界条件.限定整个流场的求解区域.

膜结构所使用的主要材料本身不具有基本刚度和形状, 即在自然状态下不具有承载能力, 只有对膜材施加预应力才能获得结构承载所必需的刚度和形状, 因此其设计过程首先是找形分析.先给定结构内预应力分布状态, 然后求满足该预应力状态和边界条件的结构形状.膜结构在荷载作用下, 尽管应变很小, 但位移较大, 所以在建立平衡方程时就必须考虑结构非线性变形的影响.

薄膜结构的非线性动力响应运动方程为

Mu¨(t)+Cu˙(t)+Ku(t)=F(t)   (6)$\mathit{{\rm M}}\ddot{\mathit{u}}(t)+\mathit{C}\dot{\mathit{u}}(t)+\mathit{{\rm K}}\mathit{u}(t)=\mathit{F}(t)(6)$

式中, M是质量矩阵;C是阻尼矩阵;K是刚度矩阵;F是由流体导致的作用在结构上的外力;u˙$\dot{\mathit{u}}$ (t) , u¨$\ddot{\mathit{u}}$ (t) 是速度和加速度.

耦合采用的是分区耦合求解的方法.流体域模拟是由基于有限体积法的CFD软件完成, 用两方程k-ε模型求解RANS (Reynolds-Averaged-Navier-Stokes) 方程[6].对于结构域的找形和求解在CSD软件ANSYS中完成, 其动力非线性响应 (大位移/小变形) 采用的是基于有限元的方法和隐式迭代求解方法.计算过程通过在CFD和CSD之间进行迭代求解, 直到在每一个荷载步当中都收敛.

在每一次迭代中, 把经过CFD计算出来的结果 (主要是压力) 作为CSD计算的输入节点荷载.CSD计算后得出结构的位移, 根据这个位移再进行流体域中结构的几何边界的修正.如果结构是大变形 (如本文中的膜结构) , 可以采用低松弛技术进行结构边界几何图形的校正, 从而提高收敛性.对于膜结构, 通过找形分析确定在预应力作用下的表面形状.结构被划分为非结构化的三角形网格, 再将其几何信息通过一定的数据格式输出, 导入到CFD中.对于CFD, 也可以通过相应的网格处理器, 将膜划分为三维的六面体结构化网格.图1显示相应的薄膜结构的结构域网格及流体域网格.流体域和结构域的计算分别在各自的处理器中进行, 用MPCCI耦合界面管理ANSYS和STAR-CD之间的数据交换.

另外, 由于CFD计算采用的是四边形的有限体网格, 而CSD计算采用的是三角形的有限元网格, 因而两种网格并不匹配, 他们之间的节点和单元数据的传递交换是通过建立一个中间的几何模型实现的, 这个中间的几何模型是以结构的表面尺寸为原形, 在流场中建立一个称作“挡板 (baffle) ”的单元.尽管两种网格都绘制出同样的表面, 但他们的节点并不一致.因此需要进行一些相互间的插值 (图2) .

需要特别指出的是, 在每一次的流体-结构耦合迭代过程中, 流体域的网格不断进行更新改变, 以适应新的“挡板”单元的边界层位置.这个过程通过用户子程序采用动网格技术实现.

为验证所阐述的耦合算法的可行性, 本文对一鞍形膜结构进行了风振响应分析.所采用的风速为v=15 m/s.

膜结构的边长为10 m, 弹性模量为900 MPa, 泊松比为0.47, 密度为1 500 kg/m3, 厚度为1 mm, 预应力为20 MPa.湍流模型采用高雷诺数的k-ε模型, 湍流密度取0.01, 进行流固耦合计算分析.限于篇幅, 取6个不同时刻膜单元表面压力分布如图3所示.

由于数值风洞方法要涉及到多学科知识的综合运用, 其理论和实际操作中的困难是显而易见的.目前国内外关于膜结构数值风洞的研究基本上是空白, 作者在这方面的研究也刚刚起步.本文采用分区耦合的方法, 应用CFD数值模拟技术来求解膜结构与风荷载的流固耦合问题, 结果证明这种方法是完全可行的, 可以得到在不同的时刻, 膜结构表面的风压分布变化值, 可以克服气动弹性模型和风洞试验方法所带来的一些缺点.需要指出的是对于膜结构风工程这样复杂的钝体绕流问题还缺乏成熟的湍流计算模型, 定量化的数据还有待于进一步研究给出.

应用CFD数值模拟技术来求解膜结构的流固耦合振动问题无疑是一种有前途的研究方法.但所面临的困难也是巨大的.本文通过对一工程实例的数值模拟, 给出了膜结构随时间风压分布的变化, 验证了这种耦合算法的可行性以及稳定的收敛特性.但是还有很多问题期待进一步的解决.例如, 流体域中网格随时间的改变即动网格技术可以进行更好的改善, 使其具有普适性.而对于结构的模拟, 目前只能应用于类似膜这样的薄壳结构, 应该将这种方法拓展成三维的实体单元.但从这里已经可以看到数值风洞技术在结构风工程领域应用的广阔前景, 为下一步进行数值模拟和实验观测数据的比较奠定了基础.