基于体积力的海水自流循环系统流速确定方法

陈康，姚志崇，周恩东，夏彬

1 深海载人装备国家重点实验室，江苏无锡214082

2 中国船舶科学研究中心，江苏无锡214082

0 引 言

在现代高速舰船和核动力潜艇中，可以采用自流循环系统为换热器提供冷却水［1］。艇体在巡航工况下时，冷却系统海水泵不工作而保持随动，利用艇体航行时与海水相对运动产生的动水头来抽吸海水，使海水能克服管路阻力自行流过换热器，从而有效降低海水泵冷却系统的负荷。艇体在高速工况下，需简化系统布置，无需另外设置海水泵；在低速工况下，采用海水泵可辅助实现海水循环功能，可降低海水循环系统对海水泵功率的需求，从而降低艇体的能耗，提高系统可靠性［2］。

自流循环系统自流能力的强弱决定了其能否满足换热需求，具有良好的自流循环能力也是设计自流循环系统的根本要求。因此，准确评估自流循环能力（自流循环系统海水流量或者流速）是设计自流循环系统的关键［3］。由于船身效应会产生压差，舰船自流循环系统进、出口缺乏有效的工程估算方法，因此无法确定舰船自流循环系统的自流能力。目前，舰船自流循环系统的自流能力一般采用CFD 方法进行计算。

裴金亮［4］在运用CFD 方法对自流循环系统及船体进行整体模拟时，由于自流循环系统中的换热器、海水泵、阀门等设备含有非常细小的部件，该部件与船体尺度相差较大，从而导致网格划分存在一定的困难。同时，模拟自流循环系统的流场需要布置大量的网格，计算量很大。可见，利用CFD 方法进行流量评估是自流循环系统设计过程中需要考虑的关键问题，由此采用用体积力法等效代替自流循环系统中部件的方法进行模拟具有优势。

采用力场模拟实体的体积力法在研究螺旋桨和船体的相互作用时运用较多，其还可以代替真实的螺旋桨参与仿真计算。Xu 等［5］和Paterson等［6］运用体积力法对船的伴流场进行了计算。傅慧萍等［7］采用螺旋桨体积力模型，实现了模型尺度下实船自航点的全粘带自由面计算。吴召华等［8］运用一种便捷的描述型体积力模型代替实体螺旋桨的作用，预报了船舶做单平面定常回转运动时的受力和力矩。力场模拟方法的计算量比船桨整体计算的小，还可以分析桨对船体的影响，为螺旋桨的设计和开发提供指导。参照体积力法代替螺旋桨的方法，可以将体积力法运用到自流循环系统的数值模拟中。

海水自流循环系统的管路阻力、关键部件（换热器、海水泵、阀门等）阻力都可以通过理论计算、仿真或实验得到，由于进、出口形状复杂，其与艇体的相互作用无法通过理论计算得到，故预报时只能通过仿真进行模拟［9］。为此，本文将仿照体积力等效代替螺旋桨的方法，采用体积力法等效代替换热器、海水泵和阀门等主要部件，对自流循环系统的流场进行数值模拟，建立一种可以确定船舶系统流速和流量的方法，从而为预报实艇海水自流循环系统的流速和流量提供参考。

1 体积力模拟原理

图1 所示为某艇的海水自流循环系统组成示意图。自流循环系统一般包括进水口部分、换热器部分、海水泵、中间管路和出水口部分［2］。艇体在水下航行时，利用航行产生的进、出口压差来驱动海水进入自流循环系统，通过换热器完成热量的交换。当艇体以低速巡航时，自流动力不足，需要开启海水泵抽吸海水。自流循环系统的水头平衡方程为

式中：vA，vB分别为进水口和出水口断面的平均速度；v0，v1，v2……为自流循环系统各断面处的平均速度；PA，PB分别为进水口和出水口断面压力；ρ为海水密度；g为重力加速度；λ为管路沿程阻力系数；L为管道长度；D为管道直径；ξ进为管路进口阻力系数；ξ出为管路出口阻力系数；ξ换为换热器阻力系数；ξ泵为海水泵阻力系数。

图1 自流循环系统组成示意图Fig.1 Composition of the artesian circulating system

式（1）左边是舰船以一定速度航行时，所能提供的总水头，右边是管路及各设备部件的水头损失，其中换热器和海水泵的水头损失为

若已知vA，vB，PA，PB以及各部件的损失系数，就可以通过式（1）直接求解出自流循环系统的流速。体积力模拟示意图如图2 所示。船体表面为复杂的三维曲面，自流循环水系统进、出口处的压力分布没有准确的理论求解方法，只能通过仿真计算方法对自流循环系统内流场和船体外流场整体进行求解来获取。自流循环系统内各设备部件可借助体积力来模拟其流动损失效应，从而避免直接对各设备部件进行数值模拟。

图2 体积力模拟示意图Fig.2 Body force simulation

