噴水推進(jìn)器進(jìn)流面獲取方法及其應(yīng)用

郭 軍， 陳作鋼， 戴原星， 陳建平

（1.上海交通大學(xué) 海洋工程國家重點(diǎn)實(shí)驗(yàn)室，上海200240；2.上海交通大學(xué) 高新船舶與深海開發(fā)裝備協(xié)同創(chuàng)新中心，上海200240；3.上海交通大學(xué) 船舶海洋與建筑工程學(xué)院，上海200240；4.中國船舶及海洋工程設(shè)計(jì)研究院，上海200011；5.噴水推進(jìn)技術(shù)重點(diǎn)實(shí)驗(yàn)室，上海200011）

噴水推進(jìn)器是利用噴出水流的反作用力來推動船舶前進(jìn)的一種特殊推進(jìn)器，具有高航速時推進(jìn)效率高、操縱性好、抗空化、低噪聲等優(yōu)勢，在高性能船舶上得到越來越多的應(yīng)用.隨著計(jì)算流體動力學(xué)（CFD）的飛速發(fā)展，CFD技術(shù)在船舶領(lǐng)域發(fā)揮著越來越重要的作用.在噴水推進(jìn)器方面，CFD技術(shù)廣泛運(yùn)用于性能研究和優(yōu)化設(shè)計(jì)中，現(xiàn)已取得相當(dāng)進(jìn)展.在國外，Bulten[1]對噴水推進(jìn)器進(jìn)行了數(shù)值研究，分析了定常多重參考系（MRF）模型和瞬時剛體運(yùn)動模型對噴水推進(jìn)器性能的影響.Takai等[2]對噴水推進(jìn)船流道形狀的優(yōu)化進(jìn)行了研究，采用體積力法來模擬泵的效果，研究表明對流道上部曲線和唇部形狀進(jìn)行優(yōu)化可以減小壓力損失，對流道進(jìn)口進(jìn)行優(yōu)化可以提高入口效率.Peri等[3]采用優(yōu)化軟件對高速噴水推進(jìn)雙體船進(jìn)行了優(yōu)化設(shè)計(jì)，以減小入口損失和減小阻力為目標(biāo)，分別對噴水推進(jìn)流道形狀和船體型線進(jìn)行了優(yōu)化，最終優(yōu)化后阻力減小了6%.Guo等[4]對噴水推進(jìn)器的進(jìn)口流動進(jìn)行了研究，分析了不同進(jìn)速比下進(jìn)口的壓力分布，研究表明進(jìn)速比越大進(jìn)口壓力損失越大.Eslamdoost等[5－6]提出了一種壓力跳躍的方法來代替泵的作用，對噴水推進(jìn)器與船體的相互作用進(jìn)行了研究，分析了噴水推進(jìn)負(fù)推力減額的原因.Lu等[7]對噴水推進(jìn)泵進(jìn)行了數(shù)值計(jì)算，研究了噴水推進(jìn)泵的水動力性能.在國內(nèi)，孫存樓等[8]進(jìn)行了噴水推進(jìn)船負(fù)推力減額機(jī)理的研究，認(rèn)為中高速時黏性阻力降低是推力減額的主要原因.龔杰等[9]進(jìn)行了噴水推進(jìn)船的自航數(shù)值模擬，比較了虛擬盤法和重疊網(wǎng)格法對船內(nèi)外流場的影響，研究表明兩種方法的外流場計(jì)算結(jié)果基本相同，內(nèi)流場采用重疊網(wǎng)格法更接近真實(shí)葉輪情況.易文彬等[10]采用試驗(yàn)和數(shù)值模擬相結(jié)合的方法分析了浸沒式噴水推進(jìn)的水動力特點(diǎn)，試驗(yàn)值和CFD預(yù)報值吻合較好.

