LS-DYNA 材料的二次开发

ANSYS/LS-DYNA专题培训 LS-DYNA材料的二次开发

内容 ● 二次开发环境 ● 主程序及入口条件 ● 开发材料的本构、子程序及求解输入文件描述 ● 编译、运行新的求解器 ● 开发Kelvin_voigt粘弹材料 ● 用新材料模式做大变形分析

二次开发环境 LS-DYNA二次开发基于FORTRAN环境在PC和UNIX平台下都需要进行连接编译，生成新的求解器 Pc平台需安装 digital visual fortran 5.0或Microsoft power station 4.0 提供的资源包括： Ls-dyna.f 用户自定义本构子程序 Ls-dyna.lib 静态连接库 Ls-dyna.dsp digital FORTRAN workspace 文件或MAKEFILE 用于（包括主程序） Readme.txt 说明文件

用户平台需安装FORTRAN77 支持的平台： 940.2 版 DEC NEC IBM HP SGI SUN 950d 版 COMPAQ CRAY SGI IBM LINUX SUN 在UNIX平台提供如下资源： Makefile 执行编译批处理文件（？） object 文件（内含主程序） dyn21.f 用户定义本构子程序

二次开发如何实现？ 用户自定义的本构代替ls-dyna.f或dyn21.f中的相关本构描述 LS-DYNA共提供10种*user_defined_material_model,由这些输入数据为自定义本构提供参数，完成分析在程序中使用的自定义subroutine要和Jobname.K中指定的相同在一次分析中，用户最多可同时使用10种自定义材料本构

主程序及入口条件 c****************************************************************** c| LS-DYNA main program entry | c****************************************************************** program lsdyna3d call dyna3d stop end c******************************************************************

入口条件 参数传递 c cm(1)=young's modulus c cm(2)=poisson's ratio C 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”壳单元 c eq.“beam”梁单元 c c 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应力 c c hisv(1)=历史变量1 c hisv(2)=历史变量2 … c hisv(n)=历史变量n

用户在点K文件中给定如下参数： 弹性模量波松比其它参数cm(n) Ls-dyna.f或dyn21.f应完成的工作：求出6个应力增量求出其它可能涉及的历史变量hisv（n）每个积分步、主程序提供如下这些已知量： 6个应变增量可能涉及的历史变量hisv（n）单元类型的字符串当前时间当前时间步长

参数说明 由主程序提供的所有参数基于单元坐标系，计算得到的应力显然如此，之后由主程序将其转换到整体坐标系 所有的历史变量在初始调用子程序时将置零 能量计算完全由主程序完成

子程序举例 c****************************************************************** c| user-defined subroutine example | c****************************************************************** subroutine umat41 (cm,eps,sig,hisv,dt1,capa,etype,time) c isotropic elastic material (sample user subroutine) c variables c cm(1)=young's modulus c cm(2)=poisson's ratio c c eps(1)=local x strain increment c eps(2)=local y strain increment c eps(3)=local z strain increment c eps(4)=local xy strain increment c eps(5)=local yz strain increment c eps(6)=local zx strain increment

c sig(1)=local x stress c sig(2)=local y stress c sig(3)=local z stress c sig(4)=local xy stress c sig(5)=local yz stress c sig(6)=local zx stress c hisv(1)=1st history variable c hisv(2)=2nd history variable c hisv(n)=nth history variable c dt1=current time step size c capa=reduction factor for transverse shear c etype: c eq."brick" for solid elements c eq."shell" for all shell elements c eq."beam" for all beam elements c time=current problem time.

character*(*) etype dimension cm(*),eps(*),sig(*),hisv(*) c c compute shear modulus, g c g2 =cm(1)/(1.+cm(2)) g =.5*g2 c 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)

elseif (etype.eq.'shell') then c 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)=sig(6)+gc*eps(6) c

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 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 end

对应的点K中的材料描述 *MAT_USER_DEFINED_MATERIAL_MODELS $ MID RO MT LMC NHV IORTHO IBULK 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+3

练习：在PC上开发并应用kelvin-voigt粘弹材料  橡胶采用kelvin-voigt模型,本构方程由下式给定： = E0 +E1（  /  t ）其中： E0 =0.6437Mpa, E1 = 0.0136Mpa·s; 密度：4000Kg/m3 应用此材料做大变形分析球直径10cm,下面由地板支撑，上部由一钢板在10ms将其到厚度为5cm的圆饼

步骤： • 得到LSTC公司提供的资源 • Ls-dyna.f • Ls-dyna.lib • Ls-dyna.dsp • 打开digit visual fortran 在此环境中打开Ls-dyna.dsp ，将

Ls-dyna.f subroutine umat41 (cm,eps,sig,hisv,dt1,capa,etype,time,d,s,t) character*(*) etype dimension cm(*),eps(*),sig(*),hisv(*),d(6),s(6),t(6) C c character*(*) etype dimension cm(*),eps(*),sig(*),hisv(*),d(6),s(6),t(6) c g2 =cm(1)/(1.+cm(2)) g =.5*g2

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)+davg) s(4)=g*eps(4) s(5)=g*eps(5) s(6)=g*eps(6) d(1)=eps(1)/dt1 d(2)=eps(2)/dt1 d(3)=eps(3)/dt1 d(4)=eps(4)/dt1 d(5)=eps(5)/dt1 d(6)=eps(6)/dt1 C for the second term in the constitutive g2 =cm(5)/(1.+cm(2)) davg=(-d(1)-d(2)-d(3))/3.

p=-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)=sig(6)+s(6)+t(6) endif c return end

Kelvin.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

应用此材料模式求解大变形问题，与Moony_rivlin计算结果比较（略）应用此材料模式求解大变形问题，与Moony_rivlin计算结果比较（略）

LS-DYNA 材料的二次开发