您所在位置：主题：Fortran语言编写金属材料的自定义本构方程，请大家帮我分析，谢谢！
lijjie
[专家分：0] 发布于 2011-07-21 21:08:00
大家好，我因为要用Fortran语言编写金属材料的自定义本构方程，需要把一个程序里面的参数给换掉就可以了，这个程序是是Advantedge软件里面的材料自定义模板。下面我把这个程序放上来，希望大家帮我分析一下这个程序是什么意思，软件应该把哪些部分给替换成自己的材料参数，不胜感激！（字母C后面都是注解，最好看附件吧，有字体颜色，看的清楚。红色字体代表程序的注解，蓝色和黑色字体才是真正的程序内容！）另外说明一下，我是机械专业的，估计有些内容涉及到材料专业的知识。分析完了最好给我发邮件注解一些程序的意思和应该替换的材料内容，494545750@126.com。
C
SUBROUTINE MAT_USER(d,time,dtime,temp,engine_s,user_s,
1 eps,deps,sig,df)
c
c AdvantEdge 2D/3D User defined material
c Total formulation example
c
c COPYRIGHT Third Wave Systems Inc. 2010
c www.thirdwavesys.com
c
#ifndef linux
!DEC$ ATTRIBUTES DLLEXPORT :: MAT_USER
#endif
c
c
implicit real*8 (a-h,o-z)
c
c
c Total formulation
c
c sig(3,3) Previous and updated stress [Pa]
c sig(1,1) = Sxx, sig(2,2)=Syy, sig(1,2)=Sxy
c dtime time step [s]
c temp temperature [c]
c
c engine_s(1) Plastic strain [1]
c engine_s(3) Plastic work rate (heat generation) [w]
c engine_s(4) Plastic strain rate [1/s]
c engine_s(5) damage
c engine_s(6:15) Engine state variables
c
c user_s(1:100) User state variable
c user_s(1:5) Tecplot output with "user"
c
c deps(3,3) Defomation tensor * dtime (strain increment)
c deps(1,1) = Dxx*dtime, deps(1,2)=Dxy*dtime, deps(2,2)=Dyy*dtime
c
c eps(3,3) Total Deformation gradient
c
c
c df(3,3) Relative Deformation gradient
c
c
c Material propeties are read from _wp.twm file
c
c Reserved material parameters
c d(5) lambda (Lame's constant)
c d(6) mu (Lame's constant)
c d(98) Conductivity
c d(99) Specific Heat * Density
c d(100) Density
c User defined material parameters
c d(25) UMATPAR01
c d(26) UMATPAR02
c d(27) UMATPAR03
c d(28) UMATPAR04
c d(29) UMATPAR05
c d(30) UMATPAR06
c d(31) UMATPAR07
c d(32) UMATPAR08
c d(33) UMATPAR09
c d(34) UMATPAR10
c d(35) UMATPAR11
c d(36) UMATPAR12
c d(37) UMATPAR13
c d(38) UMATPAR14
c d(39) UMATPAR15
c d(40) UMATPAR16
c d(41) UMATPAR17
c d(42) UMATPAR18
c d(43) UMATPAR19
c d(44) UMATPAR20
c d(45) UMATPAR21
c d(46) UMATPAR22
c d(47) UMATPAR23
c d(48) UMATPAR24
c d(49) UMATPAR25
c d(50) UMATPAR26
c d(51) UMATPAR27
c d(52) UMATPAR28
c d(53) UMATPAR29
c d(54) UMATPAR30
c d(55) UMATPAR31
c d(56) UMATPAR32
c d(57) UMATPAR33
c d(58) UMATPAR34
c d(59) UMATPAR35
c d(60) UMATPAR36
c d(61) UMATPAR37
c d(62) UMATPAR38
c d(63) UMATPAR39
c d(64) UMATPAR40
c d(65) UMATPAR41
c d(66) UMATPAR42
c d(67) UMATPAR43
c d(68) UMATPAR44
c d(69) UMATPAR45
c d(70) UMATPAR46
c d(71) UMATPAR47
c d(72) UMATPAR48
c d(73) UMATPAR49
c d(74) UMATPAR50
c
c Radial Return method elastic perfect-plastic material
c
real*8 dtime,temp,d(100),engine_s(15),user_s(100),eps(3,3),
& deps(3,3),sig(3,3)
real*8 sigtr(3,3),dsig(3,3),sigdiv(3,3),df(3,3)
real*8 df_(3,3)
real*8 bb(3,3),temp3x3(3,3)
real*8 b_trial(3,3),s_trial(3,3)
c
c
c
c Initialize state varianble
c
if ( time .le. dtime )then
engine_s(1:15)=0.0d0
eps(1:3,1:3)=0.0d0
eps(1,1)=1.0d0
eps(2,2)=1.0d0
eps(3,3)=1.0d0
user_s(1:100) = 0.0d0
end if
bb(1,1) = user_s(11) + 1.0d0
bb(2,1) = user_s(12)
bb(3,1) = user_s(13)
bb(1,2) = user_s(14)
bb(2,2) = user_s(15) + 1.0d0
bb(3,2) = user_s(16)
bb(1,3) = user_s(17)
bb(2,3) = user_s(18)
bb(3,3) = user_s(19) + 1.0d0
c Initialize plastic strain rate
engine_s(4)=0.0d0
c Initialize plastic work (heat generation)
engine_s(3)= 0.0d0
c
c Compute elastic Predictor
c
det = df(1,1)*df(2,2)*df(3,3) + df(2,1)*df(3,2)*df(1,3)
1 + df(1,2)*df(2,3)*df(3,1) - df(3,1)*df(2,2)*df(1,3)
2 - df(1,2)*df(2,1)*df(3,3) - df(3,2)*df(2,3)*df(1,1)
c
c
c
if ( det .le. 0.0d0 .or. isnan(det) ) then
c
c Failed element
c
user_s(1:15)=0.0d0
eps(1:3,1:3)=0.0d0
eps(1,1)=1.0d0
eps(2,2)=1.0d0
eps(3,3)=1.0d0
c
c sig(1:3,1:3) = 0.0d0
c
c p_ = dK*(det-1.0d0)
c sig(1,1)= sig(1,1) + p_
c sig(2,2)= sig(2,2) + p_
c sig(3,3)= sig(3,3) + p_
detF = eps(1,1)*eps(2,2)*eps(3,3) + eps(2,1)*eps(3,2)*eps(1,3)
1 + eps(1,2)*eps(2,3)*eps(3,1) - eps(3,1)*eps(2,2)*eps(1,3)
2 - eps(1,2)*eps(2,1)*eps(3,3) - eps(3,2)*eps(2,3)*eps(1,1)
write(6,*)"Negative det df",det,detF
c
return
end if
c
d_c = det**(-1.0/3.0)
c
df_(1:3,1:3)= d_c*df(1:3,1:3)
forall(i=1:3,j=1:3)
temp3x3(i,j)=df_(i,1)*bb(1,j)+df_(i,2)*bb(2,j)+
1 df_(i,3)*bb(3,j)
end forall
forall(i=1:3,j=1:3)
b_trial(i,j)=temp3x3(i,1)*df_(j,1)+temp3x3(i,2)*df_(j,2)+
1 temp3x3(i,3)*df_(j,3)
end forall
c
b_trace = (b_trial(1,1)+b_trial(2,2)+b_trial(3,3))
d_i_ = b_trace/3.0d0
s_trial(1:3,1:3) = b_trial(1:3,1:3)
s_trial(1,1) = s_trial(1,1) - d_i_ !b_trace/3.0
s_trial(2,2) = s_trial(2,2) - d_i_ !b_trace/3.0
s_trial(3,3) = s_trial(3,3) - d_i_ !b_trace/3.0
c
c
s_trial(1:3,1:3)= d(6)*s_trial(1:3,1:3)
yield_str = d(25) ! Initial yield stress
s_trial_norm = umat_div_norm(s_trial)
f_trial = s_trial_norm - sqrt(3.0/2.0)*yield_str
c
c Check for plastic loading
c
if ( f_trial .le. 0.0d0 ) then
sig(1:3,1:3) = s_trial(1:3,1:3)
bb(1:3,1:3) = b_trial(1:3,1:3)
else
c
c Return-mapping algorithm
c
dmu_ = d_i_*d(6)
c
delta_gm = f_trial/(2.0d0*dmu_)
c
c Plastic Strain
engine_s(1) = engine_s(1) + delta_gm
c
c Plastic Strain rate
engine_s(4) = delta_gm/dtime
c
c Plastic work rate (heat generation )
c
engine_s(3) = delta_gm*yield_str/dtime
c
temp3x3(1:3,1:3) = s_trial(1:3,1:3)/s_trial_norm
c Return map
sig(1:3,1:3)=
1 s_trial(1:3,1:3)-2.0d0*dmu_*delta_gm*temp3x3(1:3,1:3)
c Update of intermediate configuration
bb(1:3,1:3) = sig(1:3,1:3)/d(6)
bb(1,1) = d_i_ + bb(1,1) ! -1.0d0
bb(2,2) = d_i_ + bb(2,2) ! -1.0d0
bb(3,3) = d_i_ + bb(3,3) ! -1.0d0
end if
user_s(11) = bb(1,1) - 1.0d0
user_s(12) = bb(2,1)
user_s(13) = bb(3,1)
user_s(14) = bb(1,2)
user_s(15) = bb(2,2) - 1.0d0
user_s(16) = bb(3,2)
user_s(17) = bb(1,3)
user_s(18) = bb(2,3)
user_s(19) = bb(3,3) - 1.0d0
c Test user_s(1:5) output
user_s(1) = engine_s(1) ! plastic strain
c
c Bulk modulus
c
! dK = d(5) + 2.0d0/3.0d0*d(6)
dK = 5.0d-1 * d(5) + d(6) / 3.0d0
c
c
c Add elastic pressure term
c
det = eps(1,1)*eps(2,2)*eps(3,3) + eps(2,1)*eps(3,2)*eps(1,3)
1 + eps(1,2)*eps(2,3)*eps(3,1) - eps(3,1)*eps(2,2)*eps(1,3)
2 - eps(1,2)*eps(2,1)*eps(3,3) - eps(3,2)*eps(2,3)*eps(1,1)
c
p_ = dK * (det - 1.0d0 / det )
sig(1,1)= sig(1,1) + p_
sig(2,2)= sig(2,2) + p_
sig(3,3)= sig(3,3) + p_
END SUBROUTINE
c
c
double precision function umat_div_norm(sigdiv)
c
c Calculate diviatoric stress norm
c
real*8 sigdiv(3,3)
real*8 sigma_norm
c
sigma_norm = sigdiv(1,1)*sigdiv(1,1)
1 +sigdiv(1,2)*sigdiv(1,2)
1 +sigdiv(1,3)*sigdiv(1,3)
1 +sigdiv(2,1)*sigdiv(2,1)
1 +sigdiv(2,2)*sigdiv(2,2)
1 +sigdiv(2,3)*sigdiv(2,3)
1 +sigdiv(3,1)*sigdiv(3,1)
1 +sigdiv(3,2)*sigdiv(3,2)
1 +sigdiv(3,3)*sigdiv(3,3)
umat_div_norm = dsqrt(sigma_norm)
return
end
c
C
SUBROUTINE MAT_USER(d,time,dtime,temp,engine_s,user_s,
1 eps,deps,sig,df)
c
c AdvantEdge 2D/3D User defined material
c Total formulation example
c
c COPYRIGHT Third Wave Systems Inc. 2010
c www.thirdwavesys.com
c
#ifndef linux
!DEC$ ATTRIBUTES DLLEXPORT :: MAT_USER
#endif
c
c
implicit real*8 (a-h,o-z)
c
c
c Total formulation
c
c sig(3,3) Previous and updated stress [Pa]
c sig(1,1) = Sxx, sig(2,2)=Syy, sig(1,2)=Sxy
c dtime time step [s]
c temp temperature [c]
c
c engine_s(1) Plastic strain [1]
c engine_s(3) Plastic work rate (heat generation) [w]
c engine_s(4) Plastic strain rate [1/s]
c engine_s(5) damage
c engine_s(6:15) Engine state variables
c
c user_s(1:100) User state variable
c user_s(1:5) Tecplot output with "user"
c
c deps(3,3) Defomation tensor * dtime (strain increment)
c deps(1,1) = Dxx*dtime, deps(1,2)=Dxy*dtime, deps(2,2)=Dyy*dtime
c
c eps(3,3) Total Deformation gradient
c
c
c df(3,3) Relative Deformation gradient
c
c
c Material propeties are read from _wp.twm file
c
c Reserved material parameters
c d(5) lambda (Lame's constant)
c d(6) mu (Lame's constant)
c d(98) Conductivity
c d(99) Specific Heat * Density
c d(100) Density
c User defined material parameters
c d(25) UMATPAR01
c d(26) UMATPAR02
c d(27) UMATPAR03
c d(28) UMATPAR04
c d(29) UMATPAR05
c d(30) UMATPAR06
c d(31) UMATPAR07
c d(32) UMATPAR08
c d(33) UMATPAR09
c d(34) UMATPAR10
c d(35) UMATPAR11
c d(36) UMATPAR12
c d(37) UMATPAR13
c d(38) UMATPAR14
c d(39) UMATPAR15
c d(40) UMATPAR16
c d(41) UMATPAR17
c d(42) UMATPAR18
c d(43) UMATPAR19
c d(44) UMATPAR20
c d(45) UMATPAR21
c d(46) UMATPAR22
c d(47) UMATPAR23
c d(48) UMATPAR24
c d(49) UMATPAR25
c d(50) UMATPAR26
c d(51) UMATPAR27
c d(52) UMATPAR28
c d(53) UMATPAR29
c d(54) UMATPAR30
c d(55) UMATPAR31
c d(56) UMATPAR32
c d(57) UMATPAR33
c d(58) UMATPAR34
c d(59) UMATPAR35
c d(60) UMATPAR36
c d(61) UMATPAR37
c d(62) UMATPAR38
c d(63) UMATPAR39
c d(64) UMATPAR40
c d(65) UMATPAR41
c d(66) UMATPAR42
c d(67) UMATPAR43
c d(68) UMATPAR44
c d(69) UMATPAR45
c d(70) UMATPAR46
c d(71) UMATPAR47
c d(72) UMATPAR48
c d(73) UMATPAR49
c d(74) UMATPAR50
c
c Radial Return method elastic perfect-plastic material
c
real*8 dtime,temp,d(100),engine_s(15),user_s(100),eps(3,3),
& deps(3,3),sig(3,3)
real*8 sigtr(3,3),dsig(3,3),sigdiv(3,3),df(3,3)
real*8 df_(3,3)
real*8 bb(3,3),temp3x3(3,3)
real*8 b_trial(3,3),s_trial(3,3)
c
c
c
c Initialize state varianble
c
if ( time .le. dtime )then
engine_s(1:15)=0.0d0
eps(1:3,1:3)=0.0d0
eps(1,1)=1.0d0
eps(2,2)=1.0d0
eps(3,3)=1.0d0
user_s(1:100) = 0.0d0
end if
bb(1,1) = user_s(11) + 1.0d0
bb(2,1) = user_s(12)
bb(3,1) = user_s(13)
bb(1,2) = user_s(14)
bb(2,2) = user_s(15) + 1.0d0
bb(3,2) = user_s(16)
bb(1,3) = user_s(17)
bb(2,3) = user_s(18)
bb(3,3) = user_s(19) + 1.0d0
c Initialize plastic strain rate
engine_s(4)=0.0d0
c Initialize plastic work (heat generation)
engine_s(3)= 0.0d0
c
c Compute elastic Predictor
c
det = df(1,1)*df(2,2)*df(3,3) + df(2,1)*df(3,2)*df(1,3)
1 + df(1,2)*df(2,3)*df(3,1) - df(3,1)*df(2,2)*df(1,3)
2 - df(1,2)*df(2,1)*df(3,3) - df(3,2)*df(2,3)*df(1,1)
c
c
c
if ( det .le. 0.0d0 .or. isnan(det) ) then
c
c Failed element
c
user_s(1:15)=0.0d0
eps(1:3,1:3)=0.0d0
eps(1,1)=1.0d0
eps(2,2)=1.0d0
eps(3,3)=1.0d0
c
c sig(1:3,1:3) = 0.0d0
c
c p_ = dK*(det-1.0d0)
c sig(1,1)= sig(1,1) + p_
c sig(2,2)= sig(2,2) + p_
c sig(3,3)= sig(3,3) + p_
detF = eps(1,1)*eps(2,2)*eps(3,3) + eps(2,1)*eps(3,2)*eps(1,3)
1 + eps(1,2)*eps(2,3)*eps(3,1) - eps(3,1)*eps(2,2)*eps(1,3)
2 - eps(1,2)*eps(2,1)*eps(3,3) - eps(3,2)*eps(2,3)*eps(1,1)
write(6,*)"Negative det df",det,detF
c
return
end if
c
d_c = det**(-1.0/3.0)
c
df_(1:3,1:3)= d_c*df(1:3,1:3)
forall(i=1:3,j=1:3)
temp3x3(i,j)=df_(i,1)*bb(1,j)+df_(i,2)*bb(2,j)+
1 df_(i,3)*bb(3,j)
end forall
forall(i=1:3,j=1:3)
b_trial(i,j)=temp3x3(i,1)*df_(j,1)+temp3x3(i,2)*df_(j,2)+
1 temp3x3(i,3)*df_(j,3)
end forall
c
b_trace = (b_trial(1,1)+b_trial(2,2)+b_trial(3,3))
d_i_ = b_trace/3.0d0
s_trial(1:3,1:3) = b_trial(1:3,1:3)
s_trial(1,1) = s_trial(1,1) - d_i_ !b_trace/3.0
s_trial(2,2) = s_trial(2,2) - d_i_ !b_trace/3.0
s_trial(3,3) = s_trial(3,3) - d_i_ !b_trace/3.0
c
c
s_trial(1:3,1:3)= d(6)*s_trial(1:3,1:3)
yield_str = d(25) ! Initial yield stress
s_trial_norm = umat_div_norm(s_trial)
f_trial = s_trial_norm - sqrt(3.0/2.0)*yield_str
c
c Check for plastic loading
c
if ( f_trial .le. 0.0d0 ) then
sig(1:3,1:3) = s_trial(1:3,1:3)
bb(1:3,1:3) = b_trial(1:3,1:3)
else
c
c Return-mapping algorithm
c
dmu_ = d_i_*d(6)
c
delta_gm = f_trial/(2.0d0*dmu_)
c
c Plastic Strain
engine_s(1) = engine_s(1) + delta_gm
c
c Plastic Strain rate
engine_s(4) = delta_gm/dtime
c
c Plastic work rate (heat generation )
c
engine_s(3) = delta_gm*yield_str/dtime
c
temp3x3(1:3,1:3) = s_trial(1:3,1:3)/s_trial_norm
c Return map
sig(1:3,1:3)=
1 s_trial(1:3,1:3)-2.0d0*dmu_*delta_gm*temp3x3(1:3,1:3)
c Update of intermediate configuration
bb(1:3,1:3) = sig(1:3,1:3)/d(6)
bb(1,1) = d_i_ + bb(1,1) ! -1.0d0
bb(2,2) = d_i_ + bb(2,2) ! -1.0d0
bb(3,3) = d_i_ + bb(3,3) ! -1.0d0
end if
user_s(11) = bb(1,1) - 1.0d0
user_s(12) = bb(2,1)
user_s(13) = bb(3,1)
user_s(14) = bb(1,2)
user_s(15) = bb(2,2) - 1.0d0
user_s(16) = bb(3,2)
user_s(17) = bb(1,3)
user_s(18) = bb(2,3)
user_s(19) = bb(3,3) - 1.0d0
c Test user_s(1:5) output
user_s(1) = engine_s(1) ! plastic strain
c
c Bulk modulus
c
! dK = d(5) + 2.0d0/3.0d0*d(6)
dK = 5.0d-1 * d(5) + d(6) / 3.0d0
c
c
c Add elastic pressure term
c
det = eps(1,1)*eps(2,2)*eps(3,3) + eps(2,1)*eps(3,2)*eps(1,3)
1 + eps(1,2)*eps(2,3)*eps(3,1) - eps(3,1)*eps(2,2)*eps(1,3)
2 - eps(1,2)*eps(2,1)*eps(3,3) - eps(3,2)*eps(2,3)*eps(1,1)
c
p_ = dK * (det - 1.0d0 / det )
sig(1,1)= sig(1,1) + p_
sig(2,2)= sig(2,2) + p_
sig(3,3)= sig(3,3) + p_
END SUBROUTINE
c
c
double precision function umat_div_norm(sigdiv)
c
c Calculate diviatoric stress norm
c
real*8 sigdiv(3,3)
real*8 sigma_norm
c
sigma_norm = sigdiv(1,1)*sigdiv(1,1)
1 +sigdiv(1,2)*sigdiv(1,2)
1 +sigdiv(1,3)*sigdiv(1,3)
1 +sigdiv(2,1)*sigdiv(2,1)
1 +sigdiv(2,2)*sigdiv(2,2)
1 +sigdiv(2,3)*sigdiv(2,3)
1 +sigdiv(3,1)*sigdiv(3,1)
1 +sigdiv(3,2)*sigdiv(3,2)
1 +sigdiv(3,3)*sigdiv(3,3)
umat_div_norm = dsqrt(sigma_norm)
return
end
c