螺旋槳梢渦空化數(shù)值模擬的網(wǎng)格策略研究

張 彬，張志榮

（中國(guó)船舶科學(xué)研究中心，江蘇 無錫 214082）

0 引 言

螺旋槳梢渦是流體在壓差作用下從槳葉壓力面流向吸力面，經(jīng)過梢部時(shí)發(fā)生流動(dòng)分離而產(chǎn)生的旋渦。一般認(rèn)為當(dāng)螺旋槳梢渦渦心壓力低于臨界壓力時(shí)將出現(xiàn)梢渦空化，嚴(yán)重影響艦船振動(dòng)、噪聲等性能，因此工程上對(duì)于梢渦空化的準(zhǔn)確預(yù)報(bào)有著迫切需求。

研究者們通過PIV和LDV技術(shù)對(duì)螺旋槳梢渦空化進(jìn)行了大量試驗(yàn)觀測(cè)和流場(chǎng)測(cè)量,為梢渦空化研究提供了豐富的試驗(yàn)數(shù)據(jù)[1-3]。隨著CFD 的發(fā)展和計(jì)算能力的提高，數(shù)值模擬在梢渦空化研究中發(fā)揮著越來越重要的作用。但是由于梢渦渦心的尺寸小且壓力梯度大，目前較難精確地模擬梢渦空化。提高螺旋槳梢渦空化數(shù)值模擬的精度有許多方法，例如在渦心中設(shè)置較多的網(wǎng)格節(jié)點(diǎn)，采用合適的湍流處理方法和空化模型等。

Chow 等[4]通過計(jì)算驗(yàn)證認(rèn)為，梢渦數(shù)值模擬中至少需要15～20 個(gè)網(wǎng)格點(diǎn)覆蓋梢渦渦心。Tuomas等[5]分別使用標(biāo)準(zhǔn)k-ε湍流模型和k-ωSST湍流模型，基于三套網(wǎng)格對(duì)PPTC螺旋槳進(jìn)行了梢渦空化數(shù)值模擬，發(fā)現(xiàn)網(wǎng)格分辨率對(duì)梢渦流場(chǎng)數(shù)值模擬精度的影響顯著大于湍流模型的影響；劉芳遠(yuǎn)等[6]使用RNGk-ε湍流模型和Zwart-Gerber-Belamri 空化模型對(duì)PPTC 槳進(jìn)行了梢渦數(shù)值模擬，通過局部加密梢渦脫泄區(qū)域的網(wǎng)格，成功模擬了梢渦空化，推力、扭矩系數(shù)與試驗(yàn)接近；胡建等[7]采用大渦模擬（large eddy simulations,LES）方法，基于螺旋加密網(wǎng)格對(duì)E779A螺旋槳進(jìn)行了梢渦空化數(shù)值模擬，展示了精細(xì)的梢渦流場(chǎng)，梢渦空化形態(tài)與試驗(yàn)結(jié)果接近，但梢渦空化區(qū)域較小；Lloyd等[8]分別使用RANS方法和DDES方法對(duì)E779A螺旋槳進(jìn)行了梢渦空化數(shù)值模擬，對(duì)其中一個(gè)槳葉基于無量綱Q準(zhǔn)則判據(jù)進(jìn)行了識(shí)別并對(duì)梢渦區(qū)域進(jìn)行了加密，發(fā)現(xiàn)結(jié)合延遲分離渦模擬DDES方法和自適應(yīng)網(wǎng)格加密方法能夠明顯降低梢渦渦心的最低壓力系數(shù)，梢渦空化模擬結(jié)果與試驗(yàn)更吻合；Yilmaz等[9]使用STAR-CCM+軟件，采用LES方法和Schnerr-Sauer空化模型對(duì)E779A螺旋槳進(jìn)行了空化數(shù)值模擬，并在螺旋加密初始網(wǎng)格的基礎(chǔ)上，通過自適應(yīng)網(wǎng)格加密方法對(duì)絕對(duì)壓力接近飽和蒸氣壓的網(wǎng)格進(jìn)行了加密，計(jì)算所得梢渦空化形態(tài)與試驗(yàn)結(jié)果比較接近。分析比較文獻(xiàn)[7～9]可以得出這一結(jié)論，即不同的湍流模擬方法（RANS、DDES和LES）對(duì)螺旋槳梢渦空化數(shù)值模擬精度有較大影響，但小于網(wǎng)格分辨率的影響。

