一種高強(qiáng)度鋼的動(dòng)態(tài)力學(xué)本構(gòu)參數(shù)識(shí)別和數(shù)值模擬*

謝 恒，呂振華

(清華大學(xué)汽車工程系，北京100084)

目前有許多學(xué)者對(duì)彈道沖擊下的高強(qiáng)度鋼等金屬材料的動(dòng)態(tài)力學(xué)特性進(jìn)行了研究［1-5］，雖然已有一些成果，但這些研究多是針對(duì)壓縮屈服特性的，達(dá)到的應(yīng)變率不夠高且很難得到數(shù)值模擬所需的斷裂參數(shù)。因此對(duì)國產(chǎn)高強(qiáng)度鋼進(jìn)行動(dòng)態(tài)拉伸力學(xué)特性測(cè)試并確定本構(gòu)模型的斷裂參數(shù)具有重要的工程價(jià)值。另外，高強(qiáng)度鋼板難以加工成圓柱形啞鈴試件，因此平板試件的實(shí)驗(yàn)和數(shù)據(jù)處理顯得至關(guān)重要。

平板試件進(jìn)入頸縮階段時(shí)，最小截面等效應(yīng)力分布較均勻，而真實(shí)應(yīng)變的分布很復(fù)雜。對(duì)于平板試件真實(shí)應(yīng)變的計(jì)算，主要有2 種方法。Z.L.Zhang 等［6］采用截面積計(jì)算的平均應(yīng)變作為真實(shí)應(yīng)變(以下簡(jiǎn)稱平均法)，即ε=ln (S0/S) ，但應(yīng)變一般是不均勻的，截面中心應(yīng)變比平均應(yīng)變大得多;I.Scheider等［7］采用數(shù)字圖像相關(guān)算法(DIC)直接得到試件表面應(yīng)變(以下簡(jiǎn)稱表面法)。然而用高速攝影對(duì)試件進(jìn)行應(yīng)變測(cè)量非常困難，且表面應(yīng)變也不一定等同于中心應(yīng)變。本文中，擬結(jié)合有限元數(shù)值模擬提出一套簡(jiǎn)單易行的平板試件真實(shí)應(yīng)變識(shí)別方法。

1 SHTB 實(shí)驗(yàn)

1.1 實(shí)驗(yàn)方案

分離式Hopkinson 桿技術(shù)的理論基礎(chǔ)是一維應(yīng)力波理論。在短試件中的彈性波多次反射獲得均勻應(yīng)力。通過測(cè)定桿上的應(yīng)變來推導(dǎo)試樣的應(yīng)力應(yīng)變關(guān)系。旋轉(zhuǎn)盤式SHTB 裝置如圖1 所示。試件外形及其連接方式如圖2 所示。

利用一維應(yīng)力假定和均勻性假定，可以獲得試件的應(yīng)變、應(yīng)變率、應(yīng)力隨時(shí)間的變化關(guān)系，進(jìn)而得到試件材料在各個(gè)應(yīng)變率下的應(yīng)力應(yīng)變關(guān)系［8］

圖1 SHTB 實(shí)驗(yàn)裝置簡(jiǎn)圖Fig.1 Scheme of SHTB

1.2 實(shí)驗(yàn)結(jié)果

制備了N-1 鋼的準(zhǔn)靜態(tài)、動(dòng)態(tài)2 種平板試件，利用旋轉(zhuǎn)盤式SHTB 裝置進(jìn)行了高應(yīng)變率沖擊拉伸實(shí)驗(yàn)，得到了N-1 鋼在各種應(yīng)變率下的工程應(yīng)力應(yīng)變關(guān)系，如圖3 所示。表1 給出了屈服應(yīng)力σ0.2、σ0.5隨應(yīng)變率的變化關(guān)系。平均應(yīng)變率為約1 000 s-1時(shí)，N-1 鋼的屈服應(yīng)力比準(zhǔn)靜態(tài)值增大約10%。在高應(yīng)變率下，只有在最大載荷σp，max之前的工程應(yīng)力應(yīng)變關(guān)系可以近似等價(jià)于真實(shí)應(yīng)力真實(shí)應(yīng)變關(guān)系。之后由于頸縮和絕熱溫升，應(yīng)力應(yīng)變關(guān)系耦合了結(jié)構(gòu)響應(yīng)和熱軟化效應(yīng)，需要經(jīng)過處理才能得到真實(shí)應(yīng)力應(yīng)變關(guān)系曲線。

