diff --git a/code/config/material.config b/code/config/material.config index 73e7efcbf..8645cf7b8 100644 --- a/code/config/material.config +++ b/code/config/material.config @@ -267,7 +267,6 @@ rhoSglScrewNeg0 0.25e10 # Initial negative screw single rhoDipEdge0 1e8 # Initial edge dipole dislocation density in m/m**3 (per slip family) rhoDipScrew0 1e8 # Initial screw dipole dislocation density in m/m**3 (per slip family) rhoSglScatter 0 # standard deviation of scatter in initial single dislocation density -cutoffRadius 1e-3 # cutoff radius for dislocation stress in m minimumDipoleHeightEdge 2e-9 # minimum distance for stable edge dipoles in m (per slip family) minimumDipoleHeightScrew 2e-9 # minimum distance for stable screw dipoles in m (per slip family) lambda0 80 # prefactor for mean free path @@ -286,10 +285,12 @@ q 1 # exponent for thermal barrier attackFrequency 50e9 # attack frequency in Hz surfaceTransmissivity 1.0 # transmissivity of free surfaces for dislocation flux grainboundaryTransmissivity 0.0 # transmissivity of grain boundaries for dislocation flux (grain bundaries are identified as interfaces with different textures on both sides); if not set or set to negative number, the subroutine automatically determines the transmissivity at the grain boundary -shortRangeStressCorrection 0 # switch for use of short range correction stress interaction_SlipSlip 0 0 0.625 0.07 0.137 0.122 # Dislocation interaction coefficient +shortRangeStressCorrection 0 # switch for use of short range correction stress +cutoffRadius 1e-3 # cutoff radius for dislocation stress in m CFLfactor 2.0 # safety factor for CFL flux check (numerical parameter) -atol_rho 1e3 # dislocation density considered relevant in m/m**3 +significantRho 1e6 # minimum dislocation density considered relevant in m/m**3 +absoluteToleranceRho 1e4 # absolute tolerance for dislocation density in m/m**3 [BCC_Ferrite] diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index 08112a8f0..baa045e3f 100644 --- a/code/constitutive_nonlocal.f90 +++ b/code/constitutive_nonlocal.f90 @@ -115,6 +115,7 @@ constitutive_nonlocal_atomicVolume, & ! atomic vo constitutive_nonlocal_Dsd0, & ! prefactor for self-diffusion coefficient constitutive_nonlocal_Qsd, & ! activation enthalpy for diffusion constitutive_nonlocal_aTolRho, & ! absolute tolerance for dislocation density in state integration +constitutive_nonlocal_significantRho, & ! density considered significant constitutive_nonlocal_R, & ! cutoff radius for dislocation stress constitutive_nonlocal_doublekinkwidth, & ! width of a doubkle kink in multiples of the burgers vector length b constitutive_nonlocal_solidSolutionEnergy, & ! activation energy for solid solution in J @@ -317,6 +318,7 @@ allocate(constitutive_nonlocal_atomicVolume(maxNinstance)) allocate(constitutive_nonlocal_Dsd0(maxNinstance)) allocate(constitutive_nonlocal_Qsd(maxNinstance)) allocate(constitutive_nonlocal_aTolRho(maxNinstance)) +allocate(constitutive_nonlocal_significantRho(maxNinstance)) allocate(constitutive_nonlocal_Cslip_66(6,6,maxNinstance)) allocate(constitutive_nonlocal_Cslip_3333(3,3,3,3,maxNinstance)) allocate(constitutive_nonlocal_R(maxNinstance)) @@ -344,6 +346,7 @@ constitutive_nonlocal_atomicVolume = 0.0_pReal constitutive_nonlocal_Dsd0 = -1.0_pReal constitutive_nonlocal_Qsd = 0.0_pReal constitutive_nonlocal_aTolRho = 0.0_pReal +constitutive_nonlocal_significantRho = 0.0_pReal constitutive_nonlocal_nu = 0.0_pReal constitutive_nonlocal_Cslip_66 = 0.0_pReal constitutive_nonlocal_Cslip_3333 = 0.0_pReal @@ -469,15 +472,17 @@ do constitutive_nonlocal_Dsd0(i) = IO_floatValue(line,positions,2_pInt) case('selfdiffusionenergy','qsd') constitutive_nonlocal_Qsd(i) = IO_floatValue(line,positions,2_pInt) - case('atol_rho') + case('atol_rho','absolutetolerancerho','absolutetolerance_rho','absolutetolerancedensity','absolutetolerance_density') constitutive_nonlocal_aTolRho(i) = IO_floatValue(line,positions,2_pInt) + case('significantrho','significant_rho','significantdensity','significant_density') + constitutive_nonlocal_significantRho(i) = IO_floatValue(line,positions,2_pInt) case ('interaction_slipslip') forall (it = 1_pInt:lattice_maxNinteraction) & constitutive_nonlocal_interactionSlipSlip(it,i) = IO_floatValue(line,positions,1_pInt+it) - case('peierlsstressedge') + case('peierlsstressedge','peierlsstress_edge') forall (f = 1_pInt:lattice_maxNslipFamily) & constitutive_nonlocal_peierlsStressPerSlipFamily(f,1_pInt,i) = IO_floatValue(line,positions,1_pInt+f) - case('peierlsstressscrew') + case('peierlsstressscrew','peierlsstress_screw') forall (f = 1_pInt:lattice_maxNslipFamily) & constitutive_nonlocal_peierlsStressPerSlipFamily(f,2_pInt,i) = IO_floatValue(line,positions,1_pInt+f) case('doublekinkwidth') @@ -573,6 +578,8 @@ enddo //constitutive_nonlocal_label//')') if (constitutive_nonlocal_aTolRho(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='aTol_rho (' & //constitutive_nonlocal_label//')') + if (constitutive_nonlocal_significantRho(i) < 0.0_pReal) call IO_error(211_pInt,ext_msg='significantRho (' & + //constitutive_nonlocal_label//')') if (constitutive_nonlocal_doublekinkwidth(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='doublekinkwidth (' & //constitutive_nonlocal_label//')') if (constitutive_nonlocal_solidSolutionEnergy(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='solidSolutionEnergy (' & @@ -1129,6 +1136,14 @@ forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) & rhoDip(s,c) = max(state(g,ip,el)%p((7_pInt+c)*ns+s), 0.0_pReal) ! ensure positive dipole densities +!*** below significantRho, density quickly drops to zero + +rhoSgl = sign((rhoSgl**4.0_pReal + constitutive_nonlocal_significantRho(instance)**4.0_pReal)**0.25_pReal & + - constitutive_nonlocal_significantRho(instance), rhoSgl) +rhoDip = (rhoDip**4.0_pReal + constitutive_nonlocal_significantRho(instance)**4.0_pReal)**0.25_pReal & + - constitutive_nonlocal_significantRho(instance) + + !*** calculate the forest dislocation density !*** (= projection of screw and edge dislocations) @@ -1501,8 +1516,9 @@ 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(constitutive_nonlocal_totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),8) :: & + rhoSgl ! single dislocation densities (including used) real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),4) :: & - rhoSgl, & ! single dislocation densities (including used) v, & ! velocity dv_dtau ! velocity derivative with respect to the shear stress real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_plasticityInstance(material_phase(g,ip,el)))) :: & @@ -1526,9 +1542,17 @@ ns = constitutive_nonlocal_totalNslip(myInstance) forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) & rhoSgl(s,t) = max(state%p((t-1_pInt)*ns+s), 0.0_pReal) +forall (s = 1_pInt:ns, t = 5_pInt:8_pInt) & + rhoSgl(s,t) = state%p((t-1_pInt)*ns+s) tauBack = state%p(12_pInt*ns+1:13_pInt*ns) +!*** below significantRho, density quickly drops to zero + +rhoSgl = sign((rhoSgl**4.0_pReal + constitutive_nonlocal_significantRho(myInstance)**4.0_pReal)**0.25_pReal & + - constitutive_nonlocal_significantRho(myInstance), rhoSgl) + + !*** get effective resolved shear stress do s = 1_pInt,ns @@ -1557,14 +1581,14 @@ endif !*** Bauschinger effect -forall (s = 1_pInt:ns, t = 5_pInt:8_pInt, state%p((t-1)*ns+s) * v(s,t-4_pInt) < 0.0_pReal) & - rhoSgl(s,t-4_pInt) = rhoSgl(s,t-4_pInt) + abs(state%p((t-1_pInt)*ns+s)) +forall (s = 1_pInt:ns, t = 5_pInt:8_pInt, rhoSgl(s,t) * v(s,t-4_pInt) < 0.0_pReal) & + rhoSgl(s,t-4_pInt) = rhoSgl(s,t-4_pInt) + abs(rhoSgl(s,t)) !*** Calculation of gdot and its tangent -gdotTotal = sum(rhoSgl * v, 2) * constitutive_nonlocal_burgers(1:ns,myInstance) -dgdotTotal_dtau = sum(rhoSgl * dv_dtau, 2) * constitutive_nonlocal_burgers(1:ns,myInstance) +gdotTotal = sum(rhoSgl(1:ns,1:4) * v, 2) * constitutive_nonlocal_burgers(1:ns,myInstance) +dgdotTotal_dtau = sum(rhoSgl(1:ns,1:4) * dv_dtau, 2) * constitutive_nonlocal_burgers(1:ns,myInstance) !*** Calculation of Lp and its tangent @@ -1694,6 +1718,15 @@ forall (c = 1_pInt:2_pInt) & +!*** below significantRho, density quickly drops to zero + +rhoSgl = sign((rhoSgl**4.0_pReal + constitutive_nonlocal_significantRho(myInstance)**4.0_pReal)**0.25_pReal & + - constitutive_nonlocal_significantRho(myInstance), rhoSgl) +rhoDip = (rhoDip**4.0_pReal + constitutive_nonlocal_significantRho(myInstance)**4.0_pReal)**0.25_pReal & + - constitutive_nonlocal_significantRho(myInstance) + + + !**************************************************************************** !*** dislocation remobilization (bauschinger effect) @@ -1934,6 +1967,15 @@ forall (t = 1_pInt:4_pInt) & v(1_pInt:ns,t) = state(g,ip,el)%p((12_pInt+t)*ns+1_pInt:(13_pInt+t)*ns) +!*** below significantRho, density quickly drops to zero + +rhoSgl = sign((rhoSgl**4.0_pReal + constitutive_nonlocal_significantRho(myInstance)**4.0_pReal)**0.25_pReal & + - constitutive_nonlocal_significantRho(myInstance), rhoSgl) +rhoDip = (rhoDip**4.0_pReal + constitutive_nonlocal_significantRho(myInstance)**4.0_pReal)**0.25_pReal & + - constitutive_nonlocal_significantRho(myInstance) + + + !*** sanity check for timestep if (timestep <= 0.0_pReal) then ! if illegal timestep... @@ -2233,22 +2275,6 @@ if (numerics_integrationMode == 1_pInt) then endif -if ( any(rhoSgl(1:ns,1:4) + rhoDot(1:ns,1:4) * timestep < - constitutive_nonlocal_aTolRho(myInstance)) & - .or. any(rhoDip(1:ns,1:2) + rhoDot(1:ns,9:10) * timestep < - constitutive_nonlocal_aTolRho(myInstance))) then -#ifndef _OPENMP - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) 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 - constitutive_nonlocal_dotState = DAMASK_NaN - return -else - constitutive_nonlocal_dotState(1:10_pInt*ns) = reshape(rhoDot,(/10_pInt*ns/)) -endif - - - #ifndef _OPENMP if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt & .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)& @@ -2268,6 +2294,21 @@ endif endif #endif + +if ( any(rhoSgl(1:ns,1:4) + rhoDot(1:ns,1:4) * timestep < -constitutive_nonlocal_significantRho(myInstance)) & + .or. any(rhoDip(1:ns,1:2) + rhoDot(1:ns,9:10) * timestep < -constitutive_nonlocal_significantRho(myInstance))) then +#ifndef _OPENMP + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) 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 + constitutive_nonlocal_dotState = DAMASK_NaN + return +else + constitutive_nonlocal_dotState(1:10_pInt*ns) = reshape(rhoDot,(/10_pInt*ns/)) +endif + endfunction @@ -2590,6 +2631,12 @@ forall (t = 5_pInt:8_pInt) & rhoSgl(1:ns,t) = state(g,ip,el)%p((t-1_pInt)*ns+1_pInt:t*ns) +!*** below significantRho, density quickly drops to zero + +rhoSgl = sign((rhoSgl**4.0_pReal + constitutive_nonlocal_significantRho(instance)**4.0_pReal)**0.25_pReal & + - constitutive_nonlocal_significantRho(instance), rhoSgl) + + !*** calculate the dislocation stress of the neighboring excess dislocation densities !*** zero for material points of local plasticity @@ -2674,7 +2721,7 @@ ipLoop: do neighboring_ip = 1_pInt,FE_Nips(mesh_element(2,neighboring_el)) !* and add up the stress contributions from egde and screw excess on these slip systems (if significant) do s = 1_pInt,neighboring_ns - if (all(abs(neighboring_rhoExcess(:,:,s)) < constitutive_nonlocal_aTolRho(instance))) then + if (all(abs(neighboring_rhoExcess(:,:,s)) < constitutive_nonlocal_significantRho(instance))) then cycle ! not significant endif @@ -2696,7 +2743,7 @@ ipLoop: do neighboring_ip = 1_pInt,FE_Nips(mesh_element(2,neighboring_el)) zsquare = z * z do j = 1_pInt,2_pInt - if (abs(neighboring_rhoExcess(1,j,s)) < constitutive_nonlocal_aTolRho(instance)) then + if (abs(neighboring_rhoExcess(1,j,s)) < constitutive_nonlocal_significantRho(instance)) then cycle elseif (j > 1_pInt) then x = connection_neighboringSlip(1) + sign(0.5_pReal * segmentLength, & @@ -2736,7 +2783,7 @@ ipLoop: do neighboring_ip = 1_pInt,FE_Nips(mesh_element(2,neighboring_el)) x = connection_neighboringSlip(1) ! have to restore this value, because position might have been adapted for edge deads before do j = 1_pInt,2_pInt - if (abs(neighboring_rhoExcess(2,j,s)) < constitutive_nonlocal_aTolRho(instance)) then + if (abs(neighboring_rhoExcess(2,j,s)) < constitutive_nonlocal_significantRho(instance)) then cycle elseif (j > 1_pInt) then y = connection_neighboringSlip(2) + sign(0.5_pReal * segmentLength, & @@ -2799,7 +2846,7 @@ ipLoop: do neighboring_ip = 1_pInt,FE_Nips(mesh_element(2,neighboring_el)) + state(g,ip,el)%p((2_pInt*c+3_pInt)*ns+s) ! negative deads (here we use symmetry: if this has negative sign it is treated as positive density at positive position instead of negative density at negative position) do s = 1_pInt,ns - if (all(abs(rhoExcessDead(:,s)) < constitutive_nonlocal_aTolRho(instance))) then + if (all(abs(rhoExcessDead(:,s)) < constitutive_nonlocal_significantRho(instance))) then cycle ! not significant endif sigma = 0.0_pReal ! all components except for sigma13 are zero @@ -2922,8 +2969,12 @@ constitutive_nonlocal_postResults = 0.0_pReal !* short hand notations for state variables -forall (t = 1_pInt:8_pInt) rhoSgl(1:ns,t) = state(g,ip,el)%p((t-1_pInt)*ns+1_pInt:t*ns) -forall (c = 1_pInt:2_pInt) rhoDip(1:ns,c) = state(g,ip,el)%p((7_pInt+c)*ns+1_pInt:(8_pInt+c)*ns) +forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) & + rhoSgl(s,t) = max(state(g,ip,el)%p((t-1_pInt)*ns+s), 0.0_pReal) +forall (s = 1_pInt:ns, t = 5_pInt:8_pInt) & + rhoSgl(s,t) = state(g,ip,el)%p((t-1_pInt)*ns+s) +forall (c = 1_pInt:2_pInt) & + rhoDip(1:ns,c) = max(state(g,ip,el)%p((7_pInt+c)*ns+1_pInt:(8_pInt+c)*ns), 0.0_pReal) rhoForest = state(g,ip,el)%p(10_pInt*ns+1:11_pInt*ns) tauThreshold = state(g,ip,el)%p(11_pInt*ns+1:12_pInt*ns) tauBack = state(g,ip,el)%p(12_pInt*ns+1:13_pInt*ns) @@ -2932,6 +2983,15 @@ forall (c = 1_pInt:2_pInt) rhoDotDip(1:ns,c) = dotState%p((7_pInt+c)*ns+1_pInt:( forall (t = 1_pInt:4_pInt) v(1:ns,t) = state(g,ip,el)%p((12_pInt+t)*ns+1_pInt:(13_pInt+t)*ns) +!*** below significantRho, density quickly drops to zero + +rhoSgl = sign((rhoSgl**4.0_pReal + constitutive_nonlocal_significantRho(myInstance)**4.0_pReal)**0.25_pReal & + - constitutive_nonlocal_significantRho(myInstance), rhoSgl) +rhoDip = (rhoDip**4.0_pReal + constitutive_nonlocal_significantRho(myInstance)**4.0_pReal)**0.25_pReal & + - constitutive_nonlocal_significantRho(myInstance) + + + !* Calculate shear rate do t = 1_pInt,4_pInt @@ -3016,11 +3076,11 @@ do o = 1_pInt,phase_Noutput(material_phase(g,ip,el)) cs = cs + ns case ('rho_sgl_edge_pos_mobile') - constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,1) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = state(g,ip,el)%p(1:ns) cs = cs + ns case ('rho_sgl_edge_pos_immobile') - constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,5) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = state(g,ip,el)%p(4*ns+1:5*ns) cs = cs + ns case ('rho_sgl_edge_neg') @@ -3028,15 +3088,15 @@ do o = 1_pInt,phase_Noutput(material_phase(g,ip,el)) cs = cs + ns case ('rho_sgl_edge_neg_mobile') - constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,2) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = state(g,ip,el)%p(ns+1:2*ns) cs = cs + ns case ('rho_sgl_edge_neg_immobile') - constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,6) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = state(g,ip,el)%p(5*ns+1:6*ns) cs = cs + ns case ('rho_dip_edge') - constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDip(1:ns,1) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = state(g,ip,el)%p(8*ns+1:9*ns) cs = cs + ns case ('rho_screw') @@ -3060,11 +3120,11 @@ do o = 1_pInt,phase_Noutput(material_phase(g,ip,el)) cs = cs + ns case ('rho_sgl_screw_pos_mobile') - constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,3) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = state(g,ip,el)%p(2*ns+1:3*ns) cs = cs + ns case ('rho_sgl_screw_pos_immobile') - constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,7) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = state(g,ip,el)%p(6*ns+1:7*ns) cs = cs + ns case ('rho_sgl_screw_neg') @@ -3072,15 +3132,15 @@ do o = 1_pInt,phase_Noutput(material_phase(g,ip,el)) cs = cs + ns case ('rho_sgl_screw_neg_mobile') - constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,4) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = state(g,ip,el)%p(3*ns+1:4*ns) cs = cs + ns case ('rho_sgl_screw_neg_immobile') - constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,8) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = state(g,ip,el)%p(7*ns+1:8*ns) cs = cs + ns case ('rho_dip_screw') - constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDip(1:ns,2) + constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = state(g,ip,el)%p(9*ns+1:10*ns) cs = cs + ns case ('excess_rho')