本文在國(guó)內(nèi)外研究的基礎(chǔ)上，以E779A 槳為研究對(duì)象，基于開源CFD 軟件OpenFOAM 的inter-PhaseChangeFoam 求解器，通過局部網(wǎng)格加密和自適應(yīng)網(wǎng)格加密方法，精確加密梢渦所處的網(wǎng)格，以相對(duì)較少的計(jì)算量準(zhǔn)確模擬近尾流場(chǎng)的梢渦空化。

1 數(shù)值模擬方法

1.1 控制方程

對(duì)于空化流求解，基于Volume of Fluid（VOF）方法將兩相流可看作以汽液兩相按一定比例混合而成的單相流模型?；旌舷嗟拿芏圈裮與動(dòng)力粘性系數(shù)μm可以通過如下公式得到：

式中，αl為液相體積分?jǐn)?shù)，αv為汽相體積分?jǐn)?shù)，ρl為液相的密度，ρv為汽相的密度，μl為液相的動(dòng)力粘性系數(shù)，μv為汽相的動(dòng)力粘性系數(shù)。

假定各項(xiàng)流體不可壓縮，不考慮流體的熱傳遞效應(yīng)和體積力，則連續(xù)性方程和動(dòng)量方程為

由于直接數(shù)值模擬計(jì)算量過大，引入RANS方法得到時(shí)均流動(dòng)的連續(xù)性方程和動(dòng)量方程為

式中，ui和uj為未經(jīng)雷諾平均的速度矢量，和為雷諾平均后的速度矢量，p為未經(jīng)雷諾平均的壓力，pˉ為雷諾平均后的壓力。式（7）中最后一項(xiàng)為雷諾應(yīng)力項(xiàng)，通過引入k-ωSST湍流模型求解。

為了求解體積分?jǐn)?shù)，引入Schnerr-Sauer 空化模型[10],由式（3）可知只需要求解其中一項(xiàng)的體積分?jǐn)?shù)。汽相體積分?jǐn)?shù)輸運(yùn)方程為

聯(lián)立式（1）、（6）和（8）可得

Schnerr-Sauer空化模型基于R-P方程和球形氣泡半徑公式推導(dǎo)得到?αv/?t，最終冷凝率和汽化率的表達(dá)式為

式中，Cc和Cv均為經(jīng)驗(yàn)系數(shù)，本文都設(shè)置為默認(rèn)值1。pv為液體的飽和蒸氣壓,本文設(shè)置為2350 Pa。氣泡半徑R滿足：

式中，n0為單位液體內(nèi)的氣泡數(shù)，設(shè)置為1.6×1013。

1.2 幾何模型

本文研究對(duì)象為INSEAN E779A模型槳，擁有豐富的螺旋槳試驗(yàn)研究數(shù)據(jù)。E779A模型槳為四葉槳，直徑D為0.227 27 m。主要參數(shù)如表1所示。坐標(biāo)系原點(diǎn)設(shè)置在槳中心點(diǎn)，各軸方向如圖1所示。

表1 E779A模型槳主要參數(shù)Tab.1 Particulars of the propeller

1.3 計(jì)算域與工況設(shè)置

計(jì)算域?yàn)?.942D的圓柱體，計(jì)算域左側(cè)為速度入口邊界，到槳中心的距離為1.25D；右側(cè)為壓力出口邊界，到槳中心的距離為4D。在螺旋槳附近區(qū)域設(shè)置圓柱體旋轉(zhuǎn)域，在旋轉(zhuǎn)域內(nèi)使用滑移網(wǎng)格模擬螺旋槳的旋轉(zhuǎn)，其它區(qū)域網(wǎng)格靜止，最后在區(qū)域交界面AMI上進(jìn)行流場(chǎng)數(shù)據(jù)交互。旋轉(zhuǎn)域直徑為1.32D，其左端面和右端面到槳中心的距離分別為0.25D和0.65D。計(jì)算域示意圖如圖2所示。

為了與試驗(yàn)結(jié)果[11]比較，數(shù)值模擬在進(jìn)速系數(shù)J=0.71 和空化數(shù)σ=1.763 工況下進(jìn)行，來流速度Uinf=5.808 m/s，轉(zhuǎn)速n=36r/s，時(shí)間步長(zhǎng)設(shè)置為Δt=1.93×10-5s，對(duì)應(yīng)螺旋槳旋轉(zhuǎn)角度0.25°。空化數(shù)σ定義如下：

