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

李嘉文，李美求，李旭東，羅競(jìng)波

(長(zhǎng)江大學(xué)機(jī)械工程學(xué)院, 湖北 荊州 434023)

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

李嘉文，李美求，李旭東，羅競(jìng)波

(長(zhǎng)江大學(xué)機(jī)械工程學(xué)院, 湖北 荊州 434023)

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

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

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

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

1 用戶子程序的使用

圖1 UMAT文件工作簡(jiǎn)化流程

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

1.1UMAT文件工作流程

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

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

用戶需要在用F語(yǔ)言[7]編寫(xiě)的UMAT文件開(kāi)頭定義如下語(yǔ)言，作為數(shù)據(jù)接口，正確實(shí)現(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個(gè)UMAT同時(shí)參與計(jì)算。現(xiàn)介紹其中主要的數(shù)組，cm是用戶指定的材料參數(shù)數(shù)組；eps是當(dāng)前時(shí)步的6個(gè)應(yīng)變?cè)隽繑?shù)組；sig是本子程序計(jì)算得到的6個(gè)應(yīng)力；hisv為用戶自定的歷史變量；dt是當(dāng)前時(shí)步的時(shí)間步長(zhǎng)；capa為殼單元的橫向剪應(yīng)力系數(shù)；etype為單元類(lèi)型；time是當(dāng)前的計(jì)算時(shí)間。

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

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

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

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

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

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

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

σ=Eε

(1)

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

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

σ=σs

(2)

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

Tresca屈服條件可以表示為：

(3)

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

Mises屈服條件則為：

(4)

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

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

(5)

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

(6)

其中，σx、σy、σz、τxy、τyz、τxz為6個(gè)應(yīng)力分量，依次對(duì)應(yīng)UMAT中的sig數(shù)組中的6個(gè)值,MPa； Δεx、Δεy、Δεz、Δγxy、Δγyz、Δγxz為6個(gè)應(yīng)變?cè)隽?，依次?duì)應(yīng)UMAT中的eps數(shù)組的6個(gè)值，1； Δεm為平均應(yīng)變?cè)隽?，Δε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)前時(shí)間步LS-DYNA輸入U(xiǎn)MAT的應(yīng)力值，MPa。

σ′=σ*

(8)

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

σ′=σ0

(9)

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

3 算例驗(yàn)證

試件為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單元?jiǎng)澐謫卧Ｐ?。如圖3所示。編譯UMAT文件后導(dǎo)入求解，提取該單元的應(yīng)力-時(shí)間、應(yīng)變-時(shí)間數(shù)據(jù)，可繪制單元的應(yīng)力-應(yīng)變曲線，如圖4所示。

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

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

4 結(jié)論

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

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

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

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

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

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

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

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

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

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

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

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

[編輯] 洪云飛

TG113.25

A

1673-1409(2017)17-0049-04

2017-06-18

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

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

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