From 684f81980d8d79f77469acefc31421dae9c082c5 Mon Sep 17 00:00:00 2001 From: jgrogan Date: Tue, 4 Jun 2024 12:32:16 +0100 Subject: [PATCH] Initial hpc batch submission script and clean whitespace in umat. --- JMBBM13b/README.md | 5 +- JMBBM13b/UCrys_HCP_Only.for | 2007 +++++++++--------- JMBBM13b/batch_submission/lxp_system_test.sh | 25 + 3 files changed, 1032 insertions(+), 1005 deletions(-) create mode 100644 JMBBM13b/batch_submission/lxp_system_test.sh 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