基于黏性涡域法的二维钝体绕流模拟

2023-11-02赵国辉王治超郝键铭高广中

空气动力学学报 2023年9期

赵国辉，王治超，郝键铭，高广中

(1.长安大学公路学院，西安 710064；2.武汉市政工程设计研究院有限责任公司，武汉 430023)

0 引言

涡方法是计算流体力学中的一类无网格拉格朗日方法，该方法可以将计算资源集中到占比较小的涡流区域，因此可以提高计算效率。由于计算点的运动与流体颗粒的运动接近或重合，用涡来描述流场，使得涡域具有很好的可视化效果，便于分析旋涡的形成机理。涡方法不仅能模拟流体中的大尺度拟序结构，也能模拟小尺度涡结构，且在模拟高雷诺数时不需要湍流模型。此外，由于没有计算网格，涡方法可以被应用于更广泛的几何形状。

最初的涡方法研究未考虑黏性效应，仅适用于理想流体。但在钝体周围，不仅分离点的位置未知并且还会随时间变化，这种情况下黏性起主要作用。Chorin[1]将涡量-速度方程分解为对流项和黏性扩散项两部分，而处理黏性扩散的方法有随机涡方法和确定性涡方法。常见的随机涡方法有随机走步法与重采样方法。随机涡方法在计算过程中产生随机数，虽然原理和实现方法均较简单，但由于扩散的随机性，此类方法在求解涡量时存在较大的统计误差且收敛速度较低。常见的确定性涡方法有粒子强度交换法（particle strength exchange method,PSE）[2]、涡量重分布法[3]和扩散速度法（diffusion velocity method,DVM）[4]。粒子强度交换法中，粒子位置的拉普拉斯算子由周围粒子的积分算子代替，其精确度强烈依赖于用于离散积分的求积规则。之后进一步演化了重新网格化的PSE 方法[5-6]，该方法不使用积分算子因而得以保留无网格化，但需要在每个时间步内求解N个粒子的N个欠定系统，因而随着粒子数量增加，计算效率显著降低。这种新PSE 方法对重新划分网格的严重依赖促使了涡量重分布法的发展。扩散速度法[4]提出了扩散速度的概念，在该方法中，涡流粒子保持其环流并根据速度场运动，速度场是对流速度和扩散速度的叠加，主要计算变量为涡量。Clarke 等[7]发现，扩散速度法在距离物体较近时效果很好，但在距离物体较远时就变成了随机走步法。且由于旋涡排斥的非单调特性，即当旋涡彼此接近时排斥率趋于零，导致旋涡“黏附”，从而不能准确模拟（尤其是物体表面线附近的）涡量演变。Dynnikova 等[8]认为，流域中的涡量演变可以看作是它沿总速度场流线的转移，从而发展了扩散速度法，使得涡量演变更为准确，并且很好地描述了边界层，确保了在物体表面附近的扩散速度的正确表达，并将这种方法命名为黏性涡域法（viscous vortex domains method,VVD）。

本文拟从黏性涡域法的发展、原理、实例等各方面展开阐述，分步介绍黏性涡域法的实现以及各步的理论公式；利用该方法模拟不同雷诺数下圆柱绕流、不同无量纲频率和无量纲振幅的振动圆柱气动力时程、不同风向角下的方柱绕流，并将计算结果与现有文献的结果进行对比，验证黏性涡域法的有效性和计算精度。

1 黏性涡域法

黏性涡域法模拟过程为：

1）生成涡量。该方法主要计算变量为涡量，涡量用物体表面线上的涡片代替，涡片强度通过无滑移边界条件用边界积分方程的形式写出。

2）涡量演变。将涡量从物体表面线转移到流场中，涡流尾流演变模拟的同时得到流体速度。

3）气动计算。通过Uhlman 提出的积分方程求解压力分布，将物体表面线各点压力与摩擦力的竖向及水平分量累积求和并无因次化，可得升力系数与阻力系数。

1.1 生成涡量

物体表面的涡量可用强度待定的自由涡片（沿物体表面线的无限薄的涡量分布）表示，并将其转移至流场中，这些涡片作为流场涡源并根据控制方程发展运动。这就意味着模拟物体表面线影响的涡片应为自由涡片而非附加涡片，为此Lighthill[9]提出，由固体边界的边界层效应产生了沿固体边界连续分布的涡量，该连续涡流近似用离散的涡片模拟，然后集中在一定强度涡片中的涡量成为涡流尾迹的一部分。

涡片强度用边界积分的形式表示为[10]：

