噴水推進器進流面獲取方法及其應用

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

（1.上海交通大學 海洋工程國家重點實驗室，上海200240；2.上海交通大學 高新船舶與深海開發裝備協同創新中心，上海200240；3.上海交通大學 船舶海洋與建筑工程學院，上海200240；4.中國船舶及海洋工程設計研究院，上海200011；5.噴水推進技術重點實驗室，上海200011）

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

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

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

1 理論模型

1.1 基本方程

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

式中：ui和uj為速度時均量；xi和xj為坐標；t為時間；ρ為流體密度；p為壓力時均量；ν為流體運動黏性系數為雷諾應力.

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

式（3）和（4）中各參數的具體定義以及含義可參見文獻[13].

1.2 動量流量法

根據ITTC規程，動量流量法求解噴水推進器推力的流場控制體如圖1所示.

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

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

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

式中：S為噴水推進器的所有壁面面積；px為壁面上的壓應力；σx為壁面上的剪切應力.

本文基于流線法獲取進流面的基本思想如下：

（1）泵進口面A3設置流線起點，沿上游計算流線，流線分布區域的半徑應略大于流道半徑，適當設置流線種子的數量以使流線基本布滿控制體，這樣計算的形狀更加接近于進流面的形狀.

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

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

式中：c0，c1，c2，c3，c4分別為4次多項式擬合系數.

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

圖2 流線法Fig.2 Streamline method

圖3 進流面獲取程序Fig.3 Program of capture area

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

圖4 進流面對流線的依賴性Fig.4 Dependence of capture area on streamlines

2 計算模型

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

網格使用STAR－CCM＋軟件生成，網格類型為Trim非結構網格.在葉輪和導葉的導邊和隨邊處生成特征線，進行局部加密.在流道入口處流場變化劇烈，設置網格控制體進行局部加密.計算域總的網格數約為291萬個，葉輪旋轉區域約為94萬個，其余靜止區域約為197萬個.計算域的網格見圖5.

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

圖5 計算域的網格Fig.5 Grids of computational domain

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

計算域入口、兩側和底部設為速度進口；計算域出口和推進器噴口為壓力出口；計算域頂部、流道壁面和導葉設為無滑移壁面；將計算域分為旋轉區域和靜止區域，旋轉區域采用穩態多參考系（MFR）方法模擬泵的旋轉，在旋轉區域的速度根據旋轉角速度給出，葉輪和導流帽設為靜止壁面；旋轉區域和靜止區域的交界設為交界面.

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

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

3 計算結果及分析

對噴水推進器中的各水動力系數定義如下：

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

3.1 網格無關性分析

噴水推進器葉輪、導葉和流道以及加密區域的網格大小都與基本尺寸相關，由基本尺寸的變化生成由密到疏的5套網格，分別編號為1～5，來進行網格無關性的分析，5套網格的數目見表1.

在設計轉速n0下，進行網格無關性分析，將最密網格（1號網格）的推力和轉矩計算結果分別記為CTn1和Qshaft1，其他網格下推力和轉矩分別以該值為基準進行歸一化處理，網格無關性計算結果見圖7.

表1 網格的數目與尺寸Tab.1 Grid sizes and numbers of five grids

圖7 網格數目依賴性Fig.7 Dependence of grid number

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

3.2 噴水推進器性能分析

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

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

表2 影響噴水推進器性能的參數Tab.2 Parameters affecting the performance of waterjet

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

圖8 進流面形狀Fig.8 Shapes of capture area

表3 進流面的形狀參數Tab.3 Shape parameters of capture area

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

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

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

噴水推進器唇部的流場變化劇烈，會產生低壓區，容易發生空化現象，圖12所示為不同轉速下唇部表面壓力分布云圖和流線圖.定義無量綱的壓力系數為

式中：p為絕對壓力；p0為參考壓力（取1個標準大氣壓）.

由圖12可以發現，在唇部發生了流動分流，而發生流動分流的位置隨著轉速的增加向下移動，這是因為隨著轉速的增加，流量增加，需要進口從更深的位置抽吸水流.低壓區域與流動分流點的變化呈相反趨勢，低轉速時，低壓區域在唇部的下緣；高轉速時，低壓區域在唇部的上緣.

圖9 進流面的速度分布云圖Fig.9 Velocity contour of capture area

圖10 泵前進口的速度分布云圖Fig.10 Velocity contour of pump intake

圖11 泵前進口的速度分布云圖Fig.11 Velocity contour of pump intake

圖13所示為推力系數曲線，總推力系數曲線和凈推力系數曲線基本重合，說明ITTC推薦的動量流量法和壁面積分法都可以進行噴水推進器的推力預報，也間接說明本文對進流面獲取的準確性，動量流量計算的推力精度較高.在噴水推進計算中，如果采用壓力跳越法[5]或其他方法代替泵的作用，壁面積分法就不再適用，只能采用動量流量法來計算推力，因此動量流量法的適用性更廣，此時采用本文開發的進流面獲取程序將會大大提高效率.

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

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

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

圖14 噴水推進器性能曲線Fig.14 Performance curves of waterjet

4 結語

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

致謝 本研究得到了噴水推進技術重點實驗室的支持，得益于噴水推進技術學術研討會（WORKSHOP2018）多位與會專家的技術交流，在此深表謝意.