LS-DYNA材料的二次開發(fā)

1、ANSYS/LS-DYNA Training 2001LS-DYNA材料的材料的二次開發(fā)二次開發(fā)ANSYS/LS-DYNA專題培訓(xùn)專題培訓(xùn)ANSYS/LS-DYNA Training 2001內(nèi)內(nèi) 容容 二次開發(fā)環(huán)境 主程序及入口條件 開發(fā)材料的本構(gòu)、子程序及求解輸入文件描述 編譯、運行新的求解器 開發(fā)Kelvin_voigt粘彈材料 用新材料模式做大變形分析ANSYS/LS-DYNA Training 2001LS-DYNA二次開發(fā)基于FORTRAN環(huán)境在PC和UNIX平臺下都需要進行連接編譯，生成新的求解器Pc平臺需安裝 digital visual fortran 5.0或Microso

2、ft power station 4.0提供的資源包括：Ls-dyna.f 用戶自定義本構(gòu)子程序Ls-dyna.lib 靜態(tài)連接庫Ls-dyna.dsp digital FORTRAN workspace 文件或MAKEFILE 用于（包括主程序）Readme.txt 說明文件二次開發(fā)環(huán)境二次開發(fā)環(huán)境ANSYS/LS-DYNA Training 2001在UNIX平臺提供如下資源：Makefile 執(zhí)行編譯批處理文件（？）object 文件（內(nèi)含主程序）dyn21.f 用戶定義本構(gòu)子程序支持的平臺：940.2 版 DEC NEC IBMHP SGI SUN950d 版COMPAQ CRAY S

3、GIIBM LINUX SUN用戶平臺需安裝FORTRAN77ANSYS/LS-DYNA Training 2001二次開發(fā)如何實現(xiàn)？二次開發(fā)如何實現(xiàn)？用戶自定義的本構(gòu)代替ls-dyna.f或dyn21.f中的相關(guān)本構(gòu)描述LS-DYNA共提供10種*user_defined_material_model,由這些輸入數(shù)據(jù)為自定義本構(gòu)提供參數(shù)，完成分析在程序中使用的自定義subroutine要和Jobname.K中指定的相同在一次分析中，用戶最多可同時使用10種自定義材料本構(gòu)ANSYS/LS-DYNA Training 2001主程序及入口條件主程序及入口條件c*c| LS-DYNA main p

4、rogram entry |c* program lsdyna3d call dyna3d stop endc*ANSYS/LS-DYNA Training 2001入口條件入口條件參數(shù)傳遞c cm(1)=youngs modulusc cm(2)=poissons ratioC cm(n)=用戶在點K中給定的新本構(gòu)參數(shù)c eps(1)= x應(yīng)變增量c eps(2)= y應(yīng)變增量c eps(3)= z應(yīng)變增量c eps(4)= xy應(yīng)變增量c eps(5)= yz應(yīng)變增量c eps(6)= zx應(yīng)變增量c 單元類型 etype:c eq.“brick” 實體單元c eq.“shell” 殼單元

5、c eq.“beam” 梁單元cc time=當(dāng)前時間c dt1=當(dāng)前時間步長c capa=縱向剪切縮減因子c sig(1)= x應(yīng)力c sig(2)=y應(yīng)力c sig(3)=z應(yīng)力c sig(4)=xy應(yīng)力c sig(5)=yz應(yīng)力c sig(6)=zx應(yīng)力cc hisv(1)=歷史變量1c hisv(2)=歷史變量2c hisv(n)=歷史變量nANSYS/LS-DYNA Training 2001每個積分步、主程序提供如下這些已知量：6個應(yīng)變增量可能涉及的歷史變量hisv（n）單元類型的字符串當(dāng)前時間當(dāng)前時間步長用戶在點K文件中給定如下參數(shù)：彈性模量波松比其它參數(shù)cm(n)Ls-dyna

6、.f或dyn21.f應(yīng)完成的工作：求出6個應(yīng)力增量求出其它可能涉及的歷史變量hisv（n）ANSYS/LS-DYNA Training 2001參數(shù)說明參數(shù)說明由主程序提供的所有參數(shù)基于單元坐標(biāo)系，計算得到的應(yīng)力顯然如此，之后由主程序?qū)⑵滢D(zhuǎn)換到整體坐標(biāo)系所有的歷史變量在初始調(diào)用子程序時將置零能量計算完全由主程序完成ANSYS/LS-DYNA Training 2001子程序舉例子程序舉例c*c| user-defined subroutine example |c* subroutine umat41 (cm,eps,sig,hisv,dt1,capa,etype,time)c isotrop

7、ic elastic material (sample user subroutine)c variablesc cm(1)=youngs modulusc cm(2)=poissons ratiocc eps(1)=local x strain incrementc eps(2)=local y strain incrementc eps(3)=local z strain incrementc eps(4)=local xy strain incrementc eps(5)=local yz strain incrementc eps(6)=local zx strain incremen

8、tANSYS/LS-DYNA Training 2001c sig(1)=local x stressc sig(2)=local y stressc sig(3)=local z stressc sig(4)=local xy stressc sig(5)=local yz stressc sig(6)=local zx stressc hisv(1)=1st history variablec hisv(2)=2nd history variablec hisv(n)=nth history variablec dt1=current time step sizec capa=reduct

