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

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

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

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

白以龍等[13-14]提出了熱塑剪切帶問題，強(qiáng)調(diào)耦合熱傳導(dǎo)和塑性功效應(yīng)的重要性，并基于假設(shè)計(jì)算出一維熱塑剪切帶問題的準(zhǔn)定常近似解，將剪切帶寬度的形成看作是塑性功和熱擴(kuò)散相互平衡的產(chǎn)物。

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

1 控制方程

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

不考慮體力作用，動(dòng)量守恒方程表示為

(1)

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

(2)

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

基于小變形假設(shè)

(3)

(4)

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

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

D=Del+DT+Dp

(5)

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

(6)

剪切帶形成過程中伴隨著應(yīng)變硬化、應(yīng)變率硬化及熱軟化等作用，假設(shè)屈服應(yīng)力與應(yīng)力、溫度及塑性變形相關(guān)，本文采用修正Litonski[16-17]塑性本構(gòu)方程

(7)

(8)

(9)

上式表明塑性應(yīng)變張量與偏應(yīng)力張量方向相同。

由于

(10)

(11)

式中：“:”指張量的雙點(diǎn)積運(yùn)算。

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

(12)

邊界條件為

(13)

2 離散模型

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

(14)

(15)

(16)

3 時(shí)間離散模型

3.1 隱式算法

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

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

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

(17)

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

(18)

式中：

(19)

(20)

(21)

(22)

因此式(18)可簡化為

(23)

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

(24)

用數(shù)學(xué)方法展開

(25)

3.2 半顯式算法

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

(26)

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

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

(27)

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

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

Δtth=min(ρcph2/2κ)

(28)

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

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

(29)

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

4 數(shù)值算例

4.1 算例1

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

缺陷內(nèi)屈服應(yīng)力和屈服應(yīng)變值由二維β函數(shù)折減

(30)

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

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

圖1 方形板結(jié)構(gòu)及邊界條件(虛線區(qū)域?yàn)槿毕菸恢?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

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

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

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

4.2 算例2

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

(a) 應(yīng)力時(shí)程曲線

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

圖4 對角線處熱耗散項(xiàng)與熱源項(xiàng)分布曲線(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)格， 對應(yīng)的網(wǎng)格尺寸分別為h=1 μm、h=0.625 μm和h=0.5 μm。依據(jù)式(28)和式(29)及不同網(wǎng)格尺寸，顯式模擬的極限時(shí)間步長分別為Δtcr=3.48 ns，Δtcr=1.36 ns和Δtcr=0.87 ns。考慮非線性問題，最終顯式時(shí)間步長分別取為Δt=0.04 ns，Δt=0.02 ns和Δt=0.02 ns。隱式模擬的時(shí)間步長均取為Δt=10 ns。計(jì)算總時(shí)長為2 μs。施加的速度邊界條件與算例1相同，即vr=5 m/s，tr=1 μs。為了區(qū)分沖擊荷載對計(jì)算結(jié)果的影響，對于100×100的網(wǎng)格還分別進(jìn)行了vr=10 m/s的兩種不同時(shí)間離散的數(shù)值模擬。

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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

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

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

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

參 考 文 獻(xiàn)

[1] 徐永波，白以龍.動(dòng)態(tài)載荷下剪切變形局部化、微結(jié)構(gòu)演化與剪切斷裂研究進(jìn)展.力學(xué)進(jìn)展，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].振動(dòng)與沖擊, 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] 白以龍，鄭哲敏,俞善炳. 關(guān)于熱塑剪切帶的演變[J].力學(xué)學(xué)報(bào),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] 白以龍，鄭哲敏，俞善炳. 熱塑剪切帶后期演變的準(zhǔn)定常近似[J].固體力學(xué)學(xué)報(bào), 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.