基于响应面代理的SCR后处理系统快速预测模型及性能优化设计

2022-02-16刘严殷雷吴涛涛张锴孟境辉

中南大学学报（自然科学版） 2022年12期

刘严，殷雷，吴涛涛，张锴，孟境辉

(1.浙江大学控制科学与工程学院，浙江杭州，310027；2.上海镭融科技有限公司，上海，201800；3.天纳克(苏州)排放系统有限公司，江苏苏州，215300；4.华北电力大学热电生产过程污染物监测与控制北京市重点实验室，北京，102206)

柴油发动机具有高燃烧经济性和高可靠性，被广泛应用于车辆航海等领域，但由此造成的环境污染问题日益突出，其尾气排放物已成为空气污染的重要来源。柴油发动机排放的主要污染物为氮氧化合物(NOx)与颗粒物(PM)[1]。目前，控制柴油发动机排放尾气中污染物含量通常采用选择性催化还原(selectivecatalyticreduction，SCR)技术[2-4]。SCR后处理系统主要性能指标有排气压降、速度均匀性和氨均匀性[5-7]。国内学者通常基于SCR 后处理系统数值模型提高系统性能，帅石金等[8]对柴油机中的SCR催化器进行了优化设计，并进行了后处理过程数值计算，分析了不同方案下载体前端面的尿素浓度分布及NOx转化率；朱博文[9]设计了一种优化混合器，提高了速度均匀性，氨气均匀性、温度分布和尿素转化效率等指标；HUANG等[10]提出了一种混合器的优化方案，发现SCR 入口处氨均匀性增加，混合器压降降低，从而改善了后处理系统的性能。

通常，提升SCR 后处理系统性能的方法是根据设计人员的经验主观调整系统结构设计参数，并通过仿真计算或实验进行验证。朱博文[9]针对现有后处理系统均匀性指标较低的问题，结合后处理系统传统拓扑设计经验，提出新型优化混合器结构，对比优化前，系统性能确有提升，然而事实上，当后处理系统的几何拓扑结构发生改变时，原有几何结构下的流动传热数值模型不再适用。这种方法不仅耗时，也未考虑到SCR 后处理统几何结构与数值模型中多物理场之间的强耦合问题。因此，一般情况下很难得到SCR 后处理系统合理的最优设计。

针对传统优化方法难于实现最优化的问题，智能优化算法作为一种启发式算法，在工业优化领域有很多成功的研究。其中，基于遗传算法改进的非支配排序遗传算法[11-12](NSGA-II)收敛快，计算复杂度低，适用于复杂多目标非线性优化问题求解。WANG等[13]在考虑温度均匀性和压降的前提下，采用代理模型和多岛遗传算法对换热器的结构参数进行了优化；HAUSWIRTH等[14]提出了在可行集上的连续时间投影梯度下降算法，并将其用于电力系统在线潮流优化；MENG等[15]采用基于非支配排序遗传算法的多目标优化方法，实现对TEGs/TECs 和TEG-TEC 集成系统的自动优化；WANG等[16]利用非支配排序算法对微通道散热器进行了优化，找到了最佳热阻和抽运功率。

然而，数值模型求解复杂，将其作为正问题求解器与智能优化算法结合的优化策略计算成本过高。因此，国内外学者对基于响应面分析的快速建模方法进行了大量研究。XIA等[17]提出了一种多物理场仿真和响应面方法结合的优化方法，减少了仿真试验次数；NAYAK等[18]采用基于响应面法的统计模型，评估了控制因素对复合材料腐蚀速率的影响；REN等[19]基于整数遗传算法与响应面方法，提出一种用于可再生制冷和供暖系统的优化设计，显著降低了计算成本。

