趙國(guó)輝，王治超，郝鍵銘，高廣中

(1.長(zhǎng)安大學(xué) 公路學(xué)院，西安 710064；2.武漢市政工程設(shè)計(jì)研究院有限責(zé)任公司，武漢 430023)

0 引言

渦方法是計(jì)算流體力學(xué)中的一類(lèi)無(wú)網(wǎng)格拉格朗日方法，該方法可以將計(jì)算資源集中到占比較小的渦流區(qū)域，因此可以提高計(jì)算效率。由于計(jì)算點(diǎn)的運(yùn)動(dòng)與流體顆粒的運(yùn)動(dòng)接近或重合，用渦來(lái)描述流場(chǎng)，使得渦域具有很好的可視化效果，便于分析旋渦的形成機(jī)理。渦方法不僅能模擬流體中的大尺度擬序結(jié)構(gòu)，也能模擬小尺度渦結(jié)構(gòu)，且在模擬高雷諾數(shù)時(shí)不需要湍流模型。此外，由于沒(méi)有計(jì)算網(wǎng)格，渦方法可以被應(yīng)用于更廣泛的幾何形狀。

最初的渦方法研究未考慮黏性效應(yīng)，僅適用于理想流體。但在鈍體周?chē)?，不僅分離點(diǎn)的位置未知并且還會(huì)隨時(shí)間變化，這種情況下黏性起主要作用。Chorin[1]將渦量-速度方程分解為對(duì)流項(xiàng)和黏性擴(kuò)散項(xiàng)兩部分，而處理黏性擴(kuò)散的方法有隨機(jī)渦方法和確定性渦方法。常見(jiàn)的隨機(jī)渦方法有隨機(jī)走步法與重采樣方法。隨機(jī)渦方法在計(jì)算過(guò)程中產(chǎn)生隨機(jī)數(shù)，雖然原理和實(shí)現(xiàn)方法均較簡(jiǎn)單，但由于擴(kuò)散的隨機(jī)性，此類(lèi)方法在求解渦量時(shí)存在較大的統(tǒng)計(jì)誤差且收斂速度較低。常見(jiàn)的確定性渦方法有粒子強(qiáng)度交換法（particle strength exchange method,PSE）[2]、渦量重分布法[3]和擴(kuò)散速度法（diffusion velocity method,DVM）[4]。粒子強(qiáng)度交換法中，粒子位置的拉普拉斯算子由周?chē)Ｗ拥姆e分算子代替，其精確度強(qiáng)烈依賴(lài)于用于離散積分的求積規(guī)則。之后進(jìn)一步演化了重新網(wǎng)格化的PSE 方法[5-6]，該方法不使用積分算子因而得以保留無(wú)網(wǎng)格化，但需要在每個(gè)時(shí)間步內(nèi)求解N個(gè)粒子的N個(gè)欠定系統(tǒng)，因而隨著粒子數(shù)量增加，計(jì)算效率顯著降低。這種新PSE 方法對(duì)重新劃分網(wǎng)格的嚴(yán)重依賴(lài)促使了渦量重分布法的發(fā)展。擴(kuò)散速度法[4]提出了擴(kuò)散速度的概念，在該方法中，渦流粒子保持其環(huán)流并根據(jù)速度場(chǎng)運(yùn)動(dòng)，速度場(chǎng)是對(duì)流速度和擴(kuò)散速度的疊加，主要計(jì)算變量為渦量。Clarke 等[7]發(fā)現(xiàn)，擴(kuò)散速度法在距離物體較近時(shí)效果很好，但在距離物體較遠(yuǎn)時(shí)就變成了隨機(jī)走步法。且由于旋渦排斥的非單調(diào)特性，即當(dāng)旋渦彼此接近時(shí)排斥率趨于零，導(dǎo)致旋渦“黏附”，從而不能準(zhǔn)確模擬（尤其是物體表面線(xiàn)附近的）渦量演變。Dynnikova 等[8]認(rèn)為，流域中的渦量演變可以看作是它沿總速度場(chǎng)流線(xiàn)的轉(zhuǎn)移，從而發(fā)展了擴(kuò)散速度法，使得渦量演變更為準(zhǔn)確，并且很好地描述了邊界層，確保了在物體表面附近的擴(kuò)散速度的正確表達(dá)，并將這種方法命名為黏性渦域法（viscous vortex domains method,VVD）。

