Marketing Strategy

Size: px

Start display at page:

Download "Marketing Strategy"

始申前匡
5 years ago
Views:

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

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

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

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

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

6 主程序及入口条件 ****************************************************************** LS-DYNA main program entry ****************************************************************** program lsdyna3d all dyna3d stop end ******************************************************************

7 入口条件参数传递 m(1)=young's modulus m(2)=poisson's ratio C m(n)= 用户在点 K 中给定的新本构参数 eps(1)= x 应变增量 eps(2)= y 应变增量 eps(3)= z 应变增量 eps(4)= xy 应变增量 eps(5)= yz 应变增量 eps(6)= zx 应变增量单元类型 etype: eq. brik 实体单元 eq. shell 壳单元 eq. beam 梁单元 time= 当前时间 dt1= 当前时间步长 apa= 纵向剪切缩减因子 sig(1)= x 应力 sig(2)=y 应力 sig(3)=z 应力 sig(4)=xy 应力 sig(5)=yz 应力 sig(6)=zx 应力 hisv(1)= 历史变量 1 hisv(2)= 历史变量 2 hisv(n)= 历史变量 n

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

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

10 子程序举例 ****************************************************************** user-defined subroutine example ****************************************************************** subroutine umat41 (m,eps,sig,hisv,dt1,apa,etype,time) isotropi elasti material (sample user subroutine) variables m(1)=young's modulus m(2)=poisson's ratio eps(1)=loal x strain inrement eps(2)=loal y strain inrement eps(3)=loal z strain inrement eps(4)=loal xy strain inrement eps(5)=loal yz strain inrement eps(6)=loal zx strain inrement

11 sig(1)=loal x stress sig(2)=loal y stress sig(3)=loal z stress sig(4)=loal xy stress sig(5)=loal yz stress sig(6)=loal zx stress hisv(1)=1st history variable hisv(2)=2nd history variable hisv(n)=nth history variable dt1=urrent time step size apa=redution fator for transverse shear etype: eq."brik" for solid elements eq."shell" for all shell elements eq."beam" for all beam elements time=urrent problem time.

12 harater*(*) etype dimension m(*),eps(*),sig(*),hisv(*) ompute shear modulus, g g2 =m(1)/(1.+m(2)) g =.5*g2 if (etype.eq.'brik') then davg=(-eps(1)-eps(2)-eps(3))/3. p=-davg*m(1)/((1.-2.*m(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)

13 elseif (etype.eq.'shell') then g =apa*g q1 =m(1)*m(2)/((1.0+m(2))*( *m(2))) q3 =1./(q1+g2) eps(3)=-q1*(eps(1)+eps(2))*q3 davg =(-eps(1)-eps(2)-eps(3))/3. p =-davg*m(1)/((1.-2.*m(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)+g*eps(5) sig(6)=sig(6)+g*eps(6)

14 elseif (etype.eq.'beam' ) then q1 =m(1)*m(2)/((1.0+m(2))*( *m(2))) q3 =q1+2.0*g g =apa*g deti =1./(q3*q3-q1*q1) 22i = q3*deti; 23i =-q1*deti fa =(22i+23i)*q1 eps(2)=-eps(1)*fa-sig(2)*22i-sig(3)*23i eps(3)=-eps(1)*fa-sig(2)*23i-sig(3)*22i davg =(-eps(1)-eps(2)-eps(3))/3. p =-davg*m(1)/((1.-2.*m(2))) sig(1)=sig(1)+p+g2*(eps(1)+davg) sig(2)=0.0 sig(3)=0.0 sig(4)=sig(4)+g*eps(4) sig(5)=0.0 sig(6)=sig(6)+g*eps(6) endif return ANSYS/LS-DYNA Training end 2001

15 对应的点 K 中的材料描述 *MAT_USER_DEFINED_MATERIAL_MODELS $ MID RO MT LMC NHV IORTHO IBULK IG E $ IVECT IFAIL 0 0 $ P1(E) P2(NU) P3(G) P4(K) 2.100E E E E+3

16 练习 : 在 PC 上开发并应用 kelvin-voigt 粘弹材料橡胶采用 kelvin-voigt 模型, 本构方程由下式给定 : σ= E 0 ε+e 1 ( ε/ t ) 其中 : E 0 =0.6437Mpa, E 1 = Mpa s; 密度 :4000Kg/m 3 应用此材料做大变形分析球直径 10m, 下面由地板支撑, 上部由一钢板在 10ms 将其到厚度为 5m 的圆饼

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

18 Ls-dyna.f subroutine umat41 (m,eps,sig,hisv,dt1,apa,etype,time,d,s,t) harater*(*) etype dimension m(*),eps(*),sig(*),hisv(*),d(6),s(6),t(6) C harater*(*) etype dimension m(*),eps(*),sig(*),hisv(*),d(6),s(6),t(6) g2 =m(1)/(1.+m(2)) g =.5*g2

19 if (etype.eq.'brik') then davg=(-eps(1)-eps(2)-eps(3))/3. p=-davg*m(1)/((1.-2.*m(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 seond term in the onstitutive g2 =m(5)/(1.+m(2)) davg=(-d(1)-d(2)-d(3))/3.

20 p=-davg*m(5)/((1.-2.*m(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 return end

21 Kelvin.k *MAT_USER_DEFINED_MATERIAL_MODELS $ MID RO MT LMC NHV IORTHO IBULK IG 2 4e e5.49 8e6 2e7 1.36e4

22 应用此材料模式求解大变形问题, 与 Moony_rivlin 计算结果比较 ( 略 )

3 92 59235 89 5 60,00024,000 84,000 2028476 6

山中傳奇與水的牽絆二峰圳 4 4 0 00 922 92 92002 60 2028476 3 92 59235 89 5 60,00024,000 84,000 2028476 6 86.6 679 58,800 4 220.4 970 23,800 69 58,800 69,733 4 23,800 252,478 69 76 6 925 4 7 5 970 38 50 25