洪鋒，薛環(huán)鋮,張帆，胡濤*

(1. 三峽大學(xué)水電機械設(shè)備設(shè)計與維護湖北省重點試驗室，湖北 宜昌 443002；2. 三峽大學(xué)機械與動力學(xué)院，湖北 宜昌 443002；3.江蘇大學(xué)國家水泵及系統(tǒng)工程技術(shù)研究中心，江蘇 鎮(zhèn)江 212013)

空化是一種復(fù)雜兩相湍流流動，包含了空泡的產(chǎn)生、生長、收縮和潰滅等多種非穩(wěn)態(tài)過程，涉及可壓縮、旋渦運動和多尺度流動分離等非定常特征.空化在流體機械中廣泛存在，嚴(yán)重時會造成機組振動、噪聲及空蝕等問題[1-2].

由于空化具有高湍流特性，湍流模型在準(zhǔn)確模擬空化流中起著重要作用[3].大渦模擬(LES)已被證明能夠精確再現(xiàn)流場的不穩(wěn)定性[4]，并在水翼云空化流動問題中得到了廣泛應(yīng)用[5-6].為研究二維幾何條件下LES模型的適用性，ROOHI等[7]和謝慶墨等[8]分別在OpenFOAM和Fluent平臺，實現(xiàn)了對二維Clark-Y水翼云空化的精確計算，驗證了2D LES方法模擬翼型空化的可行性.

LES引入了一種基于空間平均的濾波方法，對大尺度流場直接求解，并通過亞格子模型對較小的渦流進行?；幚?Smagorinsky- Lilly渦黏性模型為最早的亞格子模型，該模型中的系數(shù)Cs是一個經(jīng)驗值，導(dǎo)致該模型在計算過程中常產(chǎn)生過度的耗散[9].隨著研究的深入，相繼產(chǎn)生了動態(tài)亞格子模型(KET)，以及考慮湍流壁面效應(yīng)的亞格子模型，如WALE模型.不同亞格子模型主要聚焦于單相、簡單流動的研究，如后向臺階流動問題[10].然而，空化流場的劇烈非定常特征，導(dǎo)致了流動結(jié)構(gòu)的復(fù)雜性和多變性，在LES中所選取的亞格子模型需充分捕獲大小尺度湍流結(jié)構(gòu)之間的相互關(guān)系，并最終決定了空化流動數(shù)值計算的精確度.但目前水翼空化流場數(shù)值計算較少關(guān)注大渦模擬方法中亞格子模型的適用性問題.

為了分析不同亞格子模型在空化非定常流動中的計算精度及適用性，文中分別使用4種亞格子模型，對二維Clark-Y水翼周圍云空化流動特性進行模擬，并與相應(yīng)試驗結(jié)果對比研究，從而對上述亞格子模型計算效果以及精確性進行評價.

1 數(shù)值模型

文中采用Mixture多相流模型，假設(shè)多組分具有相同的分速度.空化模型采用Zwart模型，該模型基于Rayleigh-Plesset方程，并假設(shè)系統(tǒng)中氣泡具有相同的尺度.Zwart模型在水翼、離心泵等流體機械空化流動數(shù)值計算中得到了廣泛應(yīng)用[11-12].

LES對Navier-Stokes方程進行濾波，并將湍流分解為可解大尺度和不可解小尺度2個部分.濾波后的Navier-Stokes方程[13]為

(1)

式中：τij為亞格子應(yīng)力；Sij為應(yīng)變率張量.

Smagorinsky-Lilly亞格子模型湍流黏度系數(shù)[14]計算式為

(2)

式中：Ls為亞格子尺度的混合長度；D為翼型表面第一層網(wǎng)格厚度；k為Karman常數(shù)；Cs為Smagorinsky常數(shù)，取0.1;Δ為局部網(wǎng)格尺寸.

在WALE亞格子模型中，湍流黏度系數(shù)[15]為

(3)

(4)

其中Cs取0.325.

WMLES模型會在空間和時間尺度都過小的近壁區(qū)域激活RANS模式，在網(wǎng)格大小足以解析自由流的局部大尺寸時切換到大渦模擬LES，其湍流黏度系數(shù)[15]為

