第十五讲基于UD+CFD耦合优化的复混肥造粒塔进风口结构设计探讨

徐静安 吴文华 陈明良

技术讲坛

第十五讲基于UD+CFD耦合优化的复混肥造粒塔进风口结构设计探讨

徐静安 吴文华 陈明良

0 引言

熔体高塔造粒过程主要是以熔融尿素为载体，与配比好的其他物料在热熔状态下搅拌混合，形成具有一定流动性的低温共熔料浆，再通过旋转差动喷头将料浆以液流的形式喷出。液流在塔内自然下落过程中依靠表面张力成型，并与上升的空气热交换，冷却凝固成球形颗粒。复混肥在造粒塔内的冷却效果直接影响整个造粒过程，塔底复混肥黏壁、颗粒不圆整等均与塔内复混肥冷却效果欠佳有关。而影响塔内复混肥冷却过程的因素主要有进风口气温、复混肥颗粒粒径、进出风口结构、塔高等。造粒塔的结构不仅影响着复混肥冷却过程，也与造粒塔的造价成本息息相关。因此，在造粒塔设计阶段，有必要对造粒塔结构（如进风口结构等）进行研究分析，尝试找到合适的造粒塔结构。

迄今为止，国内已建成上百套造粒塔，但关于造粒塔结构对复混肥冷却过程影响的研究甚少，尤其是对造粒塔透风窗结构方面的研究。国内不少学者对造粒塔内颗粒冷却过程进行理论分析，将造粒塔内的空气流动简化成活塞流，在计算塔内风量时往往忽视了造粒塔进风口结构对塔内风量的影响，这势必将降低计算结果的准确度。上海化工研究院长年致力于复混肥高塔造粒工艺研究工作，近年来采用计算流体动力学（CFD）软件与现场测试数据相结合的方法建立了一套数学模型，该数学模型根据造粒塔的实际几何结构进行建模，模型计算预测出的结果与现场测试数据的误差在可接受范围内。鉴于造粒塔设备庞大，塔高90～120 m，通过实验及实测来掌握其性能工作量很大，本文在己建立的流体力学数学模型基础上，通过耦合优化方法改变几何模型，初步探索研究造粒塔进风口结构对复混肥熔融造粒塔内风量的影响，以期为造粒塔几何结构设计提供技术参考。

耦合优化的逻辑框图见本刊2017年第五期出版的第十三讲“对基于‘试验设计+模拟计算’耦合优化研究、设计方法的推广”。

1 初始模型计算方法

1.1 造粒塔初始结构

造粒塔内区域如图1所示。造粒塔由混凝土筒体、进风口、出风口、旋转差动喷头和锥形接料斗等组成。造粒塔初始结构参数如表1所示。

图1 造粒塔结构示意图

1.2 网格划分

造粒塔的模拟区域比较大，采用分区结构画网格。将喷头以上区域定义为气流区，该区域网格为非结构化网格(Tet/Hyblrid)。将喷头以下区域定义为颗粒区，该区域网格为结构化网格（Hex/Wedge）。另外，根据网格尺寸及数量对造粒塔进行网格独立性检验，根据计算结果和计算时间选择恰当的网格划分方法。本文根据造粒塔尺寸共进行了三种水平的网格划分，并对不同网格水平的的造粒塔模型分别进行模拟计算。鉴于网格数较大,常规计算机无法胜任，本次模拟计算工作均在上海化工研究院IBM高性能计算机平台上完成。具体网格数、模拟计算时间及模拟结果见表2。模拟结果表明，Level 2的网格数（即7543158）具有一定的合理性，在此基础上增大网格数，模拟误差减小不明显，而计算时间增长，且残差值增大甚至可能无法收敛。考虑时间消耗，对造粒塔结构探索工作中选择网格密度水平为Level 2。

表1 造粒塔初始结构参数

表2 网格独立性检验结果

1.3 颗粒追踪方式

本文数值模拟计算中，连续相均采用稳态的计算方式，颗粒追踪方式分别采用稳态追踪和非稳态追踪。颗粒稳态追踪时，“Number of Continuous PhaseIterations Per DPM Iteration”设置为20，即分散相每迭代一步连续相迭代20步。

颗粒非稳态追踪时，时间步长设置为0.1 s，时间步数为1，即每隔0.1 s颗粒喷射一次，每0.1 s对每个颗粒进行一轮包含一步的轨迹计算，每一轮颗粒轨迹计算中连续相迭代20步。

颗粒非稳态追踪与颗粒稳态追踪的模拟计算结果对比如表3所示。

表3 颗粒非稳态追踪模拟结果对比

