趙國輝，王治超，郝鍵銘，高廣中

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

0 引言

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

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

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

1 黏性渦域法

黏性渦域法模擬過程為：

1）生成渦量。該方法主要計算變量為渦量，渦量用物體表面線上的渦片代替，渦片強度通過無滑移邊界條件用邊界積分方程的形式寫出。

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

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

1.1 生成渦量

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

渦片強度用邊界積分的形式表示為[10]：

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

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

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

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

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

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

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

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

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

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

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

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

圖2 相鄰兩個曲線面元和靠近的渦流粒子Fig.2 Two adjacent curved panels andapproaching vortex elements

1.2 渦量演變

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

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

渦方法模擬無黏、不可壓縮流的渦流尾流演化時，需要對以下常微分方程組進行求解：

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

對于擴散速度，有：

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

1.3 氣動計算

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

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

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

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

物體靜止時，有：

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

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

式（20）中，等號右端第一項表示物體表面渦元的存在對物面壓強的貢獻，第二項代表流場中自由渦的運動對物面壓強的貢獻。求解式（20）得到Bi，再將其代入式（14），可得到物面上i點的壓強：

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

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

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

2 數(shù)值模擬

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

2.1 靜止圓柱繞流模擬

分別取雷諾數(shù)Re=2 × 102、1 × 103、1 × 105進行模擬計算和驗證，根據(jù)多次試算綜合考慮計算精度與效率，面元長度取0.01 m。圖4～圖6 是基于黏性渦域法模擬的靜止圓柱渦流演變過程、壓力云圖及氣動力時程。用點模擬渦流尾流的離散渦渦流粒子位置，分別用紅色、藍色點表示正、負循環(huán)渦流。t*=Ut/D為無量綱時間。阻力系數(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ù)靜止圓柱氣動力時程Fig.6 Aerodynamic time histories of a static cylinder with different Reynolds numbersbers

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

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

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

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

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

通過氣動力時程曲線以及與各方法氣動力系數(shù)幅值對比，可得：

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

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

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

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

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 中，尾流呈現(xiàn)出經(jīng)典的卡門渦街。在圓柱單個振蕩周期內(nèi)，渦從圓柱兩側(cè)交替脫落，即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]的計算結(jié)果進行對比（圖9～圖10），可見，在不同頻率和振幅下，雖然氣動力時程曲線略有不同，但大體趨勢基本保持一致。為更加直觀地進行對比，圖11 給出了阻力系數(shù)平均值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 振動圓柱結(jié)果對比Fig.11 Oscillating cylinder results comparison

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

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

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

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

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

圖14 不同方法各風(fēng)向角方柱升、阻力系數(shù)對比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°時剪切層Fig.15 Shear layer when the wind direction angle is 15°

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

3 結(jié)論

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

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

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

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

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

4 后期研究方向展望

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

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

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

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