为简化说明，系统中的部件只选取换热器。对计算域中换热器部分施加等效换热器水头损失的体积力。单位体积的体积力大小为fb，并使

式中：Fb为施加体积力的总值；dV为单位体积。体积力的流动控制方程为

在实际工程中，需要求解的是自流循环系统的流速，而等效体积力大小是无法用式（2）和式（3）确定的。使用体积力等效代替阻力部件，在对参考泵选型时，可以通过2 种特性曲线求交点的方法求得海水自流循环系统的流速点。水泵的性能曲线如图3 所示。在水泵的性能曲线上，任意流量点都可以找到与之相对应的扬程，该组流量与扬程参数即为工况。管路的特性曲线是管路的水头损失随管路流量变化的曲线，因此，管路的特性曲线与水泵性能曲线的交点即为运行工况点（即图3 中的M 点）。

根据上述方法，通过求取所替代部件的性能曲线与系统特性曲线交点的方法，可以确定实艇在不同航速下海水自流循环系统的工况点。系统特性曲线是通过体积力代替主要部件后再与艇体联合计算而求得，部件性能曲线是通过对所替代部件进行单独计算或仿真而得到。

图3 水泵的性能曲线交点图Fig.3 Intersection of the pump performance curves

由于用体积力法等效替代阻力部件的方法可以忽略不必要的流动细节，故在将海水自流循环系统与大尺寸的艇体联合计算时，采用体积力法可对阻力部件进行单独计算。同时，在保证计算精准度和有效性的前提下，体积力法还能大幅减少计算网格，提高计算效率。

2 自流循环系统流速确定过程

借助体积力等效代替换热器等部件来确定自流循环系统流速的方法，过程如下：

第1 步，在舰船特定航速下，求解自流循环系统的内流场及船体外流场，其中换热器所在区域全部设置为流体域，并分别施加一系列大小不同的体积力fb。计算求解后，提取换热器流体域两端截面的流量和压力，绘制出流量压差曲线。

第2 步，由式（2）设定流速值，求得换热器的水头损失，也即海水流经换热器后的压差。将流速换算为流量，同样绘制出流量压差曲线。

第3 步，采用类似于图3 的方法，求解前2 步绘制的流量压差曲线的交点，该点即为某航速下自流循环系统的自流点。第1 步求出的流量压差曲线相当于管路特性曲线，第2 步求出的流量压差曲线相当于水泵性能曲线。读取自流点对应的流量，就可求出相应的流速，也即自流循环系统的自流流速。

需要进一步说明的是，式（2）中换热器的流动损失系数ξ换需根据经验或实验确定。而根据舰船自流循环系统设计的需求，换热器一般无法使用标准换热器。在设计非标准换热器时，其流动损失系数未知，这时，可以采用CFD 方法确定。使用该方法确定流动损失系数，只需要单独模拟换热器，最重要的是，还能与自流循环系统及船体设计同步进行，且不影响系统多方案设计的对比，如自流循环系统进、出口流道的设计以及位置选择等。在海水自流循环系统上采用体积力法有2个优势：一是换热器采用CFD 模拟法可用于自流循环系统多方案设计中，而不是每设计一种换热器方案都要对整艇进行数值模拟；二是换热器和自流循环系统的设计可同步进行，从而缩短设计周期。

其他设备部件也可以采用体积力法进行等效代替。为简化说明，后续实例中的部件只考虑换热器，需要加入其他部件时只需要将部件的特性曲线并入到总的部件性能曲线中即可。

3 CFD 数值模拟方法

3.1 网格划分

整个自流循环系统的研究对象不仅包括系统内部的内流场，还包括艇体外部的外流场。为了模拟整个自流循环系统的内部流动特性，参加计算的部分包括外流场、进出水管道以及换热器。海水自流循环系统模拟的内、外流场的流动速度与几何尺寸相差较大，导致其流动特性也存在很大的不同。使用GAMBIT 软件分块画网格，然后再组合到一块进行内、外流场联合计算。实际选用圆柱形流体计算域，流场半径为艇体型深的12倍；在长度方向上，计算域向艇体艏部延伸艇长的2 倍，向艇体艉部延伸艇长的3 倍。自流循环系统流体域的网格划分如图4 所示。

图4 模型网格划分Fig.4 Model mesh division

在使用体积力法等效代替换热器进行数值模拟时，代入的是体积力的总值而不是系统的流场分布。因此，体积力的分布方式不影响自流循环系统流速点的求取，实际将体积力均匀分布到圆柱形流体域上。体积力法计算示意图如图5所示。图中：A，B 分别为海水进、出口截面；C，D 为施加体积力的圆柱形区域边界的截面。体积力代替部分网格数为1 万，整个流场的网格总数为800 万。

