高雷諾數下水翼渦發放頻率預報方法

侯磊，丁云峰，王晴，宗智，孫雷，姜宜辰*

1 海軍裝備部駐武漢地區第二軍事代表室，湖北武漢430064

2 大連理工大學船舶工程學院，遼寧大連116024

3 中國艦船研究設計中心，湖北武漢430064

4 船舶振動噪聲重點實驗室，湖北武漢430064

0 引 言

渦激振動（VIV）被視為流體力學、結構力學、計算流體力學和聲學等不同領域的重點問題。這些領域中關于VIV 的研究可見于Rockwell［1］，Williamson 和Govardhan［2］，Sarpkaya［3］，Bearman［4］和Assi 等［5］的綜述中。流體與結構物相互作用會增加結構的振動，在特定工況下會造成結構破壞。此外，如果泄渦頻率與結構的某一階固有頻率相等，共振會顯著增加結構的振幅，結構的振動位移會俘獲流體激勵，造成所謂的“鎖定”現象。Ausoni 等［6］在研究二維鈍體水翼時發現，如果結構運動位于共振區間之外，泄渦頻率遵循斯特勞哈爾定律，即“脫離鎖定”。在鎖定狀態下，彎曲的渦量線變成平行于尾端的直線，渦強也會增大［7］。因此，研究水翼的繞流特性及渦發放機理，并以此為基礎，對水翼泄渦流致振動進行探究是十分必要的。

Blake 和Gershfeld［8］在20 世紀80 年代全面回顧了前人關于水翼的研究工作。這些試驗說明尾渦泄放在對稱翼尾和鈍體翼尾同時存在。Huerre與Monkewitz［9］和Oertel［10］分別給出了流動不穩定性和鈍體尾流的概述。Lotfy 和Rockwell［11］做了鈍體振動翼尾尾渦泄放研究，為之后的泄渦研究提供了參考。然而，鈍體渦發放的研究結果并不能直接轉換到流線型體，因為相比于鈍體，流線型體首部和尾部的較大距離抑制了窄渦波動和首部停滯點之間的相互作用。

Ohmi 等［12］對機翼在大攻角下作強迫旋轉和平移振動時渦的形成進行了深入的試驗探究。Morikawa 和Gr?nig［13］則研究了此時旋渦的生成及其結構。Huang 等［14］對NACA 0012 翼型進行了試驗分析，并探討了翼根連接處與翼梢對渦激振動的影響。Jung 等［15］通過試驗對低雷諾數下機翼繞流形成的渦街進行分析，發現機翼振動情況下的脫渦頻率相對于機翼靜止不動時要小。Davidson實驗室針對水翼做了大量的試驗研究，主要集中在兩自由度理想翼型方面，觀察測量水翼在水洞試驗中的振動、空泡、壓力分布等。Chae 等［16-17］對兩自由度下的二維水翼進行了系統的研究，探討了附連水質量、水動力阻尼、非線性流固耦合響應等對水翼穩定性的影響。

利用離散渦方法進行水翼渦發放頻率特性的研究，在國內外并不多見。現有的研究大部分基于傳統的經驗模型或數值方法。經驗模型，如尾流振子模型，只能用來求解結構的渦激振動響應，無法對尾流場的特征進行分析，且求解結果僅在數量級上接近真實值；而傳統的數值方法，如有限體積、有限差分、有限元、大渦模擬和直接數值模擬等計算方法和模型在求解VIV 時需要巨量網格，難以獲得足夠高的計算效率和精度［18-19］。

離散渦方法是拉格朗日體系下的一種渦動力學算法，該方法由集中渦模型發展而來，其基本思想是在物體表面直接產生離散渦元，渦元在流場中有限的區域內通過相互作用來模擬渦量場，對渦量場進行積分可以得到速度場，進而可以得到整個流場結構。該方法不依賴于流場網格，在渦量集中區渦尺度分辨率高，算法簡單，程序易于并行化處理，可以選用大時間步長，由于整個求解過程中無需計算全部流場信息，相對于傳統網格法，其計算量大大降低，極為適合高雷諾數情況下結構繞流和渦激振動問題的求解［20-22］。按照離散的渦元基本形狀，離散渦（DVM）可以分為點渦法、渦片法、渦絲法、渦桿法、渦團法和混合方法等；根據初始新生渦強度確定方式的不同，可以分為面元法和渦量—流函數法；根據粘性的處理方式，可以分為面積擴散渦方法、隨機渦方法和確定性渦方法等；根據渦元誘導方式，可以分為直接法、格子渦法和快速多極子展開法；根據邊界處理方式，可以分為鏡像渦法、渦元法、容量矩陣法等［23-26］。