圖2 動(dòng)態(tài)試件及其連接方式Fig.2 The dynamic specimen and its connection mode

表1 N-1 鋼的屈服應(yīng)力與應(yīng)變率關(guān)系Table 1 Relation between yield stress and strain rate of N-1 steel

圖3 各種應(yīng)變率下N-1 鋼的工程應(yīng)力-應(yīng)變關(guān)系Fig.3 Engineering stress-strain relation of N-1 steel at different strain rates

2 本構(gòu)模型參數(shù)識(shí)別

采用修正的Johnson-Cook(MJC)模型［9-10］，其等效應(yīng)力

式中:εeq為等效塑性應(yīng)變，A、B、n、C、m 為材料參數(shù)。量綱一塑性應(yīng)變率由確定為參考應(yīng)變率。歸一化的溫度由T*=(T－Tr)/(Tm－Tr)確定，T 為絕對(duì)溫度，Tr、Tm分別為室溫和材料的熔點(diǎn)。由絕熱引起的溫升

式中:ρ 為材料密度，cp為比定壓熱容，χ 為Taylor-Quinney 系數(shù)即塑性功轉(zhuǎn)熱比例。

MJC 模型的斷裂準(zhǔn)則部分由Cockcroft-Latham(CL)準(zhǔn)則［11］定義，模型表示為

即只有當(dāng)拉應(yīng)力產(chǎn)生時(shí)才可能斷裂。Wcr為單位體積的臨界斷裂塑性功。在有限元分析中，當(dāng)單元的塑性功W 達(dá)到臨界值Wcr時(shí)，刪除該單元。

由表1 中σ0.2、σ0.5隨應(yīng)變率的變化關(guān)系可分別擬合得出應(yīng)變率強(qiáng)化參數(shù)C 的2 個(gè)值，取其平均值C=0.009。根據(jù)文獻(xiàn)［1］，鋼的溫度軟化系數(shù)m 可近似取1.0。選擇N-1 鋼的應(yīng)變率為1 358 s-1的工程應(yīng)力應(yīng)變實(shí)驗(yàn)曲線進(jìn)行參數(shù)A、B、n 和Wcr的識(shí)別。

2.1 真實(shí)應(yīng)力的識(shí)別

由于頸縮后已經(jīng)不是單向拉伸應(yīng)力的狀態(tài)，等效應(yīng)力不再等于軸向應(yīng)力。無論將試件的原始截面積還是瞬時(shí)截面積直接代入應(yīng)力=載荷÷面積的計(jì)算公式，都是不準(zhǔn)確的。因此在識(shí)別參數(shù)A、B、n 之前必須對(duì)這個(gè)階段的真實(shí)應(yīng)力、真實(shí)應(yīng)變進(jìn)行識(shí)別。試件頸縮區(qū)外形如圖4 所示。

P.W.Bridgman［12］提出的關(guān)于平板單向拉伸時(shí)頸縮區(qū)的應(yīng)力修正公式為

式中:σeq為等效應(yīng)力為平均真實(shí)軸向應(yīng)力，σeng為工程應(yīng)力，S0、S 分別為初始截面積和當(dāng)前截面積，S0=3.85 mm2。c 為頸縮區(qū)應(yīng)力修正系數(shù)，幾何參數(shù)a、R 分別為最小截面處的厚度和表面曲率半徑。a/R 由經(jīng)驗(yàn)公式［13］計(jì)算

式中:ε=ln (S0/S) 為當(dāng)前平均對(duì)數(shù)應(yīng)變，εp，max為工程應(yīng)力應(yīng)變曲線最大載荷時(shí)的平均對(duì)數(shù)應(yīng)變。式(8)對(duì)各種鋼材均有效。

從實(shí)驗(yàn)結(jié)果可以看出，動(dòng)態(tài)試件斷裂時(shí)的截面形狀為枕形，見圖5，最終截面面積S=2.70 mm2。

圖4 試件頸縮變形示意圖Fig.4 Deformation of the necking in the specimen

圖5 斷裂截面形狀與尺寸Fig.5 Shape and dimensions of the fractured cross-section

