改进的物理融合神经网络在瑞利-泰勒不稳定性问题中的应用1)

丘润荻 王静竹 黄仁芳 杜特专 王一伟, 黄晨光

* (中国科学院力学研究所流固耦合系统力学重点实验室,北京 100190)

† (中国科学院大学未来技术学院,北京 100049)

** (中国科学院大学工程科学学院,北京 100049)

†† (中国科学院合肥物质科学研究院,合肥 230031)

引言

两相流动演化问题普遍存在于海洋工程、生物医学、工业制造等领域[1-5].数值模拟方法是工程中分析两相流动问题的最重要手段,两相界面的捕捉精度直接决定计算与分析的准确程度.因此,提高相界面捕捉精度和降低计算误差一直是数值模拟两相流动演化过程的热点问题[6-9].相场法(phase-field method)由于具有物理含义明确、算法质量守恒和界面捕捉精度高等优势,被广泛用于求解黏性与表面张力主导的界面演化、大密度比水汽界面演化、界面强非线性演化、热毛细对流等问题[10-13].其中,瑞利-泰勒(Rayleigh-Taylor,RT)不稳定性包含非定常、强非线性等因素,对数值方法的计算规模和计算精度要求很高[14-17].经典的RT 不稳定性发生在密度分层流体中,当重密度流体处于轻密度流体之上时,在重力作用下,相界面上的小扰动会逐渐演化发展,重密度流体向下延伸形成一个尖峰形状的头部,下方的轻密度流体向上移动形成气泡的形状,最终出现不断扩展的湍流混合层.

针对基于神经网络的两相流场求解问题,最近引人关注的是通过物理融合神经网络(physicsinformed neural networks,PINNs) 求解偏微分方程[18-22],通过添加控制方程和边界条件作为神经网络的约束,大大提高流场建模和反演能力.Raissi等[23]基于PINNs 框架,使用已知的浓度分布信息推断未知的速度场,并将该方法成功应用于圆柱绕流和三维血管造影中.Cai 等[24]在PINNs 框架下,使用实验获得的温度场反推出了低雷诺数范围下流场的流动情况,并将该方法推广到传热问题的流动分析中.Buhendwad 等[25]探索了PINNs 在正问题和反问题中的应用,计算了包括Couette 流动、Poiseuille 流动、气泡变形、气泡震荡和气泡上升等两相流算例,验证了PINNs 求解两相流问题的可行性.

Qiu 等[26]使用相场法作为界面捕捉方法,结合纳维-斯托克斯(Navier-Stokes,N-S)方程作为控制方程,建立了基于相场法的物理融合神经网络(PINNs for phase-field method,PF-PINNs),并通过单涡剪切流和气泡上升算例,均获得了满意的效果,证实了PF-PINNs 求解大密度比两相流的可靠性.然而,在该架构中,PF-PINNs 直接求解相场法方程中最高阶导数项,速度仍有待提升.针对PINNs 的效率问题,Lü等[27]提出了深度混合残差方法(deep mixed residual method,MIM),可以加速神经网络的训练过程,该方法通过将一个高阶导数降为多个低阶导数来提升求解速度,该架构成功用于Poisson 方程、Monge-Ampere方程、Biharmonic 方程和Kortewegde Vries 方程的求解中.

为了提升两相流PINNs 的计算效率,本文在PF-PINNs 框架基础上,利用深度混合残差方法,加入化学能作为神经网络中的辅助变量,用于求解RT 不稳定性问题,来分析改进的PF-PINNs 的计算精度和计算效率.本文的具体内容包括:第一节中介绍了两相流问题的基本控制方程;第二节中介绍物理融合神经网络方法以及计算过程中采用的加速策略;第三节中将改进的PF-PINNs 用于求解RT 不稳定性算例,并分析了改进算法的加速效果;最后一节中对本文进行总结和展望.

1 基于相场模型的二维N-S 方程

