流固耦合现象在斜拉桥、超高层建筑、塔桅结构、输电线、索膜结构等土木工程领域扮演着相当重要的角色[1]。过去由于缺乏对流固耦合现象的认识，曾发生过1879年苏格兰Tay大桥的垮塌，1948年西雅图Tacoma大桥风振倒塌，以及1965年英格兰三座冷却塔倒塌等重大灾难[2]。如今，一方面诸如弛振、颤振等风致结构气动现象已为人们所熟知。另一方面由于大量轻盈结构的使用 (例如亚特兰大七万人体育场屋顶采用膜结构，同时结构跨度达到240m) ，建筑对风的响应越来越敏感。为了准确获得具有外形复杂或是受外界风场影响显著的建筑结构的风荷载，工程界已经进行了大量费用高昂的风洞试验或是采用半经验方法对流固耦合问题进行分析计算。此外，工程界已经可以通过使用相对成熟的有限元分析方法 (FEM, Finite Element Method) 对结构在风荷载下的应力、变形等效应进行精确计算。同时，也能够通过使用计算流体动力学 (CFD, Computational Fluid Dynamics) 软件对结构所受的风荷载进行准确估计。然而，对于进行流固耦合计算的方法和工程百富策略白菜网却少之又少。

松耦合方法也称弱耦合方法，是指通过流体计算得到风压、风力并将其输入到动力计算中，求得位移、速度及加速度响应。然后，基于这些响应在流体计算中对建筑物周围的计算网格进行移动从而改变建筑物表面的边界条件。如此，流体的控制方程及建筑物的控制方程可以独立并交互求解[3]。最初也是最容易想到的流体与固体之间的交互程序是使用所谓的一致网格，即固体表面的网格划分与流体与固体交界面处的流体网格划分完全一致。这样无需通过插值计算，结构上每个节点的变形都可以完全传递给流体节点;而流体中计算得到的流固交界面之间的荷载数据也可以无损地传递到结构表面。然而，使用这种方法的缺点是需要占用较大的计算机存储资源。

事实上由于计算方法和精度要求的不同，有限元模型和流体力学模型对流固耦合面上网格数量的要求往往相差很大。流体力学模型通常需要非常精细的网格才能捕捉得到诸如边界层、激振等现象;而结构力学所要求的网格精度往往没那么高。此外，不同领域的求解器对网格的类型往往也非常不一致，这就出现了网格节点之间的匹配问题，见图1ㄢ

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

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

整理式 (1) 得到:

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

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

求解式 (4) ， (5) 得到:

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

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

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

作为2010年上海世博会立体交通组织重要载体的世博轴，由两层地下空间和一层高架步道组成。高架步道上覆连续张拉式索膜结构屋顶，总长约840m，最大跨度约97m，总面积近64 000 m2。包括支点系统和膜面系统两部分。

支点系统通过水平索将31组外桅杆及背索、19个下拉点、18个与阳光谷的拉节点连接成稳定体系。膜结构轻盈飘逸，见图4，膜面采用国家A级PTFE膜材，主要由脊索、边索、谷索和膜形成连续的三角倒锥形膜单元。

世博轴地处台风多发的上海地区 (50年重现期基本风压为0.55kN/m2) ，项目现场位于黄浦江畔，周边地形开阔无高大建筑阻挡，属C类地貌。膜结构最大标高38m，平均标高超出10m平台约15m，属于低矮房屋。对于此类建筑，我国相应规范中既无体型系数可供参考，更无风振系数的计算方法。本文计算基于通用有限元和流体计算软件，通过二次开发，采用内存为8G的16核64位工作站进行以上模型的计算。

形状确定是膜结构设计的第一步，也是决定一个膜结构设计是否成功的关键。相比于普通小品膜结构而言，世博轴索膜结构具有跨度更大、结构更柔的特点，采用有限元方法分析存在收敛困难等问题。文献[7]采用分步计数、逐步刚度逼近的方法成功解决这一问题，并进一步考察了不同结构自重、不同弹性模量、泊松比对找形及其后续计算的影响，此处不再赘述。有限元模型如图5, 6所示，取具有代表性的3片膜进行讨论，膜面总面积约为9 300m2ㄢ

