diff --git a/code/constitutive_dislotwin.f90 b/code/constitutive_dislotwin.f90 index 7c66b516c..54a5875c7 100644 --- a/code/constitutive_dislotwin.f90 +++ b/code/constitutive_dislotwin.f90 @@ -1139,14 +1139,14 @@ do f = 1_pInt,lattice_maxNslipFamily ! loop over StressRatio_pminus1*(1-StressRatio_p)**(constitutive_dislotwin_q(myInstance)-1.0_pReal) !* Plastic velocity gradient for dislocation glide - Lp = Lp + (1.0_pReal - sumf)*gdot_slip(j)*lattice_Sslip(:,:,index_myFamily+i,myStructure) + Lp = Lp + (1.0_pReal - sumf)*gdot_slip(j)*lattice_Sslip(:,:,1,index_myFamily+i,myStructure) !* Calculation of the tangent of Lp forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & dLp_dTstar3333(k,l,m,n) = & dLp_dTstar3333(k,l,m,n) + dgdot_dtauslip(j)*& - lattice_Sslip(k,l,index_myFamily+i,myStructure)*& - lattice_Sslip(m,n,index_myFamily+i,myStructure) + lattice_Sslip(k,l,1,index_myFamily+i,myStructure)*& + lattice_Sslip(m,n,1,index_myFamily+i,myStructure) enddo enddo diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index d8c3055b0..f7ae2d22a 100644 --- a/code/constitutive_nonlocal.f90 +++ b/code/constitutive_nonlocal.f90 @@ -183,7 +183,8 @@ rhoDotFluxOutput, & rhoDotMultiplicationOutput, & rhoDotSingle2DipoleGlideOutput, & rhoDotAthermalAnnihilationOutput, & -rhoDotThermalAnnihilationOutput +rhoDotThermalAnnihilationOutput, & +screwStressProjection !< combined projection of Schmid and non-Schmid stresses for resolved shear stress acting on screws real(pReal), dimension(:,:,:,:,:,:), allocatable, private :: & compatibility !< slip system compatibility between me and my neighbors @@ -414,7 +415,7 @@ allocate(peierlsStressPerSlipFamily(lattice_maxNslipFamily,2,maxNinstance)) minDipoleHeightPerSlipFamily = -1.0_pReal peierlsStressPerSlipFamily = 0.0_pReal -allocate(nonSchmidCoeff(lattice_maxNonSchmid,maxNinstance)) +allocate(nonSchmidCoeff(lattice_maxNnonSchmid,maxNinstance)) nonSchmidCoeff = 0.0_pReal !*** readout data from material.config file @@ -577,7 +578,7 @@ do while (trim(line) /= '#EOF#') case('shortrangestresscorrection') shortRangeStressCorrection(i) = IO_floatValue(line,positions,2_pInt) > 0.0_pReal case ('nonschmid_coefficients') - do f = 1_pInt, lattice_maxNonSchmid + do f = 1_pInt, lattice_maxNnonSchmid nonSchmidCoeff(f,i) = IO_floatValue(line,positions,1_pInt+f) enddo case('deadzonescaling','deadzone','deadscaling') @@ -769,6 +770,10 @@ peierlsStress = 0.0_pReal allocate(colinearSystem(maxTotalNslip,maxNinstance)) colinearSystem = 0_pInt +allocate(screwStressProjection(3,3,4,maxTotalNslip,maxNinstance)) +screwStressProjection = 0.0_pReal + + do i = 1,maxNinstance myStructure = constitutive_nonlocal_structure(i) ! lattice structure of this instance @@ -1003,6 +1008,29 @@ do i = 1,maxNinstance -lattice_st(1:3, slipSystemLattice(s1,i), myStructure), & lattice_sn(1:3, slipSystemLattice(s1,i), myStructure)], [3,3])) enddo + + + !*** combined projection of Schmid and non-Schmid stress contributions to resolved shear stress for screws + !* four types t: + !* 1) positive screw at positive resolved stress + !* 2) positive screw at negative resolved stress + !* 3) negative screw at positive resolved stress + !* 4) negative screw at negative resolved stress + + screwStressProjection = 0.0_pReal + do s = 1_pInt,ns + do l = 1_pInt,lattice_NnonSchmid(myStructure) + screwStressProjection(1:3,1:3,1,s,i) = screwStressProjection(1:3,1:3,1,s,i) & + + nonSchmidCoeff(l,i) * lattice_Sslip(1:3,1:3,2*l,slipSystemLattice(s,i),myStructure) + screwStressProjection(1:3,1:3,2,s,i) = screwStressProjection(1:3,1:3,2,s,i) & + + nonSchmidCoeff(l,i) * lattice_Sslip(1:3,1:3,2*l+1,slipSystemLattice(s,i),myStructure) + enddo + screwStressProjection(1:3,1:3,3,s,i) = -screwStressProjection(1:3,1:3,2,s,i) + screwStressProjection(1:3,1:3,4,s,i) = -screwStressProjection(1:3,1:3,1,s,i) + forall (t = 1:4) & + screwStressProjection(1:3,1:3,t,s,i) = screwStressProjection(1:3,1:3,t,s,i) & + + lattice_Sslip(1:3,1:3,1,slipSystemLattice(s,i),myStructure) + enddo enddo @@ -1708,6 +1736,7 @@ subroutine constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar99, Tstar_v, Temp use math, only: math_Plain3333to99, & math_mul6x6, & + math_mul33xx33, & math_Mandel6to33 use debug, only: debug_level, & debug_constitutive, & @@ -1722,7 +1751,7 @@ use material, only: homogenization_maxNgrains, & phase_plasticityInstance use lattice, only: lattice_Sslip, & lattice_Sslip_v, & - NnonSchmid + lattice_NnonSchmid use mesh, only: mesh_ipVolume implicit none @@ -1754,8 +1783,6 @@ integer(pInt) myInstance, & ! curren s, & ! index of my current slip system sLattice ! index of my current slip system according to lattice order real(pReal), dimension(3,3,3,3) :: dLp_dTstar3333 ! derivative of Lp with respect to Tstar (3x3x3x3 matrix) -real(pReal), dimension(3,3,2,totalNslip(phase_plasticityInstance(material_phase(g,ip,el)))) :: & - nonSchmidTensor real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),8) :: & rhoSgl ! single dislocation densities (including blocked) real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),4) :: & @@ -1773,7 +1800,6 @@ real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,e Lp = 0.0_pReal dLp_dTstar3333 = 0.0_pReal -nonSchmidTensor = 0.0_pReal myInstance = phase_plasticityInstance(material_phase(g,ip,el)) myStructure = constitutive_nonlocal_structure(myInstance) @@ -1795,31 +1821,26 @@ tauBack = state%p(iTauB(1:ns,myInstance)) !*** get effective resolved shear stress -!*** add non schmid contributions to ONLY screw components if present (i.e. if NnonSchmid(myStructure) > 0) +!*** for screws possible non-schmid contributions are also taken into account do s = 1_pInt,ns - sLattice = slipSystemLattice(s,myInstance) - tau(s,1:4) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,myStructure)) + tauBack(s) - nonSchmidTensor(1:3,1:3,1,s) = lattice_Sslip(1:3,1:3,sLattice,myStructure) - nonSchmidTensor(1:3,1:3,2,s) = nonSchmidTensor(1:3,1:3,1,s) - do k = 1_pInt, NnonSchmid(myStructure) - tau(s,3) = tau(s,3) + nonSchmidCoeff(k,myInstance) & - * math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,2*k,sLattice,myStructure)) - tau(s,4) = tau(s,4) + nonSchmidCoeff(k,myInstance) & - * math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,2*k+1,sLattice,myStructure)) - nonSchmidTensor(1:3,1:3,1,s) = nonSchmidTensor(1:3,1:3,1,s) & - + nonSchmidCoeff(k,myInstance) & - * math_Mandel6to33(lattice_Sslip_v(1:6,2*k,sLattice,myStructure)) - nonSchmidTensor(1:3,1:3,2,s) = nonSchmidTensor(1:3,1:3,2,s) & - + nonSchmidCoeff(k,myInstance) & - * math_Mandel6to33(lattice_Sslip_v(1:6,2*k+1,sLattice,myStructure)) - enddo + sLattice = slipSystemLattice(s,myInstance) + tau(s,1:2) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,myStructure)) + if (tau(s,1) > 0.0_pReal) then + tau(s,3) = math_mul33xx33(math_Mandel6to33(Tstar_v), screwStressProjection(1:3,1:3,1,s,myInstance)) + tau(s,4) = math_mul33xx33(math_Mandel6to33(Tstar_v), screwStressProjection(1:3,1:3,3,s,myInstance)) + else + tau(s,3) = math_mul33xx33(math_Mandel6to33(Tstar_v), screwStressProjection(1:3,1:3,2,s,myInstance)) + tau(s,4) = math_mul33xx33(math_Mandel6to33(Tstar_v), screwStressProjection(1:3,1:3,4,s,myInstance)) + endif + forall (t = 1_pInt:4_pInt) & + tau(s,t) = tau(s,t) + tauBack(s) ! add backstress enddo !*** get dislocation velocity and its tangent and store the velocity in the state array -if (myStructure == 1_pInt .and. NnonSchmid(myStructure) == 0_pInt) then ! for fcc all velcities are equal +if (myStructure == 1_pInt .and. lattice_NnonSchmid(myStructure) == 0_pInt) then ! for fcc all velcities are equal call constitutive_nonlocal_kinetics(v(1:ns,1), tau(1:ns,1), 1_pInt, Temperature, state, & g, ip, el, dv_dtau(1:ns,1)) do t = 1_pInt,4_pInt @@ -1860,13 +1881,22 @@ enddo do s = 1_pInt,ns sLattice = slipSystemLattice(s,myInstance) - Lp = Lp + gdotTotal(s) * lattice_Sslip(1:3,1:3,sLattice,myStructure) - forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt,k=1_pInt:3_pInt,l=1_pInt:3_pInt) & - dLp_dTstar3333(i,j,k,l) = dLp_dTstar3333(i,j,k,l) & - + dgdot_dtau(s,1) * lattice_Sslip(i,j,sLattice,myStructure) * lattice_Sslip(k,l,sLattice,myStructure) & - + dgdot_dtau(s,2) * lattice_Sslip(i,j,sLattice,myStructure) * lattice_Sslip(k,l,sLattice,myStructure) & - + dgdot_dtau(s,3) * lattice_Sslip(i,j,sLattice,myStructure) * nonSchmidTensor(k,l,1,s) & - + dgdot_dtau(s,4) * lattice_Sslip(i,j,sLattice,myStructure) * nonSchmidTensor(k,l,2,s) + Lp = Lp + gdotTotal(s) * lattice_Sslip(1:3,1:3,1,sLattice,myStructure) + if (tau(s,1) > 0.0_pReal) then + forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt,k=1_pInt:3_pInt,l=1_pInt:3_pInt) & + dLp_dTstar3333(i,j,k,l) = dLp_dTstar3333(i,j,k,l) & + + dgdot_dtau(s,1) * lattice_Sslip(i,j,1,sLattice,myStructure) * lattice_Sslip(k,l,1,sLattice,myStructure) & + + dgdot_dtau(s,2) * lattice_Sslip(i,j,1,sLattice,myStructure) * lattice_Sslip(k,l,1,sLattice,myStructure) & + + dgdot_dtau(s,3) * lattice_Sslip(i,j,1,sLattice,myStructure) * screwStressProjection(k,l,1,s,myInstance) & + + dgdot_dtau(s,4) * lattice_Sslip(i,j,1,sLattice,myStructure) * screwStressProjection(k,l,3,s,myInstance) + else + forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt,k=1_pInt:3_pInt,l=1_pInt:3_pInt) & + dLp_dTstar3333(i,j,k,l) = dLp_dTstar3333(i,j,k,l) & + + dgdot_dtau(s,1) * lattice_Sslip(i,j,1,sLattice,myStructure) * lattice_Sslip(k,l,1,sLattice,myStructure) & + + dgdot_dtau(s,2) * lattice_Sslip(i,j,1,sLattice,myStructure) * lattice_Sslip(k,l,1,sLattice,myStructure) & + + dgdot_dtau(s,3) * lattice_Sslip(i,j,1,sLattice,myStructure) * screwStressProjection(k,l,2,s,myInstance) & + + dgdot_dtau(s,4) * lattice_Sslip(i,j,1,sLattice,myStructure) * screwStressProjection(k,l,4,s,myInstance) + endif enddo dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333) diff --git a/code/constitutive_phenopowerlaw.f90 b/code/constitutive_phenopowerlaw.f90 index f6a68db38..1f42b3097 100644 --- a/code/constitutive_phenopowerlaw.f90 +++ b/code/constitutive_phenopowerlaw.f90 @@ -240,7 +240,7 @@ subroutine constitutive_phenopowerlaw_init(myFile) constitutive_phenopowerlaw_aTolShear = 0.0_pReal allocate(constitutive_phenopowerlaw_aTolTwinfrac(maxNinstance)) constitutive_phenopowerlaw_aTolTwinfrac = 0.0_pReal - allocate(constitutive_phenopowerlaw_nonSchmidCoeff(lattice_maxNonSchmid,maxNinstance)) + allocate(constitutive_phenopowerlaw_nonSchmidCoeff(lattice_maxNnonSchmid,maxNinstance)) constitutive_phenopowerlaw_nonSchmidCoeff = 0.0_pReal rewind(myFile) @@ -368,7 +368,7 @@ subroutine constitutive_phenopowerlaw_init(myFile) constitutive_phenopowerlaw_interaction_TwinTwin(j,i) = IO_floatValue(line,positions,1_pInt+j) enddo case ('nonschmid_coefficients') - do j = 1_pInt, lattice_maxNonSchmid + do j = 1_pInt, lattice_maxNnonSchmid constitutive_phenopowerlaw_nonSchmidCoeff(j,i) = IO_floatValue(line,positions,1_pInt+j) enddo case default @@ -682,7 +682,7 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar lattice_maxNtwinFamily, & lattice_NslipSystem, & lattice_NtwinSystem, & - NnonSchmid + lattice_NnonSchmid use mesh, only: & mesh_NcpElems, & mesh_maxNips @@ -745,17 +745,17 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar ! Calculation of Lp tau_slip_pos(j) = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,structID)) tau_slip_neg(j) = tau_slip_pos(j) - nonSchmid_tensor(1:3,1:3,1) = math_Mandel6to33(lattice_Sslip_v(1:6,1,index_myFamily+i,structID)) - nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,1) - do k = 1, NnonSchmid(structID) + nonSchmid_tensor(1:3,1:3,1) = lattice_Sslip(1:3,1:3,1,index_myFamily+i,structID) + nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,1) + do k = 1,lattice_NnonSchmid(structID) tau_slip_pos(j) = tau_slip_pos(j) + constitutive_phenopowerlaw_nonSchmidCoeff(k,matID)* & dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k,index_myFamily+i,structID)) tau_slip_neg(j) = tau_slip_neg(j) + constitutive_phenopowerlaw_nonSchmidCoeff(k,matID)* & dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,structID)) nonSchmid_tensor(1:3,1:3,1) = nonSchmid_tensor(1:3,1:3,1) + constitutive_phenopowerlaw_nonSchmidCoeff(k,matID)*& - math_Mandel6to33(lattice_Sslip_v(1:6,2*k,index_myFamily+i,structID)) + lattice_Sslip(1:3,1:3,2*k,index_myFamily+i,structID) nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,2) + constitutive_phenopowerlaw_nonSchmidCoeff(k,matID)*& - math_Mandel6to33(lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,structID)) + lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+i,structID) enddo gdot_slip_pos(j) = 0.5_pReal*constitutive_phenopowerlaw_gdot0_slip(matID)* & ((abs(tau_slip_pos(j))/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(matID))*& @@ -764,7 +764,7 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar ((abs(tau_slip_neg(j))/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(matID))*& sign(1.0_pReal,tau_slip_neg(j)) Lp = Lp + (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F - (gdot_slip_pos(j)+gdot_slip_neg(j))*lattice_Sslip(1:3,1:3,index_myFamily+i,structID) + (gdot_slip_pos(j)+gdot_slip_neg(j))*lattice_Sslip(1:3,1:3,1,index_myFamily+i,structID) !-------------------------------------------------------------------------------------------------- ! Calculation of the tangent of Lp @@ -772,7 +772,7 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar dgdot_dtauslip_pos(j) = gdot_slip_pos(j)*constitutive_phenopowerlaw_n_slip(matID)/tau_slip_pos(j) forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + & - dgdot_dtauslip_pos(j)*lattice_Sslip(k,l,index_myFamily+i,structID)* & + dgdot_dtauslip_pos(j)*lattice_Sslip(k,l,1,index_myFamily+i,structID)* & nonSchmid_tensor(m,n,1) endif @@ -780,7 +780,7 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar dgdot_dtauslip_neg(j) = gdot_slip_neg(j)*constitutive_phenopowerlaw_n_slip(matID)/tau_slip_neg(j) forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + & - dgdot_dtauslip_neg(j)*lattice_Sslip(k,l,index_myFamily+i,structID)* & + dgdot_dtauslip_neg(j)*lattice_Sslip(k,l,1,index_myFamily+i,structID)* & nonSchmid_tensor(m,n,2) endif enddo @@ -832,7 +832,7 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,temperature,state,ipc,ip,el lattice_NslipSystem, & lattice_NtwinSystem, & lattice_shearTwin, & - NnonSchmid + lattice_NnonSchmid use mesh, only: & mesh_NcpElems,& mesh_maxNips @@ -911,7 +911,7 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,temperature,state,ipc,ip,el ! Calculation of dot gamma tau_slip_pos(j) = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,structID)) tau_slip_neg(j) = tau_slip_pos(j) - do k = 1, NnonSchmid(structID) + do k = 1,lattice_NnonSchmid(structID) tau_slip_pos(j) = tau_slip_pos(j) + constitutive_phenopowerlaw_nonSchmidCoeff(k,matID)* & dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k,index_myFamily+i,structID)) tau_slip_neg(j) = tau_slip_neg(j) + constitutive_phenopowerlaw_nonSchmidCoeff(k,matID)* & @@ -1069,7 +1069,8 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,temperature,dt,stat lattice_maxNslipFamily, & lattice_maxNtwinFamily, & lattice_NslipSystem, & - lattice_NtwinSystem,NnonSchmid + lattice_NtwinSystem, & + lattice_NnonSchmid use mesh, only: & mesh_NcpElems, & mesh_maxNips @@ -1132,7 +1133,7 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,temperature,dt,stat j = j + 1_pInt tau_slip_pos = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,structID)) tau_slip_neg = tau_slip_pos - do k = 1, NnonSchmid(structID) + do k = 1,lattice_NnonSchmid(structID) tau_slip_pos = tau_slip_pos + constitutive_phenopowerlaw_nonSchmidCoeff(k,matID)* & dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k,index_myFamily+i,structID)) tau_slip_neg = tau_slip_neg + constitutive_phenopowerlaw_nonSchmidCoeff(k,matID)* & diff --git a/code/constitutive_titanmod.f90 b/code/constitutive_titanmod.f90 index 1dd5d8e14..6e4f3670f 100644 --- a/code/constitutive_titanmod.f90 +++ b/code/constitutive_titanmod.f90 @@ -1246,8 +1246,7 @@ subroutine constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,& lattice_maxNslipFamily, & lattice_maxNtwinFamily, & lattice_NslipSystem, & - lattice_NtwinSystem, & - NnonSchmid + lattice_NtwinSystem use mesh, only: & mesh_NcpElems, & mesh_maxNips @@ -1441,14 +1440,14 @@ subroutine constitutive_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,& !************************************************* !sumf=0.0_pReal !* Plastic velocity gradient for dislocation glide - Lp = Lp + (1.0_pReal - sumf)*gdot_slip(j)*lattice_Sslip(:,:,index_myFamily+i,myStructure) + Lp = Lp + (1.0_pReal - sumf)*gdot_slip(j)*lattice_Sslip(1:3,1:3,1,index_myFamily+i,myStructure) !* Calculation of the tangent of Lp forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & dLp_dTstar3333(k,l,m,n) = & dLp_dTstar3333(k,l,m,n) + dgdot_dtauslip(j)*& - lattice_Sslip(k,l,index_myFamily+i,myStructure)*& - lattice_Sslip(m,n,index_myFamily+i,myStructure) + lattice_Sslip(k,l,1,index_myFamily+i,myStructure)*& + lattice_Sslip(m,n,1,index_myFamily+i,myStructure) enddo enddo diff --git a/code/lattice.f90 b/code/lattice.f90 index aca62bdcd..38f6f3e94 100644 --- a/code/lattice.f90 +++ b/code/lattice.f90 @@ -38,7 +38,7 @@ module lattice lattice_maxNslip = 33_pInt, & !< max # of slip systems over lattice structures lattice_maxNtwin = 24_pInt, & !< max # of twin systems over lattice structures lattice_maxNinteraction = 42_pInt, & !< max # of interaction types (in hardening matrix part) - lattice_maxNonSchmid = 6_pInt !< max # of non schmid contributions over lattice structures + lattice_maxNnonSchmid = 6_pInt !< max # of non schmid contributions over lattice structures integer(pInt), allocatable, dimension(:,:), protected, public :: & lattice_NslipSystem, & !< total # of slip systems in each family @@ -51,14 +51,16 @@ module lattice lattice_interactionTwinTwin !< Twin--twin interaction type + real(pReal), allocatable, dimension(:,:,:,:,:), protected, public :: & + lattice_Sslip !< Schmid and non-Schmid matrices + real(pReal), allocatable, dimension(:,:,:,:), protected, public :: & - lattice_Sslip_v, & - lattice_Sslip !< Schmid matrices, normal, shear direction and d x n of slip systems + lattice_Sslip_v !< Mandel notation of lattice_Sslip real(pReal), allocatable, dimension(:,:,:), protected, public :: & - lattice_sn, & - lattice_sd, & - lattice_st + lattice_sn, & !< normal direction of slip system + lattice_sd, & !< slip direction of slip system + lattice_st !< sd x sn ! rotation and Schmid matrices, normal, shear direction and d x n of twin systems real(pReal), allocatable, dimension(:,:,:,:), protected, public :: & @@ -85,7 +87,7 @@ module lattice interactionTwinTwin integer(pInt), allocatable, dimension(:), protected, public :: & - NnonSchmid !< total # of non-Schmid contributions for each structure + lattice_NnonSchmid !< total # of non-Schmid contributions for each structure !-------------------------------------------------------------------------------------------------- ! fcc (1) @@ -97,7 +99,8 @@ module lattice integer(pInt), parameter, private :: & lattice_fcc_Nslip = 12_pInt, & ! sum(lattice_fcc_NslipSystem), & !< total # of slip systems for fcc - lattice_fcc_Ntwin = 12_pInt ! sum(lattice_fcc_NtwinSystem) !< total # of twin systems for fcc + lattice_fcc_Ntwin = 12_pInt, & ! sum(lattice_fcc_NtwinSystem) !< total # of twin systems for fcc + lattice_fcc_NnonSchmid = 0_pInt !< total # of non-Schmid contributions for fcc integer(pInt), private :: & lattice_fcc_Nstructure = 0_pInt @@ -224,10 +227,6 @@ module lattice 2,2,2,2,2,2,2,2,2,1,1,1 & ],pInt),[lattice_fcc_Ntwin,lattice_fcc_Ntwin],order=[2,1]) !< Twin--twin interaction types for fcc - integer(pInt), parameter, private :: NnonSchmid_fcc = 0_pInt !< total # of non-Schmid contributions for fcc - - real(pReal), dimension(3,3,2,NnonSchmid_fcc,lattice_fcc_Nslip), parameter, private :: & - lattice_nonSchmid_fcc = 0.0_pReal ! reshape([],[3,3,2,NnonSchmid_fcc,lattice_fcc_Nslip]) !< Tensor for each non-Schmid contribution for fcc !-------------------------------------------------------------------------------------------------- @@ -240,7 +239,8 @@ module lattice integer(pInt), parameter, private :: & lattice_bcc_Nslip = 24_pInt, & ! sum(lattice_bcc_NslipSystem), & !< total # of slip systems for bcc - lattice_bcc_Ntwin = 12_pInt ! sum(lattice_bcc_NtwinSystem) !< total # of twin systems for bcc + lattice_bcc_Ntwin = 12_pInt, & ! sum(lattice_bcc_NtwinSystem) !< total # of twin systems for bcc + lattice_bcc_NnonSchmid = 6_pInt !< # of non-Schmid contributions for bcc. 6 known non schmid contributions for BCC (A. Koester, A. Ma, A. Hartmaier 2012) integer(pInt), private :: & lattice_bcc_Nstructure = 0_pInt @@ -419,10 +419,7 @@ module lattice !< 1: self interaction !< 2: collinear interaction !< 3: other interaction - integer(pInt), parameter, private :: NnonSchmid_bcc = 0_pInt !< # of non-Schmid contributions for bcc - real(pReal), dimension(3,3,2,NnonSchmid_bcc,lattice_bcc_Nslip), parameter, private :: & - lattice_nonSchmid_bcc = 0.0_pReal ! reshape([],[3,3,2,NnonSchmid_bcc,lattice_bcc_Nslip]) !< Tensor for each non-Schmid contribution for bcc !-------------------------------------------------------------------------------------------------- ! hex (3+) @@ -434,7 +431,8 @@ module lattice integer(pInt), parameter , private :: & lattice_hex_Nslip = 33_pInt, & ! sum(lattice_hex_NslipSystem), !< total # of slip systems for hex - lattice_hex_Ntwin = 24_pInt ! sum(lattice_hex_NtwinSystem) !< total # of twin systems for hex + lattice_hex_Ntwin = 24_pInt, & ! sum(lattice_hex_NtwinSystem) !< total # of twin systems for hex + lattice_hex_NnonSchmid = 0_pInt !< # of non-Schmid contributions for hex integer(pInt), private :: & lattice_hex_Nstructure = 0_pInt @@ -689,10 +687,7 @@ module lattice 20,20,20,20,20,20, 19,19,19,19,19,19, 18,18,18,18,18,18, 17,17,17,17,17, 4 & ],pInt),[lattice_hex_Ntwin,lattice_hex_Ntwin],order=[2,1]) !< Twin--slip interaction types for hex (isotropic, 16 in total) - integer(pInt), parameter, private :: NnonSchmid_hex = 0_pInt !< # of non-Schmid contributions for hex - real(pReal), dimension(3,3,2,NnonSchmid_hex,lattice_hex_Nslip), parameter, private :: & - lattice_nonSchmid_hex = 0.0_pReal ! reshape([],[3,3,2,NnonSchmid_hex,lattice_hex_Nslip]) !< Tensor for each non-Schmid contribution for hex public :: & lattice_init, & @@ -745,9 +740,9 @@ subroutine lattice_init write(6,'(a16,1x,i5,/)') ' # structures:',lattice_Nstructure endif - allocate(NnonSchmid(lattice_Nstructure)); NnonSchmid = 0_pInt - allocate(lattice_Sslip(3,3,lattice_maxNslip,lattice_Nstructure)); lattice_Sslip = 0.0_pReal - allocate(lattice_Sslip_v(6,1+2*lattice_maxNonSchmid,lattice_maxNslip,lattice_Nstructure)); lattice_Sslip_v = 0.0_pReal + allocate(lattice_NnonSchmid(lattice_Nstructure)); lattice_NnonSchmid = 0_pInt + allocate(lattice_Sslip(3,3,1+2*lattice_maxNnonSchmid,lattice_maxNslip,lattice_Nstructure)); lattice_Sslip = 0.0_pReal + allocate(lattice_Sslip_v(6,1+2*lattice_maxNnonSchmid,lattice_maxNslip,lattice_Nstructure)); lattice_Sslip_v = 0.0_pReal allocate(lattice_sd(3,lattice_maxNslip,lattice_Nstructure)); lattice_sd = 0.0_pReal allocate(lattice_st(3,lattice_maxNslip,lattice_Nstructure)); lattice_st = 0.0_pReal allocate(lattice_sn(3,lattice_maxNslip,lattice_Nstructure)); lattice_sn = 0.0_pReal @@ -784,6 +779,7 @@ integer(pInt) function lattice_initializeStructure(struct,CoverA) math_vectorproduct, & math_tensorproduct, & math_norm3, & + math_mul33x3, & math_trace33, & math_symmetric33, & math_Mandel33to6, & @@ -795,9 +791,13 @@ integer(pInt) function lattice_initializeStructure(struct,CoverA) implicit none character(len=*) struct real(pReal) CoverA + real(pReal), dimension(3) :: sdU = 0.0_pReal, & + snU = 0.0_pReal, & + np = 0.0_pReal, & + nn = 0.0_pReal real(pReal), dimension(3,lattice_maxNslip) :: sd = 0.0_pReal, & sn = 0.0_pReal - real(pReal), dimension(12,lattice_maxNonSchmid,lattice_maxNslip) :: sns = 0.0_pReal + real(pReal), dimension(3,3,2,lattice_maxNnonSchmid,lattice_maxNslip) :: sns = 0.0_pReal real(pReal), dimension(3,lattice_maxNtwin) :: td = 0.0_pReal, & tn = 0.0_pReal real(pReal), dimension(lattice_maxNtwin) :: ts = 0.0_pReal @@ -818,13 +818,13 @@ integer(pInt) function lattice_initializeStructure(struct,CoverA) lattice_fcc_Nstructure = lattice_fcc_Nstructure + 1_pInt ! count fcc instances if (lattice_fcc_Nstructure == 1_pInt) then ! me is first fcc structure processMe = .true. - NnonSchmid(myStructure) = NnonSchmid_fcc ! Currently no known non schmid contributions for FCC (to be changed later) + lattice_NnonSchmid(myStructure) = lattice_fcc_NnonSchmid ! Currently no known non schmid contributions for FCC (to be changed later) do i = 1_pInt,myNslip ! assign slip system vectors sd(1:3,i) = lattice_fcc_systemSlip(1:3,i) sn(1:3,i) = lattice_fcc_systemSlip(4:6,i) - do j = 1_pInt, NnonSchmid_fcc - sns(1:6,j,i) = math_Mandel33to6(lattice_nonSchmid_fcc(1:3,1:3,1,j,i)) - sns(7:12,j,i) = math_Mandel33to6(lattice_nonSchmid_fcc(1:3,1:3,2,j,i)) + do j = 1_pInt,lattice_fcc_NnonSchmid + sns(1:3,1:3,1,j,i) = 0.0_pReal + sns(1:3,1:3,2,j,i) = 0.0_pReal enddo enddo do i = 1_pInt,myNtwin ! assign twin system vectors and shears @@ -847,14 +847,26 @@ integer(pInt) function lattice_initializeStructure(struct,CoverA) lattice_bcc_Nstructure = lattice_bcc_Nstructure + 1_pInt ! count bcc instances if (lattice_bcc_Nstructure == 1_pInt) then ! me is first bcc structure processMe = .true. - NnonSchmid(myStructure) = NnonSchmid_BCC ! 5 known non schmid contributions for BCC (A. Koester, A. Ma, A. Hartmaier 2012) + lattice_NnonSchmid(myStructure) = lattice_bcc_NnonSchmid do i = 1_pInt,myNslip ! assign slip system vectors sd(1:3,i) = lattice_bcc_systemSlip(1:3,i) sn(1:3,i) = lattice_bcc_systemSlip(4:6,i) - do j = 1_pInt, NnonSchmid_bcc - sns(1:6,j,i) = math_Mandel33to6(lattice_nonSchmid_bcc(1:3,1:3,1,j,i)) - sns(7:12,j,i) = math_Mandel33to6(lattice_nonSchmid_bcc(1:3,1:3,2,j,i)) - enddo + sdU = sd(1:3,i) / math_norm3(sd(1:3,i)) + snU = sn(1:3,i) / math_norm3(sn(1:3,i)) + np = math_mul33x3(math_axisAngleToR(sdU,60.0_pReal*INRAD), snU) + nn = math_mul33x3(math_axisAngleToR(-sdU,60.0_pReal*INRAD), snU) + sns(1:3,1:3,1,1,i) = math_tensorproduct(sdU, np) + sns(1:3,1:3,2,1,i) = math_tensorproduct(-sdU, nn) + sns(1:3,1:3,1,2,i) = math_tensorproduct(math_vectorproduct(snU, sdU), snU) + sns(1:3,1:3,2,2,i) = math_tensorproduct(math_vectorproduct(snU, -sdU), snU) + sns(1:3,1:3,1,3,i) = math_tensorproduct(math_vectorproduct(np, sdU), np) + sns(1:3,1:3,2,3,i) = math_tensorproduct(math_vectorproduct(nn, -sdU), nn) + sns(1:3,1:3,1,4,i) = math_tensorproduct(snU, snU) + sns(1:3,1:3,2,4,i) = math_tensorproduct(snU, snU) + sns(1:3,1:3,1,5,i) = math_tensorproduct(math_vectorproduct(snU, sdU), math_vectorproduct(snU, sdU)) + sns(1:3,1:3,2,5,i) = math_tensorproduct(math_vectorproduct(snU, -sdU), math_vectorproduct(snU, -sdU)) + sns(1:3,1:3,1,6,i) = math_tensorproduct(sdU, sdU) + sns(1:3,1:3,2,6,i) = math_tensorproduct(-sdU, -sdU) enddo do i = 1_pInt,myNtwin ! assign twin system vectors and shears td(1:3,i) = lattice_bcc_systemTwin(1:3,i) @@ -876,7 +888,7 @@ integer(pInt) function lattice_initializeStructure(struct,CoverA) myNslip = lattice_hex_Nslip ! overall number of slip systems myNtwin = lattice_hex_Ntwin ! overall number of twin systems processMe = .true. - NnonSchmid(myStructure) = NnonSchmid_hex ! Currently no known non schmid contributions for hex (to be changed later) + lattice_NnonSchmid(myStructure) = lattice_hex_NnonSchmid ! Currently no known non schmid contributions for hex (to be changed later) ! converting from 4 axes coordinate system (a1=a2=a3=c) to ortho-hexgonal system (a, b, c) do i = 1_pInt,myNslip sd(1,i) = lattice_hex_systemSlip(1,i)*1.5_pReal ! direction [uvtw]->[3u/2 (u+2v)*sqrt(3)/2 w*(c/a)] @@ -885,9 +897,9 @@ integer(pInt) function lattice_initializeStructure(struct,CoverA) sn(1,i) = lattice_hex_systemSlip(5,i) ! plane (hkil)->(h (h+2k)/sqrt(3) l/(c/a)) sn(2,i) = (lattice_hex_systemSlip(5,i)+2.0_pReal*lattice_hex_systemSlip(6,i))/sqrt(3.0_pReal) sn(3,i) = lattice_hex_systemSlip(8,i)/CoverA - do j = 1_pInt, NnonSchmid_hex - sns(1:6,j,i) = math_Mandel33to6(lattice_nonSchmid_hex(1:3,1:3,1,j,i)) - sns(7:12,j,i) = math_Mandel33to6(lattice_nonSchmid_hex(1:3,1:3,2,j,i)) + do j = 1_pInt,lattice_hex_NnonSchmid + sns(1:3,1:3,1,j,i) = 0.0_pReal + sns(1:3,1:3,2,j,i) = 0.0_pReal enddo enddo do i = 1_pInt,myNtwin @@ -925,14 +937,17 @@ integer(pInt) function lattice_initializeStructure(struct,CoverA) lattice_sn(1:3,i,myStructure) = sn(1:3,i)/math_norm3(sn(1:3,i)) ! make unit vector lattice_st(1:3,i,myStructure) = math_vectorproduct(lattice_sd(1:3,i,myStructure), & lattice_sn(1:3,i,myStructure)) - lattice_Sslip(1:3,1:3,i,myStructure) = math_tensorproduct(lattice_sd(1:3,i,myStructure), & - lattice_sn(1:3,i,myStructure)) - lattice_Sslip_v(1:6,1,i,myStructure) = math_Mandel33to6(math_symmetric33(lattice_Sslip(1:3,1:3,i,myStructure))) - do j = 1_pInt, NnonSchmid(myStructure) - lattice_Sslip_v(1:6,2*j,i,myStructure) = sns(1:6,j,i) - lattice_Sslip_v(1:6,2*j+1,i,myStructure) = sns(7:12,j,i) + lattice_Sslip(1:3,1:3,1,i,myStructure) = math_tensorproduct(lattice_sd(1:3,i,myStructure), & + lattice_sn(1:3,i,myStructure)) + do j = 1_pInt,lattice_NnonSchmid(myStructure) + lattice_Sslip(1:3,1:3,2*j ,i,myStructure) = sns(1:3,1:3,1,j,i) + lattice_Sslip(1:3,1:3,2*j+1,i,myStructure) = sns(1:3,1:3,2,j,i) enddo - if (abs(math_trace33(lattice_Sslip(1:3,1:3,i,myStructure))) > 1.0e-8_pReal) & + do j = 1_pInt,1_pInt+2_pInt*lattice_NnonSchmid(myStructure) + lattice_Sslip_v(1:6,j,i,myStructure) = & + math_Mandel33to6(math_symmetric33(lattice_Sslip(1:3,1:3,j,i,myStructure))) + enddo + if (abs(math_trace33(lattice_Sslip(1:3,1:3,1,i,myStructure))) > 1.0e-8_pReal) & call IO_error(0_pInt,myStructure,i,0_pInt,ext_msg = 'dilatational slip Schmid matrix') enddo do i = 1_pInt,myNtwin ! store twin system vectors and Schmid plus rotation matrix for my structure