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

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

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

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

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

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

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

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

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

1.1 RT方法

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

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

Stot×CNDEF(j,t)

(1)

其中：CNDEF(j,t)為t時(shí)刻尺寸為j的NDEF類型缺陷的濃度；GNDEF(j)為級(jí)聯(lián)碰撞過(guò)程產(chǎn)生尺寸為j的NDEF類型缺陷的速率；w(k,j)為缺陷尺寸由k變?yōu)閖的轉(zhuǎn)化速率；Stot為材料內(nèi)固有缺陷對(duì)尺寸為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)部固有缺陷對(duì)尺寸j的團(tuán)簇的吸收項(xiàng)。

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

1.2 ETD算法

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

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

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

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

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

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

(2)

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

(3)

(4)

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

(5)

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

(6)

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

(7)

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

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

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

(8)

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

ETD1格式：

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

(9)

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

(10)

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

(11)

其中：Cvac為單尺寸空位原子在t時(shí)刻的濃度；Cn為尺寸為n的團(tuán)簇在t時(shí)刻的濃度；n為團(tuán)簇尺寸；N為團(tuán)簇的總體尺寸；F(t)矩陣的1～N行分別對(duì)應(yīng)非線性算子對(duì)尺寸為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ù)值求解時(shí)的精度控制。

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

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

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

為驗(yàn)證算法穩(wěn)定性，對(duì)基于缺陷數(shù)密度的守恒量進(jìn)行數(shù)值驗(yàn)證。在熱老化條件下，雖然空位能聚集形成空位團(tuán)簇，但沒(méi)有外來(lái)缺陷的影響，總的空位數(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時(shí)空位團(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為初始時(shí)間；tm為第m個(gè)時(shí)間步，在演化過(guò)程中，空位數(shù)量值保持穩(wěn)定，且與初始狀態(tài)空位數(shù)相同。

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

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

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

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

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

(25)

其中：

(26)

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

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

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

(27)

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

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

(28)

對(duì)方程(27)的求解等價(jià)于求解下面的方程：

(29)

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

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

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

(30)

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

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

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

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

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

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

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

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

表3 中子輻照SA-304 SS位錯(cuò)環(huán)演化的團(tuán)簇動(dòng)力學(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)簇動(dòng)力學(xué)方法[22]的模擬結(jié)果相吻合，且與實(shí)驗(yàn)觀察數(shù)據(jù)的變化趨勢(shì)相符，在高損傷劑量下，模擬結(jié)果略高于實(shí)驗(yàn)結(jié)果，但二者沒(méi)有量級(jí)的差別。利用ETD方法求解的計(jì)算數(shù)密度結(jié)果與304鋼的實(shí)驗(yàn)值較為接近，數(shù)值求解格式解得的半徑分布與實(shí)驗(yàn)結(jié)果非常接近，從而驗(yàn)證了ETD方法在點(diǎn)缺陷可動(dòng)模型求解中的可行性與有效性。

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

3 總結(jié)與展望

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

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