基于單向拉伸仿真對各向異性薄板本構修正的研究

張東凱，李曉廣，余亞克

(河鋼集團邯鋼技術中心，河北 邯鄲 056038)

傳統板料沖壓成形主要依靠設計者的經驗，即通過不斷修模、試模、調整間距及壓邊力等試錯方式，來解決成形中出現起皺、破裂等問題。隨著汽車工業的發展，板料沖壓成形工藝趨于復雜，傳統方法導致零件生產周期長、效率低，已經不能適應成形工藝的發展。而計算機技術和有限元方法的發展給板料成形帶來了新方法，通過計算機仿真模擬沖壓工藝、指導沖壓設計、減少試制次數、縮短生產周期，從而降低制造成本，提高生產效率。要提高模擬精度，使仿真模擬更加接近實際成形過程，必須選擇合理的本構模型[1]，如屈服準則、硬化模型、損傷斷裂模型等。軋制板料帶有各向異性行為，Hill48屈服準則可以準確預測鋼材的屈服軌跡，胡成武[2]利用Hill48平面應力方程和Hill塑性應變比建立了預測任意方向屈服方程，預測DC04和DC05板料在0°、45°、90°的Hill48屈服值，并將結果與試驗數據對比，對比結果十分吻合。牛超[3]利用Hill48屈服函數和Swift等向強化本構數值模擬了0°、45°、90°真應力-應變曲線。趙坤民[4]利用ABAQUS/STANDARD進行拉伸仿真分析，得到模擬頸縮后各向異性塑性變形的合適單元。ManSoo Joun等[5]使用剛體-塑性有限元方法分析了拉伸試驗結果，迭代后的應力-應變曲線精確地預測了頸縮點，減小了拉伸試驗和分析結果的之間的載荷差異。王友德[6]將單軸拉伸試驗的載荷-位移曲線作為目標，通過試驗-數值耦合方法(NE法)，不斷迭代修正輸入的應力-應變曲線，直至載荷-位移曲線達到要求，即得到真實應力-應變曲線，并比較了NE法較傳統理論方法的優點。王玉凱[7]利用ABAQUS/EXPLICIT對Q235冷軋鋼進行拉伸斷裂數值模擬，驗證了Johnson-Cook本構模型對鋼材拉伸性能的合理性。陳輝[8]應用TF方法對單軸本構模型進行有限元迭代計算使得模擬計算結果與圓環試樣壓縮試驗結果滿足一致性，從而得到304材料的本構關系參數。姚迪[9]基于TF方法模擬的載荷-位移曲線、漏斗根部直徑-位移曲線和漏斗變形輪廓線等均與試驗結果吻合良好，試樣表面模擬應變也與DIC測試結果吻合。蔡力勛[10]基于單軸拉伸真應力-真應變曲線，采用二分法并結合有限元模擬迭代載荷-位移曲線，得到真實反映材料拉伸全程變形特征的等效真應力-真應變曲線。祁爽[11]基于有限元輔助方法得到了四種工程材料的等效全程單軸本構關系，同時分析了應力三軸度。彭云強[12]基于ANSYS APDL語言和有限元輔助方法獲得了某型鋼的本構關系，并借助ABAQUS中GTN損傷模型預測了裂紋擴展。以上仿真采用的是各向同性屈服準則，未分析各向異性屈服準則的影響。

本文對DC04進行了單向拉伸試驗，得到載荷位移曲線和真應力-塑性應變曲線。基于有限元輔助方法，在ABAQUS/EXPLICIT中合理簡化拉伸模型，耦合邊界條件，引入Hill48各向異性屈服準則，采用動力顯示算法和適當質量縮放進行準靜態單向拉伸仿真，迭代后的載荷-位移曲線與試驗結果吻合，此時得到DC04板材的合理塑性本構，可作為成形仿真的本構輸入。最后對比分析了Hill48屈服準則和Mises屈服準則下的仿真結果，結果顯示Hill48屈服準則狀態下的載荷位移曲線與試驗結果符合更好。

1 試驗材料及結果