本文將以二維NACA 66 系列水翼為例，首先通過與試驗數據對比，驗證離散渦方法對于模擬高雷諾數下渦發放頻率預報的適用性和準確性。然后，應用此方法，模擬不同流速下NACA 66（MOD）水翼的渦發放情況，并分析流速對水翼泄渦頻率的改變情況。最后，在某一相同高雷諾數下，簡要探討水翼攻角對渦發放頻率的影響規律，以便為分析螺旋槳葉剖面渦發放頻率特性打下堅實基礎。

1 基本原理

假設流體非理想不可壓，由動量守恒和質量守恒組成的控制方程可寫為：

式中：U為流體速度；t 為時間；f 為流體體積力；p為壓力；ρ為流體密度；ν為流體運動粘性系數。對于重力場f=(0，0，-g)，方程兩邊同時取旋度，由于?×f=0，?×(?p)=0，應用渦量—速度關系式ω=?×U，可得

式（3）又稱為旋渦擴散方程或渦量輸運公式。從中可看出，對于非理想不可壓流體，對流項與粘性擴散項之和控制渦量的變化。當流動為二維平面情況時，存在關系式：

將式（6）代入渦量—速度關系式，得到如式（7）所示的控制方程：

但上述方程難以直接求解，根據Chorin［27］所提出算子分裂的思想，動力學控制方程可以分裂為：

1）對流項

2）粘性擴散項

假設流場為無限域，只需要考慮物面邊界條件和遠場邊界條件，其中物面條件為法向不可穿透和切向無滑移條件：

式中：n和τ分別為物面單位法向量和切向量；U為靠近物面處流體粒子速度矢量；Ub為物體速度矢量。在粘性渦中，只要滿足無滑移和不可穿透條件中的任意一個即可，因為粘性渦的這兩個條件是等價的［28］。將渦元布置在物面，點渦由物面產生脫落，并向外運動擴展，即可使物面不可穿透條件自動得到滿足。物面點渦由于渦元間的誘導產生誘導速度，誘導速度和滑移速度相互抵消，從而滿足無滑移條件。

在無窮遠處滿足無窮遠速度條件：

根據Biot-Savart 定理，渦元對無窮遠處的誘導速度為0，于是無窮遠處的流體速度只有來流速度，所以遠場的邊界條件同樣自動滿足。

根據Chorin［27］提出的渦團速度分布公式：

式中：r 為誘導點位置與渦團中心位置的距離；θ為兩點的相位角；σ為渦團半徑；Γ為渦團的渦強，其由流函數方程確定；u 為誘導點切向速度，逆時針為正。將誘導速度代入式（6），并轉化為直角坐標形式，則流函數可表示為

由于粘性作用，流場內的旋渦表現出一定的聚集性，只在流場有限區域內產生相互作用。離散渦方法基于渦量場這一特征，把從物體表面脫落，發展形成的連續旋渦場離散成一系列的渦元，通過模擬渦量場，得到速度場，進而獲得整個流場結構［29-30］。根據這一思想，渦量場內任一點r 處的渦量為：

式中：Γi為位于ri處渦元的渦強；n 為當前流場內渦元總數；δ(r-ri)為狄拉克函數。如果認為流場是無限流場，求解式（8）對流項中的泊松方程可得到渦元誘導速度：

從式（15）可看出，渦元運動速度包括渦元之間的相互誘導部分和來流速度，而相互誘導項就是Biot-Savart定理。

計算出渦元t時刻速度U（r，t）后，經過時間步長dt，新時刻的位置矢量可以表示為

應該注意的是，上式并未考慮由于流體粘性而產生的位移，對于粘性擴散部分，類似熱擴散方程，其基本解是格林函數：

式（17）與均值為0、標準差為s的正態分布的概率密度函數的形式相似：

式中：P 為概率密度函數；η為隨機變量。這樣，Δt時間步內渦元的隨機擴散位移為

對于灰色定權聚類分析方法而言，白化權函數的表達式通常由相關行業的專家給定，決策權ηj的值可通過專家法進行確定，也可通過下述方法確定：

