0 引　言

北极地区能源储量丰富，随着全球气候变化和周边国家经济技术发展，各国对于北极自然环境、航道开发和资源分布的研究考察愈发深入。极地船舶是极地考察、航运和资源开采的基础装备。北极地区地理位置特殊，常年被冰雪覆盖，气候环境恶劣，冬季平均气温仅为−34 ℃，极端情况下最低气温可低至−67 ℃^{[1 − 3]}。极地长期的低温和冰雪环境是海洋装备结构设计的重要影响因素，因为船舶在极地作业过程中常会面临结冰问题，甲板、上建、内部管路、轮机和电气系统等能否在低温结冰情况下正常运作，是极地船舶设计制造和航运实践中需要重点关注的问题^[4]。

目前,极地船舶结冰研究方向多集中在水滴结冰或船舶甲板等船舶外部结构结冰现象研究，通常采用经验分析方法、模型试验和数值模拟等研究方法。韩端锋等^[5]运用Wheeler相场模型模拟海水六角冰晶的生长，揭示了冰晶形成机理，并对无量纲过冷度、各向异性强度和噪声强度等参数进行敏感性分析。谢强等^[6]阐述了极地海洋平台甲板结冰的各项危害并归纳总结了现有除冰技术的优缺点以及适用情况；沈杰等^{[7 − 8]}将 Fluent 软件和 FENSAPICE 软件相结合，对极地船舶杆件结构和极地海洋平台甲板进行数值模拟，得到不同风速影响下的表面结冰量变化。针对船舶低温条件下的结构响应研究，以LNG船低温货舱的温度场分布和结构热应力研究为主，比如，余祥虎^[9]以一艘薄膜型LNG船为研究对象，利用Ansys计算了1/4液舱段液舱绝热层殷钢膜主屏壁绝热膜完好和殷钢膜主屏壁绝热膜破损2种工况下，液舱温度及温度应力分布情况。李敏^[10]探讨了液货舱的传热方式，利用Ansys有限元软件对1/2液货舱进行建模，施加边界条件和载荷，对液货舱的稳态温度场进行数值模拟，得到了液货舱的温度场分布图并对结果进行了分析。目前，虽然针对船体液舱内海水结冰现象的研究工作已经陆续开展，但考虑船舶真实极地作业环境、多种温度源共同影响下的舱室结冰过程研究工作尚未广泛开展，针对液舱结冰引起的船体结构安全问题也尚未引起足够重视，可查询的相关文献资料很少。各国船级社虽然都对液舱在极地温度的技术参数以及防寒措施方面制定了详细规范，但都尚未进行明确说明。

综上所述，尽管对极地船舶在低温环境下液舱内海水结冰现象的研究工作已经初步开展，但目前研究成果较少，对于低温条件下船舶液舱结冰相变行为的评估方法明显不足。因此，为了准确评估低温环境下极地船舶液舱结冰相变行为，本文以极地船舶单侧舷侧压载舱为研究对象，基于水结冰相变计算方法，运用多物理场仿真分析软件COMSOL Multiphysics，建立含压载水的液舱传热模型，分析低温环境下压载舱冻结的水-冰相变过程，计算温度场的时空分布，研究舱内压载水结冰相变过程，进而开展舱内液位高度、舱外空气温度和室温等参数对压载舱内水结冰相变行为的敏感性分析。

1 压载舱内水冰相变行为数值模型 1.1 传热方式

本文以某极地船舶靠近船中区域的舷侧压载舱为研究对象，纵向取2个横舱壁之间，长18 m，宽度方向取内、外壳体之间的结构，高度取型深12 m。因位于舷侧，压载水舱内空气和水域温度主要受到外界冷空气和海水的影响。根据船体结构特征，压载舱被纵向平台板分为4个舱室，如图1所示。H₄为最上层舱室空气域，下方H₁～H₃为压载水舱，被H₂顶端舱内平台板和H₃顶端的内底板隔离开；图中H₁内压载水不满舱，压载水面以上为空气域。降温过程中，压载舱内存在的传热形式以图中箭头表示，加粗箭头表示热传导，曲线箭头表示热对流，传热过程涉及到热量在物质内部，如空气域和水域内部的传递，以及物体与周围环境之间，如船体板与外界空气之间的热量交换。

船舶在极地航行过程中，压载舱内主要传热方式包括：1）热量在压载水以及空气域中的热传导；2）水线以下船体与外部低温海水进行自然对流换热；3）水线以上船体与外界冷空气进行自然对流换热；4）压载舱内空气与压载水进行自然对流换热；5）压载舱内冷热空气之间自然对流换热。

船舶在极地航行过程中，热量在压载舱空气与压载水中的传递可以看作三维的流体传热问题，过程遵循热传导定律（傅里叶导热定律），传热方程为：