颗粒非稳态追踪与稳态追踪的数值模拟结果非常接近，模拟结果与实测结果的相对误差均在5%以内，其中颗粒非稳态追踪的模拟结果误差相对较小。理论上，颗粒稳态追踪仅适合于颗粒稀少且颗粒对连续相影响较小的气粒多相流的数值模拟，而颗粒非稳态追踪在模拟计算颗粒较密集、颗粒对连续相有一定影响的气粒多相流领域更胜一筹。由于造粒高塔体积庞大,颗粒下降速度快、在塔内停留时间短，塔内气粒多相流中颗粒所占体积分数较小，颗粒稳态追踪方式也适合于用来模拟计算塔内颗粒运动轨迹，最终体现在颗粒非稳态追踪与稳态追踪两种方式模拟计算塔内气粒多相流的模拟结果非常接近。考虑到非稳态计算耗时长且对计算机的性能要求高，本文在对造粒塔结构的模拟探索计算中均采用颗粒稳态追踪的方式。

1.4 物性参数及边界条件

1.4.1 气相定解条件

对于气体入口与出口边界条件，只需给定入口压强p、温度T、湍流强度I以及水力直径Di。以其中一个模拟计算为例，气相进口条件及参数值为：p= 101 325 Pa，速度方向为垂直于进风口截面，T=296.9 K，I=0.04，Di=3.4 m。气相出口条件与参数值为：p=101 325 Pa，速度方向为垂直于出风口截面，气流温度T=296.9 K，I=0.05，Di=6.8 m。

1.4.2 颗粒定解条件

颗粒边界条件主要包括颗粒进口条件和颗粒与壁面碰撞条件。

（1）颗粒进口条件

颗粒进口流量为10.56 kg/s，温度为386.85 K (113.7℃)，假设固体颗粒除了比热容外其他各种性质与料浆相同，进口物料熔液物性参数见表4。

颗粒从锥台射出，速度大小为3.31 m/s，切向分量百分数0.999 989 1。锥台注射的颗粒粒径分布参考塔底颗粒粒径分布（如表5所示）。

表4 复混肥物性参数

表5 塔底颗粒粒径分布

模拟计算过程中将颗粒粒径分布近似成R-R分布，遵从R-R分布的颗粒群中，颗粒直径d与尺寸小于d的颗粒所占的质量分数Rd有如下关系：

式中，d为特征尺寸，其含义为直径小于d的颗粒质量分数Rd=1-e-1=0.632 1；n为尺寸分布指数，该参数值越大，颗粒尺寸分布越均匀。数值模拟过程中，颗粒群的特征尺寸d为2.48 mm，分布指数n为3.5。

复混肥颗粒在塔内的降落过程还伴随着凝固放热过程。颗粒从外往里开始凝固，因此存在开始凝固温度Tf和完全凝固时的平均温度Ts。近似地认为颗粒温度高于Tf时为液态，温度低于Ts时为固态，凝固潜热采用修改比热容的方式加入。具体修改方式如下：

式中，C1为液态料浆热容，CS为固态颗粒热容，L为凝固潜热。

复混肥配方为20-5-20（硫基），该种复混肥的开始凝固温度Tf=86.1℃，完全凝固时的平均温度Ts=83.7℃，凝固潜热L=87.34 kJ/kg。

（2）颗粒出口条件

简化模型中将颗粒出口面设置为“壁面”，颗粒与出口面的碰撞条件取“逃逸”，也就是颗粒到达出口面时计算终止。

（3）颗粒碰撞条件

液滴进入造粒塔内，与空气逆流接触实现热量交换，由于现有工艺条件下料浆液滴较少黏壁，固体肥料碰壁后假设成弹性碰撞，因此颗粒与壁面的碰撞条件取“弹性碰撞”。

1.5 计算原理

采用分散相模型进行计算，控制方程包括连续相方程、分散相方程和气粒两相间热量传递方程。连续相求解过程为稳态迭代计算，颗粒追踪方式采用稳态追踪。所谓稳态方式是指每隔若干个连续相流场迭代步，在当前流场状态下，逐个对每个颗粒进行从起始位置直到其终了的轨迹积分计算及源项计算。稳态方式得到某一时刻连续相流场条件下，在一系列积分时间步的颗粒状态，一系列颗粒位置可连成运动轨迹线。

热压自然通风主要是依靠浮升力作用的一种通风方式，自然通风的数值模拟研究方法一般为扩大计算域，在建筑的外围建立一较大的模型形成无限远边界，在外围入口处设置为速度入口边界条件来考虑侧风对其的影响。不考虑侧风的影响时，可以将速度边界的速度值设得很小，或者改为压力边界条件，或者直接将建筑的进出口处均设为压力入口边界条件。综合考虑后，决定将造粒塔进出口边界条件均设为压力边界条件，内部空气密度也采用波希涅斯克(Boussinesq)假设。

