傳統的試驗方法研究材料的宏觀形變行為已非常成熟,但采用該方法存在試驗周期長、成本高等不足,并且無法直接觀測材料內部的微觀形變及力學行為。因此,為了降低材料特性研究時對試驗的過分依賴,可描述多晶相互作用的晶體塑性理論及分子動力學理論得以發展。并且,近年來計算機性能的提升也加速了基于晶體塑性理論或分子動力學理論的金屬材料形變特性的研究進程。但是,晶體塑性有限元理論多用于模擬晶粒級細觀尺度下的形變行為,無法體現材料在形變過程中的原子構型及位錯的演化,而用于描述原子在外力作用下最基本運動細節的分子動力學模擬則可以彌補這一短板。分子動力學模擬需依據材料經驗力場及原子初始位置計算其勢能,通過求導的方法分析體系內各個原子的受力情況,進而基于牛頓運動方程,在極短的時間間隔(1×10-15s)內計算原子新的位置及受力,如此反復直至計算收斂[1-2]。因此,分子動力學模擬可以實時捕捉體系的動態形變,為在原子尺度上研究微觀形變機制奠定了基礎。ChenB等[3]建立了具有多種取向位錯網的單晶鎳基高溫合金的分子動力學模型,通過對比體系能量及位錯網演化過程,為建立更加精確的界面位錯網模型提供依據。之后,LiNL和WuWP等[4]對引入界面位錯網的分子動力學模型施加單軸拉伸載荷及循環載荷,通過分析材料形變過程中的位錯演化、層錯/滑移運動及應力-應變行為,評估應變率、溫度對材料力學行為及循環形變機制的影響。LiXT等[5]揭示了位錯發射與孔洞和缺陷的相互作用機理,認為位錯發射會導致裂紋尖端鈍化,是裂紋尖端脆性擴展演化為韌性擴展的主要因素。祝占偉[6]采用分子動力學模擬分析了鎳基合金中五重孿晶的形成機制及形變機理。本文建立了LAMMPS分子動力學仿真原子結構模型,進行了室溫(300K)與高溫(600和920K)條件下的單軸拉伸仿真,討論了溫度對模型內部位錯演化及力學行為的影響。

1、分子動力學模型

1.1　相互作用勢函數

本文采用嵌入原子法(Embedded-AtomMethod,EAM)勢函數[7]來描述原子間作用力的相互作用,EAM勢是一種描述金屬鍵最常用、最精確的相互作用模型,并且常用于描述原子尺度下的擴散、變形、熔化和燒結等過程[8-10]。EAM勢中總能量Etot的計算公式如下,其中,系數1/2是為了避免一對相互作用的重復計算。

式中:vij(rij)為分離距離為rij的原子i與原子j之間的一對相互作用能;Fi為原子i的嵌入能,可以描述為局部電子密度ρi的函數;N為體系內共參與計算的原子對的數量。ρi的表達式為:

式中:ρj(rij)為將原子i賦值到分離距離為rij的原子j的電子密度函數。采用用于表示接近等原子成分組成下的Fe-NiCr-Co-Al的高度理想化面心立方FCC混合物的EAM勢函數[7],該EAM勢旨在模擬復雜系統的性質,可以較為精確地模擬系統中各類原子之間的相互作用。

1.2　原子結構模型

