首页/文章/ 详情

JCP:涡轮泵中气蚀现象的建模与仿真

2月前浏览2086
 

   

   

   

摘要:本文提出了一种模拟包括气蚀和凝结的两相流模型,重点是其在涡轮泵流动中的应用。该模型通过扩展Kapila的方法,将热和自由能弛豫项结合到旋转参考系中的可压缩多相流方程中。所采用的数值方法基于Godunov格式的有限体积法,确保压力平衡和总能量通过弛豫步骤守恒。该模型使用了分别用于液相和气相的状态方程(EOS),实验饱和曲线以psat(T)表示,这是气蚀建模唯一需要的参数。为了扩大其适用范围到轴向泵之外,开发了一种特殊的黎曼求解器,以耦合旋转和非旋转区域内的流动域。热力学相变过程通过简单案例得到验证,并进一步用于模拟涡轮泵感应器中的水流。通过多种质量流量条件下的非气蚀条件获得了泵的特性曲线,并且在严重气蚀条件下成功验证了该模型,展示了由气泡破裂引起的性能下降。建模还提供了详细描述,包括气泡和其对叶片压力影响的分析。

SIMPOP
1. 引言    

在航天发动机中,涡轮泵用于在将推进剂注入燃烧室之前对其加压。为了减轻上方储罐的重量,涡轮泵必须在低进口压力和高转速下工作。在这种条件下,常常观察到气蚀片在轴向诱导器阶段发展。如果不加控制,这些两相流结构会导致泵的致命性能故障。为了更好地理解这一现象并最终改进诱导器的设计,数值研究是一种宝贵的分析工具。不幸的是,这种流动的建模仍然是大多数数值模型的主要技术挑战。一方面,均质模型带有额外的相连续性方程,包括由Rayleigh-Plesset型闭合规律驱动的质量传输源项,这些模型对经验相变参数非常敏感,并且不总是能够根据实验数据进行校准。另一方面,虽然采用混合物等温状态方程封闭的模型在上述问题上不那么突出,但对低温敏感流体(如液氢或液氧)的使用并未允许气蚀延迟现象的描述。事实上,导致蒸气压下降的温度变化无法通过不包括能量方程的模型来捕捉。

为了克服这些问题,提出了一种两相可压缩流动方法。关注点在于单速度和压力不平衡模型及其在气蚀流动中的验证和一致性属性,特别是在超速水下火箭、液滴冲击诱导气蚀、气泡动力学、液滴雾化以及其他物理问题(如炸药流、表面张力驱动流等)中的应用。在这种框架下,每个相位由其自身的状态方程(EOS)和在气蚀条件下发生的相变描述过程,通过热力学平衡(压力、温度和吉布斯自由能的平衡)进行参数化。为了适应数值泵计算的特定要求,采用了一个移动参考框架。在这种策略中,假设叶片保持固定网格,通过导出叶片旋转参考框架中的流动模型来考虑叶片运动。这种建模整合在ECOGEN中,一个用于多相、可压缩、多物理流动的开源工具,以验证当前方法进行涡轮泵诱导器的完整性能分析的能力。 

本文的组织如下:第二部分描述了感兴趣的两相流模型及其为了考虑泵叶片运动所需的导出。第三部分展示了解决两相流模型的数值方法,特别是水动力求解器,在旋转参考系中解决的数值特性,以及通过热力学弛豫过程描述相变。最后,第四部分对水供给涡轮泵诱导器在非气蚀和气蚀条件下进行了研究。提供了实验数据的压降比较,并讨论了气蚀空泡对叶片负荷的影响。

2. 控制方程    

2.1 两相流模型

本文提出的建模基于一种两相可压缩流动方法,其中每个相具有其自身的热力学量,但共享相同的速度。流动模型源自Kapila等人提出的模型,可能的相约束以弛豫项的形式表示在控制方程中。这些项,如在数值建模部分(第3节)中所描述的,对于气蚀建模是有意义的。在这种情况下,每个相由其自身的体积分数、质量和内能演化方程控制。

为了清楚起见,系统给出如下两相:

(1a)

(1b)

1c)

(1d)

(1e)

1f)

这些项表示为弛豫项:

(2)

然而,在无限弛豫情况下,这些项以及弛豫项在数值方案中不需要(见第3节)。

