LS-DYNA材料的二次开发汇总

ANSYS/LS-DYNATraining2001LS-DYNA材料的二次开发ANSYS/LS-DYNA专题培训ANSYS/LS-DYNATraining2001内容●二次开发环境●主程序及入口条件●开发材料的本构、子程序及求解输入文件描述●编译、运行新的求解器●开发Kelvin_voigt粘弹材料●用新材料模式做大变形分析ANSYS/LS-DYNATraining2001LS-DYNA二次开发基于FORTRAN环境在PC和UNIX平台下都需要进行连接编译，生成新的求解器Pc平台需安装digitalvisualfortran5.0或Microsoftpowerstation4.0提供的资源包括：Ls-dyna.f用户自定义本构子程序Ls-dyna.lib静态连接库Ls-dyna.dspdigitalFORTRANworkspace文件或MAKEFILE用于（包括主程序）Readme.txt说明文件二次开发环境ANSYS/LS-DYNATraining2001在UNIX平台提供如下资源：Makefile执行编译批处理文件（？）object文件（内含主程序）dyn21.f用户定义本构子程序支持的平台：940.2版DECNECIBMHPSGISUN950d版COMPAQCRAYSGIIBMLINUXSUN用户平台需安装FORTRAN77ANSYS/LS-DYNATraining2001二次开发如何实现？用户自定义的本构代替ls-dyna.f或dyn21.f中的相关本构描述LS-DYNA共提供10种*user_defined_material_model,由这些输入数据为自定义本构提供参数，完成分析在程序中使用的自定义subroutine要和Jobname.K中指定的相同在一次分析中，用户最多可同时使用10种自定义材料本构ANSYS/LS-DYNATraining2001主程序及入口条件c******************************************************************c|LS-DYNAmainprogramentry|c******************************************************************programlsdyna3dcalldyna3dstopendc******************************************************************ANSYS/LS-DYNATraining2001入口条件参数传递ccm(1)=young'smodulusccm(2)=poisson'sratioCcm(n)=用户在点K中给定的新本构参数ceps(1)=x应变增量ceps(2)=y应变增量ceps(3)=z应变增量ceps(4)=xy应变增量ceps(5)=yz应变增量ceps(6)=zx应变增量c单元类型etype:ceq.“brick”实体单元ceq.“shell”壳单元ceq.“beam”梁单元cctime=当前时间cdt1=当前时间步长ccapa=纵向剪切缩减因子csig(1)=x应力csig(2)=y应力csig(3)=z应力csig(4)=xy应力csig(5)=yz应力csig(6)=zx应力cchisv(1)=历史变量1chisv(2)=历史变量2…chisv(n)=历史变量nANSYS/LS-DYNATraining2001每个积分步、主程序提供如下这些已知量：6个应变增量可能涉及的历史变量hisv（n）单元类型的字符串当前时间当前时间步长用户在点K文件中给定如下参数：弹性模量波松比其它参数cm(n)Ls-dyna.f或dyn21.f应完成的工作：求出6个应力增量求出其它可能涉及的历史变量hisv（n）ANSYS/LS-DYNATraining2001参数说明由主程序提供的所有参数基于单元坐标系，计算得到的应力显然如此，之后由主程序将其转换到整体坐标系所有的历史变量在初始调用子程序时将置零能量计算完全由主程序完成ANSYS/LS-DYNATraining2001子程序举例c******************************************************************c|user-definedsubroutineexample|c******************************************************************subroutineumat41(cm,eps,sig,hisv,dt1,capa,etype,time)cisotropicelasticmaterial(sampleusersubroutine)cvariablesccm(1)=young'smodulusccm(2)=poisson'sratiocceps(1)=localxstrainincrementceps(2)=localystrainincrementceps(3)=localzstrainincrementceps(4)=localxystrainincrementceps(5)=localyzstrainincrementceps(6)=localzxstrainincrementANSYS/LS-DYNATraining2001csig(1)=localxstresscsig(2)=localystresscsig(3)=localzstresscsig(4)=localxystresscsig(5)=localyzstresscsig(6)=localzxstresschisv(1)=1sthistoryvariablechisv(2)=2ndhistoryvariablechisv(n)=nthhistoryvariablecdt1=currenttimestepsizeccapa=reductionfactorfortransverseshearcetype:ceq.brickforsolidelementsceq.shellforallshellelementsceq.beamforallbeamelementsctime=currentproblemtime.ANSYS/LS-DYNATraining2001character*(*)etypedimensioncm(*),eps(*),sig(*),hisv(*)cccomputeshearmodulus,gcg2=cm(1)/(1.+cm(2))g=.5*g2cif(etype.eq.'brick')thendavg=(-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-DYNATraining2001elseif(etype.eq.'shell')thencgc=capa*gq1=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))*q3davg=(-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.0sig(4)=sig(4)+g*eps(4)sig(5)=sig(5)+gc*eps(5)sig(6)=sig(6)+gc*eps(6)cANSYS/LS-DYNATraining2001elseif(etype.eq.'beam')thenq1=cm(1)*cm(2)/((1.0+cm(2))*(1.0-2.0*cm(2)))q3=q1+2.0*ggc=capa*gdeti=1./(q3*q3-q1*q1)c22i=q3*deti;c23i=-q1*detifac=(c22i+c23i)*q1eps(2)=-eps(1)*fac-sig(2)*c22i-sig(3)*c23ieps(3)=-eps(1)*fac-sig(2)*c23i-sig(3)*c22idavg=(-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.0sig(3)=0.0sig(4)=sig(4)+gc*eps(4)sig(5)=0.0sig(6)=sig(6)+gc*eps(6)endifreturnendANSYS/LS-DYNATraining2001对应的点K中的材料描述*MAT_USER_DEFINED_MATERIAL_MODELS$MIDROMTLMCNHVIORTHOIBULKIG17.890E-094140043$IVECTIFAIL00$P1(E)P2(NU)P3(G)P4(K)2.100E+053.000E-0180.769E+3175.0E+3ANSYS/LS-DYNATraining2001练习：在PC上开发并应用kelvin-voigt粘弹材料橡胶采用kelvin-voigt模型,本构方程由下式给定：=E0+E1（/t）其中：E0=0.6437Mpa,E1=0.0136Mpa·s;密度：4000Kg/m3应用此材料做大变形分析球直径10cm,下面由地板支撑，上部由一钢板在10ms将其到厚度为5cm的圆饼ANSYS/LS-DYNATraining2001步骤：得到LSTC公司提供的资源Ls-dyna.fLs-dyna.libLs-dyna.dsp打开digitvisualfortran在此环境中打开Ls-dyna.dsp，将ANSYS/LS-DYNATraining2001subroutineumat41(cm,eps,sig,hisv,dt1,capa,etype,time,d,s,t)character*(*)etypedimensioncm(*),eps(*),sig(*),hisv(*),d(6),s(6),t(6)Cccharacter*(*)etypedimensioncm(*),eps(*),sig(*),hisv(*),d(6),s(6),t(6)cg2=cm(1)/(1.+cm(2))g=.5*g2Ls-dyna.fANSYS/LS-DYNATraining2001if(etype.eq.'brick')thendavg=(-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)/dt1d(2)=eps(2)/dt1d(3)=eps(3)/dt1d(4)=eps(4)/dt1d(5)=eps(5)/dt1d(6)=eps(6)/dt1Cforthesecondtermintheconstitutiveg2=cm(5)/(1.+cm(2))davg=(-d(1)-d(2)-d(3))/3.ANSYS/LS-DYNATraining2001p=-davg*cm(5)/((1.-