張玉倩，曾劍鋒，張弘強，張星燁，王志純

(1.天津天科工程監(jiān)理咨詢事務(wù)所，天津 300452；2.南開大學(xué) 軟件學(xué)院，天津 300071；3.天津大學(xué) 水利工程仿真與安全國家重點實驗室，天津 300072；4.天津大學(xué) 建筑工程學(xué)院，天津 300072；5.法國達(dá)飛海運集團(tuán)公司，馬賽 13000；6.交通運輸部天津水運工程科學(xué)研究所，天津 300456)

隨著全球航運行業(yè)的不斷發(fā)展，船舶螺旋槳射流也愈發(fā)受到國內(nèi)外研究人員的關(guān)注。無論是分析船槳、船舵、船體之間的相互干擾，或是探討附加裝置的節(jié)能效果，都要對螺旋槳射流流場作深入研究。

國內(nèi)外專家學(xué)者關(guān)于螺旋槳射流的研究對象，大部分都是固定于水體中的螺旋槳(以下簡稱為定槳)。其中，Albertson等[1]基于軸向動量理論，采用平面自由射流理論研究螺旋槳射流，得到了平面射流初速度公式，為之后的研究奠定了理論基礎(chǔ)。之后，Blaauw與Kaa[2]、Berger與Cederwall[3]、Verhey[4]、Fuehrer與Romisch[5]在Albertson的理論基礎(chǔ)上，通過建立螺旋槳物理模型，對Albertson理論作了修正，提出了一系列推求射流流場各位置流速的計算方法。Stewart[6]、Hamill與Johnston[7]、Hashmi[8]、Lam等[9]對螺旋槳射流流場的起始流速及流速衰減情況作了進(jìn)一步研究，通過物模試驗對Albertson等經(jīng)驗公式作了修正，得到了更接近于真實情況的螺旋槳射流流速分布規(guī)律。

然而，在實船航行過程中，螺旋槳是推動船體并隨之一起運動的，對于某些特定條件下的定槳物模試驗，其對實船螺旋槳射流真實情況的反映會有所欠缺。溫春鵬等[10]在分析適航資源分布的同時，發(fā)現(xiàn)船舶在航行過程中會對浮泥產(chǎn)生擾動，尾流有較為明顯的渾濁帶。張定軍等[11]通過船舶操縱模擬試驗，研究了船舶港內(nèi)作業(yè)對港池水域的要求。目前，預(yù)報螺旋槳水動力特性的計算方法主要有升力線法、升力面法、面元法、CFD(Computational Fluid Dynamics)方法等。雖然CFD方法起步較晚，但是其適用范圍更廣、計算精度更高。因此，本文在前人的定槳物模試驗基礎(chǔ)上，采用CFD軟件FLUENT對推進(jìn)過程中的螺旋槳(以下簡稱為動槳)射流進(jìn)行三維數(shù)值計算，以得到更符合實際情況的動槳射流流場分布規(guī)律。同時，采用自航船模進(jìn)行動槳射流試驗，通過高精度的PIV(Particle Image Velocimetry)系統(tǒng)進(jìn)行流場測量，以驗證數(shù)值計算方法的可行性。

常規(guī)CFD計算大多涉及的是穩(wěn)態(tài)問題，然而對于旋轉(zhuǎn)、推進(jìn)同時進(jìn)行的螺旋槳，其必然存在運動邊界的問題。因此，本文采用動網(wǎng)格技術(shù)，對存在運動邊界的非定常流動進(jìn)行數(shù)值計算。另外，在螺旋槳射流流場中，射流軸向速度占據(jù)主要地位[12]，因此本文所提及的射流均表示軸向射流。

1 計算對象

1.1 計算域

圖1 計算域Fig.1 Calculation domain

