辛之夼，聶寧明，賀新福，王彥棡，吳 石，王 玨

(1.中國科學(xué)院 計(jì)算機(jī)網(wǎng)絡(luò)信息中心，北京 100190；2.中國科學(xué)院大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院，北京 100049；3.中國原子能科學(xué)研究院，北京 102413)

材料輻照損傷一直是材料計(jì)算重要的研究領(lǐng)域。材料在反應(yīng)堆中經(jīng)受高能中子轟擊產(chǎn)生過飽和缺陷，輻照缺陷在溫度場、應(yīng)力場等作用下經(jīng)過長時間演化，導(dǎo)致材料微結(jié)構(gòu)演化乃至宏觀性能退化，是影響反應(yīng)堆安全性和經(jīng)濟(jì)性的關(guān)鍵因素之一。隨著物理理論研究的深入和高性能計(jì)算技術(shù)的發(fā)展，研究人員采用實(shí)驗(yàn)加先進(jìn)建模再加計(jì)算模擬的手段對材料輻照損傷的復(fù)雜機(jī)理開展研究。

速率理論[1](RT)是一種用于模擬材料輻照損傷微觀結(jié)構(gòu)長時間演化的重要方法，相比較分子動力學(xué)(MD)、動力學(xué)蒙特卡羅(KMC)[2]這類直接模擬演化軌跡的方法，RT方法基于平均場近似，忽略了缺陷演化的空間信息，在均勻介質(zhì)中采用反應(yīng)速率表征缺陷間相互作用的快慢，從而實(shí)現(xiàn)對微觀結(jié)構(gòu)演化行為的高效模擬。RT方法模擬相較于模擬金屬中輻照缺陷隨時間演化的常用方法——實(shí)體動力學(xué)蒙特卡羅方法[3](OKMC)而言，其能模擬約100 dpa下各種類型缺陷的演化過程，遠(yuǎn)強(qiáng)于OKMC方法，因此在模擬高損傷劑量材料微觀結(jié)構(gòu)演化時，能明顯體現(xiàn)出優(yōu)勢。而且RT在介觀層面上進(jìn)行模擬，是建立材料宏觀性能與微觀組織之間關(guān)聯(lián)關(guān)系的重要手段。近30年來，國際上采用RT方法成功模擬了反應(yīng)堆結(jié)構(gòu)材料內(nèi)輻照誘導(dǎo)各類微觀缺陷隨損傷劑量的演化行為，如位錯環(huán)[4-5]、析出物[6]等，并用于預(yù)測輻照脆化[7]、輻照腫脹[8]等宏觀性能。

RT方法的應(yīng)用難點(diǎn)主要在于對剛性極強(qiáng)的大規(guī)模RT方程組進(jìn)行數(shù)值求解。受限于計(jì)算能力，RT方程求解方法大致可歸為確定性方法和隨機(jī)性方法兩類，這兩類方法的基本思想均是采用近似方法對方程組進(jìn)行方程數(shù)量縮減，以減少計(jì)算量。對于確定性方法，一種是按團(tuán)簇尺寸進(jìn)行分組求解的分組方法(Group-Method)[9]：用“平均方程”來代替方程組，能有效地減少欲求解方程的數(shù)量，降低了計(jì)算復(fù)雜度，但分組需要技巧。另一種為 Fokker-Planck方法[10]：對于小尺寸團(tuán)簇用主方程求解獲得離散的解，對于大尺寸團(tuán)簇將其轉(zhuǎn)變?yōu)镕okker-Planck方程進(jìn)行連續(xù)化近似求解。該方法對于求解二元體系(間隙-空位)很有效，但對于有3種或以上類型缺陷體系的求解很受限制。對于隨機(jī)性方法，其中比較經(jīng)典的是隨機(jī)團(tuán)簇動力學(xué)方法(SCD)[11]。SCD方法是隨機(jī)模擬算法(stochastic simulation arithmetic,SSA)[12]在核材料模擬中的一個擴(kuò)展，利用隨機(jī)方法求解方程組，減少了欲求解方程的數(shù)量，從而降低了求解難度。兩種方法的目的均是為了在保證準(zhǔn)確度的同時，減少計(jì)算量，提高求解效率。