$ \rho {C_{{p}}}\frac{{\partial {{T}}}}{{\partial {{t}}}} + \rho {C_{{p}}}{{u}} \cdot \nabla {{T}} + \nabla \cdot {{q}} = {{Q}} + {{{Q}}_{{p}}} + {{{Q}}_{{{vd}}}}，$

(1)

$ {{q}} = - {{k}}\nabla {{T}} ，$

(2)

$ \rho =\frac{{P}_{A}}{{R}_{\text{s}}T}，\mathrm{在}\mathrm{理}\mathrm{想}\mathrm{气}\mathrm{体}\mathrm{域}\mathrm{中}。$

(3)

式中：T为任意时刻流体温度，K；t为时间，s；C_p为流体恒压热容，J/(kg·K)；u为流体速度，m/s；ρ为流体密度，kg/m³；$ \nabla T $ 为温度梯度，K；$ \nabla $为哈密顿微分算子；q为通量，W/m²；k为流体导热系数，W/(m·K)；Q为热源，W/m³；Q_P为压力功，J；Q_vd为黏性耗散，J；P_A为绝对压力，Pa；R_S为比气体常数，J/(kg·K)。

船体钢板、空气以及水中的热传导过程中往往伴随着热对流现象的发生，热对流的常见情况是流体与温度不同的固体壁面接触时发生热量交换，其热流密度q与固体表面温度T和流体温度T₀之差成比例，比例系数为对流换热系数h的表达式为：

$ q = h\left( {T - {T_0}} \right)。$

(4)

式中：h为对流换热系数，W/(m²·K)。

1.2 水冰相变函数

相变是指物质由一种相态转变为另一种相态的过程，水由液态凝固为固态冰的相变过程中，尽管温度保持不变，系统却会释放热量，即为潜热（333 J/g）。计算时使用融化/凝固模型来描述这一过程，在COMSOL热分析过程，相变潜热会在相变的界面发生。通过引入一个相变温度范围，以平滑过渡的方式处理相变的过程。当温度低于冻结点时，水相向冰相转变，温度场参数和不同相的组分变化为：

$ \rho = {\theta _1}{\rho _1} + {\theta _2}{\rho _2}，$

(5)

$ {c_p} = \frac{1}{\rho }({\theta _1}{\rho _1}{c_{p,1}} + {\theta _2}{\rho _2}{C_{p,2}}) + {L_{1 \to 2}}\frac{{\partial {\alpha _m}}}{{\partial T}}，$

(6)

$ {\alpha _m} = \frac{1}{2}\frac{{{\theta _2}{\rho _2} - {\theta _1}{\rho _1}}}{{{\theta _1}{\rho _1} + {\theta _2}{\rho _2}}}，$

(7)

$ k={\theta }_{1}{k}_{1}+{\theta }_{2}{k}_{2}\text{，}{\theta }_{1}+{\theta }_{2}=1。$

(8)

式中：θ₁、θ₂为相1、相2百分数；L_1→2为相1到相2的潜热，J/kg；ρ₁、ρ₂为相1、相2密度，kg/m³；C_p,1、C_p,2为相1、相2恒压热容，J/(kg·K)；k₁、k₂为相1、相2导热系数，W/(m·K)。

相变函数如图2所示，当实施相变函数α(T)时，在相变温度T_pc，12 的区间ΔT_1→2内，发生了固相到液相的平滑过渡。设θ₀为水相的百分数，当θ₀=1表明系统均为水相，还未发生冻结凝固；当θ₀＜1表明部分水相转变为冰相；当θ₀=0表明所有水相完成冻结相变过程，全都转化为冰相。因此，当θ₀和θ₁的值均处于0～1间，表明此时系统内为固液混合物。

1.3 压载舱传热数值模型及边界条件

压载舱传热数值模型如图3所示。图3(a)为舱室几何模型，图3(b)为压载舱数值模型网格（局部），舱室结构件包括甲板、船底板、舷侧外板、内侧纵舱壁板、纵桁、肋骨及其他扶强材网格采用二维壳单元建模，图3(c)为液舱中水以及空气网格，采用三维实体单元建模。网格特征尺寸300 mm，实体单元共33306个，壳单元共33787个。

在压载舱的舷侧外板与外界海水接触处、舷侧外板与外界冷空气接触处以及内侧纵壁板与船体内部空气接触处，采用定温边界条件。设置外界低温冷空气温度为−40 ℃，外界海水温度为−2 ℃，船体内部空气温度为5 ℃。压载水液面、舱内空气以及船体板之间的热交换通过对流边界条件描述，设置相关界面处的换热系数。

1.4 材料参数

压载舱水冰相变数值仿真模型中采用材料包含低温钢EH36、舱内外海水、空气以及冰体，热分析所需材料参数见表1。