离散方式采用有限差分法，连续相方程采用SIMPLE算法迭代求解，分散相建立在连续相初始求解值上，并与连续相进行耦合求解。

2 造粒塔初始条件模拟结果分析

本文针对两座工程塔的三次不同工艺参数下的造粒过程分别进行模拟计算，具体模拟结果及工程标定误差如表6所示。其中“模拟1”与“模拟2”针对的造粒塔进风口结构：n=9，s=74.46 m2、θ=139.81°；“模拟3”针对的造粒塔进风口结构：n=14，s=90.72 m2，θ=160.51°。表6显示，模拟结果与工程化高塔实测结果的相对误差均在10%以内，说明建立的数值模型可靠，得出的数值模拟结果具有一定的可信度。

表6 数值模拟结果对比

具体针对“模拟1”文件进行造粒塔内气速气温云图分析，如图2所示。并对造粒塔内颗粒分布、运动降温规律等进行分析，如图5与图6所示。造粒塔外面的冷空气在压差的作用下，从9个进风口进入造粒塔，9个不同气速的气流在塔中心汇流形成涡流。从图2中可以明显看出塔中心的涡流，使得粒子发生错流，轨迹发生变化，部分粒子可能从气速小的窗口逃逸，现场测试中也发现通风窗口落有不少粒子。气体汇流后沿着塔轴向上升，上升过程中与下降的复混肥颗粒进行热量交换与动量交换。

图2 塔进风口处横截面上气速矢量图

图3和图4显示的是造粒塔内xz截面温度云图、速度云图、温度等值线和速度等值线。从图3中可以看出温度等值曲线呈现规则的W形，出现两处低谷一处高峰。这与造粒塔颗粒分布规律有关，如图5所示，颗粒主要集中在塔半径4～7 m区域内，颗粒分布规律是由喷头喷出颗粒的粒径分布规律导致的，粒子的粒径主要集中在2～3.36 mm。塔半径4～7 m区域内粒子比较多，导致该区域空气温度上升得比较快，升高到相同温度所需上升的高度比起他区域要低，最终导致轴截面的温度等值线中出现两处低谷。塔中心粒子很少，热交换不明显，气温升高到相同温度所需升高的高度比其他区域要高，最终导致轴截面的温度等值线中出现高峰。气体上升过程中，温度也在上升，在喷头区域温度梯度明显变大。这是因为喷头区域下降的是高温熔融复混肥并且颗粒比较密集，复混肥与空气的温差比较大，传热推动力大。喷头区域的热空气继续向上运动，由于渐渐远离热源以及与塔壁附近的冷空气换热，气温逐渐下降。

图3 塔xz轴截面上气温云图

图5 距进风口35 m处横截面上颗粒分布图

从图4中可以看出进风口气速最大，空气进入造粒塔后由于截面迅速增大导致速度迅速减小，空气上升过程中沿轴向速度变化不是很大。进入喷头区域由于热压增大，因此气体明显提速，进入热气通道后横截面减小，气速再次提高。塔半径4～7 m处气速明显较小，这与粒子集中在该区域，空气受到较大的阻力有关。

图6 塔内颗粒降温规律图

颗粒在造粒塔内下降并与上升空气发生对流换热，颗粒温度逐渐降低，依次经历液态冷却过程、凝固过程、固态冷却过程。从图6可以看出2.0 mm粒径颗粒下降5 m左右进入凝固过程，下降35 m左右进入固态冷却过程；2.5 mm粒径颗粒下降10 m左右进入凝固过程，下降50 m左右进入固态冷却过程；而3.9 mm与4.1 mm粒径颗粒下降25 m左右进入凝固状态，落至塔底时依旧处于凝固状态，颗粒温度最终为86.2℃。不同粒径颗粒的传热系数与停留时间的差异最终导致了颗粒降温曲线之间的差异。大粒径颗粒对流传热系数与停留时间均小于小粒径颗粒的对流传热系数与停留时间，降落至塔底时凝固潜热尚未完全释放出来。

3 进风口结构的耦合优化探索

3.1 均匀试验设计

本文主要是在已有的被校核的模拟计算基础上，通过改变进风口结构分析研究进风口结构对自然通风造粒塔内风量V的影响，其中造粒塔结构参数及工艺参数如表7所示。考虑到影响塔内风量的因素较多，只选取对塔内风量影响较大的3个因素。即通风窗个数n、进风口总面积s以及通风窗对应的总圆心角θ。各因素选取水平参照均匀设计表U7*(74)，各因素具体取值如表8所示。