式中，p0為靜壓，pv為液體的飽和蒸汽壓，ρl為液體密度。

2 網(wǎng)格策略

2.1 網(wǎng)格局部加密

2.1.1 網(wǎng)格劃分

首先采用六面體網(wǎng)格劃分非結(jié)構(gòu)網(wǎng)格，其體網(wǎng)格變化率為2。然后對(duì)槳葉表面、導(dǎo)隨邊和梢部的網(wǎng)格均勻加密，其中導(dǎo)隨邊和梢部附近的網(wǎng)格尺寸為槳葉表面網(wǎng)格尺寸的1/4。接著在槳葉和槳轂表面添加邊界層網(wǎng)格，其高度變化率為1.2，第一層邊界層網(wǎng)格高度約為槳葉表面網(wǎng)格尺寸的1/10。最后對(duì)螺旋槳附近區(qū)域的全部網(wǎng)格進(jìn)行加密，網(wǎng)格尺寸為槳葉表面網(wǎng)格尺寸的2 倍。槳葉表面網(wǎng)格和螺旋槳附近網(wǎng)格如圖3所示。

若對(duì)近尾流場(chǎng)全部網(wǎng)格進(jìn)一步精細(xì)加密，將導(dǎo)致網(wǎng)格數(shù)量過大。為了在精確模擬梢渦空化的同時(shí)減少不必要的計(jì)算量，對(duì)梢渦渦管所在圓環(huán)柱內(nèi)的兩個(gè)區(qū)域進(jìn)行局部加密，如圖4（a）所示。

圖4（a）中藍(lán)色的網(wǎng)格加密區(qū)域記為A，每個(gè)槳葉附近的紅色區(qū)域?yàn)橐月菪龢行狞c(diǎn)為圓心，扇環(huán)內(nèi)徑r1=0.45D，扇環(huán)外徑R1=0.51D，旋轉(zhuǎn)角為36°，軸向從x=0拉伸到x=0.06D的扇環(huán)區(qū)域。在A區(qū)域中加密網(wǎng)格的原因是槳葉附近網(wǎng)格正交性較差而不適用于自適應(yīng)網(wǎng)格加密，因此需要提前在槳葉附近的梢渦區(qū)域布置足夠密的網(wǎng)格。

圖4（b）中紅色的網(wǎng)格加密區(qū)域記為B，其為x軸為對(duì)稱軸，圓環(huán)內(nèi)徑r2=0.42D，圓環(huán)外徑R2=0.51D，軸向從x=0 拉伸到x=0.25D的圓環(huán)柱區(qū)域。本文為了減少計(jì)算量?jī)H加密至軸向x=0.25D，若計(jì)算資源充足可盡量將該區(qū)域向后延伸。

為了比較說明A和B區(qū)域的網(wǎng)格分辨率對(duì)梢渦空化數(shù)值模擬的影響，根據(jù)上述網(wǎng)格劃分策略共劃分了四套網(wǎng)格，網(wǎng)格信息如表2所示，加密后的G3網(wǎng)格如圖4（b）所示。

表2 網(wǎng)格信息Tab.2 Information of different grids

2.1.2 計(jì)算結(jié)果與分析

基于上述四套網(wǎng)格對(duì)螺旋槳梢渦空化進(jìn)行數(shù)值模擬。表3列出了各套網(wǎng)格數(shù)值模擬結(jié)果的推力系數(shù)KT和扭矩系數(shù)KQ，結(jié)果顯示均與試驗(yàn)結(jié)果非常接近，證明本文采取的網(wǎng)格劃分方法和數(shù)值模擬方法比較可靠。

表3 各套網(wǎng)格計(jì)算結(jié)果的推力系數(shù)和扭矩系數(shù)Tab.3 Thrust coefficients and torque coefficients of calculation results of each grid

采用汽體體積分?jǐn)?shù)αv=0.2 的等值面表示梢渦空化，結(jié)果顯示數(shù)值模擬的梢渦空化形態(tài)均與試驗(yàn)比較吻合，如圖5所示。對(duì)比四套網(wǎng)格的計(jì)算結(jié)果，可以觀察到梢渦空化區(qū)域隨著網(wǎng)格的增加而不斷增大，其中G3與G4的計(jì)算結(jié)果差距相對(duì)較小。但是數(shù)值模擬的梢渦空化長(zhǎng)度與試驗(yàn)的差距均很大，仍需要進(jìn)一步提高梢渦區(qū)域的網(wǎng)格分辨率。比較G2與G3的數(shù)值模擬結(jié)果可知，對(duì)A區(qū)域的局部加密有利于提高梢渦空化數(shù)值模擬精度。

