基于體積力的海水自流循環(huán)系統(tǒng)流速確定方法

陳康，姚志崇，周恩東，夏彬

1 深海載人裝備國家重點實驗室，江蘇無錫214082

2 中國船舶科學研究中心，江蘇無錫214082

0 引 言

在現(xiàn)代高速艦船和核動力潛艇中，可以采用自流循環(huán)系統(tǒng)為換熱器提供冷卻水［1］。艇體在巡航工況下時，冷卻系統(tǒng)海水泵不工作而保持隨動，利用艇體航行時與海水相對運動產生的動水頭來抽吸海水，使海水能克服管路阻力自行流過換熱器，從而有效降低海水泵冷卻系統(tǒng)的負荷。艇體在高速工況下，需簡化系統(tǒng)布置，無需另外設置海水泵；在低速工況下，采用海水泵可輔助實現(xiàn)海水循環(huán)功能，可降低海水循環(huán)系統(tǒng)對海水泵功率的需求，從而降低艇體的能耗，提高系統(tǒng)可靠性［2］。

自流循環(huán)系統(tǒng)自流能力的強弱決定了其能否滿足換熱需求，具有良好的自流循環(huán)能力也是設計自流循環(huán)系統(tǒng)的根本要求。因此，準確評估自流循環(huán)能力（自流循環(huán)系統(tǒng)海水流量或者流速）是設計自流循環(huán)系統(tǒng)的關鍵［3］。由于船身效應會產生壓差，艦船自流循環(huán)系統(tǒng)進、出口缺乏有效的工程估算方法，因此無法確定艦船自流循環(huán)系統(tǒng)的自流能力。目前，艦船自流循環(huán)系統(tǒng)的自流能力一般采用CFD 方法進行計算。

裴金亮［4］在運用CFD 方法對自流循環(huán)系統(tǒng)及船體進行整體模擬時，由于自流循環(huán)系統(tǒng)中的換熱器、海水泵、閥門等設備含有非常細小的部件，該部件與船體尺度相差較大，從而導致網格劃分存在一定的困難。同時，模擬自流循環(huán)系統(tǒng)的流場需要布置大量的網格，計算量很大?？梢姡肅FD 方法進行流量評估是自流循環(huán)系統(tǒng)設計過程中需要考慮的關鍵問題，由此采用用體積力法等效代替自流循環(huán)系統(tǒng)中部件的方法進行模擬具有優(yōu)勢。

采用力場模擬實體的體積力法在研究螺旋槳和船體的相互作用時運用較多，其還可以代替真實的螺旋槳參與仿真計算。Xu 等［5］和Paterson等［6］運用體積力法對船的伴流場進行了計算。傅慧萍等［7］采用螺旋槳體積力模型，實現(xiàn)了模型尺度下實船自航點的全粘帶自由面計算。吳召華等［8］運用一種便捷的描述型體積力模型代替實體螺旋槳的作用，預報了船舶做單平面定?；剞D運動時的受力和力矩。力場模擬方法的計算量比船槳整體計算的小，還可以分析槳對船體的影響，為螺旋槳的設計和開發(fā)提供指導。參照體積力法代替螺旋槳的方法，可以將體積力法運用到自流循環(huán)系統(tǒng)的數值模擬中。

海水自流循環(huán)系統(tǒng)的管路阻力、關鍵部件（換熱器、海水泵、閥門等）阻力都可以通過理論計算、仿真或實驗得到，由于進、出口形狀復雜，其與艇體的相互作用無法通過理論計算得到，故預報時只能通過仿真進行模擬［9］。為此，本文將仿照體積力等效代替螺旋槳的方法，采用體積力法等效代替換熱器、海水泵和閥門等主要部件，對自流循環(huán)系統(tǒng)的流場進行數值模擬，建立一種可以確定船舶系統(tǒng)流速和流量的方法，從而為預報實艇海水自流循環(huán)系統(tǒng)的流速和流量提供參考。

1 體積力模擬原理

圖1 所示為某艇的海水自流循環(huán)系統(tǒng)組成示意圖。自流循環(huán)系統(tǒng)一般包括進水口部分、換熱器部分、海水泵、中間管路和出水口部分［2］。艇體在水下航行時，利用航行產生的進、出口壓差來驅動海水進入自流循環(huán)系統(tǒng)，通過換熱器完成熱量的交換。當艇體以低速巡航時，自流動力不足，需要開啟海水泵抽吸海水。自流循環(huán)系統(tǒng)的水頭平衡方程為

式中：vA，vB分別為進水口和出水口斷面的平均速度；v0，v1，v2……為自流循環(huán)系統(tǒng)各斷面處的平均速度；PA，PB分別為進水口和出水口斷面壓力；ρ為海水密度；g為重力加速度；λ為管路沿程阻力系數；L為管道長度；D為管道直徑；ξ進為管路進口阻力系數；ξ出為管路出口阻力系數；ξ換為換熱器阻力系數；ξ泵為海水泵阻力系數。

