非均質材料與位錯交互能的數值等效夾雜算法

李 璞，朱 凱，侯佳卉，謝東東，錢厚鵬，金曉清,

(1. 重慶大學機械傳動國家重點實驗室，重慶 400030；2. 重慶大學航空航天學院，重慶 400044)

非均質材料因其優(yōu)異的力學性能，被廣泛地應用于航空航天、軌道交通、新能源汽車等領域中，如新能源汽車儲氫載體普遍使用碳纖維增強環(huán)氧樹脂[1]，飛機起落架、機翼通常選擇鈦基復合材料[2]。然而，隨著機械零部件在高溫、重載、高速等極端環(huán)境中的運行，由位錯或裂紋引起的表界面失效越來越成為影響機械零部件服役壽命的關鍵因素[3?4]。位錯作為工程結構中的重要缺陷，其在微觀層面上的遷移變化，會極大地改變材料的宏觀力學性能[5]，在非均質機械零部件中，位錯在基體內的運動還會受到周圍夾雜物如微孔洞、微裂紋、析出相、晶界等因素的極大影響，研究位錯與周圍介質的這種相互作用是工程領域和力學領域的熱點問題之一[6?8]。

美國西北大學的Dundurs教授和Mura教授[9]首次推導出了平面問題下圓形雜質和刃型位錯相互作用能(也稱交互能)和位錯受力的解析解。隨后，Stagni和Lizzio[10]進一步考慮了自由表面對圓形雜質和刃型位錯的影響，并給出了橢圓雜質與刃型位錯之間的相互作用解[11]。他們的解[10?11]是以Laurent多項式給出的，但無窮級數的解在實際應用上則顯得較為不便[12]。Santare和Keer[12]的表達式中規(guī)避了無窮級數帶來的不便，他們利用Muskhelishvili的復變函數法[13]，對刃型位錯與橢圓形剛體雜質的相互作用問題進行了解析研究。應用該方法[13]，Chen等[14]進一步給出了刃型位錯與三層橢圓形雜質的彈性解。當雜質所在區(qū)域變成圓形孔洞時，Dai[15]進一步研究了半平面對刃型位錯的影響，他們的研究是對Dundurs和Mura[9]全平面問題下位錯與圓形孔洞相互作用解的有益補充。

需要指出的是，單個位錯與雜質的相互作用問題是理解材料形變和微結構演化機制的基礎問題，在眾多方面有著重要的應用。如借助離散位錯方法[16]，可以進一步將雜質與單個位錯的相互作用關系擴展到雜質與裂紋的相互作用問題[17?19]中去，這里的雜質-位錯解被作為了求解雜質-裂紋相互作用問題的格林函數[20?21]。另外，當考慮雜質的尺度效應時，可以分析納米夾雜物對位錯的影響[19?22]，相關的研究技術通常需要引入原子模擬[23 ? 24]。

從上述研究[9?12,14?15,19,22,25?27]中可以看出，即使是圓形(或橢圓形)雜質與刃型位錯的相互作用問題，其解析解的獲得都顯得頗為不易，而相關的表達式在數學上也顯得較為復雜。但是，在實際的工程材料中，夾雜物不僅具有很強的分布離散性[28]。而且其大小、形態(tài)往往都不一樣，這造成了夾雜物對周圍介質不同的應力擾動狀態(tài)[29?30]。Ehselby[31]所提出的等效夾雜理論為求解任意形狀夾雜物對彈性場的擾動效應提供了可能。等效夾雜法(Equivalent inclusion method，EIM)通過將非均質材料問題轉為均質材料疊加一個待定的“本征應變”問題，可以有效地統(tǒng)一解決塑性應變、殘余應力、材料相變等問題[32]。雖然Eshelby的等效夾雜法最初是用于求解三維橢球雜質問題的，但是正如Jin等[33]在其研究中所指出的那樣：Eshelby等效夾雜法同樣是求解平面應變或平面應力雜質問題的一種有效方法。

