大跨度桥梁属于轻质柔性结构，是现代桥梁中重要的建筑形式之一。由于大跨度桥梁的特点，风荷载是其主要的控制荷载。自1940年塔科马大桥遭遇风荷载破坏以来，流固耦合 (简称FSI)作用就成为了大跨度桥梁设计中不可或缺的一部分^[1]。并且随着建筑科技的不断进步，大跨度悬索桥越来越多地被应用到实际工程中，由于桥梁的跨度的不断增大，柔性也在不断增加，因此对于风与桥梁的流固耦合作用的研究就显得尤为重要，研究流固耦合作用对于准确判断桥梁的风致涡激振，以及其他气弹现象，如颤振、抖振、驰振等桥梁振动有着重要的理论意义和工程价值。

目前国内外对于大跨度桥梁的风振响应的数值模拟研究有一些，但对于流固耦合的研究还很有限。Frandsen^{[2, 3]}开发了基于有限元方法具有动网格技术的程序SPECTRUM，并对大带东桥的涡激振动与颤振进行了全耦合求解，获得了与其他研究有很好一致性的结果。Taylor^[4]釆用分离耦合方法研究了 Storeboelt大桥的颠振行为，结果与风洞试验一致。Selvam^[5]开展了基于有限元的流固耦合模型对大带东桥的颠振临界风速进行了分析。Sarwar^[6]釆用大涡模拟(LES)研究了矩形箱梁在加装不同抗风构造后对风致耦合激振动的抑制作用，最终验证了数值模拟结果的可靠性。国内学者徐讽^[7]等使用CFD方法数值模拟了低雷诺数下不同截面的柱体受到弹性支撑时的振动现象，并使用动网格技术来实现柱体与流场间的非线性耦合作用，研究了不同截面形状对柱体振动形式的影响。在风致祸激振动方面，陈文礼等^[8]将圆柱看作一个单由度的弹簧振子系统，采用基于RANS方法的剪切应力输运模型(简称SST模型)，考虑了圆柱和流场之间的双向耦合作用，观测到了锁定现象和“拍”的现象。可以看出，考虑流体和结构耦合的相关研究较为有限。而考虑流固耦合作用的风振响应分析由于引入了计算结构力学的计算，突破了传统计算流体力学只能研究截面形状与气动参数之间关系的局限，可以将桥梁结构的动力特性以及结构阻尼纳入到计算之中，从而能够更加全面充分地研究桥梁的流固耦合现象。

目前流固耦合问题的数值计算方法有弱耦合法和强耦合法^{[9, 10]}，以往的大部分流固耦合方法是首先分别对流体方程和结构方程单独求解，再交换数据进行循环求解的弱耦合方法^[9]，这种方法的缺点是流体域网格和结构域网格不一致，使流固交界面上的平衡条件不能精确满足，导致结果易发散^[10]。本文拟采用同时求解流体方程和结构方程的强耦合方法解决这一问题，以增强结果的稳定性和准确性。

本文对考虑流固耦合作用的大跨度悬索桥风振响应进行研究，介绍了一种流固耦合分析的强耦合方法。该方法引入弹性模型处理流体变形，使其在边界上与结构网格保持一致。采用同时求解流体方程、结构方程和动网格弹性模型方程的变分方程的方法得到流固耦合体系的解，并给出了与同时求解方法相应的湍流模型和边界条件。最后将该方法用于大跨度悬索桥的风振响应分析，计算了桥梁横断面的气动力系数，并对桥梁进行了颤振分析。 2 控制方程和边界条件 2.1 空气流体控制方程

空气流体控制方程由为不可压黏性Navie-Stokes方程(简称N-S方程)，即连续方程和动量方程:

点击浏览公式

点击浏览公式

边界条件为：

点击浏览公式

点击浏览公式

点击浏览公式

大跨度悬索桥属于经历大变形的柔性结构，具有较强的几何非线性，因此其控制方程采用Total Lagrangian方程和大变形理论来描述：

在T.L.方法中，差分平衡方程基于初始未变形位形表示为:

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

