SUBROUTINE HYPELA(D,G,E,DE,S,TEMP,
1 DTEMP,NGENS,N,NN,KC,MATS,NDI,NSHEAR)
c implicit none
c* * * * * *
c
c user subroutine to define young's modulus and poisson's ratio
c as function of stress in non-linear elastic small strain
c material.
c
c d stress strain law to be formed by user
c g change in stress due to temperature effects
c e total strain
c de increment of strain
c s stress - should be updated by user
c temp state variables
c dtemp increment of state variables
c ngens size of stress - strain law
c n element number
c nn integration point number
c kc layer number
c mats material i.d.
c ndi number of direct components
c nshear number of shear components
c
c* * * * * *
implicit real*8 (a-h,o-z)
INCLUDE '../common/concom' ! 通过concom模块得到当前的计算步数
integer :: ngens,nn,kc,mats,ndi,nshear
real*8 :: e(1),de(1),temp(40),dtemp(40),g(1),d(ngens,ngens),s(1)
! temp(40)需要在前面initial condition 里面设置完成
integer :: n(2)

C* * * * * *C Local Arrays
C-------------------------------------
C EELAS - Elastic Strains ! 弹性应变
C EPLAS - Plastic Strains ! 塑性应变
C ALPHA - Shift Tensor ! 硬化参数
C Flow - Plasitc Flow Directions ! 塑性流动方向
C SIG - Stress at start of increment ! 增量步开始时的应力
C EPSPL - Plastic Strains at start of increment ! 增量步开始时的塑性应变
C
real*8 EELAS(ngens),EPLAS(ngens),ALPHA(ngens)
1 ,FLOW(ngens),SIG(ngens),EPSPL(ngens)
c integer :: inc, incsub, ncycle
real*8 E0,ENU,EBULK,EG,ELAM ! 初始弹性模量，泊松比，体积模量，剪切模量
integer :: K1, K2 ! 循环变量
real*8 :: Smises, SYIELD, hard ! VM应力, 屈服应力，硬化模量
real*8 :: SIGM,DEQPL ! 平均正应力, 等效应变增量
real*8 :: effg, efflam, effhard ! 等效剪切模量，等效lame常数，等效硬化模量
C
real*8,parameter :: ENUMAX=.49999D0,TOLER=1.D-6

C
C temp (1) - E ! 弹性模量
C temp (2) - NU ! 泊松比
C temp (3) - SYIELD ! 屈服强度
C temp (4) - HARD ! 硬化模量 C ------------------------------
C 初始状态赋值
C
if(inc==0.and.incsub==0.and.ncycle==0) then
TEMP(1)=200e3; ! 弹性模量 200GPa
TEMP(2)=0.27; ! 泊松比为 0.27
TEMP(3)=210.0; ! 屈服强度 210MPa
TEMP(4)=20e3; ! 硬化模量 20GPa
end if

E0=TEMP(1)
ENU=Min(TEMP(2),ENUMAX)
EBULK=E0/(1.-2.*ENU)/3.
EG=E0/(1.+ENU)/2.
ELAM=(EBULK*3.-EG*2.)/3.

C 弹性矩阵

D=0.
do K1=1, NDI
do K2=1, NDI
D(K2,K1)=ELAM
end do
D(K1,K1)=EG*2.+ELAM
END do
do K1=NDI+1, NGENS
D(K1,K1)=EG
end do

C 数组赋值

EELAS=TEMP(4+1:4+NGENS)
EPLAS=TEMP(4+NGENS+1:4+2*NGENS)
ALPHA=TEMP(4+2*NGENS+1:4+3*NGENS)

C
C 保存初始应力并计算试算应力
C
do k1=1, NGENS
SIG(K1)=S(K1)
EPSPL(K1)=EPLAS(K1)
EELAS(K1)=EELAS(K1)+DE(K1)
do k2=1,NGENS
s(k2)=s(k2)+d(k2,k1)*de(k1)
end do
end do

C
C 计算von mises 应力
C
smises=(s(1)-alpha(1)-s(2)+alpha(2))**2
1 +(s(2)-alpha(2)-s(3)+alpha(3))**2
2 +(s(3)-alpha(3)-s(1)+alpha(1))**2
do K1=NDI+1,NGENS
smises=smises+6.*(s(k1)-alpha(k1))**2
end do
smises=sqrt(smises/2.)

C
C 得到屈服应力和硬化模量
C
SYIELD=temp(3)
hard=temp(4)

C
C 判断是否屈服
C
if(smises.gt.(1+toler)*syield) then
C
C 如果屈服，则计算硬化方向 dF/dSIG
C
SIGM=(s(1)+s(2)+s(3))/3.
do K1=1, NDI
flow(k1)=(s(k1)-alpha(k1)-SIGM)/smises
end do
do k1=NDI+1, NGENS
Flow(k1)=(s(k1)-alpha(k1))/smises
end do
C
C 计算等效应变增量
C
deqpl=(smises-syield)/(EG*3.+HARD)

C
C Update shift tensor, elastic an plastic strains and stress
C
do k1=1, ndi
alpha(k1)=alpha(k1)+hard*flow(k1)*deqpl
eplas(k1)=eplas(k1)+3./2.*flow(k1)*deqpl
eelas(k1)=eelas(k1)-3./2.*flow(k1)*deqpl
s(k1)=alpha(k1)+flow(k1)*syield+SIGM
end do
do k1=ndi+1, ngens
alpha(k1)=alpha(k1)+hard*flow(k1)*deqpl
eplas(k1)=eplas(k1)+3.*flow(k1)*deqpl
eelas(k1)=eelas(k1)-3.*flow(k1)*deqpl
s(k1)=alpha(k1)+flow(k1)*syield
end do

C
C 得到切线刚度矩阵
C
effg=eg*(syield+hard*deqpl)/smises
efflam=(ebulk*3.-effg*2.)/3.
effhard=eg*3.*hard/(eg*3.+hard)-effg*3.
do k1=1, ndi
do k2=1, ndi
d(k2,k1)=efflam
end do
d(k1,k1)=effg*2.+efflam
end do
do k1=ndi+1,ngens
d(k1,k1)=effg
end do
do k1=1,ngens
do k2=1, ngens
d(k2,k1)=d(k2,k1)+effhard*flow(k2)*flow(k1)
end do
end do
end if
C
C 保存当前弹性应变，塑性应变和硬化参数
C
do k1=1, ngens
temp(4+k1)=eelas(k1)
temp(4+k1+ngens)=eplas(k1)
temp(4+k1+2*ngens)=alpha(k1)
end do
return
end