基于SPH方法与能量最小化算法平衡态模拟分析

梁广兵,罗贤兵,孙树瑜

(1. 贵州大学 数学与统计学院,贵州 贵阳 550025; 2.阿卜杜拉国王科技大学 计算传递现象实验室,沙特阿拉伯 图瓦 23955-6900)

0 引 言

能量是物理学上的一个重要概念。一切物理过程都基于能量的变化,因此基于能量对物理过程进行描述和模拟是科学研究的热点,能够帮助人们理解物理过程最底层的原理和规律。不同形式的能量在不同学科领域中都有着重要的应用。

文献[1]介绍了双阱势(double-well potential,DWP)是量子力学中最重要的势之一,它在量子信息论中已经变得非常重要;文献[2]扩展了解析传递矩阵方法来求解任意对称双阱势中的能量分裂,同时提出了一种显式的、并与分裂能级相对应的色散方程,数值计算表明,该方法可以对对称双阱势给出极精确的结果;文献[3]在相场模型中研究了具有双阱势的艾伦-卡恩方程的数值逼近,提出了一种新的线性、能量稳定和保持最大值原理的方案,结合了最近发展的能量分解方法来半隐式处理双阱势函数,该方法确保了扩大相变域内的能量不等式;文献[4]关于高阶多项式自由能的艾伦-卡恩方程提出了一种无条件稳定的数值格式,且根据双阱势和梯度项组成的总能量泛函的梯度流推导出对流方程,该方程已被广泛用于图像处理、树突生长、平均曲率运动和多相流体流动。因此,关于双阱势函数的研究具有一定的实际意义。

在流体动力学中,常用粒子来描述某种状态方程中流体中所携带的物质,系统在达到平衡态时内部会形成一些能量最低的分布。在所有的自然定律当中,能量最小原理是最重要的定律之一[5]。能量越低,系统所处的状态越稳定。因此一个系统总是要调整自己,使系统的总能量达到最低,使自己处于稳定的平衡态。基于该思想,在一维不考虑粒子速度的情况下,以每个粒子的位置作为输入,以系统的总能量作为输出来研究系统平衡态问题。首先给定粒子初始位置,对应一个初始能量分布。然后按照一定的距离步长让粒子运动,去寻找所有粒子的最佳位置。在此运动过程中,当系统能量降低时,接受粒子位置;当系统能量升高时,拒绝粒子当前位置,然后粒子继续运动。整个系统一直向能量最低的地方变化,当系统能量达到最低值的附近时, 此时粒子就不再运动。基于光滑粒子流体动力学(smooth particle hydrodynamics,SPH)理论,当粒子所处位置使得系统能量最小时,粒子分布状态代表了系统的平衡态。

SPH是一种拉格朗日无网格粒子方法,已经成功应用于工程和科学等众多领域,它由文献[6-7]分别提出,最初用于求解三维开放空间天体物理学问题。随着文献[8-9]成功地模拟了自由表面流动和多相流,SPH方法逐渐扩展到流体的各个领域。SPH方法与有限元、有限差分法等网格法不同,它无需对计算区域进行网格划分,而是将计算区域离散为一系列有相互作用的粒子。这些粒子携带流体的各种特性,包括质量、密度、速度、加速度、能量等,粒子间通过核函数进行相互作用。从而在处理大变形、爆炸、变形边界以及界面追踪等问题具有天然的优势。基于SPH方法的思想,对双阱势能量泛函模型进行离散,得到易于求解的数值格式。最初的想法是应用蒙特卡洛分子模拟求解,其原理是让每个粒子做热运动,在此过程中,总能量低时接受粒子位置分布,总能量高时不接受。但该方法的收敛速度太慢而且计算耗时过大,因此可行性不高。 为解决此问题,出现了梯度下降法[10]、Adam算法[11]以及模拟退火算法等,这些算法已被广泛应用于解决机器学习中产生的大规模优化问题,不需要每次迭代求解线性系统,具有类似牛顿方法的优势[12]。经过对算法以及其他因素的考虑,本文选择梯度下降法作为数值模拟方法,同时对梯度下降法作详细的描述,以简单的回归方程对梯度下降法进行检验。

本文为验证该模型和算法设计了2个算例。为减少边界效应,采用周期边界条件进行模拟。计算结果表明,该双阱势能量模型结合SPH、梯度下降法的计算效率较高,结果准确。

1 SPH方法

1.1 核近似

对任意连续函数f(x),可由狄拉克函数将其表示为[9]:

(1)

其中,Ω为积分区域。

若用光滑核函数W来代替狄拉克函数,则核近似表示为:

(2)

光滑核函数W应具备以下性质[13]:

(1) 归一化。计算公式为:

(3)

(2) 紧支性。计算公式为:

W(x-x′,h)=0,|x-x′|≥0

(4)

(3) 狄拉克性。计算公式为:

(5)

(4) 偶函数性。

