李嘉文，李美求，李旭東，羅競波

(長江大學(xué)機械工程學(xué)院, 湖北 荊州 434023)

基于LS-DYNA的理想彈塑性材料本構(gòu)關(guān)系模型開發(fā)基本方法

李嘉文，李美求，李旭東，羅競波

(長江大學(xué)機械工程學(xué)院, 湖北 荊州 434023)

LS-DYNA是一款功能強大的有限元分析計算軟件，軟件自身含有強大的材料庫，更允許用戶添加自定義的材料本構(gòu)關(guān)系模型，用于得到更準(zhǔn)確的仿真結(jié)果。基于該平臺，介紹了該軟件與用戶自定義材料子程序之間的數(shù)據(jù)交互關(guān)系，并利用彈塑性力學(xué)相關(guān)知識，分別建立了理想彈塑性材料的彈性本構(gòu)關(guān)系與塑性本構(gòu)關(guān)系，最后通過算例驗證了所構(gòu)建材料本構(gòu)關(guān)系的正確性，并編寫了基于LS-DYNA平臺的理想彈塑性材料本構(gòu)關(guān)系，介紹了開發(fā)的基本方法，為應(yīng)用自編寫復(fù)雜本構(gòu)關(guān)系材料的工程仿真提供了基礎(chǔ)。

LS-DYNA；材料本構(gòu)關(guān)系；二次開發(fā)

LS-DYNA是一款功能齊全的幾何非線性、材料非線性以及摩擦和接觸分離等界面狀態(tài)非線性有限元數(shù)值計算軟件。凡是涉及接觸-碰撞、爆炸、穿甲與侵徹、應(yīng)力波傳播、金屬加工、流固耦合等問題，LS-DYNA都可以進(jìn)行求解。目前，LS-DYNA在國防軍工、巖土工程、土木工程、建筑、汽車等領(lǐng)域中均獲得廣泛應(yīng)用[1]。雖然LS-DYNA軟件自身擁有較完善的材料庫，涵蓋了彈性、彈塑性、復(fù)合材料、蜂窩材料、絲織物、膠類物、生物肌體、炸藥、推進(jìn)劑、混凝土、土壤、地質(zhì)、超彈性、橡膠、泡沫、玻璃、粘性流體、剛體等各種材料模型[2]，已經(jīng)可以滿足大部分的工程需要，但在實際工程中，針對特殊情況，用戶更希望使用自定義的材料本構(gòu)關(guān)系模型，如特殊的材料屈服條件，用以得到更準(zhǔn)確的仿真結(jié)果。

2006年，金乾坤在TCK模型的基礎(chǔ)上開發(fā)了一個用于混凝土或鋼筋混凝土碰撞或侵徹的損傷模型，并用有關(guān)混凝土侵徹試驗結(jié)果對模型參數(shù)進(jìn)行了標(biāo)定[3]。2011年，張安康等指出了采用數(shù)值算法求解本構(gòu)方程組時，半隱式的圖形返回算法具有良好的應(yīng)用性[4]。2015年，莊蔚敏等建立了高強鋼熱成形損傷-相變本構(gòu)關(guān)系模型,通過編寫該力學(xué)模型的用戶自定義材料子程序，將其應(yīng)用于LS-DYNA軟件的高強鋼熱成形數(shù)值模擬中，該試驗驗證了本構(gòu)模型的有效性[5]。2016年，胡志強等二次開發(fā)功能，將冰材料模型嵌入LS-DYNA程序，并驗證了該模型的準(zhǔn)確性和適用性[6]。下面,筆者主要從編寫材料本構(gòu)關(guān)系的基礎(chǔ)出發(fā)，介紹自定義本構(gòu)關(guān)系模型開發(fā)的基本方法。

1 用戶子程序的使用

圖1 UMAT文件工作簡化流程

LSTC公司在軟件中提供了用于用戶自定義二次開發(fā)的動態(tài)鏈接庫文件(LSDYNA.LIB)并預(yù)設(shè)了程序接口數(shù)組實現(xiàn)數(shù)據(jù)傳遞。

1.1UMAT文件工作流程

UMAT文件工作流程分為3步：①用戶需要利用以FORTRAN語言編寫的UMAT文件，結(jié)合LS-DYNA.LIB數(shù)據(jù)庫，在Fortran Compiler中編譯得到新的可執(zhí)行文件LSDYNA.EXE；②使用關(guān)鍵字*MAT_USER_DEFINED_MATERIAL_MODELS定義用戶材料參數(shù)后生成K文件；③求解前在LS-DYNA MANAGER中選擇Solver為生成的LSDYNA.EXE進(jìn)行求解，實現(xiàn)UMAT文件的使用。UMAT文件簡化流程圖如圖1所示。

1.2軟件與UMAT文件的接口參數(shù)

用戶需要在用F語言[7]編寫的UMAT文件開頭定義如下語言，作為數(shù)據(jù)接口，正確實現(xiàn)參數(shù)的相互傳遞，以LS-DYNA970版本為例，軟件預(yù)設(shè)的接口為[1]：