GH4169鎳基高溫合金粉末中存在超過10種元素,而Ni、Cr、Al與Fe元素的占比超過了元素總體含量的90%(質量百分數)。因此,為了在簡化原子結構模型建模的同時,盡可能地保證多體相互作用的準確性,建立GH4169鎳基高溫合金原子結構模型時僅考慮Ni、Cr、Al與Fe共4種主要元素[9],4種元素的含量分別為52%、10%、1%與37%,即Ni52Fe37Cr10Al1。由于GH4169鎳基高溫合金中Ni元素含量最多,因此,在后續建模過程中,首先,依據鎳基合金晶體類型建立適當尺寸的面心立方結構Ni單晶晶格模型。然后,依據Cr、Al與Fe這3種元素占比,分別在Ni單晶晶格模型中的FCC格位上隨機替換Ni原子,形成一種高度類似于GH4169鎳基高溫合金元素組成的Ni52Fe37Cr10Al1原子結構,用于后續分子動力學計算。需要說明的是,Ni元素的點陣常數aNi=3.52?作為常數用于后續建模。計算原子尺度下的GH4169鎳基高溫合金單軸拉伸應力-應變關系時,首先,建立尺寸(無量綱)為40×40×40的面心立方結構Ni單晶晶格模型,結合Ni元素的點陣常數可知,該晶格模型在x、y與z方向上的實際尺寸分別為40aNi、40aNi與40aNi。最終模型包含的原子總數為256000個,其中Ni原子有133120個、Fe原子有94720個、Cr原子有25600個、Al原子有2560個,最終模型如圖1所示。

施加單軸拉伸載荷之前,首先在NVE(原子數-體積能量,Number-Volume-Energy)系綜(恒定原子序數,恒溫,不恒壓)下對模型弛豫50ps,保證系統處于平衡狀態。然后,利用Nose-Hoover熱浴法設置環境溫度,并且通過NVT(原子數-體積-溫度,Number-Volume-Temperature)系綜用應變控制對原子結構模型施加單軸拉伸載荷,應變率為1×109s^-1。計算過程中每間隔5ps輸出環境熱力學信息,以監控仿真是否平穩運行。通過設置變量的方法計算不同時刻體系的應力與應變數值,以繪制應力-應變曲線。仿真過程時間步長為0.001ps,每隔5ps輸出體系形變結果。此外,通過改變環境溫度(300、600和920K),以探究溫度對GH4169鎳基高溫合金應力應變關系的影響。需要說明的是,整體計算過程包含minimize過程、弛豫過程、fix-deform過程,拉伸載荷加載過程中以載荷施加完畢為形變過程結束時刻,因此,并未關注整體計算時間。5000個時間步長是為了降低計算結果文件體積的同時保證輸出計算結果可以精確反應整體形變過程。

2、結果與討論

2.1　不同溫度下的微觀應力-應變關系

為了探究原子級納米尺度下溫度對GH4169鎳基高溫合金單軸拉伸力學性能及形變機制的影響,在300、600和920K條件下對原子結構施加單軸拉伸載荷,計算得到的真實應力-真實應變曲線,如圖2所示。研究表明,分子動力學模擬過程中應變率的大小對材料彈性行為的影響可以忽略,更高的應變率僅表現出更高的屈服強度與更差的塑性[4]。從圖2中可以看,不同溫度下的真實應力-真實應變曲線均可劃分為3個區域,即隨著應變的增加,依次為應力線性增加的彈性區、應力驟降的塑性激活區與應力穩定的應力釋放區[11]。與宏觀應力-應變曲線不同的是,微觀尺度下的應力-應變曲線在進入塑性激活區后應力變化呈現先驟降后波動的特征,這是因為在進入屈服階段后原子結構內部會形成大量位錯,隨著應變的增加,位錯會增殖與運動,因此,位錯與位錯之間、位錯與其他晶體學結構之間會出現相互作用,這些相互作用一定程度上會使位錯塞積,甚至消失,最終導致應力波動。應力驟降是由于Shockley位錯發生反應,演化為其他類型位錯或泯滅,此時體系內應力釋放,導致應力下降。此外,隨著溫度的升高,GH4169鎳基高溫合金的屈服強度ReL逐漸下降,分別為22.38、18.58和15.02GPa。同樣地,屈服應變也隨溫度的升高而降低,分別為16.65%、14.63%和12.67%,這表明溫度越高,納米原子GH4169鎳基高溫合金結構越容易形成位錯,導致塑性形變。利用不同溫度下彈性區的真實應力-真實應變曲線擬合得到的彈性模量分別為125.9、125.3和119.1GPa,可以看出,當溫度低于920K時,材料的彈性模量十分接近,而當溫度達到920K時,彈性模量下降了6.8GPa,這從納米尺度上證明了GH4169鎳基高溫合金具有良好的920K高溫抗性。