圖1 自流循環(huán)系統(tǒng)組成示意圖Fig.1 Composition of the artesian circulating system

式（1）左邊是艦船以一定速度航行時，所能提供的總水頭，右邊是管路及各設備部件的水頭損失，其中換熱器和海水泵的水頭損失為

若已知vA，vB，PA，PB以及各部件的損失系數，就可以通過式（1）直接求解出自流循環(huán)系統(tǒng)的流速。體積力模擬示意圖如圖2 所示。船體表面為復雜的三維曲面，自流循環(huán)水系統(tǒng)進、出口處的壓力分布沒有準確的理論求解方法，只能通過仿真計算方法對自流循環(huán)系統(tǒng)內流場和船體外流場整體進行求解來獲取。自流循環(huán)系統(tǒng)內各設備部件可借助體積力來模擬其流動損失效應，從而避免直接對各設備部件進行數值模擬。

圖2 體積力模擬示意圖Fig.2 Body force simulation

為簡化說明，系統(tǒng)中的部件只選取換熱器。對計算域中換熱器部分施加等效換熱器水頭損失的體積力。單位體積的體積力大小為fb，并使

式中：Fb為施加體積力的總值；dV為單位體積。體積力的流動控制方程為

在實際工程中，需要求解的是自流循環(huán)系統(tǒng)的流速，而等效體積力大小是無法用式（2）和式（3）確定的。使用體積力等效代替阻力部件，在對參考泵選型時，可以通過2 種特性曲線求交點的方法求得海水自流循環(huán)系統(tǒng)的流速點。水泵的性能曲線如圖3 所示。在水泵的性能曲線上，任意流量點都可以找到與之相對應的揚程，該組流量與揚程參數即為工況。管路的特性曲線是管路的水頭損失隨管路流量變化的曲線，因此，管路的特性曲線與水泵性能曲線的交點即為運行工況點（即圖3 中的M 點）。

根據上述方法，通過求取所替代部件的性能曲線與系統(tǒng)特性曲線交點的方法，可以確定實艇在不同航速下海水自流循環(huán)系統(tǒng)的工況點。系統(tǒng)特性曲線是通過體積力代替主要部件后再與艇體聯(lián)合計算而求得，部件性能曲線是通過對所替代部件進行單獨計算或仿真而得到。

圖3 水泵的性能曲線交點圖Fig.3 Intersection of the pump performance curves

由于用體積力法等效替代阻力部件的方法可以忽略不必要的流動細節(jié)，故在將海水自流循環(huán)系統(tǒng)與大尺寸的艇體聯(lián)合計算時，采用體積力法可對阻力部件進行單獨計算。同時，在保證計算精準度和有效性的前提下，體積力法還能大幅減少計算網格，提高計算效率。

2 自流循環(huán)系統(tǒng)流速確定過程

借助體積力等效代替換熱器等部件來確定自流循環(huán)系統(tǒng)流速的方法，過程如下：

第1 步，在艦船特定航速下，求解自流循環(huán)系統(tǒng)的內流場及船體外流場，其中換熱器所在區(qū)域全部設置為流體域，并分別施加一系列大小不同的體積力fb。計算求解后，提取換熱器流體域兩端截面的流量和壓力，繪制出流量壓差曲線。

第2 步，由式（2）設定流速值，求得換熱器的水頭損失，也即海水流經換熱器后的壓差。將流速換算為流量，同樣繪制出流量壓差曲線。

第3 步，采用類似于圖3 的方法，求解前2 步繪制的流量壓差曲線的交點，該點即為某航速下自流循環(huán)系統(tǒng)的自流點。第1 步求出的流量壓差曲線相當于管路特性曲線，第2 步求出的流量壓差曲線相當于水泵性能曲線。讀取自流點對應的流量，就可求出相應的流速，也即自流循環(huán)系統(tǒng)的自流流速。

需要進一步說明的是，式（2）中換熱器的流動損失系數ξ換需根據經驗或實驗確定。而根據艦船自流循環(huán)系統(tǒng)設計的需求，換熱器一般無法使用標準換熱器。在設計非標準換熱器時，其流動損失系數未知，這時，可以采用CFD 方法確定。使用該方法確定流動損失系數，只需要單獨模擬換熱器，最重要的是，還能與自流循環(huán)系統(tǒng)及船體設計同步進行，且不影響系統(tǒng)多方案設計的對比，如自流循環(huán)系統(tǒng)進、出口流道的設計以及位置選擇等。在海水自流循環(huán)系統(tǒng)上采用體積力法有2個優(yōu)勢：一是換熱器采用CFD 模擬法可用于自流循環(huán)系統(tǒng)多方案設計中，而不是每設計一種換熱器方案都要對整艇進行數值模擬；二是換熱器和自流循環(huán)系統(tǒng)的設計可同步進行，從而縮短設計周期。