(5)

式中：dw為距壁面的距離；S為應(yīng)變率；取κ= 0.41，Cs= 0.2；hmax為最大網(wǎng)格尺寸；hwn為標(biāo)準(zhǔn)網(wǎng)格間距；取Cw=0.15.

KET亞格子模型增加了對亞格子尺度湍流效應(yīng)的模擬，其湍流黏度系數(shù)[16]定義為

(6)

式中：ksgs為亞格子動能；參數(shù)Ck的值根據(jù)流場狀態(tài)調(diào)整.

2 計算域及網(wǎng)格劃分

2.1 數(shù)值計算模型

研究對象為弦長C=70 mm、翼型攻角α=8°的Clark-Y型水翼.水翼位于矩形流場中，流場域的幾何參數(shù)如圖1所示.計算域邊界條件設(shè)定與ROOHI等[7]的研究保持一致，即流場入口設(shè)置為速度入口邊界，速度Uin=10 m/s；出口設(shè)置為壓力出口邊界，壓力值由空化數(shù)決定；上壁面和下壁面以及水翼表面均設(shè)為無滑移壁面條件.

圖1 計算域Fig.1 Computational domain

2.2 網(wǎng)格無關(guān)性分析

采用WALE亞格子模型計算水翼非定常流動，來驗證網(wǎng)格無關(guān)性.驗證的網(wǎng)格數(shù)分別為80 636(Mesh1)，127 696(Mesh2)，170 066(Mesh3)及253 764(Mesh4).翼型壁面周圍采用O型網(wǎng)格，并對近壁面網(wǎng)格進行加密.文中采用壓力-速度直接耦合法(couple)求解控制方程，擴散項采用二階中心差分格式，對流項采用高階分辨率格式.瞬態(tài)項采用非穩(wěn)態(tài)二階隱式時間積分格式.為滿足湍流模型對壁面網(wǎng)格的要求，y+的值不大于1.根據(jù)Courant-Friedrichs-Lewy condition條件，確定時間步長Δt= 1.0×10-6s.

表1為σ=0.8條件下，采用4種不同網(wǎng)格計算得到的平均升、阻力系數(shù)Cl,Cd.從表中可以發(fā)現(xiàn)，Mesh2計算所得的升阻力系數(shù)與試驗值已較為吻合，且繼續(xù)增大網(wǎng)格數(shù)量，升阻力系數(shù)的變化基本穩(wěn)定，表明網(wǎng)格數(shù)量沒有導(dǎo)致求解誤差.因此，選擇Mesh2作為文中的計算網(wǎng)格.

表1 不同網(wǎng)格數(shù)下的升、阻力系數(shù)對比Tab.1 Comparison of predicted lift and drag coe-fficients based on different grids

3 結(jié)果與討論

3.1 升阻力系數(shù)變化規(guī)律

圖2為4種亞格子模型在一個空化周期內(nèi)升力系數(shù)的瞬時變化，并與試驗值進行對比.

圖2 不同亞格子模型預(yù)測的瞬時升力系數(shù)與試驗值的比較Fig.2 Comparison of transient lift coefficient obser-ved by different sub-grid scale models and experiments in cavitation cycle

由于空泡生長、脫落和潰滅，升力隨時間呈現(xiàn)振蕩變化.空泡長度達到最大時升力系數(shù)最大，而小幅度的振蕩說明云空化階段空泡的小規(guī)模分離.對比發(fā)現(xiàn)，WALE模型預(yù)測的瞬時升力系數(shù)周期性變化更接近試驗數(shù)據(jù)，KET模型及WMLES模型預(yù)測的升力系數(shù)峰值較高，升力系數(shù)變化過程較為劇烈，Smagorinsky-Lilly模型預(yù)測的升力峰值雖然與試驗較為接近，但周期性變化明顯弱于WALE模型的預(yù)測結(jié)果.