對Eshelby上述理論的數值實現稱為“數值等效夾雜法(NEIM)[34]”，即通過對雜質所在區(qū)域進行數值離散，相應單元上的等效本征應變通過迭代法來進行確定。Zhou等[34]應用數值等效夾雜法對遠端受載的平面雜質問題進行了研究并指出：對于任意Dundurs參數的材料組合，數值等效夾雜法均具有很好的數值穩(wěn)定性和魯棒性。應該指出，著名力學家Hutchinson[35]提出了一種近似法，用于研究裂紋尖端大量微裂紋的影響，尤其是微裂紋區(qū)域內應力強度的降低或屏蔽效應。由于Hutchinson的方法只考慮零階的局部擾動效應[35]而沒有進行迭代逼近，因而在文獻中[34]稱為零階迭代近似解。Hutchinson方法的這種思想也被Li等[36?37]應用于解決雜質和位錯的相互作用問題。但是，和精確解[9]相比。這種近似法具有較大的數值誤差，尤其是當計算點位于雜質附近區(qū)域時[34?38]。

位錯和夾雜物的交互能及位錯受力對于理解位錯運動及其引起的材料微結構變化具有相當重要的作用。現有文獻[9 ? 12, 14 ? 15, 19, 22, 25 ? 27,36 ? 37]對雜質-位錯交互能和位錯受力的解析解研究，主要探討特定的夾雜形狀，如圓形或橢圓形。而當用數值方法求解交互能時，需有效處理無窮域上積分，及位錯奇異性問題。近年來，Li等[38]在螺旋位錯與夾雜物的交互能研究上，取得了理論分析和數值算法的重要突破。在這項工作中，基于反平面剪切問題力學模型，他們[38]提出相應的等效夾雜算法并對圓形異質夾雜與螺型位錯的相互作用能進行了解析求解，獲得了交互能和彈性能的顯式單元解，在此基礎上，進一步提出了反平面剪切問題中求解任意形狀雜質和螺型位錯的數值計算方案。需要指出的是，這項研究首次報道了一種非多項式非指數形式的等效本征應變，并論證了等效夾雜法可以有效地用于解析求解非均勻本征應變。

然而，應用等效夾雜原理求解刃型位錯問題的有效性尚未得到驗證。一方面，刃型位錯的數學表達式相對螺型位錯而言更加復雜，加之周圍雜質對其產生的擾動效應，這使得對該問題交互能的解析求解變得異常困難；另一方面，由于控制方程中待求解的未知量陡然增多，即使是構建數值算法進行計算，也將不得不考量由位錯奇異性和材料錯配性所造成的數值波動和計算失穩(wěn)。本文所述及的夾雜物與刃型位錯的交互能研究，進一步將反平面計算模型[38]擴展到了平面問題中，通過引入Dundurs參數，推導了平面問題下的理論模型，并論證了等效夾雜法在求解刃型位錯問題中的有效性；同時，由于Dundurs參數的引入，本文數值算法相對前述論文[38]所提出的計算方法得到了較大改善，因本文算法極大地削減了迭代計算中的變量數目，這進一步提高了本文算法的魯棒性和穩(wěn)定性；此外，通過引入相對誤差的范數形式，本文定量分析了計算結果的相對誤差，并對本算法的有效性進行了評估。需要指出的是，三維問題下夾雜物與位錯的耦合計算往往較難獲得解析解，通常需要借助分子動力學模擬或者ab-initio理論。本文二維分析的計算，基于連續(xù)介質力學的理論框架，相應模型的計算結果在一定的宏觀尺度上可以揭示具體工程實際的物理本質。

本文借助等效夾雜法，建立了非均質材料中雜質與刃型位錯的應力控制方程，結合邊界條件、平衡方程和高斯定理，推導了交互能的計算公式。同時，本文結合快速傅里葉變換算法，給出了基于面單元離散的數值化計算方案，通過和相應的精確解對比，驗證了本文計算方案的有效性；并對本文計算方案和零次迭代解的無窮范數和2-范數進行了相對誤差分析，證明了本文計算方案在求解任意形狀雜質與位錯問題上的數值穩(wěn)定性和有效性。此外，本文還示例分析了任意雜質的形狀參數對交互能的影響。

1 任意形狀夾雜物與位錯交互能的數值計算方法

在無限大的平面D中，在(ξ,η)處放置有一個刃型位錯，在位錯附近分布有一個任意形狀的雜質Ω(圖1)。在本研究中，假定雜質和其周圍的基體是完美固結的，并且基體和雜質的彈性模量分別為Cijkl和kl。

圖1 任意形狀雜質與刃型位錯的相互作用示意圖Fig. 1 Schematic of an arbitrarily shaped inhomogeneity interacting with an edge dislocation

1.1 夾雜物-位錯問題中的等效夾雜法