为考虑周边建筑对膜面风荷载的影响，数值风洞模型包含了阳光谷、10m平台等其他建筑。为提高计算效率，仅对有限元模型所对应的膜面设置动网格属性，其余部分保持“不动”。非稳态计算采用LES模拟[8,9]、隐式计算。简单的弱耦合算法会在膜面上生成过大的加速度从而导致求解不稳定，通过在1个时间步上进行流体与结构间反复交替求解来得到稳定的时程解。

为了研究不同网格体系的计算精度和效率，本文设计了以下3套网格工况:工况一:有限元模型和流体力学模型均采用具有3 924个单元的一致网格体系;工况二:有限元模型和流体力学模型均采用具有11 935个单元的一致网格体系;工况三:有限元模型采用具有3 924个单元的稀疏网格，流体力学模型采用具有11 935个单元的稠密网格，两者之间采用上文所述的变形协调系统传递数据。不同密度网格系统如图7所示。

CFD计算采用通用流体力学计算软件Fluent，在来流为10m/s的均匀流条件下，A点的位移见图8。计算开始的初始1min结构并不参与计算，即网格不动，直至流体计算稳定后才释放结构使网格发生运动。因此，结构在被释放后将有一个脉冲运动，见图8 (图中起始时刻为计算第2min) 。

计算表明，当单元尺寸小于0.1m时计算结果与网格尺度无关，A点处的风压为0.03Pa。本算例中稀疏网格的单元尺度基本在5m左右，最大单元尺度接近8m;而稠密网格的单元尺寸在0.5m左右。非稳态计算稳定时刻CFD模型的稀疏网格和稠密网格在A点处的风压发现，前者在A点处的风压为0.04Pa，后者为0.032Pa，相对误差分别为33%和7%。由此可见，CFD网格密度对风压计算结果影响较大。

对于有限元模型而言，单元数量对索膜结构有限元计算速度有着非常重要的影响，随着有限元单元数量增多，非线性计算的收敛性变差，需要更多的迭代步和更小的步长，从而导致了计算速度的下降。另一方面单元数量对于计算精度提高有限。表1为具有不同单元数量的有限元模型模态计算结果和耗费的机时。因此，为了既保证数值计算的准确性，同时也为了提高计算效率，实际计算中可以采取有限元模型为稀疏网格，CFD计算采用稠密网格的策略。这也说明了采用非一致网格的必要性。本文采用稀疏一致网格、稠密一致网格和非一致网格对世博轴索膜结构算例进行时长为2min的动力时程计算，计算机实际耗费的机时见图9ㄢ

图9中列出4种工况:工况A代表FEM单元数量和CFD单元数量均为3 924的一致网格;工况B代表FEM单元数量为3 924和CFD单元数量为11 935的非一致网格;工况C代表FEM单元数量和CFD单元数量均为11 935的一致网格;工况D代表FEM单元数量为3 924和CFD单元数量为11 935的非一致网格。工况B和工况D的区别在于，工况B采用在计算初始化时就已经设定的FEM单元与CFD单元之间的对应关系，在之后的每一次网格更新计算中无需重新建立这种关系;而工况D却没有采取这一初始化操作，它需要在每一次网格更新计算时，重新建立一次非一致网格间的对应关系。

从工况A与工况C的比较发现，单元数量对索膜结构有限元计算速度有着非常重要的影响。这是由于随着有限元单元数量增多，非线性计算的收敛性变差，需要更多的迭代步和更小的步长，从而导致了计算速度的下降。从工况B与工况D的比较发现，非一致网格间对应关系的建立是一个非常耗费机时的运算。要在3 924个FEM单元和11 935个CFD单元之间建立一次对应关系就需要进行一次3 924×11 935次运算。因此，在初始化中即建立好非一致网格间的对应关系对于提高流固耦合的计算效率非常重要。