2.2 真實(shí)應(yīng)變的識(shí)別

通過對(duì)中心最大真實(shí)應(yīng)變與最小截面的厚度縮減率關(guān)系的研究，得出經(jīng)驗(yàn)公式用以計(jì)算真實(shí)應(yīng)變，方法如下。設(shè)最小截面的厚度縮減率為

式中:a0、a 分別為試件的最小截面的初始厚度、當(dāng)前厚度。εeq可寫為rt與應(yīng)變修正因子fe的乘積

為確定fe，取A=1 400 MPa，B=500，1 000，2 000 MPa，n=0.1，0.5，1.0;對(duì)這幾種硬化參數(shù)相差很大的材料試件進(jìn)行動(dòng)態(tài)單向拉伸的數(shù)值模擬。

數(shù)值模擬的試件外形、拉伸速度等具體設(shè)置將在第3 節(jié)中介紹。此處模型省略了入射桿、透射桿，直接對(duì)試件兩端施加位移。數(shù)值模擬時(shí)確保各材料試件的最小截面的厚度縮減率稍大于實(shí)驗(yàn)值0.27。

數(shù)值模擬實(shí)驗(yàn)結(jié)果如圖6 所示。盡管材料硬化參數(shù)不同，fe與rt存在幾乎相同的非線性函數(shù)關(guān)系。選取擬合函數(shù)和擬合結(jié)果分別為

圖6 不同硬化參數(shù)材料的應(yīng)變修正因子Fig.6 Strain correction factors for different material laws

以B=1 000 MPa、n=0.5 的材料為例，將本文中提出的方法與文獻(xiàn)中常用的2 種方法進(jìn)行比較。以數(shù)值模擬結(jié)果為精確值，當(dāng)變形到rt=0.272 7 時(shí)，此時(shí)中心真實(shí)應(yīng)變已經(jīng)達(dá)到0.527 1。利用本文中提出的修正因子法即式(9)～(12)，求得真實(shí)應(yīng)變?yōu)?.520 4;借鑒表面法［7］，取有限元模型試件表面的應(yīng)變求得真實(shí)應(yīng)變?yōu)?.474 4;而基于平均法［6］求得真實(shí)應(yīng)變僅為0.35。

用面積變化來表征式(10)也是可以的。但在本例中，鋼的硬化特性函數(shù)不同時(shí)，fe與面積縮減率(S0－S)/S0關(guān)系的一致性稍差。

最小截面的厚度縮減率比面積縮減率準(zhǔn)確的原因與試件外形尺寸有關(guān)。本節(jié)待識(shí)別的是最大應(yīng)變，即最小截面的中心應(yīng)變。對(duì)于本例中的平板試件，中心應(yīng)變與最小截面的厚度縮減率直接相關(guān)。而對(duì)于圓柱試件，由于試件在徑向各個(gè)方向均勻收縮，因此中心應(yīng)變與最小截面的面積縮減率直接相關(guān)。

2.3 硬化參數(shù)A、B、n 的識(shí)別

根據(jù)斷裂時(shí)的工程應(yīng)力和式(6)～(8)，可以確定斷裂時(shí)的真實(shí)應(yīng)力σeq，f。根據(jù)斷裂時(shí)的截面尺寸和式(9)～(12)，可以確定斷裂時(shí)的真實(shí)應(yīng)變?chǔ)舉q，f。由工程應(yīng)力應(yīng)變關(guān)系曲線的最大載荷之前的部分和σeq，f、εeq，f可擬合出參數(shù)A、B、n，結(jié)果如圖7 所示。

2.4 CL 斷裂準(zhǔn)則參數(shù)Wcr的確定

由式(4)可知，為確定Wcr，需知σ1和εeq的關(guān)系。而σ1在頸縮區(qū)分布不均勻，難以解析計(jì)算得到。利用以上步驟確定的本構(gòu)參數(shù)進(jìn)行試件動(dòng)態(tài)拉伸實(shí)驗(yàn)的數(shù)值模擬。先取Wcr為一個(gè)大數(shù)，即假設(shè)不發(fā)生斷裂，模型的其他設(shè)置見第3 節(jié)。模擬得到斷裂面中心的W-εeq關(guān)系，見圖8，再由εeq，f確定Wcr。

