金屬材料在沖擊荷載下局部變形的數(shù)值模擬及分析

陳小翠， 杜成斌， 江守燕(河海大學 力學與材料學院， 南京 211106)

當金屬材料受高應變率載荷作用時，如爆炸、子彈射擊、高速撞擊、切削等，會產生集中于很窄的帶狀區(qū)域內的嚴重塑性變形，通常叫做絕熱剪切帶[1]。剪切帶寬度一般為10 μm量級，而且傳播的速度很快，為微秒量級[2-3]。金屬扭轉、拉伸實驗中經常觀察到帶內產生劇烈的塑性變形，且?guī)葴囟雀哌_100 °C量級。這些因素對材料造成不可逆損傷，如孔洞和斷裂的生長[4-5]，甚至造成材料的斷裂破壞[6]。

高應變率作用下，塑性變形通常使用與應變率、應力和溫度相關的本構方程來描述，塑性變形產生的熱能具有軟化效應，促使塑性流動產生，而應變率和應力的增加則產生材料硬化效應。Marchand和Duffy的實驗指出[7]，剪切變形局部化形成過程主要分為3個階段：第一階段主要是局部化發(fā)生之前，材料從彈性階段發(fā)展到均勻塑性變形階段；第二階段材料溫度快速增長，塑性變形增大，應力應變曲線出現(xiàn)應力峰值，材料開始應變軟化；第三階段即材料失穩(wěn)階段，材料應力極限下降，產生嚴重的材料局部化，導致金屬材料承載能力下降。

剪切帶現(xiàn)象不僅是多物理場的熱力耦合問題，更具有微米級的空間尺寸和微秒級的時間尺度，準確模擬剪切帶的產生和生長一直是研究的難題[8]?？紤]到問題的復雜性，在模擬過程中人們一般采用絕熱剪切的假設來簡化模擬過程中的數(shù)值計算難度[9-11]。但該假設忽略了熱擴散效應，即不考慮局部化形成過程中熱傳導與熱能之間的相互作用，這樣通常伴隨著網(wǎng)格敏感性[12]。

白以龍等[13-14]提出了熱塑剪切帶問題，強調耦合熱傳導和塑性功效應的重要性，并基于假設計算出一維熱塑剪切帶問題的準定常近似解，將剪切帶寬度的形成看作是塑性功和熱擴散相互平衡的產物。

本文基于有限元分析軟件FEAP (Finite Element Analysis Program)[15]，用Fortran語言編寫一個新的單元。該單元在能量方程中考慮了熱擴散作用，基于混合有限元的理論框架，同時將位移、溫度、應力和等效塑性應變場作為未知量，同步進行求解適用于高應變率作用下塑性變形局部化問題的偏微分方程組，適用于模擬二維塑性變形局部化問題。文中基于該單元進行了經典的45°剪切帶問題模擬，并將計算結果與絕熱假設結果比較，觀查網(wǎng)格敏感性。文中同時考慮了使用顯式算法和隱式算法進行剪切帶數(shù)值模擬，并將兩者的計算精度及計算成本進行對比。

1 控制方程

沖擊荷載作用下的剪切帶問題的模擬用偏微分方程組(PDEs)來表示，主要為動量守恒方程、能量守恒方程、彈性本構和塑性本構關系。

不考慮體力作用，動量守恒方程表示為

(1)

式中：u為位移向量，σ為應力張量，ρ為材料密度。

(2)

式中：T為溫度場，κ和cp分別為熱導率和比熱容，Dp為塑性變形率。

基于小變形假設

(3)

(4)

式中：v為速度向量；L為速度梯度張量。

根據(jù)彈塑性力學理論，考慮溫度效應的熱變形，則在彈塑性狀態(tài)下本構方程中變形率由三部分組成：彈性應變率Del、 溫度變化引起的應變率DT和非彈性應變率Dp，即

D=Del+DT+Dp

(5)

根據(jù)胡克定律(Hooke’s Law)，彈性本構方程為

(6)

剪切帶形成過程中伴隨著應變硬化、應變率硬化及熱軟化等作用，假設屈服應力與應力、溫度及塑性變形相關，本文采用修正Litonski[16-17]塑性本構方程