由于實船螺旋槳的外觀、尺寸及空間形態(tài)極其復(fù)雜，且本文研究對象為螺旋槳槳后射流流場，并未涉及螺旋槳推力、扭矩、表面壓力分布等水動力特性。因此，在計算域的建立過程中對螺旋槳進(jìn)行概化，采用邊界類型中的Fan模型替代直徑(以下簡稱為槳徑，用D表示)0.015 m、轉(zhuǎn)速3 000 rpm的螺旋槳。計算域如圖1所示，Gx、Gy、Gz軸的交點為三維坐標(biāo)原點，螺旋槳中軸線(以下簡稱為槳軸線)坐標(biāo)為(x，0.075，-0.015)，坐標(biāo)軸定義為：x軸為縱向，與槳軸線一致，指向下游；y軸為橫向，與槳平面一致；z軸為垂向，服從右手系。

若螺旋槳射流在發(fā)展過程中過早碰到水槽邊壁，則射流流速將會發(fā)生反射，從而影響整個射流流場的計算結(jié)果。為減小邊壁反射對射流流場的不利影響，并綜合計算機的硬件配置及計算精度，計算域采用長0.9 m、寬0.15 m、高0.07 m的長方體區(qū)域。其中，z=-0.045～0的下部區(qū)域為水體部分(即研究區(qū)域)，z=0～0.025的上部區(qū)域為空氣部分。

1.2 計算網(wǎng)格

圖2 計算網(wǎng)格Fig.2 Computational grids

計算網(wǎng)格劃分是數(shù)值計算過程中最為耗時的環(huán)節(jié)之一，但也是直接影響計算精度與效率的關(guān)鍵因素之一。計算網(wǎng)格如圖2所示，本文采用局部加密方法對網(wǎng)格進(jìn)行劃分，不但可以提高計算精度，而且可以避免流場變化平緩區(qū)域的計算資源浪費。

在緊鄰螺旋槳的區(qū)域內(nèi)建立一個長0.9 m、寬0.017 m、高0.017 m的長方體，將螺旋槳的運動軌跡完全包裹在內(nèi)，并采用適應(yīng)性較好的四面體網(wǎng)格對該長方體及其相鄰區(qū)域進(jìn)行加密。在向外平緩過渡的過程中，網(wǎng)格尺寸逐漸增大、結(jié)構(gòu)逐漸稀疏，并采用計算較易收斂的六面體網(wǎng)格。完成計算網(wǎng)格劃分，網(wǎng)格總數(shù)近90萬個。

2 計算方法

2.1 控制方程

對于計算域內(nèi)形式復(fù)雜的三維非定常不可壓縮湍流流動，本文采用Reynolds平均法，即對非穩(wěn)態(tài)Navier-Stokes(N-S)方程采用時間平均法，得到時均形式的控制方程。

目前，Reynolds平均法是使用最為廣泛的湍流計算方法，其核心思路是不直接求解瞬時N-S方程，而是求解時均化Reynolds方程[13]。實際工程關(guān)注的重點是湍流所引起的平均流場變化與流場整體分布，正如本文著重研究槳后射流流場。

2.2 動網(wǎng)格技術(shù)

動網(wǎng)格技術(shù)是用于處理運動邊界所引起的非定常流動的常用方法，該方法在各個研究領(lǐng)域內(nèi)都能夠得到較為廣泛的應(yīng)用[14]。

在數(shù)值計算過程中，為了使網(wǎng)格能夠適應(yīng)運動邊界移動所引起的變化，需要對計算網(wǎng)格進(jìn)行修正。本文采用彈性光滑與局部重劃相結(jié)合的方法，對計算網(wǎng)格進(jìn)行修正。

2.3 UDF技術(shù)

常規(guī)螺旋槳射流數(shù)值計算的研究對象一般也是定槳，需要通過改變來流大小來控制進(jìn)速系數(shù)(用J表示)，以替代實際情況下的進(jìn)速系數(shù)變化。

為真實模擬動槳射流，實現(xiàn)真正意義上的“進(jìn)速系數(shù)”，本文采用可被動態(tài)連接到FLUENT求解器中的UDF(User Defined Function)技術(shù)。采用編譯器VC++對UDF進(jìn)行編譯，控制螺旋槳沿x軸正方向分別作J=0.2的勻速直線運動。計算前，對整個過程中的動網(wǎng)格進(jìn)行更新檢查，通過調(diào)整相關(guān)參數(shù)使得網(wǎng)格質(zhì)量滿足計算要求。

