分子動力學(xué)模擬在C/SiC材料虛擬試驗中的應(yīng)用與展望

李 堯，侯傳濤， 朱元夫，任 方，劉寶瑞，吳振強

(1.北京強度環(huán)境研究所 可靠性與環(huán)境技術(shù)重點實驗室, 北京 100076;2.北京強度環(huán)境研究所, 北京 100076)

0 引言

臨近空間、可重復(fù)使用等飛行器在高速飛行時面臨嚴酷的氣動力、熱和氧環(huán)境，其熱防護結(jié)構(gòu)需要同時滿足耐高溫、氣動維形、高承載等要求。近年來， 隨著具有耐高溫、高比強度和高韌性特點的陶瓷基復(fù)合材料制備工藝逐漸成熟，高速飛行器熱防護結(jié)構(gòu)研制在“冷熱相分離”結(jié)構(gòu)的基礎(chǔ)上，一些關(guān)鍵部位逐步采用“防熱承載一體化”結(jié)構(gòu)。這些新材料、新型結(jié)構(gòu)的應(yīng)用使飛行器性能指標得到了提升，但同時給地面驗證和評估試驗帶來困難。一方面，更加嚴酷高溫以及復(fù)雜的多物理場耦合環(huán)境，使試驗環(huán)境模擬難度大、成本高；另一方面，溫度、熱流、應(yīng)變以及損傷失效參數(shù)測試變得更為困難。

近年來，材料結(jié)構(gòu)性能虛擬試驗技術(shù)逐漸受到關(guān)注并得到應(yīng)用，不僅是對真實試驗的補充，而且逐步成為評估和驗證材料結(jié)構(gòu)性能的新的技術(shù)手段和發(fā)展方向。以C/SiC為代表的新一代陶瓷基復(fù)合材料優(yōu)異的綜合力學(xué)性能來自于多尺度結(jié)構(gòu)設(shè)計，其內(nèi)部包含亞微米厚度的界面層、數(shù)十微米直徑的碳纖維、毫米尺寸的纖維束以及SiC基體。此外，制備工藝中不可避免產(chǎn)生的大量裂紋、孔隙等缺陷，使其宏觀力學(xué)行為與傳統(tǒng)金屬材料顯著不同。國內(nèi)外學(xué)者針對C/SiC材料開展了大量多尺度建模與分析工作，能夠較準確預(yù)示室溫下的力學(xué)行為。然而，這些多尺度分析多基于有限元方法，其預(yù)示精度依賴于強度準則和材料參數(shù)的選取。對高溫、氧環(huán)境下C/SiC材料力學(xué)行為機理和規(guī)律的認識比較初步，材料參數(shù)難以獲取、失效準則缺失。因此，C/SiC材料和結(jié)構(gòu)的虛擬試驗仍然面臨挑戰(zhàn)。

分子動力學(xué)方法是一種通過模擬構(gòu)成材料的原子運動，自下而上地預(yù)測材料性能及力學(xué)行為的方法；但是由于考慮所有組成原子自由度，計算量十分龐大。隨著高性能計算、人工智能等技術(shù)發(fā)展，分子動力學(xué)計算能力得到顯著提升，應(yīng)用場景不斷豐富，對解決C/SiC虛擬試驗的瓶頸問題，展現(xiàn)出較大潛力。本文將首先介紹分子動力學(xué)的基本原理，然后對在C/SiC材料模擬方面的進展進行簡要綜述，最后對未來的應(yīng)用場景和發(fā)展方向進行展望。

1 分子動力學(xué)方法介紹

分子動力學(xué)有著堅實的理論基礎(chǔ)并在金屬材料體系中得到大量應(yīng)用。近年來，結(jié)合機器學(xué)習(xí)等技術(shù)，分子動力學(xué)模型和算法也在不斷發(fā)展。本章將從基本原理和發(fā)展趨勢兩個方面對分子動力學(xué)方法進行介紹。

