動槳射流數(shù)值模擬

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

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

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

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

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

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

1 計算對象

1.1 計算域

圖1 計算域Fig.1 Calculation domain

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

若螺旋槳射流在發(fā)展過程中過早碰到水槽邊壁，則射流流速將會發(fā)生反射，從而影響整個射流流場的計算結果。為減小邊壁反射對射流流場的不利影響，并綜合計算機的硬件配置及計算精度，計算域采用長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é)之一，但也是直接影響計算精度與效率的關鍵因素之一。計算網(wǎng)格如圖2所示，本文采用局部加密方法對網(wǎng)格進行劃分，不但可以提高計算精度，而且可以避免流場變化平緩區(qū)域的計算資源浪費。

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

2 計算方法

2.1 控制方程

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

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

2.2 動網(wǎng)格技術

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

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

2.3 UDF技術

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

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

2.4 湍流模式與多相流模式

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

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

2.5 邊界條件設置

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

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

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

3 計算結果分析

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

3.1.1 縱截面流場分析

計算完成之后，對進速系數(shù)J=0.2時的射流流場進行分析。從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(即槳側距離槳軸線1/4、1/2、3/4倍槳徑)等3個xoz平面?？v截面射流流速等值圖如圖3所示，通過這4幅縱截面圖，對不同橫向位置的二維射流流場進行分析。

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

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

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

3.1.2 中軸面流速變化分析

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

對于射流流場最具有代表性、流速最為顯著的中軸面，有必要對其流速變化進行深入分析。在中軸面上，沿水深方向依次提取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的位置達到峰值；峰值過后，射流流速迅速減小，衰減程度較大；逐漸的，衰減程度減小，射流流速緩慢減小。并且，槳后越遠位置，射流流速衰減程度越小。

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

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

3.1.3 三維流場分析

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

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

3.2 模型驗證

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

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

在長20.0 m、寬0.4 m、高0.5 m的試驗水槽中，水深0.1 m，靜水；自航船模螺旋槳槳徑0.015 m，轉速3 000 rpm，槳中心位于水面以下0.015 m位置，螺旋槳推動船體以進速系數(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軸的直線，得到槳后不同距離射流沿水深方向的流速分布。相應的，在數(shù)值仿真結果中提取中軸面上位于槳后0.015 m、0.030 m、0.060 m、0.150 m等4條垂直于x軸的直線，與物模試驗所得到的槳后不同距離射流沿水深方向的流速分布進行對比驗證。

物模試驗與數(shù)值計算轉速(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可得，相對于試驗結果，模擬結果中流速沿程衰減稍快、斷面流速最大值稍大，但兩者總體趨勢吻合較好。由圖7-c可得，相對于試驗結果，模擬結果中流速沿程衰減稍慢、斷面流速最大值稍大，兩者總體趨勢依舊吻合較好。由圖7-d可得，相對于試驗結果，模擬結果中僅流速沿程衰減稍慢，其余兩者均吻合較好。

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

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

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

4 結論

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

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

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

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

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