地表水源热泵作为一种清洁的可再生能源利用技术，近年来受到广泛关注[1]．传统地表水源热泵系统易受到水源温度影响，我国北方地区冬季地表水温过低，蒸发器面临结冰风险[2]．为解决这一问题，学者提出了冰源热泵技术[3]，在传统地表水源热泵蒸发器侧增加乙二醇溶液循环，利用过冷水低于0 ℃不结冰的特性[4]，使用超声波控制过冷水在制冰装置内解除过冷状态，形成含有冰晶颗粒的冰浆，再通过出冰管道排放到外界环境中．冰源热泵技术扩大了地表水源热泵的使用范围，并大幅度减少用水量，减少对环境的影响[5]．冰源热泵系统使用过冷水法制取冰浆，在制冰装置中，超声波可以使过冷水内部产生超声空化现象[6]，形成微小的气泡．气泡崩溃时瞬间在气泡周围产生较大的温度与压力梯度，降低液体的冻结温度导致成核[7]．产生的冲击波可使冰晶粒径更小、分布更均匀，能有效地防止冰晶在管道内的沉积．低浓度冰浆过滤后，经由管道直接排出到外界环境．冰浆在出冰管道内的安全运输直接影响系统性能稳定性，若管道发生冰堵，则须暂停系统融冰，因此保证管道内冰浆流速高于形成冰堵的临界流速是系统连续安全运行的关键．目前国内外学者已对管内冰浆流动进行了多项研究，主要集中于压降变化、流态变化、换热特性及临界流速等方面．文献[8]对乙二醇冰浆在过渡区的流动和传热特性进行了试验研究，发现：在低含冰率区域管道摩擦系数和传热系数是恒定的，而在高含冰率区域会出现分层现象．文献[9]对乙二醇冰浆在水平加热管中的换热情况进行了数值研究．结果表明：在层流进口条件下，径向传热是通过热传导，熔化只发生在管壁附近．在湍流进口条件下，熔化发生在垂直于管道轴线的整个横截面上．文献[10]介绍了一种计算冰浆流体沉积速度的新方法，可用于丙二醇溶液不同浓度和管道直径的冰浆沉积速度计算．前人的研究主要集中于添加乙二醇等冰点抑制剂的冰浆，对非绝热条件下纯水冰浆及其管道内沉积的临界流速研究较少，本研究结合两种临界流速判定标准，对管道内不同初始含冰率(ice packing fraction，IPF)冰浆在非绝热条件下形成冰堵的临界流速进行研究．1 冰浆管道换热模型的建立1.1　建立物理模型利用ICEM软件对出冰管道进行物理建模和网格划分，管道模型有水平直管和90°弯管两种．水平直管基准模型见图1，图中：长度L1=5 m；管径D1=0.1 m．经验证，当网格总数达到2.808×105后，模拟结果随网格数目增长变化可忽略，满足网格独立性．网格总体质量均在0.7以上，满足仿真软件Fluent的要求．90°弯管基准模型见图2，图中：弯管管道轴线长度L2=5 m；弯管管径D2=0.1 m；弯管曲率半径R=0.3 m，弯曲部分位于管道中间．图2模型同样满足网格独立性．10.13245/j.hust.221011.F001图1水平直管模型10.13245/j.hust.221011.F002图290°弯管模型模拟试验中将冰晶颗粒视为球体，结合两种临界流速判断依据，分别研究冰晶粒径、管道管径、管长及弯管曲率半径对不同IPF冰浆在管道内冰堵临界流速的影响．现有过冷法制取冰浆的冰晶粒径ds可达到0.1~0.3 mm[11]，冰浆中水和冰晶的物性参数如表1所示，表中：ρ为密度；μ为运动黏度．10.13245/j.hust.221011.T001表1冰浆物性参数材料ρ/(kg∙m-3)μ/(mPa∙s)ds/mm水999.91.72×10-2冰晶915.01.790.1，0.15，0.20，0.251.2　数值模型本研究选择RNG k-ε模型和欧拉模型[12]对水平直管内冰浆流动特性进行模拟仿真．欧拉模型是一种研究多相流动的常用模型，被广泛应用在颗粒悬浮流、流化床等领域的模拟研究中．欧拉模型将多相流描述为互相穿插的连续体，分别建立各相的Navier-Stokes方程，各相均满足质量守恒定律、动量守恒定律和能量守恒定律[13]．第q项的体积Vq=∫αqdV，式中：αq为第q相的体积分数，各项体积分数之和为1，即∑αq=1；V为多相流体积．1.2.1　连续性方程第q相的连续性方程为∂∂t(αqρq)+∇(αqρqvq)=∑p=1n(m˙pq-m˙qp)+Sq，式中：t为时间；ρq为第q相的密度；vq为第q相的速度；m˙pq为p相到q相的传质；m˙qp为q相到p相的传质；Sq为质量源项．1.2.2　动量方程第q相的动量方程可描述为[14]∂∂t(αqρqvq)+∇(αqρqvqvq)=-αq∇P+∇∙τq+αqρqg+∑p=1n(Rpq+vpqm˙pq-vqpm˙qp)+(Fq+Flq+Fvq+Ftq),式中：g为重力；Fq为体积力；Flq为浮升力；Fvq为虚拟质量力；Ftq为湍流扩散力；Rpq为第p相和第q相之间的相互作用力，P是所有相的共同压力；vpq为p相与q相之间的相对速度；τq为第q相的应力应变张量，有τq=αqμq(∇vq+vqT)+αq(λq-2μq/3)∇⋅vqI，其中，μq和λq为第q相的剪切黏度和体积黏度，I为单位矢量．1.2.3　能量方程对于液相和固相的能量守恒方程为      ∂∂t(αqρq(eq+vq2/2))+∇⋅(αqρqvq(hq+vq2/2))=∇⋅αqkeq∇Tq-∑jhj,qJj,q+τeq⋅vq+∑p=1n(Qqp+m˙pqHpq-m˙qpHqp)+P∂αq∂t+Sq，式中：eq为q相的势能；hq为q相的焓；τeq为q相的有效应力应变张量；Tq为q相温度；keq为有效导热系数；Sq为源相；Qqp为q相和p相之间的换热强度；Hpq为p相转变为q相的相变潜热；hj,q为q相中j组分的焓；Jj,q为q相中j组分的扩散通量．相之间的热交换符合局部平衡条件Qqp=-Qpq和Qqq=0．1.2.4　相间传热传质模型由于欧拉模型不能与Fluent软件内置的换热模型直接耦合，因此须要通过编写用户自定义方程(user defined function，UDF)内嵌．冰晶颗粒与水之间的相间传热系数Kls使用文献[15]模型计算，有Kls=6λlαlαsNu/ds2，式中：λl为液相导热系数；Nu为努塞尔数，αl为液相体积分数；αs为固相体积分数；ds为固相颗粒直径．又有Nu=(7-10αs+αs2)(1+0.7Re0.2Pr1/3)+(1.33-2.4αs+αs2)Re0.2Pr1/3. (1)式(1)中，Pr为普朗特数，有Pr=cp1μl/λl，式中：cp1为液相定压比热容；μl为液相动力黏度．式(1)中，Re为雷诺数，有Re=ρlds|vl-vs|/μl，式中：ρl为液相密度；vl为液相流速；vs为固相颗粒流速．冰晶由于融化向水传递的质量mls=Kls(Tl-Ts)/ΔH，式中：Tl为液相温度；Ts为固相颗粒温度；ΔH为固相冰晶的潜热．1.2.5　湍流模型本研究采用RNG湍流模型，描述各相湍流模型的湍流脉动动能k和湍流耗散率ε方程分别为∂∂t(ρk)+∂∂xi(ρkvi)=∂∂xjβkμe∂k∂xj+Gk+Gb-ρε-YM+Sk；∂∂t(ρε)+∂∂xi(ρεvi)=∂∂xjβεμe∂ε∂xj+C1ε(Gk+C3εGb)/k-C2ερε2/k-χε+Sε,式中：vi为第i相速度；μe为有效黏度；Gk为平均速度梯度产生的湍流动能；Gb为浮力产生的湍流动能；YM为可压缩湍流中脉动膨胀对总耗散率的作用；βk为k的逆有效普朗特数；βε为ε的逆有效普朗特数；Sk为k的源项；Sε为ε的源项；C1ε，C2ε和C3ε为系数，其中C1ε=1.42；C2ε=1.68；xi为第i相x方向的分量；xj为第j相x方向的分量；χε为适应应变率和流线曲率变化的附加项．1.3　边界及初始条件入口边界设置为速度入口，水和冰晶以均匀流速进入管道，设置冰晶颗粒在管道入口截面内均匀分布，管道初始含冰率IPF即指入口边界的平均含冰率，范围为1%~10%．入口边界设置湍流强度和水力直径，水力直径为管道直径．出口边界设置为压力出口．水与冰晶的温度在入口边界设置为0 ℃，管壁面处设置为恒定温度15 ℃．使用标准壁面函数处理近壁面的流动，近壁面第一层网格尺寸满足壁面率y+值约为30．1.4　临界流速判断依据管道是否发生冰堵有两种判断依据：管道内的冰晶颗粒最低速度为0 m/s，对应临界流速一(v1)；管道内局部最大含冰率大于52%，对应临界流速二(v2)[16]．第二种判断依据是由球体在三维空间的不同堆积方式确定的，冰浆在流动过程中有悬浮流、流化床和固定床[17]三种流态．当含冰率大于52%时，冰浆流态由悬浮流转变为流化床，容易造成管道冰堵，两种流态转变速度对应临界流速二．此判断依据参考了文献[16]．1.5　模型验证为了确保模拟仿真结果的可靠性，须要验证本研究中使用的数值模型．文献[16]对管道内冰浆固液两相流动进行了一系列研究，通过将模拟数据和实测数据进行对比，具有很高的一致性，其内容与本研究基本原理相同，试验条件设置具有较高相似性．表2显示了模拟结果和论文模拟数据的对比，物理参数为管径DN50，乙醇溶液体积分数为10%，ds=0.3 mm．本文模型模拟得到的结果与文献中的数据绝对误差在±0.01 m/s以内，相对误差为相对±10%以内，可以得出模拟结果与文献[16]中的数据有较高的一致性，验证了本文所用数值模型的可靠性．10.13245/j.hust.221011.T002表2模型验证模型IPF值/%510152030文献[16]0.0900.1100.1200.1200.130本文0.1000.1170.1260.1300.141m·s-12 非绝热条件冰浆在管道内流动特性分析2.1　水平直管内冰浆流动临界流速2.1.1　粒径对临界流速的影响图3和图4显示了IPF值在1%~10%的范围内，ds对管内冰浆流动两种临界流速的影响．由图3可以得出：v1随着冰晶粒径增加显著增长，增长率也不断扩大，v1的变化范围为0.95~1.35 m/s．由图4可得出：v2随着冰晶粒径增加呈近似线性增长趋势，在IPF值为7%~10%的范围内更为明显，且增长率相近，v2的变化范围为0.31~0.49 m/s．可以得出： 冰晶粒径对v1的影响显著，对v2的影响相对较小．10.13245/j.hust.221011.F003图3ds对v1的影响10.13245/j.hust.221011.F004图4ds对v2的影响2.1.2　管径对临界流速的影响图5和图6显示了IPF值在1%~10%的范围内，直管管径D1对管内冰浆流动两种临界流速的影响．由图5可以得出：v1随着直管管径的增加呈现增长趋势，增长幅度逐渐平缓，v1的变化范围为0.07~0.10 m/s．由图6可以得出：v2不断增长，在不同含冰率下增长率相近，线性特征明显，v2的变化范围为0.05~0.11 m/s．可以得出：直管管径对两种临界流速的影响较小，影响程度大致相当．10.13245/j.hust.221011.F005图5D1对v1的影响10.13245/j.hust.221011.F006图6D1对v2的影响2.1.3　管长对临界流速的影响图7和图8显示了IPF值在1%~10%的范围内，直管管长L1对管内冰浆流动两种临界流速的影响．由图7中可以得出：v1随着直管管长的增加，呈现显著上升趋势，同时增长率缓慢增加，v1的变化范围为1.04~1.62 m/s．由图8可以得出，v2呈近似线性增长趋势，增长率随IPF值增加略有上升，v2的变化范围为0.14~0.44 m/s．可以得出：直管管径对v1的影响显著，对v2的影响较弱．10.13245/j.hust.221011.F007图7L1对v1的影响10.13245/j.hust.221011.F008图8L1对v2的影响2.2　90°弯管内冰浆流动临界流速2.2.1　粒径对临界流速的影响图9和图10显示了IPF值在1%~10%的范围内，冰晶粒径对90°弯管内冰浆流动两种临界流速的影响．由图9中可以得出：v1随着冰晶粒径的增加，呈现显著增长趋势，且增长率略有上升，与直管内的现象类似，v1的变化范围为0.37~0.57 m/s．由图10可以得出：v2随冰晶粒径增加呈增长趋势，增长率逐渐下降，v2的变化范围为0.20~0.36 m/s．可以得出：冰晶粒径对v1的影响较大，对v2的影响相对较弱．10.13245/j.hust.221011.F009图9ds对90°弯管v1的影响10.13245/j.hust.221011.F010图10ds对90°弯管v2的影响2.2.2　管径对临界流速的影响图11和图12显示了IPF值在1%~10%的范围内，弯管管径D2对管内冰浆流动两种临界流速的影响．由图11中可以得出：v1随着弯管管径的增加，呈现微弱上升趋势，增长率较小且逐渐降低，v1的变化范围为0.02~0.04 m/s．由图12可以得出：v2同样呈现微弱上升趋势，增长率较小且逐渐趋于平缓，v2的变化范围为0.04~0.07 m/s．可以得出：弯管管径对临界流速的影响较小，且对v2的影响略大于对v1的影响．10.13245/j.hust.221011.F011图11D2对90°弯管v1的影响10.13245/j.hust.221011.F012图12D2对90°弯管v2的影响2.2.3　曲率半径对临界流速的影响图13和图14显示了IPF值在1% ~ 10%的范围内，弯管曲率R对管内冰浆流动两种临界流速的影响．由图13中可以得出：v1随着弯管曲率的增加，呈现微弱下降趋势，v1的变化范围为0.02~0.03 m/s．由图14可以得出：v2呈近似微弱下降趋势，v2的变化范围为0.01~0.03 m/s．可以得出：弯管曲率对临界流速的影响微弱，且对两种临界流速的影响大致相当．10.13245/j.hust.221011.F013图13R对90°弯管v1的影响10.13245/j.hust.221011.F014图14R对90°弯管v2的影响3 结论本研究利用CFD软件Fluent对水平直管和90°弯管内IPF值为1%~10%的冰浆流动情况进行了模拟，并获得冰晶粒径、管径、管长、弯管曲率等参数对两种临界流速的影响，结果如下．a. 冰晶粒径对水平直管和90°弯管的两种临界流速的影响显著．对于水平直管，v1随着冰晶粒径的增加，上升趋势明显，增长率也不断扩大，v2呈近似线性增长趋势．对于90°弯管，v1随着冰晶粒径的增加不断增长，且增长率不断上升，与直管内的现象类似；v2呈增长趋势，且增长率随IPF值增加有所下降．冰晶粒径对v1的影响远大于对v2的影响，对水平直管的影响大于对90°弯管的影响．b. 管道管径对水平直管和90°弯管的两种临界流速的影响较小．对于水平直管，v1随着直管管径的增加，呈现显著上升趋势，增长率逐渐趋于平缓；v2呈近似线性增长趋势，不同IPF值下增长率相近．对于90°弯管，两种临界流速随着弯管管径的增加，呈现微弱上升趋势，增长率较小且逐渐降低，管道管径对两种临界流速的影响程度大致相当．c. 直管管长对管内冰浆流动两种临界流速的影响较大，且对v1的影响远大于对v2的影响．v1随着直管管长的增加，呈现显著上升趋势，增长率逐渐增加．v2呈近似线性增长趋势，且增长率随IPF值增加有所上升．d. 弯管曲率对管内冰浆流动两种临界流速的影响微弱，对两种临界流速的影响大致相当，临界流速随着弯管曲率的增加呈现微弱下降趋势．