基于最佳预后元模型的颗粒污垢特性全局敏感性分析

2024-02-23谢广烁张斯亮何松肖娟王斯民

化工进展 2024年1期

谢广烁，张斯亮，何松，肖娟，王斯民

（1 西安交通大学化学工程与技术学院，陕西西安 710049；2 上海蓝滨石化设备有限责任公司，上海 201518）

换热器广泛应用在化工、石油、电力、食品等工业生产过程中，是关系到国计民生的重要设备[1-4]。换热管作为管壳式换热器的主要元件，污垢积累会恶化传热性能，增加流动阻力，降低使用寿命。Akinç 等[5]发现炼油厂内换热器结垢造成的能源损失占总能耗的2%。颗粒污垢是换热器内主要的污垢类型之一，不仅会增加析晶污垢的成核数量，还会为生物的聚集、生长、繁殖提供场所，使多种污垢并存而难以去除[6-7]。

众多学者分析了不同操作工况对颗粒沉积特性的影响[8-10]。刘洪涛等[11]研究了流动方式和流速对壁面上颗粒沉积速率的影响，发现水平流动下沉积速率随颗粒量纲为1弛豫时间的增加而增加。张宁等[12]模拟了圆管内的颗粒污垢生长特性，研究了流速、入口温度和管壁粗糙度对污垢热阻的影响。于晓燕[13]计算了圆管内浓度、流速和温度对颗粒沉积的影响，结果表明浓度和流速对污垢特性影响较为明显。以上研究分析了常见操作参数对颗粒沉积的影响，如颗粒直径、浓度、流速和温度等，但并未对各参数的影响程度进行定量描述。

敏感性分析作为定量评估模型或系统的输入变量对输出变量影响程度的重要方法，其中全局敏感性分析可分析每个变量及其交互作用对输出变量的影响，量化设计变量的不确定性并有效识别关键参数，避免优化无关参数导致有限的性能提升及计算负荷浪费[14-16]。Du 等[17]采用Taguchi 方法研究了螺旋角、管径、管中心距等几何参数对管壳式换热器流动和传热性能的影响。Fesanghary 等[18]基于Sobol'全局敏感性量化了几何参数对折流板换热器总成本的影响，结果表明密封条数量是最不重要的参数，在优化过程中将其剔除。基于元模型的敏感性分析具有计算量小、与经验无关等优点，得到广泛应用。元模型是一种表示输入与输出变量关系的近似数学算法，常用的元模型有多项式回归、克里金和遗传聚合模型等。由于工程问题中输入与输出变量关系的复杂性，元模型决定着敏感性分析的准确性，因此最佳预后元模型（metamodel of optimal prognosis，MOP）应运而生，从模型库中自动选择合适的元模型并处理变量，具有较高的通用性和准确性[19]。Xiao 等[20]采用基于MOP 的全局敏感性分析，研究了折弯角、折弯率、相对高度和壳程流速对旋梯式螺旋折流板换热器流动传热性能的影响，结果表明壳侧流速和折弯角对单位压降传热系数影响较为明显，总影响分别为70.34%、25.01%。

内插件不仅增强管内传热，还具备在线防污除垢的作用，根据不同的运动方式分为固定式和转动式。自转式内插件由于流体流过产生动力矩使其旋转，管内压降通常低于静止内插件的压降[21]，且比电机带动下的主动式旋转内插件减少了额外动力消耗。因此，本文以内置转子防结垢为研究背景，基于MOP 对传热管内颗粒污垢特性进行全局敏感性分析。颗粒直径、入口流速、颗粒浓度和入口温度为输入变量，沉积率和污垢热阻为输出参数，并且分析对扩散沉积、湍流泳沉积、热泳沉积和重力沉积四个主要沉积速率的影响。训练数据空间通过数值计算获得，采用了欧拉-欧拉模型和颗粒污垢模型，得到了输入参数对输出参数的主影响和交互影响，为传热管内防结垢找寻关键变量。

1 数值计算方法与几何模型

1.1 几何模型

如图1所示，传热管内径28mm，长度300mm；转子由一个六边形面扭转90°得到，并在径向最大处两端各倒出一个半径10mm 的圆角。转子直径26mm，节距78mm，为了方便转子在实际应用中的安装，其左右两端各加工一个直径5mm、高10mm的圆台。

图1 几何模型(单位：mm)

1.2 数学模型

1.2.1 控制方程欧拉-欧拉模型通过压力和相间交换系数实现动量方程与连续性方程的耦合。对于泥浆流动、水力输送和颗粒沉积的计算，使用欧拉-欧拉模型较为合适[22]。计算中进行了以下简化假设：①流动和传热都是稳态过程；②忽略换热管向环境的热量损失；③不考虑质量力和流体的黏性耗散。因此，欧拉-欧拉模型的控制方程如下式(1)～式(3)。

连续性方程