需要注意的是,模型可以在没有热和质量传输的情况下使用。在这种特定情况下,弛豫项将负责液-气混合物中的蒸气出现。这一力学效应已在若干研究中证明模拟了刚性弛豫系数和有限值的气蚀过程。

当进一步考虑热和质量传输时,模型将能够使蒸气在达到饱和时出现,然后通过气泡生长模型包括力、热和化学效应。

需要注意的是,该模型也可以在以相同状态方程(EOS)表示的相的混合物相态理论方程中找到。两者是等效的,使用一种形式的主要动机而不是与离散方程相关的数值方法,是为了确保混合物总能量的守恒(有关更多细节,见第3.1节)。

在这种建模方法中,每个相的热力学闭合由其自身的状态方程(EOS)保证。液相在这里通过Stiffened Gas方程表示。对于某一个相,其表示如下:

(3)

蒸气相服从理想气体规律的特定情况可以描述为Stiffened Gas EOS。在相变数值过程中,温度同样通过状态方程给出:

4)

由于流动在接近相变时的行为,采用实验数据拟合理论饱和曲线以确定相的EOS参数。当各相在饱和时达到完全热力学平衡状态时,这些情况使我们能够表示理论饱和温度:

(5)

其中:


2.2 泵叶片运动

为了描述涡轮泵中的流动,考虑叶片的运动是至关重要的。为此,建议不再在绝对参考系中描述流动,而是在以叶片速度旋转的参考系中描述流动。在这种描述中,被称为移动参考系(MRF)方法的流动模型必须被推导出来。即使这种方法在地球物理系统中被广泛使用,其对两相可压缩流动的推导并不是那么简单。特别是对于欧拉单相流模型,总能量方程在不可压缩流动的情况下没有给出,或者如文献中所见,其形式存在较大争议。因此,特别关注在旋转参考系中建立可压缩两相流动的方程。

为了在旋转参考系中获得流动模型,主要有两种方法。前者基于拉格朗日描述,后者则应用了连续的伽利略和旋转变换。在以下内容中,由于其严格定义了不同时间和空间算子,所以选择后者。 

为了便于推导,我们将绝对参考系记为R,移动参考系记为R'(见图1)。

图1 绝对参考系与移动参考系

考虑到在我们的应用案例中参考系之间的变换是纯旋转,位置矢量可以写成:

从位置矢量的定义,现在可以定义不同参考系中速度之间的联系:

由于公式6,上述关系可以重写为:

(7)

考虑到M是反对称矩阵,最终移动参考系中的速度可以简化为:

在下面的推导中,需要注意任何标量场$b$不会受到参考系变化的影响。因此,对于任何标量量,引号将被省略。此外,已定义参考系变化后,提醒位置矢量和时间是独立变量。


  • 标量方程

为了推导,考虑任何守恒方程可以用物质导数表示:

(8)

从公式(8)可以将物质导数的时间部分转换为移动参考系的计算:

公式(8)的对流部分可以展开为:

因此,移动参考系中标量量b的物质导数为:

(9)

为了展开公式(9)中的最后一项,我们注意到:

使用这些性质,可以建立标量量$b$的保守传输方程在移动参考系中具有相同形式:

(10)

这种转换可以直接应用于控制方程(1a)、(1b)、(1c)、(1e)和(1f)。


  • 混合动量方程

混合动量方程的转换需要注意进一步的计算。为推导目的,混合动量方程(1d)用物质导数表示为:

(11)

压力梯度转换如下:

速度物质导数分两步展开。首先,将物质导数算子转换为移动参考系:

其次,绝对速度被替换:

使用压力梯度和速度物质导数,混合动量方程在移动参考系中通过将其乘以旋转矩阵$O$得到:

可以写成保守形式:

(12)


  • 混合总能量

由于数值原因,混合总能量方程在叶片参考系中的形式被给出。由于混合总能量包括混合动能和混合内能,建议通过将这些分量相加来构建它们。动能可以通过动量方程的标量乘积与相对速度获得:

动量方程的积分与相对速度得到的总能量方程:

(13)


  • 在移动参考系中的控制方程

现在,整个控制方程系统在移动参考系中使用(10)和(12)写出。与原始流动模型(1)相比,除了速度定义外,在动量和能量方程右侧项的额外项是明显的。