表7 结构参数及工艺参数

表8 因素水平表

其中造粒塔具体进风口结构如图7所示，造粒塔底部径向分布着若干个矩形通风窗，通风窗个数n一般为6～18个，进风口总面积s一般为66～96m2，通风窗对应的总圆心角θ一般为100～180°。

图7 某造粒塔进风口结构布置图

3.2 均匀设计结果分析

对于高约l00 m的工程造粒塔，现场实测标定一次工作量相当大，涉及到结构尺寸的变动改造几乎不可能。而基于耦合优化方法的均匀设计方案，开展对应的数值模拟计算，以造粒塔内风量V为分析指标可以得到相应的技术信息，具体结果见表9。

表9 模拟试验结果

通过数据处理系统（DPS）软件对模拟试验结果进行二次多项式逐步回归统计建模，计算结果如下所示：

计算结果表明各项检验显著，模型拟合效果在可接受范围之内。t检验显示在实验范围内对V的影响次序：s×θ＞n×θ。数学统计模型显示造粒塔内风量随着通风窗个数增大而减小，随着通风窗总面积与总圆心角的增大而增大。此外,数学模型预测的最优化进风口结构为：n=6，s=96 m2，θ=180°。

4 耦合优化结果验证

统计模型耦合优化预测的进风口结构为n=6，s=96 m2，θ=180°，对造粒塔重新进行几何建模，并结合对应的工艺参数进行模拟计算，模拟结果显示优化结构下造粒塔风量为152.35 m3/s，模拟计算值与统计数学模型预测值的相对误差为1.05%。直接根据优化进风口结构建造造粒塔，造价上千万。鉴于直接工程验证优化结果成本太大，本文选择对已有的山东某造粒塔进行实测标定，验证耦合优化模型预报的质量。山东造粒塔进风口结构为：n=9，s=74.46 m2，θ=139.81°，生产工艺参数如表7中所示。根据数学模型V=141.649 736 6-0.003 530 702 055×n×θ+ 0.000 932 333 124 6×s×θ计算，计算结果如表10所示，结果显示模型预测值与实测标定值的相对误差为0.975 4%，误差在可接受范围内，说明回归的耦合优化模型具有显著的统计意义。

5 结论

（1）介绍了复混肥造粒塔CFD数值模拟计算原理与过程，并分析研究了造粒塔内空气温度、速度云图以及复混肥颗粒在塔内的运动规律。分析结果显示造粒塔内轴截面气温等值线呈现W形状，气速等值线呈现M形状，在进风口处造粒塔中心与通风窗的上端均出现明显的涡流。颗粒在造粒塔内降落过程中主要集中在塔半径4～7 m的范围内，依次经历液态冷却过程、凝固过程、固态冷却过程，颗粒温度随路程的变化曲线出现两个明显的转拆点。

表10 模型预报验证结果

（2）基于均匀设计法开展数值模拟的耦合优化研究，探索造粒塔进风口结构对塔内风量的影响。通过对模拟结果进行二次多项式回归分析，得出数学统计模型，并最终确定了试验范围内的最优化进风口结构：n=6，s=96 m2，θ=180°。

（3）结合工程实测标定，对耦合优化模型预报质量进行验证，结果显示模型预测值与实测标定值的相对误差在可接受范围内，回归的模型具有显著的统计意义，可作为工程设计技术基础。

后记：

最近，阅读了《中国人才》2015年9月刊登的一篇文章——各国加强科技人才队伍建设特点与趋势（上）。各国均在创新人才培养方面强化科学、技术、工程和数学(STEM)教育，借此推动国家创新和经济的发展。英国政府在2011年—2015年实施支持STEM教育计划，美国2013年公布了联邦STEM教育五年战略计划，加拿大政府2014年12月启动STEM的人才培养计划，韩国在2011年—2015年强化科学、技术、工程、艺术、数学(STEAM)人才的系统培养。

多年来，上海化工研究院支持和推动数学工具、数字化技术在应用基础研究、技术开发、工程化转化工作中的学习、探索、应用。2010级硕士吴文华在读期间及留院工作后，结合高塔造粒项目，不断探索UD+CFD耦合优化技术，并在院部配置的IBM高性能计算机上重新做了计算，不断完善。

“沉舟侧畔千帆过，病树前头万木春”，青年同仁们，努力啊！

徐静安男教授原上海化工研究院院长长期从事化工机械、化学工程、化工工艺开发、化工试验设计与数据处理方面的研究