(7)

(8)

(9)

上式表明塑性應變張量與偏應力張量方向相同。

由于

(10)

(11)

式中：“:”指張量的雙點積運算。

根據(jù)式(10)和式(11),控制方程簡化為

(12)

邊界條件為

(13)

2 離散模型

PDEs的弱形式通過將各個方程乘以各自的試探函數(shù)：ωu，ωT，ωσ和ωγp，并在求解域Ω上進行積分。再根據(jù)分部積分原理，對偏微分方程降階。則R表示為

(14)

(15)

(16)

3 時間離散模型

3.1 隱式算法

文中在建立本構方程時采用率形式，因此在計算過程中可直接采用速度場v來代替位移場u作為變量，在時間離散過程中只需要變量的一階導數(shù)項，相對于隱式積分，可采用后向歐拉迭代法(Backward Euler Method)。

若對變量α采用隱式時間離散方法，則變量α和β均采用隱式算法。對于隱式算法，在每一個時間增量步都需要對非線性方程進行線性化并迭代求解。

R表示為在xn附近的Taylor展開式，僅保留線性項，在xn增量方向通過Gateaux求導，得出Jacobian矩陣，如下式所示

(17)

式中：J為雅可比矩陣；θ為無窮小量。故離散后方程可以表示為

(18)

式中：

(19)

(20)

(21)

(22)

因此式(18)可簡化為

(23)

回顧上節(jié)的Schur補方法，上式簡化為

(24)

用數(shù)學方法展開

(25)

3.2 半顯式算法

倘若對變量α采用前向歐拉方法(Forward Euler Method)，即顯式時間離散方法，而變量β采用隱式算法，則該模型采用半顯式算法離散。對于變量α計算過程如下

(26)

一直以來，顯式算法經常被用于求解需要分許多時間增量步的高速動力學問題，如沖擊、碰撞、爆破之類高度非線性問題。但顯式算法時間步長的選擇直接關系到數(shù)值模擬的穩(wěn)定性及計算時間。對于剪切帶模擬過程中時間增量的選擇，文中假設材料變形一直處于彈性階段，則變量α表示的動量守恒方程和能量守恒方程需要滿足不同的時間步要求，分別為力學問題和熱問題所需要滿足的時間步條件。

對于力學問題，臨界時間步為膨脹波傳播過程中通過最小網(wǎng)格單元所需的時間[20]，表示如下

(27)

式中：h為單元尺寸；c為P波速；λ和μ為Lame常數(shù)。

對于熱效應問題，臨界時間步為熱波通過最小網(wǎng)格單元所需的時間，表達如下

Δtth=min(ρcph2/2κ)

(28)

因此，對于上述模型，臨界時間步長為

Δt=γmin(Δtth，Δtmech)

(29)

式中：γ為考慮非線性問題的折減系數(shù)。

4 數(shù)值算例

4.1 算例1

本算例考慮上下兩端豎向受拉的金屬材料方板。板的具體受力狀態(tài)及速度邊界條件，如圖1所示。方板中心設有一個圓形缺陷，來誘發(fā)剪切帶的產生(圖1中圓形區(qū)域內為缺陷)。模型及其受力狀態(tài)關于x軸與y軸對稱，取右上角1/4方板進行模擬，在1/4方板左側及底部施加對稱邊界約束。

缺陷內屈服應力和屈服應變值由二維β函數(shù)折減

(30)

式中：α為折減系數(shù)；(x0,y0)缺陷中心坐標；r0為折減半徑。本算例中取α=4%，x0=0，y0=0，r0=1 μm，即半徑為1 μm范圍內材料屈服應力和屈服應變均折減4%，來誘發(fā)剪切帶的產生。

數(shù)值計算時，假設平面應變狀態(tài)，板的寬度為100 μm，即H=50 μm。施加的速度邊界條件，如圖2所示。1 μs內速度從0線性增長到5 m/s，隨后保持不變，即vr=5 m/s，tr=1 μs。本算例采用的材料參數(shù)及材料相應的塑性本構模型的參數(shù)詳見表1。

