diff --git a/JMBBM13b/README.md b/JMBBM13b/README.md
index 4e53031..d27f2f2 100644
--- a/JMBBM13b/README.md
+++ b/JMBBM13b/README.md
@@ -1,3 +1,4 @@
-Journal Article: https://zenodo.org/records/11184080
+Journal Article: https://doi.org/10.1016/j.jmbbm.2014.01.007
+
+Supporting Data - including original software versions: https://zenodo.org/records/11184080
-Supporting Data - including original software versions: https://doi.org/10.1016/j.jmbbm.2014.01.007
\ No newline at end of file
diff --git a/JMBBM13b/UCrys_HCP_Only.for b/JMBBM13b/UCrys_HCP_Only.for
index fa8b4b6..5217b13 100644
--- a/JMBBM13b/UCrys_HCP_Only.for
+++ b/JMBBM13b/UCrys_HCP_Only.for
@@ -1,1026 +1,1027 @@
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)
+ 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
+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
- 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----- 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
+ 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
+ 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
+ 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
+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 #####
- 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)
+ 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
- end do
-C----- Update the normal to a slip plane and a slip direction
- idnor=3*nslptl
- iddir=6*nslptl
- do j=1,nslptl
+ 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 #####
+ 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
- statev(idnor)=statev(idnor)+dspnor(i,j)-dspnro(i,j)
iddir=iddir+1
- statev(iddir)=statev(iddir)+dspdir(i,j)-dspdro(i,j)
+ statev(idnor)=slpnor(i,j)
+ statev(iddir)=slpdir(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
+ 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
+c 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
+c 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----------------------------------------------------------------------
\ No newline at end of file
+C----------------------------------------------------------------------
diff --git a/JMBBM13b/batch_submission/lxp_system_test.sh b/JMBBM13b/batch_submission/lxp_system_test.sh
new file mode 100644
index 0000000..056a62a
--- /dev/null
+++ b/JMBBM13b/batch_submission/lxp_system_test.sh
@@ -0,0 +1,25 @@
+#!/bin/bash -l
+
+#SBATCH -J AbaqusUMatCheck
+#SBATCH -A myaccount
+#SBATCH --nodes=1
+#SBATCH --time=00:05:00
+#SBATCH -p cpu
+#SBATCH --qos test
+# disable multithreading. Please check both cases to see which gives better performance
+#SBATCH --hint=nomultithread
+#SBATCH --ntasks-per-node=128
+#SBATCH --cpus-per-task=1
+#SBATCH -o %x-%j.log
+
+### Load ABAQUS module
+module load abaqus/2023
+module load intel-compilers
+
+### Configure environment variables, need to unset SLURM's Global Task ID for ABAQUS's PlatformMPI to work
+unset SLURM_GTIDS
+
+### Set input file and job (file prefix) name here
+job_name=${SLURM_JOB_NAME}
+
+abaqus job=${job_name} input=CDIE_1E.inp user=UCrys_HCP_Only double=both interactive
\ No newline at end of file