注意速度定义后,系统可以以一般矢量形式表达(去掉引号以简化符号):

(14)

(15)

注意,如果在绝对(固定)参考系中考虑流动,S矢量消失以恢复原始系统(1)。在这种情况下,所考虑的速度矢量简化为绝对速度。此性质对于在第3.2节中描述的移动和静止参考系中耦合流动解非常有用。

3. 数值建模    

在移动参考系中的两相流模型需要包括相变的贡献,考虑以下形式:

(16)

其中不同矢量的定义早先在公式(15)中给出。

模型(14)使用ECOGEN开源CFD求解器求解,采用显式有限体积法和算子分裂策略。在时间步长上状态向量的演变由一系列算子执行:

(17)

其中每个算子独立地从右向左应用。

每个算子解决流动模型的特定部分:双曲算子使用Godunov类方案解决流动的水动力学部分;源算子执行源项的积分,在这种特定情况下,包括与叶片运动相关的附加项的积分;弛豫算子负责描述在气蚀现象的相变期间发生的热力学演化。

每个算子的数值方法将在以下章节中具体介绍。


3.1 双曲系统的求解

在此步骤中,方程组(14)简化为:

(18)

模型(18)通过类似Godunov的显式有限体积法方案求解。体积为V的单元i的状态向量U的时间演变,通过以下公式获得:

(19)

对于体积分数方程,速度被接触不连续性的速度替代,以确保正确处理此传输方程:

到目前为止,数值方案是相当保守的,基于求解保守系统,并添加非保守项的贡献。不幸的是,目前对相内能量方程中存在的非保守项的近似不能保证混合总能量的守恒以及相间能量的正确分配。因此,为了解决这个问题,使用Schmidmayer等人的方法校正相内能量通量。


3.2 旋转项

1)源项积分方案

在移动参考系方法中,叶片的运动通过虚构的科里奥利力和离心力来考虑。从数值角度来看,这些项被视为源项,例如其他人所做的(例如Biedron等人或Prachar)。按照分裂算子的过程,方程(14)简化为一阶常微分方程系统:

(20)

初始解由双曲算子提供。在ECOGEN中,时间积分通过一阶或更高阶的方法执行,例如Runge-Kutta 2-4阶方法。实际上,为了与水动力求解器保持一致,使用了一阶Euler方案进行计算。

2)旋转和静态区域之间的相互作用

虽然移动参考系方法允许在旋转参考系中描述流动,即泵叶片的参考系,但这种描述仅在叶片附近区域有效。为了在上下游区域中描述流动,建议利用静态参考系和旋转参考系中流动模型的差异仅涉及引入的额外源项和使用的流动速度。因此,一旦在网格中定义了静态和移动区域(见图2),建议仅将源项应用于需要的区域,并通过在黎曼问题(RP)期间分隔它们的界面来耦合静态/移动流动区域。

图2 轴流泵上静态和旋转流体区域的示例

图3 静态参考系中的单元和旋转参考系中单元之间的界面处的黎曼问题

考虑一个界面分隔了属于旋转区域的单元i和属于静态区域的单元j(见图3),建议为每个单元建立特定的流量:

(21)

RP的左右初始状态必须相应地准备好:旋转区域单元的流量在绝对参考系中进行转换。

现在单元i和j的状态处于同一参考系中,静态单元的流量用常规近似黎曼求解器计算,在这种情况下使用HLLC求解器。在此步骤中,静态参考系中黎曼问题解的原始变量矢量也被提取。


3.3 相变

在这种两相流模型中,导致液体/蒸气混合物的热量和质量传递被描述为回到热力学平衡的过程,即:

(22a)

(22b)

(22c)

这里,认为热力学平衡通过无限大的弛豫速率瞬间达到。对于其他替代方案,读者可参考弛豫系数的确定以及Schmidmayer等人关于有限压力弛豫的初步研究。

在下面的内容中,解释了用于描述气蚀过程中相变的瞬间弛豫过程及其在达到流动条件时的触发。

1)通过热力学弛豫的相变

热力学弛豫对应于分裂过程中的最后一个算子。在此步骤中,方程(14)简化为常微分方程(ODE):

(23)

