一种耦合“CFD+深度学习”的气固快速模拟方法
摘要：在计算流体力学领域，深度学习被用于重建流场、预测曳力、求解泊松方程、加速流体模拟等方面的研究。为了加速气固两相流的模拟计算，使用卷积长短时记忆网络对物理量进行预测，并基于LibTorch实现深度学习模型预测与OpenFOAM的耦合。通过与单纯OpenFOAM模拟结果对比，发现深度学习模型预测存在颗粒体积分数不守恒、极小数值预测不准确的问题，先后通过体积分数校正和网格数据过滤消除了前述影响。选取不同的三个物理量组合进行深度学习模型预测以加速CFD计算，在同样选取颗粒体积分数和气体速度的条件下，对比了增加预测颗粒速度和压力对耦合流程计算结果的影响，其中，增加预测颗粒速度的计算结果较为准确，经分析发现这可能与求解器对不同变量的调用顺序有关。
气固两相流广泛存在于化工过程中，为了理解反应器内的气固复杂流动反应过程，计算流体力学的应用日益增多。气固体系的CFD模拟方法主要分为连续和离散两种方法，其中连续方法以求解Navier-Stokes（NS）方程为主，离散方法以求解大量颗粒的运动方程为主。在不同的方法中，体系内网格或颗粒的数量、时间步长、物理模型的复杂度等，对于求解过程影响显著，在建立合理模型的任务之外，来自计算规模与效率的挑战也是重要的考量因素，如何提高计算速度一直是该领域重要的研究课题。近年来，以深度神经网络为代表的深度学习在诸多学科领域取得迅速进展[1-4]，CFD领域也不例外[5-7]。例如，Kim等[8]基于卷积神经网络（convolutional neural network, CNN）提出了一种生成神经网络模型，能够从一组简化参数中生成不发散的二维或三维流体速度场。Fukami等[9]开发了一种混合下采样跳跃连接多尺度模型，能够对二维圆柱尾流的粗糙流场图像进行超分辨率重建。Yang等[10]利用神经网络训练密相流化床能量最小多尺度曳力模型，实现了在多工况下曳力的合理预测。Muralidhar等[11]在使用神经网络训练颗粒曳力模型时，将网络中隐藏层的输出与压力场和速度场联系起来，旨在开发一个受物理启发的深度学习模型。对于气固两相流及更多复杂系统，Guo等[12]提出了一种人工智能的研究范式，拟将介科学原理引入深度学习模型设计中，从而提高模型的可解释性及拓展应用领域。鉴于多相流数值模拟计算速度的局限性和深度学习的巨大潜力，本研究希望借助深度学习来加速气固两相流的模拟计算。
目前，利用深度学习加速流体模拟有诸多方法。Maulik等[13]利用深度学习生成湍流黏度系数的代理模型，该模型允许稳态速度和压力求解器使用更大的松弛因子，使得稳态NS方程的求解速度提高了约5倍。泊松方程的求解是NS方程求解过程中比较耗时的步骤之一[5, 14]，一般基于预处理共轭梯度法（preconditioned conjugate gradient, PCG）进行求解，Tompson等[14]使用CNN替代PCG对二维和三维烟雾模拟进行了加速。Ajuria等[15]也开发了一种基于CNN模型的泊松方程求解器，由于该网络在高Richardson数下结果不稳定，他们还提出了一种基于神经网络与传统雅各比求解器的混合求解策略。Obiols-Sales等[16]利用CNN预测稳态不可压NS方程中部分迭代步后的速度场和压力场，实现了1.9~7.4倍的加速。Jeon等[17]以反应流为背景，模拟微分方程的求解过程，设计了有限体积法神经网络，根据当前时刻的温度场、浓度场和速度场预测下一时刻温度场、浓度场和速度场的变化量。Bazai等[18]在鼓泡床体系中，借助神经网络预测若干时间步后的颗粒体积分数。
在本研究中，提出根据模拟计算得到的前序物理场预测后序物理场的方法对整个模拟计算过程进行加速，研究长时间预测对计算结果的影响，以期较大程度地减少数值求解计算量，加速整个模拟计算过程。对物理场的时序预测可以看作一个时空序列预测问题，即神经网络通过观察大量的时空序列数据来捕捉体系内的动力学规律，并对未来某时刻的空间分布进行预测。时空序列预测问题已在多个领域获得了应用。在天气领域，Shi等[19-20]提出了卷积长短时记忆（convolutional long short-term memory, ConvLSTM）网络，并为降水预测问题提供了一个基准模型。Wang等[21-22]提出了一种新型时空长短时记忆单元，可以同时提取和记忆时间和空间上的特征。在交通领域，Yu等[23]提出一种时空图卷积神经网络对交通流进行时空序列预测，该网络与传统的CNN相比，具有更少的参数，训练速度更快。在海洋领域，Zhou等[24]将ConvLSTM用于海浪高度的预测，具有良好的精度和效率。Gonzalez等[25]提出了一种基于深度学习的非线性模型归约策略，利用卷积自动编码器获得输入数据的低维表示，再使用改进的长短时记忆（long short-term memory, LSTM）网络进行动力学建模，结果显示该网络优于传统的POD-Galerkin降阶模型。值得注意的是，对于不同的时空序列预测问题，根据所解决问题体现出的不同特点，网络设计会有不同的侧重点，例如画面中是否有确定的运动主体、运动主体是否具有刚性以及画面中哪些位置的相关度更高等，都是在网络设计中需要考虑的问题。综上，由于各类流场普遍存在的时空复杂结构，实现其准确预测具有挑战性。
在借助深度学习加速流体模拟的研究中，除了构建一个合理的神经网络架构之外，神经网络与模拟计算软件的耦合，以及网络模型的部署应用，也是研究者需要认真处理的问题。Obiols-Sales等[16]基于TensorFlow与OpenFOAM开发了一个将计算流体力学与神经网络耦合的计算框架CFDNet，利用神经网络预测替代CFD模拟中部分迭代步的计算，以加速稳态RANS模型求解。Weiner等[26]基于Python开发了flowTorch包，可以使研究者在Python环境下更方便地使用OpenFOAM模拟数据进行网络训练以及后处理。Maulik等[27]利用TensorFlow的C接口库将训练好的涡黏模型部署到OpenFOAM中，以加速稳态模拟。Jeon等[28]基于TensorFlow以及OpenFOAM建立了一种机器学习和CFD耦合计算的迁移学习策略，对二维不可压逆流模拟进行加速。Lu等[29]基于LibTorch（PyTorch C++版本）将训练好的曳力模型与MFix相结合，以对比不同曳力模型对模拟结果的影响。
根据以上结果，由于PyTorch深度学习框架易于调试，开源代码丰富，在科学研究领域应用广泛，因此本研究选用此框架将深度学习与OpenFOAM相耦合，实现气固两相流的模拟加速。首先对体系内重要物理量构建深度学习模型，再将深度学习模型预测与OpenFOAM模拟过程进行耦合。将耦合计算得到的结果与单纯OpenFOAM模拟结果进行比较分析，考察结果的可靠性及加速效果，并探究不同物理量预测以及深度学习模型预测的时间跨度对耦合计算结果的影响。
1.1 网络模型
LSTM网络[30]是用于时序预测问题的经典网络模型，该网络的提出是为了解决长期依赖问题和传统循环神经网络（recurrent neural network, RNN）的梯度消失和梯度爆炸问题。LSTM网络利用细胞状态实现信息相对长时间的传输，已被广泛用于各领域的时序预测问题[31-32]。LSTM网络存在很多变体，Cho等[33]提出了门控循环单元(gated recurrent unit, GRU），它将细胞状态与隐藏状态合并，减少了门控开关数量，具有更少的模型参数。Shi等[20]提出了ConvLSTM网络，将LSTM网络中的全连接层替换为卷积层，以提取数据空间特征并进行时序预测，其结构如图1所示，本文采用ConvLSTM网络对气固两相流的物理场进行预测。ConvLSTM网络常用计算公式如下：
式中，ft表示遗忘门，用于决定上一时刻的细胞状态ct-1中哪些信息被遗忘；it表示输入门，与候选细胞状态
共同决定将哪些新信息加入该时刻细胞状态ct；ot表示输出门，用于决定细胞状态中哪些信息将被输出；ht代表隐藏状态，由ot与ct计算得出；
是sigmoid激活函数；W代表网络卷积权重参数；b代表偏置项；下角标f、i、c、o分别对应遗忘门、输入门、候选细胞状态、输出门；
代表哈达玛乘积。
图1 ConvLSTM网络结构（修改自文献[19-20]）
1.2 深度学习与OpenFOAM耦合流程
图2 深度学习与OpenFOAM耦合流程图
T—模拟物理时长;t—当前物理时长;tOF—OpenFOAM单次运行物理时长;tDNN—深度学习模型单次预测物理时长
在线阶段使用Shell脚本将深度学习模型预测过程与OpenFOAM计算过程相耦合。首先调用OpenFOAM开启计算，获得前10帧计算数据；然后暂时挂起OpenFOAM计算，利用C++程序读取物理场数据，生成输入张量，加载深度学习模型对物理场进行预测，得到第11~20帧数据，之后将第20帧预测结果写入对应时刻的数据文件，对于未使用深度学习模型预测的物理场，则复 制第10帧数据写入对应时刻的数据文件，之后重启OpenFOAM进行下一轮计算。如此在OpenFOAM计算和深度学习模型预测之间交替循环进行，直至达到预设的物理时长。另外，为了保证最终结果满足物理定律，当剩余物理时间小于或等于一个完整的耦合计算周期（tOF+tDNN）时，全部使用OpenFOAM进行计算。
2.1 模拟方法
气固两相流的数值模拟方法主要有直接数值模拟（direct numerical simulation, DNS）、离散颗粒模拟（discrete particle method, DPM）和双流体模拟（two-fluid model, TFM）[34]。DNS通过细网格对颗粒进行几何解析来分析流场，气体对颗粒的作用力通过在颗粒表面积分获得，控制方程中不需要引入经验关联式，计算精度很高，但计算量较大。DPM方法在拉格朗日坐标系下跟踪每个颗粒的运动轨迹，气体被看作连续相，使用体积平均的NS方程来描述，计算量小于DNS。TFM将气体和颗粒看作可以相互贯穿的连续相，对每一相建立NS方程，计算量在三种方法中最小。对于工业级别设备，通常采用TFM方法进行模拟。本研究考虑到方法的扩展性及其在工业反应器中的应用前景，选择使用TFM模拟框架来构造数据集并考察深度学习的加速效果。首先，给出双流体模拟的控制方程[35]如下：
式中，εs、εg分别代表固相和气相的体积分数；τs、τg分别代表固相和气相剪切应力；Ms、Mg分别代表固相和气相所受的相间作用力；ps代表固相压力；µs代表固相黏度；λs代表固相体积黏度。固相应力采用颗粒动理论经验关联式封闭，更多模型细节见文献[36]。
2.2 数据集构造
本研究的模拟体系是一个高1.5 m、宽0.28 m的二维鼓泡床体系，如图3所示，静态床层高度为0.4 m，床层中颗粒体积分数为0.6，其他体系参数见表1。应用OpenFOAM中的multiphaseEulerFoam求解器对该体系进行模拟计算，并将模拟结果与文献[37]结果进行对比，图4展示了二者床层空隙率的径向分布对比结果，结果表明，本研究的OpenFOAM模拟结果与文献[37]给出的Fluent模拟结果基本一致，关键物理量的曲线走势与实验结果相近。
图3 二维床示意图
表1 模拟体系主要参数
图4 气体速度为0.46 m/s时床层的时均空隙率径向分布对比
为了构造足量的数据集，共设计了202个工况，在不同工况内改变体系中Ug、ρg、μg、ρs四个变量的取值，变动范围如下：Ug∈[0.26 m/s，0.78 m/s], ρg∈[1.225 kg/m3，2.45 kg/m3]，μg∈[18.4 μPa·s，36.8 μPa·s]，ρs∈[1650 kg/m3，2500 kg/m3]。构造数据集时，每个工况运行10 s，以0.01s作为帧间隔记录一次数据，每个工况生成1000帧数据。每个样本由20帧数据构成，前10帧是输入数据，后10帧是预测目标，样本间无重叠，故每个工况包含50个样本，所有工况共生成10100个样本。
气固两相流的模拟计算涉及众多变量，但是体系中部分变量可以由其他变量计算得到或者随时间变化较小，所以不必使用神经网络对所有变量进行预测。为了减轻模型构建的计算量，本文选取三个物理量进行模型预测，首先选取εs、Ug、Us考察模型预测效果，然后选取εs、Ug、p作为对比，考察预测不同物理量对CFD加速效果的影响。当使用神经网络预测εs、Ug、Us三个物理量时，将εs、Ugx、Ugy、Usx、Usy作为神经网络的5个通道，每帧数据的大小为5×300×56，构建数据集时对数据进行Min-Max归一化处理。训练集、验证集、测试集以8:1:1的比例进行划分，验证集和测试集分别根据速度分层采样，并保证尽可能多的工况被覆盖。
2.3 模型训练与结果分析
参照文献[19]构建神经网络，学习率为10-4，批大小为4，其他参数见表2。选用均方误差（MSE）作为损失函数，计算公式如式（13）所示。为了防止过拟合，当验证集损失经过20轮迭代后达到阈值时，停止训练。模型训练使用的节点包含4个Intel® Xeon Gold 5218 CPU (2.30 GHz)和1个NVIDIA® V100 GPU，单轮训练耗时约为18 min。本节以下分析以神经网络预测εs、Ug、Us为例。图5展示了模型损失函数随着训练轮次的变化，可以看出模型经过20轮以上迭代后已经基本收敛。表3给出了模型在测试集上的均方误差与峰值信噪比（PSNR）数据，据此可以对所训练的模型进行性能评估。PSNR是一种图像质量评价指标，计算公式如式（14）所示，数值越大表明预测图像与真实图像越接近，通常认为PSNR>30时图像质量较好[38-40]，所以此模型的训练效果是可以接受的。
表2 神经网络参数
图5 损失函数曲线
式中，N代表样本数目；Yi代表第i个样本的目标值；
代表第i个样本的预测值；MAX表示图像颜色的最大数值，本研究中值为1。
对于神经网络的预测效果，本研究分别从图像和数值两方面进行分析。图6给出了测试集中某算例的神经网络输入、最后一帧输出以及对应预测目标中εs物理场的图像，对比图6（b）、（c）可以看出，预测值与目标值非常接近。图7对比了测试集中某批次下神经网络预测值与目标值的频次分布直方图，结果显示，预测值与目标值的数据分布整体一致，部分区间略有差别。该结果表明，在常用的判定方式下，目前训练的网络模型具有良好的预测效果。
图6 神经网络中εs物理场的输入(a)、预测的最后一帧(b)和目标(c)图像
图7 神经网络预测值与目标值数据分布
2.4 模型部署与结果分析
根据图2所示流程，将深度学习模型预测过程与OpenFOAM模拟过程相耦合，耦合后的计算流程如图8所示。初始流场在OpenFOAM中计算500个时间步长，即0.1 s，接着利用深度学习模型预测0.1 s，之后再调用OpenFOAM计算0.1 s，此过程不断循环往复。该流程中，深度学习模型的一步预测替代了OpenFOAM 500个时间步的计算，节省了计算时长，而后续OpenFOAM计算又对计算结果起到了校正作用，保证了模拟结果满足守恒定律。在测试算例中计算10 s物理时间，然后统计εs、Ug、Us物理量的轴向分布，并与单纯使用OpenFOAM模拟结果相对比。由于速度大小主要体现在y轴分量上，因此对于Ug、Us只统计其y轴分量进行对比。
图8 耦合计算流程
（1）深度学习模型预测εs、Ug、Us
首先对εs、Ug、Us进行预测，εg根据εs计算得出。研究发现，计算过程中深度学习模型预测造成了颗粒总量的波动，在不同测试算例下，网格平均颗粒体积分数变化如图9所示。尽管每次模型预测引起的误差仅为10-3量级，但是多次预测的误差累积，使得10 s物理时间内网格平均颗粒体积分数的最大相对误差达到了12%以上。为了防止误差对最终结果产生较大影响，对模型预测结果中的εs进行校正，为了避免引起网格内原有数值的显著波动，将εs的偏差值平均分配到大于平均固含率的部分网格内，这样既保证了εs守恒，又不会导致空域中出现过高的εs值。
图9 网格平均颗粒体积分数变化
经过固含率校正后，图10（a）~（c）对比了不同物理量在耦合计算与OpenFOAM模拟中的轴向分布，εs、Ugy与模拟结果总体吻合良好，但是Usy在上方空域误差较大。究其原因，是Usy偏差大的区域对应的εs数量级极小（例如10-151），此处几乎没有颗粒，而在利用深度学习模型进行预测时，无法得到如此小数量级的数字；同时由于εs与Us的耦合作用，模型对εs的预测偏差导致了继续调用OpenFOAM计算时Us的计算偏差。由于该误差主要存在于床层界面之上，在实际应用中意义不大，所以可通过对预测结果进行识别过滤予以消除。观察数据发现，床层界面之上的颗粒体积分数几乎均在10-3之下，而床层界面之下的颗粒体积分数几乎均在10-3之上，因此将网格内εs小于10-3的数值全部替换为10-20，从而确保床层空域内的固含率足够小，同时又不会对模拟结果产生显著影响。此外，由于使用有限体积法进行数值求解时，网格内的εs数值还与邻近网格值密切相关，因此当某网格的邻近网格εs均小于10-3时，其εs值也替换为10-20。图10（d）~（f）展示了过滤后的结果，此时Usy在上方空域的计算结果与目标结果更加接近，而εs、Ugy的分布几乎不变，说明目前的过滤方法是可行的。
（2）深度学习模型预测εs、Ug、p
对εs、Ug、p进行预测，εg仍根据εs计算得出，各物理量的轴向分布如图11所示，结果显示，物理量的总体分布依旧与OpenFOAM模拟结果相近，但是效果不及图10（d）~（f）。以下从各物理量参与计算的过程分析该问题产生的原因。OpenFOAM中双流体模型的求解过程分为四个步骤：①采用上一时刻的p作为此刻p的初值；②根据简化的动量方程得到两相的速度，在该步计算中，速度对时间导数的离散格式是隐式Euler格式，深度学习模型预测的Ug、Us在该步参与计算；③使用PIMPLE算法进行压力修正，修正方程由动量方程与连续性方程得出，模型预测的Ug、Us在该步继续参与计算，该步收敛后得到满足界面通量连续性的p；④更新湍流模型，之后循环该过程直至速度收敛。由此得出，模型预测的p仅作为p迭代初值参与计算，而模型预测的Ug和Us在每个分步求解中均参与了计算。对于εs、Ug、Us的预测组合，p值由于是拷贝0.1 s之前的数据，所以其更新迟滞了0.1 s；同理，对于εs、Ug、p的预测组合，Us的更新也迟滞了0.1 s。从计算结果来看，在深度学习模型预测εs、Ug的相同条件下，增加预测Us比增加预测p在耦合计算中的效果更好。以上分析表明，当选取不同物理量进行模型预测时，需要考虑不同变量在求解器中的调用顺序与实际使用流程。
图11 深度学习模型预测εs、Ug、p时典型物理量的轴向分布对比
根据以上结果可知，深度学习模型预测εs、Ug、Us时，耦合计算流程与单纯OpenFOAM模拟结果吻合较好。为进一步提升整体的计算速度，本研究增大深度学习模型的预测跨度，将模型单次预测结果作为模型输入进行连续预测，以实现更好的加速效果。图12展示了耦合计算流程中，深度学习模型预测跨越不同物理时长时εs的轴向分布对比。其中单轮计算中，OpenFOAM均运行0.1s，深度学习模型预测分别跨越0.3、0.5、1.0s。结果显示，当深度学习模型预测与OpenFOAM以3:1或者5:1的比例运行时，结果偏差较小；以10:1的比例运行时，结果出现较大偏差，如果继续增大该比例，预测偏差可能会超出能够接受的范围。
图12 不同深度学习模型预测跨度下εs轴向分布对比
图13给出了不同深度学习模型预测跨度下耦合流程的计算时长与加速比。随着模型预测跨度的增大，计算总耗时逐渐下降且下降速率变缓，图中计算总耗时等于OpenFOAM计算耗时与DNN计算耗时之和。在图示范围内，DNN计算耗时与OpenFOAM计算耗时相比可以忽略不计，所以计算总耗时主要由OpenFOAM计算耗时所决定，而OpenFOAM计算耗时与其所承担的计算物理时间密切相关，二者变化趋势基本一致，均呈非线性变化趋势。在当前考察的1:1，3:1，5:1和10:1共4个跨度比下，耦合计算流程加速比分别为1.8、3.3、4.6和9.1。
图13 不同深度学习模型预测跨度下耦合流程的计算耗时与加速比
本文基于LibTorch与开源ConvLSTM网络，在二维气固鼓泡床中实现了深度学习与OpenFOAM耦合计算流程。该计算流程中，深度学习模型预测与OpenFOAM模拟计算交替进行，既减少了计算时长，又保证了计算结果满足物理定律。在控制体系内颗粒守恒后，εs、Ugy的轴向分布与单纯OpenFOAM模拟结果相近，Usy则由于深度学习模型对εs的极小数值预测不准确的原因在上方空域存在误差，为此提出了一种极小量修正方法消除该误差。此外，对比了深度学习模型预测不同物理量对结果的影响，发现预测εs、Ug、Us比预测εs、Ug、p在耦合计算中的效果稍好，分析后认为，深度学习模型预测对物理量的选取需要根据其在求解器内的调用情况确定。进一步地，考察了不同深度学习模型预测跨度比（DNN:OF）对加速效果的影响。当DNN:OF为5:1时，耦合计算能够加速4.6倍左右，εs轴向分布的误差范围较小；而当DNN:OF为10:1时，能够加速9.1倍，但是预测误差较大；耦合计算流程的总耗时主要由OpenFOAM的计算耗时决定。对于进一步减少整个模拟计算时长以及改善预测结果的工作，可以从以下两方面入手：（1）改进网络模型，从而增加深度学习模型预测的时间跨度；（2）将深度学习模型预测过程进一步耦合到OpenFOAM求解器中，以减少数据文件的读写时间。
转载自：《化工学报》：”一种耦合CFD与深度学习的气固快速模拟方法“，作者来自中国科学院过程工程研究所
