多晶體壓剪試樣靜態加載有限元計算

趙偉業，趙 聃，呂 品，金 濤，馬勝國

（1.太原理工大學機械與運載工程學院應用力學研究所，山西 太原 030024；2.材料強度與結構沖擊山西省重點實驗室，山西 太原 030024；3.太原理工大學力學國家級實驗教學示范中心，山西 太原 030024）

多晶體是由許多不同取向的單晶體構成，晶粒的取向對金屬材料的塑性變形特征以及損傷、斷裂等都有很大影響[1]。從微觀機制上看，晶體塑性變形方式有位錯滑移、孿生、相變等多種方式，晶體塑性理論主要分析上述機制影響下材料的變形規律[2]。由于各晶粒的取向不同，塑性變形首先發生在位錯滑移容易開動的晶粒上，且晶粒間的相互作用協調各自的變形，導致多晶體內部的應變場和應力場具有強烈的不均勻性。晶體塑性有限元法是一種重要的分析晶粒間相互作用的方法[3]。國內一些學者基于多晶體材料的特點建立了Voronoi有限元模型，針對不同工況并結合材料自身性質，對多晶體模型進行了數值分析[4-6]。

在材料成型過程中，經常出現復雜的受力狀態，須明晰復雜應力狀態下材料的強度和變形特性，其中壓剪復合加載導致的復合應力狀態是國內外學者的關注點。鄭文等[7]通過添加具有不同角度傾斜端面的墊塊，研究了靜態和動態壓剪復合加載下材料的力學響應，給出了相應的數據處理方法，并利用有限元模擬驗證了實驗的可行性。章超等[8]基于鄭文的實驗設計方法，利用材料試驗系統（MTS）靜態實驗裝置和分離式霍普金森壓桿（SHPB)動態實驗裝置，分別進行了不同角度斜端面墊塊下花崗巖的準靜態與動態加載試驗，分析了不同加載路徑下巖石的力學性能以及破壞情況。李雪艷等[9]通過添加單斜端面墊塊和長方體套筒，利用材料試驗機對閉孔泡沫鋁進行準靜態條件下的壓剪復合加載實驗，并通過改變墊塊的角度得到泡沫鋁在不同加載路徑下的力學性能。國外Rittel等[10-11]提出帶有兩個與豎直方向呈45°的對稱矩形斜槽的圓柱體壓剪試樣，用于研究較寬應變率及應變范圍內試樣的力學響應。Dorogoy等[12-15]利用準靜態和動態兩種加載方式，結合有限元模擬探究了壓剪試樣在較寬應變率范圍內發生塑性應變時的力學性能，測試了不同應變率下材料的剪切力學響應，并將壓剪試樣的矩形斜槽變為半圓形斜槽，改善了壓剪試樣矩形斜槽邊角的應力集中現象。Vural等[16]對壓剪試樣進行了數值分析和實驗測試，建立了外部施加載荷、位移與斜截面內的等效應力、等效塑性應變等之間的關系式，在此基礎上，為了增強斜截面應力場的均勻性，改變斜槽的角度并進行分析。Zhao等[17]將金屬壓剪試樣改變為斜槽部分由聚合物代替的全新壓剪試樣，通過動態及準靜態實驗結合數值模擬，進行應力、應變分析，并將這種新試樣與原來的全金屬壓剪試樣進行全面對比。

上述研究大都假設材料為各向同性，未能體現出微觀晶粒對宏觀各向異性行為的影響。本研究引入晶體塑性理論，試圖揭示多晶體材料的微觀變形機理，以期提高材料加工成型產品的質量。將多晶體材料屬性賦予壓剪試樣，結合材料的宏、微觀特性進行有限元分析，并對壓剪復合應力狀態下微觀晶粒的變形演化進行分析。

1 晶體塑性理論

1.1 運動學

晶體塑性理論的基本框架基于Perice、Asaro和Needlman[18-19]的研究工作。材料有限變形過程中，總體變形梯度F可分解為