为了求解由公式(23)表示的瞬时弛豫过程,可以认为所有导致热力学平衡的弛豫过程同时发生,也可以认为每个弛豫过程(压力、温度、化学势)依次应用。虽然后一种方法通常减少了找到最终热力学平衡状态所需的计算时间,但在模拟涡轮泵内的流动时,难以预先确定各相之间平衡的顺序。因此,我们建议使用基于Le Métayer等人工作的同时弛豫方法。

图4 压力、温度和吉布斯自由能弛豫过程的状态

这种方法包括通过代数方程(24)求解液体/蒸气混合物在饱和状态下的热力学状态。寻找从初始非平衡状态“0”(图4)松弛到的平衡热力学状态“f”基于质量守恒(24b)和能量守恒(24a)在发生相变的计算单元的体积上;并基于各相之间的热力学平衡约束(公式(24c)到(24e))。最后,为了闭合系统(24),添加体积分数约束公式(24f)。

(24a)

(24b)

(24c)

(24d)

(24e)

(24f)

各相的密度和能量可以使用其各自的状态方程表示。根据平衡约束公式(24c)和(24d),它们变为:

(25)

此外,公式(24f)和(24b)得出:

(26)

因此,公式(24a)变为:

(27)

注意,最后一个方程仅涉及弛豫结束时的压力和温度。引入最后约束公式(24e)中的状态方程允许我们获得在饱和时温度和压力之间的关系:

(28)

当饱和曲线引入到公式(27)中时,后者仅成为最终压力的函数:

(29)

最后,可以使用迭代方法如Newton-Raphson方法在公式(29)上确定最终状态的平衡压力。完整过程的示意图见图5。

图5 压力、温度和化学势松弛的求解过程示意图

一旦找到压力,其余各相的热力学量可以用公式(28)、(25)和(26)轻松更新。

2)无阈值条件下的热量和质量传递触发

由于相变仅在达到饱和条件时才发生,现在有必要建立一个标准来评估是否应应用热力学弛豫。在未达到标准的情况下,必须找到纯相形式是否存在。通常,多相流方法中很少处理此问题,这些方法中经常使用液/气相体积分数的阈值。不过,这些阈值的使用并不基于物理现象,纯数学的考虑可能导致蒸汽生产的显著不同。本文提出的方法包括验证在所考虑的计算单元内,给定的质量和能量下,是否存在纯相态的稳定形式。

实际上,初始状态中包含的质量和内能分别赋给每个液相和气相:

(30)

然后根据状态方程(3)和(4)计算相的压力和温度。此外,使用公式(5)评估该相的饱和温度。

如果相k是液体并且液相温度低于沸点,则未达到饱和条件,稳定解是过冷液体(见图6)。

图6 (P,T)图中的相变

同样地,如果相k是气体并且气体温度高于沸点,则未达到饱和条件,稳定解是过热气体。如果不存在上述稳定条件,则混合状态对应于处于热力学平衡的液/气混合物,需要应用上一节的热力学弛豫过程。完整过程在图7中总结。

图7 松弛过程评估图。只有在没有找到纯液体或纯蒸汽形式的稳定溶液的情况下,才会进行热力学弛豫 

3)相变触发的验证

为了确保无阈值条件下的相变触发正常工作,建议研究在热绝缘盒中静止流体的热力学状态,向其中注入或抽取能量。

首先,定义纯相的热力学状态。然后,将热源项应用于控制体:当为正时,向系统提供能量;当为负时,系统失去能量。由于域的体积是恒定的,流体遵循等容过程,这可能导致在达到饱和条件时发生相变。

为此测试,状态方程的参数在宽温度范围内校准,如Saurel等人所做的那样,这些参数见表1。将这种校准与实验数据的比较也显示在图8的饱和曲线上。

表1 水的状态方程参数

图8实验和理论(p,v)饱和曲线(左)和放大感兴趣的相变过程(右)

为了验证,建议研究从中移除能量的纯蒸汽组成的流体的液化过程(见图9)。一旦获得液/气混合物(状态3),进行反向操作,包括加入能量使液体蒸发,直到流体返回状态1。

图9 冷凝/汽化相变测试说明

