某二維水翼空化數(shù)值模型對比研究

蔣文青+鄭源+葛新峰

摘要：以二維水翼計算模型為研究對象，采用數(shù)值模擬的方法，通過不同的空化模型和湍流模型對計算模型進行數(shù)值模擬研究，得到了翼型在空化過程中空泡的運動規(guī)律。結(jié)果表明：該翼型在不同空化模型下，發(fā)生空化時，空化過程包含了空泡的產(chǎn)生、發(fā)展、脫落和潰滅幾個過程，空化的運動過程總具有明顯的周期性；Schnerr-Sauer模型的計算結(jié)果與實際的水翼空化流動過程相吻合，Schnerr-Sauer空化模型可以很好地模擬繞水翼的空化流動過程。

關(guān)鍵詞：Clark-Y翼型；空化；數(shù)值模擬；空泡；運動規(guī)律

中圖分類號：TK730.2 文獻標(biāo)識碼：A 文章編號：1672-1683（2017）03-0204-05

Abstract：We took a two-dimensional hydrofoil calculation model as the research object，and adopted the method of numerical simulation.The movement of the hydrofoil cavity in the cavitation process was obtained through the numerical simulation study on different cavitation models and turbulence models.The result showed：In different types of cavitation models，the cavitation process of this hydrofoil included the generation，development，fall-off，and burst of vacuoles，and the process showed obvious periodicity.The simulation result of S-S cavitation model was consistent with the actual movement.The S-S cavitation model can well simulate the flow around a hydrofoil in the cavitation process.

Key words：Clark-Y hydrofoil；cavitation；numerical simulation；vacuoles；movement pattern

水力機械轉(zhuǎn)輪作為重要的能量轉(zhuǎn)換裝置，其葉片斷面形狀主要為翼型，所以對翼型的研究和設(shè)計是提高水力機械性能的主要途徑之一。

Huang等[1]應(yīng)用完全空化模型對某NACA66（MOD）水翼在不同的攻角時的定?？栈螒B(tài)進行了研究。褚學(xué)森等[2]應(yīng)用完全空化模型研究了回轉(zhuǎn)體和圓盤空化器的自然空化流動，給出了不同氣核質(zhì)量分?jǐn)?shù)和湍流參數(shù)對空化流動的影響結(jié)果。郝宗瑞等[3]針對翼型攻角為8°的二維NACA0015水翼展開了與流體流動特性和流場結(jié)構(gòu)相關(guān)的研究。Kang等[4]數(shù)值模擬了直徑為200 mm的四葉軸流泵的三維空化流動。劉登成等[5]對直徑為227 mm的INSEAN E779A四葉模型槳的空化流動進行了研究。張博等[6]為了可以更直觀的觀察出空化流動的結(jié)果，利用實驗的方式對繞水翼空化流動展開了研究。張博等[7]應(yīng)用商用軟件CFX的Kubata 空化模型，描述了非定常云狀空化形態(tài)的發(fā)展過程。雖然目前已有較多的使用這些模型來研究非定?？栈奈墨I，但是，由于在一種軟件中僅僅會包含一種模型，缺乏對多種模型的對比，特別是對空化流場的結(jié)構(gòu)模擬的對比研究還不是很充分。

本文利用CFD數(shù)值模擬方法對翼型的外部繞流形態(tài)進行了研究。利用FLUENT[8]對翼型進行不同空化模型以及在不同湍流模型下的空化模擬計算，得到翼型空化的規(guī)律。

1 數(shù)學(xué)模型和計算方法

1.1 空化模型

Schnerr-Sauer 這一空化模型是基于氣泡動力學(xué)和汽核半徑共同得到的，式中的汽核半徑RB與體積分?jǐn)?shù)α和汽核數(shù)密度 n的關(guān)系為RB=α1-α34πn1/3，其中n為給定的經(jīng)驗常數(shù)，默認(rèn)值為n=1×1013。