Subroutine umat41 (cm,eps,sig,hisv,dt,capa,etype,time,temp,i,ixs,x,k,j)

不僅限于Subroutine UMAT41，Subroutine umat41～umat50都是用戶可以自定義的材料，在LS-DYNA中，允許最多10個UMAT同時參與計算?，F(xiàn)介紹其中主要的數(shù)組，cm是用戶指定的材料參數(shù)數(shù)組；eps是當(dāng)前時步的6個應(yīng)變增量數(shù)組；sig是本子程序計算得到的6個應(yīng)力；hisv為用戶自定的歷史變量；dt是當(dāng)前時步的時間步長；capa為殼單元的橫向剪應(yīng)力系數(shù)；etype為單元類型；time是當(dāng)前的計算時間。

在LS-DYNA中，一個正確的UMAT文件的主要任務(wù)是定義材料的應(yīng)力-應(yīng)變關(guān)系[8]，其具體的工作方式是由軟件計算得到的應(yīng)變增量Δε，在UMAT文件中計算得到應(yīng)力增量Δσ，同時更新當(dāng)前應(yīng)力σ，并返回至軟件中參與下一個時間步的應(yīng)變增量計算。用戶在編寫UMAT文件時需要合理使用以上參數(shù)，實現(xiàn)數(shù)據(jù)在軟件與UMAT文件間的正確輸入與輸出。

1.3K文件關(guān)鍵字定義

在LS-DYNA中，用戶需要通過關(guān)鍵字*MAT_USER_DEFINED_MATERIAL_MODELS來定義用戶材料參數(shù)[1]。其中的主要控制參數(shù)如表1所示。

表1主要控制參數(shù)及說明

用戶需要在K文件中對應(yīng)參數(shù)位置輸入合理的數(shù)值，完成材料參數(shù)的定義。其中部分參數(shù)需要和UMAT文件結(jié)合使用，用戶定義這些量時需要注意，避免出現(xiàn)計算錯誤。

2 理想彈塑性材料本構(gòu)關(guān)系的編寫

理想彈塑性材料在受逐漸增大的力作用時，材料的應(yīng)力(σ)-應(yīng)變(ε)行為分為彈性階段和塑性階段(屈服階段)，在材料未進(jìn)入塑性階段時，材料可視為線彈性材料，其應(yīng)力應(yīng)變關(guān)系為：

σ=Eε

(1)

式中，E為彈性模量,MPa。

理想彈塑性材料在進(jìn)入塑性階段后，其屈服應(yīng)力σs不再增加，某種程度上，理想彈塑性材料的屈服就意味著破壞。其應(yīng)力應(yīng)變關(guān)系為：

σ=σs

(2)

基于以上理論，UMAT需要對材料的狀態(tài)階段進(jìn)行判斷，從而選擇不同的真實應(yīng)力σ′計算方式，涉及到用戶編程時對材料屈服條件的編寫。用戶可以根據(jù)具體需要編寫合適的屈服條件。以金屬材料為例，最常用的屈服準(zhǔn)則有Tresca屈服條件與Mises屈服條件2種。

Tresca屈服條件可以表示為：

(3)

式中，J2為試應(yīng)力偏量的第二不變量,MPa；θσ為Lode角。

Mises屈服條件則為：

(4)

用戶在UMAT中需要按彈性階段計算試應(yīng)力σ*，并以此計算屈服條件中的相關(guān)參數(shù)，與設(shè)定的屈服強度σs進(jìn)行比較，一旦試應(yīng)力滿足以上屈服條件，則認(rèn)為材料進(jìn)入了屈服階段。下面均以Mises屈服條件為判斷條件。

為了計算彈性試應(yīng)力σ*，需要用到彈性力學(xué)中的廣義胡克定律[9]：

(5)

為方便程序編寫，可利用其增量形式進(jìn)行表達(dá)，如式(6)所示：

(6)

其中，σx、σy、σz、τxy、τyz、τxz為6個應(yīng)力分量，依次對應(yīng)UMAT中的sig數(shù)組中的6個值,MPa； Δεx、Δεy、Δεz、Δγxy、Δγyz、Δγxz為6個應(yīng)變增量，依次對應(yīng)UMAT中的eps數(shù)組的6個值，1； Δεm為平均應(yīng)變增量，Δεm=(Δεx+Δεy+Δεz)/6；σm為平均應(yīng)力，σm=3KΔεm，K=E/3(1-2ν)為體積彈性模量；ν為材料泊松比；G為材料的剪切模量，MPa。

最終的試應(yīng)力數(shù)組σ*為：

圖2 應(yīng)力更新算法

σ*=σ0+ Δσ

(7)

式中，σ0為當(dāng)前時間步LS-DYNA輸入UMAT的應(yīng)力值，MPa。

σ′=σ*

(8)