本文采用RNGk-ε模型进行湍流求解，其加入了旋流对湍流的影响，提高了预测旋流和快速应变流的精度[23]，如式(4)、式(5)。

湍动能k方程

1.2.2 颗粒污垢模型

颗粒污垢特性是沉积和剥蚀迭加后的结果，学者们采用不同的数值模型研究了颗粒沉积过程，如拉格朗日粒子跟踪法、欧拉两相法、离散粒子模型（DPM）和Mixture模型等[24-27]。文献[28]详细介绍了本文采用的颗粒污垢模型，考虑了水中微米颗粒沉积的四种有效沉积机制：扩散沉积、湍流泳沉积、热泳沉积和重力沉积。颗粒沉积率计算如式(6)和式(7)。

颗粒剥蚀受壁面剪切力和污垢黏结强度因子影响，剥蚀率计算如式(8)。

假定污垢均匀分布在换热管内壁面上，且污垢成分一致，则污垢热阻如式(9)。

1.3 边界条件

重力加速度为y轴负方向，连续相为水，分散相为氧化镁颗粒，氧化镁颗粒的直径是20～50μm，颗粒浓度为0.4～4.0kg/m3。转子和管壁的材质分别为铝和铁，水的黏度随温度变化较大（见表1），其余物性参数均为常物性（见表2）。采用速度入口和压力出口边界，入口速度为0.15～0.35m/s，介质入口温度为293～323K，管壁温度恒定为343K。采用多重参考系法模拟转子的转动，将流体域划分为转子所在的旋转域及其他静止域。

表1 不同温度下水的黏度[29]

表2 物性参数

1.4 全局敏感性分析方法

假设输入与输出之间的关系可以描述为Y=f(x)，其中Y=(y1,…,ym)是标量输出函数，x=(x1, … ,xn)是n维自由输入变量。根据Sobol'[30]提出的全局敏感性指标，一阶效应敏感性，也称为主影响，用来衡量xi对y的单独贡献，如式(10)。

式中，V(Y)为f(x)的无条件方差；V[E(Y|xi)]为fi(xi)的方差。

由于一阶敏感性仅表现为xi的变化，只有各个输入变量的解耦效应，而高阶敏感性指标通常不计算。Homma 和Saltelli[31]提出将x=(x1,…,xn)分成两部分，即xi和其互补子集x～i，互补子集包括除xi以外的所有因素，总影响如式(11)。

STi=Si+Si,～i= 1 -S～i(11)

式中，Si,～i为xi与其他所有元素的相互影响；S～i为x～i的主影响且不包括xi对应的任何影响。

最佳预后元模型的建立是一个自动选择过程，元模型库包括多项式回归、最小二乘法和克里金模型。此外，考虑Box-Cox变换以提高元模型的近似质量，该幂变换方法能够稳定方差，增强正态分布，减少非线性效应。变量筛选一方面是采用由皮尔逊积矩相关系数决定的显著性滤波器，只有当输入变量之间的相关误差在95%～99%以内时才考虑其对输出变量的显著性。另一方面采用重要性滤波器，其定义为全模型与约简模型的调整决定系数（adj.R2）之差。CoImin的最小阈值在1%～5%之间，如果CoI不小于CoImin，则将考虑该变量。在交叉验证的基础上，提出预后系数（coefficient of prognosis，CoP）衡量MOP 的近似质量。结合基于方差的敏感性指标，将CoP 与MOP 确定的总影响敏感性指标相乘，计算出单个输入变量的方差贡献，如式(12)～式(15)。

图2 展示了基于MOP 的全局敏感性分析流程。输入变量为颗粒直径d（20～50μm）、颗粒浓度c（0.4～4.0kg/m3）、入口速度v（0.15～0.35m/s）、介质入口温度T（293～323K），沉积率和污垢热阻作为输出参数，采用增强的旋转中心差分方法生成49 个实验点。

图2 基于MOP的全局敏感性分析

2 网格无关性验证与模型验证

当颗粒直径为20μm、颗粒浓度为0.4kg/m3、入口温度为303K、入口流速为0.3m/s 时，对传热管模型进行网格独立性验证。流体域采用四面体网格，并局部加密管壁边界层和转子壁面网格。通过改变整体网格、管壁第一层网格和转子壁面网格尺寸来调整网格数量，结果如图3所示。随着网格数量的增加，传热管内的沉积率和压降都是先减小后趋于稳定。当网格数量为1727959 时，沉积率为2.62×10-4kg/(m2·s)，相较于网格数为1250704 时变化了2.2%，同时压降变化了1.4%。综合考虑计算精度和计算时间，本文采用的网格数为1250704。

图3 网格独立性测试

通过复现文献[32]的实验结果，对MRF模型和RNGk-ε模型计算内置转子传热管内流场和温度场的准确性进行验证。以内管中的流体为研究对象，入口温度为313K，管壁温度为295K，传热管内径为20mm，扭带宽度和长度分别为14mm、2000mm，节距比为6.0，建立了与文献中一致的物理模型，模拟结果与实验数据对比如图4所示，模拟值与实验值的变化趋势一致，努塞尔数的偏差最小为6.1%，平均偏差为8.4%，表明内置转子传热管内流场和温度场计算的精度较高。