式中：F?為彈性變形梯度，包含彈性變形和剛性轉動；Fp是體現晶體滑移的塑性變形梯度。速度梯度的表達式為

式中：L?為速度梯度的彈性部分，Lp為速度梯度的塑性部分。

其中張量D和 ? 也可以分解為弾性變形部分和塑性變形部分（上標*和p分別表示彈性和塑性），即

并且滿足

1.2 本構模型

在晶體變形過程中，彈性變形率D?和Cauchy應力張量的Jaumann率表示為

式中：σ為Cauchy應力，I為二階單位張量，C為四階對稱彈性模量張量。

分切應力 τ(α)與Cauchy應力之間有如下關系

本研究以面心立方（FCC）晶體為研究對象。面心立方晶體包含4個滑移面，每個滑移面包含3個滑移方向，因此FCC晶體共有{111}〈110〉 12個滑移系。

1.3 率相關硬化模型

由Pierce等[18]提出的以Schmid準則為基礎的率相關剪切應變率模型為

應變硬化通過臨界分切應力的增量關系演化表示為

式中：hαβ為硬化矩陣。hαα（α = β）表示自硬化模量，hαβ（α ≠ β）表示潛硬化模量。

晶體變形時，不同硬化階段的硬化矩陣不同，準確表達的難度很大，所以通常需要對硬化矩陣進行簡化。對自硬化模量可做簡單假設[20]

式中：h0為初始硬化模量； τ0為臨界切應力；τs為飽和切應力；γ為所有滑移系上切應變的Taylor累積積分，即總切應變，表示為

潛硬化模量表示為

式中：q為硬化矩陣系數。

2 有限元數值模擬

2.1 材料初始織構對數值計算的影響

以6061鋁合金棒材為研究對象進行數值計算?？紤]到棒材軸向拉拔后會帶有絲織構[21]，將對數值結果產生影響，為此使用Abaqus有限元軟件和Python語言建立了包含54個晶粒的Voronoi代表性體積單元（Representative volume element, RVE）模型進行壓縮計算，模型分別包含了初始隨機織構和初始絲織構，其他條件完全相同。

圖1給出了多晶體的兩種初始取向，即隨機織構和絲織構在（111）面上的極圖[22]，利用表1所示的6061鋁合金的材料參數[23]進行數值計算，單向壓縮時的應力-應變（σ33-ε33）曲線計算結果如圖2所示。表1中，C11、C12、C44為彈性模量分量，m為率敏感指數，q為硬化矩陣系數。

圖1 多晶體隨機織構和絲織構的極圖Fig.1 Polar graphs of polycrystalline with random texture and fiber texture

表1 6061鋁合金材料參數[23]Table 1 Material parameters of 6061 aluminum alloy[23]

由圖2可知：當塑性變形較小時，初始織構為絲織構的模型在壓縮方向上的應力絕對值高于隨機取向對應模型，最大相對偏差約為7%；當塑性變形增加到一定程度時，兩種模型對應的曲線基本重合。這表明當塑性變形較小時，初始織構為絲織構的模型受影響較為明顯；隨著塑性變形的增大，多晶體的晶體取向發生了一定程度的旋轉，使得塑性變形較大時兩種模型的應力比較接近。在后續的有限元計算中，均以具有初始絲織構的材料為研究對象。

圖2 具有不同初始織構的單向壓縮應力-應變曲線Fig.2 Unidirectional compression stress-strain curves with different initial textures

2.2 創建多晶體壓剪模型

使用Abaqus有限元軟件和Python語言建立多晶體模型，其中1個單元代表1個晶粒，創建兩個平板剛體模型，平板剛體與多晶體模型上下底面進行接觸設置，多晶體壓剪試樣幾何尺寸及模型如圖3和圖4所示。接觸屬性為面面硬接觸，并且設置相應的摩擦系數。給予下板參考點全約束邊界條件；上板參考點處設置一個2.5 mm的壓縮位移載荷，其他邊界條件全約束。