分別截取包含A和B區(qū)域在內(nèi)的兩個(gè)截面，其無量綱軸向渦量分布如圖6 所示。圖6（a）顯示，隨著網(wǎng)格分辨率的提高，梢渦上方的軸向渦量明顯增大，而下方的渦量變化很小。這表明梢部端點(diǎn)附近的流動(dòng)狀況更加復(fù)雜，速度壓力變化較大，需要更高的網(wǎng)格分辨率減少數(shù)值耗散。圖6（b）中只有G1的軸向渦量相對(duì)較小，而其它網(wǎng)格的計(jì)算結(jié)果差距很小，這是因?yàn)镚1 網(wǎng)格的整體網(wǎng)格分辨率以及其它三套網(wǎng)格的B區(qū)域網(wǎng)格分辨率均較低。比較兩個(gè)截面的無量綱軸向渦量分布可以發(fā)現(xiàn)，隨著梢渦從A區(qū)域發(fā)展到B區(qū)域，其軸向渦量衰減十分迅速，這表明需要在B區(qū)域中布置足夠密的網(wǎng)格以提高數(shù)值模擬精度。

從圖5與圖6可以觀察到，G3與G4的計(jì)算結(jié)果差距較小，為了節(jié)約計(jì)算資源以及提高計(jì)算效率，選擇G3網(wǎng)格作為下一步自適應(yīng)網(wǎng)格加密的初始網(wǎng)格。

2.2 自適應(yīng)網(wǎng)格加密

2.2.1 網(wǎng)格劃分

若對(duì)B區(qū)域的所有網(wǎng)格進(jìn)行加密，將大幅增加網(wǎng)格數(shù)量。本文采用基于Q準(zhǔn)則的自適應(yīng)網(wǎng)格加密方法，精確提高B區(qū)域中梢渦網(wǎng)格的網(wǎng)格分辨率，從而減少不必要的計(jì)算量。Q為速度梯度張量的第二不變量，計(jì)算公式如下：

式中，Ω為渦張量，S為應(yīng)變率張量。

自適應(yīng)網(wǎng)格加密流程如圖7所示。首先在一定工況下對(duì)局部加密后的網(wǎng)格進(jìn)行螺旋槳梢渦空化數(shù)值模擬，再根據(jù)Q方法識(shí)別計(jì)算結(jié)果中的梢渦結(jié)構(gòu)，然后使用topoSet工具標(biāo)記Q取值大于一定值的網(wǎng)格，最后利用refineMesh 工具對(duì)標(biāo)記的網(wǎng)格進(jìn)行八叉樹加密?？芍貜?fù)自適應(yīng)加密步驟，直到梢渦區(qū)域網(wǎng)格分辨率滿足精細(xì)數(shù)值模擬的需求。螺旋槳近尾流場(chǎng)渦結(jié)構(gòu)復(fù)雜，為了避免加密近尾流場(chǎng)中的無關(guān)渦結(jié)構(gòu)，在讀取并標(biāo)記網(wǎng)格時(shí)Q的最小閾值Qlimit應(yīng)當(dāng)設(shè)置得較大。

基于G3的梢渦空化數(shù)值模擬結(jié)果，標(biāo)記并加密了計(jì)算域內(nèi)軸向x/D=0.06到x/D=0.25的區(qū)域內(nèi)Q＞2.0×106的網(wǎng)格，加密后的網(wǎng)格記作G5，網(wǎng)格總數(shù)為2.13×107，自適應(yīng)加密后的網(wǎng)格如圖8所示。

除Q準(zhǔn)則之外，自適應(yīng)網(wǎng)格加密方法也可以通過Δ 準(zhǔn)則、λ2準(zhǔn)則和Ω方法等識(shí)別并加密梢渦區(qū)域的網(wǎng)格，或者通過氣體體積分?jǐn)?shù)或者絕對(duì)壓力值識(shí)別并加密梢渦空化區(qū)域的網(wǎng)格。