圖1 方形板結構及邊界條件(虛線區(qū)域為缺陷位置)Fig.1 Square plate configuration and boundary conditions

圖2 速度邊界條件Fig.2 Velocity profile used as a loading function

表1 算例1采用材料參數(shù)Tab.1 Material parameters used in the first example

為顯示熱傳導效應在該模擬過程中的重要性及其對網(wǎng)格依賴性的緩解作用，本文分別進行了絕熱假設模擬(κ=0)和熱擴散模擬，將四分之一板離散成50×50，80×80和100×100的均勻網(wǎng)格，使用隱式時間算法。計算總時長2.0 μs，時間步長取5 ns。絕熱假設和熱擴散模擬的數(shù)值分析結果，如圖3所示。其中包括應力隨時間變化曲線，t=2 μs時方形板主對角線方向(圖1虛線處)的等效塑性應變(EQPS) 分布圖。

從圖1可知，考慮熱傳導效應模擬剪切帶，應力時程曲線及剪切帶寬度、等效塑性應變峰值的結果都比較穩(wěn)定，對網(wǎng)格依賴性小。反之，絕熱假設導致應力時程曲線的應力軟化階段出現(xiàn)明顯的網(wǎng)格依賴性。剪切帶內的塑性變形被放大，過高地估計了剪切帶內的塑性應變及溫度。

圖4表示出模擬終止時100×100網(wǎng)格模型熱擴散項與塑性功產生的熱源項沿對角線的分布圖。圖5為相同工況下的等效塑性應變云圖。從圖5可知，剪切局部化范圍內熱擴散項與熱能同一量級，符號相反。與文獻[13-14]結論一致，在能量方程中引入熱傳導使得熱能與熱擴散互相對抗，形成模型內部尺寸效應，從而減弱網(wǎng)格依賴性。因此，不管從物理角度還是從數(shù)學角度，在能量守恒中考慮熱耗散項均具有重要意義。

4.2 算例2

本算例采用“4.1”節(jié)的模型，材料參數(shù)及網(wǎng)格尺寸均與算例1相同，但此算例進行半顯式離散模擬及隱式算法模擬的計算結果、計算成本等比較。基于算例1，在隱式和顯式模擬過程中均考慮熱擴散效應。

(a) 應力時程曲線

(b) 方形板對角線方向的EQPS分布曲線(t =2 μs)圖3 絕熱假設與熱耗散模擬結果比較Fig.3 Comparison of adiabatic and diffusion assumptions with different meshes

圖4 對角線處熱耗散項與熱源項分布曲線(100×100網(wǎng)格，t=2 μs)Fig.4 Source and diffusive term along the diagonal of the plate plotted at the final time

采用50×50，80×80和100×100的均勻網(wǎng)格， 對應的網(wǎng)格尺寸分別為h=1 μm、h=0.625 μm和h=0.5 μm。依據(jù)式(28)和式(29)及不同網(wǎng)格尺寸，顯式模擬的極限時間步長分別為Δtcr=3.48 ns，Δtcr=1.36 ns和Δtcr=0.87 ns?？紤]非線性問題，最終顯式時間步長分別取為Δt=0.04 ns，Δt=0.02 ns和Δt=0.02 ns。隱式模擬的時間步長均取為Δt=10 ns。計算總時長為2 μs。施加的速度邊界條件與算例1相同，即vr=5 m/s，tr=1 μs。為了區(qū)分沖擊荷載對計算結果的影響，對于100×100的網(wǎng)格還分別進行了vr=10 m/s的兩種不同時間離散的數(shù)值模擬。

圖5 EQPS云圖(100×100網(wǎng)格，t=2 μs)Fig.5 EQPS field at the final time

