王姍姍，張鑒

1. 中國(guó)科學(xué)院計(jì)算機(jī)網(wǎng)絡(luò)信息中心，北京 100190

2. 中國(guó)科學(xué)院大學(xué)，北京 100049

引言

相場(chǎng)法是近年來(lái)出現(xiàn)的一種強(qiáng)有力的計(jì)算方法，用于建模和預(yù)測(cè)材料的中尺度形態(tài)和微觀結(jié)構(gòu)演變[1]。相場(chǎng)法使用一組守恒和非守恒的在區(qū)域內(nèi)連續(xù)的場(chǎng)變量來(lái)描述微觀結(jié)構(gòu)。經(jīng)典的相場(chǎng)方程——Allen-Cahn 方程[2]和Cahn-Hilliard 方程[3]，描述了場(chǎng)變量在空間和時(shí)間上的改變。微觀結(jié)構(gòu)的繁雜演化可以被相場(chǎng)法有效的描述，而無(wú)需明確跟蹤界面的位置。相場(chǎng)法已成為模擬中尺度微觀結(jié)構(gòu)演化的一種重要而通用的方法。

對(duì)于經(jīng)典的Allen-Cahn方程和Cahn-Hilliard方程，已有基于有限差分格式的時(shí)間推進(jìn)求解方案。這些求解方案普遍存在著缺陷：很難設(shè)計(jì)高階的求解格式、時(shí)間步長(zhǎng)無(wú)法取到令人滿意的長(zhǎng)度等。本文所使用的緊致指數(shù)時(shí)間差分方法（compact Exponential Time Difference, cETD）[4-5]，對(duì)于經(jīng)典的相場(chǎng)方程的高階導(dǎo)數(shù)項(xiàng)，采用積分因子法進(jìn)行準(zhǔn)確積分，而對(duì)于它們的非線性項(xiàng)，又與經(jīng)典的積分因子法不同，緊致指數(shù)時(shí)間差分方法對(duì)方程的非線性部分利用多步逼近、龍格庫(kù)塔、預(yù)估校正等方法設(shè)計(jì)了高階的求解格式，并且利用算子分裂格式增強(qiáng)了整體的穩(wěn)定性。

微觀結(jié)構(gòu)演化的動(dòng)力是系統(tǒng)的總能量下降，系統(tǒng)的總自由能包括化學(xué)自由能、界面能、彈性應(yīng)變能、電磁能、靜電能等。其中，界面能各向異性和彈性應(yīng)變能是材料產(chǎn)生各向異性的主要因素。界面能是微結(jié)構(gòu)的界面上由于組成和結(jié)構(gòu)的不均勻產(chǎn)生的額外的能量，由于固體的晶體性質(zhì)，界面能通常是各向異性的，界面能各向異性的程度可以對(duì)粒子的生長(zhǎng)形態(tài)和平衡形狀產(chǎn)生重大影響，在相場(chǎng)方程中，界面能各向異性表現(xiàn)為包含非線性系數(shù)的高階空間導(dǎo)數(shù)。彈性應(yīng)變能產(chǎn)生于微結(jié)構(gòu)中相之間的晶格失配造成彈性位移，通過(guò)Khachaturyan 理論[6]在相場(chǎng)方程中引入彈性應(yīng)變能，將彈性應(yīng)變能表達(dá)為場(chǎng)變量的函數(shù)。已有的關(guān)于彈性應(yīng)變能計(jì)算的工作[7-8]，多為使用顯格式進(jìn)行計(jì)算的。

本文研究的重點(diǎn)是緊致指數(shù)時(shí)間差分方法中界面能各向異性和彈性應(yīng)變能的計(jì)算?；瘜W(xué)自由能在相場(chǎng)方程中主要表現(xiàn)為高階空間導(dǎo)數(shù)部分，已有的緊致指數(shù)時(shí)間差分方法通過(guò)適當(dāng)算子分裂格式處理化學(xué)自由能，將其歸于非線性項(xiàng)的計(jì)算，并證明了能量穩(wěn)定性[9-11]。本文在緊致指數(shù)時(shí)間差分方法框架下引入了界面能各向異性和彈性應(yīng)變能的計(jì)算，將界面能各向異性和彈性應(yīng)變能同樣歸于緊致指數(shù)時(shí)間差分方法的非線性項(xiàng)統(tǒng)一處理，為其設(shè)計(jì)了適當(dāng)?shù)乃阕臃至迅袷?，證明了算子分裂能夠保證能量穩(wěn)定。并通過(guò)鎳基合金以及Zr 的氫化物的材料腐蝕相場(chǎng)模型的數(shù)值實(shí)驗(yàn)，驗(yàn)證了該算法的正確性和能量穩(wěn)定性。