其他設備部件也可以采用體積力法進行等效代替。為簡化說明，后續(xù)實例中的部件只考慮換熱器，需要加入其他部件時只需要將部件的特性曲線并入到總的部件性能曲線中即可。

3 CFD 數值模擬方法

3.1 網格劃分

整個自流循環(huán)系統(tǒng)的研究對象不僅包括系統(tǒng)內部的內流場，還包括艇體外部的外流場。為了模擬整個自流循環(huán)系統(tǒng)的內部流動特性，參加計算的部分包括外流場、進出水管道以及換熱器。海水自流循環(huán)系統(tǒng)模擬的內、外流場的流動速度與幾何尺寸相差較大，導致其流動特性也存在很大的不同。使用GAMBIT 軟件分塊畫網格，然后再組合到一塊進行內、外流場聯(lián)合計算。實際選用圓柱形流體計算域，流場半徑為艇體型深的12倍；在長度方向上，計算域向艇體艏部延伸艇長的2 倍，向艇體艉部延伸艇長的3 倍。自流循環(huán)系統(tǒng)流體域的網格劃分如圖4 所示。

圖4 模型網格劃分Fig.4 Model mesh division

在使用體積力法等效代替換熱器進行數值模擬時，代入的是體積力的總值而不是系統(tǒng)的流場分布。因此，體積力的分布方式不影響自流循環(huán)系統(tǒng)流速點的求取，實際將體積力均勻分布到圓柱形流體域上。體積力法計算示意圖如圖5所示。圖中：A，B 分別為海水進、出口截面；C，D 為施加體積力的圓柱形區(qū)域邊界的截面。體積力代替部分網格數為1 萬，整個流場的網格總數為800 萬。

為驗證體積力法等效代替部件方法的準確性，將體積力法與傳統(tǒng)計算方法的結果進行對比。將換熱器實體劃分網格后加到自流循環(huán)系統(tǒng)流道中，然后再進行整艇的聯(lián)合模擬計算。實體聯(lián)合計算網格如圖6 所示，最后采用的換熱器網格數為1 036 萬，整個流場網格總數為1 832 萬。為了在模擬艇體外流場的同時兼顧自流循環(huán)系統(tǒng)的內流場，于計算后再劃分網格，使換熱器的壁面y+≈4，艇體艉部和自流循環(huán)系統(tǒng)進、出口的壁面y+≈102，艇體和外流場其他區(qū)域的壁面y+≈275。這里，用y+表征網格密度。

圖5 體積力法計算示意圖Fig.5 Calculation diagram of body force method

圖6 實體聯(lián)合計算示意圖Fig.6 Combined calculation diagram of component and hull

自流循環(huán)系統(tǒng)使用的換熱器以管殼式換熱器為例。劃分管殼式換熱器網格的難點在于，換熱器的換熱管太過細長，與換熱器整體尺度相比較小。由于系統(tǒng)部件的尺寸變化大，為保證計算的精準度，需要在換熱管的內壁面劃分大量網格。在劃分換熱器網格時，為了能更好地對管壁處的流場進行模擬，需要對壁面附近的網格進行加密處理。換熱器網格如圖7 所示。

圖7 換熱器網格圖Fig.7 Heat exchanger grid

3.2 湍流模型

仿真計算涉及艇體的外部流場和內部海水管路流場。在外流場和內流場的混合計算過程中，需采用合適的湍流模型。FLUENT 軟件提供的兩方程湍流模型主要包括標準k-ε模型、RNG k-ε模型、Realizable k-ε模型、標準k-ω模型和SST k-ω模型［10］。張楠等［11］結合這幾種湍流模型，將計算得到的潛艇阻力和尾流場與模型試驗結果進行了對比分析，驗證了數值方法的可靠性。張志榮等［12］針對這幾種不同湍流模型的粘性繞流場進行了數值模擬，將計算得出的阻力系數、槳盤面處流場與試驗值進行了比較。結果表明，在復雜船型粘性流場計算的湍流模式中，和標準k-ω模型相比，使用SST k-ω模型具有更高的精準度和可靠性，其綜合了k-ω模型在近壁面區(qū)域和遠場區(qū)域計算的優(yōu)點，同時增加了橫向耗散導數項，在湍流粘度定義中考慮了湍流剪切應力的運輸過程，更加適用于廣域流場。因此，本文實際選用的是SST k-ω模型進行計算。

3.3 數值計算方法

數值模擬求解三維定常不可壓縮RANS 方程，采用有限體積法離散動量方程，擴散項使用中心差分格式，對流項采用二階迎風格式，壓力速度耦合采用SIMPLE 算法，離散得到的代數方程組用Gauss-Seidel 法迭代求解，并使用代數多重網格法加速收斂。