综上，本文基于改进的非支配排序遗传-TOPSIS决策算法(NSGA-II-TOPSIS)对SCR后处理系统结构进行优化设计。首先，利用有限元原理，建立了SCR后处理系统有限元数值预测模型，基于部分数值模型计算结果建立系统响应面数学模型，确定影响系统性能的关键结构参数；其次，通过系统的数学模型，将优化问题转化为一个在可行域内的多目标寻优问题，以系统压降最小、氨均匀性和速度均匀性最大为多目标，以关键结构参数为优化变量，采用改进的非支配排序遗传算法进行优化；最后，结合TOPSIS决策算法在互相矛盾的性能目标函数中取得折中优化设计方案，并在STAR CCM+中验证了最优设计的有效性。

1 有限元数值模型

1.1 几何模型

图1所示为SCR后处理系统结构图。由图1可见：柴油机车尾气SCR后处理系统通常为直线型，主要由入口区域、混合器区域、氧化型催化器(DOC)区域、壁流式颗粒捕集器(DPF)区域、选择性还原催化器(SCR)区域和出口区域6个部分组成。尾气由进口区域流入，首先通过DOC 区域，将尾气中的CO、碳氢化合物转化为无害的CO2和H2O；经过初步处理的尾气再流入DPF 区域，进行柴油颗粒的捕捉；经过上述处理后的尾气流入到混合器中，并与从尿素喷嘴中喷射的尿素颗粒充分混合后进入到SCR 区域进行选择性催化还原，将混合气体中的NOx还原成N2和H2O；最终，满足排放要求的尾气自出口区域排出。后处理系统的关键结构是混合器，主要由进口多孔挡板、尾气尿素混合反应管道以及出口多孔挡板3 个部分构成。表1所示为系统多孔介质区域结构参数；表2所示为系统混合器结构参数。

图1 SCR后处理系统结构图Fig.1 Structure of SCR post-processing system

表1 系统多孔介质区域结构参数Table 1 Structure parameter of system porous media area

表2 系统混合器结构参数Table 2 System mixer structure parameter

1.2 控制方程

本文建立的柴油机SCR 后处理系统模型满足一般流体的流动、传热物理规律，控制方程主要为质量守恒方程、能量守恒方程与动量守恒方程。陈莲芳等[20]不考虑结晶风险，对SCR 反应器入口烟道和导流部件进行模拟优化设计，提高了催化剂入口的速度均匀性。需要指出的是，本文优化主要针对系统典型性能指标，尚不考虑结晶风险对系统长期寿命的影响。

模型满足以下假设条件：

1) 流动具有不可压缩的稳态特性；

2) 忽略重力对流体流动的影响；

3) 流体和固体的性能参数不随其他物理量的变化而变化；

4) 不考虑壁面液滴附着和壁面液膜的形成；

5) 后处理系统外表面绝热；

6) 进口尾气采用理想干空气。

1.2.1 连续性方程

单位时间内流进、流出系统的体积流量相等。

式中：ρ为尾气密度，kg/m3；u为尾气速度矢量，m/s；t为时间，s。

1.2.2 动量方程

尾气在排气管中的流动模型为湍流，采用可实现的k-ε两层湍流模型，该模型将可实现k-ε模型与两层方法[21]相结合，主要用以求解湍动能k和耗散率ε。

式中：为流体的平均速度，m/s；μt为流体的动力黏度，Pa·s；σ为普朗特常数，本模型中σε=1.2，σk=1.0，Cε1=1.44，Cε2=1.44，Sε，Sk为均为系统定义的常数；ε0为抵消湍流衰减的环境湍流值；f2为大涡流的时间尺度；T0和Te分别为时间尺度。

式中：fc为曲率修正系数；Gk为流体流动的平均速度梯度所引起的湍动能；Gb为流体浮力引起的湍动能；γM为膨胀对总的耗散率的影响项；Gε3为模型系数。

1.2.3 能量方程

单位体积能量的增加等于流进单位体积的热量和外力对单位体积所做的功之和，其表达式如下：