根據Eshelby等效夾雜法的基本原理[29]，在刃型位錯作用下，由于雜質Ω存在所產生的應力擾動，可以通過一個分布有合適本征應變ε?ij的夾雜問題來模擬(圖2)。因此，圖1中雜質與位錯的相互作用應力解可以通過疊加一個均質材料解(即位錯解)和相應的夾雜解來得到(圖2)。

圖2 雜質-位錯問題中的Eshelby等效夾雜法Fig. 2 Schematic of Eshelby's equivalent inclusion method for the interaction between an edge dislocation and an inhomogeneity

式中，α和β為Dundurs參數，其具體表達式為：

本文以平面應力情況為例，相應解可以用來解釋平面應變的情況。上述式(4)中的κ為Kolosov常數，在平面應力問題中，κ=(3?ν)/(1+ν)；μ和ν分別表示材料的剪切模量和泊松比。需要指出的是，下標1和2分別表示基體和雜質。

1.2 夾雜物與位錯的交互能及位錯受力

在圖1所示的任意形狀雜質-位錯問題中，在x3方向上單位長度內的彈性能W可以看作是雜質對刃型位錯做功[38]：

式中：σij為單獨由雜質產生擾動場；為由位錯造成的位移梯度。

單獨由位錯在均質材料中引起的彈性能W0為：

因此，雜質與位錯之間的交互能ΔW 為：

為了進一步簡化上述式(7)，注意到有下列關系式：

式中，uk,l為雜質產生的位移梯度。

考慮自由邊界條件和平衡方程，并應用散度定理，式(7)還可以進一步推導成下述形式：

由于本征應變只作用在夾雜內部，而在夾雜外部為零。因此，式(9)將積分運算僅僅限定在了一個有限且確定的作用區(qū)域上，從而避免了如式(7)所示的在一個無窮大區(qū)域上的積分計算，這樣使得在數值計算過程中，大大降低了計算機存儲空間和不必要的時間浪費。

一旦獲得了位錯與雜質之間的交互能，可以進一步得到作用在位錯上的力[39]：

式中，i 和j為笛卡爾坐標系下的單位基矢量。需要指出的是，正向的位錯受力表示雜質和位錯之間是排斥關系，相反，位錯則受到雜質的吸引。

1.3 交互能的數值計算方法

從交互能的式(9)可以看出，求解的關鍵在于確定雜質-位錯問題中的等效本征應變以及相應的應力擾動場σij。圖3顯示了基于數值等效夾雜法求解交互能和位錯受力的計算流程圖。參數N為位錯的位置編號，K為迭代次數。當位錯的位置確定時，在整個循環(huán)流程中，位錯應力場σi保持不變，其具體表達式可由Hills等[16]專著得到。需要說明的是，在等效夾雜法的數值實現中，初始的應力擾動場輸入均為零，即K=1時，σ=0。在不斷迭代過程中，當迭代次數大于最大迭代次數Kmax，或者迭代精度小于程序設置誤差δ時，輸出最終的等效本征應變ε和相應的應力擾動場σ。

圖3 雜質-位錯問題中的交互能和位錯受力的計算流程圖Fig. 3 Flow-chart of the present computational scheme for solving the interaction energy and force on dislocation due to an edge dislocation interacting with an inhomogeneity

在上述計算過程中，夾雜應力場的耗時占據了整個流程的大部分時間。根據疊加法的基本原理，夾雜應力場可以表示為如下公式：

式中：Nx、Ny分別為夾雜區(qū)域Ω上沿x和y方向的離散單元數；(I,J)和(I′,J′)分別表示目標場和激勵場的單元編號。這里的Tijkl為影響系數 (其矩陣形式見式(3))，相關的具體解析表達式見文獻[40]。

式(11)展示了一類卷積型疊加，利用快速傅里葉變換算法的加速性質，可以顯著地提高其運算速度，具體操作過程簡述如下[41?43]：首先，將影響系數矩陣和本征應變所代表的激勵矩陣擴展到(2Nx×2Ny)的區(qū)域上；接著分別對上述兩個矩陣完成包卷循環(huán)填充(wrap-around order padding)和補零(zero-padding)操作；經過離散傅里葉變換后，將上述兩個序列在傅里葉頻域內進行點對點相乘；通過對上述乘積進行傅里葉逆變換，拋去擴展區(qū)域上的無效數值，即得到傅里葉時域內的夾雜擾動應力解。對比直接疊加法和離散傅里葉變換算法可知，式(11)所需的時間復雜度為O(9Nx2Ny2)，其中，O(·)表示乘法運算的階數，相應地，經過傅里葉變換算法加速后，所需的時間復雜度為O(108NxNylog2(4NxNy))。由此可知，當網格數取32 ×32時，直接疊加法所消耗的時間將是傅里葉變換算法所消耗時間的近24倍。隨著網格數的增多，傅里葉變換算法所取得的時間加速將更加明顯。需要指出的是，快速傅里葉變換算法是對離散卷積的精確求和[41?43]，可以顯著地提高計算速度；由于，本文采用了單元離散-疊加的方法，因此，可以用于求解任意形狀雜質，包括邊界為平滑或者非平滑截面的夾雜物[34,44]。

