基于RANSE的螺旋槳模型敞水?dāng)?shù)值模擬方法研究

金奕星，吳乘勝，王建春，王 星

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

0 引 言

螺旋槳是應(yīng)用最為廣泛的船用推進(jìn)器[1]，其水動(dòng)力性能直接影響船舶的總體性能。作為基礎(chǔ)水動(dòng)力性能，螺旋槳敞水性能預(yù)報(bào)是船舶設(shè)計(jì)和總體性能預(yù)報(bào)中不可或缺的重要環(huán)節(jié)。數(shù)值計(jì)算作為與模型試驗(yàn)并駕齊驅(qū)的主要技術(shù)手段，廣泛應(yīng)用于螺旋槳敞水性能預(yù)報(bào)[2]。

經(jīng)過幾十年的發(fā)展，螺旋槳水動(dòng)力性能數(shù)值計(jì)算經(jīng)歷了從勢(shì)流方法到粘流方法的轉(zhuǎn)變。與勢(shì)流方法相比，粘流方法的流動(dòng)控制方程考慮了流體粘性的作用、旋渦的影響以及旋渦在流動(dòng)中的傳輸，能夠模擬螺旋槳附近流場(chǎng)的精細(xì)結(jié)構(gòu)，進(jìn)而能更準(zhǔn)確地預(yù)報(bào)水動(dòng)力性能。目前，粘流CFD計(jì)算方法已經(jīng)成為螺旋槳敞水性能預(yù)報(bào)的主要手段之一。

基于粘流CFD 方法國(guó)內(nèi)外不少研究人員都開展了螺旋槳敞水性能數(shù)值計(jì)算研究。Shotatro Uto[3]使用RANS 方程和B-L 零方程湍流模型對(duì)螺旋槳繞流場(chǎng)進(jìn)行了數(shù)值模擬。Sanchez Caja 等[4]采用芬蘭赫爾辛基技術(shù)大學(xué)開發(fā)的RANS 求解器FINFLO 對(duì)單槳進(jìn)行了水動(dòng)力分析，在設(shè)計(jì)點(diǎn)工況處得到的推力和扭矩系數(shù)與試驗(yàn)值相比，誤差在1.5%以內(nèi)。Feneno[5]采用RANS方法對(duì)大側(cè)斜螺旋槳水動(dòng)力性能進(jìn)行了預(yù)報(bào)，獲得的非定常計(jì)算結(jié)果與試驗(yàn)結(jié)果比較一致。Rhee 等[6]以非結(jié)構(gòu)化網(wǎng)格為基礎(chǔ)，結(jié)合RANS方程和k-ω湍流模型對(duì)五葉螺旋槳的敞水性能進(jìn)行了計(jì)算，所得的推力和扭矩與試驗(yàn)值相差在10%以內(nèi)。國(guó)內(nèi)的黃勝等[7]研究了不同湍流模型對(duì)螺旋槳水動(dòng)力性能計(jì)算的影響。龔呂等[8]采用非結(jié)構(gòu)化網(wǎng)格和標(biāo)準(zhǔn)k-ε湍流模型對(duì)六葉槳進(jìn)行了計(jì)算，所得的推力和扭矩系數(shù)與試驗(yàn)值誤差達(dá)到8%。洪方文、張志榮等[9]研究了網(wǎng)格尺度、幾何精細(xì)度表達(dá)以及邊界層網(wǎng)格形式對(duì)螺旋槳水動(dòng)力性能計(jì)算的影響。劉志華等[10]采用RANS方程結(jié)合RNGk-ε湍流模型，運(yùn)用多塊混合網(wǎng)格對(duì)螺旋槳敞水性能進(jìn)行了預(yù)報(bào)，與試驗(yàn)結(jié)果比較，滿足工程應(yīng)用需求。鄭巢生等[11]基于OpenFOAM 進(jìn)行了螺旋槳敞水性能計(jì)算，計(jì)算結(jié)果與試驗(yàn)數(shù)據(jù)符合較好。

綜合上述，由國(guó)內(nèi)外相關(guān)研究工作可見，對(duì)于螺旋槳敞水CFD 模擬這一類常見的典型應(yīng)用，即便是學(xué)術(shù)和技術(shù)素養(yǎng)很高的科技人員，獲得的結(jié)果也存在相當(dāng)?shù)牟町悺Ｒ虼耍趯?duì)同一類對(duì)象和問題的數(shù)值模擬中，如何盡可能消除人為因素影響從而獲得不“因人而異”的結(jié)果，是船舶水動(dòng)力CFD應(yīng)用和數(shù)值水池研發(fā)必須解決的問題；而在“屬性細(xì)分”的基礎(chǔ)上進(jìn)行“知識(shí)封裝”，是解決這一問題的有效途徑[12]。

影響螺旋槳敞水CFD模擬結(jié)果的因素很多，包括湍流模型、網(wǎng)格劃分以及旋轉(zhuǎn)的模擬方法等多個(gè)方面，且相互之間還可能存在耦合影響。其中，旋轉(zhuǎn)的模擬方法以及相關(guān)的網(wǎng)格處理方式，是螺旋槳敞水計(jì)算特有的典型問題，本文的研究工作主要聚焦于此。

作為一種旋轉(zhuǎn)驅(qū)動(dòng)的流動(dòng)問題，螺旋槳敞水的數(shù)值模擬方法至少有三種——非慣性坐標(biāo)系方法、多參考坐標(biāo)系方法和滑移網(wǎng)格方法，每種方法中網(wǎng)格處理方式又有多種。而數(shù)值模擬方法的選擇和網(wǎng)格處理方式的不同，是影響CFD 模擬結(jié)果的重要人為因素；這些因素的影響研究，也是“知識(shí)封裝”中的重要“知識(shí)”來源。