本文擬從黏性渦域法的發(fā)展、原理、實(shí)例等各方面展開(kāi)闡述，分步介紹黏性渦域法的實(shí)現(xiàn)以及各步的理論公式；利用該方法模擬不同雷諾數(shù)下圓柱繞流、不同無(wú)量綱頻率和無(wú)量綱振幅的振動(dòng)圓柱氣動(dòng)力時(shí)程、不同風(fēng)向角下的方柱繞流，并將計(jì)算結(jié)果與現(xiàn)有文獻(xiàn)的結(jié)果進(jìn)行對(duì)比，驗(yàn)證黏性渦域法的有效性和計(jì)算精度。

1 黏性渦域法

黏性渦域法模擬過(guò)程為：

1）生成渦量。該方法主要計(jì)算變量為渦量，渦量用物體表面線(xiàn)上的渦片代替，渦片強(qiáng)度通過(guò)無(wú)滑移邊界條件用邊界積分方程的形式寫(xiě)出。

2）渦量演變。將渦量從物體表面線(xiàn)轉(zhuǎn)移到流場(chǎng)中，渦流尾流演變模擬的同時(shí)得到流體速度。

3）氣動(dòng)計(jì)算。通過(guò)Uhlman 提出的積分方程求解壓力分布，將物體表面線(xiàn)各點(diǎn)壓力與摩擦力的豎向及水平分量累積求和并無(wú)因次化，可得升力系數(shù)與阻力系數(shù)。

1.1 生成渦量

物體表面的渦量可用強(qiáng)度待定的自由渦片（沿物體表面線(xiàn)的無(wú)限薄的渦量分布）表示，并將其轉(zhuǎn)移至流場(chǎng)中，這些渦片作為流場(chǎng)渦源并根據(jù)控制方程發(fā)展運(yùn)動(dòng)。這就意味著模擬物體表面線(xiàn)影響的渦片應(yīng)為自由渦片而非附 加渦片，為此Lighthill[9]提出，由固體邊界的邊界層效應(yīng)產(chǎn)生了沿固體邊界連續(xù)分布的渦量，該連續(xù)渦流近似用離散的渦片模擬，然后集中在一定強(qiáng)度渦片中的渦量成為渦流尾跡的一部分。

渦片強(qiáng)度用邊界積分的形式表示為[10]：

向量核Q(r,ξ)類(lèi)似于二維格林函數(shù)的梯度：

其中：K為物體表面線(xiàn)；S為流場(chǎng)；γ(ξ)為渦片強(qiáng)度；ξ為物體表面線(xiàn)上的點(diǎn)；γ(r)為附加渦片強(qiáng)度，其等于物體表面速度的切線(xiàn)分量。

為方便表示，定義式（1）等號(hào)右邊為fτ(r)，且fτ(r)=f(r)·τ(r)。n(r)和τ(r)分別是外法線(xiàn)向量和切線(xiàn)向量，且n(r) ×τ(r)=e。由此可得到取決于來(lái)流流速、物體表面線(xiàn)速度和流動(dòng)中的渦量分布的已知函數(shù)f(r)：

其中：V∞為來(lái)流流速；N為流場(chǎng)內(nèi)渦流粒子數(shù)；ri和Гi，i=(1,2,···,N)，分別為渦流粒子的位置和環(huán)量。

為得到邊界積分方程（1）的唯一解，需要引入如下附加條件：

利用Lifanov[11]的方法，將方程（5）與邊界積分方程（1）聯(lián)合求解。求解邊界積分方程最簡(jiǎn)單、最普遍的方法是伽遼金法。但傳統(tǒng)的伽遼金法無(wú)法重建邊界積分方程的解。解的變化越明顯，渦粒子越接近物體表面線(xiàn)，顯著的局部誤差會(huì)導(dǎo)致近壁區(qū)域的速度場(chǎng)重建不正確。為解決該問(wèn)題，Soldatova 等[12]提出了在數(shù)值解中加入一項(xiàng)，該項(xiàng)考慮了由于密集放置的渦粒子的影響而引起的解在一個(gè)面元上的局部行為：

