适合HCP mg合金的黄永刚umat的晶体塑性

晶体塑性有限元

2025年5月18日 16:30

参考文献：《Computational micromechanics of bioabsorbable magnesium stents》

适合HCP mg合金的黄永刚umat的晶体塑性的图1

文章的子程序中考虑了孪晶的影响，感兴趣的可以下载调试。子程序如下（部分子函数可以参考）：

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)

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

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

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

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

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

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

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

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

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

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

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

do j = 1,6

nslip = nslip+1

do i = 1,3

slpdir(i,nslip) = rwkdir(i,j)

slpnor(i,nslip) = rwknor(i,j)

enddo

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

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

do i = 1,3

slpdir(i,j) = termd(i)

slpnor(i,j) = termn(i)

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

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

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

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

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

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

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

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

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

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

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

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

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 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

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

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

do j=1,ndi

do i=1,ndi

ddsdde(i,j)=d(i,j)

ddsdde(i,j)=ddsdde(i,j)-stress(i)

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 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

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

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 if

do j=1,ntens

do i=1,ntens

ddsdde(i,j)=ddsdde(i,j)/(1.+dev)

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

do j=1,ntens

stres1(j)=stress(j)

do i=1,ntens

ddsde1(i,j)=ddsdde(i,j)

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

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

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

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 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

C----------------------------------------------------------------------

相关子函数已经传入知识星球，可以直接加入下载

适合HCP mg合金的黄永刚umat的晶体塑性的图2

登录后免费查看全文

立即登录

技术邻APP
工程师必备

项目客服
培训客服
平台客服

TOP

适合HCP mg合金的黄永刚umat的晶体塑性

推荐阅读