diff --git a/code/constitutive.f90 b/code/constitutive.f90 index 63ae9a1a5..34328d300 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -436,7 +436,7 @@ end function constitutive_homogenizedC !-------------------------------------------------------------------------------------------------- !> @brief calls microstructure function of the different constitutive models !-------------------------------------------------------------------------------------------------- -subroutine constitutive_microstructure(orientations, Fe, Fp, ipc, ip, el) +subroutine constitutive_microstructure(orientations, Fe, Fp, ipc, ip, el, F0s,Fes,Fps) use prec, only: & pReal use material, only: & @@ -473,7 +473,10 @@ subroutine constitutive_microstructure(orientations, Fe, Fp, ipc, ip, el) ho, & !< homogenization tme !< thermal member position real(pReal), intent(in), dimension(:,:,:,:) :: & - orientations !< crystal orientations as quaternions + orientations, & + F0s, & + Fes, & + Fps !< crystal orientations as quaternions ho = material_homog(ip,el) tme = thermalMapping(ho)%p(ip,el) @@ -488,7 +491,7 @@ subroutine constitutive_microstructure(orientations, Fe, Fp, ipc, ip, el) case (PLASTICITY_NONLOCAL_ID) plasticityType call plastic_nonlocal_microstructure (Fe,Fp,ip,el) case (PLASTICITY_PHENOPLUS_ID) plasticityType - call plastic_phenoplus_microstructure(orientations,ipc,ip,el,Fe,Fp) + call plastic_phenoplus_microstructure(orientations,ipc,ip,el,F0s,Fes,Fps) end select plasticityType end subroutine constitutive_microstructure diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 71a5e4743..271fd2300 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -424,7 +424,7 @@ subroutine crystallite_init crystallite_partionedFp0 = crystallite_Fp0 crystallite_partionedFi0 = crystallite_Fi0 crystallite_partionedF0 = crystallite_F0 - crystallite_partionedF = crystallite_F0 + crystallite_partionedF = crystallite_F0 call crystallite_orientations() crystallite_orientation0 = crystallite_orientation ! store initial orientations for calculation of grain rotations @@ -437,7 +437,10 @@ subroutine crystallite_init call constitutive_microstructure(crystallite_orientation, & ! pass orientation to constitutive module crystallite_Fe(1:3,1:3,c,i,e), & crystallite_Fp(1:3,1:3,c,i,e), & - c,i,e) ! update dependent state variables to be consistent with basic states + c,i,e, + crystallite_F0, + crystallite_Fe, + crystallite_Fp) ! update dependent state variables to be consistent with basic states enddo enddo enddo @@ -1714,7 +1717,10 @@ subroutine crystallite_integrateStateRK4() call constitutive_microstructure(crystallite_orientation, & crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), & - g, i, e) ! update dependent state variables to be consistent with basic states + g, i, e, + crystallite_F0, + crystallite_Fe, + crystallite_Fp) ! update dependent state variables to be consistent with basic states enddo; enddo; enddo !$OMP ENDDO @@ -2040,7 +2046,10 @@ subroutine crystallite_integrateStateRKCK45() call constitutive_microstructure(crystallite_orientation, & crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), & - g, i, e) ! update dependent state variables to be consistent with basic states + g, i, e, + crystallite_F0, + crystallite_Fe, + crystallite_Fp) ! update dependent state variables to be consistent with basic states enddo; enddo; enddo !$OMP ENDDO @@ -2260,7 +2269,10 @@ subroutine crystallite_integrateStateRKCK45() call constitutive_microstructure(crystallite_orientation, & crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), & - g, i, e) ! update dependent state variables to be consistent with basic states + g, i, e, + crystallite_F0, + crystallite_Fe, + crystallite_Fp) ! update dependent state variables to be consistent with basic states enddo; enddo; enddo !$OMP ENDDO @@ -2495,7 +2507,10 @@ subroutine crystallite_integrateStateAdaptiveEuler() call constitutive_microstructure(crystallite_orientation, & crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), & - g, i, e) ! update dependent state variables to be consistent with basic states + g, i, e, + crystallite_F0, + crystallite_Fe, + crystallite_Fp) ! update dependent state variables to be consistent with basic states enddo; enddo; enddo !$OMP ENDDO !$OMP END PARALLEL @@ -2839,7 +2854,10 @@ eIter = FEsolving_execElem(1:2) call constitutive_microstructure(crystallite_orientation, & crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), & - g, i, e) ! update dependent state variables to be consistent with basic states + g, i, e, + crystallite_F0, + crystallite_Fe, + crystallite_Fp) ! update dependent state variables to be consistent with basic states enddo; enddo; enddo !$OMP ENDDO !$OMP END PARALLEL @@ -3084,7 +3102,10 @@ subroutine crystallite_integrateStateFPI() call constitutive_microstructure(crystallite_orientation, & crystallite_Fe(1:3,1:3,g,i,e), & crystallite_Fp(1:3,1:3,g,i,e), & - g, i, e) ! update dependent state variables to be consistent with basic states + g, i, e, + crystallite_F0, + crystallite_Fe, + crystallite_Fp) ! update dependent state variables to be consistent with basic states p = phaseAt(g,i,e) c = phasememberAt(g,i,e) plasticState(p)%previousDotState2(:,c) = plasticState(p)%previousDotState(:,c) diff --git a/code/plastic_phenoplus.f90 b/code/plastic_phenoplus.f90 index 9c125db04..b2a366dd3 100644 --- a/code/plastic_phenoplus.f90 +++ b/code/plastic_phenoplus.f90 @@ -739,9 +739,11 @@ end subroutine plastic_phenoplus_aTolState !-------------------------------------------------------------------------------------------------- !> @brief calculate push-up factors (kappa) for each voxel based on its neighbors !-------------------------------------------------------------------------------------------------- -subroutine plastic_phenoplus_microstructure(orientation,ipc,ip,el,Fe,Fp) +subroutine plastic_phenoplus_microstructure(orientation,ipc,ip,el,F0,Fe,Fp) use math, only: pi, & + math_identity2nd, & math_mul33x33, & + math_mul33xx33, & math_mul3x3, & math_transpose33, & math_qDot, & @@ -773,7 +775,8 @@ subroutine plastic_phenoplus_microstructure(orientation,ipc,ip,el,Fe,Fp) ipc, & !< component-ID of integration point ip, & !< integration point el - real(pReal), dimension(3,3), intent(in) :: & + real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + F0, & ! deformation gradient from last increment Fe, & ! elastic deformation gradient Fp ! elastic deformation gradient !< element real(pReal), dimension(4,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & @@ -813,10 +816,15 @@ subroutine plastic_phenoplus_microstructure(orientation,ipc,ip,el,Fe,Fp) tmp_acshear !temp holder for accumulative shear for m' real(pReal), dimension(3,3) :: & + F0_me, & !my deformation gradient from last converged increment Fe_me, & !my elastic deformation gradient Fp_me, & !my plastic deformation gradient + dF_me, & !my deformation gradient change (delta) + dE_me, & !my Green Lagrangian strain tensor (delta) Fe_ne, & !elastic deformation gradient of my current neighbor - Fp_ne !plastic deformation gradient of my current neighbor + Fp_ne, & !plastic deformation gradient of my current neighbor + dF_ne, & !deformation gradient of my current neighbor + dE_ne !delta Green Lagrangian strain tensor real(pReal), dimension(plastic_phenoplus_totalNslip(phase_plasticityInstance(material_phase(1,ip,el)))) :: & m_primes, & !m' between me_alpha(one) and neighbor beta(all) @@ -835,10 +843,6 @@ subroutine plastic_phenoplus_microstructure(orientation,ipc,ip,el,Fe,Fp) ne_mprimes !m' between each neighbor !***Get my properties - !@TODO - ! still need to know how to access the total strain for current material point - ! also, need to figure out an efficient way to calculate gamma_dot for the material - ! point and its neighbors Nneighbors = FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el)))) ph = phaseAt(ipc,ip,el) !get my phase of = phasememberAt(ipc,ip,el) !get my spatial location offset in memory @@ -853,13 +857,18 @@ subroutine plastic_phenoplus_microstructure(orientation,ipc,ip,el,Fe,Fp) mprime_cut = 0.7_pReal !set by Dr.Bieler dtaylor_cut = 1.0_pReal !set by Chen, quick test only - !***calculate my Taylor factor - !***gather my orientation and slip systems + !***gather my orientation, F and slip systems my_orientation = orientation(1:4, ipc, ip, el) - Fe_me = Fe(ipc,ip,el) - Fp_me = Fp(ipc,ip,el) + F0_me = F0(1:3, 1:3, ipc, ip, el) + Fe_me = Fe(1:3, 1:3, ipc, ip, el) + Fp_me = Fp(1:3, 1:3, ipc, ip, el) slipNormal(1:3, 1:ns) = lattice_sn(1:3, 1:ns, ph) slipDirect(1:3, 1:ns) = lattice_sd(1:3, 1:ns, ph) + !******calculate Taylor factor for me + !@note: we need teh + F_me = math_mul33x33(Fe_me,Fp_me) + E_me = 0.5*(math_mul33x33(math_transpose33(F_me), F_me) - math_identity2nd) !E = 0.5(F^tF-I) + vonStrain !***loop into the geometry to figure out who is my closest neighbor LOOPNEIGHBORS: DO n=1_pInt, Nneighbors