圖6為兩種離散方法模擬的材料應力隨時間變化的曲線。從圖中可以看出兩種方法的應力曲線吻合很好，且網(wǎng)格的尺寸敏感性小。圖7針對100×100網(wǎng)格，比較了不同時間迭代方法及不同沖擊速度條件下得到的等效塑性應變及溫度場沿主對角線的分布情況。顯式與隱式模擬均在Dell optiplex 755電腦上進行，4 GB內存，4核2.83 GHz。結果表明對于45°剪切帶問題，增大沖擊速度明顯增大了剪切帶內的等效塑性應變及帶內溫度，使剪切帶內的塑性變形更劇烈，但是不同速度邊界條件下，針對所選的時間步長，兩種算法的計算精度都基本一致?？紤]顯式算法的條件穩(wěn)定性，時間步長需要比隱式算法小很多。不同網(wǎng)格的總模擬cpu時長及隱式算法的單位時間步的平均迭代次數(shù)Niter列于表2中(vr=5 m/s)。雖然顯式算法的計算精度與隱式算法相同，但顯式算法的cpu時長遠遠大于隱式算法的時長(見表2)。

圖6 絕熱假設與熱耗散模擬結果比較Fig.6 Comparison of adiabatic and diffusion assumptions with different meshes

(a) 對角線方向EQPS分布曲線

(b) 對角線方向溫度分布曲線圖7 不同速度邊界條件下100×100網(wǎng)格顯式算法與隱式算法結果比較

Fig.7 Results obtained with explicit and implicit scheme plotted at the final time with the finest mesh under different velocity impact

表2 不同算法的Cpu時長比較(vr=5 m/s)Tab.2 Cpu time comparison (vr=5 m/s)

因此，盡管顯式算法經常被廣泛用于求解高速動力學問題，但對于沖擊荷載下45°剪切帶問題，為了滿足計算的精度要求，需要采用很小的時間步長，導致計算所需成本比隱式算法大很多。本文編寫的單元采用隱式算法迭代收斂穩(wěn)定，計算精度高，且因為考慮了熱傳導作用，網(wǎng)格敏感性小。

5 結 論

本文研究了在高速沖擊荷載作用下金屬材料局部剪切帶形成過程的模擬方法，并使用Schur補方法節(jié)約計算成本。在混合有限元的基礎上，基于有限元分析程序FEAP (Finite Element Analysis Program)，使用Fortran語言編寫一個新單元，該單元將位移、溫度、應力和等效塑性應變場統(tǒng)一作為未知量，同步進行求解適用于高應變率作用下塑性變形局部化問題的偏微分方程組，Jacobian矩陣通過對殘余向量在未知量增量方向進行Gateaux求導計算出理論解，避免數(shù)值微分的過程及誤差。針對該計算模型，

(1) 模擬過程中在能量方程中同時考慮了熱傳導效應及塑性功產生的熱能作用，在非線性問題的模型內部引入尺寸效應，有效地限制剪切局部化變形范圍，緩解了剪切帶模擬的網(wǎng)格敏感性現(xiàn)象。

(2) 該模型在空間上Galerkin離散，在時間上分別進行了半顯式算法(應力及塑性應變變量采用隱式時間離散，而位移和溫度變量為顯式迭代)與隱式算法的剪切帶模擬。

(3) 依據(jù)顯式算法與隱式算法模擬結果對比得出，對于剪切帶問題的模擬，顯式算法的精度并不及隱式算法，為了使顯式算法的計算精度滿足要求，需要選取很小的時間步長其計算成本比隱式算法大很多，而文中的模型采用隱式算法迭代，迭代收斂穩(wěn)定，不僅時間成本低，而且計算精度高。

致謝本論文得到國家留學基金委資助(201506710022)，感謝美國哥倫比亞大學土木工程與工程力學系Haim Waisman教授及其課題組的指導與幫助。

參 考 文 獻

[1] 徐永波，白以龍.動態(tài)載荷下剪切變形局部化、微結構演化與剪切斷裂研究進展.力學進展，2007，37 (4)：496-516.

XU Yongbo, BAI Yilong. Shear localization, microstructure evolution and fracture under high-strain rate[J]. Advances in Mechanics,2007,37(4): 496-516.

[2] WRIGHT T W, BATRA R C. The initiation and growth of adiabatic shear bands[J]. International Journal of Plasticity, 1985,1(3): 205-212.