常用核函数如下。

(1) 高斯型核函数。表达式为:

(6)

其中,对应一维、二维和三维问题,αd的值分别为1/π1/2h、1/πh2、1/π3/2h2。高斯型核函数是充分光滑的,即使对于高阶导数,它也被认为是一种极不错的选择,因其很稳定且精度高,特别对于不规则粒子分布的情况。

(2) 三次样条型核函数。表达式为:

(7)

其中:对应一维、二维和三维问题,αd的值分别为1/h、,15/7πh2、3/2πh3;R=|xi-xj|/h为两粒子之间相对距离;h为光滑长度,决定了核函数支持域的大小。

利用分部积分、高斯公式以及应用光滑函数的紧支性[14],可推导函数导数的核近似式如下:

(8)

同理,可推导函数二阶导数核近似式为:

(9)

本文选择钟形函数[6]作为核函数,其表达式为:

(10)

其中,αd在一维、二维和三维空间中的值分别为5/4h、5/πh2、105/16πh3。光滑核函数及其一阶导数图像如图1所示。

图1 光滑核函数及其一阶导数

1.2 粒子近似

可将(2)式转换为支持域内所有粒子加权求和的离散形式,即

(11)

函数导数的离散形式为:

(12)

2 能量模型

双阱势是量子力学中最重要的势之一,为理解各种量子现象提供了许多有价值的见解[15]。双阱势也称为自由能函数,同时在相场模型[16-18]的应用中扮演着至关重要的角色。本文基于双阱势函数,考虑流体性质,建立数学模型。

2.1 模型推导

本文引入相位函数来确定2种流体所占据的区域:

(13)

设Ginzburg-Landau[19]双阱势函数为:

(14)

并定义混合能量泛函，即

(15)

其中:λ为混合能量密度;η为界面厚度的毛细管宽度。

考虑一维空间中的气液相问题,因为(φ2-1)2=(φ-1)2(φ+1)2,根据(14)式,定义:

F(ρ)=(ρ-ρL)2(ρ-ρV)2

(16)

其中,ρL、ρV分别为液相和气相的密度。

最后得到修正的双阱势函数以及能量函数,其表达式如下:

(17)

其中:第1项f0为修正的双阱势能量密度;fd为密度梯度项,它对能量也有贡献,用来表征两相中界面的能量,称之为表面张力能,这个能量会导致气液相的面积越小越好。

2.2 模型离散

(18)

其中的未知量为密度ρ。根据SPH方法的近似和离散思想,由(2)式和(11)式可推导出密度的求解公式为:

(19)

(19)式被称之为“密度求和法”,它遵循物理中的质量守恒,因此被广泛应用于工程中。同时根据(12)式可得密度梯度的表达式为:

(20)

3 梯度下降法

梯度下降法,通常也称为最速下降法,是求解无约束优化问题最简单和最古老的方法之一,目前是机器学习和深度学习中最常用的优化策略。它通过寻找一个函数的参数值,使损失函数尽可能地最小化,已经成为人工智能和计算机科学应用的有力工具。

设线性回归方程为:

h(x)=θ0x0+θ1x1+…+θnxn

(21)

需要确定θ这组向量的选择方法,才能使得拟合函数与目标函数y(x)尽可能地接近。因此,选择均方误差作为目标函数的损失函数[20]，即

(22)

通过寻找最优的θ,使得损失函数J(θ)达到最小。若J(θ)的值取0,则表明找到了最优的θ,构造了拟合度极高的函数,但这基本上是达不到的,只能使损失函数无限接近0,当满足一定的精度之后停止迭代。其求解过程分为2步。

(1) 计算损失函数J(θ)的导数。计算公式为:

(23)

(2) 参数θ以负梯度方向进行更新。计算公式为:

(24)

其中,α为学习率。

若α取值过小，则会导致收敛速度过慢;若α过大，则损失函数可能不会收敛,值反而越来越大。因此合理地选择学习率,能够更快地得到收敛结果。

下面介绍模型的求解。

(25)

其中:Δx=10-6h,h为光滑长度;α=0.1h为学习率。

给定粒子初始位置x0,可计算系统初始能量E(x0),然后由梯度下降法进行位置更新得到新的粒子位置x0′,计算系统的新能量E(x0′)。若E(x0′)

4 数值算例

通过以下2个算例,以Python软件进行编程验证数值模型的可行性、SPH方法以及梯度下降法的有效性。

4.1 算例1

在一维空间上考虑气液相平衡问题。由热力学第二定律可知,当气液相平衡时,系统的能量达到最低。本文寻找系统的能量最低值,并求出所对应的参数值。该算例中所有的参数均是无量纲的。

