基于SPH方法與能量最小化算法平衡態(tài)模擬分析

梁廣兵,羅賢兵,孫樹瑜

(1. 貴州大學(xué) 數(shù)學(xué)與統(tǒng)計學(xué)院,貴州 貴陽 550025; 2.阿卜杜拉國王科技大學(xué) 計算傳遞現(xiàn)象實驗室,沙特阿拉伯 圖瓦 23955-6900)

0 引 言

能量是物理學(xué)上的一個重要概念。一切物理過程都基于能量的變化,因此基于能量對物理過程進(jìn)行描述和模擬是科學(xué)研究的熱點,能夠幫助人們理解物理過程最底層的原理和規(guī)律。不同形式的能量在不同學(xué)科領(lǐng)域中都有著重要的應(yīng)用。

文獻(xiàn)[1]介紹了雙阱勢(double-well potential,DWP)是量子力學(xué)中最重要的勢之一,它在量子信息論中已經(jīng)變得非常重要;文獻(xiàn)[2]擴展了解析傳遞矩陣方法來求解任意對稱雙阱勢中的能量分裂,同時提出了一種顯式的、并與分裂能級相對應(yīng)的色散方程,數(shù)值計算表明,該方法可以對對稱雙阱勢給出極精確的結(jié)果;文獻(xiàn)[3]在相場模型中研究了具有雙阱勢的艾倫-卡恩方程的數(shù)值逼近,提出了一種新的線性、能量穩(wěn)定和保持最大值原理的方案,結(jié)合了最近發(fā)展的能量分解方法來半隱式處理雙阱勢函數(shù),該方法確保了擴大相變域內(nèi)的能量不等式;文獻(xiàn)[4]關(guān)于高階多項式自由能的艾倫-卡恩方程提出了一種無條件穩(wěn)定的數(shù)值格式,且根據(jù)雙阱勢和梯度項組成的總能量泛函的梯度流推導(dǎo)出對流方程,該方程已被廣泛用于圖像處理、樹突生長、平均曲率運動和多相流體流動。因此,關(guān)于雙阱勢函數(shù)的研究具有一定的實際意義。

在流體動力學(xué)中,常用粒子來描述某種狀態(tài)方程中流體中所攜帶的物質(zhì),系統(tǒng)在達(dá)到平衡態(tài)時內(nèi)部會形成一些能量最低的分布。在所有的自然定律當(dāng)中,能量最小原理是最重要的定律之一[5]。能量越低,系統(tǒng)所處的狀態(tài)越穩(wěn)定。因此一個系統(tǒng)總是要調(diào)整自己,使系統(tǒng)的總能量達(dá)到最低,使自己處于穩(wěn)定的平衡態(tài)。基于該思想,在一維不考慮粒子速度的情況下,以每個粒子的位置作為輸入,以系統(tǒng)的總能量作為輸出來研究系統(tǒng)平衡態(tài)問題。首先給定粒子初始位置,對應(yīng)一個初始能量分布。然后按照一定的距離步長讓粒子運動,去尋找所有粒子的最佳位置。在此運動過程中,當(dāng)系統(tǒng)能量降低時,接受粒子位置;當(dāng)系統(tǒng)能量升高時,拒絕粒子當(dāng)前位置,然后粒子繼續(xù)運動。整個系統(tǒng)一直向能量最低的地方變化,當(dāng)系統(tǒng)能量達(dá)到最低值的附近時, 此時粒子就不再運動。基于光滑粒子流體動力學(xué)(smooth particle hydrodynamics,SPH)理論,當(dāng)粒子所處位置使得系統(tǒng)能量最小時,粒子分布狀態(tài)代表了系統(tǒng)的平衡態(tài)。