其中：γ*為修正項(xiàng)；?q為基函數(shù)（當(dāng)q=0 時(shí)為常數(shù)，q=1時(shí)為線(xiàn)性，q=2 時(shí)為二次）；γq為未知系數(shù)，可以通過(guò)伽遼金方法確定。為簡(jiǎn)便起見(jiàn)，考慮正交基函數(shù)，當(dāng)ξ在對(duì)應(yīng)面元上時(shí)，?0(ξ)≡1 并且?1和?2都等于1，在面元之外與其相關(guān)的所有基函數(shù)?q都等于零。

通過(guò)節(jié)點(diǎn)劃分對(duì)物體表面線(xiàn)進(jìn)行離散，將物體表面線(xiàn)劃分為多個(gè)節(jié)點(diǎn)，兩節(jié)點(diǎn)間的節(jié)段稱(chēng)為面元，如圖1 所示。選擇離渦流粒子相當(dāng)近的面元，設(shè)面元k∈[kb,ke]，這些面元被曲率為κk的密切圓圓弧所代替，利用共形映射理論為每個(gè)相當(dāng)接近渦流粒子的面元寫(xiě)下精確解：

圖1 物體表面節(jié)點(diǎn)劃分和面元Fig.1 Node division of an airfoil and panels

其中：Γw為位于rw點(diǎn)的渦流粒子環(huán)量；rk(ξ)為密切圓圓弧上一點(diǎn)；nk(ξ)為點(diǎn)rk(ξ)處的外單位法向量；為密切圓的圓心；；Rk為密切圓半徑。

位于點(diǎn)rw渦流粒子的環(huán)量可進(jìn)一步表示為：

其中：Гw,k為分布在第k面元上的總環(huán)量，用來(lái)考慮渦粒子的貢獻(xiàn)；dw,k為從渦流粒子到第k面元的距離；s為弧長(zhǎng)參數(shù)值，，sw對(duì)應(yīng)于渦流粒子位置在第k個(gè)面元上的投影。各物理量如圖2 所示。

圖2 相鄰兩個(gè)曲線(xiàn)面元和靠近的渦流粒子Fig.2 Two adjacent curved panels andapproaching vortex elements

1.2 渦量演變

渦量生成后，需要將渦量從物體表面線(xiàn)轉(zhuǎn)移到流域。在物體表面線(xiàn)處形成渦片的分布渦量被轉(zhuǎn)化為單獨(dú)的渦流粒子，成為旋渦尾流的一部分。不考慮流動(dòng)可壓縮性，就渦量而言，Navier-Stokes 方程可表示為：

其中：V為對(duì)流速度；W為由黏性效應(yīng)引起的擴(kuò)散速度；渦量場(chǎng)Ω只有一個(gè)非零分量，與流動(dòng)平面正交，可寫(xiě)為Ω=Ωe。方程（9）是指流動(dòng)中存在的渦量隨速度移動(dòng)，新的渦量只在流動(dòng)區(qū)域的邊界上（即物體的邊界上）產(chǎn)生。

渦方法模擬無(wú)黏、不可壓縮流的渦流尾流演化時(shí)，需要對(duì)以下常微分方程組進(jìn)行求解：

其中，ri是第i個(gè)渦流粒子的位置。由于渦量只是隨著流速在流域中傳遞，對(duì)流速度V用Biot-Savart 定律由渦量分布計(jì)算得到：

對(duì)于擴(kuò)散速度，有：

式中 υ為運(yùn)動(dòng)黏度。式（12）表明擴(kuò)散速度與運(yùn)動(dòng)黏度υ成正比，并且取決于點(diǎn)r附近流動(dòng)區(qū)域中的渦量分布和流動(dòng)區(qū)域的邊界形狀。

1.3 氣動(dòng)計(jì)算

圖3 為圓柱繞流示意圖，圖中Fp為壓力、τ為物體表面的切向應(yīng)力。物體表面受到流體的壓力和剪切應(yīng)力，由Uhlman 提出的積分方程，可以求解流場(chǎng)壓力分布，從而計(jì)算壓力：

圖3 圓柱繞流示意圖Fig.3 Flow around a cylinderer