1.1 基本原理

微觀尺度上，通常認為由原子或分子運動起主導(dǎo)作用，常用于分析固體形成的微結(jié)構(gòu)和缺陷，如：位錯、晶界、多相等，發(fā)展出晶格動力學(xué)、位錯動力學(xué)等理論分析方法，以及分子動力學(xué)、蒙特卡洛等模擬方法。分子動力學(xué)不僅模擬體系的能量，還可以獲得每一時刻原子的受力，在力學(xué)行為研究方面相較蒙特卡洛模擬具有一定優(yōu)勢。分子動力學(xué)以統(tǒng)計力學(xué)為理論基礎(chǔ)，忽略原子核和電子的量子效應(yīng)，將原子核作為經(jīng)典粒子，用經(jīng)典的牛頓運動方程描述組成系統(tǒng)的原子運動過程。式(1)中的，和分別代表時刻第個原子的位置、質(zhì)量以及受力。體系的總勢能()可以通過力場或勢函數(shù)計算獲得，利用式(2)可以計算新的位置處的原子受力。

(1)

(2)

分子動力學(xué)可以計算平衡態(tài)系統(tǒng)的結(jié)構(gòu)和能量，以及非平衡態(tài)系統(tǒng)動態(tài)演化過程中原子的力與運動，例如變形、相變、化學(xué)反應(yīng)等。對原子與原子間相互作用的描述建立在力場(Force field)或勢函數(shù)之上：在給定的初始條件下，首先利用力場或勢函數(shù)計算系統(tǒng)的勢能()；然后求解每個原子受力并計算原子的運動過程，獲得新的原子位置和原子速度矢量；重復(fù)這一過程，最終得到體系內(nèi)原子位置隨時間的演變過程。相互作用原子的種類、所處的化學(xué)環(huán)境等不同，力場或勢函數(shù)對各種原子定義也不同，并采用不同的力常數(shù)和經(jīng)驗勢參數(shù)來描述成鍵、非鍵等相互作用。針對不同體系已經(jīng)形成了Len-nard-Jones勢、Tersorff勢、SW勢、Morse勢、嵌入原子勢(Embedded Atom Method)、CHARMM勢、AMBER勢等經(jīng)驗勢函數(shù)。經(jīng)驗勢參數(shù)和力常數(shù)可以依據(jù)材料物性的實驗結(jié)果獲得，或者通過擬合量子力學(xué)計算結(jié)果獲得。

目前，根據(jù)研究領(lǐng)域和對象的不同，分子動力學(xué)模擬已開發(fā)了LAMMPS、GROMACS、NAMD等計算軟件，以及OVTIO、VMD等可視化軟件，基于高性能計算可實現(xiàn)上億原子體系的模擬。

1.2 發(fā)展趨勢