2.3 摩擦力對數值模擬結果的影響

圖3 多晶體壓剪試樣幾何尺寸Fig.3 Geometrical dimension of polycrystal shear-compression specimen

考慮到壓剪模型與剛體接觸板之間摩擦的影響，首先對摩擦系數μ進行分析。根據齊康等[24]的總結，鋁合金在良好潤滑情況下摩擦系數可達到0.02左右。于是設置3組摩擦系數進行數值模擬，分別為0.025、0.050、0.100。每組摩擦系數對應的Mises應力分布如圖5所示，可以看到試件內應力分布不均勻，隨著摩擦系數的增加，斜槽以外部分的應力明顯增加，且整個模型中間向外膨脹的變形更明顯。每組模型變形后上板都會與壓剪模型發生相對滑動，如黑色箭頭所示。相對滑動距離分別為1.282、0.688、0.214 mm，其中摩擦系數最小組的相對滑動距離幾乎是摩擦系數最大組的6倍，對斜槽及斜槽以外部分的交界面的變形產生了明顯影響。

圖4 多晶體壓剪模型Fig.4 Polycrystal shear-compression model

圖6和圖7分別為模型斜槽部分單元的平均切應力-平均切應變曲線和平均正應力-平均正應變曲線。由兩圖可知，摩擦系數對切應變和正應變有明顯影響，且最大切應變（見圖6插圖）和最大正應變（見圖7插圖）都與摩擦系數呈負相關。切應變的最大相對偏差為13.7%，正應變的最大相對偏差為56.1%。圖8為模型斜槽部分對應單元的平均Mises應力-平均等效應變曲線，其最大等效應變隨摩擦系數的變化規律與圖6、圖7相似，由插圖知不同摩擦系數下等效應變的最大相對偏差為34.5%。圖9為模型頂面的位移-載荷曲線，當模型壓縮2.5 mm時，不同摩擦系數下所需載荷的最大相對偏差為12.4%，即摩擦系數越大，壓縮模型時需要的載荷越大。對不同摩擦系數情況下開槽部位的最大剪切應變、最大正應變以及等效應變進行總結。結果列于表2。

圖5 不同摩擦系數情況下的應力云圖Fig.5 Stress nephogram with different friction coefficients

圖6 不同摩擦系數情況下模型斜槽單元的平均切應力-平均切應變曲線Fig.6 Average shear stress-average shear strain curve of the model's chute element corresponding to different friction coefficients

圖7 不同摩擦系數情況下模型斜槽單元的平均正應力-平均正應變曲線Fig.7 Average normal stress-average normal strain curve of the model's chute element corresponding to different friction coefficients

圖8 摩擦系數不同時模型斜槽單元的平均Mises應力-平均等效應變曲線Fig.8 Average Mises stress-average equivalent strain curve of the model's chute element corresponding to different friction coefficients

圖9 摩擦系數不同時模型頂面的力-位移曲線Fig.9 Force-displacement curve of the top surface of the model corresponding to different friction coefficients

表2 不同摩擦系數對應模型的數值結果比較Table 2 Numerical results of models corresponding to different friction coefficients

以上結果表明：隨著摩擦系數的增加，模型與剛體接觸板之間的滑動阻力增大；當3組模型都壓縮2.5 mm時，摩擦系數的增加使模型與剛體接觸板的相對滑動減小，斜槽部分沿切向和正向的變形都減小。因此，斜槽部分沿切向和正向的變形在豎直方向的分量之和減小，而斜槽以外部分在豎直方向上的壓縮量增加。此分析結果與Dorogoy等[13]關于摩擦力的分析相吻合。因此，在模型壓縮2.5 mm時，為了減少所需壓力，增加斜槽的壓剪變形程度，應當盡可能地減小摩擦系數。

為了進一步分析壓剪模型關鍵區域的微觀變形特征，在斜槽部分按照圖10箭頭所示方向選取了16個單元作為特征晶粒。各晶粒的初始歐拉角 (φ1,ψ,φ2)如表3 所示。