噴水推進(jìn)器推力性能是噴水推進(jìn)系統(tǒng)性能的一項(xiàng)重要內(nèi)容.試驗(yàn)時噴水推進(jìn)器推力的測量有直接測量法和動量流量法，動量流量法通過進(jìn)流面和噴口的動量差來計(jì)算推力.由于直接測量法所需的測量設(shè)備非常復(fù)雜，技術(shù)門檻較高，目前僅有國外少數(shù)幾家單位能夠直接測量；動量流量法不但可以選取滿足流量要求的任意泵，而且船體與噴水推進(jìn)之間無須復(fù)雜的水密結(jié)構(gòu)，已成為國內(nèi)外研究機(jī)構(gòu)進(jìn)行試驗(yàn)和數(shù)值計(jì)算時采用的主要方法.從第21至第24屆國際拖曳水池會議（ITTC），經(jīng)過噴水推進(jìn)專家委員會的不斷完善，最終推出了基于動量流量法的噴水推進(jìn)模型試驗(yàn)的推薦指南[11]，形成了比較嚴(yán)謹(jǐn)?shù)睦碚擉w系.

根據(jù)動量流量法，在CFD數(shù)值計(jì)算噴水推進(jìn)器推力的過程中，流量與噴口的速度很容易獲得，但進(jìn)流面的的速度不易獲得，主要是進(jìn)流面的形狀無法確定，因此如何準(zhǔn)確獲得進(jìn)流面的形狀成為準(zhǔn)確計(jì)算噴水推進(jìn)器推力的關(guān)鍵.Ding等[12]提出了自定義標(biāo)量函數(shù)的方法來確定噴水推進(jìn)系統(tǒng)的進(jìn)流面，但此種方法需要用戶在商業(yè)軟件中編寫后處理程序進(jìn)行二次開發(fā)，其要求較高，不適用于大多數(shù)研究人員.本文提出一種基于流線法的進(jìn)流面獲取方法，該方法簡單實(shí)用，適用于大多數(shù)的研究人員.使用MATLAB語言開發(fā)成程序，使過程大大簡化.在此基礎(chǔ)上，采用CFD軟件STAR－CCM＋對噴水推進(jìn)器的性能進(jìn)行數(shù)值計(jì)算.首先進(jìn)行了網(wǎng)格無關(guān)性驗(yàn)證，之后選取合適的網(wǎng)格，系統(tǒng)研究了噴水推進(jìn)器在不同轉(zhuǎn)數(shù)下的性能，為噴水推進(jìn)器系統(tǒng)性能的提高和改善奠定了基礎(chǔ).

1 理論模型

1.1 基本方程

流體力學(xué)問題要遵守基本守恒定律，包括質(zhì)量守恒、動量守恒和能量守恒.在船舶流體力學(xué)問題中，通常不考慮能量交換，故忽略能量守恒，則不可壓縮流體對應(yīng)的連續(xù)性方程和動量守恒方程分別為

式中：ui和uj為速度時均量；xi和xj為坐標(biāo)；t為時間；ρ為流體密度；p為壓力時均量；ν為流體運(yùn)動黏性系數(shù)為雷諾應(yīng)力.

由于式（2）是不封閉的，需引入湍流模型，所以本文計(jì)算中采用SSTk－ω湍流模型.該模型結(jié)合k－ω和k－ε模型，在預(yù)報具有逆壓梯度的流動和流動分離等方面具有優(yōu)勢，其表達(dá)形式為

式（3）和（4）中各參數(shù)的具體定義以及含義可參見文獻(xiàn)[13].

1.2 動量流量法

根據(jù)ITTC規(guī)程，動量流量法求解噴水推進(jìn)器推力的流場控制體如圖1所示.

圖1 流場控制體Fig.1 Control volume of waterjet

進(jìn)流面A1與來流方向垂直，位于流道斜坡與船體的切點(diǎn)A′之前，距離通常取1倍的流道直徑，以避免流道入口的影響.通常認(rèn)為收縮斷面A7與噴口面積A6的直徑相同，將噴口面積A6作為控制體的出口面.流場控制體在x方向的總推力為

流體作用在噴水推進(jìn)器壁面的力，直接傳遞在船體上，稱之為凈推力.在計(jì)算過程中可以很方便地監(jiān)測噴水推進(jìn)器各壁面上的壓力和黏性力，通過壁面積分求解出噴水推進(jìn)器的凈推力，即噴水推進(jìn)器在x方向的合力為