現(xiàn)實中很多現(xiàn)象跨越多個時間和空間尺度，例如，固體變形直至破壞的過程，跨越了從原子結(jié)構(gòu)到宏觀結(jié)構(gòu)的多個尺度，是一種復(fù)雜的多尺度現(xiàn)象。為了更大的模擬，基于分子動力學(xué)的多尺度模擬的研究也在興起，主要以粗粒化(Coarse Grained)分子動力學(xué)方法、橋域法(Bridging domain Method)等多尺度方法為代表。粗?；肿觿恿W(xué)方法是通過聯(lián)合原子來構(gòu)建模型，例如將質(zhì)量較小的原子附著在重原子上，形成具有新的范德華半徑和質(zhì)量的聯(lián)合原子，而多個聯(lián)合原子可以更進一步簡化為一個粒子以降低體系的自由度。粗?；Ｐ秃喕瞬糠纸Y(jié)構(gòu)細節(jié)，使總粒子數(shù)目得到降低，因此所能模擬的空間尺度可增加至原來的幾倍，時間尺度提高到原來的幾十倍。橋域方法代表多尺度方法中的一類重要方法，將研究區(qū)域分為不同區(qū)域，采用不同尺度方法進行計算，在區(qū)域的邊界處進行信息交換。針對固體破壞問題建立的MAAD(Macroscopic, Atomis-tic, Abinitio Dynamics)方法認為對于固體的破壞總是從局部發(fā)生，所以可以在局部建立小尺度模型，而在其他地方建立大尺度模型：例如在裂紋極端采用緊束縛模型模擬原子鍵的破壞，在周圍采用利用分子動力學(xué)模擬處理裂紋面,而在更遠的地方使用連續(xù)介質(zhì)模擬來模擬斷裂過程。耦合QM/MM(Quantum Mechnics/Molecular Mechanics)方法是一種結(jié)合量子力學(xué)的計算精度和分子動力學(xué)計算速度的原子模擬方法，將系統(tǒng)分為電子性質(zhì)重要的區(qū)域用量子力學(xué)計算而其他區(qū)域采用經(jīng)典描述。類似地還發(fā)展出了分子動力學(xué)與有限元耦合方法等。經(jīng)過多年的發(fā)展加深了金屬材料斷裂破壞的認識，并推動了制備工藝的進步和新合金材料的設(shè)計，例如高強高韌的納米孿晶、高溫超合金、高熵合金等。

2 C/SiC材料的模擬

陶瓷基復(fù)合材料是由纖維、基體、界面層等組分復(fù)合形成的多相非均質(zhì)材料，結(jié)構(gòu)本身具有多尺度特征，還包含大量的大大小小的孔隙和缺陷。目前，陶瓷基復(fù)合材料模擬研究仍以有限元方法為主，針對微觀、細觀、宏觀結(jié)構(gòu)分別建模計算，而不同尺度間信息可采用傳遞、協(xié)同、并發(fā)等耦合方式。其中，以代表性體積單元方法(Representative Volume Element)和通用單胞模型方法(Generalized Method of Cells)應(yīng)用最為廣泛。代表性體積單元方法的主要思路是選取合適區(qū)域建立細觀結(jié)構(gòu)模型作為代表性體積單元，通過有限元計算獲取應(yīng)力、應(yīng)變場量，最后通過均勻化方法對細觀場量進行處理，得到宏觀材料的應(yīng)力-應(yīng)變關(guān)系。美國國家航空航天局(NASA)針對復(fù)合材料結(jié)構(gòu)多尺度仿真分析提出的通用單胞模型方法，具有宏、細觀結(jié)構(gòu)統(tǒng)一的特點，在細觀尺度上，進行簡化近似并通過解析求解獲得細觀應(yīng)力和應(yīng)變場，而宏觀尺度上利用細觀均勻化結(jié)果進行數(shù)值求解，可用于復(fù)合材料結(jié)構(gòu)在復(fù)雜載荷下的計算分析。然而，有限元計算依賴于微/細觀材料組分物理性能的準確獲取，同時難以直接模擬氧化、燒蝕等化學(xué)反應(yīng)過程。

近年來，分子動力學(xué)模擬研究陶瓷基復(fù)合材料物性、化學(xué)反應(yīng)過程的工作不斷涌現(xiàn)，對獲取陶瓷基復(fù)合材料關(guān)鍵力學(xué)性能提供了新的途徑。下面將通過幾個例子介紹分子動力學(xué)在C/SiC材料模擬方面的近期研究進展。

2.1 界面相微結(jié)構(gòu)建模