本文引入指數(shù)時間差分(ETD)方法對大規(guī)模的RT方程組進(jìn)行并行求解。ETD方法[13]是一類用于求解剛性常微分方程組的高精度數(shù)值算法，與傳統(tǒng)的半隱式歐拉法相比，ETD方法均能達(dá)到對應(yīng)的最優(yōu)時間收斂率且表現(xiàn)優(yōu)于傳統(tǒng)的半隱式歐拉法[14]。ETD方法具有良好的穩(wěn)定性[15]，其對非線性方程的線性部分進(jìn)行精確求解，對非線性部分進(jìn)行插值近似，在保證格式精度的同時，可很好地解決由于穩(wěn)定條件的限制導(dǎo)致時間步長過小的問題。利用ETD方法的數(shù)值計(jì)算優(yōu)勢對RT主方程直接求解，可得到更符合RT方程的數(shù)值解，從而對材料輻照損傷模擬有更好的模擬效果。

本文針對模擬輻照損傷的RT方法和ETD方法進(jìn)行介紹，分析ETD方法在求解RT主方程中使用的可能性，并介紹ETD方法在空位團(tuán)簇演化模型和點(diǎn)缺陷可動模型兩類RT方程求解中的不同計(jì)算實(shí)現(xiàn)。

1 RT模型和數(shù)值求解算法介紹

輻照導(dǎo)致材料內(nèi)微觀缺陷的形核、長大行為是一個跨越多個時間(從缺陷產(chǎn)生約10-15s到微結(jié)構(gòu)演化約s/h/a)、空間尺度(點(diǎn)缺陷約0.1 nm到微結(jié)構(gòu)約μm)的演化過程，描述微觀結(jié)構(gòu)演化的MD、KMC等方法受計(jì)算資源的限制，往往無法模擬微結(jié)構(gòu)的長時間演化行為，而RT忽略了體系的隨機(jī)效應(yīng)和空間的關(guān)聯(lián)性，因此在模擬缺陷演化的過程中不會受到時空尺度的限制，其模擬結(jié)果能與微觀表征結(jié)果(如TEM、SANS等)進(jìn)行直接對比。因此RT在輻照導(dǎo)致材料內(nèi)微觀缺陷演化過程的研究中得到了非常廣泛的應(yīng)用。

1.1 RT方法

RT方法類似于化學(xué)反應(yīng)過程的求解，將缺陷之間的相互作用簡化為一組微分方程來描述。在輻照條件下，材料中同時存在多種缺陷，包括點(diǎn)缺陷、缺陷團(tuán)簇以及位錯、晶界等缺陷阱。缺陷的產(chǎn)生、聚集、釋放以及阱吸收等相互作用均會在缺陷演化RT主方程中予以考慮。RT的一般形式如式(1)所示。

CNDEF(k,t)-∑w(j,l)×CNDEF(j,t)-

Stot×CNDEF(j,t)

(1)

其中：CNDEF(j,t)為t時刻尺寸為j的NDEF類型缺陷的濃度；GNDEF(j)為級聯(lián)碰撞過程產(chǎn)生尺寸為j的NDEF類型缺陷的速率；w(k,j)為缺陷尺寸由k變?yōu)閖的轉(zhuǎn)化速率；Stot為材料內(nèi)固有缺陷對尺寸為j的缺陷的吸收速率；主方程右側(cè)的第2項(xiàng)表示由尺寸k的缺陷團(tuán)簇生成尺寸j的團(tuán)簇的反應(yīng)速率，右側(cè)第3項(xiàng)表示由尺寸j的團(tuán)簇生成其他尺寸的團(tuán)簇的生成速率，右側(cè)第4項(xiàng)表示材料內(nèi)部固有缺陷對尺寸j的團(tuán)簇的吸收項(xiàng)。

RT方法通過RT主方程(式(1))描述團(tuán)簇?cái)?shù)密度隨時間變化的演化過程。方程維數(shù)與缺陷數(shù)量相關(guān)。伴隨損傷劑量增加，位錯環(huán)、空洞尺寸逐漸增大，預(yù)測高劑量下位錯環(huán)等缺陷的演化行為，通常需要百萬量級的方程數(shù)量。伴隨著合金成分的增加，方程數(shù)量呈幾何型增長。且這一類方程通常是剛性病態(tài)的微分方程組，對數(shù)值求解提出了較高要求。因此，對該大規(guī)模微分方程組的高效精確的數(shù)值求解是其主要應(yīng)用難點(diǎn)。