两相流动中的控制方程包括连续性方程、动量方程与界面演化方程.本文的界面演化方程由基于Cahn-Hillard 模型的相场法描述,其考虑不可压缩不相溶两相流体组成的混合流体,并假设界面扩散速度和化学能的梯度成正比来精确描述界面附近的扩散效应,确保界面演化方程在满足质量守恒的同时能准确捕捉界面.关于相场法的推导和本文使用的相分数定义、范围等可参考Qiu 等[26]的工作.将相场法的界面演化方程与原始的N-S 方程结合,可导出不可压、不相溶两相流动的控制方程与表面张力的形式

其中,C为混合物的相分数,本文中相分数的范围定义在-1 到1 之间,即在重流体中C=1,在轻流体中C=-1,界面附近的区域内-1 ≤C≤1.u=(u,v) 为流场的速度矢量,p表示压强,为重力加速度矢量;M0为渗透率,φ为化学能,ε为界面宽度;为表面张力,σ0为表面张力系数.ρM和 μM为混合物的密度和黏度,其可分别表示为

其中,下标L表示重流体的物理量,下标G表示轻流体的物理量.

2 物理融合神经网络

物理融合神经网络主要用于构建流体动力学模型,即建立时空坐标与待求解物理量之间的映射关系.在原始PF-PINNs 中,该映射关系可表示为如下形式

其中,(x,y,t) 为时空坐标,(u,v,C,p) 为待求物理量,fNN为三个输入、四个输出的全连接神经网络,Θ 为全连接神经网络的权重和偏差.将式(1)～式(4)通过残差形式写入神经网络,可得到PF-PINNs 的物理约束部分.为了提升PF-PINNs 的计算效率,本文参考MIM[27]的思想,将化学能视为一个辅助变量,也加入到神经网络的输出当中,并修改了损失函数中与相分数方程有关的约束形式.本文所述普通形式神经网络与深度混合残差形式神经网络在界面演化方程和相应损失函数的区别如表1 所示.

表1 普通形式神经网络与深度混合残差形式神经网络在界面演化方程中的区别Table 1 Differences between a normal neural network and a deep mixed residual method network in interface evolution equation

由化学能的公式可知,化学能与物理问题之间并不存在直接关联;它由相场法的定义给出,且在两相流问题中仅会影响界面附近的区域.但化学能的表达式却引入了相分数的二阶导数,高阶导数的存在扩大了神经网络计算图的规模,使自动微分的计算成本上升.为了使读者大致了解MIM 策略生效的原因,本文以自动微分计算二阶导uxx为例进行说明.在大部分机器学习库中,自动微分通过函数ux=tf.gradient(u,x)实现,通过u计算ux需要调用一次计算图;而对uxx的计算则通过uxx=tf.gradient(ux,x)实现,由于调用ux的同时也需要调用ux的计算图,故通过ux计算uxx需要调用两次计算图,那么通过u计算uxx需要调用三次计算图.依此类推,通过u计算uxxx需要调用六次计算图,通过u计算uxxxx需要调用十次计算图,计算高阶导数的开销可通过降低导数阶数得到,若增加一个网络输出q,并将uxxxx改写为qxx,令q=uxx,则计算qxx共需要调用2 × (1+2)=6 次计算图.虽然将硬约束改为软约束会使结果损失一定的精度,但是计算所需的开销大幅下降,这是MIM 策略在PF-PINNs 中生效的主要原因.

为了降低计算图的规模,可在尽可能保证计算精确的情况下,由软约束条件限定化学能的值.改用软约束条件要求神经网络将化学能作为输出,故采用MIM 思想改进后的神经网络映射关系如下

根据改进后的映射关系与物理约束,可建立如图1 所示的PF-PINNs 框架,图中的 σ 表示非线性激活函数,本文选用如下形式的Swish 函数:σ(Y)=在每一个迭代步中,程序首先通过神经网络得到时空坐标到物理量的函数映射关系;随后使用自动微分计算物理量的各阶偏导数;再由物理方程和初边界条件确定总残差值,优化器根据总残差值对神经网络的参数进行优化;上述过程将被反复执行,直到每次随机采样的总残差值不再下降或迭代步达到最大,此时可认为训练完成.