式中：˙为热通量向量，W/m2；T为黏性应力张量，N；v为流体速度，m/s；fb为单位体积流体受体积合力；Su为外界提供能量，kJ；E为总能量，kJ；H为总焓，kJ/kg；h为静态焓，kJ/kg。

1.2.4 多孔介质模型

本文建立的模型通过黏性损失项与惯性损失项计算流体的流速，该部分是在动量方程基础上实现，其计算方程如下：

式中：Dij与Cij分别为黏性损失系数矩阵和惯性损失矩阵中的系数，由下式确定。

式中：Dp为多孔介质中的颗粒平均直径，m；χ为多孔介质孔隙率。

如图1所示，在本文建立的模型中，多孔区域有3个部分，其中DOC、DPF和SCR这3个多孔区域z方向黏性阻力系数分别设置为1 260.66，5 729.39 和5 729.39 kg/(m3·s)，惯性阻力系数分别设置为9.89，22.47 和30.01 kg/m4。需要指出的是，考虑到尾气沿z方向流动，以上3个区域x方向与y方向的黏性阻力系数和惯性阻力系数均设置为z方向的 1 000 倍，这样既可以保证流体在轴向方向有很好的流动性，又能保证结果的收敛性。

1.3 评价指标

SCR后处理系统的评价指标主要有氨均匀性、速度均匀性以及排气压降。

均匀性一般是以SCR 入口面作为统计面，均匀性直接反映催化反应的效果以及尾气的流动情况，氨均匀性Um与速度均匀性Uv的计算公式如下：

式中：ωm为NH3质量分数；ωi为局部NH3质量分数；vm为选定面上的速度；vi为局部速度；n为取值点的个数；Mi为单元格i的面积；M为选定面的总面积。

排气压降是评价后处理系统性能的另一重要指标，排气压降大表明在排气管路发生堵塞，易导致柴油发动机输出功率减小，增加发动机油耗；排气压降较小的SCR 后处理系统往往动力性与经济性均较好。

式中：n′为系统入口处横截面中的单元格总数；Pi,in为系统入口面各单元压力；m′为系统出口处横截面中的单元格总数；Pi,out为系统出口面各单元的压力。

1.4 边界条件与物性参数

尾气通道进口设置为质量流量入口，流量为Min=733 kg/h，入口温度Tin=608 K。尾气通道出口设置为压力出口，设置压力Pout=100 kPa，以实现流体在流道内充分自由流动，并抑制气体回流。尿素喷嘴入口流量M0设置为0.108 kg/h，尿素颗粒温度T0设置为300 K，其中尿素水溶液是由质量分数为32.5%的尿素和67.5%的水混合而成，喷射方式为贴壁射流，尿素颗粒碰壁模型采用Bai-Gosman 碰壁模型[22]，液滴破碎模型选用TAB 液滴破碎模型[23-25]，并使用泰勒模拟[26]。对于喷射器，喷射外锥角θ=9°，不考虑液膜。整体系统均采用混合钢材料，密度ρ=7 700 kg/m，比热容C=1 200 J/(kg·℃)，导热系数λ=0.12 W/(m·K)。

1.5 模型验证

为了验证本文建立的SCR 后处理系统数值模型，将其预测结果与文献[27]的计算结果进行对比。文献[27]采用计算流体力学方法分析了某国六b柴油机后处理系统排气压降与各段载体速度均匀性，本文采用与文献[27]相同的设置，文献[27]与本文数值模型的排气压降分别为23.03 kPa 与 23.87 kPa，相对误差为3.65%；速度均匀性分别为0.897 和0.884，相对误差为-1.46%。通过模型验证，验证了本文所建立的直线型SCR 后处理系统数值模型的合理性与正确性。

2 响应面数学模型

2.1 响应面分析法