拉伸試驗材料為1 mm的DC04冷軋鋼板，按照GB/T 228.1-2010[13]切取與軋制方向成0°、45°、90°的標準試樣，試樣尺寸見圖1。拉伸試驗在150 kN的電子力學拉伸機上進行，根據GB/T 5027-2007[14],夾頭移動速度定為2 mm/min。拉伸斷裂試件見圖2，由試驗所得到的載荷位移曲線見圖3，各個方向屈服應力σα和塑性應變比rα見表1。

圖1 單向拉伸試樣形狀及尺寸 (單位：mm)Fig.1 Shape and geometric dimension of uniaxial tensile specimen(unit: mm)

圖2 拉伸后的DC04試件(從上到下依次為0°、45°、90°)Fig.2 The specimens of DC04 after stretching(From top to bottom, 0°, 45°, 90°)

圖3 單向拉伸試驗載荷位移曲線Fig.3 Load-displacement curve of unidirectional tensile test

基于體積不可壓縮理論，將工程應力-工程應變按式(1)換算為真應力-應變曲線，去除彈性應變后得到真應力-塑性應變曲線，見圖4，將其作為單向拉伸仿真初始本構輸入。

圖4 真應力-塑性應變曲線 Fig.4 True stress and plastic strain

(1)

式中，σT為真應力，εT為真應變，σE為工程應力，εE為工程應變。

2 有限元數值模擬

2.1 有限元輔助方法

圖5 迭代修正方法Fig.5 Iterative correction method

2.2 幾何及邊界條件

由于考察的是拉伸件的中部區域，由圣維南原理，去除了夾持段和圓角過渡段，留下標距段，尺寸為1 mm×20 mm×80 mm，見圖6。

圖6 單向拉伸幾何模型及邊界條件Fig.6 Unidirectional stretching geometric model and boundary conditions

為了模擬邊界條件，將模型兩邊分別設置兩個參考點RP1(右端)，RP2(左端)，位置為過渡圓角與直段的相交中心，按圖中坐標為(-60,0)，(60,0)，將RP1與右端面耦合，耦合位移轉角，同理應用于RP2。

2.3 材料模型

采用動力顯示算法，需要輸入DC04的密度7.85×10-9t/mm3，彈性模量210 000 MPa，泊松比0.3。

塑性本構首次輸入試驗真應力-塑性應變曲線，這里以真應力和塑性應變表格形式，將擬合曲線輸入到ABAQUS塑性屬性設置中，需要注意的是真應力為遞增形式，在循環迭代本構后，重新輸入塑性表格數據時需去除遞減真應力段。

由ABAQUS手冊，對于板料(平面應力狀態)，Hill48各向異性屈服準則中的應變比R11，R22，R33，R12，R13，R23用各向異性系數r值表示為

R11=R13=R23=1

(2)

(3)

(4)

(5)

按公式2～5，DC04的R11，R22，R33，R12，R13，R23值見表2。由于引入各向異性，需要定義材料主方向，這里定義幾何模型縱向為材料主方向。

表2 Hill48屈服準則的應變比參數

2.4 分析步及網格設置

由于采用準靜態單向拉伸仿真，所以分析步采用動力顯示分析，時間間隔為0.1 s，對幾何模型進行適當的質量縮放，后續分析表明質量縮放因子取1000，滿足準靜態仿真要求。

由于要循環修正真應力和塑性應變，所以在ABAQUS場輸出中需勾選MISES，PEEQ作為修正結果，考慮計算效率，時間均勻間隔200。

為了分析計算結果，得到載荷位移曲線，需要在歷史輸出中勾選參考點RP1的支反力RF1及縱向位移U1，同理，時間均勻間隔200。

整體網格尺寸為1，厚度方向布置了4個種子，將厚度方向分割了3段，中間部分對縱向邊線增加種子個數，實現網格細化，見圖7。在仿真中總結經驗，調節單元類型，其中各向異性采用非協調模式。

圖7 網格模型Fig.7 Grid model

3 仿真結果

3.1 準靜態分析

在ABAQUS/EXPLICIT中，能量平衡方程可以表示為

constant=ETOT=EI+EV+EFD+EKE-EW

(6)

其中，EI為內能(包括彈性應變能、塑性應變能和與沙漏控制相關的偽能)，EV為粘性機制耗散的能量，EFD為摩擦耗散的能量，EKE為動能，EW為外力功，ETOT為系統的總能量。