其中：β為常數(shù)，二維問(wèn)題的計(jì)算域內(nèi)部β=2p，計(jì)算域邊界上β=p；B為停滯熱焓；p為壓強(qiáng)；Z為計(jì)算域內(nèi)部；Y為計(jì)算域邊界；G為標(biāo)量格林函數(shù)。

設(shè)在物面處有Nw個(gè)物面渦元，且第j個(gè)物面渦元的環(huán)量為Γj，在計(jì)算域內(nèi)部Z有Nz個(gè)自由渦，且第k個(gè)自由渦的環(huán)量為Γk，則物面i點(diǎn)處的壓力積分方程的離散形式為：

物體靜止時(shí)，有：

物體表面渦和流場(chǎng)中自由渦的渦強(qiáng)與環(huán)量關(guān)系為：

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

式（20）中，等號(hào)右端第一項(xiàng)表示物體表面渦元的存在對(duì)物面壓強(qiáng)的貢獻(xiàn)，第二項(xiàng)代表流場(chǎng)中自由渦的運(yùn)動(dòng)對(duì)物面壓強(qiáng)的貢獻(xiàn)。求解式（20）得到Bi，再將其代入式（14），可得到物面上i點(diǎn)的壓強(qiáng)：

對(duì)于物體表面的剪切應(yīng)力τw，假定邊界層內(nèi)速度是線(xiàn)性分布的，則物體所受的剪切應(yīng)力可通過(guò)下式求解：

綜上，物體所受的流體力為：

式中，en和eτ分別代表物體表面法向量和切向量。將物體所受流體力沿順流向和橫流向分解，可得到物體所受的阻力FD和升力FL，無(wú)因次化后可得阻力系數(shù)CD和升力系數(shù)CL：

2 數(shù)值模擬

基于黏性渦域法，在Linux 系統(tǒng)下利用Lua 語(yǔ)言編制圓柱繞流的計(jì)算程序，利用Gnuplot 實(shí)現(xiàn)流場(chǎng)可視化。進(jìn)行了低雷諾數(shù)和高雷諾數(shù)下靜止和振動(dòng)圓柱的模擬。

2.1 靜止圓柱繞流模擬

分別取雷諾數(shù)Re=2 × 102、1 × 103、1 × 105進(jìn)行模擬計(jì)算和驗(yàn)證，根據(jù)多次試算綜合考慮計(jì)算精度與效率，面元長(zhǎng)度取0.01 m。圖4～圖6 是基于黏性渦域法模擬的靜止圓柱渦流演變過(guò)程、壓力云圖及氣動(dòng)力時(shí)程。用點(diǎn)模擬渦流尾流的離散渦渦流粒子位置，分別用紅色、藍(lán)色點(diǎn)表示正、負(fù)循環(huán)渦流。t*=Ut/D為無(wú)量綱時(shí)間。阻力系數(shù)CD=FD/(0.5ρDU2)，升力系數(shù)CL=FL/(0.5ρDU2)，ρ為流體密度，F(xiàn)D和FL分別為柱體所受的阻力和升力。

圖4 靜止圓柱渦量生成及演變Fig.4 Generation and evolution of vorticity around a static cylinderc cylinder

圖5 靜止圓柱壓力云圖Fig.5 Static cylinder pressure contour

圖6 不同雷諾數(shù)靜止圓柱氣動(dòng)力時(shí)程Fig.6 Aerodynamic time histories of a static cylinder with different Reynolds numbersbers

由圖4 可見(jiàn)，VVD 能夠模擬渦流演變過(guò)程，且利用該方法模擬的卡門(mén)渦街與周志勇等[13]利用離散渦方法、王亞玲等[14]利用計(jì)算流體軟件CFX-4、桑文慧等[15]利用SIMPLEC 算法和張偉偉等[16]的模擬結(jié)果基本一致。時(shí)間步長(zhǎng)為0.05、計(jì)算步為20 000 時(shí)，VVD 方法的計(jì)算時(shí)間約為10 min。

由圖6 可知，在不同雷諾數(shù)下，本文方法與文獻(xiàn)[21,23]的氣動(dòng)力時(shí)程計(jì)算結(jié)果基本一致。為進(jìn)一步驗(yàn)證本文方法的準(zhǔn)確性，將氣動(dòng)計(jì)算結(jié)果與現(xiàn)有文獻(xiàn)結(jié)果進(jìn)行對(duì)比，結(jié)果見(jiàn)表1～表3，表中數(shù)據(jù)均為旋渦脫落穩(wěn)定后的數(shù)據(jù)。