本文采用基于RANSE 的船舶水動(dòng)力學(xué)CFD 求解器，開展螺旋槳模型敞水?dāng)?shù)值模擬研究。當(dāng)前求解器擁有三種可用于螺旋槳敞水模擬的計(jì)算方法：非慣性坐標(biāo)系方法、多參考坐標(biāo)系方法和滑移網(wǎng)格方法，其中后兩種方法需要將計(jì)算域拆分為旋轉(zhuǎn)域和靜止域，不同計(jì)算域上求解的控制方程也存在差異。為了實(shí)現(xiàn)不同計(jì)算域之間的耦合計(jì)算，當(dāng)前CFD求解器引入了交界面方法，通過構(gòu)造虛擬單元進(jìn)行域之間流場(chǎng)信息的插值傳遞。為了解決交界面兩側(cè)網(wǎng)格單元類型和網(wǎng)格尺度差異，以及網(wǎng)格相對(duì)位置關(guān)系變化造成的求解器計(jì)算穩(wěn)定性和計(jì)算精度下降問題，求解器采用了基于Laplace 權(quán)函數(shù)的插值計(jì)算方法，充分利用了貢獻(xiàn)單元和插值點(diǎn)的距離關(guān)系和方向關(guān)系。最后，本文以典型螺旋槳模型為對(duì)象，設(shè)計(jì)了包括三種計(jì)算方法和三種網(wǎng)格處理方式的計(jì)算方案，研究不同計(jì)算方法和網(wǎng)格形式，以及交界面位置對(duì)計(jì)算結(jié)果的影響。計(jì)算結(jié)果表明，采用基于Laplace 權(quán)函數(shù)插值計(jì)算的交界面方法，能夠提高螺旋槳敞水計(jì)算穩(wěn)定性，保證計(jì)算精度。與此同時(shí)，研究給出了采用當(dāng)前求解器進(jìn)行螺旋槳敞水性能數(shù)值計(jì)算的推薦方案，并采用推薦方案開展了典型螺旋槳模型敞水性能曲線數(shù)值計(jì)算，結(jié)果與模型試驗(yàn)符合良好。

1 數(shù)值計(jì)算方法

本文使用的船舶水動(dòng)力CFD 求解器NaViiX（Naval Hydrodynamics Oriented CFD Solvers），由中國(guó)船舶科學(xué)研究中心獨(dú)立自主研發(fā)，具有完全自主知識(shí)產(chǎn)權(quán)。求解器基于求解RANS方程（Reynolds Averaged Navier-Stokes Equations），采用有限體積法（Finite Volume Method，F(xiàn)VM）離散控制方程，其中對(duì)流項(xiàng)采用二階迎風(fēng)格式，擴(kuò)散項(xiàng)離散采用中心差分格式，速度壓力耦合采用基于SIMPLE（Semi-Implicit Method for Pressure Linked Equations）算法的同位網(wǎng)格（Collocated Grid）技術(shù)進(jìn)行解耦，具有以下功能：

（1）具備三維航行體湍流繞流CFD模擬能力；

（2）支持結(jié)構(gòu)化網(wǎng)格、非結(jié)構(gòu)化網(wǎng)格、混合網(wǎng)格、交界面網(wǎng)格和滑移網(wǎng)格；

（3）支持標(biāo)準(zhǔn)k-ε、RNGk-ε、k-ω和SSTk-ω等常用湍流模型；

（4）支持慣性坐標(biāo)系、非慣性坐標(biāo)系和多參考坐標(biāo)系求解；

（5）支持MPI并行計(jì)算。

本文所開展的螺旋槳模型敞水?dāng)?shù)值模擬研究，采用了非慣性坐標(biāo)系方法、多參考坐標(biāo)系方法、滑移網(wǎng)格方法，以及多計(jì)算域耦合計(jì)算的交界面方法，下面對(duì)其進(jìn)行簡(jiǎn)要介紹。

1.1 非慣性坐標(biāo)系方法

非慣性坐標(biāo)系方法（Non-inertial Reference Frames Method）是在非慣性參考坐標(biāo)系下求解流動(dòng)控制方程。對(duì)于軸對(duì)稱的單一旋轉(zhuǎn)部件在做周期性旋轉(zhuǎn)運(yùn)動(dòng)時(shí)，在慣性坐標(biāo)系下觀察，流動(dòng)是非定常的；但如果在轉(zhuǎn)動(dòng)坐標(biāo)系（非慣性坐標(biāo)系的一種特例）下觀察，則流動(dòng)是定常的。對(duì)于螺旋槳敞水性能計(jì)算這類問題，在非慣性系下，可以將非定常問題轉(zhuǎn)化為準(zhǔn)定常問題進(jìn)行求解，在實(shí)踐中有著較為廣泛的應(yīng)用。

流體力學(xué)問題中，非慣性坐標(biāo)系是基于地面慣性坐標(biāo)系定義的一個(gè)參考坐標(biāo)系。非慣性坐標(biāo)系相對(duì)于慣性坐標(biāo)系，通過附加離心慣性力和科氏慣性力，質(zhì)量、動(dòng)量和能量守恒定律依舊適用，只是控制方程需要進(jìn)行相應(yīng)修改，具體修改方法見文獻(xiàn)[13]。

非慣性系下絕對(duì)速度形式控制方程求解與慣性系下控制方程求解的差異主要有三點(diǎn)：控制體邊界速度修正、牽連慣性力源項(xiàng)動(dòng)量方程中增加和邊界條件的相應(yīng)修改，具體見文獻(xiàn)[13]。