2.4 湍流模式與多相流模式

大部分實際工程中，流體運動都處于湍流狀態(tài)，螺旋槳射流也不例外。相對于其它湍流模式，RNGk-ε模型更適合對旋流、射流等進(jìn)行計算[15]。因此，本文采用RNGk-ε模式作模擬。

處理多相流的計算方法，主要有歐拉-拉格朗日法、歐拉-歐拉法。在FLUENT中，有VOF模式、混合模式、歐拉模式等3種歐拉-歐拉多相流模式可供選擇。由于計算內(nèi)容包含流體與固體表面的相互作用，根據(jù)FLUENT對多相流選擇的基本原則，采用VOF模式作模擬。

2.5 邊界條件設(shè)置

在計算域中，入口為Velocity_Inlet，出口為Outflow，螺旋槳與流體的接觸面為Fan，水槽底面、邊壁均為Wall，動靜交接面?zhèn)鬟f方式采用混合面法[16]。

對流項采用二階迎風(fēng)格式進(jìn)行離散，對擴(kuò)散項采用具有二階精度的中心差分方式進(jìn)行離散。計算初期，分別選取0.002 s、0.01 s兩個時間步長進(jìn)行計算。計算結(jié)果顯示，以上兩個時間步長更新后的網(wǎng)格質(zhì)量均能滿足計算精度的要求，且計算結(jié)果較為接近，說明采用0.01 s的時間步長可以得到與時間步長無關(guān)的數(shù)值解。因此，本文選取0.01 s作為時間步長進(jìn)行計算，每一時間步長內(nèi)迭代50次。

圖3 縱截面射流流速等值圖Fig.3 Jet velocity contours in 4 longitudinal sections

3 計算結(jié)果分析

3.1 進(jìn)速系數(shù)為0.2時的流場分析

3.1.1 縱截面流場分析

計算完成之后，對進(jìn)速系數(shù)J=0.2時的射流流場進(jìn)行分析。從y=0.075的xoz平面(即包含槳軸線的xoz平面，以下簡稱為中軸面)開始，依次提取y=0.075+1/4D=0.078 75、y=0.075+1/2D=0.082 5、y=0.075+3/4D=0.086 25(即槳側(cè)距離槳軸線1/4、1/2、3/4倍槳徑)等3個xoz平面。縱截面射流流速等值圖如圖3所示，通過這4幅縱截面圖，對不同橫向位置的二維射流流場進(jìn)行分析。

由圖3-a可得，射流流場以槳軸線為中心呈對稱分布，流速分布呈雙峰型，斷面最大射流流速位于槳葉中部。在遠(yuǎn)離槳平面的擴(kuò)散過程中，射流流速不斷減小、衰減程度逐漸減小。此外，射流流速在槳后極短距離內(nèi)激增，并于0.37D位置達(dá)到最大流速，該計算結(jié)果也與“Hamill發(fā)現(xiàn)最大流速值在0.35D區(qū)域內(nèi)無衰減”[17]這一研究成果相吻合。

由圖3-b可得，縱截面橫向移動1/4D距離，射流流場仍以槳軸線為中心呈對稱分布，斷面最大射流流速位置已由槳葉中部逐漸偏向槳中心，且雙峰型已逐漸向單峰型衰減。由圖3-c可得，縱截面橫向移動1/2D距離，雙峰型已完全衰減為單峰型，斷面最大射流流速位于槳中心。由圖3-d可得，縱截面橫向移動3/4D距離，流速分布呈細(xì)長橢圓型，斷面最大射流流速位于槳中心。

縱觀4幅縱截面射流流速等值圖，隨著縱截面距中軸面橫向距離的增大，射流流速及其衰減程度不斷減小，流場影響范圍逐漸減小，雙峰型特征愈發(fā)減弱，斷面最大射流流速位置已由槳葉中部逐漸集中于槳中心。

3.1.2 中軸面流速變化分析

圖4 中軸面射流流速變化圖Fig.4 Jet velocity graph of center axial plane