式中：P，Q兩個獨立隨機變量滿足P∈（0，1）和Q∈（0，2 π）。于是t+dt 時刻的渦元位置矢量可以表示為

式中，δr=(δx，δy)。

本文采用的隨機渦方法需要將結構表面離散成節點，同時需要確定適用于該數值模型的新生渦生成位置及數目。首先，將物體看作一個np邊形，其中np為邊界單元數，單元端點為xi（i=1，2，…，n）。配置點xc，i，即滿足速度邊界條件的點。通常有兩種布置方式：一種是將其與相應的外法向量ni布置在單元中點（如圖1（a）所示），并假設物面上的新生渦強Γi集中在配置點上；另一種是將離散的渦元布置在流域內人工設置的邊界層上，即沿端點外法向，距離端點一定距離，這個距離就是當地邊界層厚度d（如圖1（b）所示），渦強同樣集中在配置點。

圖1 配置點布置方式Fig.1 Configuration of discrete nodes arrangement

以圓柱為例探究兩種布置方式對于本文使用的離散渦模型的適用性。一直徑為D=0.2 m 的固定圓柱在遭遇速度為1.0 m/s 的均勻來流時，其新生渦強分布的表達式為

式中，θ為圓柱表面離散節點徑向與流向的夾角。利用上述兩種布置方式對該圓柱繞流問題進行模擬，離散節點數均為n，第2 種布置方式中人工邊界層厚度d 根據經驗設置為D/n。兩種方式得到的新生渦強如圖2 所示。

由圖中對比可知，在物面離散節點數相同時，設置人工邊界層的配置點布置方式得到的新生渦強更接近理論解，因此，在本文的水翼渦發放模擬中，采用第2 種布置方式。

圖2 兩種配置點布置方式新生渦強對比Fig.2 Comparison of new vortex strength between two node arrangements

此外，傳統的CFD 方法通常利用網格對結構或流場進行劃分求解，模擬的結果很大程度上取決于網格劃分的合理性。本文采用的離散渦方法雖為無網格方法，但新生渦的數量同樣會影響計算精度和效率。因此，選取上述固定圓柱進行討論，同樣根據新生渦強的數值對比來衡量最佳新生渦數目。令圓柱直徑D 為0.2 m，流速為1.0 m/s，選取不同數量（n=40，80，120，160，200）的新生渦數，采用第2 種布置方式對該圓柱繞流問題進行模擬，并將新生渦強的結果與理論值進行對比。

圖3 不同新生渦數目求解新生渦強理論值誤差對比Fig.3 Tolerance comparison on vortex strength using various numbers of new-generated vortices

2 算法驗證

2.1 算例模型

采用與Leroux 等［31］試驗相同的NACA 66 水翼模型，弦長C=0.15 m，最大厚度為弦長的12%，最大厚度位于距離前緣45%弦長處。最大拱度為2°，位于弦長中心。試驗中，水翼的寬度為0.191 m，水翼由不銹鋼制作。NACA 66 水翼截面型式如圖4 所示。

圖4 NACA 66 系列水翼截面（最大厚度為弦長的12%，位于x/C=-0.05）Fig.4 Section of NACA 66 series hydrofoil（maximum thickness is 12% of chord length at x/C=-0.05）

2.2 剛性水翼渦發放

采用離散渦法對上述NACA 66 水翼在穩流下的渦發放進行模擬。設置NACA 66 水翼為剛體，攻角保持6°不變，模擬了4 種不同流速情況，對應的雷諾數分別為0.8×106，1.0×106，1.2×106和1.5×106。模擬的總時間長度為30 個周期，其中每個周期T為流體流過一個弦長所需時間（即T=C/U）。

圖5 給出了雷諾數為1.0×106情況下水翼垂向載荷Fy在30 個周期內的變化歷程。由圖可見，垂向載荷時間歷程中存在兩種不同頻率的震蕩，其中周期約為0.5T 的高頻振蕩應是由渦發放直接引發，而周期約為10T 的振蕩則是由渦團之間的相互影響產生。

圖5 垂向載荷時間歷程（Re=1.0×106）Fig.5 Time history of vertical load（Re=1.0×106）

2.3 渦發放頻率對比