1.2 多參考坐標(biāo)系方法

多參考坐標(biāo)系方法[14]（Multiple Reference Frames Method）的主要思想是，將計(jì)算域劃分為旋轉(zhuǎn)域和靜止域，用旋轉(zhuǎn)域包含旋轉(zhuǎn)部件，用靜止域包含靜止部件，旋轉(zhuǎn)域在非慣性坐標(biāo)系下求解流動(dòng)控制方程，靜止域在慣性坐標(biāo)系下求解流動(dòng)控制方程。旋轉(zhuǎn)域和靜止域之間通過交界面（Interface）連接，這里的交界面是重疊的兩個(gè)面，這兩個(gè)面上的網(wǎng)格點(diǎn)無需一一對(duì)應(yīng)，即在交界面上只要保證面搭接而不需要保證點(diǎn)搭接；域之間通過交界面方法進(jìn)行流場(chǎng)信息的傳遞，從而實(shí)現(xiàn)域之間的耦合計(jì)算。

多參考坐標(biāo)系方法的主要應(yīng)用場(chǎng)景是多體相對(duì)運(yùn)動(dòng)的流場(chǎng)數(shù)值模擬，尤其是同時(shí)存在轉(zhuǎn)動(dòng)部件和靜止部件的繞流模擬。因此，多參考坐標(biāo)系方法也可用于螺旋槳敞水性能數(shù)值計(jì)算。由于該方法所采用的計(jì)算網(wǎng)格劃分方式相當(dāng)靈活，在實(shí)際工程應(yīng)用中被廣泛采用。

1.3 滑移網(wǎng)格方法

滑移網(wǎng)格（Sliding Mesh）是動(dòng)網(wǎng)格的一種特例，與常規(guī)動(dòng)網(wǎng)格方法不同，滑移網(wǎng)格方法中的網(wǎng)格節(jié)點(diǎn)只做剛體運(yùn)動(dòng)，控制體單元形狀和網(wǎng)格點(diǎn)之間的距離都保持不變。

滑移網(wǎng)格方法和多參考坐標(biāo)系方法有許多類似之處，都是將計(jì)算域劃分為旋轉(zhuǎn)域和靜止域，域之間通過交界面連接，通過構(gòu)造插值點(diǎn)，以插值的方式進(jìn)行流場(chǎng)信息傳遞。二者最大的不同之處在于滑移網(wǎng)格是非定常計(jì)算方法，多參考坐標(biāo)系是準(zhǔn)定常計(jì)算方法；滑移網(wǎng)格的旋轉(zhuǎn)域隨著時(shí)間的推進(jìn)是在真實(shí)地轉(zhuǎn)動(dòng)，交界面兩側(cè)的位置信息也隨時(shí)間改變；而多參考系方法的旋轉(zhuǎn)域?qū)嶋H上是不轉(zhuǎn)動(dòng)的，只是在非慣性旋轉(zhuǎn)坐標(biāo)系下進(jìn)行求解。

滑移網(wǎng)格方法交界面上插值點(diǎn)構(gòu)造、宿主單元搜索和插值方法與多參考坐標(biāo)系方法是一樣的。需要注意的是，在每個(gè)時(shí)間步開始計(jì)算之前，都需要重新進(jìn)行滑移交界面的構(gòu)造、宿主單元搜索以及變量的插值計(jì)算。

1.4 交界面方法

無論是多參考坐標(biāo)系方法還是滑移網(wǎng)格方法，都面臨計(jì)算域拆分為多個(gè)子計(jì)算域進(jìn)行網(wǎng)格生成的問題，這就要求流場(chǎng)求解器具備同時(shí)讀入多套網(wǎng)格進(jìn)行耦合計(jì)算的能力。無論是多參考坐標(biāo)系方法還是滑移網(wǎng)格方法，不同計(jì)算域求解的控制方程存在差異，這是由計(jì)算對(duì)象是否運(yùn)動(dòng)決定的。盡管不同計(jì)算域中求解的控制方程可能存在差異，但必須保證不同計(jì)算域之間的流場(chǎng)是緊密耦合的，因?yàn)槟M的繞流場(chǎng)本身是一個(gè)連續(xù)流場(chǎng)，是人為地將其拆分為多個(gè)計(jì)算域進(jìn)行數(shù)值模擬，這樣做的目的是方便實(shí)現(xiàn)多體相對(duì)運(yùn)動(dòng)問題的流場(chǎng)數(shù)值模擬。為了實(shí)現(xiàn)不同計(jì)算域之間的緊密耦合計(jì)算，就不得不采用交界面方法建立計(jì)算域之間的聯(lián)系。

交界面方法是通過計(jì)算域之間的兩個(gè)面來傳遞域之間的流場(chǎng)信息。這兩個(gè)面分別歸屬于兩個(gè)不同的計(jì)算域，因此無論是交界面上的面網(wǎng)格還是交界面附近區(qū)域的體網(wǎng)格，都不存在相互的制約與限制。這樣是方便了網(wǎng)格的生成，但卻大大增加了求解器實(shí)現(xiàn)多計(jì)算域耦合計(jì)算的難度，因?yàn)榻唤缑鎯蓚?cè)的網(wǎng)格情況變得復(fù)雜。交界面兩側(cè)可能存在網(wǎng)格類型不同、網(wǎng)格尺度不同和網(wǎng)格單元的相互位置關(guān)系變化的情況，這些因素都會(huì)導(dǎo)致域之間流場(chǎng)信息傳遞不準(zhǔn)確，進(jìn)而導(dǎo)致域之間流場(chǎng)出現(xiàn)間斷，最終導(dǎo)致計(jì)算發(fā)散，或者導(dǎo)致人為拆分計(jì)算域產(chǎn)生的計(jì)算誤差太大，無法滿足應(yīng)用要求。

