多级闪蒸是海水淡化工业中最成熟，运行安全性最高，弹性大，适合于大型和超大型淡化装置的一项技术^[1]。多级闪蒸过程中的循环盐水经过预热器加热，回收闪蒸淡水蒸汽的热量，温度升高，在流经换热器的过程中，某些组分的溶解度可能达到过饱和，导致结晶析出，在管道里形成析晶污垢^[2-3]，沉积于换热表面。污垢的存在不仅增大流体阻力，增加能量消耗^[4]，降低传热效率，严重时还会堵塞管路，使设备瘫痪。据分析海水析晶污垢的主要成分是CaSO₄·2H₂O^[5]，质量分数占到60.4%，其次为Mg(OH)₂，质量分数占到39.6%。因此，探讨CaSO₄·2H₂O析晶污垢在换热管道里的形成过程很有必要。

Fahmi Brahim^[6]在分析传质和表面反应的基础上，建立了污垢沉积率的数学模型，探讨了CaCO₃在换热面上结垢过程的温度场和浓度场的分布规律，并依据温度场和浓度场计算CaCO₃污垢的生成率、剥蚀率和污垢热阻等，但由于实际上污垢同时含有H₂O和CaCO₃，该研究没有考虑H₂O分子对析晶污垢的影响。E.Nebot等^[7]通过改进Konak模型，考察了流体流速对冷凝器中污垢沉积的影响，并指出增大流速有利于降低污垢热阻。程浩明等^[8]应用FLUENT构建了Hasson^[9]模型，从污垢沉积质量出发，计算出污垢热阻随时间的变化规律，并对比了光滑管和波纹管内CaCO₃污垢沉积和剥蚀的规律，结果表明波纹管具有较好的抑垢特性。该研究将结垢模型从直管拓展到复杂管路，使模型的适用性得到加强，但Hasson模型涉及参数较多，形式复杂，求解困难，有待进一步改进。Feng Zhang等^[10]使用MATLAB和多物理场耦合分析软件，将污垢形成过程以虚拟动态的形式直观地呈现出来。

目前，对CaSO₄析晶污垢形成过程的模拟已有报道，但存在如下问题：一方面，大多研究是将主流体到热边界层的传质过程和垢层表面发生的反应结晶过程进行综合处理，很少有研究直接使用表面反应速率计算结晶率，且忽视了真实状况下CaSO₄析晶污垢以CaSO₄·2H₂O的形式存在，即析晶污垢内CaSO₄与H₂O分子的摩尔数之比为1:2，未考虑H₂O分子对结垢速率的影响。另一方面，忽视了结垢过程中流动区域、主体流速以及各组分浓度场的变化。随着垢层厚度不断增加，流体流动区域会随之发生改变，这将导致网格划分区域发生变化，在模拟计算中需要对主体流速进行修正，同时还应考虑析晶后各组分在流动场中浓度的变化。应该采用动网格进行模拟计算，但计算量大，且难以收敛。针对以上问题，本文拟通过在传质方程中添加组分源项对析晶结垢而导致的流体主体中各组分浓度的变化进行修正，建立化学反应模型。直接采用表面反应速率计算结晶率，进一步修正传统的Hasson模型，以减少计算量。同时在析晶污垢生成速率中添加H₂O分子的质量，使模拟结果更接近实际生产的情况。以CaSO₄的析晶结垢过程为例，运用提出的模型和求解方法探讨影响析晶污垢生成的流速、壁面温度因素，得到流速、温度对污垢生成过程的影响规律，为多级闪蒸海水淡化的设计及运行提供理论指导。

1 数学模型 1.1 物理模型

本文探讨光滑直管的析晶结垢过程。由于直管内流体流动、传热、结垢等特性参数在径向上的变化规律是相同的，而流体仅沿轴向方向流动，故可将圆管内流体流动简化为二维模型(见图 1)。在二维空间内建立析晶结垢的数学模型并进行求解，获得的速度、温度和浓度分布结果，既能够满足计算精度要求，又能清晰地反应问题的本质。

1.2 控制方程

CaSO₄析晶污垢的形成是一个典型的不可压缩、无内热源的二维对流换热问题。为了得到整个计算区域的速度场、温度场以及浓度场，需要将连续性方程、动量方程、能量方程以及传质方程联合求解，主要控制方程如下：

$ 连续性方程：\frac{{\partial \rho }}{{\partial t}} + \nabla \cdot \left( {\rho V} \right) = 0。$

(1)

动量方程：