1.2 ETD算法

ETD算法是由Cox和Matthews提出的一類用于數(shù)值求解剛性微分方程的時間積分方法[16-17]。由于隱式算法格式通常具有較高的計(jì)算穩(wěn)定性，剛性微分方程的數(shù)值求解通常采用隱式求解方法。但隱式算法格式在計(jì)算過程中需對方程進(jìn)行迭代求解，計(jì)算量較大，對于長時間模擬不友好。在ETD方法的求解過程中，首先將剛性方程組分解成線性和非線性部分，對于被積函數(shù)中的線性部分直接進(jìn)行精確計(jì)算，對非線性部分采用插值逼近的方式進(jìn)行求解，這樣可得到非線性部分的高精度積分。所以ETD方法可在保證精度的情況下，對剛性微分方程進(jìn)行顯式地穩(wěn)定求解，同時相較于隱式算法格式大幅降低了計(jì)算量。Cox和Matthews給出了ETD-RK格式[16]，應(yīng)用到具有快速衰變特征的剛性方程的求解上。ETD還被應(yīng)用于相場方程求解[18-19]，結(jié)果表明與傳統(tǒng)相場求解方法相比，求解效果與使用半隱式歐拉法的求解效果相比具有更好的收斂精度。

2 ETD方法應(yīng)用于RT模型求解的典型算例

基于上述優(yōu)勢，采用ETD方法對Ni中空位團(tuán)簇演化以及奧氏體不銹鋼中微結(jié)構(gòu)(位錯環(huán)以及空洞)演化過程建立RT模型并進(jìn)行數(shù)值求解，并詳細(xì)介紹整個應(yīng)用過程。

2.1 空位團(tuán)簇演化模型

1) 模型建模與算法格式推導(dǎo)

當(dāng)金屬材料內(nèi)某個原子逃離原有晶體點(diǎn)陣位置，則在該位置出現(xiàn)空缺，并導(dǎo)致空位形成，在一定溫度條件下，金屬材料內(nèi)的空位處于熱平衡狀態(tài)[20]。長時間熱老化服役條件下，材料內(nèi)的空位會在熱力學(xué)驅(qū)動下發(fā)生擴(kuò)散，并在遷移過程中不斷發(fā)生聚集、釋放反應(yīng)。空位間的聚集、俘獲反應(yīng)導(dǎo)致空位團(tuán)簇不斷長大，而空位的釋放則引起空位團(tuán)縮小，伴隨服役時間延長，材料內(nèi)的空位型缺陷不斷發(fā)生相互作用，從而引起各尺寸空位缺陷濃度隨時間發(fā)生變化，其演化方程形式如式(2)所示。

(2)

其中：Cvac為單個空位的濃度；Cn為尺寸為n的空位團(tuán)簇的數(shù)密度，n=2,…,N；βn為吸收速率；αn為發(fā)射速率，可通過下式計(jì)算：

(3)

(4)

由于模擬過程中沒有新的空位點(diǎn)缺陷產(chǎn)生，所以要求保持物質(zhì)總量Qtot守恒：

(5)

初始條件通過下式計(jì)算：

(6)

為后續(xù)ETD格式的推導(dǎo)，上述模型可表示為如下矩陣形式：

(7)

其中：C(t)為t時刻的缺陷濃度，C(t)=(Cvac,C2,…,Cn)；右邊第1項(xiàng)為線性部分，第2項(xiàng)F為非線性部分，表示非線性算子對C(t)的作用。

首先將上式乘以積分因子e-tL，然后在1個時間步(從t=tk到tk+1=tk+Δt)上對方程進(jìn)行積分，得到如下形式：

C(tk+1)=eΔt·LC(tk)+

(8)

其中：k為迭代的時間步數(shù)；Δt為時間步長。采用數(shù)值積分求解式中積分項(xiàng)。右端的積分項(xiàng)采用不同數(shù)值近似算法可得到不同的ETD格式。

ETD1格式：

C(t+Δt)=C(t)·eΔt·L+

(9)

其中，I為單位對角陣。在空位團(tuán)簇演化模型求解中，L由下式的對角陣給出：

(10)

F(t)由下式計(jì)算：

(11)