為了進一步獲得交互能的數值解，對等效夾雜所在的區(qū)域進行矩形單元離散，式(9)可以進一步：

在獲得交互能的基礎上，通過中心差分法即可得到作用在位錯上的力。

2 算例分析

2.1 數值化計算方法的驗證

為了驗證本文計算方法的有效性，考慮在無限大的平面介質中分布有一個半徑為a的圓形夾雜物(圖4)，其剪切模量和泊松比分別為μ2和ν2。在雜質外(2a, 0)處作用有一刃型位錯，相應的基體彈性模量和泊松比分別為μ1和ν1。表1列出了相關計算參數的數值大小，其中刃型位錯的Burgers矢量bi取值大小為(1, 0)，其符號的方向定義為：沿著位錯軌跡從位錯核走向無窮遠處，軌跡的右側為bi的正向，左側為bi的負向。

圖4 圓形雜質與刃型位錯的相互作用示意圖Fig. 4 Schematic of a circular inhomogeneity interacting with an edge dislocation

表1 夾雜物和相應基體的計算參數Table 1 Computational parameters of an inhomogeneity and the corresponding matrix

本文算法采用了迭代次數K=5的數值計算方案，其中，在面單元所離散的夾雜區(qū)域上，單元數取Nx=32,Ny=32。圖5、圖6中的零次迭代解表示在計算過程中忽略夾雜擾動場所得到的數值結果。本文的解析解來自于Dundurs和Mura教授[9]的計算結果，其中的交互能可以從總能量中除去位錯自身能量得到。

圖5 無量綱交互能的變化圖，其中位錯位置ξ/a在x1軸上從1.1移動到1.6Fig. 5 Variation of the normalized interaction energy with dislocation position varying along x1-axis from 1.1 to 1.6

圖6 無量綱位錯受力的變化圖，其中位錯位置ξ/a在x1軸上從1.1移動到1.6Fig. 6 Variation of the normalized force with dislocation position varying along x1-axis from 1.1 to 1.6

在圖5和圖6中，分別對夾雜物(SiC和Ti-6Al-V)和位錯之間的交互能和位錯受力進行了無量綱化研究，其中，ΔW0=μ1/[2π(k1+1)]，F0=μ1/[πa(k1+1)]。隨著刃型位錯位置不斷地靠近雜質，二者之間的交互能絕對值越來越大，相應地，位錯受到雜質的排斥效應也越來越劇烈。從上述兩幅圖中還可以明顯看出，相對于零次迭代解，本文計算方案與準確解具有更好的數值吻合度。

在文獻[45]中，給出了在刃型位錯(位置在(3.0×10?3mm, 0))周邊存在一個圓形的非金屬氧化物Al2O3(如圖4所示)，其剪切模量為152 GPa，泊松比為0.22，半徑為2.0×10?3mm，伯格斯矢量的大小為(1.0×10?7mm, 0)。通過本文的計算，此時非金屬夾雜物與刃型位錯的交互能和位錯受力如圖7和圖8所示。

圖7 刃型位錯與Al2O3交互能的變化圖，其中位錯位置ξ/a在x1軸上從2.1移動到2.6Fig. 7 Variation of the interaction energy between edge dislocation and Al2O3 where dislocation position varies along x1-axis from 2.1 to 2.6

圖8 位錯受力的變化圖，其中位錯位置ξ/a在x1軸上從2.1移動到2.6Fig. 8 Variation of the force on edge dislocation with dislocation position ξ/a varying along x1-axis from 2.1 to 2.6

2.2 本文計算方法的可靠性分析