在本文中，响应面分析法将后处理系统三大性能指标分别作为多个因素(优化变量)的函数(优化目标)，以小部分数值模型的仿真结果作为实验数据，设计响应面实验拟合出优化目标与优化变量之间的函数关系，将正问题求解器从计算复杂的数值模型转化为简单高效的显式数学代理模型。

响应面分析法的数学基础是泰勒展开，其核心思想就是用多项式来拟合任意的函数，对于多变量问题采用泰勒展开策略，当只考虑一次项与常数项时，拟合出的函数值往往与实际值偏差较大，一般采用二次项拟合就会出现较为理想的结果，即

式中：y为输出量；β为系数；x为输入变量；e为观测误差；k′为超参数的数目。

将2个变量的乘积合并用一个变量代替，写成如下形式，

多次观测后的矩阵形式为

拟合误差L计算公式为

当拟合误差最小时，拟合曲面和实际值最接近，对L的偏导数为0时，拟合误差达到最小值。

拟合的响应面即为

2.2 实验设计

将直线型SCR 后处理系统中混合器的关键结构参数即混合器管路半径(因素A)、混合器管路侧向开孔率(因素B)以及混合器后侧挡板开孔间距(因素C)确定为输入变量，将后处理系统压降、氨均匀性和速度均匀性确定为响应变量。为保证响应面数学预测模型的精度，采用设计空间内均匀撒点的自定义设计方法进行响应面分析。

输入变量取值见表3，为了兼顾数据的整体性，在每个变量的给定范围内均匀取值，共组成64组样本点。

表3 输入变量取值Table 3 Input variable value

2.3 因素显著性分析及数学模型

输入变量对拟合模型影响显著性程度主要由失拟项P来衡量，一般对于F值较大的拟合模型，P需小于10-5。因素显著性分析结果如表4 所示，由表4可见：输入变量均是显著性影响因素，进一步说明混合器管路半径、混合器管路侧向开孔率以及混合器后侧挡板开孔间距是影响SCR 后处理系统性能指标的关键结构参数。

表4 因素显著性分析Table 4 Factor significance analysis

模型拟合精度由拟合优度衡量，其为剔除回归方程中所包含项数的相关系数，该值越接近1，说明拟合效果越好。计算公式如下，

式中：S1和S2分别为残差的各项平方和与模型的各项平方和；f1和f2分别为残差的自由度与模型的自由度。通过响应面拟合，排气压降、速度均匀性与氨均匀性这3个响应面模型拟合优度均较高，分别为0.944，0.917与0.931，3个响应面数学模型均取得了较好的拟合效果。

本文最终获得响应面拟合的SCR 后处理系统数学模型，分别如下：

式中：Y1为SCR 后处理系统排气压降；Y2为速度均匀性；Y3为氨均匀性。

2.4 数学模型验证

为验证拟合的数学预测模型的精确性，以数值模型结果为实际值，以响应面数学模型输出结果为预测值，将二者进行对比。图2所示为排气压降、速度均匀性和氨均匀性实际值与预测值的比较结果。由图2可见：排气压降、速度均匀性和氨均匀性的预测值和实际值散点均匀分布在函数y=x两侧，因此，模型排气压降、速度均匀性、氨均匀性的预测值与实际值符合度较高，其中排气压降、速度均匀性和氨均匀性的预测最大误差分别为9.74%，2.05%和3.29%。

图2 响应面数学模型预测值与实际值对比Fig.2 Comparison between actual value and predicted value of response surface mathematical model

3 结构参数优化

3.1 初始结构数值仿真

对优化前SCR 后处理系统进行数值计算，得到如图3所示的后处理系统压力分布云图，SCR后处理系统各区域各分段的压降及排气总压降如表5所示。由图3 和表5 可见：混合器区域压降最大，主要因为该区域设置有结构复杂的多孔挡板结构，尾气流动阻力大、扰流强，压力损失为4.39 kPa，其余各区域压降较小，均在1～2 kPa，优化前系统总压降为11.22 kPa。

