757 lines
20 KiB
Text
757 lines
20 KiB
Text
|
subroutine vumat(nblock, ndir, nshr, nstatev, nfieldv, nprops,
|
||
|
|
* lanneal, steptime, totaltime, dt, cmname, coordmp, charlength,
|
||
|
|
* props, density, straininc, relspininc, tempold, stretchold,
|
||
|
|
* defgradold, fieldold, stressold, stateold, enerinternold,
|
||
|
|
* enerinelasold, tempnew, stretchnew, defgradnew, fieldnew,
|
||
|
|
* stressnew, statenew, enerinternnew, enerinelasnew)
|
||
|
|
c
|
||
|
|
include 'vaba_param.inc'
|
||
|
|
c modified from harewood vumat- jg:20/07/12
|
||
|
|
dimension stressold(nblock,ndir+nshr),
|
||
|
|
1 stressnew(nblock,ndir+nshr),
|
||
|
|
1 stateold(nblock,nstatev),statenew(nblock,nstatev),
|
||
|
|
1 straininc(nblock, ndir+nshr)
|
||
|
|
dimension slpdir(3,18),slpnor(3,18),slpdef(6,18),slpspn(3,18),
|
||
|
|
1 dspdir(3,18),dspnor(3,18),d(6,6),fslip(18),dfdxsp(18),
|
||
|
|
1 ddemsd(6,18),h(18,18),dgamma(18),dtausp(18),dgslip(18),
|
||
|
|
1 dstres(6),delats(6),dvgrad(3,3),dspin(3),workst(18,18),
|
||
|
|
1 indx(18),term(3,3),trm0(3,3),props(nprops),itrm(3),
|
||
|
|
1 rotate(3,3),rwkdir(3,24),rwknor(3,24),indxL(3),termd(3),
|
||
|
|
1 termn(3),tauslpA(18)
|
||
|
|
c
|
||
|
|
do km=1,nblock
|
||
|
|
C----- As the VUMAT passes in tensor shear strain and this subroutine
|
||
|
|
C----- uses engineering strain --> STRAININC(shr) x 2
|
||
|
|
do i=1,nshr
|
||
|
|
straininc(km,i+3)=straininc(km,i+3)*2.d0
|
||
|
|
end do
|
||
|
|
if (nshr.gt.1) then
|
||
|
|
save=straininc(km,5)
|
||
|
|
straininc(km,5)=straininc(km,6)
|
||
|
|
straininc(km,6)=save
|
||
|
|
save=stressold(km,5)
|
||
|
|
stressold(km,5)=stressold(km,6)
|
||
|
|
stressold(km,6)=save
|
||
|
|
end if
|
||
|
|
C Elastic Matrix {D}
|
||
|
|
gshear = props(1)/(2.d0*(1.d0+props(2)))
|
||
|
|
e11 = 2.d0*gshear*(1.d0-props(2))/(1.d0-2.d0*props(2))
|
||
|
|
e12 = 2.d0*gshear*props(2)/(1.d0-2.d0*props(2))
|
||
|
|
d = 0.d0
|
||
|
|
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 Lin Elastic Response for Exp/Packager
|
||
|
|
if(totaltime==0.)then
|
||
|
|
C Calculation of Stress Inc
|
||
|
|
dstres=0.d0
|
||
|
|
do i=1,ndir+nshr
|
||
|
|
do j=1,ndir+nshr
|
||
|
|
dstres(i)=dstres(i)+d(i,j)*straininc(km,j)
|
||
|
|
end do
|
||
|
|
end do
|
||
|
|
C Calculation of Stress New
|
||
|
|
do i=1,ndir+nshr
|
||
|
|
stressnew(km,i)=stressold(km,i)+dstres(i)
|
||
|
|
end do
|
||
|
|
cycle
|
||
|
|
endif
|
||
|
|
c------ Crystal Type:
|
||
|
|
ictype=nint(props(9))
|
||
|
|
if(ictype == 2)then
|
||
|
|
nslptl = 6
|
||
|
|
elseif(ictype == 3)then
|
||
|
|
nslptl = 12
|
||
|
|
else
|
||
|
|
nslptl = 18
|
||
|
|
endif
|
||
|
|
C----- Integration parameter: THETA
|
||
|
|
theta = 0.5d0
|
||
|
|
term = 0.d0
|
||
|
|
trm0 = 0.d0
|
||
|
|
do i=1,3
|
||
|
|
term(i,i)=2.d0
|
||
|
|
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=0.d0
|
||
|
|
do i=1,ndir
|
||
|
|
dev=dev+straininc(km,i)
|
||
|
|
end do
|
||
|
|
C----- Check whether the current stress state is the initial state
|
||
|
|
if (stateold(km,1).eq.0.) then
|
||
|
|
c ##### basal and prismatic#####
|
||
|
|
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.
|
||
|
|
do i = 1,3
|
||
|
|
do k = 1,3
|
||
|
|
slpdir(k,i) = rwkdir(k,i)
|
||
|
|
slpnor(k,i) = rwknor(k,1)
|
||
|
|
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
|
||
|
|
if(ictype >= 3)then
|
||
|
|
c ##### 2nd order pyramidal <a+c> #####
|
||
|
|
aspect = 1.623
|
||
|
|
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
|
||
|
|
rwknor(1,1) = aspect*sangle
|
||
|
|
rwknor(2,1) = 3.*aspect*cangle
|
||
|
|
rwknor(3,1) = 4.*sangle*cangle
|
||
|
|
rwknor(1,2) = aspect*sangle
|
||
|
|
rwknor(2,2) = -3.*aspect*cangle
|
||
|
|
rwknor(3,2) = 4.*sangle*cangle
|
||
|
|
rwknor(1,3) = 2.*aspect*sangle
|
||
|
|
rwknor(2,3) = 0.
|
||
|
|
rwknor(3,3) = 4.*sangle*cangle
|
||
|
|
do j = 4,6
|
||
|
|
rwknor(1,j) = rwknor(1,j-3)
|
||
|
|
rwknor(2,j) = rwknor(2,j-3)
|
||
|
|
rwknor(3,j) = -rwknor(3,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
|
||
|
|
rwknor(1,1) = 0.
|
||
|
|
rwknor(2,1) = -2.*aspect*cangle
|
||
|
|
rwknor(3,1) = -4.*sangle*cangle
|
||
|
|
rwknor(1,2) = -aspect*sangle
|
||
|
|
rwknor(2,2) = -aspect*cangle
|
||
|
|
rwknor(3,2) = -4.*sangle*cangle
|
||
|
|
rwknor(1,3) = -aspect*sangle
|
||
|
|
rwknor(2,3) = aspect*cangle
|
||
|
|
rwknor(3,3) = -4.*sangle*cangle
|
||
|
|
do j = 4,6
|
||
|
|
rwknor(1,j) = -rwknor(1,j-3)
|
||
|
|
rwknor(2,j) = -rwknor(2,j-3)
|
||
|
|
rwknor(3,j) = rwknor(3,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
|
||
|
|
rotate(j,j) = 1.
|
||
|
|
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----- 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
|
||
|
|
stateold(km,nstatev) = nslptl
|
||
|
|
idnor=3*nslptl
|
||
|
|
iddir=6*nslptl
|
||
|
|
do j=1,nslptl
|
||
|
|
do i=1,3
|
||
|
|
idnor=idnor+1
|
||
|
|
stateold(km,idnor)=slpnor(i,j)
|
||
|
|
iddir=iddir+1
|
||
|
|
stateold(km,iddir)=slpdir(i,j)
|
||
|
|
end do
|
||
|
|
end do
|
||
|
|
C----- Initial value of the current strength for all slip systems
|
||
|
|
do j=1,nslptl
|
||
|
|
if(j<=3)then
|
||
|
|
stateold(km,j)=props(10)
|
||
|
|
elseif(j<=6)then
|
||
|
|
stateold(km,j)=props(13)
|
||
|
|
elseif(j<=12)then
|
||
|
|
stateold(km,j)=props(16)
|
||
|
|
else
|
||
|
|
stateold(km,j)=props(19)
|
||
|
|
endif
|
||
|
|
enddo
|
||
|
|
C----- Initial value of shear strain in slip systems
|
||
|
|
do i=1,nslptl
|
||
|
|
stateold(km,nslptl+i)=0.d0
|
||
|
|
end do
|
||
|
|
stateold(km,9*nslptl+1)=0.d0
|
||
|
|
C----- Initial value of the resolved shear stress in slip systems
|
||
|
|
do i=1,nslptl
|
||
|
|
term1=0.
|
||
|
|
do j=1,ndir+nshr
|
||
|
|
if (j.le.ndir) then
|
||
|
|
term1=term1+slpdef(j,i)*stressold(km,j)
|
||
|
|
else
|
||
|
|
term1=term1+slpdef(j-ndir+3,i)*stressold(km,j)
|
||
|
|
end if
|
||
|
|
end do
|
||
|
|
stateold(km,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
|
||
|
|
slpnor(i,j)=stateold(km,idnor)
|
||
|
|
iddir=iddir+1
|
||
|
|
slpdir(i,j)=stateold(km,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.5d0*(slpdir(1,j)*slpnor(2,j)-
|
||
|
|
2 slpdir(2,j)*slpnor(1,j))
|
||
|
|
slpspn(2,j)=0.5d0*(slpdir(3,j)*slpnor(1,j)-
|
||
|
|
2 slpdir(1,j)*slpnor(3,j))
|
||
|
|
slpspn(3,j)=0.5d0*(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.d0
|
||
|
|
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)*stressold(km,1)
|
||
|
|
ddemsd(5,j)=ddemsd(5,j)+slpspn(2,j)*stressold(km,1)
|
||
|
|
if (ndir.gt.1) then
|
||
|
|
ddemsd(4,j)=ddemsd(4,j)+slpspn(1,j)*stressold(km,2)
|
||
|
|
ddemsd(6,j)=ddemsd(6,j)-slpspn(3,j)*stressold(km,2)
|
||
|
|
end if
|
||
|
|
if (ndir.gt.2) then
|
||
|
|
ddemsd(5,j)=ddemsd(5,j)-slpspn(2,j)*stressold(km,3)
|
||
|
|
ddemsd(6,j)=ddemsd(6,j)+slpspn(3,j)*stressold(km,3)
|
||
|
|
end if
|
||
|
|
if (nshr.ge.1) then
|
||
|
|
ddemsd(1,j)=ddemsd(1,j)+slpspn(1,j)
|
||
|
|
2 *stressold(km,ndir+1)
|
||
|
|
ddemsd(2,j)=ddemsd(2,j)-slpspn(1,j)
|
||
|
|
2 *stressold(km,ndir+1)
|
||
|
|
ddemsd(5,j)=ddemsd(5,j)-slpspn(3,j)
|
||
|
|
2 *stressold(km,ndir+1)
|
||
|
|
ddemsd(6,j)=ddemsd(6,j)+slpspn(2,j)
|
||
|
|
2 *stressold(km,ndir+1)
|
||
|
|
end if
|
||
|
|
if (nshr.ge.2) then
|
||
|
|
ddemsd(1,j)=ddemsd(1,j)-slpspn(2,j)
|
||
|
|
2 *stressold(km,ndir+2)
|
||
|
|
ddemsd(3,j)=ddemsd(3,j)+slpspn(2,j)
|
||
|
|
2 *stressold(km,ndir+2)
|
||
|
|
ddemsd(4,j)=ddemsd(4,j)+slpspn(3,j)
|
||
|
|
2 *stressold(km,ndir+2)
|
||
|
|
ddemsd(6,j)=ddemsd(6,j)-slpspn(1,j)
|
||
|
|
2 *stressold(km,ndir+2)
|
||
|
|
end if
|
||
|
|
if (nshr.eq.3) then
|
||
|
|
ddemsd(2,j)=ddemsd(2,j)+slpspn(3,j)
|
||
|
|
2 *stressold(km,ndir+3)
|
||
|
|
ddemsd(3,j)=ddemsd(3,j)-slpspn(3,j)
|
||
|
|
2 *stressold(km,ndir+3)
|
||
|
|
ddemsd(4,j)=ddemsd(4,j)-slpspn(2,j)
|
||
|
|
2 *stressold(km,ndir+3)
|
||
|
|
ddemsd(5,j)=ddemsd(5,j)+slpspn(1,j)
|
||
|
|
2 *stressold(km,ndir+3)
|
||
|
|
end if
|
||
|
|
end do
|
||
|
|
C----- Shear strain-rate in a slip system at the start of increment:
|
||
|
|
do i=1,nslptl
|
||
|
|
tauslp=stateold(km,2*nslptl+i)
|
||
|
|
if(i>=13.and.tauslp<=0)then
|
||
|
|
yield=1.e6
|
||
|
|
else
|
||
|
|
yield=stateold(km,i)
|
||
|
|
endif
|
||
|
|
x=tauslp/yield
|
||
|
|
fslip(i)=0.001d0*((abs(x))**50.)*dsign(1.d0,x)
|
||
|
|
dfdxsp(i)=50.d0*0.001d0*(abs(x))**(50.-1.0)
|
||
|
|
end do
|
||
|
|
C----- Self- and latent-hardening laws
|
||
|
|
gamtol=stateold(km,9*nslptl+1)
|
||
|
|
do iself = 1,nslptl
|
||
|
|
do latent = 1,nslptl
|
||
|
|
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
|
||
|
|
h(iself,latent) = hlatnt
|
||
|
|
enddo
|
||
|
|
end do
|
||
|
|
C----- Solve the increment of shear strain in a slip system
|
||
|
|
term1=theta*dt
|
||
|
|
do i=1,nslptl
|
||
|
|
tauslp=stateold(km,2*nslptl+i)
|
||
|
|
if(i>=13.and.tauslp<=0)then
|
||
|
|
yield=1.e6
|
||
|
|
else
|
||
|
|
yield=stateold(km,i)
|
||
|
|
endif
|
||
|
|
x=tauslp/yield
|
||
|
|
term2=term1*dfdxsp(i)/yield
|
||
|
|
term3=term1*x*dfdxsp(i)/yield
|
||
|
|
do j=1,nslptl
|
||
|
|
term4=0.d0
|
||
|
|
do k=1,6
|
||
|
|
term4=term4+ddemsd(k,i)*slpdef(k,j)
|
||
|
|
end do
|
||
|
|
workst(i,j)=term2*term4+h(i,j)*term3
|
||
|
|
2 *dsign(1.d0,fslip(j))
|
||
|
|
end do
|
||
|
|
workst(i,i)=workst(i,i)+1.d0
|
||
|
|
end do
|
||
|
|
call ludcmp (workst, nslptl, 18, indx, ddcmp)
|
||
|
|
C----- Increment of shear strain in a slip system
|
||
|
|
term1=theta*dt
|
||
|
|
do i=1,nslptl
|
||
|
|
tauslp=stateold(km,2*nslptl+i)
|
||
|
|
if(i>=13.and.tauslp<=0)then
|
||
|
|
yield=1.e6
|
||
|
|
else
|
||
|
|
yield=stateold(km,i)
|
||
|
|
endif
|
||
|
|
term2=term1*dfdxsp(i)/yield
|
||
|
|
dgamma(i)=0.
|
||
|
|
do j=1,ndir
|
||
|
|
dgamma(i)=dgamma(i)+ddemsd(j,i)*straininc(km,j)
|
||
|
|
end do
|
||
|
|
if (nshr.gt.0) then
|
||
|
|
do j=1,nshr
|
||
|
|
if (j.eq.1) then
|
||
|
|
dgamma(i)=dgamma(i)+ddemsd(4,i)
|
||
|
|
2 *straininc(km,4)
|
||
|
|
elseif (j.eq.2) then
|
||
|
|
dgamma(i)=dgamma(i)+ddemsd(6,i)
|
||
|
|
2 *straininc(km,5)
|
||
|
|
elseif (j.eq.3) then
|
||
|
|
dgamma(i)=dgamma(i)+ddemsd(5,i)
|
||
|
|
2 *straininc(km,6)
|
||
|
|
end if
|
||
|
|
end do
|
||
|
|
end if
|
||
|
|
dgamma(i)=dgamma(i)*term2+fslip(i)*dt
|
||
|
|
end do
|
||
|
|
call lubksb (workst, nslptl, 18, indx, dgamma)
|
||
|
|
C----- Update the shear strain in a slip system:
|
||
|
|
do i=1,nslptl
|
||
|
|
stateold(km,nslptl+i)=stateold(km,nslptl+i)+dgamma(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
|
||
|
|
stateold(km,i)=stateold(km,i)+dgslip(i)
|
||
|
|
end do
|
||
|
|
C----- Increment of strain associated with lattice stretching:
|
||
|
|
delats=0.
|
||
|
|
do j=1,3
|
||
|
|
if (j.le.ndir) delats(j)=straininc(km,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)=straininc(km,j+ndir)
|
||
|
|
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
|
||
|
|
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
|
||
|
|
stateold(km,2*nslptl+i)=stateold(km,2*nslptl+i)
|
||
|
|
2 +dtausp(i)
|
||
|
|
end do
|
||
|
|
C----- Increment of stress: DSTRES
|
||
|
|
do i=1,ndir+nshr
|
||
|
|
dstres(i)=-stressold(km,i)*dev
|
||
|
|
end do
|
||
|
|
do i=1,ndir
|
||
|
|
do j=1,ndir
|
||
|
|
dstres(i)=dstres(i)+d(i,j)*straininc(km,j)
|
||
|
|
end do
|
||
|
|
if (nshr.gt.0)then
|
||
|
|
do j=1,nshr
|
||
|
|
dstres(i)=dstres(i)+d(i,j+3)*
|
||
|
|
2 straininc(km,j+ndir)
|
||
|
|
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,ndir
|
||
|
|
dstres(i+ndir)=dstres(i+ndir)+d(i+3,j)
|
||
|
|
2 *straininc(km,j)
|
||
|
|
end do
|
||
|
|
do j=1,nshr
|
||
|
|
dstres(i+ndir)=dstres(i+ndir)+d(i+3,j+3)*
|
||
|
|
2 straininc(km,j+ndir)
|
||
|
|
end do
|
||
|
|
do j=1,nslptl
|
||
|
|
dstres(i+ndir)=dstres(i+ndir)-ddemsd(i+3,j)
|
||
|
|
2 *dgamma(j)
|
||
|
|
end do
|
||
|
|
end do
|
||
|
|
end if
|
||
|
|
C----- Update the stress: STRESSOLD
|
||
|
|
do i=1,ndir+nshr
|
||
|
|
stressold(km,i)=stressold(km,i)+dstres(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
|
||
|
|
stateold(km,idnor)=stateold(km,idnor)+dspnor(i,j)
|
||
|
|
iddir=iddir+1
|
||
|
|
stateold(km,iddir)=stateold(km,iddir)+dspdir(i,j)
|
||
|
|
end do
|
||
|
|
end do
|
||
|
|
C----- Total cumulative shear strains on all slip systems
|
||
|
|
do i=1,nslptl
|
||
|
|
stateold(km,9*nslptl+1)=stateold(km,9*nslptl+1)
|
||
|
|
2 +abs(dgamma(i))
|
||
|
|
end do
|
||
|
|
c----- update stressold to stressnew
|
||
|
|
do i=1,ndir+nshr
|
||
|
|
stressnew(km,i)=stressold(km,i)
|
||
|
|
end do
|
||
|
|
c----- update stateold to statenew for 1 - nstatev
|
||
|
|
do i=1,nstatev
|
||
|
|
statenew(km,i)=stateold(km,i)
|
||
|
|
end do
|
||
|
|
if (nshr.gt.1) then
|
||
|
|
save=straininc(km,5)
|
||
|
|
straininc(km,5)=straininc(km,6)
|
||
|
|
straininc(km,6)=save
|
||
|
|
save=stressnew(km,5)
|
||
|
|
stressnew(km,5)=stressnew(km,6)
|
||
|
|
stressnew(km,6)=save
|
||
|
|
end if
|
||
|
|
enddo
|
||
|
|
return
|
||
|
|
end
|
||
|
|
c----------------------------------------------------------------------
|
||
|
|
subroutine ludcmp (a, n, np, indx, d)
|
||
|
|
include 'vaba_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 'vaba_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
|
||
|
|
C----------------------------------------------------------------------
|