此空化模型假設(shè)空化氣泡不相互影響，并且隨著蒸汽體積分?jǐn)?shù)的增加，氣核密度相應(yīng)的減小。式中，氣核半徑RB取值為10-6m，形核位置體積分?jǐn)?shù)αnuc取值為5×10-4，蒸發(fā)系數(shù)Fvap取為50，冷凝系數(shù)Fcond取為0.001。

1.2 湍流模型

湍流的數(shù)值模擬方法主要采用了非直接模擬方法中的大渦模擬方法和雷諾平均方法中的RNG k-ε模型。大渦模擬[9-15]方法的基本思想大致上可以概括為：對于大尺寸的渦，采用瞬時的Navier-Stokes方程來直接模擬，并不直接模擬小尺寸的渦，對于小尺寸的渦，則需要通過建立模型考慮其影響。大渦模擬方法對計算機的要求仍然是比較高的，依舊是計算機的內(nèi)存以及計算速度方面的要求，不過，相比直接數(shù)值模擬方法而言，要求是有所降低的，因此，在目前的CFD研究和應(yīng)用中成為了一個熱點。

RNG k-ε模型被稱為重整化群，是由Yakhot及Orzag[16]在研究的過程中提出來的。在這一計算模型中，將尺寸比較大的空泡運動的影響，通過修正的方法，從而將黏度的影響體現(xiàn)在小尺寸的空泡上，這樣一來就可以將小尺寸的影響效果從整體的運動系統(tǒng)中有效的剔除出去，最終才能夠得到了計算需要的k方程和ε方程。RNG k-ε模型在對于修正湍動黏度的過程中，一般需要將平均流動中所產(chǎn)生的漩渦和旋流等流動狀況列入考慮的范圍；RNG k-ε模型中產(chǎn)生項與流體的流動有一定的相關(guān)性，不過在解決問題的過程中還可以看作是空間坐標(biāo)的函數(shù)的問題；RNG k-ε模型的優(yōu)點在于在面對應(yīng)變率較高的狀態(tài)或者流線的彎曲程度很大的情況時，模型的處理效果與實際的發(fā)展情況是很接近的。

1.3 計算區(qū)域和網(wǎng)格

本文的研究對象選取了Clark-Y翼型，翼型為非對稱的結(jié)構(gòu)，該翼型為上彎下平，其中最凸出處為弦長20%處，機翼的厚度為弦長的11.7%，弦長L為0.15m，攻角為零度，整個計算區(qū)域見圖1，上游來流進口距翼型頭部距離為4.5 L，下游去流出口距翼型尾部距離為8 L，兩側(cè)壁面距翼型中心各為8 L。在ICEM CFD中將Clark-Y翼型模型導(dǎo)入，對該模型進行網(wǎng)格劃分，最終得到適應(yīng)性強的結(jié)構(gòu)化網(wǎng)格。在翼型的表面處流體的流動并不是很規(guī)律的，流動變化較快，因此一般會對翼型的特殊區(qū)域進行了網(wǎng)格加密處理，這樣一來便可以較為真實的模擬復(fù)雜的流動狀態(tài)。對于重要的區(qū)域進行了網(wǎng)格加密，而對于普通的流動區(qū)域便可以適當(dāng)?shù)臏p少網(wǎng)格的數(shù)目，防止網(wǎng)格數(shù)太多增加計算的時間和計算量。并且經(jīng)過網(wǎng)格無關(guān)性驗證之后，最終確定的網(wǎng)格節(jié)點總數(shù)為69 969，并確定了壁面附近的y+值為30～60之間，見圖2。

圖1 計算區(qū)域

Fig.1 The computational domain

圖2 計算域和網(wǎng)格示意圖

Fig.2 Sketch of the computational domain and grids

1.4 邊界條件的設(shè)置