為了解決上述不得不面對(duì)的難題，在構(gòu)造完交界面插值點(diǎn)和完成插值點(diǎn)在相鄰計(jì)算域中宿主單元搜索之后，本文采用了一種基于Laplace 權(quán)函數(shù)的插值計(jì)算方法，通過插值的方式傳遞域之間的流場(chǎng)信息。插值計(jì)算公式如下：

式中，U是流場(chǎng)變量，如速度、壓強(qiáng)等，n是貢獻(xiàn)單元的數(shù)量，Ui是第i個(gè)貢獻(xiàn)單元的流場(chǎng)變量值，φi是第i個(gè)貢獻(xiàn)單元的插值系數(shù)，φi需要滿足插值系數(shù)和為1的條件。

采用基于Laplace權(quán)函數(shù)的插值計(jì)算方法時(shí)，插值系數(shù)φi需要滿足Laplace算子

式中，n是貢獻(xiàn)單元數(shù)量，Xi=(xi,yi,zi)是貢獻(xiàn)單元坐標(biāo)，Xp=(xp,yp,zp)是目標(biāo)插值點(diǎn)的坐標(biāo)，那么插值系數(shù)可以定義為

式中，λx、λz、λz就是拉格朗日乘數(shù)，可以根據(jù)貢獻(xiàn)單元坐標(biāo)和目標(biāo)插值點(diǎn)坐標(biāo)計(jì)算得到。通過上述方法，就可以基于Laplace權(quán)函數(shù)計(jì)算出所有貢獻(xiàn)單元的插值系數(shù)。

相比普通的插值系數(shù)構(gòu)造方式，Laplace 權(quán)函數(shù)插值計(jì)算方法用到了更多貢獻(xiàn)單元，并且同時(shí)考慮了貢獻(xiàn)單元和插值點(diǎn)的距離關(guān)系和方向關(guān)系，因此具有更高的插值精度，而且該插值系數(shù)計(jì)算方式適用于各種單元類型的網(wǎng)格。為了應(yīng)對(duì)局部網(wǎng)格尺度差異巨大的問題，本文還對(duì)插值系數(shù)進(jìn)行了限制，將插值系數(shù)限定在0～2 的范圍內(nèi)，這樣可以有效去除負(fù)值和過大的插值系數(shù)，增加流場(chǎng)求解器的穩(wěn)定性。

圖1給出了一典型算例，算例中采用滑移網(wǎng)格方法進(jìn)行非定常數(shù)值計(jì)算，交界面分別使用普通插值方法和Laplace 權(quán)函數(shù)插值方法進(jìn)行流場(chǎng)信息傳遞，圖中給出了兩種插值方法計(jì)算得到的螺旋槳推力曲線隨時(shí)間變化的對(duì)比。從圖中可以看出：使用普通插值方法，在計(jì)算的某些時(shí)刻，螺旋槳推力會(huì)出現(xiàn)脈沖式波動(dòng)，并且振蕩幅度很大，偏離正常值；使用Laplace 權(quán)函數(shù)插值方法后，在整個(gè)計(jì)算過程中，螺旋槳推力呈現(xiàn)出合理的非定常波動(dòng)，異常的脈沖式波動(dòng)現(xiàn)象消失了。

圖1 螺旋槳推力隨時(shí)間變化比較Fig.1 Comparison of propeller thrusts varying with time

可見，CFD 求解器中，交界面采用Laplace 權(quán)函數(shù)插值方法進(jìn)行域之間流場(chǎng)信息傳遞，可以提高求解器計(jì)算穩(wěn)定性，還可以提高計(jì)算結(jié)果的準(zhǔn)確性。

2 計(jì)算方案設(shè)計(jì)

2.1 研究對(duì)象與計(jì)算工況

本文針對(duì)一固定螺距四葉螺旋槳，開展螺旋槳模型敞水?dāng)?shù)值模擬研究。固定螺距四葉槳是散貨船、油船和集裝箱船三大主力船型廣泛采用的推進(jìn)器，具有較強(qiáng)的代表性；同時(shí)，本文選用的螺旋槳模型，在中國(guó)船舶科學(xué)研究中心進(jìn)行過系統(tǒng)的敞水性能基準(zhǔn)檢驗(yàn)?zāi)Ｐ驮囼?yàn)和不確定度分析，擁有高準(zhǔn)確度和高置信度的基準(zhǔn)檢驗(yàn)?zāi)Ｐ驮囼?yàn)數(shù)據(jù)[16]，適用于CFD 計(jì)算方法的考核和結(jié)果的驗(yàn)證。螺旋槳模型幾何外形如圖2所示，主尺度參數(shù)列于表1。

表1 螺旋槳模型主尺度參數(shù)Tab.1 Principal parameters of propeller model

圖2 螺旋槳模型幾何外形Fig.2 Geometry of propeller model

數(shù)值模擬工況為進(jìn)速系數(shù)J=0.4，該工況在螺旋槳設(shè)計(jì)點(diǎn)附近，開展了多輪次模型試驗(yàn)，試驗(yàn)數(shù)據(jù)準(zhǔn)確度高、置信度高。數(shù)值模擬中，流體介質(zhì)屬性及螺旋槳轉(zhuǎn)速與模型試驗(yàn)保持一致，流體密度為ρ=996.76 kg/m3，流體粘性為μ=0.000 868 kg/（m·s），螺旋槳轉(zhuǎn)速n=18 r/s，根據(jù)進(jìn)速系數(shù)公式（5）計(jì)算得到來流速度為V=1.8 m/s。