表1 不同方法氣動(dòng)力計(jì)算結(jié)果對(duì)比（Re=2 × 102）Table 1 Comparison of aerodynamics obtained by different methods (Re=2 × 102)

表2 不同方法氣動(dòng)力計(jì)算結(jié)果對(duì)比（Re=1 × 103）Table 2 Comparison of aerodynamics obtained by different methods (Re=1 × 103)

表3 不同方法氣動(dòng)力計(jì)算結(jié)果對(duì)比（Re=1 × 105）Table 3 Comparison of aerodynamics obtained by different methods (Re=1 × 105)

通過(guò)氣動(dòng)力時(shí)程曲線(xiàn)以及與各方法氣動(dòng)力系數(shù)幅值對(duì)比，可得：

1）Re=2 × 102時(shí)，本文方法與現(xiàn)有文獻(xiàn)中各方法計(jì)算結(jié)果基本一致，其中與Fluent 計(jì)算結(jié)果最為吻合，阻力系數(shù)比離散渦方法結(jié)果稍微偏大，斯特勞哈爾數(shù)與其他方法相比偏小，但與文獻(xiàn)[28]結(jié)果接近；

2）Re=1 × 103時(shí)，與Fluent 計(jì)算結(jié)果較為吻合，但比其他方法（尤其比試驗(yàn)）所得結(jié)果偏大，斯特勞哈爾數(shù)基本一致；

3）Re=1 × 105時(shí)，升力和阻力系數(shù)計(jì)算結(jié)果與其他文獻(xiàn)吻合較好，斯特勞哈爾數(shù)基本一致。

綜上，雷諾數(shù)低于1 × 105時(shí)，黏性渦域法能夠較為準(zhǔn)確地模擬靜止圓柱繞流流場(chǎng)和計(jì)算升力/阻力系數(shù)。

2.2 振動(dòng)圓柱繞流模擬

模擬圓柱沿垂直來(lái)流方向做受迫振動(dòng)，運(yùn)動(dòng)方程為y=Asin(2πfet)，fe為圓柱的振動(dòng)頻率，t為圓柱受迫運(yùn)動(dòng)的運(yùn)動(dòng)時(shí)刻。無(wú)量綱振動(dòng)頻率f*=fe/fs分別取0.5、0.7、0.75、0.95、1.0、1.2，fs為靜止圓柱自然脫落頻率。無(wú)量綱振幅A*=A/D分別取0.2、0.4，A為振蕩幅值。振動(dòng)圓柱渦流演變過(guò)程、壓力云圖及氣動(dòng)力時(shí)程見(jiàn)圖7～圖10。圖7 中，尾流呈現(xiàn)出經(jīng)典的卡門(mén)渦街。在圓柱單個(gè)振蕩周期內(nèi)，渦從圓柱兩側(cè)交替脫落，即2S 模式，且渦量場(chǎng)輪廓彼此十分相似。

圖7 振動(dòng)圓柱渦量生成及演變Fig.7 Generation and evolution of vorticity around an oscillating cylindering cylinder

圖8 振動(dòng)圓柱壓力云圖Fig.8 Pressure contours around an oscillating cylinderng cylinder

將本文方法與其他方法[29-33]的計(jì)算結(jié)果進(jìn)行對(duì)比（圖9～圖10），可見(jiàn)，在不同頻率和振幅下，雖然氣動(dòng)力時(shí)程曲線(xiàn)略有不同，但大體趨勢(shì)基本保持一致。為更加直觀地進(jìn)行對(duì)比，圖11 給出了阻力系數(shù)平均值CDM的對(duì)比情況。

圖9 振動(dòng)圓柱氣動(dòng)力時(shí)程（A＊=0.2）Fig.9 Aerodynamic timeoscilhistories of an oscillating cylinder（A＊=0.2）

圖10 振動(dòng)圓柱氣動(dòng)力時(shí)程（A＊=0.4）Fig.10 Oscillating cylinder aerodynamic time histohistories of an oscillating cylinder（A＊=0.4）er（A＊=0.4）