2 压载水结冰相变过程数值模拟 2.1 压载水满舱结冰相变分析

研究对象吃水为6.5 m，压载舱中满舱。舱内压载水和空气的初始温度设为0.15 ℃，外部空气温度为−40 ℃，外部海水温度为−2 ℃，船体内部温度为5 ℃，设置舱内压载水的相变温度为0 ℃。计算得到了压载舱降温冻结过程中，舱内温度场随时间的动态变化，分别给出冻结时间3、24、100、300、500、960 h共6个典型时刻的整舱温度场情况，以及压载水和结冰后冰水混合物的冻结状态。

压载舱内温度场分布和冻结状态随时间变化情况如图4所示。随冻结时间增加，压载舱室舷侧上部首先开始降温，最上层舱内空气域降温幅度大于舱内压载水，原因是水的比热容远大于空气。同时，舱内压载水高出吃水的部分降温幅度大于水线以下部分，原因是舷侧钢板的比热很小，导热性强，致使水线以上接近舷侧外板的压载水受到外界冷空气控制急速降温结冰，水线以下的压载水的缓慢降温则是受到外界海水影响。

冻结的起始点为压载水自由液面靠近舷侧外板处，此处距离外界冷空气的距离最近，结冰由该处向下和向内扩张。图中冻结状态数轴中冻结率为1表示完全成冰，0表示水，0～1表示冰水混合物。冻结时间3 h，处于水线以上部分的压载水首先降温结冰，冰块先在水线以上的舷侧外板内壁面形成。冻结时间24～100 h的过程中，在外界冷空气与低温海水的共同作用下，冰块沿着舷侧外板生长，并沿着液面横向向内、向下生长，向下的生长速率大于向内生长速率，100 h时压载水液面上尚未形成完整的横向贯通冰盖。冻结时间达到300 h时刻，压载水液面基本冻结形成完整的冰盖，并沿着内侧舱壁向下持续发展，在500～960 h过程中以后，舱内温度场变化已不明显，压载水冻结区域逐渐稳定，最终形成了较大的冰体。

为了进一步观测压载舱内冻结过程中的水/冰温度变化情况，取压载水自由液面靠近舷侧外板点a、内壁点c、最上层舱室内壁与自由液面顶部交点b、外板侧距船底高5 m处d和3 m处e，作为温度监测点，监测点温度随冻结时间变化情况如图5所示。

a点由于通过舷侧外板板传热直接受到外界冷空气的控制作用，温度骤降伴随着急剧冻结相变，冻结时间越长a点温度越低，200 h后降温趋于平缓，最终温度接近−30 ℃。初始时刻，b点和c点主要受到船体内部热空气的控制作用，温度−时间曲线先呈现处短暂小幅上升，再随着压载水冻结趋势由舷侧向内壁发展，b、c点温度下降，由于冻结的横向发展较缓慢，b、c点降温较为平缓。

压载舱深度方向上，由a点至d点再到e点，离外界冷空气的距离逐渐变远，降温速率越低，温度-时间曲线越平缓，发生水冰相变的时间越晚，最终温度越高，表明外界冷空气是液舱内降温速率和舱内压载水冻结相变的最重要影响因素。

2.2 数值模型可靠性验证

数值模拟中，网格的数量和质量都有可能会影响到计算结果的可靠性，因此网格划分时，应在保证质量的前提下确定合适的网格尺寸，使计算结果在满足精度要求的同时又可以有较高的计算效率。为探究网格尺寸和密度对舱室降温和结冰过程计算结果的影响，针对液舱传热几何模型划分了3种不同单元尺寸的网格模型，特征尺寸分别为1号网格450 mm、2号网格300 mm和3号网格150 mm，模型单元数量见表2。

不同特征尺寸的模型降温冻结计算结果中，取a点温度变化情况进行对比，见图6。a点的温度-时间变化曲线基本重合，表明3种尺寸网格的模型，同一测点降温速率和同一时刻的温度数值几乎一致，网格尺寸对计算结果的影响很小。本文舱内温度场和相变过程数值模拟，采用当前网格特征尺寸300 mm的数值模型计算结果具有可靠性，综合考虑计算效率和结果准确性，选择网格特征尺寸300 mm继续开展结冰行为影响因素计算分析。

3 结冰行为影响因素分析 3.1 舱内水位影响分析

船舶航行、作业过程中，为调整姿态或吃水，需合理调控压载水量，因而导致压载舱内液位高度改变。压载水液位高度改变，液舱中空气和压水占比随之增减。为探究舱内液位对冻结情况的影响，设置舱内液位高度3、5、6.5 m（满舱），热边界条件与前文计算工况一致。