2.2 計(jì)算方案設(shè)計(jì)

本文主要研究不同計(jì)算方法、網(wǎng)格形式、以及交界面位置對(duì)螺旋槳敞水計(jì)算結(jié)果的影響。計(jì)算所使用的CFD 求解器支持采用非慣性坐標(biāo)系方法、多參考坐標(biāo)系方法和滑移網(wǎng)格方法進(jìn)行螺旋槳敞水性能計(jì)算，支持結(jié)構(gòu)化網(wǎng)格、非結(jié)構(gòu)化網(wǎng)格、混合網(wǎng)格和交界面網(wǎng)格對(duì)計(jì)算域進(jìn)行劃分。基于上述三種計(jì)算方法，結(jié)合三種不同的網(wǎng)格形式，設(shè)計(jì)了五種計(jì)算方案，如表2所示。

表2 計(jì)算方案Tab.2 Introduction of calculation schemes

五種計(jì)算方案除了計(jì)算方法和網(wǎng)格形式不同外，其他計(jì)算參數(shù)設(shè)置保持一致。基于中國(guó)船舶科學(xué)研究中心關(guān)于螺旋槳敞水性能CFD計(jì)算的實(shí)踐，結(jié)合當(dāng)前CFD求解器的功能特性，數(shù)值模擬采用的設(shè)置為：對(duì)流項(xiàng)采用二階迎風(fēng)格式，擴(kuò)散項(xiàng)離散采用帶非正交修正的中心差分格式，代數(shù)方程組求解采用穩(wěn)定雙共軛梯度法（BiCGSTAB），湍流模擬采用常用的RNGk-ε湍流模型。

3 螺旋槳敞水?dāng)?shù)值模擬

3.1 網(wǎng)格收斂性分析

網(wǎng)格數(shù)量對(duì)螺旋槳敞水計(jì)算結(jié)果的影響主要來自于槳葉表面以及槳葉附近區(qū)域的網(wǎng)格，而上述計(jì)算方案中槳葉表面和槳葉附近區(qū)域都采用非結(jié)構(gòu)化網(wǎng)格，并且不同方案槳葉附近網(wǎng)格密度基本一致，因此以方案三為代表開展網(wǎng)格收斂性研究。

為了研究網(wǎng)格數(shù)量對(duì)螺旋槳敞水性能計(jì)算的影響，采用三套不同疏密程度的網(wǎng)格計(jì)算J=0.4工況下螺旋槳模型的推力系數(shù)KT和扭矩系數(shù)KQ。三套網(wǎng)格具體信息及對(duì)應(yīng)的計(jì)算結(jié)果列于表3 中，水動(dòng)力系數(shù)KT和KQ計(jì)算結(jié)果隨網(wǎng)格數(shù)變化如圖3所示。

圖3 水動(dòng)力系數(shù)計(jì)算結(jié)果隨網(wǎng)格數(shù)變化Fig.3 Computational results of hydrodynamic coefficients varying with grid number

表3 計(jì)算網(wǎng)格與計(jì)算結(jié)果Tab.3 Computational meshes and computational results

推力系數(shù)和扭矩系數(shù)計(jì)算公式分別為

式中，T和Q分別是螺旋槳推力和扭矩值，ρ為密度，n為轉(zhuǎn)速，D為槳葉直徑。

從計(jì)算結(jié)果可以看出，計(jì)算所得螺旋槳推力系數(shù)和扭矩系數(shù)隨網(wǎng)格數(shù)量增加都是收斂的，細(xì)網(wǎng)格和中網(wǎng)格計(jì)算結(jié)果之間的差別僅為0.6%左右。本文后續(xù)計(jì)算中，皆采用細(xì)網(wǎng)格。

3.2 計(jì)算網(wǎng)格與邊界條件

根據(jù)上述五種設(shè)計(jì)方案的需求，結(jié)合網(wǎng)格依賴性研究結(jié)論，劃分了三種計(jì)算網(wǎng)格；三種網(wǎng)格的計(jì)算域都為圓柱體，計(jì)算域前端位于螺旋槳槳盤面前方6.0D（D為槳盤面直徑）處，周向邊界到槳軸中心距離也為6.0D，計(jì)算域后端位于槳盤面后方12.0D處。

網(wǎng)格一采用單一計(jì)算域，全域采用非結(jié)構(gòu)化網(wǎng)格；網(wǎng)格二采用交界面方式將計(jì)算域劃分為內(nèi)域（旋轉(zhuǎn)域）和外域（靜止域），內(nèi)外域皆采用非結(jié)構(gòu)化網(wǎng)格，域之間采用交界面連接；網(wǎng)格三同樣采用交界面方式將計(jì)算域劃分為內(nèi)域和外域，內(nèi)域采用非結(jié)構(gòu)化網(wǎng)格，外域采用結(jié)構(gòu)化網(wǎng)格，內(nèi)外域以交界面連接。三種計(jì)算網(wǎng)格如圖4所示。

圖4 三種計(jì)算網(wǎng)格Fig.4 Three types of computational meshes