SPH是一種拉格朗日無網(wǎng)格粒子方法,已經(jīng)成功應(yīng)用于工程和科學(xué)等眾多領(lǐng)域,它由文獻(xiàn)[6-7]分別提出,最初用于求解三維開放空間天體物理學(xué)問題。隨著文獻(xiàn)[8-9]成功地模擬了自由表面流動和多相流,SPH方法逐漸擴展到流體的各個領(lǐng)域。SPH方法與有限元、有限差分法等網(wǎng)格法不同,它無需對計算區(qū)域進(jìn)行網(wǎng)格劃分,而是將計算區(qū)域離散為一系列有相互作用的粒子。這些粒子攜帶流體的各種特性,包括質(zhì)量、密度、速度、加速度、能量等,粒子間通過核函數(shù)進(jìn)行相互作用。從而在處理大變形、爆炸、變形邊界以及界面追蹤等問題具有天然的優(yōu)勢。基于SPH方法的思想,對雙阱勢能量泛函模型進(jìn)行離散,得到易于求解的數(shù)值格式。最初的想法是應(yīng)用蒙特卡洛分子模擬求解,其原理是讓每個粒子做熱運動,在此過程中,總能量低時接受粒子位置分布,總能量高時不接受。但該方法的收斂速度太慢而且計算耗時過大,因此可行性不高。 為解決此問題,出現(xiàn)了梯度下降法[10]、Adam算法[11]以及模擬退火算法等,這些算法已被廣泛應(yīng)用于解決機器學(xué)習(xí)中產(chǎn)生的大規(guī)模優(yōu)化問題,不需要每次迭代求解線性系統(tǒng),具有類似牛頓方法的優(yōu)勢[12]。經(jīng)過對算法以及其他因素的考慮,本文選擇梯度下降法作為數(shù)值模擬方法,同時對梯度下降法作詳細(xì)的描述,以簡單的回歸方程對梯度下降法進(jìn)行檢驗。

本文為驗證該模型和算法設(shè)計了2個算例。為減少邊界效應(yīng),采用周期邊界條件進(jìn)行模擬。計算結(jié)果表明,該雙阱勢能量模型結(jié)合SPH、梯度下降法的計算效率較高,結(jié)果準(zhǔn)確。

1 SPH方法

1.1 核近似

對任意連續(xù)函數(shù)f(x),可由狄拉克函數(shù)將其表示為[9]:

(1)

其中,Ω為積分區(qū)域。

若用光滑核函數(shù)W來代替狄拉克函數(shù),則核近似表示為:

(2)

光滑核函數(shù)W應(yīng)具備以下性質(zhì)[13]:

(1) 歸一化。計算公式為:

(3)

(2) 緊支性。計算公式為:

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

(4)

(3) 狄拉克性。計算公式為:

(5)

(4) 偶函數(shù)性。

常用核函數(shù)如下。

(1) 高斯型核函數(shù)。表達(dá)式為:

(6)

其中,對應(yīng)一維、二維和三維問題,αd的值分別為1/π1/2h、1/πh2、1/π3/2h2。高斯型核函數(shù)是充分光滑的,即使對于高階導(dǎo)數(shù),它也被認(rèn)為是一種極不錯的選擇,因其很穩(wěn)定且精度高,特別對于不規(guī)則粒子分布的情況。

(2) 三次樣條型核函數(shù)。表達(dá)式為:

(7)

其中:對應(yīng)一維、二維和三維問題,αd的值分別為1/h、,15/7πh2、3/2πh3;R=|xi-xj|/h為兩粒子之間相對距離;h為光滑長度,決定了核函數(shù)支持域的大小。

利用分部積分、高斯公式以及應(yīng)用光滑函數(shù)的緊支性[14],可推導(dǎo)函數(shù)導(dǎo)數(shù)的核近似式如下:

(8)

同理,可推導(dǎo)函數(shù)二階導(dǎo)數(shù)核近似式為:

(9)

本文選擇鐘形函數(shù)[6]作為核函數(shù),其表達(dá)式為:

(10)

其中,αd在一維、二維和三維空間中的值分別為5/4h、5/πh2、105/16πh3。光滑核函數(shù)及其一階導(dǎo)數(shù)圖像如圖1所示。

圖1 光滑核函數(shù)及其一階導(dǎo)數(shù)

1.2 粒子近似

可將(2)式轉(zhuǎn)換為支持域內(nèi)所有粒子加權(quán)求和的離散形式,即

(11)

函數(shù)導(dǎo)數(shù)的離散形式為:

(12)

2 能量模型

雙阱勢是量子力學(xué)中最重要的勢之一,為理解各種量子現(xiàn)象提供了許多有價值的見解[15]。雙阱勢也稱為自由能函數(shù),同時在相場模型[16-18]的應(yīng)用中扮演著至關(guān)重要的角色。本文基于雙阱勢函數(shù),考慮流體性質(zhì),建立數(shù)學(xué)模型。

2.1 模型推導(dǎo)

本文引入相位函數(shù)來確定2種流體所占據(jù)的區(qū)域:

(13)

設(shè)Ginzburg-Landau[19]雙阱勢函數(shù)為:

(14)

并定義混合能量泛函，即

(15)

其中:λ為混合能量密度;η為界面厚度的毛細(xì)管寬度。