初始纯蒸汽状态的压力为p = 1.5 bar,温度为T = 450 K(见图8初始状态)。在0.1秒内通过热汇以Q= -10^6 W m^-3冷却体积,然后停止模拟,并用热源Q = 10^6 W m^-3重新启动0.1秒。注意,这些值经过精心选择,以确保穿过饱和曲线并发生相变,如以下结果所示。图10显示了相质量分数随时间的演变以及混合物的压力。特别地,可以区分以下阶段:能量损失导致体积中的压力和温度下降,使饱和在状态2达到,此时流体压力等于饱和压力;能量减少促进液体形成,如状态3中质量分数增加所示;向系统引入能量,蒸汽质量分数增加,直到穿过露点,压力保持在饱和压力,见状态3到状态2的转变。当穿过饱和曲线时蒸发完成,蒸汽质量分数再次为1。

图10 相变过程中流量随时间变化的质量分数(左)和压力(右)

考虑到相变在达到液化/蒸发过程的饱和条件时正确触发,并且对于蒸发/液化过程获得了相似的结果,接下来采用相变触发方法。

4. 结果    

在航天工业中,大多数航天器使用低温推进剂,因为其推力性能。尽管本文提出的模型能够描述此类流体,但以下研究的是一个水供给的涡轮泵诱导器,以便与现有实验数据进行比较。

图11 给水诱导器的几何形状

所关注的诱导器是一个带有可变截面轮毂的三叶片诱导器(见图11)。在实验中,将常温液态水以恒定质量流量注入以5000转/分钟旋转的诱导器中。对于给定的入口压力,测量在叶片上游和下游区域之间产生的超压。最后,重复这一过程,以不同的质量流量和入口压力来表征诱导器在广泛的操作条件下的性能。

在数值模拟中,为了重现实验装置,入口处施加质量流量和温度T= 293 K,而出口处的压力p_out设置为考虑的实验测试点,如图12所示。使用基于特征的方法数值求解这些边界条件,其中考虑了特定的黎曼问题以确定边界通量。值得注意的是,边界条件故意应用在远离叶片区域的位置,以避免与旋转流的不良数值交互。为了考虑叶片的运动,在叶片区域应用移动参考系(MRF)源项,使用第3.2节中介绍的方法。

图12 给水诱导轮的数值设置。实验压力探头分别位于距离上游和下游探头诱导轮鼻毂0.4Db和1.2Db的位置(叶片直径为λ)

图12中还显示了用于测量叶片产生的超压的压力传感器(请注意使用的是实验装置的精确位置)。

在本文提出的两相流方法中,每个相必须通过其自己的状态方程进行描述。在这种情况下,液态和气态水相分别由刚性气体和理想气体状态方程描述。其参数的校准是通过对实验饱和曲线应用最小二乘拟合方法得到的。在常温T = 293 K下校准的参数如表2所示。

表2 液态/蒸汽水状态参数方程

本节的第一部分致力于无气蚀条件下的流动模拟。目标是评估网格收敛性和CPU时间消耗方面的最佳网格折衷,同时获得对流场(压力场、叶片负载)的信息,这将有助于后续与气蚀条件下模拟的比较。同时,整体诱导器性能也将与实验数据进行比较,以证明数值设置适合描述无气蚀流动。本节的第二部分涉及在两种操作质量流量下气蚀状态下的流动模拟。将显示数值模型能够描述由于气蚀空泡形成导致的性能下降,并简要分析其对叶片负载的影响。同样将给出与实验数据的比较并进行讨论。


4.1 非气蚀流动

1)网格收敛性

为了评估数值解对空间离散化的依赖性,考虑了三种四面体网格:一个约919k单元的网格,一个中间网格约1527k单元,一个稍微更精细的网格约2543k单元。所有网格均使用Pointwise软件(V18.4R4)生成,并使用Gmsh开源软件中提供的METIS算法进行分区。根据几何区域分布的单元数量见图13。从图中可以看出,特别注意了叶片区域的网格划分。为了确保较好的网格质量,检查了诸如等体积偏斜度和最大包含角等指标,分别低于0.8和135°。

图13 根据不同诱导轮网格的几何区域,元素数量的分布

为了评估网格收敛性,在叶片上下游测量压力,并与相应的实验数据进行比较(见图12)。

图14 网格收敛:在标称质量流量下,压力的相对误差是网格尺寸的函数

结果以相对误差的形式总结在图14中。从图中可以清楚地看出,网格细化有助于上游压力向实验压力收敛。叶片下游,由于出口压力被设定,因此无论网格大小都能得到良好的压力估计。剩余的细微差异是由于在出口缺少可用的实验数据,出口压力设定为在靠近外壳的下游位置测量的实验压力。