圖11 振動(dòng)圓柱結(jié)果對(duì)比Fig.11 Oscillating cylinder results comparison

2.3 不同風(fēng)向角方柱繞流模擬

目前大跨橋梁塔柱均以方形或在其基礎(chǔ)上進(jìn)行切角或凹槽處理的截面為主。與圓柱截面不同，對(duì)于方柱這種邊緣鋒利的截面，分離點(diǎn)可能固定在拐角處，邊界層在其前方兩拐角點(diǎn)分離，對(duì)雷諾數(shù)的依賴(lài)性也小得多，阻力系數(shù)也在一定雷諾數(shù)范圍內(nèi)基本不變。然而不同長(zhǎng)寬比的矩形截面以及橋塔受到不同的風(fēng)向角時(shí)都會(huì)產(chǎn)生不同的流場(chǎng)，因此有必要驗(yàn)證方柱繞流流場(chǎng)及氣動(dòng)系數(shù)變化規(guī)律。實(shí)際應(yīng)用中，橋塔斷面長(zhǎng)寬比接近于1，且研究方柱的學(xué)者較多，因此本節(jié)以方柱為例進(jìn)行繞流模擬。

利用黏性渦域法模擬計(jì)算0°～45°風(fēng)向角下方柱繞流流場(chǎng)及氣動(dòng)系數(shù)，方柱邊長(zhǎng)取L=1 m，雷諾數(shù)為1 000，時(shí)間步長(zhǎng)0.01 s。圖12 給出了0°、15°和45°風(fēng)向角下的方柱渦量圖。從圖中可以清楚觀察到，旋渦上下交替脫落形成卡門(mén)渦街。圖13 為風(fēng)向角為0°、15°和45°時(shí)的升/阻力系數(shù)時(shí)程曲線(xiàn)。圖14 為不同方法、各風(fēng)向角下方柱升/阻力系數(shù)結(jié)果對(duì)比。圖15 為風(fēng)向角15°時(shí)剪切層示意圖?？梢钥闯觯罕疚姆椒ǖ挠?jì)算結(jié)果及趨勢(shì)與其他方法得到的基本一致；對(duì)于阻力系數(shù)，楊素哲[34]和Naudascher 等[35]計(jì)算的最小值出現(xiàn)在風(fēng)向角為10°時(shí)，Norberg[36]試驗(yàn)的最小值出現(xiàn)在風(fēng)向角為12.5°時(shí)，本文和Taylor 等[37]計(jì)算的最小值出現(xiàn)在風(fēng)向角為15°時(shí)。對(duì)于升力系數(shù)，楊素哲[34]計(jì)算的最小值出現(xiàn)在風(fēng)向角為10°時(shí)，Norberg[36]試驗(yàn)的最小值出現(xiàn)在風(fēng)向角為12.5°時(shí)，Naudascher 等[35]計(jì)算的最小值出現(xiàn)在風(fēng)向角13°左右，本文和Taylor等[37]計(jì)算的最小值出現(xiàn)在風(fēng)向角為15°時(shí)。綜上，可以得出結(jié)論，本文方法能夠有效地模擬方柱繞流。

圖12 不同風(fēng)向角下方柱繞流渦量Fig.12 Vorticity around a rectangular cylinder under different wind direction angles

圖13 不同風(fēng)向角下方柱繞流升、阻力時(shí)程Fig.13 Time histories of lift and drag coefficients around a rectangular cylinder under different wind direction angles

圖14 不同方法各風(fēng)向角方柱升、阻力系數(shù)對(duì)比Fig.14 Comparison of rectangular cylinder drag and lift coefficients at different wind direction angles by on angles by different methods

圖15 風(fēng)向角15°時(shí)剪切層Fig.15 Shear layer when the wind direction angle is 15°