對于射流流場最具有代表性、流速最為顯著的中軸面，有必要對其流速變化進(jìn)行深入分析。在中軸面上，沿水深方向依次提取4條垂直于z軸的直線，往深處依次為Z1(x，0.075，-D)、Z2(x，0.075，-1.25D)、Z3(x，0.075，-1.5D)、Z4(x，0.075，-1.75D)，由此作出不同水深位置的射流流速曲線。中軸面射流流速沿水深方向的變化曲線圖如圖4所示，圖中Vx表示螺旋槳射流的軸向速度大小(下同)，x/D表示槳后若干倍槳徑的縱向位置。

由圖4可得，不同水深位置的4條射流流速曲線，在分布規(guī)律上均保持一致：從槳平面開始，射流流速激增，并于槳后0.37D的位置達(dá)到峰值；峰值過后，射流流速迅速減小，衰減程度較大；逐漸的，衰減程度減小，射流流速緩慢減小。并且，槳后越遠(yuǎn)位置，射流流速衰減程度越小。

圖5 槳后三維射流流線圖Fig.5 3D diagram of jet streamline behind propeller

由于斷面最大射流流速位于槳葉中部，因此，水深位置為槳葉中部的直線Z2，相對于水深位置為槳軸線的直線Z1，射流流速更大，且為相同x軸位置的流速最大值。并且，隨著水深不斷增大，射流流速衰減程度逐漸減小。

3.1.3 三維流場分析

為了能夠更加直觀的對槳后射流流場進(jìn)行分析，提取槳后及其相鄰區(qū)域內(nèi)的三維射流流線圖，對其總體分布趨勢進(jìn)行分析，槳后三維射流流線圖如圖5所示(圖例與圖3-e相同)。

由圖5可得，射流流線以槳軸線為中軸線，圍繞其向后呈螺旋式發(fā)展。槳平面后的射流流速均高于外流域，這也是螺旋槳向后推水、水體反作用于螺旋槳而產(chǎn)生推力的結(jié)果。較大的流速變化均高度集中于槳平面及槳后流域，隨著遠(yuǎn)離槳平面，外流域的變化與影響程度迅速減小。內(nèi)半徑處，射流在螺旋槳作用下向內(nèi)收縮，直徑也逐漸減小，反映了螺旋槳對水體的抽吸作用。此外，該圖也與Lam等[15]用CFD方法模擬出的螺旋槳射流流線圖相吻合。

3.2 模型驗證

本文采用自航船模進(jìn)行動槳射流試驗，并采用高精度的PIV系統(tǒng)進(jìn)行流場測量，以驗證數(shù)值仿真方法的可行性。本文采用由美國TSI公司生產(chǎn)的立體PIV流場測量系統(tǒng)，該系統(tǒng)由激光器、同步器、光電編碼器、CCD攝像機等精密部件組成。采用的自航船模由國內(nèi)船模試驗研究經(jīng)驗非常豐富的交通運輸部天津水運工程科學(xué)研究所研制，并在吃水深度、靜水航速、零舵角直航操縱性、“Z”型操縱性等方面與實船進(jìn)行了相似性率定。

圖6 物模試驗布置Fig.6 Layout of physical model test

在長20.0 m、寬0.4 m、高0.5 m的試驗水槽中，水深0.1 m，靜水；自航船模螺旋槳槳徑0.015 m，轉(zhuǎn)速3 000 rpm，槳中心位于水面以下0.015 m位置，螺旋槳推動船體以進(jìn)速系數(shù)0.2勻速直線通過水槽中部試驗段(即PIV系統(tǒng)可拍攝段)，物模試驗布置如圖6所示。在試驗數(shù)據(jù)中提取中軸面上位于槳后0.015 m、0.030 m、0.060 m、0.150 m(即槳后1、2、4、10倍槳徑)等4條垂直于x軸的直線，得到槳后不同距離射流沿水深方向的流速分布。相應(yīng)的，在數(shù)值仿真結(jié)果中提取中軸面上位于槳后0.015 m、0.030 m、0.060 m、0.150 m等4條垂直于x軸的直線，與物模試驗所得到的槳后不同距離射流沿水深方向的流速分布進(jìn)行對比驗證。

