适合HCP mg合金的黄永刚umat的晶体塑性
subroutine umat(stress,statev,ddsdde,sse,spd,scd,rpl,
1 ddsddt,drplde,drpldt,stran,dstran,time,dtime,
2 temp,dtemp,predef,dpred,cmname,ndi,nshr,ntens,
3 nstatv,props,nprops,coords,drot,pnewdt,celent,
4 dfgrd0,dfgrd1,noel,npt,layer,kspt,kstep,kinc)
c
include 'aba_param.inc'
c modified from huang umat- jg:20/07/12
c *******************************************************************
c - only suitable for finite deformation isotropic elasticity
c with fcc and hcp crystal slip.
c - most parameters are now hard-coded.
c - most documentation from original umat removed.
c - constants converted to double precision
c - direction vectors must be orthogonal unit vectors
c *******************************************************************
parameter (zero=0.d0, one=1.d0, two=2.d0)
character*8 cmname
dimension stress(ntens),statev(nstatv),ddsdde(ntens,ntens),
1 dstran(ntens),props(nprops),drot(3,3)
dimension ispdir(3),ispnor(3),slpdir(3,18),slpnor(3,18),
1 slpdef(6,18),slpspn(3,18),dspdir(3,18),dspnor(3,18),
2 d(6,6),rotate(3,3),fslip(18),dfdxsp(18),ddemsd(6,18),
3 h(18,18),ddgdde(18,6),dstres(6),delats(6),dspin(3),
4 dvgrad(3,3),dgamma(18),dtausp(18),dgslip(18),
5 workst(18,18),indx(18),term(3,3),trm0(3,3),itrm(3)
dimension fslip1(18),stres1(6),gamma1(18),tausp1(18),gslp1(18),
1 spnor1(3,18),spdir1(3,18),ddsde1(6,6),dsold(6),
2 dgamod(18),dtauod(18),dgspod(18),dspnro(3,18),
3 dspdro(3,18),dhdgdg(18,18),rwkdir(3,24),rwknor(3,24),
3 indxL(3),termd(3),termn(3),gamma(18),tauslpL(18)
C----- Elastic matrix in GLOBAL system
gshear = props(1)/(2.*(1.+props(2)))
e11 = 2.*gshear*(one-props(2))/(1.-2.*props(2))
e12 = 2.*gshear*props(2)/(1.-2.*props(2))
d = 0.
do j = 1,3
d(j,j) = e11
do i = 1,3
if(i.ne.j) d(i,j) = e12
end do
d(j+3,j+3) = gshear
end do
c------ Crystal Type:
ictype=nint(props(9))
if(ictype == 1)then
c FCC
nslptl = 12
elseif(ictype == 2)then
C HCP - Basal and Prism
nslptl = 6
elseif(ictype == 3)then
c HCP - Basal, Prism, Pyr
nslptl = 12
else
c HCP - Basal, Prism, Pyr, Twin
nslptl = 18
endif
adot=0.001
C----- Implicit integration parameter: THETA
theta = 0.5
C----- Iteration
itrmax = 1
gamerr = 1.e-5
nitrtn = -1
c
dsold = zero
dgamod = zero
dtauod = zero
dgspod = zero
dspnro = zero
dspdro = zero
C----- Increment of spin associated with the material element: DSPIN
do j = 1,3
do i = 1,3
term(i,j) = drot(j,i)
trm0(i,j) = drot(j,i)
end do
term(j,j) = term(j,j)+one
trm0(j,j) = trm0(j,j)-one
end do
call ludcmp(term, 3, 3, itrm, ddcmp)
do j = 1,3
call lubksb(term, 3, 3, itrm, trm0(1,j))
end do
dspin(1) = trm0(2,1)-trm0(1,2)
dspin(2) = trm0(1,3)-trm0(3,1)
dspin(3) = trm0(3,2)-trm0(2,3)
C----- Increment of dilatational strain: DEV
dev = zero
do i = 1,ndi
dev = dev+dstran(i)
end do
C----- Iteration starts
1000 continue
nitrtn = nitrtn+1
C----- Check whether the current stress state is the initial state
if (statev(1) == 0.) then
if (ictype == 1)then
c----- generating all possible slip directions for fcc
rwkdir = 1./sqrt(2.)
do j = 1,6
do i = 1,3
if (i.eq.j.or.j-i.eq.3)rwkdir(i,j) = 0.
end do
end do
rwkdir(1,6) = -rwkdir(1,6)
rwkdir(2,4) = -rwkdir(2,4)
rwkdir(3,5) = -rwkdir(3,5)
c----- generating all possible slip planes for fcc
rwknor = 1./sqrt(3.)
do i = 1,3
do j = 1,4
if (j.eq.i+1)rwknor(i,j) = -rwknor(i,j)
end do
end do
c------ Generating all slip systems for FCC
nslip = 0
do j = 1,4
do i = 1,6
rdot = 0.
do k = 1,3
rdot = rdot+rwkdir(k,i)*rwknor(k,j)
end do
if (rdot.eq.0.) then
nslip = nslip+1
do k = 1,3
slpdir(k,nslip) = rwkdir(k,i)
slpnor(k,nslip) = rwknor(k,j)
end do
end if
end do
end do
else
c----- generating slip directions and normals for hcp-basal
rwkdir = 0.
rwknor = 0.
angle = acos(-1.)/3.
cangle = cos(angle)
sangle = sin(angle)
rwkdir(1,1) = 1.
rwkdir(2,1) = 0.
rwkdir(1,2) = cangle
rwkdir(2,2) = sangle
rwkdir(1,3) = -cangle
rwkdir(2,3) = sangle
rwknor(3,1) = 1.
rwknor(3,2) = 1.
rwknor(3,3) = 1.
do i = 1,3
do k = 1,3
slpdir(k,i) = rwkdir(k,i)
slpnor(k,i) = rwknor(k,i)
enddo
enddo
c----- generating slip directions and normals for hcp-prismatic
rwknor = 0.
rwknor(1,1) = 0.
rwknor(2,1) = -1.
rwknor(1,2) = sangle
rwknor(2,2) = -cangle
rwknor(1,3) = -sangle
rwknor(2,3) = -cangle
do i = 4,6
do k = 1,3
slpdir(k,i) = rwkdir(k,i-3)
slpnor(k,i) = rwknor(k,i-3)
enddo
enddo
endif
if(ictype >= 3)then
c ##### 2nd order pyramidal <a+c> #####
aspect = 1.624
c slip directions
do j = 1,6
rwkdir(3,j) = aspect
enddo
rwkdir(1,1) = -cangle
rwkdir(2,1) = -sangle
rwkdir(1,2) = cangle
rwkdir(2,2) = -sangle
rwkdir(1,3) = -2.*cangle
rwkdir(2,3) = 0.
do j = 4,6
rwkdir(1,j) = -rwkdir(1,j-3)
rwkdir(2,j) = -rwkdir(2,j-3)
enddo
rlength=sqrt(1.+aspect*aspect)
do j = 1,6
do i = 1,3
rwkdir(i,j) = rwkdir(i,j)/rlength
enddo
enddo
c slip normals
do j = 1,6
rwknor(3,j) = 4.*sangle*cangle
enddo
rwknor(1,1) = aspect*sangle
rwknor(2,1) = 3.*aspect*cangle
rwknor(1,2) = -aspect*sangle
rwknor(2,2) = 3.*aspect*cangle
rwknor(1,3) = 2.*aspect*sangle
rwknor(2,3) = 0.
do j = 4,6
rwknor(1,j) = -rwknor(1,j-3)
rwknor(2,j) = -rwknor(2,j-3)
enddo
rlength=sqrt(3.*(1.+aspect*aspect))
do j = 1,6
do i = 1,3
rwknor(i,j) = rwknor(i,j)/rlength
enddo
enddo
nslip = 6
do j = 1,6
nslip = nslip+1
do i = 1,3
slpdir(i,nslip) = rwkdir(i,j)
slpnor(i,nslip) = rwknor(i,j)
enddo
enddo
if(ictype == 4)then
c ###### twinning planes #####
c slip directions
do j = 1,6
rwkdir(3,j) = aspect
enddo
rwkdir(1,1) = 0.
rwkdir(2,1) = -2.*sangle
rwkdir(1,2) = -3.*cangle
rwkdir(2,2) = -1.*sangle
rwkdir(1,3) = -3.*cangle
rwkdir(2,3) = 1.*sangle
do j = 4,6
rwkdir(1,j) = -rwkdir(1,j-3)
rwkdir(2,j) = -rwkdir(2,j-3)
enddo
rlength=sqrt(3.+aspect*aspect)
do j = 1,6
do i = 1,3
rwkdir(i,j) = rwkdir(i,j)/rlength
enddo
enddo
c slip normals
do j = 1,6
rwknor(3,j) = 4.*sangle*cangle
enddo
rwknor(1,1) = 0.
rwknor(2,1) = 2.*aspect*cangle
rwknor(1,2) = aspect*sangle
rwknor(2,2) = aspect*cangle
rwknor(1,3) = aspect*sangle
rwknor(2,3) = -aspect*cangle
do j = 4,6
rwknor(1,j) = -rwknor(1,j-3)
rwknor(2,j) = -rwknor(2,j-3)
enddo
do j = 1,6
do i = 1,3
rwknor(i,j) = rwknor(i,j)/rlength
enddo
enddo
do j = 1,6
nslip = nslip+1
do i = 1,3
slpdir(i,nslip) = rwkdir(i,j)
slpnor(i,nslip) = rwknor(i,j)
enddo
enddo
endif
endif
C----- Unit vectors in slip dirs and unit norms-Global system
c----- Generate rotation matrix
do i = 1,3
term(i,1) = props(i+2)
term(i,2) = props(i+5)
enddo
term(1,3) = term(2,1)*term(3,2)-term(3,1)*term(2,2)
term(2,3) = term(3,1)*term(1,2)-term(1,1)*term(3,2)
term(3,3) = term(1,1)*term(2,2)-term(2,1)*term(1,2)
call ludcmp (term, 3, 3, indxL, dcmp)
rotate = 0.
do j = 1,3
do i = 1,3
if (i.eq.j)rotate(i,j) = 1.
end do
end do
do j = 1,3
call lubksb (term, 3, 3, indxL, rotate(1,j))
end do
c--- Rotate slip normals and directions to global system
do j = 1,nslptl
do i = 1,3
termd(i) = 0.
termn(i) = 0.
do k = 1,3
termd(i) = termd(i)+rotate(i,k)*slpdir(k,j)
termn(i) = termn(i)+rotate(i,k)*slpnor(k,j)
end do
end do
do i = 1,3
slpdir(i,j) = termd(i)
slpnor(i,j) = termn(i)
end do
end do
C----- Get Slip deformation tensor: SLPDEF (Schmid factors)
do j=1,nslptl
slpdef(1,j)=slpdir(1,j)*slpnor(1,j)
slpdef(2,j)=slpdir(2,j)*slpnor(2,j)
slpdef(3,j)=slpdir(3,j)*slpnor(3,j)
slpdef(4,j)=slpdir(1,j)*slpnor(2,j)
1 +slpdir(2,j)*slpnor(1,j)
slpdef(5,j)=slpdir(1,j)*slpnor(3,j)
1 +slpdir(3,j)*slpnor(1,j)
slpdef(6,j)=slpdir(2,j)*slpnor(3,j)
1 +slpdir(3,j)*slpnor(2,j)
end do
C----- Store normals and directions
statev(nstatv)=nslptl
idnor=3*nslptl
iddir=6*nslptl
do j=1,nslptl
do i=1,3
idnor=idnor+1
iddir=iddir+1
statev(idnor)=slpnor(i,j)
statev(iddir)=slpdir(i,j)
end do
end do
C----- Initial value of the current strength for all slip systems
do j=1,nslptl
if(ictype == 1)then
statev(j)=props(10)
else
if(j<=3)then
statev(j)=props(10)
elseif(j<=6)then
statev(j)=props(13)
elseif(j<=12)then
statev(j)=props(16)
else
statev(j)=props(19)
endif
endif
enddo
C----- Initial value of shear strain in slip systems
do i=1,nslptl
statev(nslptl+i)=0.
end do
statev(9*nslptl+1)=0.
C----- Initial value of the resolved shear stress in slip systems
do i=1,nslptl
term1=0.
do j=1,ntens
if (j.le.ndi) then
term1=term1+slpdef(j,i)*stress(j)
else
term1=term1+slpdef(j-ndi+3,i)*stress(j)
end if
end do
statev(2*nslptl+i)=term1
end do
else
C----- Current stress state
idnor=3*nslptl
iddir=6*nslptl
do j=1,nslptl
do i=1,3
idnor=idnor+1
iddir=iddir+1
slpnor(i,j)=statev(idnor)
slpdir(i,j)=statev(iddir)
end do
end do
C----- Slip deformation tensor: SLPDEF (Schmid factors)
do j=1,nslptl
slpdef(1,j)=slpdir(1,j)*slpnor(1,j)
slpdef(2,j)=slpdir(2,j)*slpnor(2,j)
slpdef(3,j)=slpdir(3,j)*slpnor(3,j)
slpdef(4,j)=slpdir(1,j)*slpnor(2,j)
1 +slpdir(2,j)*slpnor(1,j)
slpdef(5,j)=slpdir(1,j)*slpnor(3,j)
1 +slpdir(3,j)*slpnor(1,j)
slpdef(6,j)=slpdir(2,j)*slpnor(3,j)
1 +slpdir(3,j)*slpnor(2,j)
end do
end if
C----- Slip spin tensor: SLPSPN
do j=1,nslptl
slpspn(1,j)=0.5*(slpdir(1,j)*slpnor(2,j)-
2 slpdir(2,j)*slpnor(1,j))
slpspn(2,j)=0.5*(slpdir(3,j)*slpnor(1,j)-
2 slpdir(1,j)*slpnor(3,j))
slpspn(3,j)=0.5*(slpdir(2,j)*slpnor(3,j)-
2 slpdir(3,j)*slpnor(2,j))
end do
C----- Double dot product of elastic moduli tensor with the slip
C deformation tensor
do j=1,nslptl
do i=1,6
ddemsd(i,j)=0.
do k=1,6
ddemsd(i,j)=ddemsd(i,j)+d(k,i)*slpdef(k,j)
end do
end do
end do
do j=1,nslptl
ddemsd(4,j)=ddemsd(4,j)-slpspn(1,j)*stress(1)
ddemsd(5,j)=ddemsd(5,j)+slpspn(2,j)*stress(1)
if (ndi.gt.1) then
ddemsd(4,j)=ddemsd(4,j)+slpspn(1,j)*stress(2)
ddemsd(6,j)=ddemsd(6,j)-slpspn(3,j)*stress(2)
end if
if (ndi.gt.2) then
ddemsd(5,j)=ddemsd(5,j)-slpspn(2,j)*stress(3)
ddemsd(6,j)=ddemsd(6,j)+slpspn(3,j)*stress(3)
end if
if (nshr.ge.1) then
ddemsd(1,j)=ddemsd(1,j)+slpspn(1,j)*stress(ndi+1)
ddemsd(2,j)=ddemsd(2,j)-slpspn(1,j)*stress(ndi+1)
ddemsd(5,j)=ddemsd(5,j)-slpspn(3,j)*stress(ndi+1)
ddemsd(6,j)=ddemsd(6,j)+slpspn(2,j)*stress(ndi+1)
end if
if (nshr.ge.2) then
ddemsd(1,j)=ddemsd(1,j)-slpspn(2,j)*stress(ndi+2)
ddemsd(3,j)=ddemsd(3,j)+slpspn(2,j)*stress(ndi+2)
ddemsd(4,j)=ddemsd(4,j)+slpspn(3,j)*stress(ndi+2)
ddemsd(6,j)=ddemsd(6,j)-slpspn(1,j)*stress(ndi+2)
end if
if (nshr.eq.3) then
ddemsd(2,j)=ddemsd(2,j)+slpspn(3,j)*stress(ndi+3)
ddemsd(3,j)=ddemsd(3,j)-slpspn(3,j)*stress(ndi+3)
ddemsd(4,j)=ddemsd(4,j)-slpspn(2,j)*stress(ndi+3)
ddemsd(5,j)=ddemsd(5,j)+slpspn(1,j)*stress(ndi+3)
end if
end do
C----- Shear strain-rate in a slip system at the start of increment:
do i=1,nslptl
tauslp=statev(2*nslptl+i)
if(i>=13.and.tauslp<=0)then
yield=1.e6
else
yield=statev(i)
endif
x=tauslp/yield
fslip(i)=adot*((abs(x))**50.)*dsign(1.d0,x)
dfdxsp(i)=50.*adot*(abs(x))**(50.-1.0)
end do
C----- Self- and latent-hardening
do i=1,nslptl
gamma(i)=statev(nslptl+1)
enddo
gamtol=statev(9*nslptl+1)
do iself = 1,nslptl
do latent = 1,nslptl
if(ictype == 1)then
c FCC
term1 = props(12)*gamtol/(props(11)-props(10))
term2 = 2.*exp(-term1)/
* (1.+exp(-2.*term1))
hlatnt = props(12)*term2**2
else
C BASAL
if(iself <= 3)then
if(latent <= 3)then
q = 0.2
else
q = 0.5
endif
if(iself == latent)q = 1.
hlatnt = q*props(12)
C PRISM
elseif(iself <= 6)then
if(latent <= 12)then
q = 0.2
else
q = 0.5
endif
if(iself == latent)q = 1.
hlatnt = q*props(15)*(1.d0-(props(13)/
* props(14)))*exp(-props(15)*(gamtol/
* props(14)))
C PYRM
elseif(iself <= 12)then
if(latent <= 6)then
q = 1.
elseif(latent <= 12)then
q = 0.2
else
q = 0.25
endif
if(iself == latent)q = 1.
hlatnt = q*props(18)*(1.d0-props(16)/
* props(17))*exp(-props(18)*gamtol/
* props(17))
C TWIN
else
if(latent <= 6)then
q = 1.
else
q = 0.2
endif
if(iself == latent)q = 1.
if(gamtol <= props(21))then
hlatnt = q*props(20)
else
hlatnt = q*props(20)*((gamtol/props(21))
* **(props(22)-1.))
endif
endif
endif
h(iself,latent) = hlatnt
enddo
end do
C----- Solve the increment of shear strain in a slip system
term1=theta*dtime
do i=1,nslptl
tauslp=statev(2*nslptl+i)
if(i>=13.and.tauslp<=0)then
yield=1.e6
else
yield=statev(i)
endif
gslip=yield
x=tauslp/gslip
term2=term1*dfdxsp(i)/gslip
term3=term1*x*dfdxsp(i)/gslip
do j=1,nslptl
term4=0.
do k=1,6
term4=term4+ddemsd(k,i)*slpdef(k,j)
end do
workst(i,j)=term2*term4+h(i,j)*term3*dsign(1.d0,fslip(j))
if(nitrtn.gt.0)workst(i,j)=workst(i,j)+term3*dhdgdg(i,j)
end do
workst(i,i)=workst(i,i)+1.
end do
call ludcmp (workst, nslptl, 18, indx, ddcmp)
c----- increment of shear strain in a slip system: dgamma
term1=theta*dtime
do i=1,nslptl
if (nitrtn.eq.0) then
tauslp=statev(2*nslptl+i)
if(i>=13.and.tauslp<=0)then
yield=1.e6
else
yield=statev(i)
endif
gslip=yield
x=tauslp/gslip
term2=term1*dfdxsp(i)/gslip
dgamma(i)=0.
do j=1,ndi
dgamma(i)=dgamma(i)+ddemsd(j,i)*dstran(j)
end do
if (nshr.gt.0) then
do j=1,nshr
dgamma(i)=dgamma(i)+ddemsd(j+3,i)*dstran(j+ndi)
end do
end if
dgamma(i)=dgamma(i)*term2+fslip(i)*dtime
else
dgamma(i)=term1*(fslip(i)-fslip1(i))+fslip1(i)*dtime
2 -dgamod(i)
end if
end do
call lubksb (workst, nslptl, 18, indx, dgamma)
do i=1,nslptl
dgamma(i)=dgamma(i)+dgamod(i)
end do
c----- update the shear strain in a slip system:
do i=1,nslptl
statev(nslptl+i)=statev(nslptl+i)+dgamma(i)-dgamod(i)
end do
C----- Increment of current strength in a slip system: DGSLIP
do i=1,nslptl
dgslip(i)=0.
do j=1,nslptl
dgslip(i)=dgslip(i)+h(i,j)*abs(dgamma(j))
end do
end do
C----- Update the current strength in a slip system:
do i=1,nslptl
statev(i)=statev(i)+dgslip(i)-dgspod(i)
end do
C----- Increment of strain associated with lattice stretching: DELATS
do j=1,6
delats(j)=0.
end do
do j=1,3
if (j.le.ndi) delats(j)=dstran(j)
do i=1,nslptl
delats(j)=delats(j)-slpdef(j,i)*dgamma(i)
end do
end do
do j=1,3
if (j.le.nshr) delats(j+3)=dstran(j+ndi)
do i=1,nslptl
delats(j+3)=delats(j+3)-slpdef(j+3,i)*dgamma(i)
end do
end do
C----- Increment of deformation gradient associated with lattice stretching
do j=1,3
do i=1,3
if (i.eq.j) then
dvgrad(i,j)=delats(i)
else
dvgrad(i,j)=delats(i+j+1)
end if
end do
end do
do j=1,3
do i=1,j
if (j.gt.i) then
ij2=i+j-2
if (mod(ij2,2).eq.1) then
term1=1.
else
term1=-1.
end if
dvgrad(i,j)=dvgrad(i,j)+term1*dspin(ij2)
dvgrad(j,i)=dvgrad(j,i)-term1*dspin(ij2)
do k=1,nslptl
dvgrad(i,j)=dvgrad(i,j)-term1*dgamma(k)*
2 slpspn(ij2,k)
dvgrad(j,i)=dvgrad(j,i)+term1*dgamma(k)*
2 slpspn(ij2,k)
end do
end if
end do
end do
C----- Increment of resolved shear stress in a slip system: DTAUSP
do i=1,nslptl
dtausp(i)=0.
do j=1,6
dtausp(i)=dtausp(i)+ddemsd(j,i)*delats(j)
end do
end do
C----- Update the resolved shear stress in a slip system:
do i=1,nslptl
statev(2*nslptl+i)=statev(2*nslptl+i)+dtausp(i)-dtauod(i)
end do
C----- Increment of stress: DSTRES
do i=1,ntens
dstres(i)=-stress(i)*dev
end do
do i=1,ndi
do j=1,ndi
dstres(i)=dstres(i)+d(i,j)*dstran(j)
end do
if (nshr.gt.0) then
do j=1,nshr
dstres(i)=dstres(i)+d(i,j+3)*dstran(j+ndi)
end do
end if
do j=1,nslptl
dstres(i)=dstres(i)-ddemsd(i,j)*dgamma(j)
end do
end do
if (nshr.gt.0) then
do i=1,nshr
do j=1,ndi
dstres(i+ndi)=dstres(i+ndi)+d(i+3,j)*dstran(j)
end do
do j=1,nshr
dstres(i+ndi)=dstres(i+ndi)+d(i+3,j+3)*dstran(j+ndi)
end do
do j=1,nslptl
dstres(i+ndi)=dstres(i+ndi)-ddemsd(i+3,j)*dgamma(j)
end do
end do
end if
C----- Update the stress: STRESS
do i=1,ntens
stress(i)=stress(i)+dstres(i)-dsold(i)
end do
C----- Increment of normal to a slip plane and a slip direction
do j=1,nslptl
do i=1,3
dspnor(i,j)=0.
dspdir(i,j)=0.
do k=1,3
dspnor(i,j)=dspnor(i,j)-slpnor(k,j)*dvgrad(k,i)
dspdir(i,j)=dspdir(i,j)+slpdir(k,j)*dvgrad(i,k)
end do
end do
end do
C----- Update the normal to a slip plane and a slip direction
idnor=3*nslptl
iddir=6*nslptl
do j=1,nslptl
do i=1,3
idnor=idnor+1
statev(idnor)=statev(idnor)+dspnor(i,j)-dspnro(i,j)
iddir=iddir+1
statev(iddir)=statev(iddir)+dspdir(i,j)-dspdro(i,j)
end do
end do
C----- Derivative of shear strain increment in a slip system w.r.t.
C strain increment: DDGDDE
term1=theta*dtime
do i=1,ntens
do j=1,nslptl
tauslp=statev(2*nslptl+j)
if(j>=13.and.tauslp<=0)then
yield=1.e6
else
yield=statev(j)
endif
gslip=yield
x=tauslp/gslip
term2=term1*dfdxsp(j)/gslip
if (i.le.ndi) then
ddgdde(j,i)=term2*ddemsd(i,j)
else
ddgdde(j,i)=term2*ddemsd(i-ndi+3,j)
end if
end do
call lubksb (workst, nslptl, 18, indx, ddgdde(1,i))
end do
C----- Derivative of stress increment w.r.t. strain increment
C----- Jacobian matrix: elastic part
do j=1,ntens
do i=1,ntens
ddsdde(i,j)=0.
end do
end do
do j=1,ndi
do i=1,ndi
ddsdde(i,j)=d(i,j)
ddsdde(i,j)=ddsdde(i,j)-stress(i)
end do
end do
if (nshr.gt.0) then
do j=1,nshr
do i=1,nshr
ddsdde(i+ndi,j+ndi)=d(i+3,j+3)
end do
do i=1,ndi
ddsdde(i,j+ndi)=d(i,j+3)
ddsdde(j+ndi,i)=d(j+3,i)
ddsdde(j+ndi,i)=ddsdde(j+ndi,i)-stress(j+ndi)
end do
end do
end if
C----- Jacobian matrix: plastic part (slip)
do j=1,ndi
do i=1,ndi
do k=1,nslptl
ddsdde(i,j)=ddsdde(i,j)-ddemsd(i,k)*ddgdde(k,j)
end do
end do
end do
if (nshr.gt.0) then
do j=1,nshr
do i=1,nshr
do k=1,nslptl
ddsdde(i+ndi,j+ndi)=ddsdde(i+ndi,j+ndi)-
2 ddemsd(i+3,k)*ddgdde(k,j+ndi)
end do
end do
do i=1,ndi
do k=1,nslptl
ddsdde(i,j+ndi)=ddsdde(i,j+ndi)-
2 ddemsd(i,k)*ddgdde(k,j+ndi)
ddsdde(j+ndi,i)=ddsdde(j+ndi,i)-
2 ddemsd(j+3,k)*ddgdde(k,i)
end do
end do
end do
end if
do j=1,ntens
do i=1,ntens
ddsdde(i,j)=ddsdde(i,j)/(1.+dev)
end do
end do
C----- Save solutions (without iteration):
if (nitrtn.eq.0) then
idnor=3*nslptl
iddir=6*nslptl
do j=1,nslptl
tauslp=statev(2*nslptl+j)
if(j>=13.and.tauslp<=0)then
yield=1.e6
else
yield=statev(j)
endif
fslip1(j)=fslip(j)
gslp1(j)=yield
gamma1(j)=statev(nslptl+j)
tausp1(j)=statev(2*nslptl+j)
do i=1,3
idnor=idnor+1
spnor1(i,j)=statev(idnor)
iddir=iddir+1
spdir1(i,j)=statev(iddir)
end do
end do
do j=1,ntens
stres1(j)=stress(j)
do i=1,ntens
ddsde1(i,j)=ddsdde(i,j)
end do
end do
end if
C----- Increments of stress DSOLD, and solution dependent state
C variables DGAMOD, DTAUOD, DGSPOD, DSPNRO, DSPDRO (for the next
C iteration)
do i=1,ntens
dsold(i)=dstres(i)
end do
do j=1,nslptl
dgamod(j)=dgamma(j)
dtauod(j)=dtausp(j)
dgspod(j)=dgslip(j)
do i=1,3
dspnro(i,j)=dspnor(i,j)
dspdro(i,j)=dspdir(i,j)
end do
end do
C----- Check if the iteration solution converges
idback=0
do j=1,nslptl
tauslp=statev(2*nslptl+j)
if(j>=13.and.tauslp<=0)then
yield=1.e6
else
yield=statev(j)
endif
x=statev(2*nslptl+j)/yield
temp=(abs(x))**(50)
f=adot*temp*dsign(1.d0,x)
residu=theta*dtime*f+dtime*(1.0-theta)*
2 fslip1(j)-dgamma(j)
if (abs(residu).gt.gamerr) idback=1
end do
if (idback.ne.0.and.nitrtn.lt.itrmax) then
do i=1,nslptl
gamma(i)=statev(nslptl+1)
enddo
gamtol=statev(9*nslptl+1)
do iself = 1,nslptl
do kderiv = 1,nslptl
dhdgdg(iself,kderiv) = 0.
do latent = 1,nslptl
if(ictype == 1)then
c FCC
term1 = props(12)*gamtol/
* (props(11)-props(10))
term2 = 2.*exp(-term1)/
* (1.+exp(-2.*term1))
term3 = props(12)/(props(11)-props(10))
* *dsign(1.d0,gamma(kderiv))
dhlatn = -2.*props(12)*term2**2
* *tanh(term1)*term3
else
C BASAL
if(iself <= 3)then
dhlatn = 0.
C PRISM
elseif(iself <= 6)then
if(latent <= 12)then
q = 0.2
else
q = 0.5
endif
if(iself == latent)q = 1.
hlatnt = q*props(15)*(1.d0-(props(13)/
* props(14)))*exp(-props(15)*(gamtol/
* props(14)))
dhlatn = q*(-props(15)/props(14))*
* dsign(1.d0,gamma(kderiv))*hlatnt
C PYRM
elseif(iself <= 12)then
if(latent <= 6)then
q = 1.
elseif(latent <= 12)then
q = 0.2
else
q = 0.25
endif
if(iself == latent)q = 1.
hlatnt = q*props(18)*(1.d0-(props(16)/
* props(17)))*exp(-props(18)*(gamtol/
* props(17)))
dhlatn = q*(-props(18)/props(17))*
* dsign(1.d0,gamma(kderiv))*hlatnt
C TWIN
else
if(latent <= 6)then
q = 1.
else
q = 0.2
endif
if(iself == latent)q = 1.
if(gamtol <= props(21))then
dhlatn = 0.
else
dhlatn = q*dsign(1.d0,gamma(kderiv))*
* (props(20)/(props(21)**(props(22)-1.)))
* *(props(22)-1.)*(gamtol**(props(22)-2.))
endif
endif
endif
dhdgdg(iself,kderiv) =
* dhdgdg(iself,kderiv)+dhlatn*
* abs(dgamod(latent))
end do
enddo
end do
go to 1000
elseif (nitrtn.ge.itrmax) then
C----- Solution not converge within maximum number of iteration (the
C solution without iteration will be used)
do j=1,ntens
stress(j)=stres1(j)
do i=1,ntens
ddsdde(i,j)=ddsde1(i,j)
end do
end do
idnor=3*nslptl
iddir=6*nslptl
do j=1,nslptl
statev(j)=gslp1(j)
statev(nslptl+j)=gamma1(j)
statev(2*nslptl+j)=tausp1(j)
do i=1,3
idnor=idnor+1
statev(idnor)=spnor1(i,j)
iddir=iddir+1
statev(iddir)=spdir1(i,j)
end do
end do
end if
C----- Total cumulative shear strains on all slip systems (sum of the
C absolute values of shear strains in all slip systems)
do i=1,nslptl
statev(9*nslptl+1)=statev(9*nslptl+1)+abs(dgamma(i))
end do
return
end
c----------------------------------------------------------------------
subroutine ludcmp (a, n, np, indx, d)
include 'aba_param.inc'
parameter (nmax=200, tiny=1.0e-20)
dimension a(np,np), indx(n), vv(nmax)
d = 1.d0
do i = 1,n
aamax = 0.
do j = 1,n
if (abs(a(i,j)).gt.aamax) aamax = abs(a(i,j))
end do
if (aamax.eq.0.) pause 'singular matrix.'
vv(i) = 1./aamax
end do
do j = 1,n
do i = 1,j-1
sum = a(i,j)
do k = 1,i-1
sum = sum-a(i,k)*a(k,j)
end do
a(i,j) = sum
end do
aamax = 0.
do i = j,n
sum = a(i,j)
do k = 1,j-1
sum = sum-a(i,k)*a(k,j)
end do
a(i,j) = sum
dum = vv(i)*abs(sum)
if (dum.ge.aamax) then
imax = i
aamax = dum
end if
end do
if (j.ne.imax) then
do k = 1,n
dum = a(imax,k)
a(imax,k) = a(j,k)
a(j,k) = dum
end do
d = -d
vv(imax) = vv(j)
end if
indx(j) = imax
if (a(j,j).eq.0.) a(j,j) = tiny
if (j.ne.n) then
dum = 1./a(j,j)
do i = j+1,n
a(i,j) = a(i,j)*dum
end do
end if
end do
return
end
C----------------------------------------------------------------------
subroutine lubksb (a, n, np, indx, b)
include 'aba_param.inc'
dimension a(np,np), indx(n), b(n)
ii = 0
do i = 1,n
ll = indx(i)
sum = b(ll)
b(ll) = b(i)
if (ii.ne.0) then
do j = ii,i-1
sum = sum-a(i,j)*b(j)
end do
else if (sum.ne.0.) then
ii = i
end if
b(i) = sum
end do
do i = n,1,-1
sum = b(i)
if (i.lt.n) then
do j = i+1,n
sum = sum-a(i,j)*b(j)
end do
end if
b(i) = sum/a(i,i)
end do
return
end
