指數時間差分方法在材料輻照損傷模擬計算中的應用

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

(1.中國科學院 計算機網絡信息中心，北京 100190；2.中國科學院大學 計算機科學與技術學院，北京 100049；3.中國原子能科學研究院，北京 102413)

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

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

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

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

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

1 RT模型和數值求解算法介紹

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

1.1 RT方法

RT方法類似于化學反應過程的求解，將缺陷之間的相互作用簡化為一組微分方程來描述。在輻照條件下，材料中同時存在多種缺陷，包括點缺陷、缺陷團簇以及位錯、晶界等缺陷阱。缺陷的產生、聚集、釋放以及阱吸收等相互作用均會在缺陷演化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)為級聯碰撞過程產生尺寸為j的NDEF類型缺陷的速率；w(k,j)為缺陷尺寸由k變為j的轉化速率；Stot為材料內固有缺陷對尺寸為j的缺陷的吸收速率；主方程右側的第2項表示由尺寸k的缺陷團簇生成尺寸j的團簇的反應速率，右側第3項表示由尺寸j的團簇生成其他尺寸的團簇的生成速率，右側第4項表示材料內部固有缺陷對尺寸j的團簇的吸收項。

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

1.2 ETD算法

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

2 ETD方法應用于RT模型求解的典型算例

基于上述優勢，采用ETD方法對Ni中空位團簇演化以及奧氏體不銹鋼中微結構(位錯環以及空洞)演化過程建立RT模型并進行數值求解，并詳細介紹整個應用過程。

2.1 空位團簇演化模型

1) 模型建模與算法格式推導

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

(2)

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

(3)

(4)

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

(5)

初始條件通過下式計算：

(6)

為后續ETD格式的推導，上述模型可表示為如下矩陣形式：

(7)

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

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

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

(8)

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

ETD1格式：

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

(9)

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

(10)

F(t)由下式計算：

(11)

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

ETD4RK格式：

(12)

L由下式計算：

(13)

poly1、poly2、poly3由下式計算：

(14)

(15)

(16)

F(t)由下式計算：

(17)

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

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

(18)

(19)

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

(20)

φ1、φ2、φ3由下式計算：

(21)

(22)

(23)

其中，γ為調整系數，用于數值求解時的精度控制。

2) 空位團簇演化模型算例結果分析

基于文獻[9]中的物理模型與參數，本文嘗試采用上一節中給出的ETD算法對速率方程進行求解驗證。模型參數列于表1。

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

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

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

圖1 金屬鎳及鎳基合金823 K熱老化103 s時空位團簇尺寸分布采用不同算法的計算模擬結果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為缺陷數密度；n為缺陷尺寸；t0為初始時間；tm為第m個時間步，在演化過程中，空位數量值保持穩定，且與初始狀態空位數相同。

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

2.2 點缺陷可動模型

1) 模型建模與算法格式推導

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

圖2 守恒量驗證結果Fig.2 Verification result of conserved quantity

(25)

其中：

(26)

方程中所涉及的缺陷反應類型包括缺陷的產生、缺陷之間的復合、吸收、釋放等，點缺陷可動速率方程模型涉及的缺陷反應列于表2。

表2 點缺陷可動速率方程模型涉及的缺陷反應Table 2 Defect reaction involved in point defect moveable rate theory model

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

(27)

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

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

(28)

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

(29)

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

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

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

(30)

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

并行劃分后，模擬區域的邊界是原始方程對應的邊界，而子區域劃分產生的相鄰子區域之間的邊界為內部邊界。考慮邊界條件的計算，求解時需分為0號進程邊界和其他進程邊界。

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

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

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

2)點缺陷可動模型算例結果分析

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

采用預估校正ETD方法的計算結果、團簇動力學程序RADIEFF模擬結果以及實驗結果對比如圖5、6所示。

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

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

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

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

3 總結與展望

本文對RT方法的原理和建模以及ETD方法進行了介紹。通過對空位團簇演化模型和點缺陷可動模型兩類RT模型求解的計算實現，對ETD方法在求解模擬輻照損傷模型中的應用進行了詳細介紹，計算結果顯示了ETD方法在求解RT方程中的有效性。總之，使用ETD方法求解RT模型的優勢，主要體現在，ETD方法能對RT主方程進行高精度的數值求解，減少了額外簡化方程所帶來的誤差，而并行求解的引入使得對于大規模RT方程組的直接數值求解成為可能。可通過直接對主方程進行求解的方式，減少額外簡化方程所帶來的誤差，從而更好地對輻照損傷進行模擬。

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