流体和结构的耦合最重要的问题之一就是恰当解决流体和结构交界面处的数据信息，包括位移和牵引力等。这里引入弹性模型来处理流体变形，达到流体网格和结构网格在边界保持一致的目的，弹性模型的控制方程为：

点击浏览公式

流体和结构在交界面上的耦合条件为：

点击浏览公式

点击浏览公式

为增强解的稳定性，本文采用对流体方程、结构方程和网格中的弹性模型同时求解的方法求耦合体系的解。为了降低对解连续性的要求，这里采用Galerkin有限元法生成流体方程、结构方程和弹性模型方程的变分形式。

流体N-S方程的变分形式为:

点击浏览公式

点击浏览公式

结构方程的变分形式可以表示为:

点击浏览公式

弹性模型的变分形式为：

点击浏览公式

此时整理式(15)~式(18)，可以将流固耦合的强耦合方程写为如下形式：

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

图 1 强耦合方程的求解流程 Fig. 1Flowchart of solving strongly coupling equations

图选项

本文湍流的模拟采用两方程的SST k-ω模型。k-ω模型和k-ε模型间是通过引入调配函数来建立联系的。SST k-ω模型考虑了湍流的剪切应力并对逆压梯度下气流分离的发生和数量都具较高的预测精度。因此本文采用两方程的SST k-ω模型来模拟湍流：

点击浏览公式

点击浏览公式

式中，k是湍动能； x_j为坐标分量；w为耗散率；μ为液体黏性系数；u_i，u_j为流体的速度分量；P_k为平均应变率的张量模量。

模型中的系数是与下列模型相应系数的线性组合：

点击浏览公式

模型中的无量纲参数C_μ是平方摩擦速度与湍动能的比。为了反映紊流风场中的高频成分效应，C_μ取值为0.04，模型其他的相关参数也要随之改变。内层Wilcox k-ω模型的参数为α₁=0.413,β₁=0.033 3,σ_k1=1.176,σ_ω1=2,外层k-ω模型参数为α₂=0.20,β₂=0.036 8,σ_k2=1.0,σ_ω2=1.168^{[11, 12]}。

在强耦合整体方法中，湍流模型也采用Galerkin有限元法进行离散。湍流边界层采用滑移边界条件，垂直于壁面的流体速度由固体壁面的法线方向速度给定，并结合壁面摩擦边界条件：

点击浏览公式

点击浏览公式

点击浏览公式

点击浏览公式

这里采用上述方法对典型大跨度悬索桥—丹麦大带东桥主梁截面进行考虑流固耦合的风致响应计算，并与已由风洞试验结果进行结果对比。大带东桥有东西个引桥，跨度均为535 m，主跨跨径为1 624 m。其横断面尺寸如图 2所示，关于该桥的更详细介绍可参考文献^[14]。这里对图 2所示断面进行风振响应分析。结构密度为1.2 kg/m³，空气的运动黏性系数取为1.460 5×10^-5 m²/s,密度取为1.225 05 kg/m³，入口风速为40 m/s。雷诺数以断面宽度计算都为Re=10⁵，入口条件根据风速和攻角确定，攻角在-10~10°之间，时间步长为1.3×10^-3。这里采用FORTRAN编写网格划分和方程求解的结构化程序代码，电脑硬件条件为：Intel i5 3450S；CPU频率2800 MHz；内存容量：4 GB DDR3，1600 MHz。

图 2 大带东桥横断面尺寸(单位：mm) Fig. 2Cross-section dimensions of Dadai East Bridge(unit:mm)

图选项

计算区域外边界为矩形，计算域尺寸为279 m×93m(9B×3B,B=31 m,B为断面宽度)，上游距断面为1.5B(B为断面宽度)，两侧距断面为1B,下游距断面为6B。网格划分如图 3所示。计算网格采用分块结构化网格,在靠近物面处加密，单元数约为14万，离结构表面最近的网格间距小于0.003B，网格划分如图 3所示。

图 3 桥梁断面网格划分 Fig. 3Meshing of bridge cross-section

图选项

