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

文章编号：0258-1825(2023)09-0107-12基于黏性涡域法的二维钝体绕流模拟赵国辉1，*，王治超1，2，郝键铭1，高广中1(1.长安大学公路学院，西安710064；2.武汉市政工程设计研究院有限责任公司，武汉430023)摘要：针对传统计算流体力学方法需要划分网格且计算耗时长等问题，一种确定性涡方法（黏性涡域法）逐渐发展并日臻成熟。本文利用该方法计算模拟了二维钝体的流场、气动系数及斯特劳哈尔数，并与其他成熟 CFD 软件计算结果对比来验证模拟精度。黏性涡域法的模拟过程为生成涡量、涡量演变及气动计算。文章首先推导了该方法各过程的基本方程，然后基于 Lua 语言编制了圆柱和方柱绕流计算程序，最后利用 Gnuplot 实现流场可视化。利用黏性涡域法模拟不同雷诺数下圆柱绕流、不同无量纲频率和无量纲振幅的振动圆柱气动力时程、不同风向角下的方柱绕流，结果表明：雷诺数为 21021105时，基于黏性涡域法的圆柱绕流气动计算结果与其他方法的计算结果基本一致；基于该方法的不同风向角下方柱绕流气动计算结果与其他文献结果趋势相同。关键词：黏性涡域法；涡方法；气动力；圆柱绕流；雷诺数中图分类号：O357.5文献标识码：Adoi：10.7638/kqdlxxb-2022.0084Simulation of flow around two-dimensional blunt bodybased on the viscous vortex domain methodZHAOGuohui1，*，WANGZhichao1，2，HAOJianming1，GAOGuangzhong1(1.School of Highway,Changan University,Xian710064,China；2.Wuhan Municipal Engineering Design&Research Institute Co.,Ltd,Wuhan430023,China)Abstract:ThemeshgenerationisindispensablefortraditionalCFDmethods,butitusuallytakesalongtime,especiallyforcomplexgeometries.Bycomparison,themesh-freeviscousvortexdomain(VVD)methodisgettingincreasinglymoreattentionforflowfieldsimulationandhasundergonesignificantdevelopmentrecently.ThispaperfirstverifiestheaccuracyofVVDbysimulatingtheflowfieldaroundabluffbodysectionandcomparingtheresultswiththoseofotherCFDmethods.Then,thegoverningequationsforthethreestepsofVVD,vorticitygeneration,vorticityevolution,andaerodynamiccalculation,arederived,basedonwhichaflowfieldcalculationprogramisdevelopedusingLualanguage.Finally,useGnuplottovisualizetheflowfield.VVDis used to simulate flow fields around a circular cylinder at different Reynolds numbers,the aerodynamichistoriesofavibratingcircularcylinderatdifferentdimensionlessfrequenciesandamplitudes,andtheflowfieldsaroundarectangularcylinderatdifferentwindangles.ResultsshowthatwhentheReynoldsnumberisbetween2102and1105,theflowfieldsaroundacircularcylinderobtainedontheVVDmethodareconsistentwiththoseofothermethods.Moreover,theaerodynamicsofarectangularcylinderwithdifferentwindanglesbasedonthismethodhavethesametrendasthosereportedinpreviousstudies.Keywords:viscous vortex domain method；vortex method；aerodynamic force；flow around cylinder；Reynoldsnumber收稿日期：2022-05-29；修订日期：2022-07-26；录用日期：2022-08-18；网络出版时间：2022-09-17基金项目：陕西省自然科学基础研究计划（2022JM-198）；中央高校基本科研业务费专项资金（300102210208）作者简介：赵国辉*（1978），男，河北邯郸人，副教授，研究方向：桥梁抗风抗震及振动控制.E-mail：引用格式:赵国辉,王治超,郝键铭,等.基于黏性涡域法的二维钝体绕流模拟J.空气动力学学报,2023,41(9):107118.ZHAOGH,WANGZC,HAOJM,etal.Simulationofflowaroundtwo-dimensionalbluntbodybasedontheviscousvortexdomainmethodJ.ActaAerodynamicaSinica,2023,41(9):107118(inChinese).doi:10.7638/kqdlxxb-2022.0084第41卷第9期空气动力学学报Vol.41，No.92023年9月ACTA AERODYNAMICA SINICASept.，2023 0 引言涡方法是计算流体力学中的一类无网格拉格朗日方法，该方法可以将计算资源集中到占比较小的涡流区域，因此可以提高计算效率。由于计算点的运动与流体颗粒的运动接近或重合，用涡来描述流场，使得涡域具有很好的可视化效果，便于分析旋涡的形成机理。涡方法不仅能模拟流体中的大尺度拟序结构，也能模拟小尺度涡结构，且在模拟高雷诺数时不需要湍流模型。此外，由于没有计算网格，涡方法可以被应用于更广泛的几何形状。最初的涡方法研究未考虑黏性效应，仅适用于理想流体。但在钝体周围，不仅分离点的位置未知并且还会随时间变化，这种情况下黏性起主要作用。Chorin1将涡量-速度方程分解为对流项和黏性扩散项两部分，而处理黏性扩散的方法有随机涡方法和确定性涡方法。常见的随机涡方法有随机走步法与重采样方法。随机涡方法在计算过程中产生随机数，虽然原理和实现方法均较简单，但由于扩散的随机性，此类方法在求解涡量时存在较大的统计误差且收敛速度较低。常见的确定性涡方法有粒子强度交换法（particlestrengthexchangemethod,PSE）2、涡量重分布 法3和 扩 散 速 度 法（diffusion velocity method,DVM）4。粒子强度交换法中，粒子位置的拉普拉斯算子由周围粒子的积分算子代替，其精确度强烈依赖于用于离散积分的求积规则。之后进一步演化了重新网格化的 PSE 方法5-6，该方法不使用积分算子因而得以保留无网格化，但需要在每个时间步内求解N 个粒子的 N 个欠定系统，因而随着粒子数量增加，计算效率显著降低。这种新 PSE 方法对重新划分网格的严重依赖促使了涡量重分布法的发展。扩散速度法4提出了扩散速度的概念，在该方法中，涡流粒子保持其环流并根据速度场运动，速度场是对流速度和扩散速度的叠加，主要计算变量为涡量。Clarke 等7发现，扩散速度法在距离物体较近时效果很好，但在距离物体较远时就变成了随机走步法。且由于旋涡排斥的非单调特性，即当旋涡彼此接近时排斥率趋于零，导致旋涡“黏附”，从而不能准确模拟（尤其是物体表面线附近的）涡量演变。Dynnikova 等8认为，流域中的涡量演变可以看作是它沿总速度场流线的转移，从而发展了扩散速度法，使得涡量演变更为准确，并且很好地描述了边界层，确保了在物体表面附近的扩散速度的正确表达，并将这种方法命名为黏性涡域法（viscousvortexdomainsmethod,VVD）。本文拟从黏性涡域法的发展、原理、实例等各方面展开阐述，分步介绍黏性涡域法的实现以及各步的理论公式；利用该方法模拟不同雷诺数下圆柱绕流、不同无量纲频率和无量纲振幅的振动圆柱气动力时程、不同风向角下的方柱绕流，并将计算结果与现有文献的结果进行对比，验证黏性涡域法的有效性和计算精度。1 黏性涡域法黏性涡域法模拟过程为：1）生成涡量。该方法主要计算变量为涡量，涡量用物体表面线上的涡片代替，涡片强度通过无滑移边界条件用边界积分方程的形式写出。2）涡量演变。将涡量从物体表面线转移到流场中，涡流尾流演变模拟的同时得到流体速度。3）气动计算。通过 Uhlman 提出的积分方程求解压力分布，将物体表面线各点压力与摩擦力的竖向及水平分量累积求和并无因次化，可得升力系数与阻力系数。1.1 生成涡量物体表面的涡量可用强度待定的自由涡片（沿物体表面线的无限薄的涡量分布）表示，并将其转移至流场中，这些涡片作为流场涡源并根据控制方程发展运动。这就意味着模拟物体表面线影响的涡片应为自由涡片而非附加涡片，为此 Lighthill9提出，由固体边界的边界层效应产生了沿固体边界连续分布的涡量，该连续涡流近似用离散的涡片模拟，然后集中在一定强度涡片中的涡量成为涡流尾迹的一部分。涡片强度用边界积分的形式表示为10：zKQ(r,)()dl12(r)=wS(r)n(r)2|r|2()dSV(r)(r)(1)Q(r,)=Q(r,)n(r)(2)向量核 Q(r,)类似于二维格林函数的梯度：Q(r,)=e(r)2|r|2(3)其中：K 为物体表面线；S 为流场；()为涡片强度；为物体表面线上的点；(r)为附加涡片强度，其等于物体表面速度的切线分量。108空气动力学学报第41卷为方便表示，定义式（1）等号右边为 f(r)，且f(r)=f(r)(r)。n(r)和(r)分别是外法线向量和切线向量，且 n(r)(r)=e。由此可得到取决于来流流速、物体表面线速度和流动中的涡量分布的已知函数 f(r)：f(r)=V+Ni=1Q(r,ri)i(4)其中：V为来流流速；N 为流场内涡流粒子数；ri和i，i=(1,2,N)，分别为涡流粒子的位置和环量。为得到边界积分方程（1）的唯一解，需要引入如下附加条件：zK()dl=(5)利用 Lifanov11的方法，将方程（5）与边界积分方程（1）联合求解。求解边界积分方程最简单、最普遍的方法是伽辽金法。但传统的伽辽金法无法重建边界积分方程的解。解的变化越明显，涡粒子越接近物体表面线，显著的局部误差会导致近壁区域的速度场重建不正确。为解决该问题，Soldatova 等12提出了在数值解中加入一项，该项考虑了由于密集放置的涡粒子的影响而引起的解在一个面元上的局部行为：()=kk()0()+qqq()(6)其中：*为修正项；q为基函数（当 q=0 时为常数，q=1时为线性，q=2 时为二次）；q为未知系数，可以通过伽辽金方法确定。为简便起见，考虑正交基函数，当 在对应面元上时，0()1 并且 1和 2都等于 1，在面元之外与其相关的所有基函数 q都等于零。通过节点划分对物体表面线进行离散，将物体表面线划分为多个节点，两节点间的节段称为面元，如图 1 所示。选择离涡流粒子相当近的面元，设面元kkb,ke，这些面元被曲率为 k的密切圆圆弧所代替，利用共形映射理论为每个相当接近涡流粒子的面元写下精确解：k()=w2rk()rw|rk()rw|2rk()rmk?rk()rmk?2+rk()rck?rk()rck?2nk()(7)rckrmk=rck+(rwrck)R2k/?rwrck?2其中：w为位于 rw点的涡流粒子环量；rk()为密切圆圆弧上一点；nk()为点 rk()处的外单位法向量；为密切圆的圆心；Rk为密切圆半径。物体表面线物体表面线面元k面元k+1面元k1涡流粒子图 1 物体表面节点划分和面元Fig.1 Node division of an airfoil and panels位于点 rw涡流粒子的环量可进一步表示为：w=kek=kbw,kw,k=wKkk(s)dsw2kLk2(3)w,kdw,kk2(sin2(1)w,k+sin2(2)w,k)(8)(1)w,k=arctandw,k/(Lk2uw,k)(2)w,k=arctandw,k/(Lk2+uw,k)(3)w,k=(1)w,k(2)w,kuw,k=swsi12sk12=(sk1+sk)/2其中：w,k为分布在第 k 面元上的总环量，用来考虑涡粒子的贡献；dw,k为从涡流粒子到第 k 面元的距离；s 为弧长参数值，sw对应于涡流粒子位置在第 k 个面元上的投影。各物理量如图 2 所示。w面元k面元k1uw,kdw,kswsk12sk1sk2图 2 相邻两个曲线面元和靠近的涡流粒子Fig.2 Two adjacent curved panels andapproaching vortex elements 1.2 涡量演变涡量生成后，需要将涡量从物体表面线转移到流域。在物体表面线处形成涡片的分布涡量被转化为单独的涡流粒子，成为旋涡尾流的一部分。不考虑流动可压缩性，就涡量而言，Navier-Stokes 方程可表示为：t+(V+W)=0(9)其中：V 为对流速度；W 为由黏性效应引起的扩散速第9期赵国辉等：基于黏性涡域法的二维钝体绕流模拟109度；涡量场 只有一个非零分量，与流动平面正交，可写为=e。方程（9）是指流动中存在的涡量随速度移动，新的涡量只在流动区域的边界上（即物体的边界上）产生。涡方法模拟无黏、不可压缩流的涡流尾流演化时，需要对以下常微分方程组进行求解：dridt=V(ri)+W(ri),i=1,N(10)其中，ri是第 i 个涡流粒子的位置。由于涡量只是随着流速在流域中传递，对流速度 V 用 Biot-Savart 定律由涡量分布计算得到：V(ri)=V+ws(r)(r)2|r|2dS+zKk(r)2|r|2()dl=V+Ni=1k(rri)2|rri|2idS+zKk(r)2|r|2()dl(11)对于扩散速度，有：W(ri)=ln|(r)|=(r)(r)(12)式中 为运动黏度。式（12）表明扩散速度与运动黏度成正比，并且取决于点 r 附近流动区域中的涡量分布和流动区域的边界形状。1.3 气动计算图 3 为圆柱绕流示意图，图中 Fp为压力、为物体表面的切向应力。物体表面受到流体的压力和剪切应力，由 Uhlman 提出的积分方程，可以求解流场压力分布，从而计算压力：B+zYBGndl=zYnutG+G(n)dl+xZGudY(13)其中：为常数，二维问题的计算域内部=2p，计算域边界上=p；B 为停滞热焓；p 为压强；Z 为计算域内部；Y 为计算域边界；G 为标量格林函数。B=p+u22(14)G=ln(1r)(15)设在物面处有 Nw个物面涡元，且第 j 个物面涡元的环量为 j，在计算域内部 Z 有 Nz个自由涡，且第k 个自由涡的环量为 k，则物面 i 点处的压力积分方程的离散形式为：Bi+Nwj=1,j,iBiGijnlj=Nwj=1,j,injujtGijljNwj=1,j,iGij(njj)lj+Nzk=1Gik(ukk)Yk(16)物体静止时，有：ujt=0(17)物体表面涡和流场中自由涡的涡强与环量关系为：jlj=jnj(18)kYk=knk(19)将式（1719）代入式（16）可得：Bi+Nwj=1,j,iBiGijnlj=Nwj=1,j,iGij(njjnj)+Nzk=1Gik(ukknk)(20)式（20）中，等号右端第一项表示物体表面涡元的存在对物面压强的贡献，第二项代表流场中自由涡的运动对物面压强的贡献。求解式（20）得到 Bi，再将其代入式（14），可得到物面上 i 点的压强：pi=Biu2i2(21)对于物体表面的剪切应力 w，假定边界层内速度是线性分布的，则物体所受的剪切应力可通过下式求解：w=uy=(22)综上，物体所受的流体力为：F=zYpends+zYeds(23)yABCxDFp图 3 圆柱绕流示意图Fig.3 Flow around a cylinderer110空气动力学学报第41卷式中，en和 e分别代表物体表面法向量和切向量。将物体所受流体力沿顺流向和横流向分解，可得到物体所受的阻力 FD和升力 FL，无因次化后可得阻力系数CD和升力系数 CL：CD=2FDU2D,CL=2FLU2D(24)2 数值模拟基于黏性涡域法，在 Linux 系统下利用 Lua 语言编制圆柱绕流的计算程序，利用 Gnuplot 实现流场可视化。进行了低雷诺数和高雷诺数下静止和振动圆柱的模拟。2.1 静止圆柱绕流模拟分别取雷诺数 Re=2102、1103、1105进行模拟计算和验证，根据多次试算综合考虑计算精度与效率，面元长度取 0.01m。图 4图 6 是基于黏性涡域法模拟的静止圆柱涡流演变过程、压力云图及气动力时程。用点模拟涡流尾流的离散涡涡流粒子位置，分 别 用 红 色、蓝 色 点 表 示 正、负 循 环 涡 流。t*=Ut/D 为无量纲时间。阻力系数 CD=FD/(0.5DU2)，升力系数 CL=FL/(0.5DU2)，为流体密度，FD和 FL分别为柱体所受的阻力和升力。由图 4 可见，VVD 能够模拟涡流演变过程，且利用该方法模拟的卡门涡街与周志勇等13利用离散涡方法、王亚玲等14利用计算流体软件 CFX-4、桑文慧等15利用 SIMPLEC 算法和张伟伟等16的模拟结果基本一致。时间步长为 0.05、计算步为 20000 时，VVD 方法的计算时间约为 10min。44t*=044t*=50t*=8044t*=1t*=2t*=3t*=4t*=5t*=6(a)涡量生成(b)尾流未完全发展时涡量演变(c)尾流完全发展时涡量演变图 4 静止圆柱涡量生成及演变Fig.4 Generation and evolution of vorticity arounda static cylinderc cylinder(a)0时刻(b)尾流未完全发展(c)尾流完全发展1.51.0t*=01.51.0t*=1t*=2t*=3t*=4t*=5t*=61.51.0t*=50t*=80图 5 静止圆柱压力云图Fig.5 Static cylinder pressure contour第9期赵国辉等：基于黏性涡域法的二维钝体绕流模拟111由图 6 可知，在不同雷诺数下，本文方法与文献 21,23 的气动力时程计算结果基本一致。为进一步验证本文方法的准确性，将气动计算结果与现有文献结果进行对比，结果见表 1表 3，表中数据均为旋涡脱落稳定后的数据。表 1 不同方法气动力计算结果对比（Re=2 102）Table 1 Comparison of aerodynamics obtainedby different methods(Re=2 102)阻力系数CD升力系数CL斯特劳哈尔数St最大值 最小值最大值 最小值本文，VVD1.551.410.720.720.178丁代伟，CFX171.561.490.590.590.2Lecointe，Fluent181.501.420.700.700.23苏波，Fluent191.621.310.730.730.2周志勇，离散涡131.32（平均值）0.720.780.194Roshko，试验200.170.19Wille，试验211.30（平均值）表 2 不同方法气动力计算结果对比（Re=1 103）Table 2 Comparison of aerodynamics obtainedby different methods(Re=1 103)阻力系数CD升力系数CL斯特劳哈尔数St最大值最小值最大值最小值本文，VVD1.651.391.141.100.200苏波，Fluent191.701.331.311.33Jordon，无穷大法221.24（平均值）0.206Patel，试验230.980.005(平均值)0.21表 3 不同方法气动力计算结果对比（Re=1 105）Table 3 Comparison of aerodynamics obtainedby different methods(Re=1 105)阻力系数CD升力系数CL斯特劳哈尔数St最大值 最小值最大值 最小值本文，VVD1.251.070.410.400.200李燕玲，DES241.21（平均值）0.410.410.196余化军，CFX251.301.290.510.510.2Schewe261.20（平均值）Schlichting，试验271.25（平均值）0.19通过气动力时程曲线以及与各方法气动力系数幅值对比，可得：1）Re=2102时，本文方法与现有文献中各方法计算结果基本一致，其中与 Fluent 计算结果最为吻合，阻力系数比离散涡方法结果稍微偏大，斯特劳哈尔数与其他方法相比偏小，但与文献 28 结果接近；2）Re=1103时，与 Fluent 计算结果较为吻合，但比其他方法（尤其比试验）所得结果偏大，斯特劳哈尔数基本一致；3）Re=1105时，升力和阻力系数计算结果与其他文献吻合较好，斯特劳哈尔数基本一致。综上，雷诺数低于 1105时，黏性涡域法能够较0201030405060708090100Ut/D(b)Re=11031.51.00.500.51.01.52.02.5CD,CLCD(本文)CL(本文)1.00.500201030405060708090100Ut/D(a)Re=21020.51.01.52.02.5CD,CLCD(本文)CL(本文)0201030405060708090100Ut/D(c)Re=11051.51.00.500.51.01.52.02.5CD,CLCDCLCD Wille试验21CL Wille试验21CD Patel试验23CL Patel试验23图 6 不同雷诺数静止圆柱气动力时程Fig.6 Aerodynamic time histories of a static cylinderwith different Reynolds numbersbers112空气动力学学报第41卷为准确地模拟静止圆柱绕流流场和计算升力/阻力系数。2.2 振动圆柱绕流模拟模拟圆柱沿垂直来流方向做受迫振动，运动方程为 y=A sin(2fet)，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 模式，且涡量场轮廓彼此十分相似。44t*=0(a)涡量生成44t*=1t*=2t*=3t*=4t*=5t*=6(b)尾流未完全发展时涡量演变44t*=50t*=85(c)尾流完全发展时涡量演变图 7 振动圆柱涡量生成及演变Fig.7 Generation and evolution of vorticity aroundan oscillating cylindering cylinder将本文方法与其他方法29-33的计算结果进行对比（图 9图 10），可见，在不同频率和振幅下，虽然气动力时程曲线略有不同，但大体趋势基本保持一致。为更加直观地进行对比，图 11 给出了阻力系数平均值 CDM的对比情况。2.3 不同风向角方柱绕流模拟目前大跨桥梁塔柱均以方形或在其基础上进行切角或凹槽处理的截面为主。与圆柱截面不同，对于方柱这种边缘锋利的截面，分离点可能固定在拐角处，边界层在其前方两拐角点分离，对雷诺数的依赖性也小得多，阻力系数也在一定雷诺数范围内基本不变。然而不同长宽比的矩形截面以及桥塔受到不同的风向角时都会产生不同的流场，因此有必要验证方柱绕流流场及气动系数变化规律。实际应用中，桥塔断面长宽比接近于 1，且研究方柱的学者较多，因此本节以方柱为例进行绕流模拟。(a)0时刻1.51.0(b)尾流未完全发展1.51.0t*=1t*=2t*=3t*=4t*=5t*=6(c)尾流完全发展1.51.0t*=50t*=85图 8 振动圆柱压力云图Fig.8 Pressure contours aroundan oscillating cylinderng cylinder第9期赵国辉等：基于黏性涡域法的二维钝体绕流模拟113利用黏性涡域法模拟计算 045风向角下方柱绕流流场及气动系数，方柱边长取 L=1m，雷诺数为1000，时间步长 0.01s。图 12 给出了 0、15和 45风向角下的方柱涡量图。从图中可以清楚观察到，旋涡上下交替脱落形成卡门涡街。图 13 为风向角为 0、15和 45时的升/阻力系数时程曲线。图 14 为不同方法、各风向角下方柱升/阻力系数结果对比。图 15 为风向角 15时剪切层示意图。可以看出：本文方法的计算结果及趋势与其他方法得到的基本一致；对于阻力系数，杨素哲34和 Naudascher 等35计算的最小值1.00.500.51.01.52.02.53.00201030405060708090100Ut/D(a)f*=0.75CD,CLCD(本文)CL(本文)CD 梁亮文30CL 梁亮文3032101234CD,CL0201030405060708090100Ut/D(b)f*=1.0CD(本文)CL(本文)CD 邓跃29CL 邓跃2932101234CD,CL0201030405060708090100Ut/D(c)f*=1.2CD(本文)CL(本文)CD 梁亮文30CL 梁亮文30图 9 振动圆柱气动力时程（A*=0.2）Fig.9 Aerodynamic timeoscilhistories ofan oscillating cylinder（A*=0.2）32101234CD,CL0201030405060708090100Ut/D(a)f*=0.70CD(本文)CL(本文)CD 梁亮文30CL 梁亮文3032101234CD,CL0201030405060708090100Ut/D(b)f*=0.95CD(本文)CL(本文)CD 梁亮文30CL 梁亮文3032101234CD,CL0201030405060708090100Ut/D(c)f*=1.20CD(本文)CL(本文)CD 梁亮文30CL 梁亮文30图 10 振动圆柱气动力时程（A*=0.4）Fig.10 Oscillating cylinder aerodynamic time histohistories ofan oscillating cylinder（A*=0.4）er（A*=0.4）114空气动力学学报第41卷出现在风向角为 10时，Norberg36试验的最小值出现在风向角为 12.5时，本文和 Taylor 等37计算的最小值出现在风向角为 15时。对于升力系数，杨素哲34计算的最小值出现在风向角为 10时，Norberg36试验的最小值出现在风向角为 12.5时，Naudascher 等35计算的最小值出现在风向角 13左右，本文和 Taylor等37计算的最小值出现在风向角为 15时。综上，可以得出结论，本文方法能够有效地模拟方柱绕流。图 14 还显示，升力系数和阻力系数均随风向角的增大先减小后增大，尝试分析原因如下。当风向角0.70.80.91.01.11.21.00.500.51.01.52.02.53.0Tang33邓跃29梁亮文30王凯鹏31本文f*(a)A*=0.20.70.80.91.01.11.2Tang33Nobari32梁亮文30王凯鹏31本文1.00.500.51.01.52.02.53.0f*(b)A*=0.4CDMCDM图 11 振动圆柱结果对比Fig.11 Oscillating cylinder results comparison(c)45风向角图 12 不同风向角下方柱绕流涡量Fig.12 Vorticity around a rectangular cylinder underdifferent wind direction angles(a)0风向角(b)15风向角32101230201030405060708090100Ut/D(a)0风向角CD,CLCDCL32101230201030405060708090100Ut/D(b)15风向角CD,CLCDCL第9期赵国辉等：基于黏性涡域法的二维钝体绕流模拟115较小时，随着风向角增大，平均阻力系数减小，主要是从 B 点脱落的剪切层间接接触 C 点；当风向角约为 15时，则该剪切层完全附着，最初与 B 分离的剪切层现在与 C 分离，因而产生更窄的尾迹和更低的平均阻力系数；当风向角继续增大时，尾流宽度又会重新变宽，从而平均阻力系数又会变大。虽然不同方法和试验得到的最小平均阻力系数风向角不相同（相差在 5之内），但黏性涡域法很好地再现了这种趋势。升力系数随风向角变化也有相同的趋势。分离的涡泡产生比 AD 面更大的局部吸力，因而产生负升力系数；当风向角为 15时剪切层完全附着，此时负升力系数最大；随着风向角进一步增大，BC 面上的局部吸力减小，升力系数逐渐增加。3 结论黏性涡域法具备涡方法无网格化的优点，同时又考虑了流体的黏性效应，因此是目前较为准确且高效的离散涡方法。通过黏性涡域法模拟不同雷诺数下圆柱绕流、不同无量纲频率和无量纲振幅的振动圆柱气动力时程、不同风向角下的方柱绕流，结果表明：1）计算不同雷诺数下（21021105）静止圆柱的气动力并模拟其流场，得到的升/阻力时程曲线和斯特劳哈尔数与其他文献计算结果基本一致；2）对不同无量纲频率无量纲振幅的振动圆柱进行模拟，得到了经典的卡门涡街，黏性涡域法计算得到的阻力系数平均值与其他文献结果基本一致，趋势相同；3）不同风向角下的方柱绕流模拟结果：当风向角大于 0、小于 15时，升力系数和阻力系数均减小；当风向角大于 15小于 45时，升力系数和阻力系数则均增大。这主要与此时剪切层的脱落与再附着有关。4）黏性涡域法在雷诺数为 21021105条件下对圆柱绕流及不同风向角下的方柱绕流模拟具备足够的计算精度。32101230201030405060708090100Ut/D(c)45风向角CD,CLCDCL图 13 不同风向角下方柱绕流升、阻力时程Fig.13 Time histories of lift and drag coefficients arounda rectangular cylinder under different wind direction angles010203040501.41.61.82.02.22.4本文结果杨素哲34Taylor37Naudascher35Norberg36风向角/()(a)阻力系数对比CD01020304050本文结果杨素哲34Taylor37Naudascher35Norberg36风向角/()(b)升力系数对比1.00.80.60.40.200.20.4CL图 14 不同方法各风向角方柱升、阻力系数对比Fig.14 Comparison of rectangular cylinder drag andlift coefficients at different wind direction anglesby on angles by different methods(a)t*=0.5(b)t*=1.5ABCD图 15 风向角 15时剪切层Fig.15 Shear layer when the wind direction angle is 15116空气动力学学报第41卷 4 后期研究方向展望黏性涡域法虽然具有很多优势，但也存在一定的局限性。由此考虑继续进一步开展如下研究工作：1）黏性涡域法目前主要应用于二维计算。若要实现三维模拟，则需要借助其他有限元软件来实现分离系统迭代方法模拟。在未来的研究中，将进一步进行纯拉格朗日三维涡方法的研究。2）随着雷诺数的增加，流场中的粒子数会相应增加，粒子数的增加将导致计算量增大，需要进一步研究高雷诺数下面元划分尺度与计算效率之间的关系38，以及针对流体的对流扩散研究基于伽辽金法等39更快速的算法。3）目前该方法仅适用均匀来流的模拟，尚且未考虑来流为湍流的工况。未来仍需进一步研究实现湍流的模拟。参 考 文 献：CHORINAJ.NumericalstudyofslightlyviscousflowJ.JournalofFluidMechanics,1973,57(4):785796.doi:10.1017/s00221120730020161DEGOND P,MAS-GALLIC S.The weighted particle method forconvection-diffusion equations.Part 1:The case of an isotropicviscosityJ.MathematicsofComputation,1989,53(188):485507.doi:10.1090/s0025-5718-1989-0983559-92SHANKARS,DOMMELENLVAN.Anewdiffusionprocedureforvortex methodsJ.Journal of Computational Physics,1996,127(1):88109.doi:10.1006/jcph.1996.01603OGAMIY,AKAMATSUT.ViscousflowsimulationusingthediscretevortexmodelthediffusionvelocitymethodJ.Computers&Fluids,1991,19(3-4):433441.doi:10.1016/0045-7930(91)90068-S4PLOUMHANS P,WINCKELMANS G S.Vortex methods for high-resolution simulations of viscous flow past bluff bodies of generalgeometryJ.JournalofComputationalPhysics,2000,165(2):354406.doi:10.1006/jcph.2000.66145PLOUMHANSP,WINCKELMANSGS,SALMONJK,etal.Vortexmethodsfordirectnumericalsimulationofthree-dimensionalbluffbodyflows:applicationtothesphereatRe=300,500,and1000J.JournalofComputationalPhysics,2002,178(2):427463.doi:10.1006/jcph.2002.70356CLARKENR,TUTTYOR.Constructionandvalidationofadiscretevortex method f