式中：S為噴水推進(jìn)器的所有壁面面積；px為壁面上的壓應(yīng)力；σx為壁面上的剪切應(yīng)力.

本文基于流線法獲取進(jìn)流面的基本思想如下：

（1）泵進(jìn)口面A3設(shè)置流線起點(diǎn)，沿上游計(jì)算流線，流線分布區(qū)域的半徑應(yīng)略大于流道半徑，適當(dāng)設(shè)置流線種子的數(shù)量以使流線基本布滿控制體，這樣計(jì)算的形狀更加接近于進(jìn)流面的形狀.

（2）在進(jìn)流面的位置A1處建立與來流方向相垂直的平面截面，平面截面的對象為流線，如圖2所示.

（3）對進(jìn)流面上的流線截面數(shù)據(jù)（圖3中藍(lán)點(diǎn)）進(jìn)行平面點(diǎn)集的凸包計(jì)算，從而得到凸包點(diǎn)集最外層的邊界點(diǎn)，然后采用孔為平等[14]提出的4次多項(xiàng)式對邊界點(diǎn)進(jìn)行擬合，最終得到進(jìn)流面的形狀曲線（圖3中紅線）：

式中：c0，c1，c2，c3，c4分別為4次多項(xiàng)式擬合系數(shù).

采用MATLAB語言開發(fā)的進(jìn)流面獲取程序如圖3所示，進(jìn)流面獲取的原則是進(jìn)流面和噴口截面的流量相等.程序有6個參數(shù)：進(jìn)流面位置船尾的斜升角；進(jìn)流面形狀中心的y坐標(biāo)和z坐標(biāo)；邊界曲線的延伸系數(shù)，以確保進(jìn)流面曲線的左右邊界超過控制體；縮放系數(shù)，由于流線法確定的進(jìn)流面完全在控制體內(nèi)部，進(jìn)流面的流量會略小于噴口流量，而縮放系數(shù)可以確保進(jìn)流面的流量與噴口的流量相等.

圖2 流線法Fig.2 Streamline method

圖3 進(jìn)流面獲取程序Fig.3 Program of capture area

當(dāng)從泵進(jìn)口面A3處計(jì)算的流線足夠多時，流線所組成的包絡(luò)體即為流場控制體，流線在位置A1的截面即為理論上的進(jìn)流面.進(jìn)流面形狀對流線數(shù)量（m）的依賴性見圖4，以流道直徑D對進(jìn)流面形狀進(jìn)行無量綱處理，Y＝y(tǒng)／D表示無量綱的y坐標(biāo)，Z＝z／D表示無量綱的z坐標(biāo).可見：當(dāng)m＝500時，得到的進(jìn)流面形狀有微小偏差；當(dāng)m＞1 000后，得到的進(jìn)流面形狀不再變化，說明當(dāng)流線數(shù)大于1 000后，本文方法可準(zhǔn)確得到進(jìn)流面的形狀，也可說明本文方法的合理性.

圖4 進(jìn)流面對流線的依賴性Fig.4 Dependence of capture area on streamlines

2 計(jì)算模型

計(jì)算域長度為30D，流道入口前為25D，流道入口后為5D，寬度為10D，深度為8D.

網(wǎng)格使用STAR－CCM＋軟件生成，網(wǎng)格類型為Trim非結(jié)構(gòu)網(wǎng)格.在葉輪和導(dǎo)葉的導(dǎo)邊和隨邊處生成特征線，進(jìn)行局部加密.在流道入口處流場變化劇烈，設(shè)置網(wǎng)格控制體進(jìn)行局部加密.計(jì)算域總的網(wǎng)格數(shù)約為291萬個，葉輪旋轉(zhuǎn)區(qū)域約為94萬個，其余靜止區(qū)域約為197萬個.計(jì)算域的網(wǎng)格見圖5.

圖6所示為噴水推進(jìn)器進(jìn)口流道和葉輪表面第一層網(wǎng)格節(jié)點(diǎn)與壁面的無量綱距離y＋的分布，噴水推進(jìn)器流道的y＋為90～160，葉輪和導(dǎo)葉的y＋在90左右，整體上噴水推進(jìn)器的y＋大部分在30～300.第一層網(wǎng)格節(jié)點(diǎn)布置在對數(shù)律成立的范圍內(nèi)，以滿足壁面函數(shù)的要求.