$ \frac{{\partial \left( {\rho V} \right)}}{{\partial t}} + \nabla \cdot \left( {\rho VV} \right) = - \nabla P + \nabla \cdot \left( \tau \right) + \rho g。$

(2)

能量方程：

$ \begin{array}{l} \;\frac{{\partial \left( {\rho E} \right)}}{{\partial t}} + \nabla \cdot ((V(\rho E + p)) = \\ \nabla \cdot (k\nabla T - \sum {_{_j}} {h_j}{J_j} + \left( {\tau \cdot V} \right))。\end{array} $

(3)

溶液中发生的析晶反应：

$ {\rm{C}}{{\rm{a}}^{2 + }} + {\rm{SO}}_4^{2 - } \to {\rm{CaS}}{{\rm{O}}_4} \downarrow 。$

(4)

溶液中每种组分的浓度随着析晶结垢过程的进行将发生变化，反应过程的速率可以用反应物Ca²⁺和SO₄^2-减少的速率或生成物CaSO₄增加的速率来表示，反应速率常数^[11]可以通过阿伦尼乌斯定律计算得到：

$ {k_r} = {k_{r0}}\exp ( - \frac{E}{{R{T_f}}})。$

(5)

其中：k_r0=7.07 m⁴·kg^-1·s^-1；活化能E=37 143 J·mol^-1；T_f为垢层表面温度，K。

因此传质方程可表示为：

$ \frac{{\partial (\rho {Y_i})}}{{\partial t}} + \nabla \cdot (\rho V{Y_i}) = - \nabla \cdot {J_j} + {S_i}。$

(6)

其中Y_i表示溶质的质量分数。

针对CaSO₄组分，传质方程中的组分源项可表示为^[12-13]：

$ {S_i} = {k_r} \cdot {({c_f} - {c_s})^2}。$

(7)

其中：c_f为CaSO₄组分的浓度，kg·m^-3；c_s为CaSO₄饱和浓度，kg·m^-3，且c_s是垢层表面温度T_f的函数，可以通过溶解度曲线的拟合公式^[14]得到：

$ {c_s} = 1.46 \cdot {10^{ - 4}}{T^2} - 0.12T + 26.60。$

(8)

由于垢层成分为CaSO₄·2H₂O，即生成一个CaSO₄分子的同时需要消耗两个H₂O分子，CaSO₄分子和H₂O分子的消耗速率之比为1:2。因此，修正H₂O分子组分的源项为：

$ {{\rm{S}}_{{{\rm{H}}_{\rm{2}}}}}_{\rm{O}}{\rm{ = 2}} \cdot \frac{{{{\rm{M}}_{{{\rm{H}}_{\rm{2}}}{\rm{O}}}}}}{{{{\rm{M}}_{{\rm{CaS}}{{\rm{O}}_{\rm{4}}}}}}} \cdot {{\rm{S}}_{{\rm{CaS}}{{\rm{O}}_{\rm{4}}}}}。$

(9)

M_H2O为H₂O的摩尔分子质量，M_CaSO4为CaSO₄的摩尔分子质量。

1.3 结垢模型

CaSO₄·2H₂O析晶污垢在换热面上的结垢过程主要包括：壁面结晶和颗粒剥蚀。其生成速率取决于表面结晶反应的速率，产生的析晶污垢量可由下式计算：

$ \frac{{{\rm{d}}m}}{{{\rm{d}}t}} = {(1 + 2 \cdot \frac{{{{\rm{M}}_{{{\rm{H}}_{\rm{2}}}{\rm{O}}}}}}{{{{\rm{M}}_{{\rm{CaS}}{{\rm{O}}_{\rm{4}}}}}}} \cdot {k_r} \cdot ({c_F} - {c_S}))^2}。$

(10)

其中：m为单位面积污垢的生成量，单位为kg·m^-2；c_F为CaSO₄组分在垢层表面处的浓度，单位为kg·m^-3。k_r与c_s都是垢层表面温度T_f的函数，假设热流仅沿y轴方向一维传递，则垢层表面温度可表示为：

$ {T_f} = {T_b} + \frac{{{T_w} - {T_b}}}{{1 + \alpha \cdot {R_f}}}, $

(11)

$ \alpha = 0.023 \cdot \frac{\lambda }{d} \cdot {R_e}^{0.8} \cdot {R_r}^{0.4}。$

(12)

式中：T_b为流体主体温度，K；T_w为管壁温度，K；α为流体与垢层表面间的换热系数，W·m^-2·K^-1；R_f为污垢热阻，m²·K·W^-1；λ为溶液导热系数；d为管道直径，m；R_e为流体雷诺数；R_r为流体普朗特数。

污垢剥蚀率采用Brahim模型^[15]求解：

$ \begin{array}{l} \frac{{{\rm{d}}{m_r}}}{{{\rm{d}}t}} = \frac{K}{P} \cdot {\rho _f} \cdot \left( {1 + \beta \cdot \Delta T} \right) \cdot \\ \;\;{d_p} \cdot {({p^2} \cdot \eta \cdot g)^{1/3}} \cdot {x_f} \cdot {u^2}。\end{array} $

(13)

式(13)中：ρ_f为垢层密度，kg·m^-3；β为线胀系数，K^-1；d_p为溶液中悬浮的CaSO₄·2H₂O晶体的平均粒径，m；η为流体的动力粘度，Pa·s；x_f为垢层厚度，m；u为流体主体流速，m·s^-1；P为晶间附着力，K表征垢层内缺陷点数目，根据文献[16]P/K采用下述公式：

$ \frac{P}{K} = 83.2 \cdot {u^{0.54}}。$

(14)

污垢净生成速率为析晶污垢生成速率与剥蚀速率的差，可由下式计算：

$ \frac{{{\rm{d}}{m_f}}}{{{\rm{d}}t}} = \frac{{{\rm{d}}m}}{{{\rm{d}}t}} - \frac{{{\rm{d}}{m_r}}}{{{\rm{d}}t}}。$

(15)

单位面积上总污垢生成量m_f为每个时间步内污垢净生成量的累积：

$ {m_f} = {m_{f, t + \mathit{\Delta }t}} = {m_{f, t}} + \frac{{{\rm{d}}{m_f}}}{{{\rm{d}}t}}\cdot\Delta t。$

(16)

单位面积总的垢层厚度为x_{f, t+Δt}每个时间步内污垢累积的厚度：

$ {x_{f, t + \mathit{\Delta }t}} = {x_{f, t}} + \frac{{{m_f} \cdot \Delta t}}{{{\rho _f}}}。$

(17)

单位面积总污垢热阻R_{f, t}为总垢层厚度与垢层导热系数的比值：

$ {R_{f, t}} = \frac{{{x_{f, t}}}}{{{\lambda _f}}}。$

(18)

λ_f为CaSO₄·2H₂O析晶污垢的导热系数，W·m^-1·K^-1。

1.4 流速的修正以及变量的初始化

由1.3小节所给出的结垢模型可知，当CaSO₄·2H₂O析晶污垢在管道内生成时，管内垢层厚度随时间逐渐增加，流体的流道会减小，流速随之增加，即在轴向方向上，管内的流速是变化的。因此，计算过程中必须对流速进行修正。本文提出如下关联式进行流速的修正，以更接近管内的实际流动状况：

$ u = {u_0} \cdot \frac{{{d^2}}}{{{{(d - 2{x_{f, t}})}^2}}}。$

(19)

计算时将3 600 s作为非稳态计算的时间步，即Δt取3 600，在t=0时刻，垢层表面温度取壁面温度，污垢剥蚀率、垢层厚度、污垢热阻均为零。

2 数值计算结果 2.1 模拟条件

入流口的速度充分发展，平均速度为v，流体温度为主体温度。壁面无滑移且温度恒定。出流口的法向速度分布满足整体质量守恒。实际流速和垢层表面温度取修正过的值，污垢结晶率、剥蚀率和净生成率取处监测面上的平均值。在相同模型与边界条件下，对流域进行不同尺寸的网格划分，得到CaSO₄组分的质量分布基本一致，证明模型结果与网格数无关，最终选网格尺寸为0.5 mm。

本文采用标准的湍流模型对CaSO₄-H₂O体系进行计算。由于流体浓度很低，所以流体物性近似取壁面温度条件下水的物性^[11]，分别讨论流速、壁面温度等因素对析晶污垢的影响。

2.2 流速对CaSO₄·2H₂O污垢生成速率的影响

以CaSO₄水溶液在换热器内的流动为例，流速的变化将影响CaSO₄·2H₂O在换热器管壁的生成率和剥蚀率。模拟条件为CaSO₄的质量分数ω=0.000 1，流体温度T_b=300 K，壁面温度a320 K。分别取流速的值为：0.3，0.5，1.0，1.5，2.0 m/s，模拟结果如下：

由图 2和图 4可得，流速越大垢层表面温度越低，垢层热阻越大；随着运行时间的延长，垢层表面温度逐渐降低，垢层热阻逐渐升高，在大约150 h以后，由于污垢净生成率降为零，污垢热阻不再增加，垢层表面温度趋于稳定。其原因是在流动开始时，壁面处于洁净状态，垢层表面温度等于壁面温度，流体主体与壁面的温差大，溶液的过饱和度增大，净生成率较大；随着时间的推移，污垢逐渐累积，污垢热阻增大，垢层表面温度逐渐降低，净生成率逐渐减小；时间进行越久，污垢结晶速率和剥蚀速率达到动态平衡，最终净生成速率趋于零，垢层厚度不再增大，污垢热阻(见图 4)趋于稳定。由图 3可知，流速增大有利于降低CaSO₄·2H₂O析晶污垢的净生成率。在流速较低的条件下，对流传热系数较小，传热热阻大，因此越靠近壁面处，温度越高，CaSO₄·2H₂O污垢的结晶速率越快。随着流速越大，污垢净生成率降低的幅度逐渐减小，当流速超过1.5 m/s后，增大流速对降低污垢净生成率的效果不再明显，而流速过高会导致较高的能耗，综合考虑，流速的适宜取值范围为1.0~1.5 m/s。

通过将本文模拟结果与Brahim等^[6]的实验结果进行对比可知，在模型条件与实验条件相近的情况下，污垢热阻的计算数值与实验数值吻合较好且变化趋势相同，证明了本文建立的模型及模拟结果的可靠性。

2.3 换热管壁面温度对CaSO₄·2H₂O污垢生成速率的影响

模拟条件为：CaSO₄的质量分数ω=0.000 3，流体温度T_b=300 K，流速u=1.0 m/s。壁面温度T_W分别取：310，320，330，340，350 K，模拟结果如下：

由图 6可看出，在污垢析出的过程中，温度越高，越有利于污垢的生成，且温度对于污垢析出的影响较流速而言更为显著(见图 7)，随时间的推移，污垢净生成率逐渐减小。由图 8可知，壁面温度越高，污垢热阻增加的速率越快，污垢生成总量随时间逐渐增加，导致主体流速增大，流体对污垢的剥蚀速率增大，且壁面温度越高，剥蚀速率越大。由图 9可看出，对于不同温度，管道内主体流速都呈现逐渐增大最终趋于恒定的变化趋势，壁面温度越高，流速达到恒定时所用的时间越短，且流速越大。

2.4 CaSO₄·2H₂O污垢的生成速率在管道不同内径处的变化

模拟条件为：流速u=1.0 m/s，CaSO₄的质量分数ω=0.000 3，流体温度T_b=300 K, 壁面温度T_W=320 K。在x=0.5 m处，取不同半径处的圆环面为假想壁面，得到不同假想壁面上污垢的净生成率和剥蚀率，结果如下图所示。

由图 10、11、12可知，越靠近管中心，即r越小，流速越大，析晶污垢净生成率越小，污垢热阻越小。该生成规律符合析晶污垢沉积扩散的方向，即析出晶体从流体中不断向壁面聚集，最终导致壁面垢层的生长，垢层热阻的增大。

2.5 CaSO₄·2H₂O污垢的结晶速率和剥蚀速率随位置的变化

模拟条件为：流速u=1.0 m/s，CaSO₄的质量分数ω=0.000 3，流体温度T_b=300 K, 壁面温度T_W=320 K。在沿管道轴向的不同位置处，即x取不同值，得到不同位置上的流速、垢层厚度、污垢热阻和污垢净生成率，结果如下图所示。

由图 13(a)(b)看出，随着时间的延长，管道轴线方向不同位置处的流速都呈现出先增大后逐渐趋于平稳的变化趋势，且离入口处越远流速越大；在管道轴向的同一位置，随着时间的延长，流速也呈现出逐渐增大的趋势，且初始时间段内变化较大(t=40~120 h)，而后逐渐趋于稳定(t=160~200 h)。

由图 14(a)看出，沿管道轴线方向的不同位置，垢层的厚度都呈现出先增大，而后趋于稳定的变化趋势；在管道的同一位置，随着时间的延长，垢层厚度有一个先快速增大，而后增加速率逐渐减小，以后，垢层厚度不再增加。垢层热阻的变化趋势与垢层厚度的变化趋势相同(见图 15)，也呈现出先增大后趋于稳定的变化过程，在时以后，垢层热阻不再增加。

由图 16可知，管道内不同位置的污垢净生成率在起始时刻最大，但随时间逐渐减小，150 h后，垢层不再生长，即沉积速率与剥蚀速率达到了平衡。管子入口处净生成速率最小，沿管长方向净生成速率逐渐增大并趋于稳定。

3 结语

本文运用FLUENT软件对多级闪蒸海水淡化系统中析晶污垢形成过程进行了数值模拟。通过添加组分源项，修正流体中组分浓度的分布，采用修正后的流体主体浓度替代垢层表面处的浓度，对传统的Hasson模型进行改进；同时考虑H₂O分子对析晶污垢生成速率的贡献，使模拟结果更符合实际情况。以CaSO₄·2H₂O污垢形成过程为例，分析了流速、壁面温度对污垢形成速率的影响以及管道径向和轴向方向的析晶规律，结果表明：管内流速越大，污垢的净生成率和热阻越小，但流速越高流体流动造成的能耗越大，综合考虑结垢和能耗因素，适宜的流速范围为1.0~1.5m/s；污垢生成速率随壁面温度的增大而显著增大；污垢净生成速率沿径向方向逐渐增大，而沿轴向方向，管子入口处污垢净生成率最小。

符合说明

c_f  溶质的组份浓度(kg·m^-3)

c_F  CaSO₄组分在垢层表面上的浓度(kg·m^-2)

c_s  溶质的饱和浓度(kg·m^-3)

d  管径(m)

d_p  溶液中悬浮CaSO₄·2H₂O晶体的平均粒子直径(m)

E  活化能(J·mol^-1)

K  垢层内缺陷点的数目

k_r  反应速率常数(m⁴·kg^-1·s^-1)

k_r0  指前因子(kg·m^-2)

m  单位面积污垢的生成量(kg·m^-2)

m_r  单位面积污垢的剥蚀量(kg·m^-2)

m_f  单位面积污垢的净生成量(kg·m^-2)

P  晶间附着力

R_r  流体普朗特常数

R  通用气体常数(J·mol^-1·K^-1)

R_e  流体雷诺数

R_f  垢层热阻(m²·K·W^-1)

R_{f, t}  污垢总热阻(m²·K·W^-1)

T_b  流体主体温度(K)

T_f  垢层表面温度(K)

T_w  管壁温度(K)

u  流体主体流速(m·s^-1)

x_f  垢层厚度(m)

x_{f, t+Δt}  垢层总厚度(m)

α  传热系数(W·m^-2·K^-1)

β  线性膨胀系数(K^-1)

λ  溶液热导率(W·m^-2·K^-1)

λ_f  CaSO₄·2H₂O结晶垢的热导率(W·m^-2·K^-1)

η  动力粘度(Pa·s)

ρ_f  垢层密度(kg·m^-3)

Nomenclature

c_f  The component concentration(kg·m^-3)

c_F   The concentration of CaSO₄ component at fouling surface(kg·m^-2)

c_s  The component saturation concentration (kg·m^-3)

d  The diameter of tube (m)

d_p  The average particle size of crystal CaSO₄·2H₂O suspended in solution (m)

E  The activation energy(J·mol^-1)

K  The number of defect points inner fouling layer

k_r  The reaction rate constant (m⁴·kg^-1·s^-1)

k_r0  The pre-exponential factor(m⁴·kg^-1·s^-1)

m  The total generation quantity of crystallization fouling(kg·m^-2)

m_r  The erosion rate of crystallization fouling(kg·m^-2)

m_f  The net generation rate(kg·m^-2)

P  The adhesive force between crystals

P_r  The Prandtl number of fluid

R  The universal gas constant(J·mol^-1·K^-1)

R_e  The Reynolds number of fluid

R_f  The fouling resistance (m²·K·W^-1)

R_{f, t} The total fouling resistance(m²·K·W^-1)

T_b  The bulk temperature of fluid (K)

T_f  The surface temperature of fouling layer (K)

T_w  The tube wall temperature (K)

u  The bulk velocity of fluid(m·s^-1)

x_f  The thickness of fouling layer (m)

x_{f, t+Δt}  The total thickness of fouling layer (m)

α  The coefficient of heat transfer (W·m^-2·K^-1)

β  The coefficient of linear expansion (K^-1)

λ  The thermal conductivity of solution(W·m^-2·K^-1)

λ_f  The heat conductivity of CaSO₄·2H₂O crystallization fouling (W·m^-2·K^-1)

η  The dynamic viscosity of fluid(Pa·s)

ρ_f  The density of fouling layer(kg·m^-3)