diff --git a/kMaterialParam.f b/kMaterialParam.f index a82d285..333fffe 100644 --- a/kMaterialParam.f +++ b/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,10 +112,12 @@ 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 + XTAUC5 = 2000.0 ! Pyramidal 2nd + ! thermal expansion coefficients alpha1 = 9.5D-6 alpha2 = alpha1 @@ -123,6 +125,42 @@ SUBROUTINE kMaterialParam(iphase,caratio,compliance,G12,thermat, ! 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 *reduce CRSS by 8 % for fitting Farrell et al. ! + XTAUC4 = XTAUC2 ! Pyramidal + XTAUC1 = 1.3333D0*XTAUC2 ! Basal + XTAUC3 = crssratio*XTAUC2 ! Pyramidal 1st + XTAUC5 = XTAUC3*5.00D0 ! Pyramidal 2nd + + ! 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 type slip systems + burgerv(13:nSys) = burger2 ! all types ! assign CRSS - tauc(1:3) = XTAUC1 - tauc(4:6) = XTAUC2 - tauc(7:12) = XTAUC4 + tauc(1:3) = XTAUC1 ! Basal + tauc(4:6) = XTAUC2 ! Prismatic + tauc(7:12) = XTAUC4 ! Pyramidal + tauc(13:24) = XTAUC3 ! Pyramidal 1st + tauc(25:30) = XTAUC5 ! Pyramidal 2nd case(1) !bcc diff --git a/kdirns.f b/kdirns.f index 8434b5c..43bb0f3 100644 --- a/kdirns.f +++ b/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' diff --git a/kmat.f b/kmat.f index 3df9ed7..89aabba 100644 --- a/kmat.f +++ b/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 diff --git a/kslip6ET.f b/kslip6ET.f index 65ee8b7..2a348e4 100644 --- a/kslip6ET.f +++ b/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 diff --git a/psiAnn.f b/psiAnn.f new file mode 100644 index 0000000..939a9d8 --- /dev/null +++ b/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/) diff --git a/umat.for b/umat.for index 48589e8..cace3f4 100644 --- a/umat.for +++ b/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' diff --git a/xAlpha.f b/xAlpha.f new file mode 100644 index 0000000..46564d4 --- /dev/null +++ b/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/) \ No newline at end of file diff --git a/xDir0CH.f b/xDir0CH.f new file mode 100644 index 0000000..ecabd31 --- /dev/null +++ b/xDir0CH.f @@ -0,0 +1,54 @@ +! The following are HCP slip system directions for a c/a ratio of 1 +! 3 basal + ddir(1,:) = (/1.000000, 0.000000, 0.000000/) + ddir(2,:) = (/-0.500000, -0.866025, 0.000000/) + ddir(3,:) = (/-0.500000, 0.866025, 0.000000/) + +! 3 prismatic + ddir(4,:) = (/-0.500000, 0.866025, 0.000000/) + ddir(5,:) = (/1.000000, 0.000000, 0.000000/) + ddir(6,:) = (/-0.500000, -0.866025, 0.000000/) + +! 6 slip directions for 1st order pyramidal planes + ddir(7,:) = (/ -1.000000, 0.000000, 0.000000/) + ddir(8,:) = (/-1.000000, 0.000000, 0.000000/) + ddir(9,:) = (/-0.500000, -0.866025, 0.000000/) + ddir(10,:) = (/-0.500000, -0.866025, 0.000000/) + ddir(11,:) = (/-0.500000, 0.866025, 0.000000/) + ddir(12,:) = (/-0.500000, 0.866025, 0.000000/) + +! 12 slip directions for 1st order pyramidal planes + ddir(13,:) = (/-0.707107, 0.000000, 0.707107/) + ddir(14,:) = (/-0.353553, 0.612372, 0.707107/) + ddir(15,:) = (/-0.353553, -0.612372, 0.707107/) + ddir(16,:) = (/-0.707107, 0.000000, 0.707107/) + ddir(17,:) = (/-0.353553, -0.612372, 0.707107/) + ddir(18,:) = (/0.353553, -0.612372, 0.707107/) + ddir(19,:) = (/0.707107, 0.000000, 0.707107/) + ddir(20,:) = (/0.353553, -0.612372, 0.707107/) + ddir(21,:) = (/0.707107, 0.000000, 0.707107/) + ddir(22,:) = (/0.353553, 0.612372, 0.707107/) + ddir(23,:) = (/0.353553, 0.612372, 0.707107/) + ddir(24,:) = (/-0.353553, 0.612372, 0.707107/) + + +! 6 slip directions for 2nd order pyramidal slip planes + ddir(25,:) = (/-0.353553, 0.612372, 0.707107/) + ddir(26,:) = (/-0.707107, 0.000000, 0.707107/) + ddir(27,:) = (/ -0.353553, -0.612372, 0.707107/) + ddir(28,:) =(/ 0.353553, -0.612372, 0.707107/) + ddir(29,:) =(/0.707107, 0.000000, 0.707107/) + ddir(30,:) =(/0.353553, 0.612372, 0.707107/) + + + + + + + + + + + + + diff --git a/xDir0ET.f b/xDir0ET.f deleted file mode 100644 index 18e5825..0000000 --- a/xDir0ET.f +++ /dev/null @@ -1,23 +0,0 @@ -! 3 basal - ddir(1,:) = (/-0.5000, 0.8660, 0.0000/) - ddir(2,:) = (/1.0000, 0.0000, 0.0000/) - ddir(3,:) = (/-0.5000, -0.8660, 0.0000/) -! 3 prismatic - ddir(4,:) = (/-0.5000, 0.8660, 0.0000/) - ddir(5,:) = (/1.0000, 0.0000, 0.0000/) - ddir(6,:) = (/-0.5000, -0.8660, 0.0000/) -! 6 slip direc for 1st order pyramidal planes - ! ddir(7,:) = (/ -0.5000, 0.8660, 0./) - !ddir(8,:) = (/-0.5000, 0.8660, 0./) - !ddir(9,:) = (/1.0000, 0., 0./) - !ddir(10,:) = (/1.0000, 0., 0./) - !ddir(11,:) = (/-0.5000, -0.8660, 0./) - !ddir(12,:) = (/-0.5000, -0.8660, 0./) -! 6 slip directions for 2nd order pyramidal slip planes - ddir(7,:) = (/-0.3536, 0.6124, 0.7071/) - ddir(8,:) = (/-0.3536, 0.6124, -0.7071/) - ddir(9,:) = (/ 0.7071, 0.0000, 0.7071/) - ddir(10,:) =(/ 0.7071, 0.0000, -0.7071/) - ddir(11,:) =(/-0.3536,-0.6124, 0.7071/) - ddir(12,:) =(/-0.3536,-0.6124, -0.7071/) - diff --git a/xNorm0CH.f b/xNorm0CH.f new file mode 100644 index 0000000..ad246ce --- /dev/null +++ b/xNorm0CH.f @@ -0,0 +1,40 @@ +! The following are HCP slip system normals for a c/a ratio of 1 +! basal plane for 3 slip + dnor(1,:) = (/0.000000, 0.000000, 1.000000/) + dnor(2,:) = (/0.000000, 0.000000, 1.000000/) + dnor(3,:) = (/0.000000, 0.000000, 1.000000/) + + !prismatic planes for 3 slip + dnor(4,:) = (/0.866025, 0.500000, 0.000000/) + dnor(5,:) = (/ 0.000000, 1.000000, 0.000000/) + dnor(6,:) = (/ 0.866025, -0.500000, 0.000000/) + +! 6 1st order pyramidal slip planes for slip + dnor(7,:) = (/0.000000, -0.755929, 0.654654/) + dnor(8,:) = (/ 0.000000, 0.755929, 0.654654/) + dnor(9,:) = (/ 0.654654, -0.377964, 0.654654/) + dnor(10,:) = (/ -0.654654, 0.377964, 0.654654/) + dnor(11,:) = (/-0.654654, -0.377964, 0.654654/) + dnor(12,:) = (/0.654654, 0.377964, 0.654654/) + +! 12 1st order pyramidal slip planes for slip + dnor(13,:) = (/0.654654, -0.377964, 0.654654/) + dnor(14,:) = (/0.654654, -0.377964, 0.654654/) + dnor(15,:) = (/0.654654, 0.377964, 0.654654/) + dnor(16,:) = (/0.654654, 0.377964, 0.654654/) + dnor(17,:) = (/0.000000, 0.755929, 0.654654/) + dnor(18,:) = (/0.000000, 0.755929, 0.654654/) + dnor(19,:) = (/-0.654654, 0.377964, 0.654654/) + dnor(20,:) = (/-0.654654, 0.377964, 0.654654/) + dnor(21,:) = (/-0.654654, -0.377964, 0.654654/) + dnor(22,:) = (/-0.654654, -0.377964, 0.654654/) + dnor(23,:) = (/0.000000, -0.755929, 0.654654/) + dnor(24,:) = (/0.000000, -0.755929, 0.654654/) + +! 6 2nd order pyramidalslip planes for slip + dnor(25,:) = (/0.353553, -0.612372, 0.707107/) + dnor(26,:) = (/0.707107, 0.000000, 0.707107/) + dnor(27,:) = (/0.353553, 0.612372, 0.707107/) + dnor(28,:) =(/-0.353553, 0.612372, 0.707107/) + dnor(29,:) =(/-0.707107, 0.000000, 0.707107/) + dnor(30,:) =(/-0.353553, -0.612372, 0.707107/) \ No newline at end of file diff --git a/xNorm0ET.f b/xNorm0ET.f deleted file mode 100644 index b012f75..0000000 --- a/xNorm0ET.f +++ /dev/null @@ -1,23 +0,0 @@ -! basal plane for 3 slip - dnor(1,:) = (/0.0000,0.0000,1.0000/) - dnor(2,:) = (/0.0000,0.0000,1.0000/) - dnor(3,:) = (/0.0000,0.0000,1.0000/) - !prismatic planes for 3 slip - dnor(4,:) = (/-0.8660, -0.5000, 0.0000/) - dnor(5,:) = (/ 0.0000, -1.0000, 0.0000/) - dnor(6,:) = (/ 0.8660, -0.5000, 0.0000/) -! 6 1st order pyramidal slip planes for slip - ! dnor(7,:) = (/-0.7500, -0.4330, 0.5000/) - !dnor(8,:) = (/ 0.7500, 0.4330, 0.5000/) - !dnor(9,:) = (/ 0., 0.8660, 0.5000/) - !dnor(10,:) = (/ 0., -0.8660, 0.5000/) - !dnor(11,:) = (/0.7500, -0.4330, 0.5000/) - !dnor(12,:) = (/-0.7500, 0.4330, 0.5000/) -! 6 2nd order pyramidalslip planes for slip - dnor(7,:) = (/ 0.3536, -0.6124, 0.7071/) - dnor(8,:) = (/-0.3536, 0.6124, 0.7071/) - dnor(9,:) = (/-0.7071, 0.0000, 0.7071/) - dnor(10,:) =(/ 0.7071, 0.0000, 0.7071/) - dnor(11,:) =(/ 0.3536, 0.6124, 0.7071/) - dnor(12,:) =(/-0.3536, -0.6124, 0.7071/) -