向量核Q(r,ξ)类似于二维格林函数的梯度：

其中：K为物体表面线；S为流场；γ(ξ)为涡片强度；ξ为物体表面线上的点；γ(r)为附加涡片强度，其等于物体表面速度的切线分量。

为方便表示，定义式（1）等号右边为fτ(r)，且fτ(r)=f(r)·τ(r)。n(r)和τ(r)分别是外法线向量和切线向量，且n(r) ×τ(r)=e。由此可得到取决于来流流速、物体表面线速度和流动中的涡量分布的已知函数f(r)：

其中：V∞为来流流速；N为流场内涡流粒子数；ri和Гi，i=(1,2,···,N)，分别为涡流粒子的位置和环量。

为得到边界积分方程（1）的唯一解，需要引入如下附加条件：

利用Lifanov[11]的方法，将方程（5）与边界积分方程（1）联合求解。求解边界积分方程最简单、最普遍的方法是伽辽金法。但传统的伽辽金法无法重建边界积分方程的解。解的变化越明显，涡粒子越接近物体表面线，显著的局部误差会导致近壁区域的速度场重建不正确。为解决该问题，Soldatova 等[12]提出了在数值解中加入一项，该项考虑了由于密集放置的涡粒子的影响而引起的解在一个面元上的局部行为：

其中：γ*为修正项；ϕq为基函数（当q=0 时为常数，q=1时为线性，q=2 时为二次）；γq为未知系数，可以通过伽辽金方法确定。为简便起见，考虑正交基函数，当ξ在对应面元上时，ϕ0(ξ)≡1 并且ϕ1和ϕ2都等于1，在面元之外与其相关的所有基函数ϕq都等于零。

通过节点划分对物体表面线进行离散，将物体表面线划分为多个节点，两节点间的节段称为面元，如图1 所示。选择离涡流粒子相当近的面元，设面元k∈[kb,ke]，这些面元被曲率为κk的密切圆圆弧所代替，利用共形映射理论为每个相当接近涡流粒子的面元写下精确解：

图1 物体表面节点划分和面元Fig.1 Node division of an airfoil and panels

其中：Γw为位于rw点的涡流粒子环量；rk(ξ)为密切圆圆弧上一点；nk(ξ)为点rk(ξ)处的外单位法向量；为密切圆的圆心；；Rk为密切圆半径。

位于点rw涡流粒子的环量可进一步表示为：

其中：Гw,k为分布在第k面元上的总环量，用来考虑涡粒子的贡献；dw,k为从涡流粒子到第k面元的距离；s为弧长参数值，，sw对应于涡流粒子位置在第k个面元上的投影。各物理量如图2 所示。

图2 相邻两个曲线面元和靠近的涡流粒子Fig.2 Two adjacent curved panels andapproaching vortex elements

1.2 涡量演变

涡量生成后，需要将涡量从物体表面线转移到流域。在物体表面线处形成涡片的分布涡量被转化为单独的涡流粒子，成为旋涡尾流的一部分。不考虑流动可压缩性，就涡量而言，Navier-Stokes 方程可表示为：

其中：V为对流速度；W为由黏性效应引起的扩散速度；涡量场Ω只有一个非零分量，与流动平面正交，可写为Ω=Ωe。方程（9）是指流动中存在的涡量随速度移动，新的涡量只在流动区域的边界上（即物体的边界上）产生。

涡方法模拟无黏、不可压缩流的涡流尾流演化时，需要对以下常微分方程组进行求解：

其中，ri是第i个涡流粒子的位置。由于涡量只是随着流速在流域中传递，对流速度V用Biot-Savart 定律由涡量分布计算得到：

对于扩散速度，有：

式中 υ为运动黏度。式（12）表明扩散速度与运动黏度υ成正比，并且取决于点r附近流动区域中的涡量分布和流动区域的边界形状。

1.3 气动计算

图3 为圆柱绕流示意图，图中Fp为压力、τ为物体表面的切向应力。物体表面受到流体的压力和剪切应力，由Uhlman 提出的积分方程，可以求解流场压力分布，从而计算压力：

图3 圆柱绕流示意图Fig.3 Flow around a cylinderer

其中：β为常数，二维问题的计算域内部β=2p，计算域边界上β=p；B为停滞热焓；p为压强；Z为计算域内部；Y为计算域边界；G为标量格林函数。

设在物面处有Nw个物面涡元，且第j个物面涡元的环量为Γj，在计算域内部Z有Nz个自由涡，且第k个自由涡的环量为Γk，则物面i点处的压力积分方程的离散形式为：