表2為不同亞格子模型模擬得到的平均升、阻力系數(shù)與相應(yīng)試驗值的對比.從表中可以看到，文中采用的WALE模型模擬得到的平均升力系數(shù)與試驗值的相對誤差δ在1.00%以內(nèi)，而其余3種模型的計算值及文獻預(yù)測值均大于1.00%；就平均阻力系數(shù)而言，文中WALE與Smagorinsky-Lilly模型預(yù)測的結(jié)果與試驗值的相對誤差均在12.00%以內(nèi)，且遠遠低于ROOHI等[7]的預(yù)測值.綜上，根據(jù)平均升阻力系數(shù)及瞬時升力系數(shù)預(yù)測精度,文中采用的WALE具有較好的適用性與可行性.

表2 4種亞格子模型計算得到的平均升力、阻力系數(shù)對比Tab.2 Comparison of average lift and drag coefficients calculated by four sub-grid scale models

3.2 時均速度分布規(guī)律

圖3為σ=0.8時，沿水翼吸力面垂直方向5個不同位置(x/c=0.2，0.4，0.6，0.8，1.0)，不同亞格子模型預(yù)測的平均速度分量與相應(yīng)試驗值[13]的對比.

圖3 時均速度數(shù)值結(jié)果與試驗數(shù)據(jù)的比較Fig.3 Comparison of numerical results of time-averaged velocity with experimental data

從圖3中可以看到，當(dāng)x/c小于等于0.4時，水翼前緣處速度平均值大于0，而當(dāng)x/c大于等于0.6時，水翼后緣處由于逆壓差的影響，速度出現(xiàn)負值，表明產(chǎn)生了回射流.對比發(fā)現(xiàn)，4種亞格子模型預(yù)測的平均速度與試驗結(jié)果均吻合良好.在空化和非空化區(qū)域以及邊界層內(nèi)部，不同亞格子模型預(yù)測的結(jié)果差異也較小，這是由于4種亞格子模型在壁面及主流區(qū)域預(yù)測到的空化程度比較接近.

3.3 云狀空泡周期性演變

圖4為不同亞格子模型模擬得到的空泡形態(tài)周期性變化和試驗結(jié)果對比.觀察發(fā)現(xiàn)，不同亞格子模型計算得到的空泡演變基本符合試驗結(jié)果.在5%T至20%T，上一個空化周期的云空泡團脫離翼型表面，向下游低壓區(qū)移動并完全潰滅，此時附著型空泡在前緣產(chǎn)生并向后緣生長，空穴的厚度和長度隨之變大.隨著時間推移，空泡幾乎完全覆蓋整個翼型表面.

圖4 不同亞格子模型模擬得到的空泡形態(tài)的周期性變化和試驗值的比較Fig.4 Comparison of periodic changes of cavity shape simulated by different sub-grid scale models and experimental data

在60%T至70%T，空穴從尾緣開始破裂并向前緣推進，回射流進一步向上游移動，導(dǎo)致空穴斷裂；在81%T至97%T，前緣的附著型空泡基本潰滅，尾緣處形成大尺度空泡團.隨時間推移，前緣開始產(chǎn)生較薄的附著型空泡，尾緣大體積空泡團脫落并向下游運動.

