diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index 153e12dd1..e80f2ceab 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -1347,16 +1347,13 @@ end subroutine plastic_nonlocal_dotState !--------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !--------------------------------------------------------------------------------------------------- -subroutine plastic_nonlocal_dotState2(Mp, F, Fp, Temperature,timestep, & - instance,of,ip,el) +subroutine plastic_nonlocal_dotState2(F, Fp,timestep, & + instance,of,ip,el) - real(pReal), dimension(3,3), intent(in) :: & - Mp !< MandelStress real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: & F, & !< elastic deformation gradient Fp !< plastic deformation gradient real(pReal), intent(in) :: & - Temperature, & !< temperature timestep !< substepped crystallite time increment integer, intent(in) :: & instance, & @@ -1385,12 +1382,7 @@ subroutine plastic_nonlocal_dotState2(Mp, F, Fp, Temperature,timestep, & real(pReal), dimension(param(instance)%sum_N_sl,10) :: & rho, & rho0, & !< dislocation density at beginning of time step - rhoDot, & !< density evolution - rhoDotMultiplication, & !< density evolution by multiplication - rhoDotFlux, & !< density evolution by flux - rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide) - rhoDotAthermalAnnihilation, & !< density evolution by athermal annihilation - rhoDotThermalAnnihilation !< density evolution by thermal annihilation + rhoDotFlux !< density evolution by flux real(pReal), dimension(param(instance)%sum_N_sl,8) :: & rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles) neighbor_rhoSgl0, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles) @@ -1400,13 +1392,6 @@ subroutine plastic_nonlocal_dotState2(Mp, F, Fp, Temperature,timestep, & v0, & neighbor_v0, & !< dislocation glide velocity of enighboring ip gdot !< shear rates - real(pReal), dimension(param(instance)%sum_N_sl) :: & - tau, & !< current resolved shear stress - vClimb !< climb velocity of edge dipoles - real(pReal), dimension(param(instance)%sum_N_sl,2) :: & - rhoDip, & !< current dipole dislocation densities (screw and edge dipoles) - dLower, & !< minimum stable dipole distance for edges and screws - dUpper !< current maximum stable dipole distance for edges and screws real(pReal), dimension(3,param(instance)%sum_N_sl,4) :: & m !< direction of dislocation motion real(pReal), dimension(3,3) :: & @@ -1423,8 +1408,7 @@ subroutine plastic_nonlocal_dotState2(Mp, F, Fp, Temperature,timestep, & real(pReal) :: & area, & !< area of the current interface transmissivity, & !< overall transmissivity of dislocation flux to neighboring material point - lineLength, & !< dislocation line length leaving the current interface - selfDiffusion !< self diffusion + lineLength !< dislocation line length leaving the current interface ph = material_phaseAt(1,el) if (timestep <= 0.0_pReal) then @@ -1438,63 +1422,16 @@ subroutine plastic_nonlocal_dotState2(Mp, F, Fp, Temperature,timestep, & stt => state(instance)) ns = prm%sum_N_sl - tau = 0.0_pReal gdot = 0.0_pReal rho = getRho(instance,of,ip,el) rhoSgl = rho(:,sgl) - rhoDip = rho(:,dip) rho0 = getRho0(instance,of,ip,el) my_rhoSgl0 = rho0(:,sgl) forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,instance),of) gdot = rhoSgl(:,1:4) * v * spread(prm%burgers,2,4) -#ifdef DEBUG - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0 & - .and. ((debug_e == el .and. debug_i == ip)& - .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0 )) then - write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> rho / 1/m^2', rhoSgl, rhoDip - write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> gdot / 1/s',gdot - endif -#endif - - !**************************************************************************** - !*** limits for stable dipole height - do s = 1,ns - tau(s) = math_tensordot(Mp, prm%Schmid(1:3,1:3,s)) + dst%tau_back(s,of) - if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal - enddo - - dLower = prm%minDipoleHeight - dUpper(:,1) = prm%mu * prm%burgers/(8.0_pReal * PI * (1.0_pReal - prm%nu) * abs(tau)) - dUpper(:,2) = prm%mu * prm%burgers/(4.0_pReal * PI * abs(tau)) - - where(dNeq0(sqrt(sum(abs(rho(:,edg)),2)))) & - dUpper(:,1) = min(1.0_pReal/sqrt(sum(abs(rho(:,edg)),2)),dUpper(:,1)) - where(dNeq0(sqrt(sum(abs(rho(:,scr)),2)))) & - dUpper(:,2) = min(1.0_pReal/sqrt(sum(abs(rho(:,scr)),2)),dUpper(:,2)) - - dUpper = max(dUpper,dLower) - - !**************************************************************************** - !*** dislocation multiplication - rhoDotMultiplication = 0.0_pReal - isBCC: if (lattice_structure(ph) == LATTICE_bcc_ID) then - forall (s = 1:ns, sum(abs(v(s,1:4))) > 0.0_pReal) - rhoDotMultiplication(s,1:2) = sum(abs(gdot(s,3:4))) / prm%burgers(s) & ! assuming double-cross-slip of screws to be decisive for multiplication - * sqrt(stt%rho_forest(s,of)) / prm%lambda0(s) ! & ! mean free path - ! * 2.0_pReal * sum(abs(v(s,3:4))) / sum(abs(v(s,1:4))) ! ratio of screw to overall velocity determines edge generation - rhoDotMultiplication(s,3:4) = sum(abs(gdot(s,3:4))) /prm%burgers(s) & ! assuming double-cross-slip of screws to be decisive for multiplication - * sqrt(stt%rho_forest(s,of)) / prm%lambda0(s) ! & ! mean free path - ! * 2.0_pReal * sum(abs(v(s,1:2))) / sum(abs(v(s,1:4))) ! ratio of edge to overall velocity determines screw generation - endforall - - else isBCC - rhoDotMultiplication(:,1:4) = spread( & - (sum(abs(gdot(:,1:2)),2) * prm%fEdgeMultiplication + sum(abs(gdot(:,3:4)),2)) & - * sqrt(stt%rho_forest(:,of)) / prm%lambda0 / prm%burgers, 2, 4) - endif isBCC forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,instance),of) @@ -1644,105 +1581,6 @@ subroutine plastic_nonlocal_dotState2(Mp, F, Fp, Temperature,timestep, & enddo neighbors endif - - - !**************************************************************************** - !*** calculate dipole formation and annihilation - - !*** formation by glide - do c = 1,2 - rhoDotSingle2DipoleGlide(:,2*c-1) = -2.0_pReal * dUpper(:,c) / prm%burgers & - * ( rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) & ! negative mobile --> positive mobile - + rhoSgl(:,2*c) * abs(gdot(:,2*c-1)) & ! positive mobile --> negative mobile - + abs(rhoSgl(:,2*c+4)) * abs(gdot(:,2*c-1))) ! positive mobile --> negative immobile - - rhoDotSingle2DipoleGlide(:,2*c) = -2.0_pReal * dUpper(:,c) / prm%burgers & - * ( rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) & ! negative mobile --> positive mobile - + rhoSgl(:,2*c) * abs(gdot(:,2*c-1)) & ! positive mobile --> negative mobile - + abs(rhoSgl(:,2*c+3)) * abs(gdot(:,2*c))) ! negative mobile --> positive immobile - - rhoDotSingle2DipoleGlide(:,2*c+3) = -2.0_pReal * dUpper(:,c) / prm%burgers & - * rhoSgl(:,2*c+3) * abs(gdot(:,2*c)) ! negative mobile --> positive immobile - - rhoDotSingle2DipoleGlide(:,2*c+4) = -2.0_pReal * dUpper(:,c) / prm%burgers & - * rhoSgl(:,2*c+4) * abs(gdot(:,2*c-1)) ! positive mobile --> negative immobile - - rhoDotSingle2DipoleGlide(:,c+8) = abs(rhoDotSingle2DipoleGlide(:,2*c+3)) & - + abs(rhoDotSingle2DipoleGlide(:,2*c+4)) & - - rhoDotSingle2DipoleGlide(:,2*c-1) & - - rhoDotSingle2DipoleGlide(:,2*c) - enddo - - - !*** athermal annihilation - rhoDotAthermalAnnihilation = 0.0_pReal - forall (c=1:2) & - rhoDotAthermalAnnihilation(:,c+8) = -2.0_pReal * dLower(:,c) / prm%burgers & - * ( 2.0_pReal * (rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) + rhoSgl(:,2*c) * abs(gdot(:,2*c-1))) & ! was single hitting single - + 2.0_pReal * (abs(rhoSgl(:,2*c+3)) * abs(gdot(:,2*c)) + abs(rhoSgl(:,2*c+4)) * abs(gdot(:,2*c-1))) & ! was single hitting immobile single or was immobile single hit by single - + rhoDip(:,c) * (abs(gdot(:,2*c-1)) + abs(gdot(:,2*c)))) ! single knocks dipole constituent - - ! annihilated screw dipoles leave edge jogs behind on the colinear system - if (lattice_structure(ph) == LATTICE_fcc_ID) & - forall (s = 1:ns, prm%colinearSystem(s) > 0) & - rhoDotAthermalAnnihilation(prm%colinearSystem(s),1:2) = - rhoDotAthermalAnnihilation(s,10) & - * 0.25_pReal * sqrt(stt%rho_forest(s,of)) * (dUpper(s,2) + dLower(s,2)) * prm%edgeJogFactor - - - !*** thermally activated annihilation of edge dipoles by climb - rhoDotThermalAnnihilation = 0.0_pReal - selfDiffusion = prm%Dsd0 * exp(-prm%selfDiffusionEnergy / (kB * Temperature)) - vClimb = prm%atomicVolume * selfDiffusion * prm%mu & - / ( kB * Temperature * PI * (1.0_pReal-prm%nu) * (dUpper(:,1) + dLower(:,1))) - forall (s = 1:ns, dUpper(s,1) > dLower(s,1)) & - rhoDotThermalAnnihilation(s,9) = max(- 4.0_pReal * rhoDip(s,1) * vClimb(s) / (dUpper(s,1) - dLower(s,1)), & - - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & - - rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have - - rhoDot = rhoDotFlux & - + rhoDotMultiplication & - + rhoDotSingle2DipoleGlide & - + rhoDotAthermalAnnihilation & - + rhoDotThermalAnnihilation - -#ifdef DEBUG - if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0 & - .and. ((debug_e == el .and. debug_i == ip)& - .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0 )) then - write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> dislocation multiplication', & - rhoDotMultiplication(:,1:4) * timestep - write(6,'(a,/,8(12x,12(e12.5,1x),/))') '<< CONST >> dislocation flux', & - rhoDotFlux(:,1:8) * timestep - write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> dipole formation by glide', & - rhoDotSingle2DipoleGlide * timestep - write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> athermal dipole annihilation', & - rhoDotAthermalAnnihilation * timestep - write(6,'(a,/,2(12x,12(e12.5,1x),/))') '<< CONST >> thermally activated dipole annihilation', & - rhoDotThermalAnnihilation(:,9:10) * timestep - write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> total density change', & - rhoDot * timestep - write(6,'(a,/,10(12x,12(f12.5,1x),/))') '<< CONST >> relative density change', & - rhoDot(:,1:8) * timestep / (abs(stt%rho(:,sgl))+1.0e-10), & - rhoDot(:,9:10) * timestep / (stt%rho(:,dip)+1.0e-10) - write(6,*) - endif -#endif - - - if ( any(rho(:,mob) + rhoDot(:,1:4) * timestep < -prm%atol_rho) & - .or. any(rho(:,dip) + rhoDot(:,9:10) * timestep < -prm%atol_rho)) then -#ifdef DEBUG - if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0) then - write(6,'(a,i5,a,i2)') '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip - write(6,'(a)') '<< CONST >> enforcing cutback !!!' - endif -#endif - plasticState(ph)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN) - else - dot%rho(:,of) = pack(rhoDot,.true.) - dot%gamma(:,of) = sum(gdot,2) - endif - end associate end subroutine plastic_nonlocal_dotState2