考慮一維空間中的氣液相問題,因為(φ2-1)2=(φ-1)2(φ+1)2,根據(jù)(14)式,定義:

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

(16)

其中,ρL、ρV分別為液相和氣相的密度。

最后得到修正的雙阱勢函數(shù)以及能量函數(shù),其表達(dá)式如下:

(17)

其中:第1項f0為修正的雙阱勢能量密度;fd為密度梯度項,它對能量也有貢獻(xiàn),用來表征兩相中界面的能量,稱之為表面張力能,這個能量會導(dǎo)致氣液相的面積越小越好。

2.2 模型離散

(18)

其中的未知量為密度ρ。根據(jù)SPH方法的近似和離散思想,由(2)式和(11)式可推導(dǎo)出密度的求解公式為:

(19)

(19)式被稱之為“密度求和法”,它遵循物理中的質(zhì)量守恒,因此被廣泛應(yīng)用于工程中。同時根據(jù)(12)式可得密度梯度的表達(dá)式為:

(20)

3 梯度下降法

梯度下降法,通常也稱為最速下降法,是求解無約束優(yōu)化問題最簡單和最古老的方法之一,目前是機器學(xué)習(xí)和深度學(xué)習(xí)中最常用的優(yōu)化策略。它通過尋找一個函數(shù)的參數(shù)值,使損失函數(shù)盡可能地最小化,已經(jīng)成為人工智能和計算機科學(xué)應(yīng)用的有力工具。

設(shè)線性回歸方程為:

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

(21)

需要確定θ這組向量的選擇方法,才能使得擬合函數(shù)與目標(biāo)函數(shù)y(x)盡可能地接近。因此,選擇均方誤差作為目標(biāo)函數(shù)的損失函數(shù)[20]，即

(22)

通過尋找最優(yōu)的θ,使得損失函數(shù)J(θ)達(dá)到最小。若J(θ)的值取0,則表明找到了最優(yōu)的θ,構(gòu)造了擬合度極高的函數(shù),但這基本上是達(dá)不到的,只能使損失函數(shù)無限接近0,當(dāng)滿足一定的精度之后停止迭代。其求解過程分為2步。

(1) 計算損失函數(shù)J(θ)的導(dǎo)數(shù)。計算公式為:

(23)

(2) 參數(shù)θ以負(fù)梯度方向進(jìn)行更新。計算公式為:

(24)

其中,α為學(xué)習(xí)率。

若α取值過小，則會導(dǎo)致收斂速度過慢;若α過大，則損失函數(shù)可能不會收斂,值反而越來越大。因此合理地選擇學(xué)習(xí)率,能夠更快地得到收斂結(jié)果。

下面介紹模型的求解。

(25)

其中:Δx=10-6h,h為光滑長度;α=0.1h為學(xué)習(xí)率。

給定粒子初始位置x0,可計算系統(tǒng)初始能量E(x0),然后由梯度下降法進(jìn)行位置更新得到新的粒子位置x0′,計算系統(tǒng)的新能量E(x0′)。若E(x0′)

4 數(shù)值算例

通過以下2個算例,以Python軟件進(jìn)行編程驗證數(shù)值模型的可行性、SPH方法以及梯度下降法的有效性。

4.1 算例1

在一維空間上考慮氣液相平衡問題。由熱力學(xué)第二定律可知,當(dāng)氣液相平衡時,系統(tǒng)的能量達(dá)到最低。本文尋找系統(tǒng)的能量最低值,并求出所對應(yīng)的參數(shù)值。該算例中所有的參數(shù)均是無量綱的。

初始時,將氣液相的密度分別設(shè)置為ρV=0.5、ρL=1.0,氣相粒子數(shù)為nV=10顆,液相粒子數(shù)為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段組成。在后面計算密度以及密度梯度項時均考慮了周期邊界條件。在原來粒子位置分布的基礎(chǔ)上添加了左端和右端以后,粒子的位置分布變成了周期分布,如圖2b所示。

(26)

其中:每個粒子的質(zhì)量均設(shè)為1;光滑函數(shù)為鐘形函數(shù);光滑長度取h=5;總的循環(huán)次數(shù)設(shè)置為300 000次。同時設(shè)置ρL=0.9和ρV=0.6為氣液兩相平衡時各相密度,因此在求解的過程中系統(tǒng)中的氣液相密度接近ρL=0.9和ρV=0.6時,總能量也變成最低,從而達(dá)到系統(tǒng)的平衡態(tài)。

4.1.1 數(shù)值結(jié)果