不同液位高度下，冻结时间960 h的冻结状态计算结果如图7所示。压载水液位越低，相同冻结时间下的结冰程度越低。当舱内水位高度为3 m时，仅液面接近舷侧外板上部发生冻结，未形成完全覆盖压载水液面的冰盖。液位高度升至为5～6.5 m时，结冰区域体积显著增大，冻结率增加。原因是液位越低，压载水液面越远离外界冷空气，对流传热先导致舱内空气域温度降低，其次在降温后的舱内空气和外界海水共同影响下压载水缓慢降温结冰。与此同时，舱内水位低导致内壁侧热空气对流传热面积较大，室内热空气对舱内降温结冰的延缓作用增加。

提取自由液面靠近舷侧外板侧点a和靠近内侧纵舱壁板侧点c的温度−时间变化曲线。压载水液位高度越高，自由液面温度监测点降温速率越快。压载水液位高度3 m，a点温度曲线在前10 h为平直段，表明初始时段a点附近水域尚未冻结，结冰后温度下降至−3℃左右，但由于远离外界冷空气，最终温度接近外界海水。压载液位高度增至5 m和满舱6.5 m时，a点降温速率越快，最终温度越低。

c点受到船体内部热空气控制作用明显，压载水液位高度3 m时，c点附近水域始终没有冻结凝固，温度接近室内空气温度。液位高度增至5 m和满舱后，c点温度经历初始上升段后出现明显下降，最终温度远高于左侧测点，分别为0 ℃和−2 ℃以下。

3.2 舱内外空气温度影响

外界冷空气温度是影响结冰过程最重要的因素之一。为了探究外界冷空气和舱室内热空气对结冰过程的影响效果，在保持舱内液位高度5 m、外界海水温度−2 ℃、液舱内初始温度0.15 ℃、冻结时长为960 h的条件不变下，设定了15组不同的舱室内外环境温度，如表3所示。每3个工况为一组，室内温度相同，舱外空气温度为−20、−40、−60 ℃。每间隔3个工况为外界低温冷空气温度相同，室温温度分别为1、3、5、7、10 ℃。

各工况冻结状态计算结果见图8。分析表3可以看出，船体内部温度相同，外界空气温度的变化对压载舱内降温冻结有明显的影响，外界冷空气温度越低，结冰速率越快，冻结率越高；外界冷空气温度越高，结冰向下和向内扩展速率越慢。对比工况1、工况4、工况7、工况10、工况13这5组实验可以看出，当外界冷空气温度同为−20 ℃时，室内温度越高，结冰速率越低，越难以形成完整的冰盖，即内侧c点附近水域越难结冰，表明船体内部气温升高可以延缓结冰速率。

工况5、工况6、工况9沿船长方向在液面c点附近的水域未结冰，而在c点下方靠近H₁和H₂舱室隔板位置的水域发生冻结，其原因为c点位于压载水液面最接近内侧纵壁的位置，其温度变化主要受到舱室内部气温影响。压载水上下舱室间的平台隔板也是钢板，导热性较水和空气更强，平台隔板在横向沟通了舷侧外板和内侧纵壁，c点下方水域越靠近平台隔板的位置，其热量越容易通过钢制隔板向舷外低温环境传递，因此降温速率反而要比远离隔板的c点液面附近水域更快。因此，c点处液面水域结冰会相对于液面以下尤其是靠近平台隔板的水域更慢，呈现出未完全结冰的状态。

监测点a、c点温度-时间变化曲线如图9所示，图9(a)和图9(b)为舱外空气温度−40 ℃情况下，船体内部不同温度工况a、c点温度随时间变化曲线；图9(c)和图9(d)为船体内部温度5 ℃情况下，舱外气温不同工况a、c点温度随时间变化曲线。由图可见，舱外低温冷空气温度越低，a点、c点的整体降温速率越快，最终温度越低。船体内部的室温越高，a点、c点整体降温速率越慢，最终温度越高。但由于a点位于舷侧，主要受到外界低温冷空气得控制作用，船体内部气温度变化对于a点初始降温速率的影响较小，同理，c点主要受到船体内部室温控制作用，外界冷空气温度的变化，对于c点初始时刻降温结冰速率影响较小。

4 结　语

本文以极地船舶舷侧压载舱为研究对象，采用水冰相变计算方法和多物理场模拟方法，建立了舷侧压载舱结冰数值仿真模型，开展了压载水低温冻结过程研究，得到以下主要结论：

1）在极地低温环境下，压载水舱冻结区域首先出现在靠近舷侧外板的液面，随冻结时间增加，冰体发生横向和垂向生长，向内侧纵舱壁逐渐发展。

2）船体外界低温冷空气是影响压载水舱内结冰速率、温度场分布和结冰形态的最主要因素；舱内液位高度的变化影响压载水的降温冻结速率，合理调控压载水量是减轻压载水舱冻结情况的有效措施。