圖14 還顯示，升力系數(shù)和阻力系數(shù)均隨風(fēng)向角的增大先減小后增大，嘗試分析原因如下。當(dāng)風(fēng)向角較小時(shí)，隨著風(fēng)向角增大，平均阻力系數(shù)減小，主要是從B 點(diǎn)脫落的剪切層間接接觸C 點(diǎn)；當(dāng)風(fēng)向角約為15°時(shí)，則該剪切層完全附著，最初與B 分離的剪切層現(xiàn)在與C 分離，因而產(chǎn)生更窄的尾跡和更低的平均阻力系數(shù)；當(dāng)風(fēng)向角繼續(xù)增大時(shí)，尾流寬度又會(huì)重新變寬，從而平均阻力系數(shù)又會(huì)變大。雖然不同方法和試驗(yàn)得到的最小平均阻力系數(shù)風(fēng)向角不相同（相差在5°之內(nèi)），但黏性渦域法很好地再現(xiàn)了這種趨勢(shì)。升力系數(shù)隨風(fēng)向角變化也有相同的趨勢(shì)。分離的渦泡產(chǎn)生比AD 面更大的局部吸力，因而產(chǎn)生負(fù)升力系數(shù)；當(dāng)風(fēng)向角為15°時(shí)剪切層完全附著，此時(shí)負(fù)升力系數(shù)最大；隨著風(fēng)向角進(jìn)一步增大，BC 面上的局部吸力減小，升力系數(shù)逐漸增加。

3 結(jié)論

黏性渦域法具備渦方法無(wú)網(wǎng)格化的優(yōu)點(diǎn)，同時(shí)又考慮了流體的黏性效應(yīng)，因此是目前較為準(zhǔn)確且高效的離散渦方法。通過(guò)黏性渦域法模擬不同雷諾數(shù)下圓柱繞流、不同無(wú)量綱頻率和無(wú)量綱振幅的振動(dòng)圓柱氣動(dòng)力時(shí)程、不同風(fēng)向角下的方柱繞流，結(jié)果表明：

1）計(jì)算不同雷諾數(shù)下（2× 102～1 × 105）靜止圓柱的氣動(dòng)力并模擬其流場(chǎng)，得到的升/阻力時(shí)程曲線(xiàn)和斯特勞哈爾數(shù)與其他文獻(xiàn)計(jì)算結(jié)果基本一致；

2）對(duì)不同無(wú)量綱頻率無(wú)量綱振幅的振動(dòng)圓柱進(jìn)行模擬，得到了經(jīng)典的卡門(mén)渦街，黏性渦域法計(jì)算得到的阻力系數(shù)平均值與其他文獻(xiàn)結(jié)果基本一致，趨勢(shì)相同；

3）不同風(fēng)向角下的方柱繞流模擬結(jié)果：當(dāng)風(fēng)向角大于0°、小于15°時(shí)，升力系數(shù)和阻力系數(shù)均減小；當(dāng)風(fēng)向角大于15°小于45°時(shí)，升力系數(shù)和阻力系數(shù)則均增大。這主要與此時(shí)剪切層的脫落與再附著有關(guān)。

4）黏性渦域法在雷諾數(shù)為2 × 102～1 × 105條件下對(duì)圓柱繞流及不同風(fēng)向角下的方柱繞流模擬具備足夠的計(jì)算精度。

4 后期研究方向展望

黏性渦域法雖然具有很多優(yōu)勢(shì)，但也存在一定的局限性。由此考慮繼續(xù)進(jìn)一步開(kāi)展如下研究工作：

1）黏性渦域法目前主要應(yīng)用于二維計(jì)算。若要實(shí)現(xiàn)三維模擬，則需要借助其他有限元軟件來(lái)實(shí)現(xiàn)分離系統(tǒng)迭代方法模擬。在未來(lái)的研究中，將進(jìn)一步進(jìn)行純拉格朗日三維渦方法的研究。

2）隨著雷諾數(shù)的增加，流場(chǎng)中的粒子數(shù)會(huì)相應(yīng)增加，粒子數(shù)的增加將導(dǎo)致計(jì)算量增大，需要進(jìn)一步研究高雷諾數(shù)下面元?jiǎng)澐殖叨扰c計(jì)算效率之間的關(guān)系[38]，以及針對(duì)流體的對(duì)流擴(kuò)散研究基于伽遼金法等[39]更快速的算法。

3）目前該方法僅適用均勻來(lái)流的模擬，尚且未考慮來(lái)流為湍流的工況。未來(lái)仍需進(jìn)一步研究實(shí)現(xiàn)湍流的模擬。