為了能更準(zhǔn)確地反映出網(wǎng)格形式對(duì)計(jì)算結(jié)果的影響，三種網(wǎng)格槳葉表面網(wǎng)格保持一致，槳葉附近區(qū)域網(wǎng)格密度基本一致，網(wǎng)格單元總數(shù)也基本一致。為此，雖然網(wǎng)格一不存在交界面，但也在相應(yīng)區(qū)域進(jìn)行了局部網(wǎng)格密度控制。網(wǎng)格一單元總數(shù)為89 萬左右；網(wǎng)格二內(nèi)域（旋轉(zhuǎn)域）單元數(shù)為68 萬左右，外域（靜止域）單元數(shù)為20萬左右，網(wǎng)格單元總數(shù)為88萬左右；網(wǎng)格三內(nèi)域（旋轉(zhuǎn)域）網(wǎng)格單元與網(wǎng)格二相同，外域（靜止域）采用結(jié)構(gòu)化網(wǎng)格，單元數(shù)也為20萬左右，網(wǎng)格單元總數(shù)為88萬左右。

由于多參考坐標(biāo)系方法和滑移網(wǎng)格方法對(duì)計(jì)算域劃分要求一致，兩種方法可以使用同樣的計(jì)算網(wǎng)格。因此，給出計(jì)算方案和計(jì)算網(wǎng)格對(duì)應(yīng)關(guān)系為：方案一采用網(wǎng)格一進(jìn)行計(jì)算，方案二和方案四采用網(wǎng)格二進(jìn)行計(jì)算，方案三和方案五采用網(wǎng)格三進(jìn)行計(jì)算。

計(jì)算中，三種網(wǎng)格都涉及的邊界條件包括：入口邊界條件、出口邊界條件和物面邊界條件，其中網(wǎng)格二和網(wǎng)格三還涉及交界面邊界條件。入口邊界給定來流速度；出口邊界給定環(huán)境壓力；物面邊界采用粘性無滑移邊界條件，引入標(biāo)準(zhǔn)壁面函數(shù)；交界面邊界通過構(gòu)造虛擬單元來封閉邊界單元，插值獲取流場(chǎng)信息。

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

使用自主CFD 求解器NaViiX，對(duì)前述五種計(jì)算方案結(jié)合對(duì)應(yīng)計(jì)算網(wǎng)格開展螺旋槳模型敞水?dāng)?shù)值模擬，從計(jì)算結(jié)果和計(jì)算耗時(shí)兩個(gè)方面，對(duì)不同方案的計(jì)算結(jié)果進(jìn)行對(duì)比分析，從中選出推薦計(jì)算方案。

非慣性坐標(biāo)系方法和多參考坐標(biāo)系方法都是準(zhǔn)定常計(jì)算方法，收斂判據(jù)為通量殘值下降6 個(gè)量級(jí)，計(jì)算最大迭代步設(shè)置為1 000；滑移網(wǎng)格方法是非定常方法，設(shè)定每個(gè)真實(shí)時(shí)間步轉(zhuǎn)過3°，根據(jù)轉(zhuǎn)速確定真實(shí)時(shí)間步長(zhǎng)Δt，設(shè)定計(jì)算540個(gè)狀態(tài)（即螺旋槳轉(zhuǎn)過4.5圈）。

為了更好地比較五種計(jì)算方案，將計(jì)算得到的推力系數(shù)、扭矩系數(shù)和效率與模型試驗(yàn)數(shù)據(jù)進(jìn)行定量比較。對(duì)于方案四和方案五，由于結(jié)果存在波動(dòng)，采用最后30 個(gè)狀態(tài)的結(jié)果取平均值作為計(jì)算結(jié)果。與此同時(shí)，還考察了各方案的計(jì)算時(shí)長(zhǎng)，因?yàn)橛?jì)算耗時(shí)也是實(shí)際應(yīng)用中比較關(guān)注的因素。推力和扭矩系數(shù)計(jì)算見式（6）和式（7），敞水效率計(jì)算公式為

在進(jìn)速系數(shù)J=0.4 的計(jì)算工況下，模型試驗(yàn)給出的螺旋槳推力系數(shù)KT=0.210 0，扭矩系數(shù)10KQ=0.259 9，敞水效率η=0.514 4。表4給出了五種方案的螺旋槳推力系數(shù)、扭矩系數(shù)和敞水效率計(jì)算結(jié)果及其與模型試驗(yàn)結(jié)果之間的誤差以及各方案的計(jì)算耗時(shí)。

表4 計(jì)算結(jié)果與計(jì)算耗時(shí)Tab.4 Calculation results and calculation elapsed time

從表4 中可以看到，當(dāng)前工況下，五種計(jì)算方案的計(jì)算結(jié)果都與模型試驗(yàn)結(jié)果相當(dāng)接近，五種方案計(jì)算結(jié)果之間的螺旋槳推力差異在0.8%以內(nèi)，扭矩差異在0.3%以內(nèi)。從五個(gè)方案計(jì)算結(jié)果的對(duì)比可以看出，采用Laplace 權(quán)函數(shù)插值方法進(jìn)行域之間流場(chǎng)信息的插值傳遞，求解器能夠處理交界面兩側(cè)不同網(wǎng)格類型、不同網(wǎng)格尺度，以及交界面相互位置關(guān)系變化的復(fù)雜情況，并且計(jì)算穩(wěn)定，計(jì)算精度與單計(jì)算域計(jì)算結(jié)果精度相當(dāng)。

對(duì)于螺旋槳敞水這類周期性旋轉(zhuǎn)運(yùn)動(dòng)，準(zhǔn)定常計(jì)算方法與非定常計(jì)算方法結(jié)果基本一致，而采用準(zhǔn)定常計(jì)算方法可以大大減少計(jì)算耗時(shí)。同為準(zhǔn)定常方法，多參考坐標(biāo)系方法相比非慣性坐標(biāo)系方法還能節(jié)約20%左右的計(jì)算耗時(shí)，主要是因?yàn)榫W(wǎng)格拆分后，控制方程系數(shù)矩陣規(guī)模變小，迭代計(jì)算耗時(shí)減少。此外，因?yàn)槎鄥⒖甲鴺?biāo)系方法采用多計(jì)算域的緣故，在網(wǎng)格劃分過程中的自由度更大，更容易生成高質(zhì)量計(jì)算網(wǎng)格。