物体静止时，有：

物体表面涡和流场中自由涡的涡强与环量关系为：

将式（17～19）代入式（16）可得：

式（20）中，等号右端第一项表示物体表面涡元的存在对物面压强的贡献，第二项代表流场中自由涡的运动对物面压强的贡献。求解式（20）得到Bi，再将其代入式（14），可得到物面上i点的压强：

对于物体表面的剪切应力τw，假定边界层内速度是线性分布的，则物体所受的剪切应力可通过下式求解：

综上，物体所受的流体力为：

式中，en和eτ分别代表物体表面法向量和切向量。将物体所受流体力沿顺流向和横流向分解，可得到物体所受的阻力FD和升力FL，无因次化后可得阻力系数CD和升力系数CL：

2 数值模拟

基于黏性涡域法，在Linux 系统下利用Lua 语言编制圆柱绕流的计算程序，利用Gnuplot 实现流场可视化。进行了低雷诺数和高雷诺数下静止和振动圆柱的模拟。

2.1 静止圆柱绕流模拟

分别取雷诺数Re=2 × 102、1 × 103、1 × 105进行模拟计算和验证，根据多次试算综合考虑计算精度与效率，面元长度取0.01 m。图4～图6 是基于黏性涡域法模拟的静止圆柱涡流演变过程、压力云图及气动力时程。用点模拟涡流尾流的离散涡涡流粒子位置，分别用红色、蓝色点表示正、负循环涡流。t*=Ut/D为无量纲时间。阻力系数CD=FD/(0.5ρDU2)，升力系数CL=FL/(0.5ρDU2)，ρ为流体密度，FD和FL分别为柱体所受的阻力和升力。

图4 静止圆柱涡量生成及演变Fig.4 Generation and evolution of vorticity around a static cylinderc cylinder

图5 静止圆柱压力云图Fig.5 Static cylinder pressure contour

图6 不同雷诺数静止圆柱气动力时程Fig.6 Aerodynamic time histories of a static cylinder with different Reynolds numbersbers

由图4 可见，VVD 能够模拟涡流演变过程，且利用该方法模拟的卡门涡街与周志勇等[13]利用离散涡方法、王亚玲等[14]利用计算流体软件CFX-4、桑文慧等[15]利用SIMPLEC 算法和张伟伟等[16]的模拟结果基本一致。时间步长为0.05、计算步为20 000 时，VVD 方法的计算时间约为10 min。

由图6 可知，在不同雷诺数下，本文方法与文献[21,23]的气动力时程计算结果基本一致。为进一步验证本文方法的准确性，将气动计算结果与现有文献结果进行对比，结果见表1～表3，表中数据均为旋涡脱落稳定后的数据。

表1 不同方法气动力计算结果对比（Re=2 × 102）Table 1 Comparison of aerodynamics obtained by different methods (Re=2 × 102)

表2 不同方法气动力计算结果对比（Re=1 × 103）Table 2 Comparison of aerodynamics obtained by different methods (Re=1 × 103)

表3 不同方法气动力计算结果对比（Re=1 × 105）Table 3 Comparison of aerodynamics obtained by different methods (Re=1 × 105)

通过气动力时程曲线以及与各方法气动力系数幅值对比，可得：

1）Re=2 × 102时，本文方法与现有文献中各方法计算结果基本一致，其中与Fluent 计算结果最为吻合，阻力系数比离散涡方法结果稍微偏大，斯特劳哈尔数与其他方法相比偏小，但与文献[28]结果接近；

2）Re=1 × 103时，与Fluent 计算结果较为吻合，但比其他方法（尤其比试验）所得结果偏大，斯特劳哈尔数基本一致；

3）Re=1 × 105时，升力和阻力系数计算结果与其他文献吻合较好，斯特劳哈尔数基本一致。

综上，雷诺数低于1 × 105时，黏性涡域法能够较为准确地模拟静止圆柱绕流流场和计算升力/阻力系数。

2.2 振动圆柱绕流模拟

模拟圆柱沿垂直来流方向做受迫振动，运动方程为y=Asin(2πfet)，fe为圆柱的振动频率，t为圆柱受迫运动的运动时刻。无量纲振动频率f*=fe/fs分别取0.5、0.7、0.75、0.95、1.0、1.2，fs为静止圆柱自然脱落频率。无量纲振幅A*=A/D分别取0.2、0.4，A为振荡幅值。振动圆柱涡流演变过程、压力云图及气动力时程见图7～图10。图7 中，尾流呈现出经典的卡门涡街。在圆柱单个振荡周期内，涡从圆柱两侧交替脱落，即2S 模式，且涡量场轮廓彼此十分相似。