由于气蚀触发由上游压力驱动,因此应选择网格以最小化其误差。尽管2543k单元网格完全满足这一标准,但计算成本显著更高,如表3所示。因此,为了在模拟成本和压力估计之间实现合理的折衷,选择了1527k单元的网格大小。

表3 比较两个网格的模拟时间和消耗的资源数量

2)特性曲线 

在没有气蚀的情况下,诱导器或泵的性能通过测量在给定流量下获得的超压来确定。使用压力表在叶片的上下游区域测量外壳上的超压(见图12)。超压以无量纲形式表示为:

(31)

为了评估诱导器的特性曲线,将出口压力设置为给定质量流量的实验操作点的值,并对不同流量重复这一过程。

图15中提出了与实验数据的比较。从左图可以看出,除最低流量外,超压的一致性良好。在右图中可以看出上游压力估计不佳。这可能是由于叶片尖端和外壳之间间隙区域的数值分辨率不足。确实,在低流量时超压较高,可能会在间隙区域发生回流,导致实验中的上游压力较低,而在这些模拟中由于该区域的网格粗糙而未观察到。还要注意,数值下游压力略高于实验压力,原因与上一节提到的相同。

图15 诱导轮特性曲线:作为流量函数的超压(左)和无量纲压力(右)

3)压力场分析

虽然没有达到气蚀状态,但叶片周围的压力场允许确定可能发生气蚀的低压区域。图16给出了研究的四种质量流量的压力场的叶片间视图。图中显示,在低流量时,压力下降主要位于叶片的前缘,位于凸面侧,且流量越低,压力下降越大。然而,在最高流量时,压力下降发生在叶片通道中:进入流动已经很快,进入叶片通道后加速,导致压力下降不再被叶片负载所补偿。

图16 不同半径下叶片间压力场情况

4)叶片负载

为了进一步表征诱导器的性能,我们建议通过沿叶片上下表面测量静压来评估叶片负载,如图17所示。为了了解压力如何随叶片跨度变化,这些测量在两个半径处进行。

图17 静压传感器在下表面和上表面上的位置,用于测量叶片的工作

图18 不同质量流量下叶片上下表面的静压分布

图18展示了静压分布。可以看出,在名义流量下,上下表面之间的超压沿整个叶片长度发生。从外壳和中脉的压力测量中,还可以注意到,从弦长的50%开始,跨度对超压有所贡献。在最大流量下,可以观察到,超压主要沿叶片长度产生,跨度方向的超压较小。需要注意的是,从前缘开始的弦长20%以上,上下表面之间的超压区域为负:叶片没有将能量传递给流体,而是流体推动叶片。这解释了为什么诱导器在高流量下的性能较低。


4.2 气蚀状态研究

本节的目的是确定气蚀发生时导致诱导器性能下降的流动条件。为了确定这种性能下降,在给定流量下,将无量纲超压评估为气蚀数的函数,定义为:

(32)

过对不同出口压力进行模拟构建性能曲线。必然地,随着出口压力降低,为确保诱导器超压在正常操作条件下保持恒定,叶片上游的压力也会降低。然而,如果上游压力降至蒸汽压力以下,则会发生气蚀并形成蒸汽。如果流动受到气蚀空泡的严重影响,诱导器将无法提供预期的超压,进而导致性能下降。为了研究这种诱导器的性能,考虑了两个在无气蚀状态下与实验数据良好一致的质量流量。

1)名义质量流量下的气蚀

为了概览气蚀的发展,提出了不同气蚀数下流动的定性分析。为了帮助完成这一任务,图19表示了从0.152到0.004范围内的气蚀数的稳态蒸汽体积分数。


图19 当𝜏=0.152(左上)、𝜏=0.051(右上)、𝜏=0.009(左下)和𝜏=0.004(右下)时,标称质量流量下的蒸汽袋𝛼𝜐>5%

对于最大的气蚀数0.152,没有蒸汽。在中间范围[0.051; 0.009],一些蒸汽空泡开始在叶片尖端间隙附近的前缘周围发展。注意这些空泡首先发展的区域与无气蚀状态下的观察结果一致(见图16)。对于<0.009,气蚀空泡达到叶片间通道,导致液体流动阻塞。这种分析可以与图20中描绘的性能曲线联系起来。