由于泄渦誘發水翼的垂向載荷發生波動，因此，能量譜中的最大譜峰所對應的頻率等同于渦發放的頻率。對垂向載荷進行傅里葉變換，為了達到精度和計算時間上的平衡，時間步長選為0.05T，經過傅里葉變換后，得到渦發放頻率區間如圖6 所示，其中包括預報值上限和預報值下限，兩者的平均值視為渦發放頻率的數值模擬近似值。

圖6 垂向載荷能量譜（Re=0.8×106）Fig.6 Power spectral density of Fy（Re=0.8×106）

圖6～圖9 給出了NACA 66 水翼在4 種不同流速下的垂向載荷的能量譜，渦發放頻率值見表1。將所得頻率值與Leroux 等［31］發表的試驗測量結果進行對比，可見在4 種流速下，測量值均落于預報的下限和上限范圍內。圖10 給出了數值模擬和試驗測量的渦發放頻率對比。由圖中可以看出，離散渦法預報的頻率平均值與測量結果吻合良好，且在較高雷諾數下，模擬結果更為準確，在低雷諾數時，預報值略高于試驗結果。這可能是由于低雷諾數時水翼繞流邊界層處于層流向湍流的過渡階段，本文的數值模型未建立湍流模型，對邊界層發展和分離采用隨機走步進行近似處理，無法實現層流到湍流這一過渡狀態的精確捕捉。

圖7 垂向載荷能量譜（Re=1.0×106）Fig.7 Power spectral density of Fy（Re=1.0×106）

圖8 垂向載荷能量譜（Re=1.2×106）Fig.8 Power spectral density of Fy（Re=1.2×106）

圖9 垂向載荷能量譜（Re=1.5×106）Fig.9 Power spectral density of Fy（Re=1.5×106）

圖10 渦發放譜峰頻率試驗與數值模擬對比Fig.10 Vortex shedding peak frequency（experiment&numerical simulation）

表1 渦發放譜峰頻率對比Table 1 Comparison of vortex shedding peak frequency

通過實測數據與數值模擬數據，均可看出泄渦的頻率隨著雷諾數的增加而增加。為進一步分析流速對渦發放頻率的影響，圖11 對比了3 種流速下在t/T = 30 時刻的渦點分布圖，其中藍色虛線代表Re=1.5×106時的渦發放位置。對比可見，在水翼后緣1.5 個弦長內，Re=1.5×106時釋放的渦對最多，其次為Re=1.2×106時，再次為Re=0.8×106時，這表明流速的增加會加快泄渦的頻率。

圖11 t/T=30 時刻NACA 66 水翼渦發放頻率對比Fig.11 Vortex shedding frequency of NACA 66 series hydrofoil at t/T=30

3 流速與攻角對水翼渦發放頻率的影響

螺旋槳是艦船三大噪聲源之一，對其進行噪聲預報是艦船聲學設計過程中的一個重要環節。國內外學者開展了大量相關研究，Seol等［32-33］采用面元法得到螺旋槳的流場信息并將其代入聲類比方程中，從理論上求解了流場任一點的輻射噪聲。在數值模擬方面，Pantle 等［34］測量了自制螺旋槳在空氣中工作時的噪聲數據，之后對該槳噪聲進行計算，發現計算結果在葉頻和倍葉頻上與試驗數據吻合較好，從而驗證了數值方法的可行性。然而，由于螺旋槳在水下工作時運動復雜，現有的測量和試驗技術無法準確捕捉噪聲信號，尤其是缺乏對于較高流速（如雷諾數Re ＞107），以及有初始攻角情況下螺旋槳不同葉剖面處渦發放頻率及螺旋槳整體噪聲的試驗及數值預報研究。

無論采用理論方法還是數值模擬方法預報螺旋槳流噪聲，都要先求解螺旋槳的湍流脈動，而在數值模擬方法中，湍流脈動的準確性又取決于湍流模擬方法。魯利等［35］以DTMB 4119 螺旋槳為對象，采用雷諾平均N-S 方程、分離渦模擬法（DES）和大渦模擬法（LES）這3 種湍流模擬方法計算湍流脈動，并通過聲學邊界元預報螺旋槳流噪聲，結果發現這3 種湍流模擬方法在一階葉頻處的預報結果相差不大，隨著階數的升高，采用LES 方法預報的結果與理論方法的結果吻合更好。然而LES 方法對網格要求高、計算量較大、計算效率偏低。本文采用的離散渦數值方法擺脫了傳統CFD 方法對網格的依賴性，大大減少了計算量，同時可采用大時間步長，解析度高，十分適合解決高雷諾數下的渦發放和渦致運動問題。

