分子動力學模擬在C/SiC材料虛擬試驗中的應用與展望

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

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

0 引言

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

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

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

1 分子動力學方法介紹

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

1.1 基本原理

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

(1)

(2)

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

目前，根據研究領域和對象的不同，分子動力學模擬已開發了LAMMPS、GROMACS、NAMD等計算軟件，以及OVTIO、VMD等可視化軟件，基于高性能計算可實現上億原子體系的模擬。

1.2 發展趨勢

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

2 C/SiC材料的模擬

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

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

2.1 界面相微結構建模

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

(a) HRTEM實驗數據三維重構圖像

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

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

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

2.2 氧化反應過程模擬

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

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

(b) 處于氧氣環境中的碳纖維原子結構模型

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

2.3 力學性能仿真預示

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

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

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

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

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

3 總論與展望

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

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

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

(2)支撐極端復雜環境下結構失效分析

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

(3)支撐結構性能多尺度虛擬試驗技術發展

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

綜上所述，隨著分子動力學計算能力提升以及新模型新算法的不斷提出，必將對材料結構虛擬試驗技術的發展帶來新的機遇。