本文的組織結(jié)構(gòu)如下，第一章介紹了相場(chǎng)模型中界面能各向異性和彈性應(yīng)變能的引入；第二章主要闡述了含有界面能各向異性和彈性應(yīng)變能的相場(chǎng)模型的緊致指數(shù)時(shí)間差分方法，進(jìn)行了算子分裂格式的能量穩(wěn)定性以及求解格式的時(shí)間精度的證明；第三章給出了鎳基合金和Zr 的氫化物材料腐蝕相場(chǎng)模型的數(shù)值實(shí)驗(yàn)，并驗(yàn)證了求解格式的時(shí)間精度；第四章對(duì)本文所采用的含彈性應(yīng)變能的各向異性相場(chǎng)模型的指數(shù)時(shí)間差分方法進(jìn)行了總結(jié)。

1 含彈性應(yīng)變能的各向異性相場(chǎng)模型

在本文中，我們考慮含有界面能各向異性和彈性應(yīng)變能的耦合模型能量方程如式(1.1)所示。

其中，E是系統(tǒng)的總能量，這里考慮了總能量中的化學(xué)自由能Ec和彈性應(yīng)變能Eel?；瘜W(xué)自由能Ec可以表達(dá)為界面梯度能量項(xiàng)和局部自由能量密度的積分，如式（1.2）所示，其中界面梯度能量中含有各向異性。式(1.2)中，c是成分場(chǎng)變量，而為序參數(shù)，和是界面梯度能量項(xiàng)的系數(shù)，是3x3 的界面能各向異性矩陣，是關(guān)于場(chǎng)變量的多項(xiàng)式形式的局部能量密度，因此局部能量密度公式可導(dǎo)，局部能量密度公式關(guān)于兩種場(chǎng)變量的偏導(dǎo)數(shù)如式(1.3)所示。

其中，

根據(jù) Khachaturyan 理論，在傅里葉空間中，彈性應(yīng)變能的計(jì)算如式(1.4)的積分所示，相場(chǎng)方程中的彈性應(yīng)變能部分如式(1.5)所示，該項(xiàng)被表達(dá)為場(chǎng)變量c或的函數(shù)，引入相場(chǎng)方程中，這里以彈性應(yīng)變能為場(chǎng)變量c的函數(shù)進(jìn)行說(shuō)明，彈性應(yīng)變能為場(chǎng)變量的函數(shù)的情況與之類似。

其中，

B(n)是傅里葉空間的彈性相互作用能，是場(chǎng)變量c的結(jié)構(gòu)函數(shù)，與分子結(jié)構(gòu)相關(guān)。k是傅里葉空間的向量，單位向量表示傅里葉變換，表示傅里葉逆變換，在三維空間中，是一個(gè)由四維張量根據(jù)物質(zhì)結(jié)構(gòu)對(duì)稱性得到的二維剛度矩陣，i、j、k、l分別取1-3 的數(shù)值，是應(yīng)變矩陣，剛度矩陣和應(yīng)變矩陣決定了彈性應(yīng)變能的大小和類型。應(yīng)力矩陣

對(duì)總能量方程(1.1)進(jìn)行變分求解，可以得到描述守恒場(chǎng)變量c的Cahn-Hilliard 方程以及描述結(jié)構(gòu)場(chǎng)變量的Allen-Cahn 方程兩種方程組成的耦合方程。Cahn-Hilliard 方程和Allen-Cahn 方程的耦合方程如式(1.6)所示，其中Cahn-Hilliard 方程中的M是遷移率，而Allen-Cahn 方程中的L是擴(kuò)散系數(shù)。

2 指數(shù)時(shí)間差分方法的求解格式與能量穩(wěn)定的證明

2.1 指數(shù)時(shí)間差分方法的求解格式

下面對(duì)式(1.6)所示含彈性應(yīng)變能的各向異性相場(chǎng)方程的緊致指時(shí)間差分算法進(jìn)行推導(dǎo)。假設(shè)，模擬區(qū)域?yàn)榭臻g：

其中，

且

則有：