(1) 粒子位置變化。粒子初始位置、中間變化位置以及平衡后的位置如圖3所示。由圖3可知,粒子由最初位置開始不斷地左右運動,尋找粒子的最佳位置,最后達(dá)到穩(wěn)定的分布,系統(tǒng)的能量達(dá)到最低。

圖3 不同迭代次數(shù)的粒子位置分布

(2) 氣液相密度變化。初始時氣液相密度分布、中間變化的密度分布以及平衡時的氣液相密度分布如圖4所示。

因為粒子的位置決定了氣液相密度,所以當(dāng)粒子找到了最佳位置,氣液相密度也接近所設(shè)置的兩相平衡時的值。由數(shù)值結(jié)果可知,ρmax=0.899 93,ρmin=0.599 63,相比于氣液兩相平衡時的密度值,其數(shù)值誤差非常小。

(3) 能量變化。由于能量是由密度及密度梯度項共同確定,能量圖像反映的是:由初始粒子位置對應(yīng)初始?xì)庖合嗝芏?初始密度計算初始能量,當(dāng)粒子開始運動,對應(yīng)的密度也將發(fā)生變化,導(dǎo)致能量也變化。因為目標(biāo)是尋找系統(tǒng)的平衡態(tài),所以系統(tǒng)的能量一定會不斷降低。

圖4 算例1兩相密度的變化

能量變化如圖5所示。由圖5可知, 能量值在迭代次數(shù)T=100之前快速地衰減,隨后進(jìn)行緩慢的變化,最后穩(wěn)定在一個范圍,不再有較大的變化,即系統(tǒng)的能量從E=0.153 04衰減到E=0.012 88,認(rèn)為其達(dá)到了系統(tǒng)的平衡態(tài),即達(dá)到能量最小值。

圖5 算例1能量變化趨勢

4.2 算例2

初始時,氣液兩相的密度分別設(shè)置為ρV=0.5、ρL=1.0。每個粒子的質(zhì)量取值為1,一共有60顆粒子,總質(zhì)量mt=60。氣液相的體積Vt=80,位置分布以及其他參數(shù)均與算例1保持一致,同時將氣液相平衡時密度值設(shè)置為ρV=ρL=3/4。該算例主要是通過描述兩相變一相的物理過程來驗證系統(tǒng)的質(zhì)量守恒。因為總的密度ρt=mt/Vt=3/4,所以系統(tǒng)趨于平衡時的密度值應(yīng)為ρV=0.5,從而滿足質(zhì)量守恒定律。

4.2.1 數(shù)值結(jié)果

(1) 粒子位置變化。粒子的位置分布如圖6所示,從圖6可以看出，系統(tǒng)平衡時所有粒子均呈均勻分布,每個粒子之間的距離相等,該位置分布也體現(xiàn)了系統(tǒng)由兩相變成了一相的原理。

圖6 算例2不同迭代次數(shù)的粒子位置分布

(2) 氣液相密度變化。氣液相密度變化如圖7所示,從圖7可以觀察密度曲線的變化趨勢,系統(tǒng)由兩相變成了一相。系統(tǒng)由最初的密度ρV=0.5,ρL=1.0變成ρ=0.749 98,因此由密度的變化過程可以驗證物理上的質(zhì)量守恒定律,從而表明SPH方法中“密度求和法”的有效性。

(3) 能量變化。能量變化如圖8所示。系統(tǒng)初始總能量為E=0.315 95,由于粒子位置不斷變化,使系統(tǒng)總能量不斷降低,最終達(dá)到能量值E=1.06E-07之后就變化得非常緩慢,以至于將其認(rèn)為達(dá)到了系統(tǒng)的平衡態(tài)。此時,系統(tǒng)的粒子將呈現(xiàn)出均勻分布,密度值為ρ=0.749 98。

圖7 算例2兩相密度的變化

圖8 算例2能量變化趨勢

5 結(jié) 論

本文基于SPH方法的思想和離散原理,將其成功地應(yīng)用到本文的模型中， 最后應(yīng)用梯度下降法對能量模型進(jìn)行了數(shù)值模擬。數(shù)值實驗結(jié)果表明,SPH方法和梯度下降法具有很高的可行性和準(zhǔn)確性。

基于相場模型,以SPH方法對模型進(jìn)行離散,最后應(yīng)用梯度下降法進(jìn)行數(shù)值模擬,從而將三者聯(lián)系在一起,為后續(xù)拓展到高維空間的工作以及基于范德瓦爾斯方程或彭-羅賓遜方程的氣液相平衡問題的研究奠定基礎(chǔ)。