圖5 計(jì)算域的網(wǎng)格Fig.5 Grids of computational domain

圖6 y＋分布云圖Fig.6 y＋contour

計(jì)算域入口、兩側(cè)和底部設(shè)為速度進(jìn)口；計(jì)算域出口和推進(jìn)器噴口為壓力出口；計(jì)算域頂部、流道壁面和導(dǎo)葉設(shè)為無滑移壁面；將計(jì)算域分為旋轉(zhuǎn)區(qū)域和靜止區(qū)域，旋轉(zhuǎn)區(qū)域采用穩(wěn)態(tài)多參考系（MFR）方法模擬泵的旋轉(zhuǎn)，在旋轉(zhuǎn)區(qū)域的速度根據(jù)旋轉(zhuǎn)角速度給出，葉輪和導(dǎo)流帽設(shè)為靜止壁面；旋轉(zhuǎn)區(qū)域和靜止區(qū)域的交界設(shè)為交界面.

考慮船體邊界層的影響，入口的速度分布采用平板邊界層速度分布：式中：v0為設(shè)計(jì)航速；z為距計(jì)算域頂部邊界的距離；δ為船模邊界層厚度；L為計(jì)算域入口至船首的距離；Re為雷諾數(shù).

采用基于分離流的黏性求解器，壓力與速度的耦合采用SIMPLE算法計(jì)算，對流項(xiàng)離散格式為二階迎風(fēng)格式.

3 計(jì)算結(jié)果及分析

對噴水推進(jìn)器中的各水動力系數(shù)定義如下：

式中：CTg和CTn分別為噴水推進(jìn)器的總推力系數(shù)和凈推力系數(shù)；KQJ、KH、KQ和η分別為噴水推進(jìn)器的流量系數(shù)、揚(yáng)程系數(shù)、轉(zhuǎn)矩系數(shù)和效率；QJ、H和Qshaft分別為流量、揚(yáng)程和轉(zhuǎn)軸上的轉(zhuǎn)矩；n為轉(zhuǎn)速；RIV和RNV分別為噴水推進(jìn)器的進(jìn)速比和噴速比；v3為泵前進(jìn)口x方向的平均速度；v6為噴口x方向的平均速度；α和β為分別為噴水推進(jìn)器進(jìn)流面的動量影響系數(shù)和能量影響系數(shù)[15]；vx為x方向的速度.

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

噴水推進(jìn)器葉輪、導(dǎo)葉和流道以及加密區(qū)域的網(wǎng)格大小都與基本尺寸相關(guān)，由基本尺寸的變化生成由密到疏的5套網(wǎng)格，分別編號為1～5，來進(jìn)行網(wǎng)格無關(guān)性的分析，5套網(wǎng)格的數(shù)目見表1.

在設(shè)計(jì)轉(zhuǎn)速n0下，進(jìn)行網(wǎng)格無關(guān)性分析，將最密網(wǎng)格（1號網(wǎng)格）的推力和轉(zhuǎn)矩計(jì)算結(jié)果分別記為CTn1和Qshaft1，其他網(wǎng)格下推力和轉(zhuǎn)矩分別以該值為基準(zhǔn)進(jìn)行歸一化處理，網(wǎng)格無關(guān)性計(jì)算結(jié)果見圖7.

表1 網(wǎng)格的數(shù)目與尺寸Tab.1 Grid sizes and numbers of five grids

圖7 網(wǎng)格數(shù)目依賴性Fig.7 Dependence of grid number

由圖7可以看出，轉(zhuǎn)矩對網(wǎng)格的依賴性不明顯，推力對網(wǎng)格的依賴性較明顯.當(dāng)網(wǎng)格數(shù)增加至2.91×106時，推力和轉(zhuǎn)矩都基本接近各自的參考基準(zhǔn)值.綜合考慮計(jì)算精度和計(jì)算資源，最終采用2號網(wǎng)格方案來進(jìn)行噴水推進(jìn)器性能的數(shù)值計(jì)算.

3.2 噴水推進(jìn)器性能分析