比较了考虑与不考虑流固耦合作用，对世博轴索膜结构风振响应的影响。不考虑流固耦合作用风荷载的取值是指:1) 利用依据达文波特谱模拟得到的风速时程作为CFD计算的入口条件;2) 采用非稳态计算获得膜结构表面的风压时程;3) 将第2步中获得的风压时程直接加载到结构有限元模型上进行结构的动力时程分析，从而获得结构的位移、应力等响应时程。图10为不考虑流固耦合作用膜面上A点自由衰减曲线。从图中可知，振动受到结构非线性阻尼的影响。在振动的前60s，由于振幅较大，非线性阻尼也较大，计算阻尼比约为0.15%。随着运动幅值的减小，60s之后的阻尼比也变得非常小。振动频率约为0.6，与膜结构的基频接近。

采用同样的方法可以计算出图8中考虑流固耦合的膜结构在自由衰减过程中的总阻尼比为4%，远大于不考虑流固耦合时的阻尼比。振动频率约0.14Hz。由此可见，流固耦合作用使膜结构的动力特性发生了较大的变化。

采用来流为10m/s的紊流计算膜面的风振响应。图11为不考虑流固耦合和考虑流固耦合作用膜面上A点的位移曲线。可以看出，索膜结构的风致响应可以分解为3个部分:静态成分、准静态成分和瞬态成分[10]。所谓静态成分是指可以直接利用CFD稳态模型求解的、不受脉动风场影响的结构表面平均风压;准静态成分是指主要考虑脉动风的空间相关性对结构整体振动的影响，可以利用CFD技术和多阶模态力法求得与结构最不利变形模态对应的风荷载分布模式;瞬态成分是指需由CFD非稳态计算得到的，考虑实际风场脉动影响及结构本身动力特性的风荷载成分。其中静态耦合对应于风荷载的平均响应，准静态耦合对应风荷载的背景响应，瞬态耦合对应风荷载的共振响应。表2列出了考虑与不考虑流固耦合作用的膜结构风振响应中各成分的分量。

从表2数据可以看出:1) 静态耦合反映了几何形状的变化对结构分析结果的影响，本算例中结构的几何变形使得A点位移趋于增大;2) 准静态耦合反映了脉动风中大尺度涡对结构响应的影响，从本例中可以发现，流固耦合基本不影响大尺度涡对结构的作用;3) 瞬态耦合反映了脉动风的高频部分 (小尺度涡) 对结构响应的影响。从计算结果来看，风荷载的高频部分对膜结构的风振响应影响很弱以至于可以忽略不计。

实际工程中，设计人员往往喜欢用风振系数的概念来描述结构的风致动力效应。如果采用位移风振响应的定义来比较A点的风振系数，则可得到:

式中:β0和βc分别代表A点不考虑和考虑流固耦合作用的风振系数;Dmax表示位移最大值;Dstead表示平均风荷载下的位移。

可知，对世博轴索膜结构而言，如果不考虑流固耦合效应，而采用普通的动力时程计算，结果将偏于不安全。

(1) 基于三角形重心坐标系原理实现非一致网格的变形协调，该方法由于仅包含少量的四则运算，因此计算速度非常快。此外，通过初始化非一致网格间的对应关系能进一步提高该方法的计算效率。

(2) 静态耦合和准静态耦合对世博轴索膜结构的流固耦合效应影响较大，而瞬态耦合的影响非常小以至于可以忽略不计。

(3) 相比于不考虑流固耦合作用的动力时程分析，静态耦合使结构的响应增强，而准静态响应与不考虑流固耦合作用的数值相同。

(4) 对于索膜结构等大跨度柔性结构而言，流固耦合作用不容忽视，采用常规计算有可能导致计算结果偏于不安全。