图3 SCR后处理系统压力分布云图Fig.3 Cloud map of pressure distribution of SCR after treatment system

表5 后处理系统各部分压降分布Table 5 Pressure drop distribution of post-processing system title kPa

图4 所示为初始解SCR 均匀性分布云图。由图4(a)可见：NH3主要分布在SCR 入口面中下区域，中部偏上区域NH3体积分数较低，氨均匀性为0.78，说明SCR 区域入口面的NH3分布不均，系统氨均匀性指标有待提高。由图4(b)可见：尾气与尿素颗粒在环形挡板处混合后流出后，SCR 入口面的中部出现流速较高的区域，速度分布不均，该表面速度均匀性为0.82，需进一步优化。

图4 初始解SCR均匀性分布云图Fig.4 Velocity distribution cloud map of initial solution SCR inlet surface

3.2 多目标函数与约束条件

传统的多目标优化是通过引入权重因子，将多目标优化问题转化为单目标寻优问题。例如，一个多目标的问题，它通常有一个复合函数f=ω1f1/f1,max+ω2f2/f2,max，权重因子ω1与ω2均是人为指定的，并且伴随主观经验而存在误差，因此这个复合函数的解往往不是最优。NSGA-II不需要扩展一个复合目标函数，而是直接以响应面分析法建立的数学模型为优化目标，同时能考虑所有的设计要求。

本研究中，优化变量设置为混合器管路侧向开孔率A、混合器后侧挡板开孔间距B和混合器管路半径C，上述变量均在各自取值范围内连续。SCR后处理系统主要性能指标包括排气压降P、速度均匀性Uv和氨均匀性Um，为兼顾以上性能指标并实现同时最优，本文定义的多目标函数如下，

3.3 NSGA-II-TOPSIS算法优化

NSGA-II-TOPSIS 优化流程图如图5 所示，具体如下：

图5 NSGA-II-TOPSIS优化流程图Fig.5 Optimization flow chart of NSGA-II-TOPSIS

1) 建立SCR后处理系统的多物理场数值模型，计算SCR后处理系统性能指标；

2) 通过响应面分析法，建立SCR 后处理系统数学模型，确定混合器主要结构参数(混合器管路半径、混合器管路侧向开孔率与混合器后侧挡板开孔间距)，并分析上述结构参数对后处理系统性能的影响程度，构建多目标优化命题；

3) 采用非支配排序遗传算法进行多目标优化，获得Pareto 最优解集后通过TOPSIS 评价方法按权重选择，获得决策解；

4) 将决策解代入数值模型中进行验证。

非支配排序多目标遗传算法种群规模设置为80，迭代次数设置为100 次，算法运行到第80 代停止迭代，三维Pareto 解集如图6 所示，表6 所示为最优解及对应优化变量取值。从表6 可见：第4代排气压降最小，第72代的氨气均匀性最高，第3代的速度均匀性最高；最终三者指标综合最优的TOPSIS 决策解(第34 代)。上述4 组解对应的优化变量也列于表6中。

表6 最优解及对应优化变量取值Table 6 Optimal solution and corresponding optimization variable value title

图6 三维Pareto前沿Fig.6 3D Pareto frontier

需强调的是，每组最优解均是在一次有效的优化运行中获得的，因此对于每1个个体(每组解)来说均是非支配的。每组最优解对应3 个目标指标，同时考虑到排气压降、速度均匀性和氨均匀性均是后处理系统重要的性能评价指标，故引入等权重评价因子，TOPSIS决策解为第34组解。