對(duì)于網(wǎng)格形式，多參考坐標(biāo)系內(nèi)部旋轉(zhuǎn)域都采用非結(jié)構(gòu)化網(wǎng)格，差別在于外部靜止域。由于槳轂表面為圓柱面，外域如果采用結(jié)構(gòu)化網(wǎng)格，可以在保證物面網(wǎng)格足夠貼合槳轂的同時(shí)，減少外域網(wǎng)格量；此外，結(jié)構(gòu)化網(wǎng)格計(jì)算結(jié)果在后處理時(shí)，流場(chǎng)的可視化顯示效果一般更佳。

綜合上述分析，以上方案中，方案三為最優(yōu)計(jì)算方案，即針對(duì)當(dāng)前CFD求解器，推薦采用多參考坐標(biāo)系方法，結(jié)合內(nèi)域采用非結(jié)構(gòu)化網(wǎng)格，外域采用結(jié)構(gòu)化網(wǎng)格的網(wǎng)格劃分形式進(jìn)行螺旋槳敞水性能數(shù)值計(jì)算。

3.4 交界面位置對(duì)計(jì)算結(jié)果影響研究

由于多參考坐標(biāo)系方法涉及計(jì)算域拆分和交界面處理方法，并且交界面方法是多計(jì)算域耦合計(jì)算的關(guān)鍵，交界面的位置會(huì)直接影響交界面兩側(cè)的網(wǎng)格單元尺度和交界面兩側(cè)流場(chǎng)變化劇烈程度，因此需要進(jìn)一步研究交界面距離螺旋槳的位置對(duì)計(jì)算結(jié)果的影響。基于推薦方案所采用的網(wǎng)格三，通過減小和增大交界面與螺旋槳的距離，得到了兩套新的計(jì)算網(wǎng)格，分別是網(wǎng)格四和網(wǎng)格五，如圖5 所示。

圖5 網(wǎng)格交界面位置對(duì)比Fig.5 Comparison of the locations of grid interfaces

采用上面三套計(jì)算網(wǎng)格開展交界面位置對(duì)計(jì)算結(jié)果影響的研究，三套網(wǎng)格的內(nèi)域都是用圓柱體將槳葉包裹在里面，內(nèi)域采用非結(jié)構(gòu)化網(wǎng)格，外域用一個(gè)更大的圓柱體代表計(jì)算域，采用結(jié)構(gòu)化網(wǎng)格，槳葉表面網(wǎng)格基本一致。三套網(wǎng)格的主要差異在于交界面位置：網(wǎng)格三交界面前端距離槳盤面距離是0.4D（D為槳盤面直徑），交界面后端距離槳盤面距離是0.2D，交界面周向距離槳軸中心為0.7D；網(wǎng)格四交界面前端和后端到槳盤面距離同為0.15D，交界面周向距離槳軸中心0.6D；網(wǎng)格五交界面前端距離槳盤面距離是1.0D，交界面后端距離槳盤面距離是0.8D，交界面周向距離槳軸中心為0.8D。

計(jì)算工況和求解參數(shù)設(shè)置與前面的完全一致，采用推薦計(jì)算方案，使用三套不同交界面位置的網(wǎng)格進(jìn)行螺旋槳敞水性能計(jì)算，將計(jì)算結(jié)果與模型試驗(yàn)結(jié)果進(jìn)行對(duì)比，如表5所示。

表5 不同交界面位置的計(jì)算結(jié)果Tab.5 Computational results of different interface locations

由表5 可見，交界面的位置會(huì)對(duì)螺旋槳敞水性能計(jì)算結(jié)果產(chǎn)生一定的影響：交界面離槳盤面越近，螺旋槳推力和扭矩計(jì)算結(jié)果相對(duì)試驗(yàn)結(jié)果較小；交界面離槳盤面越遠(yuǎn)，螺旋槳推力和扭矩計(jì)算結(jié)果相比試驗(yàn)結(jié)果較大。經(jīng)分析主要原因在于：由于螺旋槳附近流場(chǎng)變化劇烈，交界面離螺旋槳過近時(shí)，內(nèi)域和外域之間流場(chǎng)信息傳遞過程中存在難以避免的誤差，會(huì)對(duì)螺旋槳水動(dòng)力計(jì)算結(jié)果產(chǎn)生一定影響；交界面離螺旋槳過遠(yuǎn)時(shí)，螺旋槳附近網(wǎng)格密度又難以高質(zhì)量地控制，影響了槳葉附近的網(wǎng)格空間分辨率，進(jìn)而會(huì)影響螺旋槳水動(dòng)力計(jì)算結(jié)果。

由此可見，選取合適的網(wǎng)格交界面位置對(duì)于螺旋槳敞水性能的準(zhǔn)確計(jì)算也相當(dāng)重要。本文推薦采用中等距離交界面位置進(jìn)行計(jì)算網(wǎng)格劃分。

3.5 推薦計(jì)算方案應(yīng)用測(cè)試

以上開展了計(jì)算方法、網(wǎng)格形式和交界面位置對(duì)螺旋槳敞水性能數(shù)值計(jì)算結(jié)果的影響研究，并獲得了推薦數(shù)值模擬方法，即采用多參考坐標(biāo)系方法，網(wǎng)格劃分采用內(nèi)域非結(jié)構(gòu)化網(wǎng)格、外域結(jié)構(gòu)化網(wǎng)格的形式，網(wǎng)格交界面采用中等距離位置。