本文對翼型模型進行空化的數(shù)值模擬研究，并通過研究結(jié)果對空化模型進行比較和分析。選定入口邊界條件為速度入口，假設(shè)流動的過程為不可壓縮流體流動，選定入口速度為10 m/s，方向為水平方向；出口邊界條件設(shè)置為壓力出口；壁面設(shè)置為無滑移壁面剪切條件，即貼在壁面的流體的運動狀況與該壁面相同；湍流模型采用與實際情況較為吻合的RNG k-ε模型；進口與回流的湍流邊界條件均給定湍流強度與特征尺寸；求解方法采用壓力基求解器，SIMPLEC壓速耦合算法；選擇的多相流模型為和混合物（Mixture）模型，本文描述了多相流動各相之間的相互作用和相互轉(zhuǎn)化情況，具體說來就是在連續(xù)的液相流動中貫穿有汽相組分。

2 計算結(jié)果與分析

2.1 數(shù)值計算結(jié)果

通過本次空化數(shù)值模擬，最終得到了在兩種不同的空化模型以及湍流模型，不同時間下的空泡體積分?jǐn)?shù)。圖3和圖4分別表示了在RNG k-ε模型下，圖5和圖6分別表示了在大渦模擬方法下，不同空化模型以及不同時刻的水蒸氣含量隨時間變化的分布云圖。

在圖3所示的Schnerr-Sauer空化模型下，較為明顯的出現(xiàn)了氣泡的成長、附著、脫落和潰滅的過程。開始時，水翼的前端會產(chǎn)生較為穩(wěn)定的空化，不過隨著時間的發(fā)展，從第1 ms開始，空化的氣泡在不斷長大的過程中也在往翼型的尾部移動，整個過程發(fā)生的非常明顯，并且因為來流具有一定的合適的速度，所以空泡在第13 ms的時候便開始發(fā)生斷裂的情況，隨著空化的進行，空化氣泡便逐漸的縮小，到達了下一個周期的發(fā)展與潰滅的過程。在15 ms的時候便進入到了第二個周期的空泡的分離與潰滅的過程。在整個周期中，發(fā)生的時間并不是很長，在十幾毫秒的時候，整個空化便差不多已經(jīng)完全了，在這其中空泡的成長占據(jù)了較多的時間。

圖4所示為在Z-G-B空化模型下的空化結(jié)果，與圖3中的結(jié)果相比較可知，圖3中的空化結(jié)果較好，無論是成長、分離還是潰滅的過程都比較明顯，雖然，圖4空泡的數(shù)量較多，不過整體的變化過程并不是很好，空泡長大的過程中，生長速度較為緩慢，但附著在翼型表面的空泡較為穩(wěn)定；空泡并沒有一個明顯的分離的過程，并且在縮小的過程中空泡的體積分布圖有著明顯的變化，在接近翼型的尾部處，空泡的體積有一段突變的狀態(tài)，導(dǎo)致整體的形態(tài)不好。從第1 ms開始，便從空化的前端發(fā)生了空化，并且前端的空泡越來越多，發(fā)生的空化現(xiàn)象很嚴(yán)重，并且會不斷地向前方進行擴散，與此同時，在翼型后側(cè)產(chǎn)生的氣泡也在不斷的向后部移動。從5 ms開始到10.5 ms這段時間內(nèi)，主要是空泡長大的過程，這一過程在這一時間段內(nèi)的變化并不是很明顯，在10.5 ms到32 ms的時間段內(nèi)，空泡便在發(fā)生潰滅，這一階段在整個的時間段中占據(jù)了很長一部分的時間，不過在潰滅的過程中，附著在翼型的壁面上的空泡體積分?jǐn)?shù)較大的一部分的形狀并不是很規(guī)整。與Schnerr-Sauer空化模型計算結(jié)果相比較而言，Z-G-B空化模型的計算結(jié)果顯示空化的形態(tài)并不是很好，高壓區(qū)總是附著在翼型的表面，并且持續(xù)的時間長，因此空化的結(jié)果并不好。

圖5 大渦模擬Schnerr-Sauer空化模型

Fig.5 Large eddy simulation of Schnerr-Sauer cavitation model