其中：Cvac為單尺寸空位原子在t時刻的濃度；Cn為尺寸為n的團(tuán)簇在t時刻的濃度；n為團(tuán)簇尺寸；N為團(tuán)簇的總體尺寸；F(t)矩陣的1～N行分別對應(yīng)非線性算子對尺寸為1～N的團(tuán)簇的作用。

ETD4RK格式：

(12)

L由下式計(jì)算：

(13)

poly1、poly2、poly3由下式計(jì)算：

(14)

(15)

(16)

F(t)由下式計(jì)算：

(17)

a(t)、b(t)、c(t)由下式計(jì)算：

a(t)=φ1Ct+φ2F(C(t),t)

(18)

(19)

C(t)=φ1a(t)+

(20)

φ1、φ2、φ3由下式計(jì)算：

(21)

(22)

(23)

其中，γ為調(diào)整系數(shù)，用于數(shù)值求解時的精度控制。

2) 空位團(tuán)簇演化模型算例結(jié)果分析

基于文獻(xiàn)[9]中的物理模型與參數(shù)，本文嘗試采用上一節(jié)中給出的ETD算法對速率方程進(jìn)行求解驗(yàn)證。模型參數(shù)列于表1。

采用ETD算法對空位及其團(tuán)簇演化的RT主方程進(jìn)行求解，并將1階的ETD1和4階ETD4RK格式計(jì)算結(jié)果與Group-Method方法的計(jì)算結(jié)果[9]進(jìn)行對比，兩種方法的對比結(jié)果如圖1所示，圖中示出了金屬鎳及鎳基合金熱老化103s條件下的空位團(tuán)簇尺寸分布。采用1階、4階ETD算法的計(jì)算結(jié)果與Group-Method方法的計(jì)算結(jié)果基本相近，證明了ETD算法在求解空位團(tuán)簇演化RT模型方程中的可行性與有效性。

為驗(yàn)證算法穩(wěn)定性，對基于缺陷數(shù)密度的守恒量進(jìn)行數(shù)值驗(yàn)證。在熱老化條件下，雖然空位能聚集形成空位團(tuán)簇，但沒有外來缺陷的影響，總的空位數(shù)量始終保持守恒，其守恒量的形式可表達(dá)為：

表1 金屬鎳及鎳基合金熱老化速率方程模型參數(shù)Table 1 Model parameters of thermal aging rate equation for nickel metal and nickel-based alloys

圖1 金屬鎳及鎳基合金823 K熱老化103 s時空位團(tuán)簇尺寸分布采用不同算法的計(jì)算模擬結(jié)果Fig.1 Simulated results of vacancy cluster size distribution of nickel metal and nickel-based alloys during thermal aging at 823 K for 103 s using different algorithms

(24)

其中：Cn為缺陷數(shù)密度；n為缺陷尺寸；t0為初始時間；tm為第m個時間步，在演化過程中，空位數(shù)量值保持穩(wěn)定，且與初始狀態(tài)空位數(shù)相同。

經(jīng)計(jì)算，初始濃度設(shè)置為1.0×10-7的條件下，守恒量為1.054 7×10-7。在106次迭代后，守恒量如圖2所示，仍然保持在這個值附近。在迭代計(jì)算的前100步左右，守恒量未立即達(dá)到1.054 7×10-7，這是由于初期計(jì)算缺陷數(shù)密度時，缺陷數(shù)密度計(jì)算結(jié)果基本均小于機(jī)器精度。隨著計(jì)算的進(jìn)行，缺陷數(shù)密度計(jì)算結(jié)果逐漸增大，最后守恒量數(shù)值逐漸增加最后趨于穩(wěn)定。

2.2 點(diǎn)缺陷可動模型

1) 模型建模與算法格式推導(dǎo)