由于不存在粘彈性材料、離散的減震器及材料阻尼，所以粘性機制耗散的能量EV很小。單向拉伸仿真沒有摩擦，故EFD為零。只剩下EI+EKE-EW，因此對于單軸拉伸仿真，若變形材料的動能EKE占內能EI的比例很小(一般為1%～5%)，那么拉伸試件的外力功EW約等于試件的內能EI，那么仿真就是準靜態的。

從圖8中可以看出，整個模型的動能ALLKE(紅色曲線)對時間的積分只占整個模型內能ALLIE(藍色曲線,基本與X軸重合)對時間積分的很小部分，所以是準靜態的。

圖8 準靜態拉伸能量歷程 Fig.8 Quasi-static tensile energy history

3.2 Hill48屈服準則下的載荷位移曲線分析

單向拉伸仿真得到等效Mises應力、等效塑性應變PEEQ云圖見圖9、10，局部出現頸縮，頸縮區域中間單元的等效Mises應力與等效塑性應變PEEQ即為新的塑性本構。按照圖5的修正方法，對材料本構進行了兩次迭代，每次迭代仿真后的載荷位移曲線及與試驗數據的標準差見圖11。

圖9 等效應力云圖(單位：Mpa) Fig.9 Equivalent stress cloud(unit: MPa)

從圖11中可以看出，將單向拉伸試驗數據以真應力-塑性應變的形式輸入，經過單向拉伸仿真后的載荷位移曲線與試驗數據的標準差為959.01，差異主要在損傷斷裂階段，說明單向拉伸試驗中頸縮部分的本構數據無法準確描述損傷斷裂段，因此需要對本構進行修正。采用迭代修正方法對塑性本構進行第一次修正后的仿真載荷位移曲線與試驗數據標準差為215.82，誤差下降約3.4倍，經過仿真后的載荷位移曲線在彈性和塑性強化階段與試驗數據吻合，損傷斷裂段開始接近試驗曲線，損傷斷裂的本構得到有效修正。經過第二次迭代修正后的仿真載荷位移曲線與試驗數據的標準差為208.89，與試驗數據的差距進一步變小，標準差相對于第一次修正下降明顯變緩，此時的真應力-塑性應變曲線(見表3)即為目標本構，可應用于汽車零部件的沖壓成形仿真。

圖10 等效塑性應變云圖 Fig.10 Equivalent plastic strain cloud

圖11 迭代仿真后的載荷位移曲線 Fig.11 Load-displacement curve after iterative simulation based

表3 迭代后的真應力-塑性應變本構

3.3 兩種屈服準則下的載荷位移曲線對比

從圖12中可以看出，Hill48屈服準則狀態下的載荷位移曲線和Mises屈服準則狀態下的結果在彈性階段和塑性強化階段一致，在損傷斷裂段出現差異,且與試驗數據的標準差分別為959.01和527.87，說明屈服準則對仿真損傷斷裂段影響較大。

圖12 不同屈服準則下的仿真結果Fig.12 Simulation results under different yield criteria

從圖13中可以看出，塑性本構經過第一次迭代后，Hill48屈服準則狀態下的載荷位移曲線和Mises屈服準則狀態下的仿真結果與試驗數據的標準差分別為215.82和330.00，Hill48屈服準則下的仿真結果標準差更小。

圖13 第一次迭代仿真后的結果Fig.13 Results after the first iteration of the simulation

圖14 第二次迭代仿真后的結果Fig.14 Results after the second iteration of the simulation

4 結論

1)按照單向拉伸試驗標準，得到DC04(厚度1 mm)薄板的0°、45°、90°載荷位移曲線，真應力-塑性應變曲線，各向屈服應力和塑性應變比。

2)合理簡化了單向拉伸仿真模型及顯示動力計算，仿真結果表明模型有效，且經分析達到了準靜態要求。

4)通過對比Hill48屈服準則和Mises屈服準狀態下的仿真結果，發現不同的屈服準則對損傷斷裂階段影響不同，Hill48各向異性屈服下的載荷位移曲線與試驗結果的一致性較各向同性Mises屈服準則更好。