图1 改进的基于相场法的物理融合神经网络(PF-PINNs)简图Fig.1 Sketch of modified physics-informed neural networks for phase-field (PF-PINNs) method

根据控制方程式(1)～式(4),可给出各方程残差项的表达形式

其中,∂a表示而 ∂ab表示LM,LU,LV,LC,Lφ为各控制方程的残差项,其分别对应式(1)～式(4)所对应的五个标量形式方程式,包括连续性方程、两个方向的动量方程、界面演化方程与化学能方程.神经网络的训练过程主要是在离散采样点上完成的,故式(10)～式(14)均需表示为离散形式的误差,以离散采样点取值表示的总误差、方程误差、初始条件误差和边界条件误差可写为

其中,Ltotal,LEqn,LICs,LBCs分别为总误差、方程误差、初始条件误差和边界条件误差.下标“Eqn”、“ICs”和“BCs”分别代表物理量或坐标在内部区域、初始条件和边界条件上的取值,上标“i”代表的是第i个采样点上的信息.NEqn,NICs,NBCs分别为内场区域、初值条件和边值条件上采样点的个数,物理量在第i个初边界采样点上的取值则根据物理问题的初边值条件确定.

离散形式的总误差Ltotal确定以后,神经网络中权重和偏差的取值 Θ 可通过随机梯度下降算法更新,该算法调节 Θ 的取值使总误差下降.本文中采用Adam 方法[28]作为随机梯度下降的优化算法,并在训练开始时使用Xavier 初始化方法[29]随机初始化神经网络的权重和偏差.完成一轮优化后,神经网络将重新在采样空间中进行采样,并重复上述流程,直到神经网络达到最大的训练次数或是总误差降低到容许值内.整个算法训练的流程示意图如图2 所示.

图2 基于相场法的物理融合神经网络(PF-PINNs)训练架构图Fig.2 Framework of physics-informed neural networks for phase-field(PF-PINNs) method

3 针对RT 不稳定性问题的神经网络训练结果

3.1 算例设置

RT 不稳定性问题的计算域和边界条件如图3 所示,其中整个区域的空间大小为Ω ∈[0,0.5] m×[0,4] m,计算总时长为t∈[0,3] s.计算区域上方为重流体L,下方为轻流体G,两流体间的初始分界线可定义为一正弦函数形式的扰动,其扰动函数为h=-0.1cos(2πx).根据扰动函数的物理含义,可定义初始时刻相分数场满足的关系式为

图3 RT 不稳定性问题的计算区域和初始条件Fig.3 Computational domain and initial condition of RT instability

假设初始时刻全场的速度为 0 m/s,则所有的初始条件已完备,扰动将在重力的影响下随时间发展.计算区域的左侧和右侧设置为滑移边界条件,上端和下端则设置为无滑移边界条件.RT 不稳定性问题的物理参数以列表的形式展示于表2 中,本问题的关键无量纲数为雷诺数Re=ρLDU/μL和阿特伍德数其中D为特征长度,为本问题的特征速度;本问题中特征长度为1 m,特征速度为1 m/s.重流体密度为 ρL=3 kg/m3,轻流体密 度则 为 ρG=1 kg/m3,两种流体的黏度均为μL=μG=0.001 N·s/m2,重力加速度为g=(0,1) m/s2,本问题中不考虑表面张力的作用,可认为 σ0=0.界面相关的参数包括界面厚度和界面渗透率,参考之前的文献研究,界面厚度可设为 ε=0.01 m,界面渗透率则设为M0=10 μm·(N·s)-1.

表2 RT 不稳定性算例中的物理参数和无量纲值Table 2 Physical parameters and non-dimensional value in RT instability

