LS-DYNA材料的二次开发

nullnullLS-DYNA材料的 二次开发ANSYS/LS-DYNA专题 培训 焊锡培训资料ppt免费下载焊接培训教程 ppt 下载特设培训下载班长管理培训下载培训时间表下载 null内 容● 二次开发环境 ● 主程序及入口条件 ● 开发材料的本构、子程序及求解输入文件描述 ● 编译、运行新的求解器 ● 开发Kelvin_voigt粘弹材料 ● 用新材料模式做大变形分析 nullLS-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 说明文件二次开发环境null在UNIX平台提供如下资源： Makefile 执行编译批处理文件（？） object 文件（内含主程序） dyn21.f 用户定义本构子程序 支持的平台： 940.2 版 DEC NEC IBM HP SGI SUN 950d 版 COMPAQ CRAY SGI IBM LINUX SUN用户平台需安装FORTRAN77 null二次开发如何实现？ 用户自定义的本构代替ls-dyna.f或dyn21.f中的相关本构描述 LS-DYNA共提供10种*user_defined_material_model,由这些输入数据为自定义本构提供参数，完成分析 在程序中使用的自定义subroutine要和Jobname.K中指定的相同 在一次分析中，用户最多可同时使用10种自定义材料本构null主程序及入口条件 c****************************************************************** c| LS-DYNA main program entry | c****************************************************************** program lsdyna3d call dyna3d stop end c****************************************************************** null入口条件参数传递 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 null每个积分步、主程序提供如下这些已知量：6个应变增量 可能涉及的历史变量hisv（n） 单元类型的字符串 当前时间 当前时间步长用户在点K文件中给定如下参数： 弹性模量 波松比 其它参数cm(n) Ls-dyna.f或dyn21.f应完成的工作： 求出6个应力增量 求出其它可能涉及的历史变量hisv（n） null参数说明由主程序提供的所有参数基于单元坐标系，计算得到的应力显然如此，之后由主程序将其转换到整体坐标系 所有的历史变量在初始调用子程序 时将置零 能量计算完全由主程序完成null子程序举例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 incrementnullc 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. null 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) null 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) cnull 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 endnull对应的点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 null练习：在PC上开发并应用kelvin-voigt粘弹材料  橡胶采用kelvin-voigt模型,本构方程由下式给定： = E0 +E1（  /  t ） 其中： E0 =0.6437Mpa, E1 = 0.0136Mpa·s; 密度：4000Kg/m3 应用此材料做大变形分析 球直径10cm,下面由地板支撑，上部由一钢板在10ms将其到厚度为5cm的圆饼 null步骤：得到LSTC公司提供的资源 Ls-dyna.f Ls-dyna.lib Ls-dyna.dsp 打开digit visual fortran 在此环境中打开Ls-dyna.dsp ，将null 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 Ls-dyna.fnull 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. nullp=-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 endnullKelvin.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 null应用此材料模式求解大变形问题，与Moony_rivlin计算结果比较（略）