物体振动的初始条件为：先固定物体计算绕流直到稳定的旋涡脱落状态，然后让物体做强迫振动(即纯竖向或纯扭转运动)，计算振动物体和流场的相互作用。绕流计算初始条件为：在整个计算域内和边界上压力为0，在物面上沿水平x轴方向和竖向y轴方向均取为0，在流场内和其他边界上沿水平x轴方向取为入流速度，竖向y轴方向取为0。模型在平衡位置处的来流攻角为0°，纯竖弯运动的振幅是0.025B，纯扭转运动的幅值角是3°。

图 4给出了0°攻角时大桥断面周围的压力场分布情况。

图 4 0°攻角时大桥断面周围的压力场分布情况(单位：Pa) Fig. 4Pressure field distribution around bridge cross-section at 0° attack angle (unit:Pa)

图选项

接着采用本文方法对该桥梁进行了颤振分析，以确定临界颤振风速。此处攻角保持为0°，稳态分析的入口风速V₀分别取16.86，33.73，50.59，67.44，84.32 m/s，颤振分析的入口风速将前述风速分别采用如下公式折算：V_r=V₀/fB，式中B是桥面的特征长度=31 m，f是桥梁的自振频率。计算过程中雷诺数保持在10⁵，其他流体参数均与前述相同。这里只考虑桥梁的扭转运动，因为有研究表明^[14]耦合的竖向位移和转动对颤振临界风速影响不大。

这里采用Scanlan^[16]提出的颤振导数方法进行颤振分析，进而计算得到颤振临界风速。本文给出了采用本文方法计算得到的气动导数(具体公式可参加文献^[16])与风洞试验和已有数值模拟结果进行了比较，结果见图 5。

图 5 桥梁断面的气动导数结果比较 Fig. 5Comparison of aerodynamic derivatives of bridge cross-section

图选项

用最小二乘法确定Scanlan提出的气动力的表达式中的气动导数，并进一步用半逆解法确定颤振临界风速。

断面的气动导数如图 5所示。可以看出，本文方法得到的气动导数与Reinhold等的风洞试验结果^[17]吻合度较高，且与Walther的数值计算结果^[17]相比，更接近试验结果。本文方法更接近于试验结果的原因在于数值模拟计算气动导数的关键在于计算运动物体所承受的气动力，难点在于流体运动和结构运动的描述方法(欧拉法和拉格朗日法)的不同，在直接求解流体和结构的耦合运动方程时,在每一计算时步上要移动靠近结构表面的计算网格，而本文方法采用弹性模型处理流体网格变形，通过引入类似结构方程组来处理流体域连续变形的方法。它将流体网格看成是弹性模型，因而将网格移动问题转化为了线性固体力学问题，此时将实际的结构网格位移看成是边界条件，这种方法避开了通常方法更新结构网格的繁琐过程，也就提高了计算精度。

为进一步验证计算结果的可靠程度，采用本文计算得到的气动导数计算各断面的颤振临界风速，并与已有不同文献的试验^[17]和数值计算结果^[17]以及ANSYS结果做了比较，如表 1所示，其中ANSYS中采用MATRIX27单元来模拟风-桥梁形成的流固耦合体系，MATRIX27单元的刚度矩阵或阻尼矩阵是风速和振动频率的导数，整个过程就是对风速搜索和对频率迭代的过程。

从表 1中可以看出，采用本文方法计算得到的颤振临界风速与风洞试验值非常接近，而小于与其他数值计算方法得到的结果，可能的原因当结构发生振颤时，是结构空气动力失稳的表现，而此时描述风与结构的耦合效应的气动阻尼应为负值，因此考虑耦合效应的颤振临界风速要小于不考虑耦合效应的数值方法计算得到的结果。

为了进一步分析流固耦合作用对桥梁断面风振响应的影响，这里也研究了未发生颤振(稳定状态)下和发生颤振(失稳状态下)桥梁断面的风振响应(角位移)时程，对采用本文方法考虑流固耦合作用和不考虑流固耦合作用(不考虑流固耦合作用时采用ANSYS计算)的计算结果进行了比较，如图 6所示。