2.2　形變過程中原子構型演化

圖3為不同溫度下隨應變增加原子結構模型中3種不同晶體結構類型的原子占比變化情況。從圖3中可以看出,形變初期3種晶體結構所含原子數保持穩定,且均為面心立方FCC結構,這與彈性階段對應。隨著應變的增加,面心立方FCC結構原子占比降低,與此同時,體心立方BCC與六方密排HCP結構原子占比增多,此時對應應力-應變曲線的塑性激活區。整體上,300、600和900K的溫度下應變分別達到17%、13%和8%時,體系內部FCC結構原子占比開始降低,300K下這種變化是突然發生的,但溫度越高,變化發生的越平緩。此外,無論此時體系溫度高低,FCC結構原子占比最多,BCC與HCP結構原子占比較少。當應變繼續增加,進入應力釋放階段后,在此階段3種晶體結構的占比隨著應變增加而趨于穩定。由此可見,原子級微觀尺度下原子排布的變化直接影響著材料的應力-應變關系,這是因為FCC結構原子占比減小會導致體系內部出現HCP與BCC結構,因此,引起位錯、層錯、滑移等結構出現。此外,溫度越高,原子運動更活躍,從而降低了原子間結合力,更容易誘發體系內部原子排序紊亂,更早地誘發不同類型位錯的形核,最終體現在彈性模量與應力-應變曲線形貌的差別上。經過以上分析可知,FCC結構原子的演化會誘發晶體學缺陷,因此,在分析材料單軸拉伸形變機制時應將完美的FCC結構原子去除。借助位錯提取算法(DislocationExtractionAlgorithm,DXA)對體系內的位錯類型及分布進行分析,將分析結果與去除FCC結構原子的體系相結合,以分析GH4169鎳基高溫合金的形變機制。圖4為3種溫度下不同應變時的原子排列結構及位錯結構特征,其中,片層狀原子為HCP結構,零星分布狀的單個原子為BCC結構,彌散團簇狀與無序狀原子為無序結構。同樣,以深淺不一的實線表征不同類型的位錯。

從圖4中可以看出,位錯的形核與運動使得體系內部原子結構重排,引發明顯的HCP結構層錯。層錯與層錯之間往往以一系列無序排列的原子為界。對于具有FCC結構的GH4169鎳基高溫合金而言,內部存在的堆垛層錯可以等同于HCP結構,宏觀下結構的演變表現為塑性形變。隨著應變的增加,層錯由塑性階段時體系兩端的多形核逐步向體系中間運動,這證明此時會發生頸縮現象。位錯大多出現于層錯周圍,這是因為在外力的作用下,體系內部異常排列的原子極易誘發位錯形核,從而降低材料的力學性能。隨著應變增加,位錯的運動引導層錯與層錯之間的相互融合趨于穩定,形成沿某個45°方向貫穿整個體系的層錯,即形成滑移。此外,位錯運動會使異常排序的原子逐步演變成FCC結構,并且吸引異常排序的原子團簇聚集于融合后的層錯周圍,宏觀表現為塑性形變。需要注意的是,在形變過程中體系內會出現由無序原子包裹的空位缺陷,這些空位可以被位錯吸收,從而牽引位錯與層錯運動,這與體系內部不同結構原子占比的變化相一致。另一方面,高溫下原子的劇烈運動會導致異常排序的原子增多,為位錯形核提供了有利條件的同時增加了體系內部空位的數量。這是因為對于晶體材料而言,空位的形成會受到溫度的影響,溫度越高,空位形成的概率越大,在此情況下,位錯及層錯更容易發生運動,與之對應的是材料力學性能的下降。

2.3　形變過程中位錯演化