本文采用NACA 66（MOD）水翼模型［36］，水翼弦長為0.72 m，最大厚度為弦長的10%，最大厚度位于距離前緣50%弦長處。NACA 66（MOD）水翼截面型式如圖12 所示。

圖12 NACA 系列水翼截面（最大厚度為弦長的10%）Fig.12 Section of NACA series hydrofoil（maximum thickness is 10% of chord length）

3.1 來流速度對渦發放頻率的影響

應用離散渦法模擬NACA 66（MOD）水翼在不同流速中的渦發放情況，水翼的攻角為10°。通過對垂向載荷進行傅里葉變換，得到了雷諾數與渦發放頻率的關系，如圖13 所示。由圖可見，隨著雷諾數的增加，泄渦頻率不斷增加。然而，在不同的雷諾數區間，泄渦頻率的增加速率不同，兩者呈現非線性關系。在雷諾數為2.5×107～3.5×107區間，泄渦頻率增加緩慢，這表明水翼在這一流速范圍內會穩定發出低頻噪聲，此時泄渦頻率在水翼結構固有頻率（60 Hz）附近變化。

圖13 雷諾數與渦發放頻率關系Fig.13 Relationship between Reynolds number and vortex shedding frequency

對于不同雷諾數，同一時刻（t=0.4×C/U）下水翼繞流的尾流場渦量云圖如圖14 所示。從圖中可以看出，雷諾數越大，新生渦運動至尾流場越快，與水翼后緣的距離越遠，這是由于流體速度的增加會加快邊界層附近流體的運動和分離，新生渦團在來流作用下快速發展并運動至尾流場，因此泄渦頻率變大。

3.2 攻角對渦發放頻率的影響

對上述NACA 66（MOD）水翼在5 種不同初始攻角情況下的泄渦過程進行模擬。攻角分別取為5°，7.5°，10°，12.5°和15°，來流速度約為14.0 m/s，對應的雷諾數為1.0×107。攻角與渦發放頻率的關系如圖15 所示。從圖中可以看出，泄渦頻率隨攻角的增大而減小。當攻角小于10°時，泄渦頻率在30～35 Hz 區間；當攻角大于10°之后，泄渦頻率下降顯著；當攻角為15°時，泄渦頻率已降低到13.5 Hz。

在Re=1.0×107條件下，同一時刻（t=0.8×C/U）不同攻角下水翼繞流的渦強對比如圖16 所示。從圖中可以看出：隨著水翼攻角的增加，泄渦尺寸和渦強不斷增加，水翼縱搖運動的能量耗散越來越多；在相同流速（相同能量輸入）情況下，泄渦變得更加緩慢，因此泄渦頻率變小。

圖14 流速與渦強的關系Fig.14 Relationship between flow velocity and vortex strength

圖15 攻角與渦發放頻率的關系（Re=1.0×107）Fig.15 Relationship between angle of attack and vortex shedding frequency（Re=1.0×107）

圖16 攻角與渦強的關系（Re=1.0×107）Fig.16 Relationship between angle of attack and vortex strength（Re=1.0×107）

對5 個不同流速下，對應的雷諾數為0.5×107，1.0×107，2.0×107，3.0×107和4.0×107，各初始攻角的泄渦過程也同樣進行了模擬，攻角與渦發放頻率的關系如圖17 所示。從圖中可以看出：在同一流速下，泄渦頻率隨著攻角的增大而減??；在某一初始攻角下，隨著流速的增加，泄渦頻率同時增大，與上面得到的結論相符。

圖17 攻角與渦發放頻率的關系Fig.17 Relationship between angle of attack and vortex shedding frequency

4 結 論

本文應用基于渦量—流函數方程的離散渦模型，對高雷諾數下二維剛性水翼渦發放問題進行了數值模擬，得到如下主要結論：

1）離散渦法能夠有效模擬NACA 66 系列水翼渦發放現象，并準確預報高雷諾數下水翼渦發放的頻率。

2）來流速度對水翼泄渦影響顯著，渦發放頻率會隨著來流速度的增大呈非線性增加。

3）隨著水翼初始攻角的增加，水翼釋放的旋渦尺寸和渦強不斷增大，泄渦頻率逐漸降低。