RESEARCH ON HETEROGENEOUS SOLID STRESS MODEL BASED ON ARTIFICIAL NEURAL NETWORK
-
摘要: 最小多尺度理论EMMS已经被引入多相质点网格法MP-PIC中, 建立了非均匀EMMS固相应力模型. 但现有的非均匀固相应力模型计算中, 中间步骤繁琐且花费时间长. 采用人工拟合的方式能获得非均匀固相应力表达式, 但需要人为确定拟合变量和拟合函数, 且针对于非均匀固相应力这种高度非线性函数所得到的拟合精度不高、用于MP-PIC模拟的结果相比原EMMS固相应力模型结果存在偏差. 针对上述问题, 本文提出通过机器学习的方法, 规避对固相体积分数的局部分布情况的表征, 并提出和建立能考虑颗粒浓度详细分布的人工神经网络ANN固相应力模型. 首先, 基于局部颗粒浓度和颗粒非均匀分布指数建立了双变量的ANN固相应力模型; 进一步将当前网格及其周边网格颗粒浓度组成的序列来详细表征颗粒浓度分布, 并建立相应的ANN固相应力模型. 然后, 将两种模型与EMMS固相应力模型进行了对比并测试了网格分辨率和粗化率对模型的影响. 研究表明: 基于ANN固相应力模型的模拟结果对EMMS固相应力模型结果有较高的还原度, 同时具有一定的网格分辨率无关性和粗化率无关性.Abstract: The Energy-Minimization Multi-Scale (EMMS) theory has been introduced into the multiphase particle-in-cell (MP-PIC) method to establish the heterogeneous EMMS solid stress model to account for the effect of non-uniform solid distribution. However, the calculation process is very complex and also very time consuming for this heterogeneous solid stress model. The expression of the heterogeneous EMMS solid stress can be obtained by manual fitting method. However, the fitting variable describes heterogeneous solid distribution as well as the fitting function describe the shape of solid stress are required for manually fitting. Since the heterogeneous solid stress function is highly nonlinear in nature, the fitting precision is not high enough for the manually fitting model. And there is an obvious deviation between the fitting correlation and the original EMMS solid stress, because it is hard to find out an appropriate parameter to characterize the heterogeneous solid concentration distribution as well as to find out an appropriate fitting function. In order to solve the above problems, an artificial neutral network (ANN) based machine learning method was proposed to avoid the characterization of the local distribution of solid volume fraction. Subsequently an ANN solid stress model which accounts for the detailed distribution of particle concentration was proposed to improve the fitting accuracy. Firstly, a two-marker based ANN solid stress model was established based on local particle concentration and particle non-uniform distribution index. Further, particle concentrations in the current cell and its neighboring cells were arrayed to represent the particle concentration distribution, thus to establish the ANN solid stress model based on particle concentration distribution. Then, the two models are compared with the EMMS solid stress model, and the effects of grid resolution and coarse-graining ratio on the model are also tested. The simulation results predicted with ANN model agreed well with that of the EMMS solid stress model, and the dependence of simulation results on grid resolution and coarse-graining ratio was also reduced.
-
引 言
早在《天工开物》(宋应星, 1637)中就记载了关于流态化的理论, 通过反复颠簸的方式[1], 使得谷物具有一定的流动性, 从而达到谷物和石子分离的效果. 在现代工业中, 尤其是化工、材料、环保等领域, 小到固体物料的加工、运输, 大到工业过程中的干燥、燃烧、合成等过程都用到了流态化理论[2-6]. 近些年, 随着计算技术的迅猛发展, 数值模拟成为研究流态化反应器中流动反应的重要方法[7], 而准确的模型是成功进行流态化数值模拟的关键[8-9].
流化模拟中最重要的两个模型分别是表征气固之间相互作用的曳力模型和表征颗粒间相互作用的固相应力模型. 近20年来的研究多关注于曳力模型的研究[10-14], 少有针对固相应力模型的研究. 固相应力模型也是MP-PIC方法中的主要模型之一. 以往的模拟实践中通常采用传统的均匀固相应力模型, 研究发现均匀模型对模拟结果有显著影响, 并且具有较大的经验性[15]. 基于均匀模型存在的问题, 姜勇[16-17]提出将EMMS (energy-minimization multi-scale)理论引入固相应力模型建模过程, 他假定(1)在气固流态化中介尺度结构普遍存在; (2)网格中的局部流场可以划分为多个密相与一个稀相; (3)每一个计算粒子可以看做一个密相; (4)每一个计算粒子周围都有稀相绕流; (5)在每一个子相中, 颗粒呈均匀分布. 在此基础上, 通过使用EMMS方法对颗粒流体系统进行分析, 对气固流动进行多尺度分解, 获取每个子相的体积分数与固相应力, 根据子相体积分数进行加权加和就得到带有介尺度结构信息的EMMS固相应力, 从而发展了考虑流化床中非均匀结构的EMMS固相应力模型. 但在实际模拟中对EMMS固相应力模型进行完整求解所需时间较长, 又提出对EMMS固相应力进行了拟合获得固相应力关联式. 将网格的固含率和表征颗粒浓度分布[18-20]的参数作为自变量, 以EMMS固相应力作为函数值, 进行非线性拟合. 但这样的拟合效果较差, 其模拟结果相比原模型对网格和粗化率存在显著的依赖性. 究其原因, 从现象上看在中等颗粒浓度范围内, 对于相同的固含率和颗粒浓度分布参数, 原始EMMS固相应力数据不唯一, 而是具有一定的“厚度”, 这就预示着采用当前二元变量进行拟合难以达到理想精度; 从原理上分析, 拟合所采用的颗粒浓度分布参数和拟合函数难以准确描述实际的非均匀颗粒分布.
所以本文拟通过构建人工神经网络模型对EMMS固相应力进行拟合, 发挥人工神经网络模型无需预先给定拟合函数[21-25], 且便于进行高维非线性拟合的优点[26], 发展具有高精度的非均匀固相应力ANN(artificial neural network)模型. 为此, 首先利用以往研究提出的表征颗粒浓度分布的参数作为输入变量, 通过神经网络模型提高拟合精度; 进一步, 放弃采用某个参数来表征颗粒分布的方法, 而直接以颗粒浓度分布序列为自变量, 通过神经网络模型对EMMS固相应力进行高维非线性建模, 以实现更为精确的EMMS固相应力拟合; 最终将训练所得的非均匀固相应力ANN模型用于循环流化床提升管的MP-PIC模拟, 对其进行验证. 本研究将人工神经网络引入到非均匀固相应力模型中, 为今后的非均匀结构建模和模拟计算提供了一种新的思路和方向.
1. 非均匀固相应力模型中颗粒分布的表征方式
为了便于模型的描述、发展和验证, 本文研究均基于二维规则网格进行. 但需要说明的是, 根据本文的建模思想, 针对颗粒浓度分布的表征方法同样适用于三维情形, 相关的非均匀固相应力AI模型也可以很方便地拓展到三维条件下. 这里首先给出了本文用到的两种表征颗粒浓度分布的方法.
1.1 颗粒非均匀分布指数
以往研究中提出利用表征颗粒浓度非均匀分布的指数δ和固相体积分数εs与非均匀固相应力ps进行关联, 其中非均匀分布指数δ的计算示意图如图1所示. 以网格(i, j)表示当前网格, 颗粒浓度非均匀分布指数的基本思想是用相邻网格(i + ii, j + jj)的颗粒浓度(此处ii, jj为−1或1)和反距离权重法来估计当前网格的颗粒浓度值, 用估计值与真实值的偏差程度来表征颗粒浓度分布非均匀性的大小. 最后可算得当前网格的颗粒浓度非均匀分布指数δ, 即
$$ \left.\begin{split} &\delta = \frac{{\varepsilon }_{s\text{,estimate}}-{\varepsilon }_{s}}{{\varepsilon }_{s}} = \frac{\displaystyle\sum _{ii = -1}^{1}\displaystyle\sum _{jj = -1}^{1} \frac{{\varepsilon }_{s,i + ii,j + jj}}{w{d}_{(i + ii,j + jj)-(i,j)}^{2}}-{\varepsilon }_{s}}{{\varepsilon }_{s}}\\ &w = \sum _{ii = -1}^{1} \sum _{jj = -1}^{1} \frac{1}{{d}_{(i + ii,j + jj)-(i,j)}^{2}} \end{split}\right\}$$ (1) 1.2 颗粒浓度分布的直接表征
本文提出的颗粒浓度分布的直接表征方法也可以由图1所示的网格结构表示. 假定中心网格(i, j)为当前网格, 将其标记为cell1; 依据从左到右, 从下到上的规则, 依次标记cell1的周边网格分别为cell2, cell3, ···, cell9 (图中的数字即代表相应的网格下标), 网格所对应的矩阵坐标分别为(i−1, j−1), (i, j−1),···, (i + 1, j + 1), 其所对应的网格颗粒浓度也分别为εs1, εs2, ···, εs9. 由于定义了1 ~ 9下标所对应的具体位置模板, 那么一组颗粒浓度数列就唯一定义了该局部网格结构的颗粒浓度分布. 采用分布(εs1, εs2,···, εs9)来近似当前网格cell1中的颗粒浓度分布. 相比上节中的颗粒浓度非均匀分布指数, 这种颗粒浓度分布的直接表征方式显然更为准确, 但是也在拟合模型中引入了更多的自变量. 形式上看是引入了更多自变量; 从物理上来看, 则是引入了一种更为准确的描述颗粒浓度分布的方法, 可望提高模型精度. 对于传统的拟合方法来说, 更多自变量的引入需要进行高维非线性拟合, 无论是拟合函数的构造还是拟合过程都是难以进行的, 而这一难点对于人工神经网络模型来说是非常容易克服的.
2. 基于人工神经网络的非均匀固相应力模型
BP (back propagation)神经网络[27-29]是当前应用最为广泛的神经网络模型, 本文工作也基于该网络模型开展. 基于人工神经网络的机器学习方法通常包含: 准备数据集、构建人工神经网络模型、在数据集上训练人工神经网络模型以获得满足训练精度要求的网络参数.
2.1 数据集
首先是准备用于训练人工神经网络模型的数据集. 利用MP-PIC方法耦合文献中的EMMS固相应力模型对经典提升管实验进行模拟, 待模拟稳定后每个时间步输出一个两相流场和固相应力场. 从上述流场和固相应力场中随机采样得到时刻t网格i内颗粒的信息, 包括固相应力ps和颗粒浓度εs, 以及周边网格颗粒浓度. 重复上述过程能够获得相应的数据集. 其中, 将固相应力ps、颗粒浓度εs以及通过式(1)计算得到的颗粒非均匀分布指数δ组成的数据集1用于基于颗粒非均匀分布指数的(双变量) ANN固相应力模型建模, 简称ANN model 1; 而将固相应力和颗粒浓度分布序列(εs1, εs2, ···, εs9) 组成的数据集2用于基于颗粒浓度分布的(多变量)ANN固相应力模型建模, 简称ANN model 2.
2.2 人工神经网络模型
本文构建的神经网络模型由输入层、隐藏层和输出层组成, 对于双变量的ANN固相应力模型来说, 输入层由颗粒浓度εs和颗粒非均匀分布指数δ两个变量组成, 输出变量为固相应力ps. 先设定一个较为简单的中间结构: 有两个隐藏层, 每层由两个神经元组成, 利用Relu函数作为激活函数、Adam作为优化器、平均绝对值误差(mean absolute error, MAE)作为损失函数; 经过训练后发现预测值与真实值存在一定的误差. 进一步, 尝试增加每个隐藏层的神经元数继续进行训练, 发现随着神经元数的不断增加, 拟合精度也不断提高, 如图2所示. 其中决定系数R2和平均绝对值误差的定义分别为
$$ {R}^{2} = \frac{\displaystyle\sum _{i = 1}^{n} {\left({\widehat{y}}_{i}-{\bar {y}}_{i}\right)}^{2}}{\displaystyle\sum _{i = 1}^{n} {\left({y}_{i}-{\bar {y}}_{i}\right)}^{2}} $$ (2) $$ L_1 = \frac{1}{n}\displaystyle\sum _{i = 1}^{n} \left|{y}_{i}-{\widehat{y}}_{i}\right| $$ (3) 其中,
$ {y}_{i} $ ,$ {\widehat{y}}_{i} $ 和${\bar {y}}_{i}$ 分别为第i个样本EMMS模型得到的固相应力值, 第i个样本使用人工神经网络模型预测的固相应力值, EMMS模型得到的所有固相应力值的平均值. 根据图2可确定每个隐藏层包含最佳的神经元数为8, 由此得到图3所示的人工神经网络模型图以及ANN固相应力模型的关联式, 即$$ \left.\begin{split} &{y}_{j} = f\left({\sum }_{i = 1}^{{n}_{X}}{w}_{i,j}{x}_{i} + {b}_{j}\right) \\ &{z}_{k} = f\left({\sum }_{j = 1}^{{n}_{y}}{w}_{j,k}{y}_{j} + {b}_{k}\right) \\ &{p}_{s,{\rm{ANN}}} = {\sum }_{k = 1}^{{n}_{z}}{w}_{k}{x}_{k} + b \\ &{n}_{X} = \left\{\begin{array}{l}2,\\ 9,\end{array}\right.\begin{array}{c}{\rm{ANN}\;{{\rm{model}}}}\;1\\ {\rm{ANN}\;{{\rm{model}}}}\;2\end{array} \\ &n_{y} = 8\text{, }n_{z} = 8 \end{split}\right\}$$ (4) 最终, 应用EMMS固相应力数据集1对人工神经网络模型进行训练后, 可得到基于双变量的ANN固相应力模型(X = εs, δ). 应用EMMS固相应力数据集2, 对该人工神经网络模型进行训练后, 可得到基于颗粒浓度分布的ANN固相应力模型(X = εs1, εs2, ···, εs9).
图4比较了人工拟合模型、人工神经网络模型计算的固相应力与完整求解EMMS固相应力模型算得的固相应力. 图4(a) ~ 图4(c)的横坐标均是EMMS固相应力,纵坐标分别是人工拟合模型、基于双变量的ANN固相应力模型和基于颗粒浓度分布的ANN固相应力模型算得的固相应力, 上虚线表示有99.9%的数据在此线之下, 下虚线表示有99.9%的数据在此线之上. 图4(d)为3个模型的误差比较. 整体看来, 图4(a) ~ 图4(c)两条虚线表示的上下偏差从人工拟合模型到ANN model 1再到ANN model 2呈现出逐渐缩小的趋势. 人工拟合模型上偏差大于下偏差, 即该模型更倾向于高估固相应力. 而两种神经网络模型下偏差大于上偏差, 倾向于低估固相应力. 固相应力位于2 ~ 8 kPa的范围时, 两种ANN固相应力模型误差均明显小于人工拟合模型, 其数据点更加均匀地分布在对角线的左右两侧. 当固相应力小于2 kPa时, 3种模型的误差相似且都较小, 人工拟合模型也具有较高的精度. 在9 kPa以上的范围3种模型误差值逐渐增大, 这种情况可能是由于这部分数据集占整个数据集的比重较小, 使得损失函数对此误差不够敏感, 从而导致拟合误差值相对较大. 整体看来, ANN model 2精度高于ANN model 1, 且两者精度均高于人工拟合模型. 综合上述分析, 可以说明放弃特定拟合函数而使用人工神经网络模型对EMMS固相应力进行拟合可以获得更高的拟合精度; 在此基础上, 采用颗粒浓度序列来直接表征颗粒浓度分布的方法则能进一步提升模型的拟合精度.
图 4 (a)人工拟合模型与EMMS固相应力模型的对比; (b) ANN model 1与EMMS固相应力模型的对比; (c) ANN model 2与EMMS固相应力模型的对比; (d)模型误差Figure 4. (a) Comparison between the manually fitted model and EMMS solid stress model; (b) Comparison between the ANN model 1 and EMMS solid stress model; (c) Comparison between the ANN model 2 and model EMMS solid stress model; (d) Model errors3. ANN固相应力模型验证
本节将以上发展的ANN固相应力模型耦合到MP-PIC方法中对经典提升管实验进行模拟. MP-PIC方法的控制方程[30-34]见表1.
表 1 MP-PIC方法的控制方程和本构模型Table 1. Governing equations and constitutive models of MP-PIC methodDescription Equation governing equations for gas phase $\begin{gathered}\dfrac{\partial }{\partial t}\left({\varepsilon }_{g}{\rho }_{g}\right) + \nabla \cdot \left({\varepsilon }_{g}{\rho }_{g}{u}_{g}\right) = 0\\ \dfrac{\partial }{\partial t}\left({\varepsilon }_{g}{\rho }_{g}{u}_{g}\right) + \nabla \cdot \left({\varepsilon }_{g}{\rho }_{g}{u}_{g}{u}_{g}\right) = -{\varepsilon }_{g}\nabla {p}_{g} + \nabla \cdot \left({\varepsilon }_{g}{\tau }_{g}\right) + {\varepsilon }_{g}{\rho }_{g}g-F\\ {\tau }_{g} = {\mu }_{g}\left[\nabla {u}_{g} + {\left(\nabla {u}_{g}\right)}^{{\rm{T}}}\right]-\dfrac{2}{3}{\mu }_{g}\left(\nabla \cdot {u}_{g}\right)I\end{gathered}$ particle distribution function $\dfrac{\partial \phi }{\partial t} + \nabla \cdot \left({u}_{p}\varphi \right) + {\nabla }_{ {u}_{p} }\cdot \left({a}_{p}\phi \right) = 0$ particle speed ${u}_{p} = \dfrac{{\rm{d}}{x}_{p} }{{\rm{d}}t}$ particle acceleration ${a}_{p} = {\beta }_{p}\left({u}_{g,p}-{u}_{p}\right)-\dfrac{\nabla {p}_{gp} }{ {\rho }_{s} } + g-\dfrac{\nabla {p}_{s,p} }{ {\varepsilon }_{s}{\rho }_{s} }$ gas-solid inter-phase momentum exchange rate $F = \displaystyle\sum _{k = 1}^{ {n}_{T} } \frac{ {n}_{p}^{k}{V}_{p}^{k} }{ {V}_{c} }{\beta }_{p}^{k}\left({u}_{g}^{k}-{u}_{p}^{k}\right)$ homogeneous drag model ${\beta }_{0} = \dfrac{3}{4}{C}_{d0}\dfrac{ {\varepsilon }_{g}{\rho }_{g}\left|{u}_{gp}-{u}_{p}\right|}{ {d}_{p} }{\varepsilon }_{g}^{-2.7}$ heterogeneous drag model[35] $\begin{gathered}{H}_{D} = a{\left(R{e}_{p} + b\right)}^{c}\\ {\beta }_{p} = {\beta }_{p0}{H}_{D}\end{gathered}$ homogeneous solid force model ${p}_{s0} = \dfrac{ {p}_{s}^{*}{\varepsilon }_{s}^{\alpha } }{{\rm{max}}\left[{\varepsilon }_{s,cp}-{\varepsilon }_{s},\delta \left(1-{\varepsilon }_{s}\right)\right]}$ EMMS solid stress model ${p}_{s,{\rm{EMMS}}} = \displaystyle\sum _{k = 1}^{ {n}_{T} } \frac{ {f}_{k}{p}_{s}^{*}{\left({\varepsilon }_{sc}^{k}\right)}^{\alpha } }{ {\rm{max} }\left[{\varepsilon }_{s,cp}-{\varepsilon }_{sc}^{k},\delta \left(1-{\varepsilon }_{sc}^{k}\right)\right]}$ artificial fitting model of EMMS solid stress $f\left({\varepsilon }_{s},{\delta }_{s}\right) = -3.273{\varepsilon }_{s}^{0.711\;5} + 1.370-{\delta }_{s}$
${ {p}_{s,{\rm{Fit}}} }\left( { {\varepsilon }_{s} },{ {\delta }_{s} } \right)\text{=}\left\{\begin{aligned} & 1.164\varepsilon _{s}^{0.986\;2}\Biggr/\left[ { {\left( \frac{ { {\varepsilon }_{s} }+1}{0.272\;1} \right)}^{-0.778\;7} }-{ {\varepsilon }_{s} } \right],f\left( { {\varepsilon }_{s} },{ {\delta }_{s} } \right)\leqslant 0 \\ & { { {\rm{e} } }^{\left( -9.373\delta _{s}^{3}-0.454\;5{ {\delta }_{s} }-7.014 \right){ {\left( { {\delta }_{s} }+3.319{ {\varepsilon }_{s} } \right)}^{-3.61} }+9.459} },f\left( { {\varepsilon }_{s} },{ {\delta }_{s} } \right) > 0 \\ \end{aligned} \right.$模拟中采用的提升管的示意图如表2, 床高为10.50 m, 直径为0.09 m, 参考Li等[36]的实验装置建立. 气相和固相从提升管底端进入其内部, 从提升管顶部流出, 颗粒流出后又被送回底部入口用以模拟颗粒循环. 提升管的两侧设定为壁面. 气固相的边界条件与文献[2]中提到的边界条件设置保持一致, 即气相的边界条件设置为速度入口、压力出口及壁面无滑移; 固相的边界条件设置为从出口到入口的单周期边界、壁面及入口处为碰撞反弹边界, 用颗粒−壁面碰撞恢复系数来计算颗粒速度. 初始在提升管内均布一定量颗粒, 全床平均固相体积分数εs0为0.058. 在所有模拟中, 均使用考虑介尺度结构的非均匀曳力模型, 相关模型如表1所示. 使用本文发展的ANN固相应力模型封闭固相应力, 作为对比, 原EMMS固相应力模型和采用人工拟合得到的固相应力模型也分别用于相同工况的模拟. 模拟所用到的控制方程和本构关系汇总于表1中. 模拟中用到的计算参数汇总于表2中. 根据模拟结果可知(如图5(b)所示), 大约经过6 s, 提升管内计算结果达到稳定. 因此, 流场中所有的平均量的统计时间都从第6 s开始, 到25 s结束, 可以确保得到数值稳定的平均性质. 参考文献中提及的粗化率和网格分辨率参数[16], 本文选择同样的参数进行比较分析. 已训练好的ANN model 1和ANN model 2模拟提升管内的气固两相流25 s所耗费的时间分别为468 s和474 s, EMMS固相应力模拟的时间524 s. 不难发现基于神经网络的固相应力模型相比EMMS固相应力模型有更高的计算效率.
表 2 算例参数Table 2. Simulation parametersSchematic drawing Parameter Value diameter, D/m 0.09 height, H/m 10.50 operating temperature, Top/K 298 operating pressure, Pop/Pa 1.01 × 105 particle diameter, dp/μm 54 particle density, ρs/(kg·m−3) 970 gas density, ρg/(kg·m−3) 1.185 gas viscosity, μg/(Pa·s) 1.90 × 10−5 superficial gas velocity, Ugo/(m·s−1) 1.52 incipient average solid volume fractions, εs0 0.058 solid volume fraction at close pack, εs,cp 0.60 solid stress parameter, ps*/Pa 100 solid stress parameter, α 2 particle–wall restitution coefficient, ew 0.90 coarse-grained ratio, np 100, 200, 500, 1000 number of grids, I × J 20 × 200, 20 × 400, 30 × 500, 40 × 600 time step, Δt/s 5 × 10−4 simulation time, t/s 25 3.1 不同固相应力模型的比较
图5对比了本文发展的两种ANN固相应力模型与EMMS固相应力模型和人工拟合模型的模拟结果, 同时给出了实验测量值. 需要说明的是, 两种ANN固相应力模型和人工拟合模型都是通过对EMMS固相应力模型进行拟合得到, 因此, 其目标是能最大程度地复现EMMS固相应力; 与之相应, 基于ANN固相应力的模拟结果也应当尽量复现EMMS固相应力的模拟结果. 图5(a)为不同固相应力模型预测的轴向孔隙率分布. 整体来看, 所有固相应力模型均能预测实验测量值所展现的上稀下浓的提升管轴向分布特征. 人工拟合模型算得的孔隙率在3 m以上区域明显高于EMMS固相应力的计算结果, 而在3 m以下区域明显低于EMMS固相应力的计算结果. 基于双变量的ANN模型与人工拟合模型的计算结果相近, 但是能显著改善3 m以上区域的预测结果, 预测的孔隙率更接近EMMS固相应力模型. 基于颗粒浓度分布的ANN模型则在基于双变量的ANN模型基础上进一步改善了模拟结果, 在0 ~ 8 m的位置更接近EMMS 模型的计算结果, 而在8 ~ 10 m的位置处尽管其不如ANN model 1准确, 但差距不大. 相比较而言, 基于颗粒浓度分布的ANN固相应力模型与EMMS固相应力模型的计算结果吻合最好. 总的来说, 放弃了用某个特定参数来表征颗粒分布, 而直接采用颗粒浓度分布作为非均匀固相应力的自变量后, ANN固相应力模型的预测值与原始EMMS固相应力模型更接近了, 表现为图中所示的基于颗粒浓度分布的ANN模型曲线相比基于双变量的ANN模型曲线更接近EMMS模型曲线. 这也证明了本文提出的采用颗粒浓度分布替代以往的颗粒非均匀分布指数δ是合理的, 而且也有助于非均匀固相应力模型精度的提升. 图6(b)对比了瞬时固相流率的模拟结果, 由图可见各固相应力模型预测的固相通量均在5 s以后达到稳定, 稳定后各固相通量曲线均围绕实验值波动, 表明各模型均能预测出合理的固相通量. 稍有不同的是在0 ~ 5 s的范围内, 相比基于双变量的ANN模型预测结果, 基于颗粒浓度分布的ANN模型得到的瞬时固相通量变化趋势与EMMS模型的结果更为接近. 表明基于颗粒浓度分布的ANN固相应力模型对EMMS固相应力模型有较好的还原度.
3.2 网格分辨率对孔隙率分布的影响
这里研究了使用两种ANN固相应力模型时, 不同网格分辨率对模拟预测的轴向和径向孔隙率分布的影响. 图6对比了当粗化率np = 1000时, 网格分辨率在20 × 200 ~ 40 × 600之间, ANN固相应力模型模拟所得提升管内孔隙率轴向分布. 图中每个算例除网格分辨率外, 其他设置均相同. 如图6(a)所示, 随着网格分辨率的变化, 双变量ANN模型模拟结果变化比较明显; 而基于颗粒浓度分布的ANN模型预测得到的轴向空隙率分布变化不明显, 表明该模型相比于双变量ANN模型对网格分辨率不敏感, 如图6(b)所示. 在当前的网格分辨率范围内, 其模拟结果能够保持一定的网格无关性, 即使使用较低的网格分辨率如20 × 200时, 基于颗粒浓度分布的ANN模型依然能够复现高分辨率网格预测的轴向空隙率分布. 当然, 图6(b)结果显示不同网格结果仍存在微小偏差, 这与以往报道的流化床模拟研究类似[37], 可能是由流化床内存在随时空演化的非均匀流动结构所致.
图7对比了不同网格分辨率下两种ANN固相应力模型算得的提升管1.5 m高处的径向孔隙率的分布, 作为参考, 同时给出了经验模型算得的径向分布[38]. 总体来看, 使用两种ANN固相应力模型时, 网格分辨率对孔隙率径向分布影响不大, 不同分辨率下的模拟结果均呈中间颗粒浓度稀、壁面颗粒浓度稠密的“环−核结构”分布. 图7(a)和图7(b)相比较而言, 双变量ANN模型的预测结果随网格分辨率变化较大, 基于颗粒浓度分布的ANN模型预测的孔隙率径向分布变化更小也更接近经验模型计算值. 尽管不同网格分辨率的结果存在一些差别, 但总体变化不大, 即使采用较粗的网格也能获得与经验模型计算值相近的结果. 上述结果表明, ANN固相应力模型, 特别是基于颗粒浓度分布的ANN模型在较低网格分辨率(也即较粗的网格)下也能够保证计算的准确性和稳定性.
3.3 粗化率对孔隙率分布的影响
MP-PIC方法中, 粗化率np代表了一个计算粒子中包含的真实颗粒数, 它是粗粒化方法中的一个重要参数. 通常, 它与网格分辨率类似, 会影响计算结果. 而实际计算中, 人们希望采用更大的np来提升计算效率. 因此, 这里研究了使用两种ANN固相应力模型时, 不同粗化率对模拟预测的轴向和径向孔隙率分布的影响.
图8对比了网格分辨率为20 × 200, 粗化率为100 ~ 2000时, ANN固相应力模型模拟所得提升管内孔隙率轴向分布. 图8(a)为双变量ANN固相应力模型预测的轴向孔隙率分布, 如图所示, 在当前粗化率范围内, 时均轴向孔隙率的分布基本保持一致, 说明在一定范围内当前模型对粗化率np的变化不敏感; 图8(b)为基于颗粒浓度分布的ANN固相应力模型预测的模拟结果, 如图所示, 孔隙率分布曲线几乎不随粗化率变化, 说明该模型计算结果对np不敏感.
图9对比了不同粗化率下ANN固相应力模型模拟所得提升管1.5 m高处的径向孔隙率的分布. 由图可见, 两个ANN模型在不同粗化率下的模拟结果均能捕捉到中间颗粒稀少、边壁颗粒稠密的“环−核结构”, 其趋势与经验模型的计算结果一致. 对比图9(a)和图9(b)可见, 双变量ANN模型在不同粗化率下预测的径向分布变化相对较大, 而基于颗粒浓度分布的ANN模型预测结果随粗化率的变化相对较小也更接近经验模型计算值, 这说明基于颗粒浓度分布的ANN模型对粗化率的选取更不敏感. 因此可以在保证结果准确性的同时, 选用较大的粗化率来降低计算成本. 这也体现了基于颗粒浓度分布的ANN固相应力模型的一个优点, 也即将来用于大规模工业装置模拟时, 可以选择较大的粗化率来提升计算效率、节省计算时间和资源.
4. 结 论
本文在非均匀固相应力模型建模中引入基于人工神经网络的机器学习方法. 通过构建恰当的人工神经网络模型, 将EMMS固相应力与非均匀的颗粒浓度分布建立关联, 获得了ANN固相应力模型. 并将所得到的ANN固相应力模型应用于MP-PIC方法中对经典提升管进行模拟, 将模拟结果与实验值和原始EMMS固相应力模型的预测值进行了比较, 同时通过对比不同网格分辨率和粗化率工况的模拟, 测试了ANN固相应力模型的性能. 为今后的非均匀模型建模开辟了新的思路, 获得主要结论如下.
(1) 基于局部颗粒浓度和颗粒非均匀分布指数建立了双变量ANN固相应力模型, 证明了人工神经网络模型用于非均匀固相应力模型建模是可行的.
(2) 用当前网格及其周边网格的颗粒浓度按一定规则组成颗粒浓度序列来直接表征颗粒浓度分布, 根据此分布建立了考虑颗粒浓度分布的ANN固相应力模型.
(3) 双变量ANN固相应力模型对EMMS固相应力模型的拟合精度高于人工拟合模型; 基于颗粒浓度分布的ANN固相应力模型的拟合精度又高于双变量ANN固相应力模型.
(4) 双变量ANN固相应力模型的模拟结果相比人工拟合模型的更接近原始EMMS固相应力模型的模拟结果, 能提高模拟精度; 基于颗粒浓度分布的ANN固相应力模型的模拟结果比双变量ANN固相应力模型的更接近原始EMMS固相应力模型的模拟结果, 模拟精度进一步提高.
(5) 在本文测试的网格分辨率和粗化率范围内, 基于颗粒浓度分布的ANN固相应力模型的预测结果具备较好网格无关性和粗化率无关性.
还需说明的是, 本文用于获取ANN固相应力模型的固相应力源数据由EMMS固相应力模型产生, 所以本文提出的基于神经网络的固相应力模型主要适用快速流态化模拟, 对于其在其他气固系统中的适用性, 尚需开展更多的模拟和验证研究. 此外, CFD-DEM (computational fluid dynamics-discrete element method)方法直接追踪颗粒碰撞过程, 相比MP-PIC方法在模拟颗粒间受力时精度更高, 进一步研究还可考虑由CFD-DEM模拟产生固相应力源数据用于本文的ANN固相应力建模, 有望进一步提高模拟精度.
-
图 4 (a)人工拟合模型与EMMS固相应力模型的对比; (b) ANN model 1与EMMS固相应力模型的对比; (c) ANN model 2与EMMS固相应力模型的对比; (d)模型误差
Figure 4. (a) Comparison between the manually fitted model and EMMS solid stress model; (b) Comparison between the ANN model 1 and EMMS solid stress model; (c) Comparison between the ANN model 2 and model EMMS solid stress model; (d) Model errors
表 1 MP-PIC方法的控制方程和本构模型
Table 1 Governing equations and constitutive models of MP-PIC method
Description Equation governing equations for gas phase $\begin{gathered}\dfrac{\partial }{\partial t}\left({\varepsilon }_{g}{\rho }_{g}\right) + \nabla \cdot \left({\varepsilon }_{g}{\rho }_{g}{u}_{g}\right) = 0\\ \dfrac{\partial }{\partial t}\left({\varepsilon }_{g}{\rho }_{g}{u}_{g}\right) + \nabla \cdot \left({\varepsilon }_{g}{\rho }_{g}{u}_{g}{u}_{g}\right) = -{\varepsilon }_{g}\nabla {p}_{g} + \nabla \cdot \left({\varepsilon }_{g}{\tau }_{g}\right) + {\varepsilon }_{g}{\rho }_{g}g-F\\ {\tau }_{g} = {\mu }_{g}\left[\nabla {u}_{g} + {\left(\nabla {u}_{g}\right)}^{{\rm{T}}}\right]-\dfrac{2}{3}{\mu }_{g}\left(\nabla \cdot {u}_{g}\right)I\end{gathered}$ particle distribution function $\dfrac{\partial \phi }{\partial t} + \nabla \cdot \left({u}_{p}\varphi \right) + {\nabla }_{ {u}_{p} }\cdot \left({a}_{p}\phi \right) = 0$ particle speed ${u}_{p} = \dfrac{{\rm{d}}{x}_{p} }{{\rm{d}}t}$ particle acceleration ${a}_{p} = {\beta }_{p}\left({u}_{g,p}-{u}_{p}\right)-\dfrac{\nabla {p}_{gp} }{ {\rho }_{s} } + g-\dfrac{\nabla {p}_{s,p} }{ {\varepsilon }_{s}{\rho }_{s} }$ gas-solid inter-phase momentum exchange rate $F = \displaystyle\sum _{k = 1}^{ {n}_{T} } \frac{ {n}_{p}^{k}{V}_{p}^{k} }{ {V}_{c} }{\beta }_{p}^{k}\left({u}_{g}^{k}-{u}_{p}^{k}\right)$ homogeneous drag model ${\beta }_{0} = \dfrac{3}{4}{C}_{d0}\dfrac{ {\varepsilon }_{g}{\rho }_{g}\left|{u}_{gp}-{u}_{p}\right|}{ {d}_{p} }{\varepsilon }_{g}^{-2.7}$ heterogeneous drag model[35] $\begin{gathered}{H}_{D} = a{\left(R{e}_{p} + b\right)}^{c}\\ {\beta }_{p} = {\beta }_{p0}{H}_{D}\end{gathered}$ homogeneous solid force model ${p}_{s0} = \dfrac{ {p}_{s}^{*}{\varepsilon }_{s}^{\alpha } }{{\rm{max}}\left[{\varepsilon }_{s,cp}-{\varepsilon }_{s},\delta \left(1-{\varepsilon }_{s}\right)\right]}$ EMMS solid stress model ${p}_{s,{\rm{EMMS}}} = \displaystyle\sum _{k = 1}^{ {n}_{T} } \frac{ {f}_{k}{p}_{s}^{*}{\left({\varepsilon }_{sc}^{k}\right)}^{\alpha } }{ {\rm{max} }\left[{\varepsilon }_{s,cp}-{\varepsilon }_{sc}^{k},\delta \left(1-{\varepsilon }_{sc}^{k}\right)\right]}$ artificial fitting model of EMMS solid stress $f\left({\varepsilon }_{s},{\delta }_{s}\right) = -3.273{\varepsilon }_{s}^{0.711\;5} + 1.370-{\delta }_{s}$
${ {p}_{s,{\rm{Fit}}} }\left( { {\varepsilon }_{s} },{ {\delta }_{s} } \right)\text{=}\left\{\begin{aligned} & 1.164\varepsilon _{s}^{0.986\;2}\Biggr/\left[ { {\left( \frac{ { {\varepsilon }_{s} }+1}{0.272\;1} \right)}^{-0.778\;7} }-{ {\varepsilon }_{s} } \right],f\left( { {\varepsilon }_{s} },{ {\delta }_{s} } \right)\leqslant 0 \\ & { { {\rm{e} } }^{\left( -9.373\delta _{s}^{3}-0.454\;5{ {\delta }_{s} }-7.014 \right){ {\left( { {\delta }_{s} }+3.319{ {\varepsilon }_{s} } \right)}^{-3.61} }+9.459} },f\left( { {\varepsilon }_{s} },{ {\delta }_{s} } \right) > 0 \\ \end{aligned} \right.$表 2 算例参数
Table 2 Simulation parameters
Schematic drawing Parameter Value diameter, D/m 0.09 height, H/m 10.50 operating temperature, Top/K 298 operating pressure, Pop/Pa 1.01 × 105 particle diameter, dp/μm 54 particle density, ρs/(kg·m−3) 970 gas density, ρg/(kg·m−3) 1.185 gas viscosity, μg/(Pa·s) 1.90 × 10−5 superficial gas velocity, Ugo/(m·s−1) 1.52 incipient average solid volume fractions, εs0 0.058 solid volume fraction at close pack, εs,cp 0.60 solid stress parameter, ps*/Pa 100 solid stress parameter, α 2 particle–wall restitution coefficient, ew 0.90 coarse-grained ratio, np 100, 200, 500, 1000 number of grids, I × J 20 × 200, 20 × 400, 30 × 500, 40 × 600 time step, Δt/s 5 × 10−4 simulation time, t/s 25 -
[1] 李洪钟, 郭慕孙. 回眸与展望流态化科学与技术. 化工学报, 2013, 64(1): 52-62 (Li Hongzhong, Guo Musun. Review and prospect of fluidization science and technology. CIESC Journal, 2013, 64(1): 52-62 (in Chinese) [2] Zhao X, Jiang Y, Li F, et al. A scaled MP-PIC method for bubbling fluidized beds. Powder Technol, 2022, 404: 117501 doi: 10.1016/j.powtec.2022.117501
[3] Chen XZ, Wang JW. A comparison of two-fluid model, dense discrete particle model and CFD-DEM method for modeling impinging gas-solid flows. Powder Technol, 2014, 254: 94-102 doi: 10.1016/j.powtec.2013.12.056
[4] Cai W, Kong X, Ye Q, et al. Numerical modelling of hydrodynamics of molten salt fluid-particles fluidized beds using CFD-DEM and TFM approaches. Powder Technol, 2022, 410: 117882
[5] Leckner B. Hundred years of fluidization for the conversion of solid fuels. Powder Technol, 2022, 411: 117935
[6] Wang J, Ku X, Lin J. Numerical investigation of biomass pyrolysis performance in a fluidized-bed reactor by a TFM-DEM hybrid model. Chemical Engineering Science, 2022, 260: 117922 doi: 10.1016/j.ces.2022.117922
[7] Chu KW, Chen YX, Ji L, et al. Coarse-grained CFD-DEM study of Gas-solid flow in gas cyclone. Chemical Engineering Science, 2022, 260: 117906 doi: 10.1016/j.ces.2022.117906
[8] Gao X, Li TW, Sarkar A, et al. Development and validation of an enhanced filtered drag model for simulating gas-solid fluidization of Geldart A particles in all flow regimes. Chemical Engineering Science, 2018, 184: 33-51 doi: 10.1016/j.ces.2018.03.038
[9] Sarkar A, Sun X, Sundaresan S. Verification of sub-grid filtered drag models for gas-particle fluidized beds with immersed cylinder arrays. Chemical Engineering Science, 2014, 114: 144-154 doi: 10.1016/j.ces.2014.04.018
[10] Ebrahimi M, Crapper M, Ooi JY. Numerical and experimental study of horizontal pneumatic transportation of spherical and low-aspect-ratio cylindrical particles. Powder Technol, 2016, 293: 48-59 doi: 10.1016/j.powtec.2015.12.019
[11] Miao Z, Kuang SB, Zughbi H, et al. CFD simulation of dilute-phase pneumatic conveying of powders. Powder Technol, 2019, 349: 70-83 doi: 10.1016/j.powtec.2019.03.031
[12] Radl S, Sundaresan S. A drag model for filtered Euler-Lagrange simulations of clustered gas-particle suspensions. Chemical Engineering Science, 2014, 117: 416-425 doi: 10.1016/j.ces.2014.07.011
[13] Yuan ZH, Wang SY, Shao BL, et al. Investigation on effect of drag models on flow behavior of power-law fluid-solid two-phase flow in fluidized bed. Particuology, 2022, 70: 43-54 doi: 10.1016/j.partic.2022.01.008
[14] Yang N, Chen JH, Ge W, et al. A conceptual model for analyzing the stability condition and regime transition in bubble columns. Chemical Engineering Science, 2010, 65(1): 517-526 doi: 10.1016/j.ces.2009.06.014
[15] Kadyrov T, Li F, Wang W. Impacts of solid stress model on MP-PIC simulation of a CFB riser with EMMS drag. Powder Technol, 2019, 354: 517-528 doi: 10.1016/j.powtec.2019.06.018
[16] 姜勇. 基于MP-PIC方法的流态化反应器快速模拟研究. [博士论文]. 北京: 中国科学院大学, 2020 Jiang Yong. Rapid simulation of fluidized reactor based on MP-PIC method. [PhD Thesis]. Beijing: University of Chinese Academy of Sciences, 2020 (in Chinese))
[17] Jiang Y, Li F, Ge W, et al. EMMS-based solid stress model for the multiphase particle-in-cell method. Powder Technol, 2020, 360: 1377-1387 doi: 10.1016/j.powtec.2019.09.031
[18] Harris A, Davidson J, Thorpe R. The prediction of particle cluster properties in the near wall region of a vertical riser (200157). Powder Technol, 2002, 127(2): 128-143 doi: 10.1016/S0032-5910(02)00114-6
[19] Wang J, Ge W, Li J. Eulerian simulation of heterogeneous gas–solid flows in CFB risers: EMMS-based sub-grid scale model with a revised cluster description. Chemical Engineering Science, 2008, 63(6): 1553-1571 doi: 10.1016/j.ces.2007.11.023
[20] Zou B, Li H, Xia Y, et al. Cluster structure in a circulating fluidized bed. Powder Technol, 1994, 78(2): 173-178 doi: 10.1016/0032-5910(93)02786-A
[21] Aguila-Leon J, Vargas-Salgado C, Chiñas-Palacios C, et al. Energy management model for a standalone hybrid microgrid through a particle Swarm optimization and artificial neural networks approach. Energy Conversion and Management, 2022, 267: 115920 doi: 10.1016/j.enconman.2022.115920
[22] Calisir T, Çolak AB, Aydin D, et al. Artificial neural network approach for investigating the impact of convector design parameters on the heat transfer and total weight of panel radiators. International Journal of Thermal Sciences, 2023, 183: 107845 doi: 10.1016/j.ijthermalsci.2022.107845
[23] Skrypnik A, Shchelchkov A, Gortyshov YF, et al. Artificial neural networks application on friction factor and heat transfer coefficients prediction in tubes with inner helical-finning. Applied Thermal Engineering, 2022, 206: 118049 doi: 10.1016/j.applthermaleng.2022.118049
[24] Kalay E, Boğoçlu ME, Bolat B. Mass flow rate prediction of screw conveyor using artificial neural network method. Powder Technol, 2022, 408: 117757 doi: 10.1016/j.powtec.2022.117757
[25] Zhang K, Zhang Z, Han Y, et al. Artificial neural network modeling for steam ejector design. Applied Thermal Engineering, 2022, 204: 117939 doi: 10.1016/j.applthermaleng.2021.117939
[26] Afandi A, Lusi N, Catrawedarma I, et al. Prediction of temperature in 2 meters temperature probe survey in Blawan geothermal field using artificial neural network (ANN) method. Case Studies in Thermal Engineering, 2022, 38: 102309 doi: 10.1016/j.csite.2022.102309
[27] Pappachan BK, Tjahjowidodo T. Parameter Prediction Using Machine Learning in Robot-Assisted Finishing Process. International Journal of Mechanical Engineering and Robotics Research, 2020, 9(3): 435-440
[28] Gao G, Li Y, Li J, et al. A hybrid model for short-term rainstorm forecasting based on a back-propagation neural network and synoptic diagnosis. Atmospheric and Oceanic Science Letters, 2021, 14(5): 100053 doi: 10.1016/j.aosl.2021.100053
[29] Chen H, Wang Y, Zuo MS, et al. A new prediction model of CO2 diffusion coefficient in crude oil under reservoir conditions based on BP neural network. Energy, 2022, 239: 122286 doi: 10.1016/j.energy.2021.122286
[30] Snider DM. An incompressible three-dimensional multiphase particle-in-cell model for dense particle flows. J Comput Phys, 2001, 170(2): 523-549 doi: 10.1006/jcph.2001.6747
[31] Li F, Song FF, Benyahia S, et al. MP-PIC simulation of CFB riser with EMMS-based drag model. Chemical Engineering Science, 2012, 82: 104-113 doi: 10.1016/j.ces.2012.07.020
[32] Ariyaratne WKH, Ratnayake C, Melaaen MC. Application of the MP-PIC method for predicting pneumatic conveying characteristics of dilute phase flows. Powder Technol, 2017, 310: 318-328 doi: 10.1016/j.powtec.2017.01.048
[33] Kadyrov T, Li F, Wang W. Bubble-based EMMS/DP drag model for MP-PIC simulation. Powder Technol, 2020, 372: 611-624 doi: 10.1016/j.powtec.2020.06.023
[34] Dymala T, Wytrwat T, Heinrich S. MP-PIC simulation of circulating fluidized beds using an EMMS based drag model for Geldart B particles. Particuology, 2021, 59: 76-90 doi: 10.1016/j.partic.2021.07.002
[35] 鲁波娜. 基于EMMS的介尺度模型及其在气固两相流模拟中的应用. [博士论文]. 北京: 中国科学院过程工程研究所, 2009 Lu Bona. EMMS-based meso-scale model and its application in simulating gas-solid two-phase flows. [PhD Thesis]. Beijing: Institute of Process Engineering, Chinese Academy of Sciences, 2009 (in Chinese)
[36] Li J, Kwauk M. Particle-fluid Two-phase Flow: The Energy-minimization Multi-scale Method. Beijing: Metallurgical Industry Press, 1994
[37] Lu B, Wang W, Li J. Eulerian simulation of gas–solid flows with particles of Geldart groups A, B and D using EMMS-based meso-scale model. Chemical Engineering Science, 2011, 66(20): 4624-4635 doi: 10.1016/j.ces.2011.06.026
[38] Wang W, Li J. Simulation of gas–solid two-phase flow by a multi-scale CFD approach —of the EMMS model to the sub-grid level. Chemical Engineering Science, 2007, 62(1-2): 208-231 doi: 10.1016/j.ces.2006.08.017