為了更加準確地分析本文計算方案與準確解之間的相對誤差，引入無窮范數和2-范數，對零次迭代解和本文數值方法所得到的交互能和位錯受力進行誤差分析，相關的對比結果見表2。需要說明的是，計算中選取了位錯位置ξ/a從1.1到1.5上均勻分布的50個點進行了計算。從表格中可以看出，零次迭代解在無窮范數和2-范數上均具有較大的相對誤差，而采用本文計算方案后，相對誤差的范數均得到極大地減小，說明本文計算方案在求解過程中具有很好的數值收斂性和穩(wěn)定性。結合圖5、圖6及表2可以進一步看出，對于只考慮材料解的零次迭代法，當計算點位于雜質附近區(qū)域時，所得到的計算誤差尤為明顯。

表2 本文數值計算方案和零次迭代解之間的相對誤差結果對比Table 2 Comparisons on relative error of results between the present numerical scheme and 0th solution

2.3 SiC或Ti-6Al-V雜質對位錯能量及受力的影響

在航空結構中，通常添加碳化硅或鈦合金顆粒來增強重要部件的機械強度，但是相關研究也表明，微裂紋的萌生通常會出現在夾雜物的周邊。為了進一步說明本文計算方法對于求解任意形狀雜質與刃型位錯交互能的適用性，對圖9所示的不同形狀的橢圓SiC雜質進行了分析，其中橢圓的兩個短半軸比a2/a1分別取0.25、0.5、1.0、2.0、4.0，其中，a1和a2的方向分別沿x1和x2。橢圓形SiC雜質和基體的材料常數采用如表1所示的計算參數。從圖9可以看出，隨著刃型位錯沿x1軸移動并不斷靠近雜質，二者之間的交互能逐漸增大，且曲線的斜率在靠近雜質處漸漸變陡；當位錯位置固定時，隨著a2/a1的比值變大，交互能則逐步增加，但增加趨勢則逐漸變小。

圖9 橢圓形SiC雜質與刃型位錯的交互能Fig. 9 Schematic of an elliptical SiC inhomogeneity interacting with an edge dislocation

已有的解析解研究通常只能針對特定的夾雜物形狀展開分析，如標準的橢圓形或圓形，但是實際工程中的夾雜物更多是以任意形狀方式進行分布的，此時的解析解在工程應用中則受到了極大地限制。在圖10和圖11中，本文選擇了由雜質形狀函數r(m)(ψ)=1+2ancos(mψ)控制的參數方程，其中，an為形狀系數，m為粗糙度系數。在本算例中，形狀系數取0.03，粗糙度系數分別取5和8。本文圖10和圖11則對這種形狀的SiC雜質和Ti-6Al-V雜質與刃型位錯之間的交互能和位錯受力進行了研究。

圖10 復雜形狀雜質與刃型位錯的交互能Fig. 10 Schematic of the interaction energy between an inhomogeneity with complex boundary and an edge dislocation

圖11 復雜形狀雜質與刃型位錯作用下的位錯受力Fig. 11 Schematic of the force due to an inhomogeneity with complex boundary and an edge dislocation

從圖中可以看出，比基體較硬的SiC雜質對刃型位錯產生的是排斥力，而較軟的鈦鋁合金Ti-6Al-V則會吸引刃型位錯。更為一般的，當夾雜物區(qū)域變?yōu)榭锥磿r，即彈性模量為零，此時孔洞對其周邊位錯產生的吸引力會進一步導致位錯在材料中的湮滅，這種理論所預測的現象已在實驗中被證實[46]。此外，對于同一雜質與刃型位錯而言，不同形狀雜質下產生的交互能和位錯受力大為不同，且這種影響在雜質附近處表現尤為明顯。

3 結論

本文主要研究了任意形狀夾雜物與刃型位錯之間的相互作用機理及其數值化計算方法，并得出以下主要結論：

(1) 通過數值等效夾雜法可以建立求解任意形狀雜質與位錯之間的控制方程，將非均質材料的非線性問題轉化為線性問題。

(2) 利用能量法、高斯定理和等效夾雜原理，論證了無窮體內夾雜與刃型位錯的交互能可以用夾雜區(qū)域上的積分來表示，避免了在無窮區(qū)域上的數值計算。

(3) 相比于零次迭代解，本文計算方案具有較好的數值收斂性和穩(wěn)定性。

(4) 本文關注夾雜物與刃型位錯的交互能及位錯受力，文中圖5和圖6分別定量地展示了它們隨夾雜物位置變化規(guī)律。