C =====================================================
C 清华大学研究生精品课程《钢筋混凝土有限元》教学程序
C 等强双线性硬化本构模型
C Isotropic harding plastic subroutine for MSC.MARC
C 陆新征 江见鲸
C 清华大学土木工程系，北京，100084
C Lu Xinzheng Jiang Jianjing, Dept. Civil Engrg. of Tsinghua University
C last revised: Sep. 2003.
C 主要参考文献：
C 1. MARC Volumn D user subroutine and Special Routines"
C 2. 江见鲸 "钢筋混凝土结构非线性有限元分析"
C
C =====================================================
C 子程序说明：
C zero: 得到等效应力
C yiel: 得到当前屈服应力
C wkslp: 得到强化模量
C assoc: 得到流动方向

subroutine assoc(stot,sinc,sc,t,ngens,ndi,nshear,n,nnn,layer)
c* * * * * *
c
c provides flow rule for general plasticity option.
C 流动法则
c
c stoc stress array !应力矩阵
c sinc flow direction array !流动向量
c sc hydrostatic stress !静水压力
c t not used
c ngens number of stress components
c ndi (3) number of direct stress components
c nshear (3) number of shear stress components
c n element number
c nnn integration point number
c layer layer number
c
c* * * * * *
implicit real*8 (a-h,o-z) dp
dimension stot(*),sinc(*)
do 26 i=1,ndi
26 sinc(i)=0.5*(3.*stot(i)-sc)
if(nshear.eq.0)go to 30
f3=3.
do 28 i=1,nshear
j=ndi+i
28 sinc(j)=f3*stot(j)
30 continue
return
end

subroutine wkslp(m,nn,kc,mats,slope,ebarp,eqrate,stryt,dt,ifirst)
implicit real*8 (a-h,o-z) dp
c* * * * * *
c
c user subroutine to define work hardening slope.
c
c 硬化子程序
c m is the current user element number
c nn is the integration point number
c kc is the layer number (=1 for non-layered elements)
c mats is the current material id
c slope work hardening slope to be defined in this routine
c ebarp is the total equivalent plastic strain
c eqrate is the equivalent plastic strain rate
c stryt current yield stress that optionally can be defined in this routine
c dt is the current total temperature
c ifirst flag distinguishing tenth cycle properties for
c ornl option
c
c the internal element number mint can be obtained with
c mint=ielint(m)
c
c* * * * * *
SLOPE=20e3 !强化模量为Es=2GPa
STRYT=210+ebarp*20e3 !强化后的屈服应力为210MPa+eps_eq_pl*Es
return
end

real*8 function yiel(m,nn,kc,yld,ifirst,dt,eplas,erate,matz,
* jprops)
c* * * * * *
C 得到屈服应力
c
c find current value of yield stress
c
c n element index
c yld initial yield stress
c ifirst flag for 10th cycle yield
c dt temperature
c eplas equivalent plastic strain
c erate equivalent plastic strain rate
c jprops table asscociated with the yield
c
c* * * * * *
implicit real*8 (a-h,o-z)
c* * * * * *
c
c
yiel=yld
c
c
c user subroutine wkslp
c 调用做功强化子程序得到强化模量和屈服应力
c
stryt=0.d0
call wkslp(m,nn,kc,mats,slope,eplas,erate,stryt,dt,ifirst)
yiel=stryt

return
end

real*8 function zero(ndi,nshear,t,iort,ianiso,yrdir,yrshr,amm,a0)
c* * * * * *
c
c 得到等效应力
c find equivalent stress for a stress tensor.
c
c ndi number of direct components of stress
c nshear number of shear components of stress
c t(i) component i of stress
c iort indicating if curvilinear coordinates are used
c iort = 0 - no curvilinear coordinates are used
c iort = 1 - get mixed components in t
c iort = 2 - assumes mixed compcomponents already in t
c ianiso flag indicating if anisotropy is used
c yrdir direct components for hill's anisotropic plasticity
c yrshr shear components for hill's anisotropic plasticity
c amm the metric if curvilinear coordinates are used
c a0 the metric scale factor if curvilinear coordinates are used
c
c* * * * * *
implicit real*8 (a-h,o-z)
dimension t(*)
dimension s(6),yrdir(3),yrshr(3),amm(3)
zero=0.d0
c
c 等效应力即为von mises 应力
c
if(ndi.eq.1) then
zero=t(1)*t(1)*2.d0
else if(ndi.eq.2) then
zero=2.d0*(t(1)*t(1)+t(2)*t(2)-t(1)*t(2))
else if(ndi.eq.3) then
zero=t(1)*t(1)+t(2)*t(2)+t(3)*t(3)
* -t(1)*t(2)-t(2)*t(3)-t(3)*t(1)
zero=2.d0*zero
endif
if(nshear.eq.1) then
zero=zero+6.d0*t(ndi+1)*t(ndi+1)
else if(nshear.eq.2) then
zero=zero+6.d0*(t(ndi+1)*t(ndi+1)+t(ndi+2)*t(ndi+2))
else if(nshear.eq.3) then
zero=zero+6.d0*(t(ndi+1)*t(ndi+1)
* +t(ndi+2)*t(ndi+2)+t(ndi+3)*t(ndi+3))
endif
c
if(zero.lt.0.d0) zero=0.d0
zero=sqrt(0.5d0*zero)
c
return
end