图7 振动圆柱涡量生成及演变Fig.7 Generation and evolution of vorticity around an oscillating cylindering cylinder

图8 振动圆柱压力云图Fig.8 Pressure contours around an oscillating cylinderng cylinder

将本文方法与其他方法[29-33]的计算结果进行对比（图9～图10），可见，在不同频率和振幅下，虽然气动力时程曲线略有不同，但大体趋势基本保持一致。为更加直观地进行对比，图11 给出了阻力系数平均值CDM的对比情况。

图9 振动圆柱气动力时程（A＊=0.2）Fig.9 Aerodynamic timeoscilhistories of an oscillating cylinder（A＊=0.2）

图10 振动圆柱气动力时程（A＊=0.4）Fig.10 Oscillating cylinder aerodynamic time histohistories of an oscillating cylinder（A＊=0.4）er（A＊=0.4）

图11 振动圆柱结果对比Fig.11 Oscillating cylinder results comparison

2.3 不同风向角方柱绕流模拟

目前大跨桥梁塔柱均以方形或在其基础上进行切角或凹槽处理的截面为主。与圆柱截面不同，对于方柱这种边缘锋利的截面，分离点可能固定在拐角处，边界层在其前方两拐角点分离，对雷诺数的依赖性也小得多，阻力系数也在一定雷诺数范围内基本不变。然而不同长宽比的矩形截面以及桥塔受到不同的风向角时都会产生不同的流场，因此有必要验证方柱绕流流场及气动系数变化规律。实际应用中，桥塔断面长宽比接近于1，且研究方柱的学者较多，因此本节以方柱为例进行绕流模拟。

利用黏性涡域法模拟计算0°～45°风向角下方柱绕流流场及气动系数，方柱边长取L=1 m，雷诺数为1 000，时间步长0.01 s。图12 给出了0°、15°和45°风向角下的方柱涡量图。从图中可以清楚观察到，旋涡上下交替脱落形成卡门涡街。图13 为风向角为0°、15°和45°时的升/阻力系数时程曲线。图14 为不同方法、各风向角下方柱升/阻力系数结果对比。图15 为风向角15°时剪切层示意图。可以看出：本文方法的计算结果及趋势与其他方法得到的基本一致；对于阻力系数，杨素哲[34]和Naudascher 等[35]计算的最小值出现在风向角为10°时，Norberg[36]试验的最小值出现在风向角为12.5°时，本文和Taylor 等[37]计算的最小值出现在风向角为15°时。对于升力系数，杨素哲[34]计算的最小值出现在风向角为10°时，Norberg[36]试验的最小值出现在风向角为12.5°时，Naudascher 等[35]计算的最小值出现在风向角13°左右，本文和Taylor等[37]计算的最小值出现在风向角为15°时。综上，可以得出结论，本文方法能够有效地模拟方柱绕流。

图12 不同风向角下方柱绕流涡量Fig.12 Vorticity around a rectangular cylinder under different wind direction angles

图13 不同风向角下方柱绕流升、阻力时程Fig.13 Time histories of lift and drag coefficients around a rectangular cylinder under different wind direction angles

图14 不同方法各风向角方柱升、阻力系数对比Fig.14 Comparison of rectangular cylinder drag and lift coefficients at different wind direction angles by on angles by different methods

图15 风向角15°时剪切层Fig.15 Shear layer when the wind direction angle is 15°

图14 还显示，升力系数和阻力系数均随风向角的增大先减小后增大，尝试分析原因如下。当风向角较小时，随着风向角增大，平均阻力系数减小，主要是从B 点脱落的剪切层间接接触C 点；当风向角约为15°时，则该剪切层完全附着，最初与B 分离的剪切层现在与C 分离，因而产生更窄的尾迹和更低的平均阻力系数；当风向角继续增大时，尾流宽度又会重新变宽，从而平均阻力系数又会变大。虽然不同方法和试验得到的最小平均阻力系数风向角不相同（相差在5°之内），但黏性涡域法很好地再现了这种趋势。升力系数随风向角变化也有相同的趋势。分离的涡泡产生比AD 面更大的局部吸力，因而产生负升力系数；当风向角为15°时剪切层完全附着，此时负升力系数最大；随着风向角进一步增大，BC 面上的局部吸力减小，升力系数逐渐增加。