本文所用的自適應(yīng)網(wǎng)格加密并非嚴(yán)格意義上的自適應(yīng)網(wǎng)格加密，不是在每個(gè)時(shí)間步后都根據(jù)某一個(gè)物理量進(jìn)行網(wǎng)格加密。這是因?yàn)闆]有一個(gè)物理量能僅僅表示梢渦結(jié)構(gòu)，例如各種渦判定方法等識(shí)別的區(qū)域都會(huì)包括槳葉表面附近的網(wǎng)格，但靠近槳葉的網(wǎng)格正交性較差而不適用于OpenFOAM 的八叉樹加密。因此舍棄了基于dynamicMesh 工具實(shí)現(xiàn)的嚴(yán)格意義上的自適應(yīng)網(wǎng)格加密，換而使用基于topoSet 工具和refineMesh工具實(shí)現(xiàn)的“靜態(tài)”的自適應(yīng)網(wǎng)格加密。

2.2.2 計(jì)算結(jié)果與分析

基于自適應(yīng)網(wǎng)格加密后的G5 網(wǎng)格進(jìn)行了梢渦空化數(shù)值模擬。表4 列出了計(jì)算結(jié)果的推力系數(shù)KT和扭矩系數(shù)KQ，與試驗(yàn)結(jié)果基本符合。由于自適應(yīng)網(wǎng)格加密只針對(duì)梢渦區(qū)域的網(wǎng)格進(jìn)行加密，對(duì)螺旋槳推力和扭矩計(jì)算結(jié)果的影響非常有限。

表4 G5計(jì)算結(jié)果的推力系數(shù)和扭矩系數(shù)Tab.4 Thrust coefficient and torque coefficient of the calculation results with G5

計(jì)算結(jié)果的梢渦空化形態(tài)由汽體體積分?jǐn)?shù)αv的等值面表示，不同的αv取值對(duì)梢渦空化形態(tài)的影響較大。在此比較了兩套網(wǎng)格計(jì)算結(jié)果中的αv=0.2、0.5 和0.8 三個(gè)等值面，如圖9 所示。在三種氣體體積分?jǐn)?shù)的取值下，自適應(yīng)加密后的網(wǎng)格的計(jì)算結(jié)果均優(yōu)于未加密的網(wǎng)格，并且隨著氣體體積分?jǐn)?shù)的增大，梢渦空化形態(tài)差距也越來越大。結(jié)果表明自適應(yīng)網(wǎng)格加密方法通過準(zhǔn)確加密梢渦空化區(qū)域的網(wǎng)格，提高了梢渦空化數(shù)值模擬精度。

將G5計(jì)算結(jié)果的αv=0.2的等值面與試驗(yàn)進(jìn)行比較，如圖10所示。數(shù)值模擬的梢渦空化形態(tài)與試驗(yàn)非常吻合，可以模擬出梢渦空化較為明顯的卷起過程和梢渦空化規(guī)律性收縮的節(jié)點(diǎn)（紅色虛線框），說明本文建立的網(wǎng)格策略能夠較好地模擬梢渦空化。但是梢渦空化僅僅在軸向上發(fā)展至x=0.16D，在長(zhǎng)度上與試驗(yàn)差距較大。為了得到更符合試驗(yàn)的計(jì)算結(jié)果，需要進(jìn)一步提高網(wǎng)格分辨率或采取更合適的空化模型和湍流模擬方法。此外，數(shù)值模擬過高預(yù)報(bào)了導(dǎo)邊附近的片空化區(qū)域。

3 結(jié) 語(yǔ)

本文結(jié)合網(wǎng)格局部加密和自適應(yīng)網(wǎng)格加密方法，對(duì)E779A 螺旋槳進(jìn)行了螺旋槳梢渦空化數(shù)值模擬。數(shù)值模擬結(jié)果顯示，在近尾流場(chǎng)內(nèi)的梢渦空化形態(tài)與試驗(yàn)非常接近，成功模擬出梢渦空化的卷起現(xiàn)象和規(guī)律性收縮的節(jié)點(diǎn)。若需提高更遠(yuǎn)區(qū)域的梢渦數(shù)值模擬精度，則應(yīng)進(jìn)一步提高網(wǎng)格分辨率或者采用更合適的湍流處理方法和空化模型。研究建立了基于OpenFOAM 的螺旋槳梢渦空化數(shù)值模擬的網(wǎng)格策略，為將來提升螺旋槳梢渦空化預(yù)報(bào)準(zhǔn)確性奠定了基礎(chǔ)。