HCP improvements

kMaterialParam.f

@@ -63,7 +63,7 @@ SUBROUTINE kMaterialParam(iphase,caratio,compliance,G12,thermat,
       ! materials available are the following:
       ! 'zirconium', 'alphauranium', 'tungsten', 'copper', 'carbide',
       ! 'olivine'
-      character(len=*), parameter :: matname = 'alphauranium' 
+      character(len=*), parameter :: matname = 'zircaloy4' 
 
 **       End of parameters to set       **
 ******************************************
@@ -112,17 +112,55 @@ SUBROUTINE kMaterialParam(iphase,caratio,compliance,G12,thermat,
         v13 = 0.04
 
         ! CRSS (MPa units)
-        XTAUC1 = 15.2 ! basal
-        XTAUC2 = 67.7 ! prismatic
+        XTAUC1 = 15.2   ! basal
+        XTAUC2 = 67.7   ! prismatic
         XTAUC4 = 2000.0 ! pyramidal
-
+        XTAUC3 = 2000.0 ! Pyramidal 1st <c+a>
+        XTAUC5 = 2000.0 ! Pyramidal 2nd <c+a>
+
         ! thermal expansion coefficients
         alpha1 = 9.5D-6
         alpha2 = alpha1
         alpha3 = 0.5895*alpha1
 
         ! prefactor for SSD evolution equation
        	gammast = 0.0
+
+        case('zircaloy4')
+
+        crssratio = 3.48      
+        caratio = 1.593 !Updated for Zr			
+        ! Burgers vectors (um)
+        burger1 = 3.20E-4 ! Metallic atomic radius of 160pm - 
+                        ! converted to microns and DOubled.
+        burger2 = 3.48*burger1 !1.88*burger1
+
+        ! elastic constants (MPa units)
+        TEMP=0
+        ! Calibrated temperature constants for Zircaloy-4 (Yang Liu)
+	  E1 = (-0.0755*TEMP + 120.4411)*1.0E3
+	  E3 = (-0.0327*TEMP + 132.8621)*1.0E3
+	  v12 = 3.4273E-4*TEMP + 0.3002
+	  v13 = -9.1182E-5*TEMP + 0.2642 
+	  G13 = (-0.0233*TEMP + 38.8367)*1.0E3
+	  G12 = E1/(2.0*(1.0+v12))
+
+        ! CRSS (MPa units)
+        TAUCSSD=-0.2576*TEMP+228.4697
+        XTAUC2 = TAUCSSD!*0.92     ! Prismatic <a>  *reduce CRSS by 8 % for fitting Farrell et al. !
+        XTAUC4 = XTAUC2           ! Pyramidal <a>
+        XTAUC1 = 1.3333D0*XTAUC2  ! Basal
+        XTAUC3 = crssratio*XTAUC2 ! Pyramidal 1st <c+a>
+        XTAUC5 = XTAUC3*5.00D0    ! Pyramidal 2nd <c+a>
+
+        ! thermal expansion coefficients
+        ! For Zr4 from Abdolvand et al., 2015 Acta Mater 93, 235-245
+        alpha1 = 1.01D-5 
+        alpha2 = alpha1 
+        alpha3 = 0.52475*alpha1
+
+        ! prefactor for SSD evolution equation
+        gammast = 130.0
 
         case default
         WRITE(*,*)"Not sure what material"
@@ -134,13 +172,15 @@ SUBROUTINE kMaterialParam(iphase,caratio,compliance,G12,thermat,
       v23 = v13
 
       ! assign Burgers vector scalars
-      burgerv(1:6) = burger1
-      burgerv(7:12) = burger2
+      burgerv(1:12) = burger1 ! all <a> type slip systems
+      burgerv(13:nSys) = burger2 ! all <a+c> types
 
       ! assign CRSS
-      tauc(1:3) = XTAUC1
-      tauc(4:6) = XTAUC2
-      tauc(7:12) = XTAUC4
+      tauc(1:3) = XTAUC1     ! Basal <a>
+      tauc(4:6) = XTAUC2     ! Prismatic <a>
+      tauc(7:12) = XTAUC4    ! Pyramidal <a>
+      tauc(13:24) = XTAUC3   ! Pyramidal 1st <c+a>
+      tauc(25:30) = XTAUC5   ! Pyramidal 2nd <c+a>       
 
       case(1) !bcc
 

kdirns.f

@@ -17,8 +17,8 @@ subroutine kdirns(xrot,TwinRot,iphase,L,nTwin,ddir1,dnor1,dtwindir1,dtwinnor1,
          ddir1 = 0.0; dnor1 = 0.0 
 
       if(iphase .eq. 0) then !hcp
-         include 'xDir0ET.f'
-         include 'xNorm0ET.f'         
+         include 'xDir0CH.f'
+         include 'xNorm0CH.f'         
       else if(iphase .eq. 1) then !bcc (24/48)
          include 'xDir1.f'
          include 'xNorm1.f'   

kmat.f

@@ -96,7 +96,7 @@ SUBROUTINE kmat(dtime,nsvars,usvars,xI,jelem,kint,time,F,
       ! 5 = Original slip rule with GND coupling 
       ! 6 = Slip rule with constants alpha and beta: alpha.sinh[ beta(tau-tauc)sgn(tau) ]
       ! 7 = Powerlaw plasticity
-      integer, parameter :: kslip = 7
+      integer, parameter :: kslip = 6
 
       ! initial temperature and temperature rate
       real*8, parameter :: Temperature = 293.0
@@ -125,7 +125,7 @@ SUBROUTINE kmat(dtime,nsvars,usvars,xI,jelem,kint,time,F,
       INTEGER, parameter :: KM=6,KN=6
 
       ! number of active slip systems considered
-      integer, parameter :: L0=12 ! HCP
+      integer, parameter :: L0=30 ! HCP
       integer, parameter :: L1=12 ! BCC
       integer, parameter :: L2=12 ! FCC
       integer, parameter :: L4=7  ! Olivine
@@ -676,7 +676,7 @@ SUBROUTINE kmat(dtime,nsvars,usvars,xI,jelem,kint,time,F,
 
       ELSE IF (kslip == 6) THEN ! updated slip rule with GND coupling
 
-      CALL  kslip6ET(xNorm,xDir,tau,signtau,tauc,burgerv,dtime,nSys,
+      CALL  kslip6ET(xNorm,xDir,tau,signtau,tauc,burgerv,dtime, nSys,
      +        iphase,irradiate,gndcut,gndtot,rhossd,Lp,tmat,gammaDot)
 
       ELSE IF (kslip == 7) THEN ! Powerlaw plasticity, orthorombic alpha-Uranium

kslip6ET.f

@@ -12,7 +12,7 @@ subroutine kslip6ET(xNorm,xDir,tau,signtau,tauc,burgerv,dtime,
       real*8, intent(out) :: Lp(3,3),tmat(6,6), gammaDot(nSys)
 
       integer :: i
-      real*8  :: alpha,beta,result1, rhom,dF,f,T,k,gamma0,b,psi,V,
+      real*8  :: alpha,beta,xtemp,result1, rhom,dF,f,T,k,gamma0,b,psi,V,
      + SNij(3,3),sni(6),nsi(6),SNNS(6,6),NSij(3,3),result4(6,6),
      + tempNorm(3), tempDir(3)
 
@@ -41,16 +41,40 @@ subroutine kslip6ET(xNorm,xDir,tau,signtau,tauc,burgerv,dtime,
           V = b*b/sqrt(psi*rhossd) ! activation volume /micron^3
 
           beta = gamma0*V/(k*T) ! rate sensitivity 1/MPa
+C
+      tmat=0.; Lp = 0.;result4=0.
+
+      !rhogndold=sum(gndold)
+      Do I=1,nSys
+
+         if (tau(I) >= tauc(I)) THEN				
+
+             gammaDot(I)=alpha*sinh(beta*signtau(I)*(tau(I) - tauc(I)) )
+
+              tempNorm = xNorm(I,:); 
+              tempDir = xDir(I,:)
+              SNij = spread(tempDir,2,3)*spread(tempNorm,1,3)
+              NSij = spread(tempNorm,2,3)*spread(tempDir,1,3)
+              CALL KGMATVEC6(SNij,sni)         
+              CALL KGMATVEC6(NSij,nsi) 
+              SNNS = spread(sni,2,6)*spread(nsi,1,6)
+              result1 = cosh(beta*signtau(I)*(tau(I) - tauc(I)))
+
+              result4 = result4 + alpha*beta*dtime*result1*SNNS          
+              Lp = Lp + gammaDot(I)*SNij
+          else
+              gammaDot(I)=0.0
+          end if
+C
+      END DO
+
+      tmat = 0.5*(result4+transpose(result4))       
 
       elseif (iphase == 2) then
 
           alpha = 0.02
-          beta = 0.1
-
-      endif
+          beta = 0.1         
 
-      !ne = microns, stress = MPa, F = microN, therefore E = pJ
-
 C
       tmat=0.; Lp = 0.;result4=0.
 
@@ -79,6 +103,53 @@ subroutine kslip6ET(xNorm,xDir,tau,signtau,tauc,burgerv,dtime,
       END DO
 
       tmat = 0.5*(result4+transpose(result4))
+
+      elseif (iphase == 0) then
+
+      rhom = 0.01                                 
+      dF = 5E-8    ! for current stage                        
+      f = 1.0E+11 ! This is in s^-1 TIME IS IN SECONDS            
+      k = 1.381E-11
+      xtemp = 293.0      
+
+C
+      tmat=0.; Lp = 0.;result4=0.
+
+      !rhogndold=sum(gndold)
+      Do I=1,nSys     
+         V = 1.0/sqrt(rhom*2.7778E+6)*burgerv(I)*burgerv(I) ! the factor 2.78e+6 originally came from a parameter called gammazero but the meaning of this is unknown atm.
+         beta = V/(k*xtemp) 
+         alpha = rhom*burgerv(I)*burgerv(I)*f*
+     +  exp(-dF/(k*xtemp))
+
+         if (tau(I) >= tauc(I)) THEN				
+
+             gammaDot(I)=alpha*sinh(beta*signtau(I)*(tau(I) - tauc(I)) )
+
+              tempNorm = xNorm(I,:); 
+              tempDir = xDir(I,:)
+              SNij = spread(tempDir,2,3)*spread(tempNorm,1,3)
+              NSij = spread(tempNorm,2,3)*spread(tempDir,1,3)
+              CALL KGMATVEC6(SNij,sni)         
+              CALL KGMATVEC6(NSij,nsi) 
+              SNNS = spread(sni,2,6)*spread(nsi,1,6)
+              result1 = cosh(beta*signtau(I)*(tau(I) - tauc(I)))
+
+              result4 = result4 + alpha*beta*dtime*result1*SNNS          
+              Lp = Lp + gammaDot(I)*SNij
+          else
+              gammaDot(I)=0.0
+          end if
+C
+      END DO
+
+      tmat = 0.5*(result4+transpose(result4))
+
+      endif
+
+      !ne = microns, stress = MPa, F = microN, therefore E = pJ
+
+
 
       return
       end 

psiAnn.f

@@ -0,0 +1,37 @@
+! Annihilation Constants for defect-dislocation interactions in HCP
+
+      psiAnn(:,1) = (/0.5,0.1,0.1/)
+	psiAnn(:,2) = (/0.1,0.5,0.1/)
+	psiAnn(:,3) = (/0.1,0.1,0.5/)
+
+	psiAnn(:,4) = (/0.5,0.1,0.1/)
+	psiAnn(:,5) = (/0.1,0.5,0.1/)
+	psiAnn(:,6) = (/0.1,0.1,0.5/)
+
+	psiAnn(:,7) = (/0.5,0.1,0.1/)
+	psiAnn(:,8) = (/0.1,0.5,0.1/)
+	psiAnn(:,9) = (/0.1,0.1,0.5/)
+	psiAnn(:,10) = (/0.5,0.1,0.1/)
+	psiAnn(:,11) = (/0.1,0.5,0.1/)
+	psiAnn(:,12) = (/0.1,0.1,0.5/)
+
+	psiAnn(:,13) = (/0.5,0.1,0.1/)
+	psiAnn(:,14) = (/0.1,0.5,0.1/)
+	psiAnn(:,15) = (/0.1,0.1,0.5/)
+	psiAnn(:,16) = (/0.5,0.1,0.1/)
+	psiAnn(:,17) = (/0.1,0.5,0.1/)
+	psiAnn(:,18) = (/0.1,0.1,0.5/)
+
+	psiAnn(:,19) = (/0.5,0.1,0.1/)
+	psiAnn(:,20) = (/0.1,0.5,0.1/)
+	psiAnn(:,21) = (/0.1,0.1,0.5/)
+	psiAnn(:,22) = (/0.5,0.1,0.1/)
+	psiAnn(:,23) = (/0.1,0.5,0.1/)
+	psiAnn(:,24) = (/0.1,0.1,0.5/)
+
+      psiAnn(:,25) = (/0.5,0.1,0.1/)
+      psiAnn(:,26) = (/0.1,0.5,0.1/)
+      psiAnn(:,27) = (/0.1,0.1,0.5/) 
+	psiAnn(:,28) = (/0.5,0.1,0.1/)
+	psiAnn(:,29) = (/0.1,0.5,0.1/)
+	psiAnn(:,30) = (/0.1,0.1,0.5/)

umat.for

@@ -1,7 +1,7 @@
 ***********************************************************
 ** Crystal Plasticity UMAT 
-**  University of Oxford   
-** Developed by Nicolo Grilli 2020,
+**  University of Oxford   
+** Developed by Nicolo Grilli 2020,
 ** rewrite of UMAT by Ed Tarleton which was based on UEL by Fionn Dunne 2007
 **********************************************************
 
@@ -113,7 +113,7 @@
       integer :: ns
 
       ! number of active slip systems considered
-      integer :: L0=12 ! HCP
+      integer :: L0=30 ! HCP
       integer :: L1=12 ! BCC
       integer :: L2=12 ! FCC
       integer :: L4=7  ! Olivine
@@ -437,12 +437,13 @@ C   A FULL INTEGRATION GRADIENT SCHEME
       include 'kcurlET.f'
       include 'kshapes.f'
       include 'utils.f'
-      include 'UEL.for'
+      !include 'UEL.for'
       include 'ktwinneighbourhood.f'
       include 'kRhoTwinInit.f'
       include 'kMaterialParam.f'
       include 'kCRSS.f'
       include 'kHardening.f'
-
+      !include 'psiAnn.f'
+      !include 'xAlpha.f'
 
 

xAlpha.f

@@ -0,0 +1,37 @@
+! Hardening Coefficients for dislocation-defect interactions in HCP
+
+           xAlpha(1,:) = (/0.1,0.2,0.2/)
+		 xAlpha(2,:) = (/0.2,0.1,0.2/)
+		 xAlpha(3,:) = (/0.2,0.2,0.1/)
+
+		 xAlpha(4,:) = (/0.1,0.2,0.2/)
+		 xAlpha(5,:) = (/0.2,0.1,0.2/)
+		 xAlpha(6,:) = (/0.2,0.2,0.1/)
+
+		 xAlpha(7,:) = (/0.1,0.2,0.2/)
+		 xAlpha(8,:) = (/0.2,0.1,0.2/)
+		 xAlpha(9,:) = (/0.2,0.2,0.1/)
+		 xAlpha(10,:) = (/0.1,0.2,0.2/)
+		 xAlpha(11,:) = (/0.2,0.1,0.2/)
+		 xAlpha(12,:) = (/0.2,0.2,0.1/)
+
+		 xAlpha(13,:) = (/0.1,0.2,0.2/)
+		 xAlpha(14,:) = (/0.2,0.1,0.2/)
+		 xAlpha(15,:) = (/0.2,0.2,0.1/)
+		 xAlpha(16,:) = (/0.1,0.2,0.2/)
+		 xAlpha(17,:) = (/0.2,0.1,0.2/)
+		 xAlpha(18,:) = (/0.2,0.2,0.1/)
+
+		 xAlpha(19,:) = (/0.1,0.2,0.2/)
+		 xAlpha(20,:) = (/0.2,0.1,0.2/)
+		 xAlpha(21,:) = (/0.2,0.2,0.1/)
+		 xAlpha(22,:) = (/0.1,0.2,0.2/)
+		 xAlpha(23,:) = (/0.2,0.1,0.2/)
+		 xAlpha(24,:) = (/0.2,0.2,0.1/)
+
+		 xAlpha(25,:) = (/0.1,0.2,0.2/)
+		 xAlpha(26,:) = (/0.2,0.1,0.2/)
+		 xAlpha(27,:) = (/0.2,0.2,0.1/)
+		 xAlpha(28,:) = (/0.1,0.2,0.2/)
+		 xAlpha(29,:) = (/0.2,0.1,0.2/)
+		 xAlpha(30,:) = (/0.2,0.2,0.1/)