9、ion factor for transverse shearc etype:c eq.brick for solid elementsc eq.shell for all shell elementsc eq.beam for all beam elementsc time=current problem time.ANSYS/LS-DYNA Training 2001 character*(*) etype dimension cm(*),eps(*),sig(*),hisv(*)cc compute shear modulus, gc g2 =cm(1)/(1.+cm(2) g =.5*

10、g2c if (etype.eq.brick) then davg=(-eps(1)-eps(2)-eps(3)/3. p=-davg*cm(1)/(1.-2.*cm(2) sig(1)=sig(1)+p+g2*(eps(1)+davg) sig(2)=sig(2)+p+g2*(eps(2)+davg) sig(3)=sig(3)+p+g2*(eps(3)+davg) sig(4)=sig(4)+g*eps(4) sig(5)=sig(5)+g*eps(5) sig(6)=sig(6)+g*eps(6)ANSYS/LS-DYNA Training 2001 elseif (etype.eq.s

11、hell) thenc gc =capa*g q1 =cm(1)*cm(2)/(1.0+cm(2)*(1.0-2.0*cm(2) q3 =1./(q1+g2) eps(3)=-q1*(eps(1)+eps(2)*q3 davg =(-eps(1)-eps(2)-eps(3)/3. p =-davg*cm(1)/(1.-2.*cm(2) sig(1)=sig(1)+p+g2*(eps(1)+davg) sig(2)=sig(2)+p+g2*(eps(2)+davg) sig(3)=0.0 sig(4)=sig(4)+g *eps(4) sig(5)=sig(5)+gc*eps(5) sig(6)

12、=sig(6)+gc*eps(6)cANSYS/LS-DYNA Training 2001 elseif (etype.eq.beam ) then q1 =cm(1)*cm(2)/(1.0+cm(2)*(1.0-2.0*cm(2) q3 =q1+2.0*g gc =capa*g deti =1./(q3*q3-q1*q1) c22i = q3*deti; c23i =-q1*deti fac =(c22i+c23i)*q1 eps(2)=-eps(1)*fac-sig(2)*c22i-sig(3)*c23i eps(3)=-eps(1)*fac-sig(2)*c23i-sig(3)*c22i

13、 davg =(-eps(1)-eps(2)-eps(3)/3. p =-davg*cm(1)/(1.-2.*cm(2) sig(1)=sig(1)+p+g2*(eps(1)+davg) sig(2)=0.0 sig(3)=0.0 sig(4)=sig(4)+gc*eps(4) sig(5)=0.0 sig(6)=sig(6)+gc*eps(6) endif return endANSYS/LS-DYNA Training 2001對應(yīng)的點對應(yīng)的點K中的材料描述中的材料描述*MAT_USER_DEFINED_MATERIAL_MODELS$ MID RO MT LMC NHV IORTHO I

14、BULK IG 1 7.890E-09 41 4 0 0 4 3$ IVECT IFAIL 0 0$ P1(E) P2(NU) P3(G) P4(K) 2.100E+05 3.000E-01 80.769E+3 175.0E+3ANSYS/LS-DYNA Training 2001練習(xí)：在練習(xí)：在PC上開發(fā)并應(yīng)用上開發(fā)并應(yīng)用kelvin-voigt粘粘彈材料彈材料 橡膠采用kelvin-voigt模型,本構(gòu)方程由下式給定： = E0 +E1（ / t ）其中： E0 =0.6437Mpa, E1 = 0.0136Mpas;密度：4000Kg/m3應(yīng)用此材料做大變形分析 球直徑10cm,下面由地

15、板支撐，上部由一鋼板在10ms將其到厚度為5cm的圓餅ANSYS/LS-DYNA Training 2001步驟：得到LSTC公司提供的資源 Ls-dyna.f Ls-dyna.lib Ls-dyna.dsp 打開digit visual fortran 在此環(huán)境中打開Ls-dyna.dsp ，將ANSYS/LS-DYNA Training 2001 subroutine umat41 (cm,eps,sig,hisv,dt1,capa,etype,time,d,s,t) character*(*) etype dimension cm(*),eps(*),sig(*),hisv(*),d(6

16、),s(6),t(6)Cc character*(*) etype dimension cm(*),eps(*),sig(*),hisv(*),d(6),s(6),t(6)c g2 =cm(1)/(1.+cm(2) g =.5*g2Ls-dyna.fANSYS/LS-DYNA Training 2001 if (etype.eq.brick) then davg=(-eps(1)-eps(2)-eps(3)/3. p=-davg*cm(1)/(1.-2.*cm(2) s(1)=p+g2*(eps(1)+davg) s(2)=p+g2*(eps(2)+davg) s(3)=p+g2*(eps(3

17、)+davg) s(4)=g*eps(4) s(5)=g*eps(5) s(6)=g*eps(6)d(1)=eps(1)/dt1d(2)=eps(2)/dt1d(3)=eps(3)/dt1d(4)=eps(4)/dt1 d(5)=eps(5)/dt1 d(6)=eps(6)/dt1C for the second term in the constitutive g2 =cm(5)/(1.+cm(2) davg=(-d(1)-d(2)-d(3)/3.ANSYS/LS-DYNA Training 2001p=-davg*cm(5)/(1.-2.*cm(2) t(1)=p+g2*(d(1)+davg) t(2)=p+g2*(d(2)+davg) t(3)=p+g2*(d(3)+davg) t(4)=g*d(4) t(5)=g*d(5) t(6)=g*d(6) sig(1)=sig(1)+s(1)+t(1) sig(2)=sig(2)+s(2)+t(2) sig(3)=sig(3)+s(3)+t(3) sig(4)=sig(4)+s(4)+t(4) sig(5)=sig(5)+s(5)+t(5) sig(6