本文用于计算RT 不稳定性的神经网络包含15 个隐藏层,每层由100 个神经元组成.该神经网络的初始学习率设置为10-3,并随着循环次数的增加而下降至10-5,总循环次数为20 万次.本问题的均匀网格在时空采样区域 (x,y,t) 中等间距生成,为了保证该问题中采样尺寸足够合适,在正式开始计算之前本文对不同尺度的均匀网格工况进行了验证.本文在验证时分别选取了细、中、粗三套均匀网格,每套网格在x和y方向的间隔如表3 所示.本文以高精度谱方法在t=2.75 s 的速度场和相场作为初始值,随后用神经网络计算t=2.75～ 3 s 的流场结果,并于图4 中给出了三种方法在t=3 s 时的界面位置,以证明三种方法在复杂流场下最终计算结果基本一致.结果表明,三个网格数量下,神经网络对相分数场的计算结果相同.由于三个网格数量下,计算结果没有显著的差别,为了提升收敛速度,本文采用细网格作为所选网格,选择的均匀网格在空间两个方向上的数量分别为nx,ny=(65,512),即网格尺寸分别为 Δx=7.7 mm,Δy=7.7 mm;均匀网格在时间上的时间步长为Δt=2 ms.

图4 t=3 s 时,不同网格密度下PF-PINNs 得到的界面外形.其中红色实细线为粗网格的结果,蓝色虚细线为中等网格的结果,黑色散点为细网格的结果Fig.4 The shape of interface in various intervals when t=3 s.The red solid line denotes the result from coarse mesh,the blue dashed line denotes the result from medium mesh,and the black scattered line denotes the result from fine mesh

表3 均匀网格间隔对结果的影响Table 3 Influence of the interval in x direction and y direction on a uniform grid

本问题的离散采样点则从均匀网格中随机抽取,在每次循环开始时,PF-PINNs 从采样区域内部取60000 个点作为内部点,从t=0 时刻抽取30000个点作为初值条件采样点,从每个边界各取2000 个点作为边界条件采样点,共计8000 个边界采样点,即=(60000,30000,2000×4).

在神经网络训练过程中,本文使用了时间分割策略.使用多个神经网络依次训练的方法借鉴了Wight 等[30]提出的“time marching strategy”的思想,第一个网络训练时间区间t∈[0,0.5] s 的结果,待第一个网络训练收敛后,将第一个网络在t=0.5 s 时的预测值传出至第二个网络,作为第二个网络的初值条件.第二个网络负责训练时间区间t∈[0.5,1] s 的结果,训练完成后将结果传给第三个网络,依此类推,直到所有时间区间上的网络均被训练.根据RT 不稳定性问题的特性,本文采用了如下的分割策略:当t∈[0,1.5] s 时,以0.5 s 作为时间间隔分割计算区域;当t∈[1.5,3] s 时,以0.25 s 作为时间间隔分割计算区域.综上,本问题需要训练9 个神经网络以顺利求解,其占据的时间区间分别为 [0,0.5] s,[0.5,1] s,[1,1.5]s 以及[1.5,1.75]s,[1.75,2]s,[2,2.25]s,[2.25,2.5]s,[2.5,2.75]s和 [2.75,3] s.这样的设置是由于RT 不稳定性问题的物理特性所导致的,当t∈[0,1.5] s 时界面形态较为稳定,流场特征并不丰富,采用较大的时间间隔可更高效地进行计算.但当t∈[1.5,3] s 时界面迅速变形,速度场变化较大,为保证神经网络结果收敛,需要缩小每个网络的时间区间范围.时间区间的长度需要根据物理问题的特征合理控制,以平衡神经网络训练的效率与精度.

3.2 训练结果

