diff --git a/src/constitutive_plastic_disloUCLA.f90 b/src/constitutive_plastic_disloUCLA.f90 index fe69ae64d..05fd7bad8 100644 --- a/src/constitutive_plastic_disloUCLA.f90 +++ b/src/constitutive_plastic_disloUCLA.f90 @@ -12,15 +12,12 @@ submodule(constitutive) plastic_disloUCLA type :: tParameters real(pReal) :: & - D, & !< grain size - mu, & !< equivalent shear modulus - D_0, & !< prefactor for self-diffusion coefficient - Q_cl !< activation energy for dislocation climb + D = 1.0_pReal, & !< grain size + mu = 1.0_pReal, & !< equivalent shear modulus + D_0 = 1.0_pReal, & !< prefactor for self-diffusion coefficient + Q_cl = 1.0_pReal !< activation energy for dislocation climb real(pReal), allocatable, dimension(:) :: & - rho_mob_0, & !< initial dislocation density - rho_dip_0, & !< initial dipole density b_sl, & !< magnitude of burgers vector [m] - nonSchmidCoeff, & D_a, & i_sl, & !< Adj. parameter for distance between 2 forest dislocations atomicVolume, & @@ -36,7 +33,7 @@ submodule(constitutive) plastic_disloUCLA omega !< attempt frequency for kink pair nucleation real(pReal), allocatable, dimension(:,:) :: & h_sl_sl, & !< slip resistance from slip activity - forestProjectionEdge + forestProjection real(pReal), allocatable, dimension(:,:,:) :: & Schmid, & nonSchmid_pos, & @@ -87,6 +84,10 @@ module subroutine plastic_disloUCLA_init startIndex, endIndex integer, dimension(:), allocatable :: & N_sl + real, dimension(:), allocatable :: & + rho_mob_0, & !< initial dislocation density + rho_dip_0, & !< initial dipole density + a !< non-Schmid coefficients character(len=pStringLen) :: & extmsg = '' @@ -104,7 +105,6 @@ module subroutine plastic_disloUCLA_init allocate(dotState(Ninstance)) allocate(dependentState(Ninstance)) - do p = 1, size(phase_plasticity) if (phase_plasticity(p) /= PLASTICITY_DISLOUCLA_ID) cycle associate(prm => param(phase_plasticityInstance(p)), & @@ -127,25 +127,22 @@ module subroutine plastic_disloUCLA_init config%getFloat('c/a',defaultVal=0.0_pReal)) if(trim(config%getString('lattice_structure')) == 'bcc') then - prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',& - defaultVal = emptyRealArray) - prm%nonSchmid_pos = lattice_nonSchmidMatrix(N_sl,prm%nonSchmidCoeff,+1) - prm%nonSchmid_neg = lattice_nonSchmidMatrix(N_sl,prm%nonSchmidCoeff,-1) + a = config%getFloats('nonschmid_coefficients',defaultVal = emptyRealArray) + prm%nonSchmid_pos = lattice_nonSchmidMatrix(N_sl,a,+1) + prm%nonSchmid_neg = lattice_nonSchmidMatrix(N_sl,a,-1) else - allocate(prm%nonSchmidCoeff(0)) - prm%nonSchmid_pos = prm%Schmid - prm%nonSchmid_neg = prm%Schmid + prm%nonSchmid_pos = prm%Schmid + prm%nonSchmid_neg = prm%Schmid endif - prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl, & - config%getFloats('interaction_slipslip'), & - config%getString('lattice_structure')) - prm%forestProjectionEdge = lattice_forestProjection_edge(N_sl,config%getString('lattice_structure'),& - config%getFloat('c/a',defaultVal=0.0_pReal)) - prm%forestProjectionEdge = transpose(prm%forestProjectionEdge) + prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,config%getFloats('interaction_slipslip'), & + config%getString('lattice_structure')) + prm%forestProjection = lattice_forestProjection_edge(N_sl,config%getString('lattice_structure'),& + config%getFloat('c/a',defaultVal=0.0_pReal)) + prm%forestProjection = transpose(prm%forestProjection) - prm%rho_mob_0 = config%getFloats('rhoedge0', requiredSize=size(N_sl)) - prm%rho_dip_0 = config%getFloats('rhoedgedip0', requiredSize=size(N_sl)) + rho_mob_0 = config%getFloats('rhoedge0', requiredSize=size(N_sl)) + rho_dip_0 = config%getFloats('rhoedgedip0', requiredSize=size(N_sl)) prm%v0 = config%getFloats('v0', requiredSize=size(N_sl)) prm%b_sl = config%getFloats('slipburgers', requiredSize=size(N_sl)) prm%delta_F = config%getFloats('qedge', requiredSize=size(N_sl)) @@ -161,16 +158,16 @@ module subroutine plastic_disloUCLA_init prm%omega = config%getFloats('omega', requiredSize=size(N_sl)) prm%B = config%getFloats('friction_coeff', requiredSize=size(N_sl)) - prm%D = config%getFloat('grainsize') - prm%D_0 = config%getFloat('d0') - prm%Q_cl = config%getFloat('qsd') - prm%atomicVolume = config%getFloat('catomicvolume') * prm%b_sl**3.0_pReal - prm%D_a = config%getFloat('cedgedipmindistance') * prm%b_sl - prm%dipoleformation = config%getFloat('dipoleformationfactor') > 0.0_pReal !should be on by default, ToDo: change to /key/-type key + prm%D = config%getFloat('grainsize') + prm%D_0 = config%getFloat('d0') + prm%Q_cl = config%getFloat('qsd') + prm%atomicVolume = config%getFloat('catomicvolume') * prm%b_sl**3.0_pReal + prm%D_a = config%getFloat('cedgedipmindistance') * prm%b_sl + prm%dipoleformation = config%getFloat('dipoleformationfactor') > 0.0_pReal !should be on by default, ToDo: change to /key/-type key ! expand: family => system - prm%rho_mob_0 = math_expand(prm%rho_mob_0, N_sl) - prm%rho_dip_0 = math_expand(prm%rho_dip_0, N_sl) + rho_mob_0 = math_expand(rho_mob_0, N_sl) + rho_dip_0 = math_expand(rho_dip_0, N_sl) prm%q = math_expand(prm%q, N_sl) prm%p = math_expand(prm%p, N_sl) prm%delta_F = math_expand(prm%delta_F, N_sl) @@ -185,22 +182,25 @@ module subroutine plastic_disloUCLA_init prm%atomicVolume = math_expand(prm%atomicVolume, N_sl) prm%D_a = math_expand(prm%D_a, N_sl) - ! sanity checks - if ( prm%D_0 <= 0.0_pReal) extmsg = trim(extmsg)//' D_0' - if ( prm%Q_cl <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl' - if (any(prm%rho_mob_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoedge0' - if (any(prm%rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoedgedip0' - if (any(prm%v0 < 0.0_pReal)) extmsg = trim(extmsg)//' v0' - if (any(prm%b_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' slipb_sl' - if (any(prm%delta_F <= 0.0_pReal)) extmsg = trim(extmsg)//' qedge' - if (any(prm%tau_0 < 0.0_pReal)) extmsg = trim(extmsg)//' tau_0' - if (any(prm%D_a <= 0.0_pReal)) extmsg = trim(extmsg)//' cedgedipmindistance or slipb_sl' - if (any(prm%atomicVolume <= 0.0_pReal)) extmsg = trim(extmsg)//' catomicvolume or slipb_sl' + if ( prm%D_0 <= 0.0_pReal) extmsg = trim(extmsg)//' D_0' + if ( prm%Q_cl <= 0.0_pReal) extmsg = trim(extmsg)//' Q_cl' + if (any(rho_mob_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoedge0' + if (any(rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoedgedip0' + if (any(prm%v0 < 0.0_pReal)) extmsg = trim(extmsg)//' v0' + if (any(prm%b_sl <= 0.0_pReal)) extmsg = trim(extmsg)//' b_sl' + if (any(prm%delta_F <= 0.0_pReal)) extmsg = trim(extmsg)//' qedge' + if (any(prm%tau_0 < 0.0_pReal)) extmsg = trim(extmsg)//' tau_0' + if (any(prm%D_a <= 0.0_pReal)) extmsg = trim(extmsg)//' cedgedipmindistance or b_sl' + if (any(prm%atomicVolume <= 0.0_pReal)) extmsg = trim(extmsg)//' catomicvolume or b_sl' else slipActive - allocate(prm%rho_mob_0(0)) - allocate(prm%rho_dip_0(0)) + rho_mob_0= emptyRealArray; rho_dip_0 = emptyRealArray + allocate(prm%b_sl,prm%D_a,prm%i_sl,prm%atomicVolume,prm%tau_0, & + prm%delta_F,prm%v0,prm%p,prm%q,prm%B,prm%kink_height,prm%w,prm%omega, & + source = emptyRealArray) + allocate(prm%forestProjection(0,0)) + allocate(prm%h_sl_sl (0,0)) endif slipActive !-------------------------------------------------------------------------------------------------- @@ -215,23 +215,23 @@ module subroutine plastic_disloUCLA_init ! state aliases and initialization startIndex = 1 endIndex = prm%sum_N_sl - stt%rho_mob=>plasticState(p)%state(startIndex:endIndex,:) - stt%rho_mob= spread(prm%rho_mob_0,2,NipcMyPhase) - dot%rho_mob=>plasticState(p)%dotState(startIndex:endIndex,:) + stt%rho_mob => plasticState(p)%state(startIndex:endIndex,:) + stt%rho_mob = spread(rho_mob_0,2,NipcMyPhase) + dot%rho_mob => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_rho',defaultVal=1.0_pReal) if (any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho' startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl - stt%rho_dip=>plasticState(p)%state(startIndex:endIndex,:) - stt%rho_dip= spread(prm%rho_dip_0,2,NipcMyPhase) - dot%rho_dip=>plasticState(p)%dotState(startIndex:endIndex,:) + stt%rho_dip => plasticState(p)%state(startIndex:endIndex,:) + stt%rho_dip = spread(rho_dip_0,2,NipcMyPhase) + dot%rho_dip => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_rho',defaultVal=1.0_pReal) startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl - stt%gamma_sl=>plasticState(p)%state(startIndex:endIndex,:) - dot%gamma_sl=>plasticState(p)%dotState(startIndex:endIndex,:) + stt%gamma_sl => plasticState(p)%state(startIndex:endIndex,:) + dot%gamma_sl => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%atol(startIndex:endIndex) = 1.0e6_pReal ! ARRG ! global alias plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) @@ -256,7 +256,7 @@ end subroutine plastic_disloUCLA_init !> @brief Calculate plastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp, & - Mp,T,instance,of) + Mp,T,instance,of) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient real(pReal), dimension(3,3,3,3), intent(out) :: & @@ -369,8 +369,7 @@ module subroutine plastic_disloUCLA_dependentState(instance,of) associate(prm => param(instance), stt => state(instance),dst => dependentState(instance)) - dislocationSpacing = sqrt(matmul(prm%forestProjectionEdge, & - stt%rho_mob(:,of)+stt%rho_dip(:,of))) + dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,of)+stt%rho_dip(:,of))) dst%threshold_stress(:,of) = prm%mu*prm%b_sl & * sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,of)+stt%rho_dip(:,of))) diff --git a/src/constitutive_plastic_phenopowerlaw.f90 b/src/constitutive_plastic_phenopowerlaw.f90 index c96e1ea9a..a772d27ed 100644 --- a/src/constitutive_plastic_phenopowerlaw.f90 +++ b/src/constitutive_plastic_phenopowerlaw.f90 @@ -8,24 +8,21 @@ submodule(constitutive) plastic_phenopowerlaw type :: tParameters real(pReal) :: & - gdot0_slip, & !< reference shear strain rate for slip - gdot0_twin, & !< reference shear strain rate for twin - n_slip, & !< stress exponent for slip - n_twin, & !< stress exponent for twin - spr, & !< push-up factor for slip saturation due to twinning - c_1, & - c_2, & - c_3, & - c_4, & - h0_SlipSlip, & !< reference hardening slip - slip - h0_TwinSlip, & !< reference hardening twin - slip - h0_TwinTwin, & !< reference hardening twin - twin - a_slip + gdot0_slip = 1.0_pReal, & !< reference shear strain rate for slip + gdot0_twin = 1.0_pReal, & !< reference shear strain rate for twin + n_slip = 1.0_pReal, & !< stress exponent for slip + n_twin = 1.0_pReal, & !< stress exponent for twin + spr = 1.0_pReal, & !< push-up factor for slip saturation due to twinning + c_1 = 1.0_pReal, & + c_2 = 1.0_pReal, & + c_3 = 1.0_pReal, & + c_4 = 1.0_pReal, & + h0_SlipSlip = 1.0_pReal, & !< reference hardening slip - slip + h0_TwinSlip = 1.0_pReal, & !< reference hardening twin - slip + h0_TwinTwin = 1.0_pReal, & !< reference hardening twin - twin + a_slip = 1.0_pReal real(pReal), allocatable, dimension(:) :: & - xi_slip_0, & !< initial critical shear stress for slip - xi_twin_0, & !< initial critical shear stress for twin xi_slip_sat, & !< maximum critical shear stress for slip - nonSchmidCoeff, & H_int, & !< per family hardening activity (optional) gamma_twin_char !< characteristic shear for twins real(pReal), allocatable, dimension(:,:) :: & @@ -41,6 +38,8 @@ submodule(constitutive) plastic_phenopowerlaw integer :: & sum_N_sl, & !< total number of active slip system sum_N_tw !< total number of active twin systems + logical :: & + nonSchmidActive = .false. character(len=pStringLen), allocatable, dimension(:) :: & output end type tParameters @@ -77,6 +76,10 @@ module subroutine plastic_phenopowerlaw_init startIndex, endIndex integer, dimension(:), allocatable :: & N_sl, N_tw + real, dimension(:), allocatable :: & + xi_slip_0, & !< initial critical shear stress for slip + xi_twin_0, & !< initial critical shear stress for twin + a !< non-Schmid coefficients character(len=pStringLen) :: & extmsg = '' @@ -97,13 +100,6 @@ module subroutine plastic_phenopowerlaw_init stt => state(phase_plasticityInstance(p)), & config => config_phase(p)) -!-------------------------------------------------------------------------------------------------- -! optional parameters that need to be defined - prm%c_1 = config%getFloat('twin_c',defaultVal=0.0_pReal) - prm%c_2 = config%getFloat('twin_b',defaultVal=1.0_pReal) - prm%c_3 = config%getFloat('twin_e',defaultVal=0.0_pReal) - prm%c_4 = config%getFloat('twin_d',defaultVal=0.0_pReal) - !-------------------------------------------------------------------------------------------------- ! slip related parameters N_sl = config%getInts('nslip',defaultVal=emptyIntArray) @@ -113,12 +109,11 @@ module subroutine plastic_phenopowerlaw_init config%getFloat('c/a',defaultVal=0.0_pReal)) if(trim(config%getString('lattice_structure')) == 'bcc') then - prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',& - defaultVal = emptyRealArray) - prm%nonSchmid_pos = lattice_nonSchmidMatrix(N_sl,prm%nonSchmidCoeff,+1) - prm%nonSchmid_neg = lattice_nonSchmidMatrix(N_sl,prm%nonSchmidCoeff,-1) + a = config%getFloats('nonschmid_coefficients',defaultVal = emptyRealArray) + if(size(a) > 0) prm%nonSchmidActive = .true. + prm%nonSchmid_pos = lattice_nonSchmidMatrix(N_sl,a,+1) + prm%nonSchmid_neg = lattice_nonSchmidMatrix(N_sl,a,-1) else - allocate(prm%nonSchmidCoeff(0)) prm%nonSchmid_pos = prm%P_sl prm%nonSchmid_neg = prm%P_sl endif @@ -126,7 +121,7 @@ module subroutine plastic_phenopowerlaw_init config%getFloats('interaction_slipslip'), & config%getString('lattice_structure')) - prm%xi_slip_0 = config%getFloats('tau0_slip', requiredSize=size(N_sl)) + xi_slip_0 = config%getFloats('tau0_slip', requiredSize=size(N_sl)) prm%xi_slip_sat = config%getFloats('tausat_slip', requiredSize=size(N_sl)) prm%H_int = config%getFloats('h_int', requiredSize=size(N_sl), & defaultVal=[(0.0_pReal,i=1,size(N_sl))]) @@ -137,7 +132,7 @@ module subroutine plastic_phenopowerlaw_init prm%h0_SlipSlip = config%getFloat('h0_slipslip') ! expand: family => system - prm%xi_slip_0 = math_expand(prm%xi_slip_0, N_sl) + xi_slip_0 = math_expand(xi_slip_0, N_sl) prm%xi_slip_sat = math_expand(prm%xi_slip_sat,N_sl) prm%H_int = math_expand(prm%H_int, N_sl) @@ -145,11 +140,13 @@ module subroutine plastic_phenopowerlaw_init if ( prm%gdot0_slip <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0_slip' if ( prm%a_slip <= 0.0_pReal) extmsg = trim(extmsg)//' a_slip' if ( prm%n_slip <= 0.0_pReal) extmsg = trim(extmsg)//' n_slip' - if (any(prm%xi_slip_0 <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_slip_0' + if (any(xi_slip_0 <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_slip_0' if (any(prm%xi_slip_sat <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_slip_sat' + else slipActive + xi_slip_0 = emptyRealArray + allocate(prm%xi_slip_sat,prm%H_int,source=emptyRealArray) allocate(prm%interaction_SlipSlip(0,0)) - allocate(prm%xi_slip_0(0)) endif slipActive !-------------------------------------------------------------------------------------------------- @@ -165,23 +162,28 @@ module subroutine plastic_phenopowerlaw_init prm%gamma_twin_char = lattice_characteristicShear_twin(N_tw,config%getString('lattice_structure'),& config%getFloat('c/a')) - prm%xi_twin_0 = config%getFloats('tau0_twin',requiredSize=size(N_tw)) + xi_twin_0 = config%getFloats('tau0_twin',requiredSize=size(N_tw)) - prm%gdot0_twin = config%getFloat('gdot0_twin') - prm%n_twin = config%getFloat('n_twin') - prm%spr = config%getFloat('s_pr') - prm%h0_TwinTwin = config%getFloat('h0_twintwin') + prm%c_1 = config%getFloat('twin_c',defaultVal=0.0_pReal) + prm%c_2 = config%getFloat('twin_b',defaultVal=1.0_pReal) + prm%c_3 = config%getFloat('twin_e',defaultVal=0.0_pReal) + prm%c_4 = config%getFloat('twin_d',defaultVal=0.0_pReal) + prm%gdot0_twin = config%getFloat('gdot0_twin') + prm%n_twin = config%getFloat('n_twin') + prm%spr = config%getFloat('s_pr') + prm%h0_TwinTwin = config%getFloat('h0_twintwin') ! expand: family => system - prm%xi_twin_0 = math_expand(prm%xi_twin_0, N_tw) + xi_twin_0 = math_expand(xi_twin_0,N_tw) ! sanity checks if (prm%gdot0_twin <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0_twin' if (prm%n_twin <= 0.0_pReal) extmsg = trim(extmsg)//' n_twin' + else twinActive + xi_twin_0 = emptyRealArray + allocate(prm%gamma_twin_char,source=emptyRealArray) allocate(prm%interaction_TwinTwin(0,0)) - allocate(prm%xi_twin_0(0)) - allocate(prm%gamma_twin_char(0)) endif twinActive !-------------------------------------------------------------------------------------------------- @@ -218,7 +220,7 @@ module subroutine plastic_phenopowerlaw_init startIndex = 1 endIndex = prm%sum_N_sl stt%xi_slip => plasticState(p)%state (startIndex:endIndex,:) - stt%xi_slip = spread(prm%xi_slip_0, 2, NipcMyPhase) + stt%xi_slip = spread(xi_slip_0, 2, NipcMyPhase) dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_xi',defaultVal=1.0_pReal) if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' @@ -226,7 +228,7 @@ module subroutine plastic_phenopowerlaw_init startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_tw stt%xi_twin => plasticState(p)%state (startIndex:endIndex,:) - stt%xi_twin = spread(prm%xi_twin_0, 2, NipcMyPhase) + stt%xi_twin = spread(xi_twin_0, 2, NipcMyPhase) dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%atol(startIndex:endIndex) = config%getFloat('atol_xi',defaultVal=1.0_pReal) if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' @@ -431,20 +433,17 @@ pure subroutine kinetics_slip(Mp,instance,of, & tau_slip_pos, & tau_slip_neg integer :: i - logical :: nonSchmidActive associate(prm => param(instance), stt => state(instance)) - nonSchmidActive = size(prm%nonSchmidCoeff) > 0 - do i = 1, prm%sum_N_sl tau_slip_pos(i) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i)) tau_slip_neg(i) = merge(math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)), & - 0.0_pReal, nonSchmidActive) + 0.0_pReal, prm%nonSchmidActive) enddo where(dNeq0(tau_slip_pos)) - gdot_slip_pos = prm%gdot0_slip * merge(0.5_pReal,1.0_pReal, nonSchmidActive) & ! 1/2 if non-Schmid active + gdot_slip_pos = prm%gdot0_slip * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active * sign(abs(tau_slip_pos/stt%xi_slip(:,of))**prm%n_slip, tau_slip_pos) else where gdot_slip_pos = 0.0_pReal