模擬區(qū)域采用周期邊界條件[9]。

為了保證能量穩(wěn)定需要進(jìn)行算子分裂，后文中對(duì)能量穩(wěn)定的算子分裂格式取值進(jìn)行了證明。這里直接引入化學(xué)自由能分裂算子和，彈性應(yīng)變能分裂算子，界面能各向異性分裂算子矩陣。算子分裂后的方程如式(2.1)所示。

其中，各項(xiàng)異性算子分裂的矩陣如下：

定義二維的方陣G 如下：

則，矩陣A,B,C可以表示為：

定義運(yùn)算：

空間離散形式的拉普拉斯算子矩陣可以寫(xiě)成[9]:

將算子代入式(2.1)，得到方程的空間離散形式如式(2.2)所示 ：

根據(jù)緊致指數(shù)時(shí)間差分方法的描述，這里將方程中含有高階導(dǎo)數(shù)的線性項(xiàng)與含有局部能量密度、界面能各向異性、彈性應(yīng)變能的非線性項(xiàng)g分離，得到式（2.3）。

其中，

在緊致指數(shù)時(shí)間差分方法中，通過(guò)離散傅里葉變換[12]來(lái)降低計(jì)算復(fù)雜度。將矩陣A，B，C根據(jù)特征值分解為如下形式：

其中，

定義算子：

代入式(2.3)，可得式（2.4），經(jīng)推導(dǎo)后可得式(2.5)。

其中，

在緊致指數(shù)時(shí)間差分方法中，式(2.5)需要與定義的如下指數(shù)時(shí)間因子相乘，

定義兩個(gè)矩陣的點(diǎn)乘運(yùn)算：

可得式(2.6)：

式(2.6)兩側(cè)均對(duì)時(shí)間積分可得到如下公式：

積分后，得到的場(chǎng)變量在第n+1 時(shí)間步和第n 時(shí)間步的迭代關(guān)系如式(2.7)所示。

根據(jù)拉格朗日多項(xiàng)式近似估計(jì)式(2.7)中關(guān)于非線性項(xiàng)的積分，得到一階緊致指數(shù)時(shí)間差分方法的求解格式(2.8)。

其中，

預(yù)估校正的緊致指數(shù)時(shí)間差分方法二階求解格式如下：

2.2 能量穩(wěn)定的證明

下面證明2.1 節(jié)中得到的含彈性應(yīng)變能的各項(xiàng)異性相場(chǎng)模型的緊致指數(shù)時(shí)間差分方法一階求解格式（2.8）和二階求解格式（2.9）的能量穩(wěn)定性。在證明過(guò)程中引用了如下引理：

Lemma 2.1對(duì)于任意給定的的三維矩陣f和g，有：

Lemma 2.2對(duì)于任意給定的的三維矩陣f和g， 如果則：

Lemma 2.3（Young’s Inequality）

下面是含彈性應(yīng)變能的各向異性耦合方程一階指數(shù)時(shí)間差分方法能量穩(wěn)定性定理：

Theorem 2.1對(duì)于含彈性應(yīng)變能的各向異性的Allen-Cahn 方程和Cahn-Hilliard 方程的耦合模型（1.6），如果滿足條件：

則一階緊致指數(shù)時(shí)間差分方法求解格式(2.8)滿足能量下降準(zhǔn)則。

Theorem 2.1的證明如下：

根據(jù)一階緊致指數(shù)時(shí)間差分方法公式（2.8），可以得到：

在方程兩側(cè)加上：

并將方程中的S 替換后可以得到：

其中，

方程兩側(cè)進(jìn)行逆變換可以得到：

方程兩側(cè)同時(shí)乘上：

并根據(jù)Lemma 2.1進(jìn)行離散積分，將兩個(gè)變量的方程合并整理后可得如下方程：

根據(jù)局部能量密度和彈性應(yīng)變能的二階泰勒展開(kāi)，得到系統(tǒng)第n個(gè)時(shí)間步和第n+1 個(gè)時(shí)間步的能量差的式兩種等價(jià)形式，式（2.10）和式（2.11）。

根據(jù)Lemma 2.2可知，式（2.11）中：

根據(jù)Lemma 2.3可得：

因此，若滿足Theorem2.1中描述的條件，則式（2.11）的右側(cè)大于0，即：

滿足能量下降準(zhǔn)則，Theorem 2.1得證。

