LS-DYNA材料的二次开发
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
ANSYS/LS-DYNA专题培训
LS-DYNA材料的 二次开发
ANSYS/LS-DYNA Training 2001
内容
二次开发环境
主程序及入口条件 开发材料的本构、子程序及求解输入文件描述 编译、运行新的求解器 开发Kelvin_voigt粘弹材料 用新材料模式做大变形分析
● ● ● ● ● ●
c****************************************************************** program lsdyna3d call dyna3d stop end c******************************************************************
ANSYS/LS-DYNA Training 2001
elseif (etype.eq.'beam' ) then q1 q3 gc =cm(1)*cm(2)/((1.0+cm(2))*(1.0-2.0*cm(2))) =q1+2.0*g =capa*g
deti =1./(q3*q3-q1*q1) c22i = q3*deti; fac =(c22i+c23i)*q1 c23i =-q1*deti
ANSYS/LS-DYNA Training 2001
二次开发环境 LS-DYNA二次开发基于FORTRAN环境 在PC和UNIX平台下都需要进行连接编译， 生成新的求解器 Pc平台需安装 digital visual fortran 5.0或 Microsoft power station 4.0 提供的资源包括：
1 7.890E-09 $ IVECT IFAIL
0
$ P1(E)
0
P2(NU) P3(G) P4(K)
2.100E+05 3.000E-01 80.769E+3 175.0E+3
ANSYS/LS-DYNA Training 2001
练习：在PC上开发并应用kelvin-voigt粘 弹材料
橡胶采用kelvin-voigt模型,本构方程由下式给定：
ANSYS/LS-DYNA Training 2001
入口条件 参数传递
c c C c c c c c c cm(1)=young's modulus cm(2)=poisson's ratio cm(n)=用户在点K中给定的新本构参数 eps(1)= x应变增量 eps(2)= y应变增量 eps(3)= z应变增量 eps(4)= xy应变增量 eps(5)= yz应变增量 eps(6)= zx应变增量
ANSYS/LS-DYNA Training 2001
c c c c c
sig(1)=local x sig(2)=local y sig(3)=local z
stress stress stress
sig(4)=local xy stress sig(5)=local yz stress
c
c c c
sig(6)=local zx stress
hisv(1)=1st history variable hisv(2)=2nd history variable hisv(n)=nth history variable
c c c c c c c
dt1=current time step size capa=reduction factor for transverse shear etype: eq."brick" for solid elements eq."shell" for all shell elements eq."beam" for all beam elements time=current problem time.
ANSYS/LS-DYNA Training 2001
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
ANSYS/LS-DYNA Training 2001
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)
支持的平台： 940.2 版
DEC NEC HP SGI IBMFra Baidu bibliotekSUN
在UNIX平台提供如下资源：
Makefile 执行编译批处理文件（？）
object 文件（内含主程序） dyn21.f 用户定义本构子程序
950d 版
COMPAQ CRAY SGI IBM LINUX SUN
ANSYS/LS-DYNA Training 2001
ANSYS/LS-DYNA Training 2001
子程序举例
c****************************************************************** c| user-defined subroutine example |
c******************************************************************
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
= E0 +E1（ / t ） 其中： E0 =0.6437Mpa, E1 = 0.0136Mpa· s; 密度：4000Kg/m3
应用此材料做大变形分析
球直径10cm,下面由地板支撑，上部由一钢板在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
Ls-dyna.f
subroutine umat41 (cm,eps,sig,hisv,dt1,capa,etype,time,d,s,t)
ANSYS/LS-DYNA Training 2001
elseif (etype.eq.'shell') then c gc =capa*g
q1
q3
=cm(1)*cm(2)/((1.0+cm(2))*(1.0-2.0*cm(2)))
=1./(q1+g2)
eps(3)=-q1*(eps(1)+eps(2))*q3 davg =(-eps(1)-eps(2)-eps(3))/3.
二次开发如何实现？
用户自定义的本构代替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 program entry |
Ls-dyna.f 用户自定义本构子程序 Ls-dyna.lib 静态连接库 Ls-dyna.dsp digital FORTRAN workspace 文件或MAKEFILE 用于 （包括主程序） Readme.txt 说明文件 ANSYS/LS-DYNA Training 2001
用户平台需安装 FORTRAN77
subroutine umat41 (cm,eps,sig,hisv,dt1,capa,etype,time) c isotropic elastic material (sample user subroutine)
c variables
c
c c c c c c c c
cm(1)=young's modulus
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 ANSYS/LS-DYNA Training 2001
对应的点K中的材料描述
*MAT_USER_DEFINED_MATERIAL_MODELS $ MID RO MT 41 LMC 4 NHV 0 IORTHO 0 IBULK 4 IG 3
c c c c c c c c c … c
sig(1)= x应力 sig(2)=y应力 sig(3)=z应力 sig(4)=xy应力 sig(5)=yz应力 sig(6)=zx应力
hisv(1)=历史变量1 hisv(2)=历史变量2
c c c c c c c
单元类型 etype: eq.“brick” 实体单元 eq.“shell” 壳单元 eq.“beam” 梁单元
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)
hisv(n)=历史变量n
time=当前时间 dt1=当前时间步长
c
capa=纵向剪切缩减因子
ANSYS/LS-DYNA Training 2001
用户在点K文件中 给定如下参数：
弹性模量 波松比 其它参数cm(n)
Ls-dyna.f或dyn21.f应 完成的工作：
求出6个应力增量
求出其它可能涉及的历 史变量hisv（n）
cm(2)=poisson's ratio
eps(1)=local x eps(2)=local y eps(3)=local z
strain increment strain increment strain increment
eps(4)=local xy strain increment eps(5)=local yz strain increment eps(6)=local zx strain increment
每个积分步、主程序提供 如下这些已知量：
6个应变增量 可能涉及的历史变量hisv（n）
单元类型的字符串
当前时间 当前时间步长
ANSYS/LS-DYNA Training 2001
参数说明
由主程序提供的所有参数基于单元坐标 系，计算得到的应力显然如此，之后由 主程序将其转换到整体坐标系 所有的历史变量在初始调用子程序 时将置零 能量计算完全由主程序完成
LS-DYNA材料的 二次开发
ANSYS/LS-DYNA Training 2001
内容
二次开发环境
主程序及入口条件 开发材料的本构、子程序及求解输入文件描述 编译、运行新的求解器 开发Kelvin_voigt粘弹材料 用新材料模式做大变形分析
● ● ● ● ● ●
c****************************************************************** program lsdyna3d call dyna3d stop end c******************************************************************
ANSYS/LS-DYNA Training 2001
elseif (etype.eq.'beam' ) then q1 q3 gc =cm(1)*cm(2)/((1.0+cm(2))*(1.0-2.0*cm(2))) =q1+2.0*g =capa*g
deti =1./(q3*q3-q1*q1) c22i = q3*deti; fac =(c22i+c23i)*q1 c23i =-q1*deti
ANSYS/LS-DYNA Training 2001
二次开发环境 LS-DYNA二次开发基于FORTRAN环境 在PC和UNIX平台下都需要进行连接编译， 生成新的求解器 Pc平台需安装 digital visual fortran 5.0或 Microsoft power station 4.0 提供的资源包括：
1 7.890E-09 $ IVECT IFAIL
0
$ P1(E)
0
P2(NU) P3(G) P4(K)
2.100E+05 3.000E-01 80.769E+3 175.0E+3
ANSYS/LS-DYNA Training 2001
练习：在PC上开发并应用kelvin-voigt粘 弹材料
橡胶采用kelvin-voigt模型,本构方程由下式给定：
ANSYS/LS-DYNA Training 2001
入口条件 参数传递
c c C c c c c c c cm(1)=young's modulus cm(2)=poisson's ratio cm(n)=用户在点K中给定的新本构参数 eps(1)= x应变增量 eps(2)= y应变增量 eps(3)= z应变增量 eps(4)= xy应变增量 eps(5)= yz应变增量 eps(6)= zx应变增量
ANSYS/LS-DYNA Training 2001
c c c c c
sig(1)=local x sig(2)=local y sig(3)=local z
stress stress stress
sig(4)=local xy stress sig(5)=local yz stress
c
c c c
sig(6)=local zx stress
hisv(1)=1st history variable hisv(2)=2nd history variable hisv(n)=nth history variable
c c c c c c c
dt1=current time step size capa=reduction factor for transverse shear etype: eq."brick" for solid elements eq."shell" for all shell elements eq."beam" for all beam elements time=current problem time.
ANSYS/LS-DYNA Training 2001
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
ANSYS/LS-DYNA Training 2001
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)
支持的平台： 940.2 版
DEC NEC HP SGI IBMFra Baidu bibliotekSUN
在UNIX平台提供如下资源：
Makefile 执行编译批处理文件（？）
object 文件（内含主程序） dyn21.f 用户定义本构子程序
950d 版
COMPAQ CRAY SGI IBM LINUX SUN
ANSYS/LS-DYNA Training 2001
ANSYS/LS-DYNA Training 2001
子程序举例
c****************************************************************** c| user-defined subroutine example |
c******************************************************************
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
= E0 +E1（ / t ） 其中： E0 =0.6437Mpa, E1 = 0.0136Mpa· s; 密度：4000Kg/m3
应用此材料做大变形分析
球直径10cm,下面由地板支撑，上部由一钢板在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
Ls-dyna.f
subroutine umat41 (cm,eps,sig,hisv,dt1,capa,etype,time,d,s,t)
ANSYS/LS-DYNA Training 2001
elseif (etype.eq.'shell') then c gc =capa*g
q1
q3
=cm(1)*cm(2)/((1.0+cm(2))*(1.0-2.0*cm(2)))
=1./(q1+g2)
eps(3)=-q1*(eps(1)+eps(2))*q3 davg =(-eps(1)-eps(2)-eps(3))/3.
二次开发如何实现？
用户自定义的本构代替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 program entry |
Ls-dyna.f 用户自定义本构子程序 Ls-dyna.lib 静态连接库 Ls-dyna.dsp digital FORTRAN workspace 文件或MAKEFILE 用于 （包括主程序） Readme.txt 说明文件 ANSYS/LS-DYNA Training 2001
用户平台需安装 FORTRAN77
subroutine umat41 (cm,eps,sig,hisv,dt1,capa,etype,time) c isotropic elastic material (sample user subroutine)
c variables
c
c c c c c c c c
cm(1)=young's modulus
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 ANSYS/LS-DYNA Training 2001
对应的点K中的材料描述
*MAT_USER_DEFINED_MATERIAL_MODELS $ MID RO MT 41 LMC 4 NHV 0 IORTHO 0 IBULK 4 IG 3
c c c c c c c c c … c
sig(1)= x应力 sig(2)=y应力 sig(3)=z应力 sig(4)=xy应力 sig(5)=yz应力 sig(6)=zx应力
hisv(1)=历史变量1 hisv(2)=历史变量2
c c c c c c c
单元类型 etype: eq.“brick” 实体单元 eq.“shell” 壳单元 eq.“beam” 梁单元
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)
hisv(n)=历史变量n
time=当前时间 dt1=当前时间步长
c
capa=纵向剪切缩减因子
ANSYS/LS-DYNA Training 2001
用户在点K文件中 给定如下参数：
弹性模量 波松比 其它参数cm(n)
Ls-dyna.f或dyn21.f应 完成的工作：
求出6个应力增量
求出其它可能涉及的历 史变量hisv（n）
cm(2)=poisson's ratio
eps(1)=local x eps(2)=local y eps(3)=local z
strain increment strain increment strain increment
eps(4)=local xy strain increment eps(5)=local yz strain increment eps(6)=local zx strain increment
每个积分步、主程序提供 如下这些已知量：
6个应变增量 可能涉及的历史变量hisv（n）
单元类型的字符串
当前时间 当前时间步长
ANSYS/LS-DYNA Training 2001
参数说明
由主程序提供的所有参数基于单元坐标 系，计算得到的应力显然如此，之后由 主程序将其转换到整体坐标系 所有的历史变量在初始调用子程序 时将置零 能量计算完全由主程序完成