图 6 桥梁在不同风速下的风振响应时程 Fig. 6Time histories of wind-indcued vibration responses at different wind speeds

图选项

分析图 6，可以得到如下结论：

(1)从图 6(a)中可以看出，未发生颤振(稳定状态)时，不论是否考虑耦合作用，结构的响应都逐渐减小，且考虑耦合作用的风振响应要小于不考虑耦合作用的风振响应，说明在气动弹性稳定的前提下，耦合作用对桥梁结构起抑制作用，可以引起能量的耗散。根据文献^[18]的研究，在气动弹性稳定的前提下，描述风与结构的耦合效应的气动阻尼一般为正值，因此考虑耦合效应的风振响应等值小于不考虑耦合效应的相应值。

(2)从图 6(b)中可以看出，发生颤振(失稳状态)时，不论是否考虑耦合作用，随着风速增大，结构的响应都呈放大的趋势，说明了颤振现象的发生。且考虑耦合作用的风振响应要大于不考虑耦合作用的风振响应，说明在气动弹失稳的情况下，流固耦合作用进一步加深了结构的不稳定性，这与前述对颤振临界风速计算得到的结论一致，因此流固耦合效应在抗风设计时是必须考虑的重要因素。

图 7(a)给出了结构的最后稳定风振响应(颤振发生前，其风速为34.5 m/s)的自谱图，图 7(b) 给出了颤振刚发生时(风速为35.1 m/s)风振响应的自谱图。

分析图 7可以看出:

(1)考虑耦合效应时桥梁结构响应的大多数能量集中在低频段，说明低阶扭转频率与悬索桥的颤振临界风速有很大关系，扭转振型在悬索桥的颤振中占主要成分。事实上，弯扭频率在桥梁抗风实验中也是十分重要的衡量桥梁动力性能的指标，这也反映了本文方法计算的正确性。

(2)悬索桥的弯扭频率较低，这是由于悬索桥主梁的跨径较长，与斜拉桥相比,扭转振型出现较晚的原因。

另外作者还计算了更精细网格划分时的情况，结果表明当网格精度提高约15%时，计算精度仅提高约2%，稳定性基本不受影响，而耗费机时却增加了约30%。因此计算结果表明网格划分的精细程度对于计算结果的精度和稳定性影响是不大的。

图 7 不同风速下桥梁风振响应自谱图 Fig. 7Spectra of bridge wind-indcued vibration response at different wind speeds

图选项

本文对考虑流固耦合作用的大跨度悬索桥风振响应进行研究，介绍了一种流固耦合分析的强耦合方法，即在每一时间步内对流体控制方程和结构控制方程同时联立求解，计算出全场变量的值。采用该方法对大跨度悬索桥进行了气动力分析和颤振分析，得到的结果证明了本文的方法是成功的，主要结论包括：

(1) 考虑流固耦合效应时的阻力系数、升力系数和扭矩系数均比不考虑耦合效应时的相对误差要小，计算精度平均可达到50%，说明耦合效应对于大跨度悬索桥计算有重要影响，也证明了本文方法的正确性。

(2) 未发生颤振(稳定状态)时，不论是否考虑耦合作用，结构的响应都逐渐减小，且考虑耦合作用的风振响应要大于不考虑耦合作用的风振响应，说明在气动弹性稳定的前提下，耦合作用对索膜结构起抑制作用，可以引起能量的耗散。

(3) 发生颤振(失稳状态)时，不论是否考虑耦合作用，随着风速增大，结构的响应都呈放大的趋势，说明了颤振现象的发生。且考虑耦合作用的风振响应要大于不考虑耦合作用的风振响应，说明在气动弹失稳的情况下，流固耦合作用进一步加深了结构的不稳定性。

(4) 考虑耦合效应的颤振临界风速要小于不考虑耦合效应的数值方法计算得到的临界风速，可能的原因当结构发生振颤时，是结构空气动力失稳的表现，而此时描述风与结构的耦合效应的气动阻尼应为负值，因此流固耦合效应在桥梁抗风设计时是必须考虑的重要因素。