圖7 修正后擬合的真實(shí)應(yīng)力應(yīng)變曲線Fig.7 Fitted true stress-strain curves after correction

圖8 斷裂面中心的W-εeq關(guān)系Fig.8 W-εeq relation at the center of the fractured cross-section

2.5 MJC 本構(gòu)模型參數(shù)識(shí)別結(jié)果

由以上步驟擬合出的MJC 模型參數(shù)分別為:A=1 400 MPa，B=712 MPa，n=0.285 8，ε·0=10-3s-1，C=0.009，Wcr=1 200 MPa。數(shù)值模擬中涉及的彈性參數(shù)、溫度效應(yīng)參數(shù)取自文獻(xiàn)［11］。

3 數(shù)值模擬驗(yàn)證

3.1 有限元模型的建立

使用非線性有限元分析軟件LS-DYNA 進(jìn)行SHTB 實(shí)驗(yàn)的數(shù)值模擬，以驗(yàn)證采用的實(shí)驗(yàn)方法、提出的本構(gòu)模型參數(shù)識(shí)別方法以及數(shù)值模擬方法的正確性。

對(duì)平板試件及桿-桿型沖擊拉伸實(shí)驗(yàn)系統(tǒng)做1/4 對(duì)稱簡(jiǎn)化建模。輸入桿長360 mm，應(yīng)變片距試件端300 mm;輸出桿長750 mm，應(yīng)變片距試件端150 mm。試件與輸入桿和輸出桿固結(jié)。輸入桿一端結(jié)點(diǎn)施加速度脈沖載荷。速度脈沖曲線由實(shí)驗(yàn)時(shí)試件中的應(yīng)變率-時(shí)間曲線估算，即v=ε·l0，試件單元尺度約0.18 mm，單元總數(shù)約500 000。有限元模型的試件段網(wǎng)格見圖9。

圖9 有限元模型試件段Fig.9 The specimen in the finite element model

3.2 數(shù)值模擬結(jié)果與討論

由數(shù)值模擬的應(yīng)變片輸出的波形如圖10 所示。由應(yīng)變波形和式(1)計(jì)算得到的工程應(yīng)力應(yīng)變關(guān)系曲線如圖11 所示。初始屈服時(shí)數(shù)值模擬的應(yīng)力與實(shí)驗(yàn)結(jié)果存在一定偏差的原因有:(1)MJC 模型的函數(shù)形式限制了擬合精度;(2)實(shí)驗(yàn)中機(jī)械系統(tǒng)多處存在濾波作用。數(shù)值模擬的斷裂應(yīng)變值與實(shí)驗(yàn)值的誤差為約10%，可能是由于N-1 鋼沖擊斷裂試件的局部頸縮不明顯，本文中利用了彌散頸縮變形信息。但這個(gè)結(jié)果遠(yuǎn)優(yōu)于基于平均法的結(jié)果，也優(yōu)于基于表面法的結(jié)果，對(duì)于工程應(yīng)用已經(jīng)足夠。

數(shù)值模擬的斷裂試件側(cè)面和截面外形見圖12 ～13?？梢?，數(shù)值模擬結(jié)果與實(shí)驗(yàn)結(jié)果較吻合。

圖11 數(shù)值模擬的工程應(yīng)力應(yīng)變曲線與實(shí)驗(yàn)結(jié)果的比較Fig.11 Engineering stress-strain curves by different methods

圖12 斷裂試件的外形Fig.12 Shapes of the fractured specimen

圖13 斷裂截面的外形Fig.13 Shapes of the fractured cross-section

4 結(jié) 論

平均應(yīng)變率為約1 000 s-1時(shí)，N-1 鋼的屈服應(yīng)力比準(zhǔn)靜態(tài)值增大約10%。

通過數(shù)值模擬研究了真實(shí)應(yīng)變與最小截面的厚度縮減率的關(guān)系，由應(yīng)變修正因子可以確定真實(shí)應(yīng)變。對(duì)于本文的試件，該方法比已有文獻(xiàn)中常用的平均法更準(zhǔn)確，比表面法更方便。根據(jù)所得到的斷裂面中心的真實(shí)應(yīng)力、真實(shí)應(yīng)變，結(jié)合頸縮前的工程應(yīng)力應(yīng)變關(guān)系，識(shí)別出了MJC 模型參數(shù)。