为了更直观地展现优化过程，随机选择进化过程中的个体与初始随机解进行对比。图7所示为优化后处理系统综合性能分析。从图7 可见：第3代解的速度均匀性指标在所有最优解中最高，为0.958 2，对比初始解中速度均匀性指标0.820 2，速度均匀性提升了14.40%；氨均匀性指标最高值为0.911 7，出现在第72 代解，相比较初始解氨均匀性，性能提升了14.31%；排气压降最优值7 kPa出现在第4 代解，与初始解相比，压降下降了37.64%。最终解为第34 代解，对比初始解，系统整体压降由11.22 kPa 下降至8.97 kPa，下降20.09%；速度均匀性由0.820 2 提升至0.915 3，速度均匀性提升10.39%；氨均匀性由0.781 2 提升至0.891 1，提升12.33%。

图7 优化后处理系统综合性能分析Fig.7 Comprehensive performance analysis of postprocessing system after optimization

由于每个最优解互相之间的关系是非支配的，上述3 组最优解虽分别在单个性能指标上达到最优，但另2个性能指标却不是较优的水平。通过分析性能指标之间的矛盾性，当SCR 后处理系统混合器结构复杂时，尿素与尾气在混合器区域混合充分，均匀性指标较高，然而复杂的混合器结构在一定程度会造成尾气流道的堵塞，引起排气压降升高。因此，系统排气压降与均匀性之间的优化关系是互相矛盾，即系统均匀性的提升往往是以升高系统排气压降为代价。为此，引入TOPSIS决策方法权衡3个指标，在所有最优解中获得了最终决策解(第34代解)，决策解的性能指标是综合最优的。

利用STAR CCM+软件对决策解进行仿真计算验证，验证结果如表7所示。排气压降的相对误差为2.99%，速度均匀性与氨均匀性指标的相对误差均小于2.00%。为了提高优化效率，基于响应面分析法建立的快速代理的数学模型进行反问题的优化，该数学模型本质上是原数值模型的降阶代理，两者并不完全相同。因此，对决策解进行仿真验证时会造成一定的仿真误差。在误差允许的范围内，决策解是正确有效的。

表7 决策解验证Table 7 Decision solution verification

决策解的系统排气背压、SCR 入口面处的氨气质量分数以及速度分布云图分别如图8和图9所示，各区域压降情况如表8所示。由图8可见：相比较于优化前，系统压降整体上均有所下降，对比各区域压降分布，发现混合器区域压降下降2.02 kPa，最显著，后处理系统未经优化的其他部分，压降也发生一定的变化，这是由于SCR 系统的后处理是一个流动传热多物理场复杂耦合的过程，并与系统的混合器几何结构关系紧密。因此，当混合器结构发生变化，一定程度上会影响整个SCR 后处理系统中的流体流场分布，进而导致系统未进行优化的其他部分压降发生变化。

图8 决策解SCR后处理系统压力分布云图Fig.8 Decision solution SCR post-processing system pressure distribution cloud map

图9 决策解SCR均匀性分布云图Fig.9 Decision solution SCR inlet surface velocity distribution cloud map

从图9可见：优化前，SCR入口面速度分布不均，入口面中部垂直区域流速大，四周流速较小；优化后，SCR 入口面速度分布不均情况明显得到改善。SCR 入口面NH3质量分数优化前区域性差异明显，中部偏下区域NH3质量分数普遍偏高，中部偏上区域NH3质量分数普遍较低，SCR 入口面氨均匀性分布不均；优化后，SCR 入口面NH3质量分数分布均匀，仅在中部偏下很小的范围内出现NH3质量分数偏高的情况。

柴油机的SCR 后处理系统的压降一般为10～30 kPa[28]，速度均匀性与氨气均匀性一般达到0.9左右[29]即能满足后处理系统对催化效率的要求。本文所得到决策解排气压降指标为8.70 kPa，远小于一般后处理系统排气压降的最小值，同时，系统均匀性指标也保持在较高的水平。考虑到理论方法与实际工业生产的区别，对于氨均匀性指标略低的情况在实际工业生产过程中是允许存在的。优化后的SCR后处理系统不仅排气压降显著降低，均匀性指标也同时保持在较高的水平，系统尾气处理性能得到大幅提升。