圖5僅體現了體系內部位錯的類型及特征,以清楚地分析形變過程中位錯的特征及演化。圖5a、圖5d與圖5g分別為300、600與920K溫度下原子結構體系進入塑性階段時的位錯形貌、類型及分布。從圖5中可以看出,材料進入塑性階段時,位錯主要以1/6<112>Shockley不全位錯為主,與此同時存在1/2<110>Perfect全位錯、1/6<110>Stair-rod壓桿位錯、1/3<111>Frank分位錯、1/3<100>Hirth分位錯及其他位錯,不同類型的位錯具有不同的柏氏矢量,并且1/3<111>Frank分位錯及1/3<100>Hirth分位錯的含量較小。從圖5a、圖5d與圖5g可以看出,材料剛進入塑性階段時內部位錯長度小、數量多,并且呈現點發射特征,這是因為1/6<112>Shockley不全位錯是從空位表面發射,并且隨著應變的增加,位錯相互結合從而形成位錯環,此時部分1/6<112>Shockley不全位錯在滑移面上以銳角相交時會形成1/6<110>Stair-rod壓桿位錯,而以鈍角相交時則會形成1/3<100>Hirth分位錯。

Shockley不全位錯通常以刃形位錯形式出現,也可以是螺位錯或混合型位錯,在外力作用下可以在{111}滑移面上運動,從而牽引層錯區擴大或縮小,當兩個層錯區之間夾帶了一個定位錯時,會形成一個Stair-rod壓桿位錯,如圖5b、圖5e與圖5h所示。此外,還可以觀察到3種溫度下均會形成由多個Stair-rod壓桿位錯包圍形成的位錯環,結合圖4可知,所圍區域為層錯區域。由于Stair-rod壓桿位錯屬于不可動位錯,因此,一定程度上可以在塑性變形階段使材料發生硬化。應變逐漸增大的同時位錯也在不斷演化,在此過程中新的Stair-rod壓桿位錯會形成,并且部分Stair-rod壓桿位錯會發生解體與擴散,演化為兩個可以進行塑性形變的跟蹤位錯,即Shockley不全位錯,如圖5c、圖5f與圖5i所示。在整個拉伸形變過程中,體系內部發生位錯形核的同時也會出現空位,空位會引導位錯運動,位錯運動又會形成新的空位,如此反復,材料內部的層錯發生移動和擴展,最終導致滑移。對比不同溫度下的位錯演化過程可知,溫度越高,位錯形核位置越均勻,這是因為高溫環境下體系內部異常排布的原子的數量多,增大了位錯形核的概率。此外,溫度越高,位錯排布方向的一致性更高,且位錯長度更長,這是位錯之間更加頻繁的相互作用導致的。對不同溫度和應變時的位錯密度進行評估,結果如圖6所示。位錯密度[12]計算公式如下:

式中:ρ為位錯密度;ldis為體系中位錯總長度;V為原子結構體積。從圖6可以看出,3種溫度下的總位錯密度與應力-應變曲線具有一致的變化規律。彈性范圍內,隨著應變的增加,體系內部各類位錯的數量均為0,當材料進入塑性階段時,Shockley不全位錯急劇增加,300、600與920K條件下Shockley不全位錯密度的峰值分別為0.0057、0.0048與0.0033nm-2,位錯密度的急劇升高證明了材料內部出現了大量層錯。隨后進入應力釋放階段,在此過程中Shockley位錯密度持續降低,而其他類型的位錯密度則有所升高,這是因為形變過程中Shockley位錯之間不同的相互作用形式使其演化為其他類型位錯。此時Stair-rod壓桿位錯密度在相對穩定的范圍內波動,分別為0.0006、0.0005與0.0004nm-2,Stair-rod壓桿位錯密度的相對穩定證明了材料具有一定的抗變形能力。而造成這種波動的原因是結構為了可以在外力的作用下持續發生變形,迫使Stair-rod壓桿位錯分解,與此同時又有新的Stair-rod壓桿位錯形成。總體而言,相對較高的位錯密度水平表明,形變過程中材料內部位錯反應頻繁,并且滑移也較為活躍,尤其是在高溫環境下,位錯的相互反應更容易發生。