为验证体积力法等效代替部件方法的准确性，将体积力法与传统计算方法的结果进行对比。将换热器实体划分网格后加到自流循环系统流道中，然后再进行整艇的联合模拟计算。实体联合计算网格如图6 所示，最后采用的换热器网格数为1 036 万，整个流场网格总数为1 832 万。为了在模拟艇体外流场的同时兼顾自流循环系统的内流场，于计算后再划分网格，使换热器的壁面y+≈4，艇体艉部和自流循环系统进、出口的壁面y+≈102，艇体和外流场其他区域的壁面y+≈275。这里，用y+表征网格密度。

图5 体积力法计算示意图Fig.5 Calculation diagram of body force method

图6 实体联合计算示意图Fig.6 Combined calculation diagram of component and hull

自流循环系统使用的换热器以管壳式换热器为例。划分管壳式换热器网格的难点在于，换热器的换热管太过细长，与换热器整体尺度相比较小。由于系统部件的尺寸变化大，为保证计算的精准度，需要在换热管的内壁面划分大量网格。在划分换热器网格时，为了能更好地对管壁处的流场进行模拟，需要对壁面附近的网格进行加密处理。换热器网格如图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-ω模型在近壁面区域和远场区域计算的优点，同时增加了横向耗散导数项，在湍流粘度定义中考虑了湍流剪切应力的运输过程，更加适用于广域流场。因此，本文实际选用的是SST k-ω模型进行计算。

3.3 数值计算方法

数值模拟求解三维定常不可压缩RANS 方程，采用有限体积法离散动量方程，扩散项使用中心差分格式，对流项采用二阶迎风格式，压力速度耦合采用SIMPLE 算法，离散得到的代数方程组用Gauss-Seidel 法迭代求解，并使用代数多重网格法加速收敛。

3.4 体积力加载

将体积力均匀施加到以换热器直径、长度所定义的圆柱形流体域内，来实现实际部件对海水的阻力影响。按照表1 所示数值施加体积力，选定某特定航速后，通过施加对应的体积力，可以模拟得到该工况下系统的流速和流量。体积力施加方式如图2 所示，体积力施加方向与海水流速方向相反。

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

4 自流循环流速点的确定

4.1 系统特性曲线

图5 中，C，D 截面的压差随流速变化的曲线即为体积力代替部件后的系统特性曲线。选定航速后，施加表1 中航速对应的体积力进行模拟计算，3 种航速下不同体积力对应的海水管道流速如图8 所示。由图8 可以看出，在一定航速下，随着体积力的增加，自流循环的流速是下降的；在相同体积力下，随着航速的增加，自流循环的流速和流量是不断增加的。系统特性曲线如图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 部件性能曲线

自流循环系统的部件性能曲线为不同管道流速下流体域等效代替部件的压差。此处流体域等效代替的部件为管壳式换热器，换热器的管程压降可以通过工程理论［13］或单独仿真计算得到。在需要添加其他阻力部件时，可通过相同的方法计算或仿真以后，再并入到总的部件性能曲线当中。自流循环系统的阻力计算如表2 所示。

表2 部件压差计算Table 2 Differential pressure of the component

4.3 求取曲线交点

将系统特性曲线和部件性能曲线绘制到同一张图中，曲线的交点即为海水自流循环系统在特定航速下的流速点。图中两条曲线存在交点说明在该航速下可以实现海水的自流循环，若没有交点，则说明此种工况下的自流循环系统流速过低，不能满足冷却水量的需求，需要开启海水泵主动循环系统，以增加吸水扬程。部件性能曲线与系统特性曲线的交点图如图10 所示。

图10 系统特性曲线与部件性能曲线交点图Fig.10 Intersection of system characteristic curves and component performance curve

5 联合计算验证

为了验证体积力法的准确性，将联合计算法与体积力法计算得到的自流循环工况点进行了比较。单独划分换热器的网格后，将换热器实体加入到自流循环管路中与艇体一起进行数值模拟，仿真计算的参数设置与体积力法的参数设置相同。选取2，4，6 kn 航速下的工况，计算得到自流循环系统的流速和部件两端的压差，体积力法与联合计算法的计算结果如表3 所示。

将体积力法与实体联合计算法得到的自流循环工况点进行比较，结果如图11 所示。由图可看出，2 种方法下的流速值与压差值相差均在5%以内。对比2 种方法的网格总数，发现体积力法具有显著优势。

表3 自流循环工况点比较Table 3 Comparison of the operating point

图11 体积力法与联合计算法工况点比较Fig.11 Comparison of operating point between body force method and combined calculation

6 结 语

本文基于体积力法建立了确定自流循环系统流速的分析方法，运用体积力法对自流循环系统进行数值模拟，类似于确定水泵工况点的方法，通过求取系统特性曲线与部件性能曲线交点的方法得到海水自循环系统的流速点。通过比较实体联合计算法和体积力法这2 种方法下的自流循环工况点，表明体积力模拟方法可以确定海水自循环系统的流速，而且计算量小。后续若需添加其他阻力部件，可以直接将对应的部件性能曲线并入到总的部件性能曲线当中，通过求交点的方法来得出实际海水流量。有关仿真结果的准确性还有待后续的试验验证。