比較4個模型模擬結(jié)果可以發(fā)現(xiàn)，WALE模型模擬得到的空泡體積分?jǐn)?shù)云圖與試驗結(jié)果更接近.具體體現(xiàn)在：① 根據(jù)試驗描述，在10%T時刻，翼型尾緣處大體積云空泡團還在翼型表面，WALE模型的結(jié)果與試驗相符，其他3個模型的云空泡團已脫離尾翼表面；② 在20%T時刻，試驗結(jié)果中，翼型表面上的片狀空泡長度與厚度達到最大，幾乎覆蓋整個吸力側(cè)表面，WALE模型較好地滿足這一特性，而KET及WMLES模型所預(yù)測片狀空泡厚度及長度都過低；③ 根據(jù)文獻[13] Clark-Y水翼試驗結(jié)果得出，云空化的不穩(wěn)定性主要歸結(jié)于回射流的作用.在60%T時刻，回射流開始作用于尾緣處并向前緣推進，空腔閉合處的反向壓力梯度有利于回射流的發(fā)展.KET與WMLES模型其后緣處有小規(guī)?？张輸嗔?，其回射流作用效果過大，相較之下，Smagorinsky-Lilly與WALE模型更符合試驗結(jié)果;④ 片狀空化在翼型中部斷裂，其空泡長度減少至上一時刻的2/3左右，WALE模型的結(jié)果與試驗數(shù)據(jù)較吻合，尤其是翼型中后部空穴特性.此外，圖5為20%T時刻，水翼尾緣處的速度矢量分布圖，對比發(fā)現(xiàn)WALE模型尾緣處空泡受到的回射流影響很小，逆壓力梯度作用不明顯，而Smagorinsky-Lilly模型此時受到回射流影響，片狀空穴尾端發(fā)生明顯回縮，回流長度增大，表明除WALE模型以外的亞格子模型過早地預(yù)測到翼型尾緣的回射流結(jié)構(gòu).因此，總體上，WALE模型能較好地預(yù)測云狀空化周期性演變過程及其流動細節(jié).

圖5 20%T時水翼尾部速度矢量分布圖Fig.5 Velocity vector distribution diagram of hydro-foil tail at 20%T

為進一步理解空化流動機理，研究空化與旋渦相互作用.利用Q準(zhǔn)則對非定常渦結(jié)構(gòu)進行可視化，Q準(zhǔn)則定義為

(7)

式中：ω為旋轉(zhuǎn)率張量；s為應(yīng)變率張量；Q準(zhǔn)則表明旋轉(zhuǎn)與變形之差，其值表示旋轉(zhuǎn)強度的大小，用正值來識別渦流，負值則說明在該區(qū)域作為主導(dǎo)的為剪切作用.

圖6為利用WALE模型預(yù)測的一個空化發(fā)展內(nèi)某一時刻Q值分布情況.觀察得到，5%T時刻，大尺度云狀空泡團運動到水翼尾緣，尾緣附近Q值分布較大，表明此處存在較強的旋轉(zhuǎn)和剪切作用；20%T時刻，片狀附著性空穴向后緣生長，在空穴內(nèi)部的Q值較為穩(wěn)定，表明在空穴內(nèi)部流體的旋轉(zhuǎn)效應(yīng)和剪切效應(yīng)相對平衡，沒有表現(xiàn)出明顯的渦流；70%T時刻，在回射流的作用下，片狀空穴發(fā)生斷裂，且斷裂空泡周圍Q值為負，說明回射流的軸向變形效應(yīng)較強；90%T時刻大尺度空泡團生成，空泡團邊界及尾緣附近的Q值為正且絕對值極大，說明這部分區(qū)域有較強的旋渦作用，旋轉(zhuǎn)與剪切效應(yīng)相互制衡和轉(zhuǎn)化，致使該區(qū)域流動出現(xiàn)不穩(wěn)定.

圖6 WALE模型計算得到的Q值分布Fig.6 Q value distribution calculated by WALE model

4 結(jié) 論

1) 與其他3種亞格子模型相比，WALE模型預(yù)測的平均升阻力系數(shù)與試驗結(jié)果最為符合，且在1個典型空化周期內(nèi)，WALE模型預(yù)測的瞬時升力變化趨勢與試驗也最接近.

2) 4種亞格子模型預(yù)測的空泡演變過程基本與試驗相符，但WALE模型模擬得到的空泡在云空化周期內(nèi)的形態(tài)、大小與位置更符合試驗結(jié)果，能夠比較準(zhǔn)確地捕捉到二維Clark-Y水翼云空化時翼型尾部回射流作用所導(dǎo)致的空泡斷裂效果，并精確預(yù)測到大型空泡團的脫落及其體積變化特征.

3) 使用Q準(zhǔn)則，描述了基于WALE模型下翼型空化流動中的旋渦結(jié)構(gòu)，得到了片狀空化脫落所產(chǎn)生的強烈旋轉(zhuǎn)與剪切作用對流場不穩(wěn)定性的影響.