RT 不稳定性问题的界面变化情况如图5 所示.其中图5(a)为改进的PF-PINNs 计算结果界面变化,界面处的初始扰动随着时间的增加而扩大,重流体从左侧边界的凹陷处下沉,轻流体则从右侧的凸起处上浮,分别形成了尖钉和气泡.在t∈[1.5,2] s 时,界面处的流体发生剧烈卷曲,在尖钉附近形成了复杂的涡结构,随轻重流体的反向运动被拉伸并向外扩张,重流体旋流结构非常复杂.图5(b)给出了高精度谱元法[32]的计算结果用作对比,其与图5(a)形态非常接近,表明改进的PF-PINNs 方法可接近传统算法的计算精度.图6 还给出了神经网络结果与高精度谱元法的计算结果在t=3.0 s 时 界面附近((x,y)∈[0,0.5]×[-1,0])速度场矢量的对比以及t=3.0 s 时两种算法的界面位置对比,结果表明,改进的PF-PINNs 计算的涡结构与传统算法得到的涡结构形状吻合,且流场中速度矢量方向与大小基本一致,矢量场各处疏密程度无明显差别.改进的PF-PINNs 的计算结果与高精度谱元法相比,在卷曲处的界面演化稍有差别,但在流场中大部分区域内界面形状保持一致,这说明了改进的PF-PINNs 的计算结果符合物理规律.

图5 RT 不稳定性问题的相分数演化和文献[33]结果的对比Fig.5 Evolution of volume fraction of RT instability and comparison with Ref.[33]

图6 在 t=3.0 s 时RT 不稳定性问题中:(a) 改进的PF-PINNs,(b) 高精度谱元法的速度矢量对比和(c) 改进的PF-PINNs 与高精度谱元法的界面形态对比.轻重流体交界面由 C=-0.8 时的等值线定义,速度矢量的大小在图中由颜色条表示.图6(c)中红线为改进的PF-PINNs 的界面,蓝线为高精度谱元法的界面Fig.6 Comparison of the velocity vector when t=3.0 s between (a) modified PF-PINNs,(b) high-order spectral element method and (c) interface evolution in RT instability.The interface between heavier fluid and lighter fluid are defined using isoline C=-0.8,and the magnitude of velocity vector is shown with color bar.The red line is the interface from the modified PF-PINNs,while the blue line is the interface from the high-order spectral element method in Fig.6(c)

为了定量分析神经网络的计算结果,本文以高精度谱元法的结果作为精确解,计算了神经网络结果相较于高精度谱元法在各个典型时刻的相对误差,以及各时刻计算结束时的总误差Ltotal,并列于表4 中.本文的相对误差定义为其中qnn为神经网络中物理量q的结果,qaccurate为精确解中物理量q的结果.通过表4 的结果可以看出,神经网络训练结果在0～ 2.5 s 内均维持着较好的精度,在2.5～ 3 s 间由于速度场的非线性演化导致误差快速上升.总体而言,随着训练时间的推移,神经网络的总误差值在不断上升,这也导致了速度相对误差的不断上升,当误差过大时计算结果将偏离真实解.本文中训练次数、学习率等参数是一直不变的,当界面流动趋于复杂时,采用原有的超参数应对当前问题可能无法保证结果有较好的收敛性.适当调节每个子网络在不同时刻的超参数可以缓解这种现象,提升神经网络的精度.图7 展示了尖钉和气泡在y方向上的位置随时间的变化曲线,并与文献[31-33]给出的结果进行了对比.其中红色线条为本文神经网络计算所得结果,其他形状的散点为参考文献给出的结果.对比结果表明,改进的PF-PINNs的模拟结果与参考文献相吻合,体现了改进的PFPINNs 在宏观结果上的计算精度.

表4 各个典型时刻神经网络与高精度谱元法的相对误差和训练结束时的总误差Table 4 Relative errors in various typical times between neural network and accurate spectral element method and the total loss when training ends

图7 气泡与尖钉位置随时间演化过程与参考文献结果对比Fig.7 Comparison of the position of bubble and spike between references and modified PF-PINNs

3.3 改进的PF-PINNs 加速效果分析