點(diǎn)缺陷可動模型采用Poker于2004年針對中子輻照奧氏體不銹鋼位錯環(huán)演化行為建立的團(tuán)簇動力學(xué)模型。在壓水堆服役環(huán)境下，堆內(nèi)構(gòu)件用奧氏體不銹鋼內(nèi)的微觀缺陷以間隙型位錯環(huán)為主，這也成為其輻照促進(jìn)腐蝕開裂的1個主要誘因。為模擬預(yù)測奧氏體不銹鋼內(nèi)位錯環(huán)的演化行為，國內(nèi)外先后建立了描述輻照誘導(dǎo)奧氏體不銹鋼內(nèi)位錯環(huán)演化行為的團(tuán)簇動力學(xué)模型，中國原子能科學(xué)研究院開發(fā)的團(tuán)簇動力學(xué)RADIEFF軟件[21]，已廣泛應(yīng)用于模擬預(yù)測不同輻照條件下奧氏體不銹鋼內(nèi)位錯環(huán)的演化行為。與上述空位團(tuán)簇演化模型不同，點(diǎn)缺陷可動模型需考慮間隙型缺陷的演化行為以及間隙、空位型缺陷之間的相互作用，其具體方程形式如式(25)所示。

圖2 守恒量驗(yàn)證結(jié)果Fig.2 Verification result of conserved quantity

(25)

其中：

(26)

方程中所涉及的缺陷反應(yīng)類型包括缺陷的產(chǎn)生、缺陷之間的復(fù)合、吸收、釋放等，點(diǎn)缺陷可動速率方程模型涉及的缺陷反應(yīng)列于表2。

表2 點(diǎn)缺陷可動速率方程模型涉及的缺陷反應(yīng)Table 2 Defect reaction involved in point defect moveable rate theory model

對于點(diǎn)缺陷模型方程(式(25))這樣復(fù)雜的非線性方程組，采用ETD方法進(jìn)行數(shù)值求解，對方程組進(jìn)行時間上的離散處理以獲得顯式的時間步進(jìn)迭代公式。首先對方程(25)進(jìn)行處理。為表述簡便，采用矩陣形式(式(27))進(jìn)行ETD算法格式的描述。

(27)

將方程組(25)的右端項(xiàng)劃分為非線性項(xiàng)和線性項(xiàng)兩個部分。線性項(xiàng)矩陣L選取形式如式(28)所示，對于其余的部分組合為非線性項(xiàng)部分Fc。

L矩陣僅選取了方程組(25)右端項(xiàng)線性部分系數(shù)的主對角線元素和副對角線元素，而未選取其他線性項(xiàng)系數(shù)。是因?yàn)榻?jīng)過實(shí)驗(yàn)發(fā)現(xiàn)，不論是在L矩陣中添加其他線性項(xiàng)的系數(shù)還是在L矩陣直接取為方程的Jacobi矩陣，均未對計(jì)算結(jié)果精度有更明顯的改進(jìn)。綜合考慮計(jì)算精度和計(jì)算效率，本文選取了形式如式(28)的L矩陣進(jìn)行ETD算法格式的構(gòu)造。

(28)

對方程(27)的求解等價于求解下面的方程：

(29)

對式(29)進(jìn)行積分后再進(jìn)行數(shù)值逼近，推導(dǎo)過程與空位團(tuán)簇模型中1階算法格式推導(dǎo)類似，得到1階的ETD格式的迭代式(式(9))。

為提高計(jì)算結(jié)果的準(zhǔn)確度，采用預(yù)估校正法對每步時間步所得的結(jié)果進(jìn)行修正，設(shè)計(jì)2階的預(yù)估校正的ETD算法格式。使用式(9)求得初步的近似值Fc0，被稱為預(yù)測值。此時預(yù)測的精度可能很差，所以再使用1次式(9)，對其進(jìn)行1次矯正，得到校正值Fc1。

Cn+1=Cn·eΔt·L+(Fc1-Fc0)·

(30)

由于團(tuán)簇尺寸在百萬量級，即方程組數(shù)量N=106，直接的數(shù)值求解在單進(jìn)程上計(jì)算效率很低甚至無法計(jì)算。本文對上述求解格式進(jìn)行并行求解，將整個方程分解為M個部分，利用M個進(jìn)程，每個進(jìn)程負(fù)責(zé)1個子區(qū)域的計(jì)算，達(dá)到子區(qū)域之間并行計(jì)算的目的，則有N=n1,n2,n3,…,nM，如圖3所示。

并行劃分后，模擬區(qū)域的邊界是原始方程對應(yīng)的邊界，而子區(qū)域劃分產(chǎn)生的相鄰子區(qū)域之間的邊界為內(nèi)部邊界。考慮邊界條件的計(jì)算，求解時需分為0號進(jìn)程邊界和其他進(jìn)程邊界。