初始时,将气液相的密度分别设置为ρV=0.5、ρL=1.0,气相粒子数为nV=10颗,液相粒子数为nL=40颗,从而求出气相粒子所占空间为vV=nV/ρV=20,液相粒子所占空间为vL=nL/ρL=40。同时单个气液相粒子体积为sV=vV/ρV=2,sL=vL/nL=1,可以分别给出气相和液相粒子的初始分布如下:气相1为(1-3-5-7-9-11-13-15-17-19)共 10颗粒子;液相(20.5-21.5-22.5-…-57.5-58.5-59.5)共40颗粒子;气相2为(61-63-65-67-69-71-73-75-77-79) 共计10颗粒子。粒子位置分布如图2所示。

图2 算例1粒子位置分布

同时,为了避免边界的影响,考虑周期边界条件,即在本身一维片段的左边和右边加上相同的片段,因此气液相粒子所占空间为-160～160,共由3段组成。在后面计算密度以及密度梯度项时均考虑了周期边界条件。在原来粒子位置分布的基础上添加了左端和右端以后,粒子的位置分布变成了周期分布,如图2b所示。

(26)

其中:每个粒子的质量均设为1;光滑函数为钟形函数;光滑长度取h=5;总的循环次数设置为300 000次。同时设置ρL=0.9和ρV=0.6为气液两相平衡时各相密度,因此在求解的过程中系统中的气液相密度接近ρL=0.9和ρV=0.6时,总能量也变成最低,从而达到系统的平衡态。

4.1.1 数值结果

(1) 粒子位置变化。粒子初始位置、中间变化位置以及平衡后的位置如图3所示。由图3可知,粒子由最初位置开始不断地左右运动,寻找粒子的最佳位置,最后达到稳定的分布,系统的能量达到最低。

图3 不同迭代次数的粒子位置分布

(2) 气液相密度变化。初始时气液相密度分布、中间变化的密度分布以及平衡时的气液相密度分布如图4所示。

因为粒子的位置决定了气液相密度,所以当粒子找到了最佳位置,气液相密度也接近所设置的两相平衡时的值。由数值结果可知,ρmax=0.899 93,ρmin=0.599 63,相比于气液两相平衡时的密度值,其数值误差非常小。

(3) 能量变化。由于能量是由密度及密度梯度项共同确定,能量图像反映的是:由初始粒子位置对应初始气液相密度,初始密度计算初始能量,当粒子开始运动,对应的密度也将发生变化,导致能量也变化。因为目标是寻找系统的平衡态,所以系统的能量一定会不断降低。

图4 算例1两相密度的变化

能量变化如图5所示。由图5可知, 能量值在迭代次数T=100之前快速地衰减,随后进行缓慢的变化,最后稳定在一个范围,不再有较大的变化,即系统的能量从E=0.153 04衰减到E=0.012 88,认为其达到了系统的平衡态,即达到能量最小值。

图5 算例1能量变化趋势

4.2 算例2

初始时,气液两相的密度分别设置为ρV=0.5、ρL=1.0。每个粒子的质量取值为1,一共有60颗粒子,总质量mt=60。气液相的体积Vt=80,位置分布以及其他参数均与算例1保持一致,同时将气液相平衡时密度值设置为ρV=ρL=3/4。该算例主要是通过描述两相变一相的物理过程来验证系统的质量守恒。因为总的密度ρt=mt/Vt=3/4,所以系统趋于平衡时的密度值应为ρV=0.5,从而满足质量守恒定律。

4.2.1 数值结果

(1) 粒子位置变化。粒子的位置分布如图6所示,从图6可以看出，系统平衡时所有粒子均呈均匀分布,每个粒子之间的距离相等,该位置分布也体现了系统由两相变成了一相的原理。

图6 算例2不同迭代次数的粒子位置分布

(2) 气液相密度变化。气液相密度变化如图7所示,从图7可以观察密度曲线的变化趋势,系统由两相变成了一相。系统由最初的密度ρV=0.5,ρL=1.0变成ρ=0.749 98,因此由密度的变化过程可以验证物理上的质量守恒定律,从而表明SPH方法中“密度求和法”的有效性。

(3) 能量变化。能量变化如图8所示。系统初始总能量为E=0.315 95,由于粒子位置不断变化,使系统总能量不断降低,最终达到能量值E=1.06E-07之后就变化得非常缓慢,以至于将其认为达到了系统的平衡态。此时,系统的粒子将呈现出均匀分布,密度值为ρ=0.749 98。

图7 算例2两相密度的变化

图8 算例2能量变化趋势

5 结 论

本文基于SPH方法的思想和离散原理,将其成功地应用到本文的模型中， 最后应用梯度下降法对能量模型进行了数值模拟。数值实验结果表明,SPH方法和梯度下降法具有很高的可行性和准确性。

基于相场模型,以SPH方法对模型进行离散,最后应用梯度下降法进行数值模拟,从而将三者联系在一起,为后续拓展到高维空间的工作以及基于范德瓦尔斯方程或彭-罗宾逊方程的气液相平衡问题的研究奠定基础。