图20 标准质量流量下的性能曲线:超压(左)和无量纲压力(右)

与实验数据相比,当气蚀尚未发生或发生适度时,压力估计与实验数据很好地一致。然而,数值模型往往预测性能下降发生在稍低的气蚀数。这如同在无气蚀结果部分所讨论的,可能是由于当前在间隙区域的网格不允许捕捉促进更快气蚀发展的回流。由于模拟成本原因进一步的网格改进超出了本次研究的范围,并且这一现象在更高流量下应较不重要,因此下一部分的研究将在更高的质量流量下进行。

2)更高质量流量下的气蚀

在本节中,气蚀状态的研究在质量流量=1.15m_n下进行,在无气蚀状态下观察到上游压力估计较好。如前所述,评估该流量下的诱导器性能,并在图21中给出相应的性能曲线:实验和数值超压均与气蚀数对比,以及上下游压力(分别为左图和右图)。可以清楚地看到,在性能下降前压力估计良好,相对误差约为5%。如在名义质量流量下所见,数值超压下降仍发生在较低的气蚀数。然而,使用该质量流量导致的临界气蚀数相比名义质量流量结果提高了约20%。

图21 质量流量𝑚=1.15×m_n时的性能曲线:超压(左)和无量纲压力(右)

为了更好地理解气蚀空泡如何影响诱导器的超压性能,建议如第4.1.4节所述沿叶片上下表面测量静压。在图22中给出了轻微气蚀和遇到性能下降时的静压分布比较。

图22 质量流量𝑚=1.15m_n,𝜏=0.011(左)和𝜐=0.009(右)时叶片上下表面的静压分布

与无气蚀状态下的叶片负载相比,可以明显注意到在弦长55%左右的凸面侧压力突然下降,当发生轻微气蚀时。这种行为是蒸汽空泡存在的特征,如图23所确认,并且可以通过气蚀空泡内压力受液/气混合物的饱和约束解释。

图23 𝜏=0.011(左)和𝜏=0.009(右)时蒸汽体积分数的叶片间视图

对最低气蚀数=0.009的压力分布检查显示,整个上表面和下表面40%的弦长上压力等于蒸汽压力。因此,叶片传递给流体的能量仅在剩余的60%弦长上完成。然而,请注意,这种有效面积的减少仅导致整体超压下降2.5%。这是因为诱导器主要在叶片尖端的跨度方向上工作,并且上表面的压力下降增加了在后缘前提供给流体的工作量。

5. 结论    

作为使用两相流方法模拟涡轮泵中空化流动的第一步,Baer–Nunziato型的单速和压力不平衡模型已被调整以满足此类计算的要求。这种模型的选择是由于它考虑了流体的可压缩性、是双曲的、符合热力学第二定律,并且在建模许多气蚀流动方面有坚实的基础。

为了考虑叶片运动对流动的影响,模型在移动参考系中进行了推导。此外,提出了一种数值方案,允许使用几何的静态和移动区域。

在这类两相流模型的框架下,与气蚀相关的相变被描述为返回热力学平衡的弛豫过程,这是一个无参数的过程。按照这种无参数相变方法,过程的触发现在由纯相状态的搜索及其与饱和条件的兼容性驱动。该方法在一个学术测试案例中进行了测试,证明了这种方法能够仅基于热力学考虑来描述相变。

最后,当前的建模在一个水供给的透平泵诱导器研究中进行了测试。首先,在无气蚀操作条件下对诱导器性能进行了表征。发现扬程下降估计与实验数据非常一致,特别是在质量流量等于或高于名义质量流量的情况下。基于这些令人鼓舞的结果,进行了气蚀操作条件下的研究。模型显示能够描述由于气蚀空穴显著发展导致的性能下降。然而,观察到下降发生在低于实验结果的入口压力下,但在较高质量流量下获得了更好的匹配。

这可能表明在外壳和叶片尖端之间间隙区域的回流捕捉不佳,特别是在名义质量流量下,需要进一步努力改善该区域流动的分辨率。在这方面,未来的工作可以包括使用总变差减小技术将ECOGEN求解器扩展到非结构化网格的二阶。