3.4 體積力加載

將體積力均勻施加到以換熱器直徑、長度所定義的圓柱形流體域內，來實現(xiàn)實際部件對海水的阻力影響。按照表1 所示數值施加體積力，選定某特定航速后，通過施加對應的體積力，可以模擬得到該工況下系統(tǒng)的流速和流量。體積力施加方式如圖2 所示，體積力施加方向與海水流速方向相反。

表1 體積力施加值Table 1 Values of the applied body force at 2，4，6 knots

4 自流循環(huán)流速點的確定

4.1 系統(tǒng)特性曲線

圖5 中，C，D 截面的壓差隨流速變化的曲線即為體積力代替部件后的系統(tǒng)特性曲線。選定航速后，施加表1 中航速對應的體積力進行模擬計算，3 種航速下不同體積力對應的海水管道流速如圖8 所示。由圖8 可以看出，在一定航速下，隨著體積力的增加，自流循環(huán)的流速是下降的；在相同體積力下，隨著航速的增加，自流循環(huán)的流速和流量是不斷增加的。系統(tǒng)特性曲線如圖9 所示，從圖9 中可以看出，在相同海水管道流速和流量下，隨著航速的增加，C，D 截面的壓差逐漸增大，符合實際趨勢。

圖8 在2，4，6 kn 航速下不同體積力對應海水管道流速Fig.8 Velocities of sea water pipe at 2，4，6 knots in different body forces

圖9 在2，4，6 kn 航速下不同流速對應C，D 截面的壓差Fig.9 Differential pressure between section C and D at 2，4，6 knots

4.2 部件性能曲線

自流循環(huán)系統(tǒng)的部件性能曲線為不同管道流速下流體域等效代替部件的壓差。此處流體域等效代替的部件為管殼式換熱器，換熱器的管程壓降可以通過工程理論［13］或單獨仿真計算得到。在需要添加其他阻力部件時，可通過相同的方法計算或仿真以后，再并入到總的部件性能曲線當中。自流循環(huán)系統(tǒng)的阻力計算如表2 所示。

表2 部件壓差計算Table 2 Differential pressure of the component

4.3 求取曲線交點

將系統(tǒng)特性曲線和部件性能曲線繪制到同一張圖中，曲線的交點即為海水自流循環(huán)系統(tǒng)在特定航速下的流速點。圖中兩條曲線存在交點說明在該航速下可以實現(xiàn)海水的自流循環(huán)，若沒有交點，則說明此種工況下的自流循環(huán)系統(tǒng)流速過低，不能滿足冷卻水量的需求，需要開啟海水泵主動循環(huán)系統(tǒng)，以增加吸水揚程。部件性能曲線與系統(tǒng)特性曲線的交點圖如圖10 所示。

圖10 系統(tǒng)特性曲線與部件性能曲線交點圖Fig.10 Intersection of system characteristic curves and component performance curve

5 聯(lián)合計算驗證

為了驗證體積力法的準確性，將聯(lián)合計算法與體積力法計算得到的自流循環(huán)工況點進行了比較。單獨劃分換熱器的網格后，將換熱器實體加入到自流循環(huán)管路中與艇體一起進行數值模擬，仿真計算的參數設置與體積力法的參數設置相同。選取2，4，6 kn 航速下的工況，計算得到自流循環(huán)系統(tǒng)的流速和部件兩端的壓差，體積力法與聯(lián)合計算法的計算結果如表3 所示。

將體積力法與實體聯(lián)合計算法得到的自流循環(huán)工況點進行比較，結果如圖11 所示。由圖可看出，2 種方法下的流速值與壓差值相差均在5%以內。對比2 種方法的網格總數，發(fā)現(xiàn)體積力法具有顯著優(yōu)勢。

表3 自流循環(huán)工況點比較Table 3 Comparison of the operating point

圖11 體積力法與聯(lián)合計算法工況點比較Fig.11 Comparison of operating point between body force method and combined calculation

6 結 語

本文基于體積力法建立了確定自流循環(huán)系統(tǒng)流速的分析方法，運用體積力法對自流循環(huán)系統(tǒng)進行數值模擬，類似于確定水泵工況點的方法，通過求取系統(tǒng)特性曲線與部件性能曲線交點的方法得到海水自循環(huán)系統(tǒng)的流速點。通過比較實體聯(lián)合計算法和體積力法這2 種方法下的自流循環(huán)工況點，表明體積力模擬方法可以確定海水自循環(huán)系統(tǒng)的流速，而且計算量小。后續(xù)若需添加其他阻力部件，可以直接將對應的部件性能曲線并入到總的部件性能曲線當中，通過求交點的方法來得出實際海水流量。有關仿真結果的準確性還有待后續(xù)的試驗驗證。