1、ANSYS/LS-DYNA Training 2001LS-DYNA材料的材料的二次开发二次开发ANSYS/LS-DYNA专题培训专题培训ANSYS/LS-DYNA Training 2001内内 容容 二次开发环境 主程序及入口条件 开发材料的本构、子程序及求解输入文件描述 编译、运行新的求解器 开发Kelvin_voigt粘弹材料 用新材料模式做大变形分析ANSYS/LS-DYNA Training 2001LS-DYNA二次开发基于FORTRAN环境在PC和UNIX平台下都需要进行连接编译，生成新的求解器Pc平台需安装 digital visual fortran 5.0或Microso
2、ft power station 4.0提供的资源包括：Ls-dyna.f 用户自定义本构子程序Ls-dyna.lib 静态连接库Ls-dyna.dsp digital FORTRAN workspace 文件或MAKEFILE 用于（包括主程序）Readme.txt 说明文件二次开发环境二次开发环境ANSYS/LS-DYNA Training 2001在UNIX平台提供如下资源：Makefile 执行编译批处理文件（？）object 文件（内含主程序）dyn21.f 用户定义本构子程序支持的平台：940.2 版 DEC NEC IBMHP SGI SUN950d 版COMPAQ CRAY S
3、GIIBM LINUX SUN用户平台需安装FORTRAN77ANSYS/LS-DYNA Training 2001二次开发如何实现？二次开发如何实现？用户自定义的本构代替ls-dyna.f或dyn21.f中的相关本构描述LS-DYNA共提供10种*user_defined_material_model,由这些输入数据为自定义本构提供参数，完成分析在程序中使用的自定义subroutine要和Jobname.K中指定的相同在一次分析中，用户最多可同时使用10种自定义材料本构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入口条件入口条件参数传递c cm(1)=youngs modulusc cm(2)=poissons ratioC cm(n)=用户在点K中给定的新本构参数c eps(1)= x应变增量c eps(2)= y应变增量c eps(3)= z应变增量c eps(4)= xy应变增量c eps(5)= yz应变增量c eps(6)= zx应变增量c 单元类型 etype:c eq.“brick” 实体单元c eq.“shell” 壳单元
5、c eq.“beam” 梁单元cc time=当前时间c dt1=当前时间步长c capa=纵向剪切缩减因子c sig(1)= x应力c sig(2)=y应力c sig(3)=z应力c sig(4)=xy应力c sig(5)=yz应力c sig(6)=zx应力cc hisv(1)=历史变量1c hisv(2)=历史变量2c hisv(n)=历史变量nANSYS/LS-DYNA Training 2001每个积分步、主程序提供如下这些已知量：6个应变增量可能涉及的历史变量hisv（n）单元类型的字符串当前时间当前时间步长用户在点K文件中给定如下参数：弹性模量波松比其它参数cm(n)Ls-dyna
6、.f或dyn21.f应完成的工作：求出6个应力增量求出其它可能涉及的历史变量hisv（n）ANSYS/LS-DYNA Training 2001参数说明参数说明由主程序提供的所有参数基于单元坐标系，计算得到的应力显然如此，之后由主程序将其转换到整体坐标系所有的历史变量在初始调用子程序时将置零能量计算完全由主程序完成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对应的点对应的点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练习：在练习：在PC上开发并应用上开发并应用kelvin-voigt粘粘弹材料弹材料 橡胶采用kelvin-voigt模型,本构方程由下式给定： = E0 +E1（ / t ）其中： E0 =0.6437Mpa, E1 = 0.0136Mpas;密度：4000Kg/m3应用此材料做大变形分析 球直径10cm,下面由地
15、板支撑，上部由一钢板在10ms将其到厚度为5cm的圆饼ANSYS/LS-DYNA Training 2001步骤：得到LSTC公司提供的资源 Ls-dyna.f Ls-dyna.lib Ls-dyna.dsp 打开digit visual fortran 在此环境中打开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)+dav
18、g) 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)=sig(6)+s(6)+t(6) endifc return endANSYS/LS-DYNA Training 2001Kelvin.k*MAT_USER_DEFINED_MATERIAL_MODELS$ MID RO MT LMC NHV IORTHO IBULK IG 2 4e3 41 5 0 0 4 3 0 0 6.437e5 .49 8e6 2e7 1.36e4 ANSYS/LS-DYNA Training 2001应用此材料模式求解大变形问题，与应用此材料模式求解大变形问题，与Moony_rivlin计算结果比较（略）计算结果比较（略）