[3] DODD B, BAI Y. Adiabatic shear localization[M]. 2nd ed. Oxford: Elsevier, 2012.

[4] OSOVSKI S, RITTEL D, LANDAU P, et al. Microstructural effects on adiabatic shear band formation[J]. Scripta Materialia, 2012, 66(1): 9-12.

[5] OSOVSKI S, RITTEL D, VENKERT A. The respective influence of microstructural and thermal softening on adiabatic shear localization[J]. Mechanics of Materials, 2013, 56:11-22.

[6] WRIGHT T W. The physics and mathematics of adiabatic shear bands[M]. Cambridge: Cambridge University Press, 2002.

[7] MARCHAND A, DUFFY J. An experimental study of the formation process of adiabatic shear bands in a structural steel[J]. Journal of the Mechanics and Physics of Solids,1988, 36(3): 251-283.

[8] JOHNSON G R, COOK W H. A constitutive model and data for metals subjected to large strains, high strain rates and high temperatures[C]∥Proceedings of the 7th International Symposium on Ballistics. Hague, Netherlands,1983.

[9] LI S, LIU W K, ROSAKIS A J, et al. Mesh free galerkin simulations of dynamic shear band propagation and failure mode transition[J]. International Journal of Solids and Structures, 2002, 39(5): 1213-1240.

[10] BELYTSCHKO T, CHIANG H Y, PLASKACZ E. High resolution two dimensional shear band computations: imperfections and mesh dependence[J]. Computer Methods in Applied Mechanics and Engineering, 1994, 119(1): 1-15.

[11] 王猛,楊明川,羅榮梅,等. 鎢合金桿式彈穿甲侵徹開坑階段絕熱剪切失效的數(shù)值模擬[J].振動與沖擊, 2016, 35(18): 111-116.

WANG Meng, YANG Mingchuan, LUO Rongmei, et al. Numerical simulation on the adiabatic shear failure at cratering stage for tungsten alloy long rod penetrator piercing into armor[J]. Journal of Vibration and Shock, 2016, 35(18): 111-116.

[12] NEEDLEMAN A. Material rate dependence and mesh sensitivity in localization problems[J]. Computer Methods in Applied Mechanics and Engineering, 1988, 67(1): 69-85.

[13] 白以龍，鄭哲敏,俞善炳. 關于熱塑剪切帶的演變[J].力學學報,1986, 18(增刊2): 377-383.

BAI Yilong, ZHENG Zhemin, YU Shanbing. On evolution of thermo-plastic shear band[J]. Chinese Journal of Theoretical and Applied Mechani, 1986, 18(Sup2): 377-383.

[14] 白以龍，鄭哲敏，俞善炳. 熱塑剪切帶后期演變的準定常近似[J].固體力學學報, 1987(2):153-158.

BAI Yilong, ZHENG Zhemin， YU Shanbing. Quasi steady approximation to the late stage evolution of thermo-plastic shear band[J]. Chinese Journal of Solid Mechanics, 1987(2): 153-158.

[15] TAYLOR R L. FEAP—A finite element analysis program[D]. Berkeley, CA:University of California, 2011.

[16] ZHOU M, RAVICHANDRAN G, ROSAKIS A J. Dynamically propagating shear bands in impact-loaded pre-notched plates—II. numerical simulations[J]. Journal of the Mechanics and Physics of Solids, 1996, 44(6): 1007-1032.

[17] MCAULIFFE C, WAISMAN H. A pian-sumihara type element for modeling shear bands at finite deformation[J]. Computational Mechanics, 2013, 53(5): 925-940.

[18] SCHUR J. über potenzreihen, die im innern des einheitskreises beschra?nkt sind[J]. Journal für die Reine und Angewandte Mathematik, 1917, 147:205-232.

[19] BERGER-VERGIAT L, MCAULIFFE C, WAISMAN L. Isogeometric analysis of shear bands[J]. Computational Mechanics, 2014, 54(2): 503-521.

[20] HUGHES T J R. The finite element method: linear static and dynamic finite element analysis[M]. 2nd ed. Mineola, NY: Dover Publications, 2000.