Schnerr-Sauer空化模型計算得出的結(jié)果較好（圖5），從產(chǎn)生、成長、飄移到潰滅的過程都是合理的，在翼型的中部往后開始產(chǎn)生空泡，隨著時間的發(fā)展，附著在翼型表面的空泡的數(shù)量越來越多，這也是空化產(chǎn)生的根本原因。于此同時，空泡的產(chǎn)生位置也在開始往后移動。從0.1 ms開始，空泡便在不斷地形成過程中，在6 ms到8 ms之間空泡由于來流具有一定的合適的流速，所以空泡有著向后移動的趨勢，并且表現(xiàn)的結(jié)果較為顯著，在6.5 ms到8 ms之間空泡是發(fā)生了斷裂的過程的，然后在8 ms時刻，空泡在不斷的潰滅、縮小，在8 ms到9 ms之間空泡又進入到了一個新的周期的循環(huán)之中，不過這一循環(huán)與上一循環(huán)相比較，整個過程不僅空泡的尺寸減小了，發(fā)生空化的時間周期也明顯的縮短了，綜上所述，在圖5中，較完整地展現(xiàn)了出了一個完整的周期以及下一個即將開始的周期，整體來說，應(yīng)用大渦模擬這一算法對于Schnerr-Sauer空化模型的計算結(jié)果還是比較可靠的，計算的結(jié)果也比較穩(wěn)定。圖5的計算結(jié)果與圖3相比較之后可知，兩者之間的空化時間是比較有些差距的，之間總共差了大約有5 ms的時間，而空泡在第一階段發(fā)生分離的時間都是在接近空化結(jié)束時候發(fā)生的，所以，在計算模型發(fā)生變化時，空化的時間是有一定的變化的，這一點在以后的研究中需要多加注意。

Z-G-B空化模型的計算結(jié)果見圖6，雖然產(chǎn)生了較多的空泡，在翼型的表面也附著了較多的氣泡，在翼型前端的空泡較為穩(wěn)定，不過隨著時間的發(fā)展，附著在翼型表面的氣泡并沒有向翼型尾部移動。在空化發(fā)生之后的潰滅過程中，效果也不是很明顯，潰滅的過程中，附著在翼型表面的空泡體積含量不斷地增加了。從第1 ms開始，空泡接近翼型的前端的位置開始發(fā)生空化現(xiàn)象，在從4 ms開始的時刻，空泡明顯的開始向翼型的尾部一側(cè)平移，但并沒有發(fā)生空泡斷裂的狀況，可能是由于計算時的速度沒有足夠的大，使得空化產(chǎn)生的氣泡并沒有能夠與翼型的壁面分離，不過隨著時間的發(fā)展，附著在翼型表面的空泡的體積分?jǐn)?shù)越來越大。將圖6與圖4對比之后發(fā)現(xiàn)，在大渦模擬的湍流計算方程下，空化發(fā)生的時間明顯變得更加短了，時間變?yōu)榱薘NG k-ε湍流計算方程計算結(jié)果的一半左右，而空化的結(jié)果是相似的。

綜上所述，在大渦模擬湍流模型的計算結(jié)果之下，空化發(fā)生的時間明顯的縮短了，但是，空化的形態(tài)并沒有發(fā)生太大的改變，整體的空化過程還是一致的，包括了空泡的產(chǎn)生、成長、分離以及潰滅，不過對于Z-G-B空化模型來說，分離和潰滅并不是一個很明顯的過程，所以，在空化數(shù)值模擬時，選擇的空化模型不一致，得到的計算結(jié)果從整體上而言是沒有太大的差別的。不過需要注意的是，空化的時間周期有了較為明顯的縮短，但是相對于整體的時間周期而言，空化發(fā)生的相對時間沒有太大的變化。

2.2 實驗對比

王國玉團隊設(shè)置的邊界條件分別是：10 m/s的速度入口，出口采用的是壓力出口，翼型表面采用的是絕熱、無滑移固壁條件[7]，與本文的設(shè)置相類似。而實驗以及數(shù)值計算的結(jié)果都只是基于修正后的RNG k-ε湍流方程，所以，在本文中就只是將兩種空化模型在湍流方程RNG k-ε下的計算結(jié)果分別于王國玉團隊的實驗以及模擬結(jié)果進行對比，并分析了產(chǎn)生結(jié)果異同的原因。