圖3 點(diǎn)缺陷可動模型的并行劃分Fig.3 Parallel partitioning of point defect movable model

將劃分到每個進(jìn)程上的部分方程分為控制域和計(jì)算域兩個部分，控制域組合成原始方程的解，計(jì)算域中通過重疊計(jì)算的方式通信數(shù)據(jù)，重疊方式如圖4所示。由于并行劃分，在劃分處會導(dǎo)致數(shù)值的丟失。實(shí)現(xiàn)過程中采用重疊計(jì)算的方式來降低并行劃分時數(shù)值丟失對結(jié)果的影響，通過通信邊界值來盡可能地減少數(shù)值的丟失。

圖4 重疊方式示意Fig.4 Overlapping method illustration

2)點(diǎn)缺陷可動模型算例結(jié)果分析

對點(diǎn)缺陷可動模型進(jìn)行ETD算法求解的算例計(jì)算。針對輻照誘導(dǎo)奧氏體不銹鋼內(nèi)位錯環(huán)演化行為建立團(tuán)簇動力學(xué)模型涉及大量模型參數(shù)，包括輻照參數(shù)、材料微觀參數(shù)等，這些參數(shù)通常來源于低尺度的模擬計(jì)算或由實(shí)驗(yàn)數(shù)據(jù)獲得的結(jié)果。為對比ETD算法的有效性，本模型所選參數(shù)與Poker等[22]所給參數(shù)一致(表3)。

采用預(yù)估校正ETD方法的計(jì)算結(jié)果、團(tuán)簇動力學(xué)程序RADIEFF模擬結(jié)果以及實(shí)驗(yàn)結(jié)果對比如圖5、6所示。

表3 中子輻照SA-304 SS位錯環(huán)演化的團(tuán)簇動力學(xué)模型所需參數(shù)Table 3 Parameter required for cluster dynamics model of neutron-irradiated SA-304 SS dislocation ring evolution

圖5 損傷劑量與直徑之間的關(guān)系Fig.5 Relationship between irradiation damage dose and diameter

由圖5、6可知，采用ETD方法獲取的計(jì)算結(jié)果與經(jīng)典團(tuán)簇動力學(xué)方法[22]的模擬結(jié)果相吻合，且與實(shí)驗(yàn)觀察數(shù)據(jù)的變化趨勢相符，在高損傷劑量下，模擬結(jié)果略高于實(shí)驗(yàn)結(jié)果，但二者沒有量級的差別。利用ETD方法求解的計(jì)算數(shù)密度結(jié)果與304鋼的實(shí)驗(yàn)值較為接近，數(shù)值求解格式解得的半徑分布與實(shí)驗(yàn)結(jié)果非常接近，從而驗(yàn)證了ETD方法在點(diǎn)缺陷可動模型求解中的可行性與有效性。

圖6 位錯環(huán)數(shù)密度隨輻照損傷劑量的變化 Fig.6 Variation of dislocation ring number density with irradiation damage dose

3 總結(jié)與展望

本文對RT方法的原理和建模以及ETD方法進(jìn)行了介紹。通過對空位團(tuán)簇演化模型和點(diǎn)缺陷可動模型兩類RT模型求解的計(jì)算實(shí)現(xiàn)，對ETD方法在求解模擬輻照損傷模型中的應(yīng)用進(jìn)行了詳細(xì)介紹，計(jì)算結(jié)果顯示了ETD方法在求解RT方程中的有效性?？傊褂肊TD方法求解RT模型的優(yōu)勢，主要體現(xiàn)在，ETD方法能對RT主方程進(jìn)行高精度的數(shù)值求解，減少了額外簡化方程所帶來的誤差，而并行求解的引入使得對于大規(guī)模RT方程組的直接數(shù)值求解成為可能?？赏ㄟ^直接對主方程進(jìn)行求解的方式，減少額外簡化方程所帶來的誤差，從而更好地對輻照損傷進(jìn)行模擬。

但ETD方法在應(yīng)用實(shí)現(xiàn)中仍存在不足。ETD計(jì)算過程中時間步長和不同精度算法的選擇需隨模型方程剛性變化而變化。為更加高效地實(shí)現(xiàn)ETD算法在速率方程高精度穩(wěn)定求解，自適應(yīng)的變階變步長ETD算法的實(shí)現(xiàn)需在下一步工作中繼續(xù)開展研究。