提高该方法的计算性能也是扩展其适用性的一种方式。一个可能的未来发展是使用隐式数值方案,如Le Martelot等人所示。

尽管仍有许多重要方面需要改进,但这些结果非常令人鼓舞,展示了将当前两相建模应用于涡轮泵中气蚀流动研究的可行性。


翻译转载自:《Journal of Computational Physics》:“Modeling and simulation of the cavitation phenomenon in turbopumps” 

来源:多相流在线
多相流燃烧化学航天UM理论Pointwise控制
著作权归作者所有,欢迎分享,未经许可,不得转载
首次发布时间:2024-09-15
最近编辑:2月前
积鼎科技
联系我们13162025768
获赞 109粉丝 112文章 302课程 0
点赞
收藏
作者推荐

一种耦合“CFD+深度学习”的气固快速模拟方法

摘要:在计算流体力学领域,深度学习被用于重建流场、预测曳力、求解泊松方程、加速流体模拟等方面的研究。为了加速气固两相流的模拟计算,使用卷积长短时记忆网络对物理量进行预测,并基于LibTorch实现深度学习模型预测与OpenFOAM的耦合。通过与单纯OpenFOAM模拟结果对比,发现深度学习模型预测存在颗粒体积分数不守恒、极小数值预测不准确的问题,先后通过体积分数校正和网格数据过滤消除了前述影响。选取不同的三个物理量组合进行深度学习模型预测以加速CFD计算,在同样选取颗粒体积分数和气体速度的条件下,对比了增加预测颗粒速度和压力对耦合流程计算结果的影响,其中,增加预测颗粒速度的计算结果较为准确,经分析发现这可能与求解器对不同变量的调用顺序有关。SIMPOP引言 气固两相流广泛存在于化工过程中,为了理解反应器内的气固复杂流动反应过程,计算流体力学的应用日益增多。气固体系的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.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耦合计算的流程,该流程分为离线阶段和在线阶段。离线阶段包含三部分工作:数据集构造、模型训练与优化和模型转换。在数据集构造阶段,主要利用数据预处理Python程序对OpenFOAM模拟计算结果数据进行标准化处理,并将其划分为训练集、验证集和测试集;在模型训练与优化阶段,利用训练集训练模型、验证集调整模型超参数、测试集测试模型效果,若预测结果符合精度要求,则保存模型参数;在模型转换阶段,通过跟踪输入在模型中的运算将PyTorch模型转换为TorchScript,转换后的模型可在C++环境中直接调用。 图2 深度学习与OpenFOAM耦合流程图T—模拟物理时长;t—当前物理时长;tOF—OpenFOAM单次运行物理时长;tDNN—深度学习模型单次预测物理时长在线阶段使用Shell脚本将深度学习模型预测过程与OpenFOAM计算过程相耦合。首先调用OpenFOAM开启计算,获得前10帧计算数据;然后暂时挂起OpenFOAM计算,利用C++程序读取物理场数据,生成输入张量,加载深度学习模型对物理场进行预测,得到第11~20帧数据,之后将第20帧预测结果写入对应时刻的数据文件,对于未使用深度学习模型预测的物理场,则复 制第10帧数据写入对应时刻的数据文件,之后重启OpenFOAM进行下一轮计算。如此在OpenFOAM计算和深度学习模型预测之间交替循环进行,直至达到预设的物理时长。另外,为了保证最终结果满足物理定律,当剩余物理时间小于或等于一个完整的耦合计算周期(tOF+tDNN)时,全部使用OpenFOAM进行计算。2. 方法验证与实例分析 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&gt;30时图像质量较好[38-40],所以此模型的训练效果是可以接受的。表2 神经网络参数 图5 损失函数曲线表3 测试集上的模型性能 式中,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的分布几乎不变,说明目前的过滤方法是可行的。图10 深度学习模型预测εs、Ug、Us时固含率校正后的典型物理量轴向分布对比[(a)~(c)]及继续处理小量后的轴向分布对比[(d)~(f)] (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 不同深度学习模型预测跨度下耦合流程的计算耗时与加速比3. 结论 本文基于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与深度学习的气固快速模拟方法“,作者来自中国科学院过程工程研究所 来源:多相流在线

未登录
还没有评论
课程
培训
服务
行家
VIP会员 学习 福利任务 兑换礼品
下载APP
联系我们
帮助与反馈