图4 内置转子传热管Nu数计算验证

通过复现文献[33]中不同入口流速下圆管内颗粒污垢热阻测试实验，对本文所采用的颗粒污垢数学模型进行验证，计算结果与实验结果对比如图5 所示。可以发现，随着时间的增加污垢热阻逐渐增加，并在一段时间后达到渐近值。模拟值与实验值增长速度相近，三种工况下渐近值的相对偏差均在12.0%以内，最小偏差为6.2%，符合两相流计算要求。渐近值最能反映换热管长期运行时所产生的污垢热阻，因此本文后面提及的污垢热阻均是渐近值。

图5 污垢热阻模拟值与实验数据对比曲线

图6 四个沉积速率的总影响CoP矩阵

3 结果与讨论

3.1 沉积速率全局敏感性分析

通过构建MOP 实现最重要参数和最佳元模型的自动检测，并验证了预后质量。由于对沉积速率进行量纲为1化时消除了颗粒浓度的影响，所以本节主要研究入口温度、入口流速和颗粒直径三个因素对各量纲为1沉积速率的影响。

3.2 沉积率与污垢热阻全局敏感性分析

图7 示出了四个输入变量对沉积率和污垢热阻影响的CoP 矩阵。从图中可以看出，md和Rf的全模型CoP 分别为97.8%、96.2%，预后模型都是基于BoX-CoX 变换的一阶多项式回归。单个CoP的最大值均为d所对应数值，说明d是影响颗粒沉积特性的最显著因素。对于md，d和v的总影响较为明显，c和T的总影响接近且较小。对于Rf，除d外其余三个变量的总影响相近，且远小于d的总影响。

图7 沉积率和污垢热阻的总影响CoP矩阵

图8 和图9 示出了四个变量对md、Rf的主影响和交互影响柱状图。图8 可以看出，四个变量对md的主影响和总影响的排列顺序一致，d对md的主影响比交互影响高出了71.6%，而v、c和T对md的主影响和交互影响十分接近。主影响和交互影响的总和分别为73.0%、59.7%，说明各变量之间的交互作用对md有着重要的影响。图9 可以发现，d对Rf的主影响明显高于交互影响，主影响为47.3%，而交互影响为12.8%，v和c对Rf的主影响均略高于交互影响，T对Rf的主影响比交互影响低48.1%。主影响和交互影响的总和分别为75.7%、42.5%，各变量之间的交互作用依然不可忽视。

图8 md的主影响和交互影响柱状图

图9 Rf的主影响和交互影响柱状图

图10和图11分别示出了混合项对md和Rf的预测响应，每个云图均在其他参数取变化范围中间值时得到。由图10 可知，最大的沉积率都出现在云图的右上角，说明随着四个输入变量的增加，沉积率增大。随着c的提升，更多的颗粒与壁面发生碰撞，单位时间内沉积到壁面上的颗粒增加。热泳运动中颗粒由温度较高向温度较低的地方移动，随着T的增加，壁面与流体的温差减小，热泳运动对颗粒沉积的抑制作用减弱；同时T提高使流体的黏度减小，湍流泳沉积迅速增加。d越大，重力的作用越明显，重力沉积快速增加。v越大，单位时间内流体携带更多的颗粒流向壁面，沉积率增加。

图10 预测混合项对md的响应

图10(a)示出了c和T在v=0.25m/s、d=35μm 时对md的影响，可以看出，在参数变化范围内沉积率一直偏低，且沉积率等值线与竖直方向呈45°夹角，说明c和T对md的综合影响较小，且两者对md的影响效果相近。在图10(b)、(c)、(f)中，当d=20μm 时，c、T、v对md的影响很小；随着d的增加，c、T、v对md的影响逐渐增加，说明c、T、v对md的局部敏感度逐渐增加。同理，在图10(a)、(d)、(e)中，等值线从图中左下角到右上角越来越密集，说明随着各变量的增加，输入变量对md的局部敏感度越来越大。

图11 表明污垢热阻的最小值都出现在每个云图的左下角，随着c、T、d、v的增加污垢热阻逐渐增大，与沉积率的变化规律一致，且图11(a)与图10(a)类似，表明c和T对Rf的影响效果相近。在图11(b)、(c)、(f)中出现了较高污垢热阻（红色区域），说明d和c、d和T、d和v对Rf的综合影响十分显著。保持c=2.2kg/m3、T=308K 不变，当d=20～26μm，Rf几乎不随v变化。此外，等值线从左下角到右上角逐渐变得密集，表示输入变量对Rf的影响越来越明显，表明随输入变量的增加，各输入变量对Rf的局部敏感度越来越高。