陶瓷基復(fù)合材料中亞微米厚度的界面層對其宏觀性能影響巨大。熱解碳(Pyrolytic Carbon，PyC)是C/C、C/SiC、SiC/SiC等最為常用的界面相材料之一，但由于難以獲得較大塊體狀的熱解碳材料，目前對于其力學(xué)性能的實驗研究較少，且結(jié)果分散性較大。為此，許多研究者從原子結(jié)構(gòu)角度出發(fā)，對熱解碳的微觀模型進行了深入研究。2003年，西北工業(yè)大學(xué)Ye等基于對C/C復(fù)合材料中的熱解碳結(jié)構(gòu)、模量、強度等性能開展了分子動力學(xué)模擬研究。2012年，Leyssale等采用反向蒙特卡洛(RMC)方法通過對HRTEM圖像數(shù)據(jù)處理、三維重構(gòu)，建立了熱解碳的原子構(gòu)型，如圖1 (a)和(b)所示；并且模擬得到納米尺度熱解碳在的彈性常數(shù)，雖然與微米級界面相的力學(xué)性能相比有一定差別，但是對于認知熱解碳微觀結(jié)構(gòu)與力學(xué)性能之間的關(guān)系具有重要意義。2019年，中國科技大學(xué)Chen等利用反應(yīng)力場模擬了熱解碳在碳纖維表面的初始沉積生長過程，如圖1(c)所示。2020年，Chen等根據(jù)沉積的結(jié)構(gòu)設(shè)計了代表性的孔洞結(jié)構(gòu)，圖1(d)所示，并提出一種便捷的熱解碳建模方法：首先在多層石墨烯片層中隨機填充該孔洞結(jié)構(gòu)，然后弛豫得到熱解碳原子構(gòu)型；利用獲得的模型，研究了熱解碳在拉、壓和剪切等載荷作用下的微觀力學(xué)行為。這些工作從原子尺度加深了研究人員對熱解碳結(jié)構(gòu)和物性的認識，對于進一步理解陶瓷基復(fù)合材料的界面性能以及斷裂、脫粘等失效過程具有重要啟發(fā)。

(a) HRTEM實驗數(shù)據(jù)三維重構(gòu)圖像

(b) 基于HRTEM圖像構(gòu)建的熱解碳原子構(gòu)型

(c) 熱解碳沉積過程模擬

(d) 熱解碳中的代表性孔洞結(jié)構(gòu)圖1 熱解碳原子建模[8-10]Fig.1 Atomic structures of pyrolytic carbon[8-10]

2.2 氧化反應(yīng)過程模擬

C/SiC復(fù)合材料結(jié)構(gòu)在高溫環(huán)境下長時間使用過程中，氧化是結(jié)構(gòu)損傷失效的重要因素。C/SiC的氧化過程極為復(fù)雜，伴隨著纖維氧化燒蝕、基體裂紋閉合等多種過程，并且這些過程受溫度控制和影響。目前，基于反應(yīng)和擴展理論，對于C/SiC的氧化反應(yīng)機理有了定性的認識，但對氧化反應(yīng)過程定量分析研究仍然處于初步階段，氧氣擴散速率、氧化反應(yīng)速率等關(guān)鍵參數(shù)仍然缺失。

(a) 在1 100 ℃下SiC(0001)表面的氧化過程

(b) 處于氧氣環(huán)境中的碳纖維原子結(jié)構(gòu)模型

(c) 碳纖維氧化過程中CO分子、O原子和CO2分子數(shù)目隨時間的變化圖2 氧化過程模擬 [14,16]Fig.2 Simulation of oxidation processes [14,16]

2.3 力學(xué)性能仿真預(yù)示