若為假，則表明此時材料已經(jīng)進(jìn)入塑性階段，此時應(yīng)退回到初始應(yīng)力狀態(tài)σ0，即：

σ′=σ0

(9)

隨后，用戶需要將得到的應(yīng)力值依次返回至sig數(shù)組中，更新應(yīng)力狀態(tài)，為下一時間的計算做準(zhǔn)備。其實現(xiàn)應(yīng)力更新過程的算法如圖2所示。

3 算例驗證

試件為5mm×5mm×5mm立方體，加載方式為下端面全約束，上端面施加速率為1mm/s、方向沿Y軸負(fù)方向的速度載荷。材料彈性模量E=200GPa，泊松比ν=0.3，密度ρ=6000kg/m3，屈服應(yīng)力設(shè)置為σs=300MPa。使用ANSYS/LS-DYNA模塊進(jìn)行前處理[10]，使用Solid164單元劃分單元模型。如圖3所示。編譯UMAT文件后導(dǎo)入求解，提取該單元的應(yīng)力-時間、應(yīng)變-時間數(shù)據(jù)，可繪制單元的應(yīng)力-應(yīng)變曲線，如圖4所示。

圖3 試件有限元模型 圖4 單元應(yīng)力-應(yīng)變曲線

由圖4可以看出，當(dāng)應(yīng)變ε≤1.5×10-3時，應(yīng)力從零開始隨著應(yīng)變的增加，應(yīng)力-應(yīng)變關(guān)系為一定斜率直線，說明在該應(yīng)變下材料處于彈性階段，符合式(1)的表達(dá)形式。當(dāng)應(yīng)變ε>1.5×10-3時，此時的等效應(yīng)力σ=300MPa為一定值，不隨應(yīng)變的增加而變化，且該定值等于屈服應(yīng)力，符合式(2)的表達(dá)形式。材料的2個階段的表現(xiàn)均符合上述理論分析結(jié)果，由此推斷該理想彈塑性材料本構(gòu)關(guān)系的UMAT完全正確。

4 結(jié)論

1)介紹了LS-DYNA的用戶子程序的文件使用流程，并說明了UMAT文件與LS-DYNA軟件間是通過應(yīng)變增量-應(yīng)力關(guān)系來傳遞更新數(shù)據(jù)。

2)從理論上分析并建立了理想彈塑性材料本構(gòu)關(guān)系，材料在不同階段采用不同的應(yīng)力算法，給出了構(gòu)建本構(gòu)關(guān)系的基本的流程，用算例驗證了本構(gòu)關(guān)系的正確性，為復(fù)雜材料本構(gòu)關(guān)系的編寫提供了參考。

[1]白金澤.LS-DYNA3D理論基礎(chǔ)與實例分析(附光盤)[M].北京：科學(xué)出版社, 2005.

[2]劉成清, 何斌, 陳馳.ANSYS/LS-DYNA工程結(jié)構(gòu)抗震、抗撞擊與抗連續(xù)倒塌分析[M].北京：中國建筑工業(yè)出版社, 2014.

[3]金乾坤.混凝土動態(tài)損傷與失效模型[J].兵工學(xué)報, 2006, 27(1):10～14.

[4]張安康, 陳士海.LS-DYNA用戶自定義材料模型開發(fā)與驗證[J].計算機應(yīng)用與軟件, 2011, 28(4):71～73.

[5]莊蔚敏, 解東旋, 余天明,等.基于損傷-相變本構(gòu)模型的高強鋼熱成形數(shù)值模擬分析[J].吉林大學(xué)學(xué)報(工學(xué)版), 2015, 45(4):1206～1212.

[6]胡志強, 高巖, 姚琪.一種理想彈塑性模擬的冰材料本構(gòu)模型[J].船舶與海洋工程, 2016(1):65～73.

[7]黃曉梅, 張偉林.FORTRAN90程序設(shè)計[M].合肥：安徽大學(xué)出版社, 2016.

[8]師訪.ANSYS二次開發(fā)及應(yīng)用實例詳解[M].北京： 中國水利水電出版社, 2012.

[9]戴宏亮.彈塑性力學(xué)[M].長沙：湖南大學(xué)出版社, 2016.

[10]熊令芳, 胡凡金.ANSYS/LS-DYNA非線性動力分析方法與工程應(yīng)用[M].北京：中國鐵道出版社, 2015.

[編輯] 洪云飛

TG113.25

A

1673-1409(2017)17-0049-04

2017-06-18

中國石油科技創(chuàng)新基金項目計劃(2015D-5006-0310)。

李嘉文(1992-)，男，碩士生，現(xiàn)主要從事機械結(jié)構(gòu)強度方面的研究工作，519702602@qq.com。

[引著格式]李嘉文，李美求，李旭東，等.基于LS-DYNA的理想彈塑性材料本構(gòu)關(guān)系模型開發(fā)基本方法[J].長江大學(xué)學(xué)報(自科版)，2017,14(17)：49～52.