下面是含彈性應(yīng)變能的各向異性耦合方程二階指數(shù)時(shí)間差分方法能量穩(wěn)定性定理：

Theorem 2.2對(duì)于含彈性應(yīng)變能的各向異性的Allen-Cahn 方程和Cahn-Hilliard 方程的耦合模型（1.6），如果滿足條件

則二階緊致指數(shù)時(shí)間差分方法求解格式(2.9)滿足能量下降準(zhǔn)則。

Theorem 2.2的證明如下：

根據(jù)二階緊致指數(shù)時(shí)間差分方法公式（2.9）可得：

在方程兩側(cè)加上：

并替換方程中的S可得：

其中r的取值和Theorem 2.1的證明中相同。

方程兩側(cè)進(jìn)行逆變換得到：

方程兩側(cè)乘

根據(jù)Lemma 2.1進(jìn)行離散積分，并根據(jù)Theorem 2.1的證明可得如下等式：

其中，

將式中的局部能量密度和彈性應(yīng)變能進(jìn)行二階泰勒展開(kāi)，整理后得到系統(tǒng)第n個(gè)時(shí)間步和第n+1個(gè)時(shí)間步的能量之差的兩種等價(jià)形式，如式（2.13）和式（2.14）所示：

根據(jù)Lemma 2.2可知，式（2.12）和式（2.14）中：

根據(jù)Lemma 2.3可得：

因此，若滿足Theorem2.2中描述的條件，則式（2.14）的右側(cè)大于0，即：

滿足能量下降準(zhǔn)則，Theorem 2.2得證。

2.3 指數(shù)時(shí)間差分方法的時(shí)間精度

下面通過(guò)局部截?cái)嗾`差分析指數(shù)時(shí)間差分方法一階和二階求解格式的精度。

為了求解第n 個(gè)時(shí)間步與第n+1 個(gè)時(shí)間步的局部截?cái)嗾`差，首先假設(shè)第n 個(gè)時(shí)間步的數(shù)值解與精確解相等：

一階求解格式的局部截?cái)嗾`差定義為：

將一階求解格式的數(shù)值解式(2.8)與精確解（2.7）作差并將積分展開(kāi)可得：

將S代入并進(jìn)行泰勒展開(kāi)可得：

因此，局部截?cái)嗾`差為：

將g代入得到：

指數(shù)時(shí)間差分方法一階求解格式的局部截?cái)嗾`差的主項(xiàng)為二階，因此該格式的時(shí)間精度為一階。

二階求解格式的局部截?cái)嗾`差定義為：

將二階求解格式（2.9）中的預(yù)估與校正公式作差可得：

根據(jù)式（2.15）可得預(yù)估步與精確解的差值：

將式（2.16）與式（2.17）相加并整理可得，二階格式的局部截?cái)嗾`差為：

將g和h代入公式可得：

指數(shù)時(shí)間差分方法二階求解格式的局部截?cái)嗾`差主項(xiàng)為三階，因此該格式的時(shí)間精度為二階。

3 數(shù)值實(shí)驗(yàn)

通過(guò)觀察鎳基合金數(shù)值實(shí)驗(yàn)和Zr 的氫化物堆疊數(shù)值實(shí)驗(yàn)兩個(gè)相場(chǎng)模型的數(shù)值模擬結(jié)果，對(duì)含有界面能各向異性和彈性應(yīng)變能的緊致指數(shù)時(shí)間差分方法的正確性以及算子分裂格式的能量穩(wěn)定性進(jìn)行驗(yàn)證。本文中的數(shù)值實(shí)驗(yàn)是在中國(guó)科學(xué)院計(jì)算機(jī)網(wǎng)絡(luò)信息中心的超級(jí)計(jì)算機(jī)“元”二期計(jì)算系統(tǒng)的 CPU節(jié)點(diǎn)進(jìn)行。以下數(shù)值實(shí)驗(yàn)均為三維模擬，為方便觀察，實(shí)驗(yàn)結(jié)果為模擬區(qū)域中心截面截圖。

3.1 鎳基合金的數(shù)值實(shí)驗(yàn)

鎳基合金的數(shù)值實(shí)驗(yàn)[13]模擬了高溫下的單一無(wú)序相分解成低溫下的富溶質(zhì)的有序相和溶質(zhì)貧化的無(wú)序相的過(guò)程。該數(shù)值實(shí)驗(yàn)的相場(chǎng)模型不含界面能各向異性，因此我們關(guān)注的重點(diǎn)是彈性應(yīng)變能的作用。其相場(chǎng)耦合方程如式(3.1)所示。