在設(shè)計(jì)航速v0下，研究噴水推進(jìn)器在不同轉(zhuǎn)速（n／n0＝0.4，0.6，0.8，1.0，1.2，1.4，1.6，1.8，2.0）時的水動力性能.

噴水推進(jìn)器的進(jìn)流面動量影響系數(shù)α、進(jìn)流面能量影響系數(shù)β、進(jìn)速比RIV和噴速比RNV是影響噴水推進(jìn)器性能的重要參數(shù)；α、RIV和RNV分別為表征A1、A3和A6平均速度的無量綱量.在不同轉(zhuǎn)速下，α、β、RIV和RNV的結(jié)果見表2.

表2 影響噴水推進(jìn)器性能的參數(shù)Tab.2 Parameters affecting the performance of waterjet

根據(jù)本文開發(fā)的程序可以很方便地得到進(jìn)流面邊界的曲線形狀，圖8所示為不同轉(zhuǎn)速下進(jìn)流面邊界的曲線形狀.由于進(jìn)流面形狀左右對稱，參數(shù)的奇數(shù)項(xiàng)系數(shù)為0，表3所示為進(jìn)流面形狀在不同轉(zhuǎn)速下的相關(guān)參數(shù)化表達(dá).其中：w為進(jìn)流面寬度；c0的絕對值為無量綱的進(jìn)流面高度.在低轉(zhuǎn)速時，進(jìn)流面邊界的曲線形狀不是橢圓形，中間稍微內(nèi)凹，說明在經(jīng)驗(yàn)設(shè)計(jì)中，如果采用橢圓方程計(jì)算的進(jìn)流面和流量會產(chǎn)生較大偏差；隨著轉(zhuǎn)速的增加，中間形狀由內(nèi)凹向直線再向外凸過渡，進(jìn)流面的邊界形狀逐漸趨向橢圓形；隨著轉(zhuǎn)速的增加，進(jìn)流面的高度和寬度都增加，高度的變化大于寬度的變化.在進(jìn)行噴水推進(jìn)試驗(yàn)時，本文計(jì)算的進(jìn)流面形狀參數(shù)可為獲取進(jìn)流面提供依據(jù).

圖8 進(jìn)流面形狀Fig.8 Shapes of capture area

表3 進(jìn)流面的形狀參數(shù)Tab.3 Shape parameters of capture area

圖9所示為不同轉(zhuǎn)速下進(jìn)流面A1的速度分布，v為流體的速度.在低于設(shè)計(jì)轉(zhuǎn)速時，進(jìn)流幾乎全部取自邊界層，隨著轉(zhuǎn)速的增加流量也相應(yīng)增加，進(jìn)流的高度逐漸大于邊界層厚度，進(jìn)流有相當(dāng)一部分取自邊界層以外的流場.

圖10所示為不同轉(zhuǎn)速下泵前進(jìn)口處A3面上速度分布云圖.可以看出：泵前進(jìn)口處速度分布并不均勻，下部速度大，上部速度小；在低轉(zhuǎn)速時，下部速度較均勻，上部速度梯度變化大；在轉(zhuǎn)速較高時，則上部速度較均勻，下部速度梯度變化大；受到泵旋轉(zhuǎn)的影響，速度向旋轉(zhuǎn)方向偏移.

圖11所示為不同轉(zhuǎn)速下噴口處A6面上速度分布云圖.可以看出：噴口處速度受到泵旋轉(zhuǎn)的影響，分布并不均勻；在噴口中心處有速度較低的區(qū)域，低速區(qū)的直徑大約為噴口直徑的20%，低速區(qū)的速度大約為噴口最大速度的50%，在低轉(zhuǎn)速時，低速區(qū)稍向下偏離噴口的幾何中心.

噴水推進(jìn)器唇部的流場變化劇烈，會產(chǎn)生低壓區(qū)，容易發(fā)生空化現(xiàn)象，圖12所示為不同轉(zhuǎn)速下唇部表面壓力分布云圖和流線圖.定義無量綱的壓力系數(shù)為

式中：p為絕對壓力；p0為參考壓力（取1個標(biāo)準(zhǔn)大氣壓）.