圖10 模型斜槽部分一條棱上所選取的16個單元Fig.10 16 selected elements on an edge of the chute part of the model

表3 特征晶粒歐拉角Table 3 Euler angle of characteristic grains

圖11顯示了不同摩擦系數下16個晶粒在z方向上的應變變化?？梢钥闯?，大多數晶粒即使位置接近，但其應變值差異也較為明顯，這是因為初始晶粒取向差較明顯。圖11中不同摩擦系數下晶粒1和晶粒16的應變絕對值較其他晶粒偏小，而且隨著摩擦系數的增加而減小，原因是摩擦系數增加使模型中間部分產生更明顯的漲型，斜槽兩端在z方向上的變形減小，1號和16號晶粒位于斜槽端部，變形過程中會受到端末效應的影響，影響其變形。

圖11 不同摩擦系數下特征晶粒在壓縮方向的應變Fig.11 Strain of characteristic grain in compression direction under different friction coefficients

2.4 晶粒尺寸對數值結果的影響

多晶體有限元計算時，為了考慮晶粒尺寸對數值結果造成的影響，在單元數目不變的情況下建立了如圖12所示的包含不同晶粒數目的多晶體計算模型，通過晶粒數目的差異體現晶粒尺寸的影響，模型上下剛體接觸板與模型之間的摩擦系數設置為0.025，其他條件不變。

圖13為不同晶粒數目模型的位移-載荷曲線，同時與1個單元表示1個晶粒的模型計算結果（Number of grains 20 000）進行對比。由圖13可知，雖然晶粒數目不同，但其體現的初始拉伸織構基本一致，且幾種模型計算的載荷結果在變化趨勢上基本一致，各曲線間相差不足1%。

圖12 包含不同晶粒數目的壓剪有限元模型Fig.12 Finite element models of compression shear with different grain numbers

圖13 不同晶粒數目對應模型的載荷-位移曲線Fig.13 Force-displacement curves of models corresponding to different numbers of grains

2.5 單元數目對數值結果的影響

為了探究單元數對多晶體壓剪模型斜槽區域力學性能的影響，控制模型的晶粒數不變（固定為500），建立了3組單元數分別為8 026、14 604、23 835的模型進行數值分析，摩擦系數設置為0.025，單元類型為C3D8，其他條件不變，見圖14。

圖15為不同單元數目的3組模型的變形Mises應力云圖，3組模型的變形情況基本一致。由于單元數目的變化，單元尺寸也隨之改變，導致建模過程中，3組模型處于同一位置上的晶粒根據臨近原則所包含的單元數目存在一定的差異，所以變形云圖也產生了輕微的差異。

圖14 單元數目不同的模型Fig.14 Models with different number of elements

圖16為3組模型同時壓縮2.5 mm時剛體接觸板上參考點的位移-載荷曲線?？梢?，3組曲線軌跡基本相似，但是在相同位移下，載荷大小關于單元數目沒有明顯的規律，當模型壓縮至2.5 mm時，3組模型載荷的最大相對偏差約6.97%，如表4所示。圖17為3組模型斜槽部分相同位置處晶粒的Mises應力-等效應變曲線。由圖17插圖可見，單元尺寸及數目的變化導致3個模型的晶粒形狀有一定差異，Mises應力云圖的差異比較明顯；3條曲線的彈性部分基本一致，塑性部分的差異較明顯。以上表明，當單元數目或尺寸發生變化時，3組模型相同位置處晶粒的形狀受到一定程度的影響，原因是指定位置的晶粒與周圍晶粒間的相互作用發生了改變，周圍晶粒的取向和形貌發生變化，進而影響該位置晶粒的應力和應變。

圖15 單元數目不同時模型的Mises變形云圖Fig.15 Mises deformation nephograms of models with different numbers of elements

圖16 單元數目不同時模型的載荷-位移曲線Fig.16 Force-displacement curves of models with different numbers of elements