2.2.1 Schnerr-Sauer空化模型對比

如圖7所示，將Schnerr-Sauer空化模型在RNG k-ε湍流模型下計算得到的結(jié)果與王國玉團隊的繞Clark-y水翼云狀空化流動的數(shù)值計算與實驗對比得到的結(jié)果進行比較，計算結(jié)果與實驗結(jié)果的云圖相似。在空化發(fā)生的時間上也進行了對比，實驗中的空化發(fā)生時間明顯多于數(shù)值模擬的時間。由于攻角的原因，空泡產(chǎn)生的位置兩者也是有區(qū)別的，在實驗中，空化氣泡從翼型前側(cè)邊緣便開始發(fā)生空化了，隨著時間增加，空泡向翼型尾部迅速移動，脫落現(xiàn)象比較明顯。由此可知，攻角的存在會影響空化的起始位置以及空化過程中空泡的脫落程度。

2.2.2 Z-G-B空化模型對比

將Z-G-B空化模型在RNG k-ε湍流模型下的模擬計算結(jié)果與繞水翼云狀空化流動結(jié)構(gòu)的數(shù)值與實驗研究得到的結(jié)果進行比較，得到的結(jié)果差別相對較大。除了上文提到的攻角的影響，使得數(shù)值模擬中得到的空泡基本附著在了翼型的中間部位，還有在Z-G-B空化模型下的計算結(jié)果即所產(chǎn)生的空化形態(tài)與實驗中給出的空化形態(tài)隨時間的變化區(qū)別較大，所以，選用的空化模型與王國玉團隊的選擇可能并不一致，才導(dǎo)致了模擬計算的結(jié)果不一致。除此之外，實驗中的湍流模型是經(jīng)過修正的，這樣可以改善空化流動過程中的計算精度，并且引入了一種與空化區(qū)域水汽相密度相關(guān)的系數(shù)，對渦黏模型進行了修正。實驗的空化過程大約進行了40 ms左右，而數(shù)值模擬所經(jīng)歷的時間為30 ms多，時間上兩者較接近。綜上所述，最終王國玉團隊與Z-G-B空化模型的數(shù)值模擬方法得到的結(jié)果存在一些差距。

3 結(jié)論

通過對二維Clark-Y翼型的非穩(wěn)態(tài)空化數(shù)值計算，研究得到的成果如下。

（1）選擇不同的空化模型和湍流模型，對翼型進行數(shù)值模擬計算，得到了翼型在空化過程中空泡的運動規(guī)律。分別設(shè)置了RNG k-ε和LES兩種湍流模型，并在Z-G-B和Schnerr-Sauer兩種空化模型下，分別進行數(shù)值模擬，最終通過空泡的體積分布云圖，得到了翼型在空化過程空泡的運動規(guī)律：對于任意的空化模型，水翼空化的非定常過程總是包括空泡的產(chǎn)生、發(fā)展、脫落和潰滅幾個過程，并且空化的運動過程具有明顯的周期性。

（2）提出了Schnerr-Sauer空化模型可以很好地模擬繞水翼的空化流動過程。選取目前較為常用的Z-G-B和Schnerr-Sauer空化模型，在相同的工況下進行空化計算，得出Schnerr-Sauer空化模型的計算結(jié)果與王國玉團隊的實驗研究結(jié)果接近，由此可知，對于一般的數(shù)值模擬可以直接采用Schnerr-Sauer空化模型。

參考文獻（References）：

[1] Huang S，He M，Wang C，et al.Simulation of cavitating flow around a 2-D hydrofoil[J].Journal of Marine Science and Application，2010，9（1）：63-68.

[2] 褚學(xué)森，王志，顏開.自然空化流動數(shù)值模擬中參數(shù)取值影響的研究[J].船舶力學(xué)，2007，11（1）：32-39.（CHU Xue-sen，WANG Zhi，YAN Kai.Parametric study on numerical simulation of natural cavitation flow[J].Journal of Ship Mechanics，2007，11（1）：32-39.（in Chinese））