本文第2 节指出,深度混合残差方法[27]能够显著提升高阶微分方程的计算速度,同时保证神经网络的精度.为了验证改进方法在简单流场与复杂流场中的效果,本文对比了原始PF-PINNs 与改进的PF-PINNs 在t∈[0,0.5] s 和t∈[2.75,3.0] s 时总误差随计算时间的变化曲线以及总误差随迭代次数的变化曲线,相关数据列于图8 与图9 中.用于对比的神经网络大小为10 层100 节点,学习率为10-3,迭代次数为80000 次,其余设置与3.1 节中给出的设置完全一致.

由图8(a)与图9(a)的结果可以看出,在其他所有条件相同的情况下,使用MIM 与不使用MIM 得到的总误差没有太大的区别,两者之间的总误差均在同一量级且差别不大.但在图8(b)的结果中,在其他所有条件相同的情况下,采用MIM 策略训练用时为40720 s,平均单步用时为0.509 s;而未采用MIM 策略训练用时为115652 s,平均单步用时为1.446 s;采用MIM 策略使平均单步用时减少了64.80%.在图9(b)的结果中,在其他所有条件相同的情况下,采用MIM 策略训练用时为34281 s,平均单步用时为0.429 s;而未采用MIM 策略训练用时为113640 s,平均单步用时为1.420 s;采用MIM 策略使平均单步用时减少了69.79%.以上结果证实了深度混合残差方法在计算初始阶段和流场复杂演化阶段的有效性与可靠性.

图8 t ∈[0,0.5] s 时采用了混合残差方法(红线)与未采用混合残差方法(蓝线)的PF-PINNs:(a) 总残差随迭代次数的变化,(b) 总残差随时间的变化.其中图8 (b)中的时间指的是GPU 计算用时Fig.8 Variations of total loss with (a) iterations and (b) times using MIM (red line) and no MIM (blue line) in t ∈[0,0.5] s.The label“time”in Fig.8(b) means the GPU computational time of PF-PINNs

图9 t ∈[2.75,3.0] s 时采用了混合残差方法(红线)与未采用混合残差方法(蓝线)的PF-PINNs:(a) 总残差随迭代次数的变化,(b) 总残差随时间的变化.其中(b)中的时间指的是GPU 计算用时Fig.9 Variations of total loss with (a) iterations and (b) times using MIM (red line) and no MIM (blue line) in t ∈[2.75,3.0] s.The label“time”in (b) means the GPU computational time of PF-PINNs

4 结论

为了实现两相流场界面的高精度和高效捕捉,本文进一步改进了基于相场法的物理融合神经网络(PF-PINNs),参考深度混合残差方法,将辅助变量化学能作为神经网络的独立输出,从而在相分数输运方程中拆分出关于化学能的额外控制方程,同时考虑两个方程的损失函数..

本文基于改进的PF-PINNs,求解了经典的Rayleigh-Taylor 不稳定问题,从而考核当存在强非线性因素条件下这一方法的求解能力、精度和效率.与高精度谱元法结果对比表明,本文的改进的PF-PINNs 有能力捕捉到以RT 不稳定为代表的两相界面的强非线性演化过程,计算结果符合物理规律.尖钉和气泡位移结果与参考文献结果定量吻合良好.与此同时,所采用的深度混合残差方法能够有效地提升计算速度,RT 不稳定算例中,在其他所有条件相同的情况下,采用MIM 策略可使训练用时减少60%以上,在保证精度的前提下带来了较明显的效率提升.

上述结果表明,PF-PINNs 方法以及深度混合残差改进方法可用于精细求解两相流界面问题.考虑到相比CFD 方法,PINNs 方法以及神经网络方法在并行效率、软硬件结合等方面具有较强的优势,在反问题求解和数据融合等方面均具有较好的应用灵活性,因此PINNs 相关方法在未来应用的潜力很大.本文工作可为进一步提升神经网络的训练速度、探索高精度数据同化方法等提供参考.