其中，

數(shù)值實(shí)驗(yàn)中的參數(shù)均采用無(wú)量綱化后的取值。其中，A1= 62.5，A2= 25，A3= 12，A4= 6.25，C'= 0.2，C''= 0.6，遷移率M為1.0，擴(kuò)散系數(shù)L為1.0，界面梯度能量的系數(shù)的取值為：

彈性應(yīng)變能計(jì)算相關(guān)的彈性應(yīng)變矩陣和剛度矩陣如下：

模擬的時(shí)間步長(zhǎng)為1e-2，模擬的總時(shí)間為20.0，網(wǎng)格尺寸為1.0，三維的模擬區(qū)域大小為128×128×128。在模擬的初始狀態(tài)，模擬區(qū)域中央給定一個(gè)半徑為8 個(gè)網(wǎng)格寬度的新相粒子，場(chǎng)變量c的初始值為0.4 和0.65，場(chǎng)變量的取值為0 和1.4。[14]算子分裂參數(shù)取值如下：

圖1 不含彈性應(yīng)變能的鎳基合金數(shù)值實(shí)驗(yàn)的模擬結(jié)果Fig.1 Simulation results of numerical experiment of Ni-based alloy without elastic strain energy

圖2 含有彈性應(yīng)變能的鎳基合金數(shù)值實(shí)驗(yàn)?zāi)M結(jié)果Fig.2 Simulation results of numerical experiment of Ni-based alloy with elastic strain energy

以下分別進(jìn)行了無(wú)彈性應(yīng)變能和有彈性應(yīng)變能的模擬。無(wú)彈性應(yīng)變能的場(chǎng)變量c在迭代步數(shù)t=0,1000,2000 的模擬結(jié)果如（圖1）所示，沒(méi)有彈性應(yīng)變能的作用，粒子隨著時(shí)間生長(zhǎng)，演化成了球形；有彈性應(yīng)變能的場(chǎng)變量c在迭代步數(shù) t=0,1000,2000的模擬結(jié)果如（圖2）所示。在彈性應(yīng)變能的作用下，粒子隨時(shí)間生長(zhǎng)，演化成了方形。其能量下降曲線如（圖3）所示，能量穩(wěn)定。

圖3 鎳基合金數(shù)值實(shí)驗(yàn)的能量下降曲線Fig.3 Energy decline curve of numerical experiment of Nibased alloy

3.2 Zr 的氫化物堆疊數(shù)值實(shí)驗(yàn)

Zr 的氫化物堆疊數(shù)值實(shí)驗(yàn)[8]是描述氫化物在Zr 中的堆積結(jié)構(gòu)形成和轉(zhuǎn)變的三維相場(chǎng)模型。該數(shù)值實(shí)驗(yàn)既含有界面能各向異性，又含有彈性應(yīng)變能。其耦合方程如式(3.2)所示：

其中，

w為1.848e9，遷移率M為5e-4，擴(kuò)散系數(shù)D為1.2302e-10，界面梯度能量的系數(shù)的取值為：

界面能各向異性矩陣取值如下：

彈性應(yīng)變能計(jì)算相關(guān)的彈性應(yīng)變矩陣和剛度矩陣為：

模擬的時(shí)間步長(zhǎng)為1e-7，模擬的總時(shí)間為2e-4，網(wǎng)格尺寸為1e-9，三維的模擬區(qū)域大小為128×128×128。在模擬的初始狀態(tài)，模擬區(qū)域中央給定一個(gè)半徑為10 個(gè)網(wǎng)格寬度的新相粒子，場(chǎng)變量c的初始值為0.1 和0.5982，場(chǎng)變量的取值為0 和1.0。

緊致指數(shù)時(shí)間差分方法算子分裂參數(shù)取值為：