3、結論

(1)室溫(300K)、高溫(600和920K)的環境下,GH4169鎳基高溫合金的微觀彈性模量分別為125.9、125.3和119.1GPa,且彈性變形范圍內的應力-應變曲線形貌差別不大,證明了GH4169鎳基高溫合金具有良好的高溫(920K)抗性。

(2)300、600和900K的溫度下應變分別達到17%、13%及8%時體系內部FCC結構原子占比開始降低,溫度越高,FCC結構數量降低至最小值時的速率越低。依據體系內部原子結構的演化可將整體形變過程分為彈性區、塑性激活區及應力釋放區。

(3)進入塑性激活區后,大量Shockley不全位錯形成,并且在隨后的形變過程中,Shockley不全位錯密度逐漸降低,演化為其他類型位錯。進入應力釋放區后,各類位錯密度保持恒定,滑移的增多及運動使材料出現了頸縮現象。

參考文獻:

[1]Plimpton S. Fast parallel algorithms for short-range molecular dynamics [J]. Journal of Computational Physics, 1995, 117: 1-19.

[2]Thompson A P, Aktulga H M, Berger R, et al. LAMMP—A flexiblesimulation tool for particle-based materials modeling at the atomic, meso, and continuum scales [ J ]. Computer Physics Communications, 2022, 271: 10817.

[3]Chen B, Wu W P, Chen M X, et al. Molecular dynamics study of fatigue mechanical properties and microstructural evolution of Nibasedsingle crystal superalloys under cyclic loading [J]. ComputationalMaterials Science, 2020, 185: 109954. [4] 　 Li N L, Wu W P, Nie K. Molecular dynamics study on the evolutionof interfacial dislocation network and mechanical properties of Ni-based single crystal superalloys [J]. Physics Letters A, 2018, 382: 1361-1367.

[5]Li X T, Peng S Y, Zhang X, et al. Microscopic and macroscopic analyses of the interaction mechanism between defect growth and dislocationemission in single-crystal aluminum [J]. Fatigue & Fractureof Engineering Materials & Structures, 2021, 44: 3008-3022.

[6]祝占偉. 鎳基合金中五重孿晶形成及變形機理的分子動力學研究 [D]. 大連: 大連理工大學, 2018. Zhun Z W. Study on the Formation and Deformation Mechanisms of Fivefold Twins in Nickel Base Alloy [D]. Dalian: Dalian Universityof Technology, 2018.

[7]Frakas D, Caro A. Model interatomic potentials for Fe-Ni-Cr-CoAlhigh-entropy alloys [J]. Journal of Materials Research, 2020, 35: 3031-3040.

[8]Nandy J, Yedla N, Gupta P, et al. Sintering of AlSi10Mg particlesin direct metal laser sintering process: A molecular dynamics simulation study [ J]. Materials Chemistry and Physics, 2019, 236: 121803.

[9]Wang J Q, Shin S, Hu A M, et al. Diffusion kinetics of transient liquid phase bonding of Ni-based superalloy with Ni nanoparticles: A molecular dynamics perspective [ J]. Computational Materials Science, 2018, 152: 228-235.

[10]Wu W P, Yao Z Z. Molecular dynamics simulation of stress distributionand microstructure evolution ahead of a growing crack in singlecrystal nickel [J]. Theoretical and Applied Fracture Mechanics,2012, 62: 67-75.

[11]楊攀峰. Ni-Co 合金拉伸力學性能和變形機制的模擬研究[D]. 蘭州: 蘭州理工大學, 2020.

Yang P F. The Simulation Research of Tensile Mechanical Propertiesand Deformation Mechanism of Ni-Co Alloys [D]. Lanzhou: Lanzhou University of Technology, 2020.

[12] Li L, Shen L M, Proust G. Fatigue crack initiation life prediction for aluminium alloy 7075 using crystal plasticity finite element simulations [J]. Mechanics of Materials, 2015, 81: 84-93.

相關鏈接