圖17 斜槽相同位置處晶粒的Mises應力-等效應變曲線Fig.17 Mises stress-equivalent strain curve of grains at the same position in the chute

表4 單元數目不同時模型在不同壓縮距離下的載荷及其最大相對偏差Table 4 Loads of the model with different numbers of elements at different compression distance and their maximum relative differences

2.6 單元類型對模擬結果的影響

為了探究單元類型的影響并提高模型計算效率，將斜槽區域單元數為2 093的四面體單元模型與斜槽區域單元數為2 040的六面體單元模型進行對比分析，兩種模型除了單元類型外其他條件與前述相同。斜槽局部變形如圖18所示，兩種網格模型的斜槽整體變形相似，應力分布都不均勻，從微觀的角度看六面體單元變形程度更大，且由云圖顏色分布可知六面體單元的計算精度更高，表明六面體單元比四面體單元的剛度更小，更適合分析模型局部變形情況。

圖18 兩種單元類型對應的模型斜槽部分Mises變形云圖Fig.18 Mises deformation nephogram of the model with two element types

圖19、圖20分別是兩種不同類型單元的斜槽的位移-載荷曲線和Mises應力-平均等效應變曲線?？梢?，當模型被壓縮到一定程度時，四面體單元的應力更大，位移為2.5 mm時，兩條曲線的應力偏差為3.4%。六面體單元模型的最大平均等效應變比四面體單元模型大4.2%。這表明當整個壓剪模型壓縮2.5 mm時，六面體單元模型在承受更小壓力的同時，模型斜槽部分單元發生了更大的變形。

圖19 兩種單元類型對應的模型頂面力-位移曲線Fig.19 Top surface force-displacement curve of model with two element types

圖20 兩種單元類型對應的模型斜槽部分單元的平均Mises應力-平均等效應變曲線Fig.20 Average Mises stress-strain curves of the sloped part of the model with two element types

3 結 論

(1)塑性變形較小時，初始織構為絲織構的RVE單向壓縮模型相對于晶粒為隨機取向的RVE單向壓縮模型的應力值更高，當塑性變形增加到一定程度時，兩種模型的應力值接近。這表明模型變形時晶粒取向會發生一定程度的轉動，從而對數值結果造成影響。

(2)數值計算了3種不同摩擦系數對多晶壓剪模型變形的影響。結果表明，隨著摩擦系數的增加，模型與剛體接觸板相對滑動減小，模型中間部分發生明顯漲型，模型斜槽部分的壓縮與剪切變形都減小，使得在壓縮相同高度時，摩擦系數大的模型的斜槽部分沿豎直方向的壓縮量減小，而斜槽以外部分沿豎直方向的壓縮量增加，同時壓縮模型所需要的力也增加。數值模擬結果顯示，選取摩擦系數為0.025進行計算能夠得到較好的壓剪應力狀態。從微觀角度分析，由于晶粒初始取向不同導致Mises應力分布不均勻，且對于位置接近的晶粒，晶粒初始取向差的存在會導致其應變值差異明顯。

(3)在相同摩擦條件下，分析固定單元數目而改變晶粒數目（尺寸）、固定晶粒數目而改變單元數目兩種情況，并討論了單元類型對數值結果的影響。結果表明：在單元數目確定的情況下，晶粒數目增加到一定程度時（1個單元代表1個晶粒），數值結果并未出現明顯差異；而固定晶粒數目、改變單元數目時，由于單元尺寸的變化導致不同組模型、相同位置的晶粒形狀有所差異，晶粒間的相互作用出現差異，從而對數值結果產生了一些影響。不同單元類型的計算結果表明，六面體單元的變形程度大于四面體單元，計算精度更高。在相同壓縮量下六面體單元模型所需壓力比四面體單元小，但計算時間較長。若考慮計算成本，可選擇單元數目相對較少的四面體單元對多晶體宏觀力學響應進行預估。