下面分別進(jìn)行了有無(wú)界面能各向異性和有無(wú)彈性應(yīng)變能的四組數(shù)值實(shí)驗(yàn)。無(wú)界面能各向異性和彈性應(yīng)變能，模擬結(jié)果如（圖4）所示，氫化物逐漸演化為球體；有彈性應(yīng)變能，無(wú)界面能各向異性，如（圖5）所示，氫化物發(fā)生旋轉(zhuǎn)，成為橢球體；無(wú)彈性應(yīng)變能，有界面能各向異性，模擬結(jié)果如（圖6）所示，氫化物在界面能各向異性的作用下變成圓盤(pán)；有彈性應(yīng)變能和界面能各向異性，模擬結(jié)果如（圖7）所示，氫化物在界面能各向異性的作用下變成圓盤(pán)，在彈性應(yīng)變能的作用下發(fā)生旋轉(zhuǎn)，但是在兩者的共同作用下，氫化物旋轉(zhuǎn)的角度小于只有彈性應(yīng)變能作用的情況下的旋轉(zhuǎn)角度。其能量下降曲線如（圖8）所示，能量穩(wěn)定。

圖4 無(wú)界面能各向異性、無(wú)彈性應(yīng)變能的Zr 的氫化物數(shù)值實(shí)驗(yàn)的演化過(guò)程Fig.4 Evolution process of numerical experiment of Zr-hydride without interface energy anisotropy and elastic strain energy

圖6 有界面能各向異性、無(wú)彈性應(yīng)變能的Zr 的氫化物數(shù)值實(shí)驗(yàn)的演化過(guò)程Fig.6 Evolution process of numerical experiment of Zr-hydride with interface energy anisotropy and without elastic strain energy

圖7 有界面能各向異性，有彈性應(yīng)變能的Zr 的氫化物數(shù)值實(shí)驗(yàn)的演化過(guò)程Fig.7 Evolution process of numerical experiment of Zr- hydride with interface energy anisotropy and elastic strain energy

圖8 Zr 的氫化物數(shù)值實(shí)驗(yàn)的能量下降曲線Fig.8 Energy decline curve of numerical experiment of Zrhydride

下面進(jìn)行了一系列時(shí)間步長(zhǎng)倍增的實(shí)驗(yàn)用以驗(yàn)證指數(shù)時(shí)間差分方法的時(shí)間精度。時(shí)間步長(zhǎng)分別取（其中δ =4e-9），以相同初始值運(yùn)行至同一時(shí)刻T= 2.56e-5，將時(shí)間步長(zhǎng)dt = 1e-9 作為基準(zhǔn)解。計(jì)算了指數(shù)時(shí)間差分方法一階格式的L2誤差和收斂率如表1所示，收斂率的值接近1；二階格式的L2誤差和收斂率如表2所示，收斂率的值接近2，驗(yàn)證了求解格式的時(shí)間精度。

表1 指數(shù)時(shí)間差分方法一階求解格式的誤差Table 1 The error of the first order scheme of the Exponential Time Difference method

表2 指數(shù)時(shí)間差分方法二階求解格式的誤差Table 2 The error of the second order scheme of the Exponential Time Difference method

4 結(jié)論與展望

本文在緊致指數(shù)時(shí)間差分方法框架下引入了界面能各向異性和彈性應(yīng)變能的計(jì)算，將界面能各向異性和彈性應(yīng)變能歸于緊致指數(shù)時(shí)間差分方法的非線性項(xiàng)統(tǒng)一處理，通過(guò)積分和預(yù)估校正的二階近似進(jìn)行計(jì)算，為保證能量穩(wěn)定，設(shè)計(jì)了界面能各向異性和彈性應(yīng)變能的算子分裂格式，并對(duì)算子分裂格式的能量穩(wěn)定性進(jìn)行了數(shù)學(xué)證明。通過(guò)材料腐蝕相場(chǎng)模型的數(shù)值實(shí)驗(yàn)，觀察了界面能各向異性和彈性應(yīng)變能對(duì)場(chǎng)變量演化的影響，驗(yàn)證了含彈性應(yīng)變能的各向異性相場(chǎng)模型的指數(shù)時(shí)間差分方法算子分裂格式的能量穩(wěn)定性。本文使用超級(jí)計(jì)算系統(tǒng)“元”對(duì)指數(shù)時(shí)間差分方法的格式精度和能量穩(wěn)定性進(jìn)行了驗(yàn)證，在后續(xù)的工作中，擬使用“天河二號(hào)”超級(jí)計(jì)算機(jī)進(jìn)行并行擴(kuò)展性的測(cè)試。[15-16]本文僅得到了指數(shù)時(shí)間差分方法的一階和二階求解格式，更高階的求解格式及其穩(wěn)定性有待進(jìn)一步探索。

利益沖突聲明

所有作者聲明不存在利益沖突關(guān)系。