物模試驗與數(shù)值計算轉(zhuǎn)速(n)比尺為1:1、槳徑(d)比尺為1:1，且射流流速V～nd[18]，因此射流流速(V)比尺為1:1。物模試驗與數(shù)值計算所得到的射流流速分布比較如圖7所示，圖中z/D表示水面以下若干倍槳徑的垂向位置。

7-a 槳后1倍槳徑位置射流流速分布7-b 槳后2倍槳徑位置射流流速分布7-c 槳后4倍槳徑位置射流流速分布7-d 槳后10倍槳徑位置射流流速分布圖7 射流流速分布比較Fig.7 Comparisons of jet velocity distributions

由圖7-a、7-b可得，相對于試驗結(jié)果，模擬結(jié)果中流速沿程衰減稍快、斷面流速最大值稍大，但兩者總體趨勢吻合較好。由圖7-c可得，相對于試驗結(jié)果，模擬結(jié)果中流速沿程衰減稍慢、斷面流速最大值稍大，兩者總體趨勢依舊吻合較好。由圖7-d可得，相對于試驗結(jié)果，模擬結(jié)果中僅流速沿程衰減稍慢，其余兩者均吻合較好。

縱觀4幅射流流速分布對比圖，達(dá)到第一個流速峰值之前，兩條曲線的斜率基本一致，流速及其衰減程度均基本相同；經(jīng)過第一個流速峰值之后大約0.6D的垂向深度范圍內(nèi)，數(shù)值計算所得到的流速略大于物模試驗結(jié)果，但偏差均在5%以內(nèi)，且二者流速衰減程度基本相同；再往下的垂向范圍內(nèi)，兩條曲線的斜率又基本趨于一致，流速及其衰減程度均基本相同。

數(shù)值計算與物模試驗中的射流流速分布總體趨勢保持一致，最大射流流速近似相等，射流流速衰減程度基本相同。至于二者在部分范圍內(nèi)的射流流速、最大射流流速出現(xiàn)位置、軸對稱情況等方面還存在些許差異，主要有兩個方面的原因：一方面，物模試驗過程中會出現(xiàn)偶然性誤差，因此會有部分?jǐn)?shù)據(jù)偏離整體趨勢，甚至不符合實際規(guī)律；另一方面，數(shù)值計算嚴(yán)格按照選定公式進(jìn)行迭代計算，因此結(jié)果較物模試驗更為理想化，數(shù)值亦會偏大。

對于水流流向復(fù)雜、流速多變的螺旋槳射流流場來說，以上些許差異，均可允許。綜上所述，數(shù)值計算結(jié)果與物模試驗結(jié)果基本吻合，驗證了數(shù)值計算方法的可行性。

4 結(jié)論

實船航行過程中，由于螺旋槳運動使得槳平面、槳后流場不斷發(fā)生變化，常規(guī)CFD穩(wěn)態(tài)分析方法的計算精度較低，不能真實有效地模擬動槳射流情況。因此，本文采用UDF及動網(wǎng)格技術(shù)，對動槳射流進(jìn)行三維數(shù)值計算，得出的主要結(jié)論如下：

(1)射流中軸面上，流場以槳軸線為中心呈對稱分布，流速分布呈雙峰型，且斷面最大流速位于槳葉中部。射流流速在槳后極短距離內(nèi)激增，并于0.37D位置達(dá)到最大流速。

(2)偏離槳軸線的橫向距離越大，射流流速及其衰減程度越小，流場影響范圍越小，雙峰型特征愈發(fā)減弱，斷面最大流速由槳葉中部逐漸集中于槳中心。偏離槳軸線的垂向距離增大，射流流速先增后減，衰減程度逐漸減小。

(3)采用自航船模進(jìn)行動槳射流試驗，并采用高精度的PIV系統(tǒng)進(jìn)行流場測量，驗證了數(shù)值計算方法的可行性，說明本文提出的方法適用于動槳射流研究，并能為較復(fù)雜工況下的動槳射流研究提供一定的理論依據(jù)。

致謝：參與本研究工作的還有長沙理工大學(xué)水利工程學(xué)院胡旭躍教授、沈小雄教授，在此并致謝意。