由圖12可以發(fā)現(xiàn)，在唇部發(fā)生了流動分流，而發(fā)生流動分流的位置隨著轉(zhuǎn)速的增加向下移動，這是因?yàn)殡S著轉(zhuǎn)速的增加，流量增加，需要進(jìn)口從更深的位置抽吸水流.低壓區(qū)域與流動分流點(diǎn)的變化呈相反趨勢，低轉(zhuǎn)速時，低壓區(qū)域在唇部的下緣；高轉(zhuǎn)速時，低壓區(qū)域在唇部的上緣.

圖9 進(jìn)流面的速度分布云圖Fig.9 Velocity contour of capture area

圖10 泵前進(jìn)口的速度分布云圖Fig.10 Velocity contour of pump intake

圖11 泵前進(jìn)口的速度分布云圖Fig.11 Velocity contour of pump intake

圖13所示為推力系數(shù)曲線，總推力系數(shù)曲線和凈推力系數(shù)曲線基本重合，說明ITTC推薦的動量流量法和壁面積分法都可以進(jìn)行噴水推進(jìn)器的推力預(yù)報，也間接說明本文對進(jìn)流面獲取的準(zhǔn)確性，動量流量計(jì)算的推力精度較高.在噴水推進(jìn)計(jì)算中，如果采用壓力跳越法[5]或其他方法代替泵的作用，壁面積分法就不再適用，只能采用動量流量法來計(jì)算推力，因此動量流量法的適用性更廣，此時采用本文開發(fā)的進(jìn)流面獲取程序?qū)蟠筇岣咝?

圖14所示為噴水推進(jìn)器的性能曲線，噴水推進(jìn)當(dāng)中的流量－揚(yáng)程曲線在物理含義上相當(dāng)于敞水螺旋槳中的進(jìn)速－推力曲線.可以看出：隨著流量系數(shù)KQJ的增加，揚(yáng)程系數(shù)KH和轉(zhuǎn)矩系數(shù)KQ不斷下降；效率η在流量系數(shù)KQJ＝0.87時達(dá)到最大，在最高效率之前效率變化緩慢，此時轉(zhuǎn)速比n／n0＝0.8～2.0，說明噴水推進(jìn)器在較大的轉(zhuǎn)速范圍內(nèi)都有較高的效率；在最高效率之后效率η下降較快，主要原因是低轉(zhuǎn)速時揚(yáng)程系數(shù)KH下降較快.

圖12 唇部的壓力系數(shù)云圖與流線Fig.12 Pressure coefficient contour and streamline in lip

圖13 總推力和凈推力曲線Fig.13 Gross thrust curve and net thrust curve

圖14 噴水推進(jìn)器性能曲線Fig.14 Performance curves of waterjet

4 結(jié)語

本文對噴水推進(jìn)器進(jìn)流面的獲取方法進(jìn)行研究，采用基于流線法的進(jìn)流面獲取方法可以準(zhǔn)確得到進(jìn)流面，該方法操作簡單方便.開發(fā)的進(jìn)流面獲取程序使后處理過程大大簡化，提高了效率.在此基礎(chǔ)上，對噴水推進(jìn)器在不同轉(zhuǎn)速下的性能進(jìn)行了數(shù)值預(yù)報和分析，給出了進(jìn)流面形狀的參數(shù)化表達(dá)，為噴水推進(jìn)試驗(yàn)提供指導(dǎo).在低轉(zhuǎn)速時，進(jìn)流面形狀不是橢圓形，隨著轉(zhuǎn)速增加，進(jìn)流面形狀逐漸趨向橢圓形.進(jìn)流面形狀、泵前進(jìn)口面速度分布以及唇部流動分離點(diǎn)位置都隨轉(zhuǎn)速變化而變化.采用動量流量法計(jì)算噴水推進(jìn)器的推力，不僅精度高，而且應(yīng)用更廣.

致謝 本研究得到了噴水推進(jìn)技術(shù)重點(diǎn)實(shí)驗(yàn)室的支持，得益于噴水推進(jìn)技術(shù)學(xué)術(shù)研討會（WORKSHOP2018）多位與會專家的技術(shù)交流，在此深表謝意.