近年來，學(xué)者們針對陶瓷、界面相等組分力學(xué)性能以及復(fù)合材料微觀力學(xué)行為開展了相關(guān)研究。2019年，Sun等通過分子動力學(xué)模擬研究了高溫濕氧環(huán)境對于Si薄膜中裂紋擴展的影響，如圖3(a)所示。2015年，李麗麗等基于Tersoff作用勢，通過分子動力學(xué)計算研究了無定型碳界面層厚度對SiC納米纖維/SiC復(fù)合材料的斷裂力學(xué)行為的影響，發(fā)現(xiàn)界面層厚度增加會降低纖維的應(yīng)力集中系數(shù),提高斷裂能，使裂紋穿透纖維的脆斷模式轉(zhuǎn)變?yōu)槔w維拔出失效模式，起到增強補韌的效果，如圖3(b)所示。圖3(c)展示了 Zhang 等理論和實驗結(jié)合研究了熱解碳微柱在拉伸和壓縮下的強度性能。2019年，Zhou等通過分子動力學(xué)模擬研究了CVI工藝C/C復(fù)合材料中熱解碳/碳纖維在拉伸、剪切等載荷下的界面力學(xué)性質(zhì)，將獲得的界面剪切模量、強度等參數(shù)帶入有限元模型的內(nèi)聚力單元中，進行了纖維拔出的有限元仿真分析，如圖3(d)所示。

(a) 高溫下Si薄膜在真空、氧氣和水氧環(huán)境中的裂紋擴展

(b) 預(yù)置基體裂紋的SiC纖維/無定型C過渡層/SiC復(fù)合材料模型

(c) 熱解碳微柱單軸拉伸實驗和模擬

(d) 熱解碳/碳纖維界面的一種典型失效過程圖3 力學(xué)行為及性能研究[17-20]Fig.3 Molecular dynamics studies on mechanical behaviors and properties [17-20]

3 總論與展望

分子動力學(xué)方法在C/SiC材料微結(jié)構(gòu)預(yù)測、氧化反應(yīng)模擬，損傷失效機制、化學(xué)反應(yīng)等方面已經(jīng)開展了一些探索性研究。本文認為分子動力學(xué)在以下幾個方面仍具有較大的應(yīng)用潛力：

(1)支撐材料C/SiC工藝改進以及新材料的研制

分子動力學(xué)模擬結(jié)合電鏡、Raman、CT等測試手段，可以加深對復(fù)合材料關(guān)鍵組分，例如界面層的微結(jié)構(gòu)形成及演化規(guī)律的認識，幫助建立其微觀和宏觀性能間關(guān)系，實現(xiàn)對材料設(shè)計和性能調(diào)控。此外，基于分子動力學(xué)的反應(yīng)過程模擬可以幫助優(yōu)化材料制備工藝參數(shù)。對于新材料的研制，可以通過分子動力學(xué)進行材料性能預(yù)測、先進行初步的篩選，排除錯誤選項，減少實驗工作量，加速研制過程。

(2)支撐極端復(fù)雜環(huán)境下結(jié)構(gòu)失效分析

對于難以開展試驗或測試的極端復(fù)雜環(huán)境，可以通過分子動力學(xué)模擬獲得給定載荷環(huán)境下材料的應(yīng)力應(yīng)變行為和破壞強度，擬合應(yīng)力應(yīng)變曲線便可以獲得本構(gòu)模型參數(shù)。對于氧化過程，可以類似地擬合出氧化反應(yīng)速率、氧化損傷因子等參數(shù)。這些模型和參數(shù)可以直接用于結(jié)構(gòu)的變形及失效分析，幫助厘清失效機理及傳播過程。

(3)支撐結(jié)構(gòu)性能多尺度虛擬試驗技術(shù)發(fā)展

以C/SiC為代表的復(fù)合材料結(jié)構(gòu)性能不僅取決于其結(jié)構(gòu)參數(shù)，還依賴于材料的微細觀特征參數(shù)以及在結(jié)構(gòu)上的取向和分布。針對這類材料結(jié)構(gòu)，需要考慮其多尺度的特征，在虛擬試驗中建立多尺度分析模型。分子動力學(xué)除了可以提供材料和模型參數(shù)外，還可以對基于有限元方法多尺度模型進行補充和發(fā)展，例如分子動力學(xué)和有限元耦合模型的構(gòu)建。

綜上所述，隨著分子動力學(xué)計算能力提升以及新模型新算法的不斷提出，必將對材料結(jié)構(gòu)虛擬試驗技術(shù)的發(fā)展帶來新的機遇。