一種高強度鋼的動態力學本構參數識別和數值模擬*

謝 恒，呂振華

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

目前有許多學者對彈道沖擊下的高強度鋼等金屬材料的動態力學特性進行了研究［1-5］，雖然已有一些成果，但這些研究多是針對壓縮屈服特性的，達到的應變率不夠高且很難得到數值模擬所需的斷裂參數。因此對國產高強度鋼進行動態拉伸力學特性測試并確定本構模型的斷裂參數具有重要的工程價值。另外，高強度鋼板難以加工成圓柱形啞鈴試件，因此平板試件的實驗和數據處理顯得至關重要。

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

1 SHTB 實驗

1.1 實驗方案

分離式Hopkinson 桿技術的理論基礎是一維應力波理論。在短試件中的彈性波多次反射獲得均勻應力。通過測定桿上的應變來推導試樣的應力應變關系。旋轉盤式SHTB 裝置如圖1 所示。試件外形及其連接方式如圖2 所示。

利用一維應力假定和均勻性假定，可以獲得試件的應變、應變率、應力隨時間的變化關系，進而得到試件材料在各個應變率下的應力應變關系［8］

圖1 SHTB 實驗裝置簡圖Fig.1 Scheme of SHTB

1.2 實驗結果

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

圖2 動態試件及其連接方式Fig.2 The dynamic specimen and its connection mode

表1 N-1 鋼的屈服應力與應變率關系Table 1 Relation between yield stress and strain rate of N-1 steel

圖3 各種應變率下N-1 鋼的工程應力-應變關系Fig.3 Engineering stress-strain relation of N-1 steel at different strain rates

2 本構模型參數識別

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

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

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

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

即只有當拉應力產生時才可能斷裂。Wcr為單位體積的臨界斷裂塑性功。在有限元分析中，當單元的塑性功W 達到臨界值Wcr時，刪除該單元。

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

2.1 真實應力的識別

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

P.W.Bridgman［12］提出的關于平板單向拉伸時頸縮區的應力修正公式為

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

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

從實驗結果可以看出，動態試件斷裂時的截面形狀為枕形，見圖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 真實應變的識別

通過對中心最大真實應變與最小截面的厚度縮減率關系的研究，得出經驗公式用以計算真實應變，方法如下。設最小截面的厚度縮減率為

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

為確定fe，取A=1 400 MPa，B=500，1 000，2 000 MPa，n=0.1，0.5，1.0;對這幾種硬化參數相差很大的材料試件進行動態單向拉伸的數值模擬。

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

數值模擬實驗結果如圖6 所示。盡管材料硬化參數不同，fe與rt存在幾乎相同的非線性函數關系。選取擬合函數和擬合結果分別為

圖6 不同硬化參數材料的應變修正因子Fig.6 Strain correction factors for different material laws

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

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

最小截面的厚度縮減率比面積縮減率準確的原因與試件外形尺寸有關。本節待識別的是最大應變，即最小截面的中心應變。對于本例中的平板試件，中心應變與最小截面的厚度縮減率直接相關。而對于圓柱試件，由于試件在徑向各個方向均勻收縮，因此中心應變與最小截面的面積縮減率直接相關。

2.3 硬化參數A、B、n 的識別

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

2.4 CL 斷裂準則參數Wcr的確定

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

圖7 修正后擬合的真實應力應變曲線Fig.7 Fitted true stress-strain curves after correction

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

2.5 MJC 本構模型參數識別結果

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

3 數值模擬驗證

3.1 有限元模型的建立

使用非線性有限元分析軟件LS-DYNA 進行SHTB 實驗的數值模擬，以驗證采用的實驗方法、提出的本構模型參數識別方法以及數值模擬方法的正確性。

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

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

3.2 數值模擬結果與討論

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

數值模擬的斷裂試件側面和截面外形見圖12 ～13。可見，數值模擬結果與實驗結果較吻合。

圖11 數值模擬的工程應力應變曲線與實驗結果的比較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 結 論

平均應變率為約1 000 s-1時，N-1 鋼的屈服應力比準靜態值增大約10%。

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

高應變率下，試件外形、工程應力應變關系的實驗結果與模擬結果均較一致，表明所識別的模型參數可較好地描述材料在高速變形時的力學行為，所用方法也可用于其他類似材料本構模型參數的確定。

對高強度鋼的動態本構模型參數擬合仍是初步的，需要對其他外形尺寸的試件開展進一步的研究工作。如果實驗中可以連續測量最小截面的厚度縮減率并利用局部頸縮的信息，擬合結果應會更準確。

［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］ 夏源明，周元鑫，楊報昌.馬氏體含量對雙相鋼力學性能的影響［J］.材料科學與工藝，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］ 胡昌明，賀紅亮，胡時勝.45 號鋼的動態力學性能研究［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］ 胡昌明，賀紅亮，胡時勝.D6AC 鋼的動態力學性能研究［J］.兵器材料科學與工程，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］ 田杰，胡時勝.G50 鋼動態力學性能的實驗研究［J］.工程力學，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］ 王禮立.應力波基礎［M］.北京:國防工業出版社，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.