高應(yīng)變率下，試件外形、工程應(yīng)力應(yīng)變關(guān)系的實(shí)驗(yàn)結(jié)果與模擬結(jié)果均較一致，表明所識(shí)別的模型參數(shù)可較好地描述材料在高速變形時(shí)的力學(xué)行為，所用方法也可用于其他類似材料本構(gòu)模型參數(shù)的確定。

對(duì)高強(qiáng)度鋼的動(dòng)態(tài)本構(gòu)模型參數(shù)擬合仍是初步的，需要對(duì)其他外形尺寸的試件開展進(jìn)一步的研究工作。如果實(shí)驗(yàn)中可以連續(xù)測(cè)量最小截面的厚度縮減率并利用局部頸縮的信息，擬合結(jié)果應(yīng)會(huì)更準(zhǔn)確。

［1］ B?rvik T，Dey S，Clausen A H.Perforation resistance of five different high-strength steel plates subjected to small-arms projectiles［J］.International Journal of Impact Engineering，2009，36(7):948-964.

［2］ 夏源明，周元鑫，楊報(bào)昌.馬氏體含量對(duì)雙相鋼力學(xué)性能的影響［J］.材料科學(xué)與工藝，1997，5(1):60-64.XIA Yuan-ming，ZHOU Yuan-xin，YANG Bao-chang.The influence of martensite fraction on the mechanical properties of dual phase steels［J］.Material Science and Technology，1997，5(1):60-64.

［3］ 胡昌明，賀紅亮，胡時(shí)勝.45 號(hào)鋼的動(dòng)態(tài)力學(xué)性能研究［J］.爆炸與沖擊，2003，23(2):188-192.HU Chang-ming，HE Hong-liang，HU Shi-sheng.A study on dynamic mechanical behaviors of 45 steel［J］.Explosion and Shock Waves，2003，23(2):188-192.

［4］ 胡昌明，賀紅亮，胡時(shí)勝.D6AC 鋼的動(dòng)態(tài)力學(xué)性能研究［J］.兵器材料科學(xué)與工程，2002，26(2):26-29.HU Chang-ming，HE Hong-liang，HU Shi-sheng.Research of dynamic mechanical behavior of D6AC steel［J］.Ordnance Material Science and Engineering，2002，26(2):26-29.

［5］ 田杰，胡時(shí)勝.G50 鋼動(dòng)態(tài)力學(xué)性能的實(shí)驗(yàn)研究［J］.工程力學(xué)，2006，23(6):107-109.TIAN Jie，HU Shi-sheng.Research of dynamic mechanical behaviors of G50 steel［J］.Engineering Mechanics，2006，23(6):107-109.

［6］ Zhang Z L，Hauge M，?deg?rd J，et al.Determining material true stress-strain curve from tensile specimens with rectangular cross-section［J］.International Journal of Solids and Structures，1999，36(23):3497-3516.

［7］ Scheider I，Brocks W，Cornec A.Procedure for the determination of true stress-strain curves from tensile tests with rectangular cross-section specimens［J］.Journal of Engineering Materials and Technology，2004，126:70-76.

［8］ 王禮立.應(yīng)力波基礎(chǔ)［M］.北京:國防工業(yè)出版社，2005:52-60.

［9］ B?rvik T，Langseth M，Hopperstad O S，et al.Ballistic penetration of steel plates［J］.International Journal of Impact Engineering，1999，22(9/10):855-886.

［10］ B?rvik T，Hopperstad O S，Berstad T，et al.Numerical simulation of plugging failure in ballistic penetration［J］.International Journal of Solids and Structures，2001，38(34/35):6241-6264.

［11］ Dey S，B?rvik T，Hopperstad O S，et al.On the influence of fracture criterion in projectile impact of steel plates［J］.Computational Materials Science，2006，38(1):176-191.

［12］ Bridgman P W.Studies in Large Plastic Flow and Fracture［M］.New York:McGraw-Hill，1952:32-37.

［13］ Le Roy G，Embury J D，Edwards G，et al.A model of ductile fracture based on the nucleation and growth of voids［J］.Acta Metallurgica，1981，29(8):1509-1522.