為了進(jìn)一步驗(yàn)證推薦數(shù)值模擬方法的可靠性，本節(jié)開展了更多工況的計(jì)算，獲得了螺旋槳敞水性能曲線，并與模型試驗(yàn)結(jié)果進(jìn)行了對(duì)比。

表6給出了螺旋槳模型敞水性能計(jì)算結(jié)果與試驗(yàn)結(jié)果的比較。

表6 螺旋槳敞水性能比較Tab.6 Comparison of open-water performance of propeller

從表中可以看出：在CFD模擬的工況范圍內(nèi)，推力系數(shù)和扭矩系數(shù)計(jì)算結(jié)果與模型試驗(yàn)結(jié)果的趨勢(shì)是一致的，都在設(shè)計(jì)點(diǎn)附近工況符合較好，隨著進(jìn)速系數(shù)越往兩端差異較大；在設(shè)計(jì)點(diǎn)附近，推力系數(shù)與試驗(yàn)結(jié)果偏差3%左右，扭矩系數(shù)與試驗(yàn)結(jié)果偏差2%左右；螺旋槳敞水效率計(jì)算結(jié)果與試驗(yàn)結(jié)果符合較好，大部分工況的偏差都在2%以內(nèi)。圖6給出了0.1～0.8 進(jìn)速系數(shù)下槳葉表面壓強(qiáng)分布云圖，圖7 給出了螺旋槳模型敞水性能曲線。

圖6 槳葉表面壓強(qiáng)分布云圖Fig.6 Pressure contour of up suface of the blades

圖7 螺旋槳敞水性能曲線比較Fig.7 Comparison of the open-water performance curves of propeller

從以上螺旋槳模型敞水性能曲線對(duì)比計(jì)算結(jié)果可以看出，基于Laplace 權(quán)函數(shù)插值方法開發(fā)的CFD 求解器，可用于螺旋槳模型敞水性能的數(shù)值計(jì)算。同時(shí)也驗(yàn)證了本文研究給出的螺旋槳敞水性能計(jì)算的推薦方法是可靠的，計(jì)算結(jié)果的精度基本滿足工程應(yīng)用需求。但由于本文測(cè)試的樣本有限，更確定性的結(jié)論還需要開展更廣泛的研究與測(cè)試。

4 結(jié) 論

本文基于中國(guó)船舶科學(xué)研究中心自主研發(fā)的CFD 流場(chǎng)求解器——NaViiX，開展了螺旋槳敞水性能數(shù)值模擬研究。以典型四葉螺旋槳為對(duì)象，研究了計(jì)算方法、網(wǎng)格劃分形式以及交界面位置對(duì)螺旋槳敞水性能計(jì)算結(jié)果的影響，并開展了相應(yīng)測(cè)試，得出如下結(jié)論：

（1）對(duì)于螺旋槳敞水這類只存在周期性旋轉(zhuǎn)部件的繞流問題，研究發(fā)現(xiàn)三種計(jì)算方法結(jié)合三種網(wǎng)格劃分形式的計(jì)算結(jié)果基本一致，多參考坐標(biāo)系方法耗時(shí)最短，非慣性坐標(biāo)系方法次之，滑移網(wǎng)格方法耗時(shí)遠(yuǎn)大于前述兩種方法。相同網(wǎng)格數(shù)量前提下，多參考坐標(biāo)系方法采用內(nèi)域非結(jié)構(gòu)化網(wǎng)格、外域結(jié)構(gòu)化網(wǎng)格的網(wǎng)格劃分形式時(shí)，網(wǎng)格質(zhì)量更好，計(jì)算結(jié)果更優(yōu)。

（2）交界面位置會(huì)對(duì)螺旋槳敞水性能計(jì)算結(jié)果產(chǎn)生影響，研究發(fā)現(xiàn)交界面距離槳盤面中等位置（即交界面前端距離槳盤面距離是0.4D（D為槳盤面直徑），交界面后端距離槳盤面距離是0.2D，交界面周向距離槳軸中心0.7D）時(shí)，計(jì)算結(jié)果與試驗(yàn)結(jié)果符合更好。

（3）綜合考慮計(jì)算精度、計(jì)算周期和網(wǎng)格質(zhì)量之后，針對(duì)當(dāng)前CFD 流場(chǎng)求解，初步給出螺旋槳敞水性能數(shù)值計(jì)算的推薦方案：采用多參考坐標(biāo)系方法，結(jié)合內(nèi)域非結(jié)構(gòu)化網(wǎng)格、外域結(jié)構(gòu)化網(wǎng)格的計(jì)算網(wǎng)格劃分策略，交界面距離槳盤面取中等距離位置。初步的應(yīng)用測(cè)試驗(yàn)證了上述推薦計(jì)算方案的可靠性。

通過本文的研究工作可見，自主研發(fā)的CFD 求解器通過在交界面處采用Laplace權(quán)函數(shù)插值方法進(jìn)行流場(chǎng)信息傳遞，使其具備處理交界面兩側(cè)為不同網(wǎng)格類型的能力，同時(shí)能夠解決交界面離模型很近時(shí)，交界面附近流場(chǎng)劇烈變化導(dǎo)致的計(jì)算不穩(wěn)定問題，可以服務(wù)于螺旋槳模型敞水性能的研究和預(yù)報(bào)。同時(shí)，也可為后續(xù)船艇自航CFD計(jì)算等APP的開發(fā)提供參考和技術(shù)支撐。