[3] 郝宗瑞，王樂勤，吳大轉(zhuǎn).水翼非定?？栈鲌龅臄?shù)值模擬[J].浙江大學(xué)學(xué)報：工學(xué)版，2010，44（5）：1043-1048.（HAO Zong-rui，WANG Le-qin，WU Da-zhuan.Numerical simulation of unsteady cavitating flow on hydrofoil[J].Journal of Zhejiang University （Engineering Science），2010，44（5）：1043-1048.（in Chinese））

[4] Kang C，Yang M G，Wu G Y，et al.Cavitation analysis near blade leading edge of an axial-flow pump[C].International Conference on Measuring Technology and Mechatronics Automation.Los Alamitos：IEEE Computer Society Press，2009：767-770.

[5] 劉登成，洪方文，張志榮等.螺旋槳片狀空泡CFD的分析[J].艦船科學(xué)技術(shù)，2009，24（6）：740-746.（LIU Deng-cheng，HONG Fang-wen，ZHANG Zhi-rong，et al.The CFD analysis of propeller sheet cavitation[J].Ship Science and Technology，2009，24（6）：740-746.（in Chinese））

[6] 張博，王國玉，黃彪，等.繞水翼空化非定常動力特性的時頻分析[J].實驗流體力學(xué)，2009，23（3）：44-49.（ZHANG Bo，WANG Guo-yu，HUANG Biao，et al，Unsteady dynamics of cavitating flows around a hydrofoil[J].Journal of Experiments in Fluid Mechanics，2009，23（3）：44-49.（in Chinese））

[7] 張博，王國玉等.修正的RNG k-ε模型在云狀空化流動計算中的應(yīng)用評價[J].北京理工大學(xué)學(xué)報，2008，28（12）：1065-1069.（ZHANG Bo，WANG Guo-yu，ZHANG Shu-li，et al.Evaluation of a modified RNG k-ε model for computations of cloud cavitating flows[J].Transactions of Beijing Institute of Technology，2008，28（12）：1065-1069.（in Chinese））

[8] 江帆，黃鵬.Fluent高級應(yīng)用與實例分析[M].北京：清華大學(xué)出版社，2010.（JIANG Fan，HUANG Peng.Fluent：advanced application and case analysis[M].Beijing：Tsinghua University Press.（in Chinese））

[9] Feiz A A，Ould-Rouis M，Lauriat G.Large eddy simulation of turbulent flow in a rotating pipe[J].International Journal of Heat and Fluid Flow，2003，24（3）：412-420.

[10] Mary I，Sagaut P.Large eddy simulation of flow around an airfoil near stall[J].AIAA Journal，2002，40（6）：1139-1145.

[11] Grigoriadis D G E，Bartzis J G，Goulas A.Efficient treatment of complex geometries for large eddy simulations of turbulent flow[J].Computers and Fluids，2004，3（2）：201-222.

[12] Shen L，Yue D K P.Large-eddy simulation of free-surface turbulence[J].Journal of Fluid Mechanics，2001，440：75-116.

[13] Julian R E，Smolarkiewiez K.Eddy resolving simulations of turbulent solar convection[J].International Journal for Numerical Methods in Fluids，2002，39（9）：855-864.

[14] Li J C.Large eddy simulation of complex turbulent flows：physical aspects and research trends[J].Acta Mechanica Sinica，2001，17 （4）：289-301.

[15] Sotiropoulos F.Progress in modeling 3D shear flows using RANS equations and advanced turbulence models[J].Calculation of Complex Turbulent Flows，WIT Press，2001，31-53.

[16] 王福軍.計算流體動力學(xué)分析-CFD軟件原理與應(yīng)用[M].北京：清華大學(xué)出版社，2004.（WANG Fu-jun.Computational fluid dynamics analysis：CFD software principles and applications[M].Beijing：Tsinghua University Press，2004.（in Chinese））