diff --git a/src/constitutive_plastic_dislotwin.f90 b/src/constitutive_plastic_dislotwin.f90 index 94e0177a0..152c630ae 100644 --- a/src/constitutive_plastic_dislotwin.f90 +++ b/src/constitutive_plastic_dislotwin.f90 @@ -8,35 +8,17 @@ !> @details to be done !-------------------------------------------------------------------------------------------------- submodule(constitutive) plastic_dislotwin - - real(pReal), parameter :: & + + real(pReal), parameter :: & kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin - - enum, bind(c) - enumerator :: & - undefined_ID, & - rho_mob_ID, & - rho_dip_ID, & - dot_gamma_sl_ID, & - gamma_sl_ID, & - Lambda_sl_ID, & - resolved_stress_slip_ID, & - tau_pass_ID, & - edge_dipole_distance_ID, & - f_tw_ID, & - Lambda_tw_ID, & - resolved_stress_twin_ID, & - tau_hat_tw_ID, & - f_tr_ID - end enum - + type :: tParameters real(pReal) :: & mu, & nu, & D0, & !< prefactor for self-diffusion coefficient Qsd, & !< activation energy for dislocation climb - omega, & !< frequency factor for dislocation climb + omega, & !< frequency factor for dislocation climb D, & !< grain size p_sb, & !< p-exponent in shear band velocity q_sb, & !< q-exponent in shear band velocity @@ -58,8 +40,8 @@ submodule(constitutive) plastic_dislotwin aTol_f_tr, & !< absolute tolerance for integration of trans volume fraction gamma_fcc_hex, & !< Free energy difference between austensite and martensite i_tr, & !< - h !< Stack height of hex nucleus - real(pReal), dimension(:), allocatable :: & + h !< Stack height of hex nucleus + real(pReal), allocatable, dimension(:) :: & rho_mob_0, & !< initial unipolar dislocation density per slip system rho_dip_0, & !< initial dipole dislocation density per slip system b_sl, & !< absolute length of burgers vector [m] for each slip system @@ -79,40 +61,39 @@ submodule(constitutive) plastic_dislotwin s, & !< s-exponent in trans nucleation rate gamma_char, & !< characteristic shear for twins B !< drag coefficient - real(pReal), dimension(:,:), allocatable :: & - h_sl_sl, & !< + real(pReal), allocatable, dimension(:,:) :: & + h_sl_sl, & !< h_sl_tw, & !< - h_tw_tw, & !< - h_sl_tr, & !< - h_tr_tr !< - integer, dimension(:,:), allocatable :: & - fcc_twinNucleationSlipPair ! ToDo: Better name? Is also use for trans - real(pReal), dimension(:,:), allocatable :: & + h_tw_tw, & !< + h_sl_tr, & !< + h_tr_tr, & !< n0_sl, & !< slip system normal forestProjection, & C66 - real(pReal), dimension(:,:,:), allocatable :: & + real(pReal), allocatable, dimension(:,:,:) :: & P_tr, & P_sl, & P_tw, & C66_tw, & C66_tr - integer :: & + integer :: & sum_N_sl, & !< total number of active slip system sum_N_tw, & !< total number of active twin system - sum_N_tr !< total number of active transformation system - integer, dimension(:), allocatable :: & + sum_N_tr !< total number of active transformation system + integer, allocatable, dimension(:) :: & N_sl, & !< number of active slip systems for each family N_tw, & !< number of active twin systems for each family N_tr !< number of active transformation systems for each family - integer(kind(undefined_ID)), dimension(:), allocatable :: & - outputID !< ID of each post result output + integer, allocatable, dimension(:,:) :: & + fcc_twinNucleationSlipPair ! ToDo: Better name? Is also use for trans + character(len=pStringLen), allocatable, dimension(:) :: & + output logical :: & ExtendedDislocations, & !< consider split into partials for climb calculation fccTwinTransNucleation, & !< twinning and transformation models are for fcc dipoleFormation !< flag indicating consideration of dipole formation end type !< container type for internal constitutive parameters - + type :: tDislotwinState real(pReal), dimension(:,:), pointer :: & rho_mob, & @@ -121,7 +102,7 @@ submodule(constitutive) plastic_dislotwin f_tw, & f_tr end type tDislotwinState - + type :: tDislotwinMicrostructure real(pReal), dimension(:,:), allocatable :: & Lambda_sl, & !< mean free path between 2 obstacles seen by a moving dislocation @@ -135,7 +116,7 @@ submodule(constitutive) plastic_dislotwin tau_r_tw, & !< stress to bring partials close together (twin) tau_r_tr !< stress to bring partials close together (trans) end type tDislotwinMicrostructure - + !-------------------------------------------------------------------------------------------------- ! containers for parameters and state type(tParameters), allocatable, dimension(:) :: param @@ -148,7 +129,7 @@ contains !-------------------------------------------------------------------------------------------------- -!> @brief module initialization +!> @brief Perform module initialization. !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- module subroutine plastic_dislotwin_init @@ -156,39 +137,34 @@ module subroutine plastic_dislotwin_init integer :: & Ninstance, & p, i, & - NipcMyPhase, outputSize, & + NipcMyPhase, & sizeState, sizeDotState, & startIndex, endIndex - - integer(kind(undefined_ID)) :: & - outputID - + character(len=pStringLen) :: & extmsg = '' - character(len=pStringLen), dimension(:), allocatable :: & - outputs - - write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_DISLOTWIN_label//' init -+>>>' - - write(6,'(/,a)') ' Ma and Roters, Acta Materialia 52(12):3603–3612, 2004' - write(6,'(a)') ' https://doi.org/10.1016/j.actamat.2004.04.012' - - write(6,'(/,a)') ' Roters et al., Computational Materials Science 39:91–95, 2007' - write(6,'(a)') ' https://doi.org/10.1016/j.commatsci.2006.04.014' - - write(6,'(/,a)') ' Wong et al., Acta Materialia 118:140–151, 2016' - write(6,'(a,/)') ' https://doi.org/10.1016/j.actamat.2016.07.032' - + + write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_DISLOTWIN_label//' init -+>>>'; flush(6) + + write(6,'(/,a)') ' Ma and Roters, Acta Materialia 52(12):3603–3612, 2004' + write(6,'(a)') ' https://doi.org/10.1016/j.actamat.2004.04.012' + + write(6,'(/,a)') ' Roters et al., Computational Materials Science 39:91–95, 2007' + write(6,'(a)') ' https://doi.org/10.1016/j.commatsci.2006.04.014' + + write(6,'(/,a)') ' Wong et al., Acta Materialia 118:140–151, 2016' + write(6,'(a,/)') ' https://doi.org/10.1016/j.actamat.2016.07.032' + Ninstance = count(phase_plasticity == PLASTICITY_DISLOTWIN_ID) - + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - + allocate(param(Ninstance)) allocate(state(Ninstance)) allocate(dotState(Ninstance)) allocate(dependentState(Ninstance)) - + do p = 1, size(phase_plasticity) if (phase_plasticity(p) /= PLASTICITY_DISLOTWIN_ID) cycle associate(prm => param(phase_plasticityInstance(p)), & @@ -196,17 +172,16 @@ module subroutine plastic_dislotwin_init stt => state(phase_plasticityInstance(p)), & dst => dependentState(phase_plasticityInstance(p)), & config => config_phase(p)) - + prm%aTol_rho = config%getFloat('atol_rho', defaultVal=0.0_pReal) prm%aTol_f_tw = config%getFloat('atol_twinfrac', defaultVal=0.0_pReal) prm%aTol_f_tr = config%getFloat('atol_transfrac', defaultVal=0.0_pReal) - + ! This data is read in already in lattice prm%mu = lattice_mu(p) prm%nu = lattice_nu(p) prm%C66 = lattice_C66(1:6,1:6,p) - - + !-------------------------------------------------------------------------------------------------- ! slip related parameters prm%N_sl = config%getInts('nslip',defaultVal=emptyIntArray) @@ -220,14 +195,14 @@ module subroutine plastic_dislotwin_init prm%forestProjection = lattice_forestProjection_edge(prm%N_sl,config%getString('lattice_structure'),& config%getFloat('c/a',defaultVal=0.0_pReal)) prm%forestProjection = transpose(prm%forestProjection) - + prm%n0_sl = lattice_slip_normal(prm%N_sl,config%getString('lattice_structure'),& - config%getFloat('c/a',defaultVal=0.0_pReal)) + config%getFloat('c/a',defaultVal=0.0_pReal)) prm%fccTwinTransNucleation = merge(.true., .false., lattice_structure(p) == LATTICE_FCC_ID) & .and. (prm%N_sl(1) == 12) if(prm%fccTwinTransNucleation) & prm%fcc_twinNucleationSlipPair = lattice_fcc_twinNucleationSlipPair - + prm%rho_mob_0 = config%getFloats('rhoedge0', requiredSize=size(prm%N_sl)) prm%rho_dip_0 = config%getFloats('rhoedgedip0',requiredSize=size(prm%N_sl)) prm%v0 = config%getFloats('v0', requiredSize=size(prm%N_sl)) @@ -238,7 +213,7 @@ module subroutine plastic_dislotwin_init prm%q = config%getFloats('q_slip', requiredSize=size(prm%N_sl)) prm%B = config%getFloats('b', requiredSize=size(prm%N_sl), & defaultVal=[(0.0_pReal, i=1,size(prm%N_sl))]) - + prm%tau_0 = config%getFloat('solidsolutionstrength') prm%CEdgeDipMinDistance = config%getFloat('cedgedipmindistance') prm%D0 = config%getFloat('d0') @@ -249,15 +224,15 @@ module subroutine plastic_dislotwin_init prm%SFE_0K = config%getFloat('sfe_0k') prm%dSFE_dT = config%getFloat('dsfe_dt') endif - + ! multiplication factor according to crystal structure (nearest neighbors bcc vs fcc/hex) !@details: Refer: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 prm%omega = config%getFloat('omega', defaultVal = 1000.0_pReal) & * merge(12.0_pReal, & 8.0_pReal, & lattice_structure(p) == LATTICE_FCC_ID .or. lattice_structure(p) == LATTICE_HEX_ID) - - + + ! expand: family => system prm%rho_mob_0 = math_expand(prm%rho_mob_0, prm%N_sl) prm%rho_dip_0 = math_expand(prm%rho_dip_0, prm%N_sl) @@ -268,8 +243,8 @@ module subroutine plastic_dislotwin_init prm%p = math_expand(prm%p, prm%N_sl) prm%q = math_expand(prm%q, prm%N_sl) prm%B = math_expand(prm%B, prm%N_sl) - prm%atomicVolume = math_expand(prm%atomicVolume,prm%N_sl) - + prm%atomicVolume = math_expand(prm%atomicVolume,prm%N_sl) + ! sanity checks if ( prm%D0 <= 0.0_pReal) extmsg = trim(extmsg)//' D0' if ( prm%Qsd <= 0.0_pReal) extmsg = trim(extmsg)//' Qsd' @@ -282,11 +257,11 @@ module subroutine plastic_dislotwin_init if (any(prm%B < 0.0_pReal)) extmsg = trim(extmsg)//' B' if (any(prm%p<=0.0_pReal .or. prm%p>1.0_pReal)) extmsg = trim(extmsg)//' p' if (any(prm%q< 1.0_pReal .or. prm%q>2.0_pReal)) extmsg = trim(extmsg)//' q' - + else slipActive allocate(prm%b_sl(0)) endif slipActive - + !-------------------------------------------------------------------------------------------------- ! twin related parameters prm%N_tw = config%getInts('ntwin', defaultVal=emptyIntArray) @@ -297,31 +272,31 @@ module subroutine plastic_dislotwin_init prm%h_tw_tw = lattice_interaction_TwinByTwin(prm%N_tw,& config%getFloats('interaction_twintwin'), & config%getString('lattice_structure')) - + prm%b_tw = config%getFloats('twinburgers', requiredSize=size(prm%N_tw)) prm%t_tw = config%getFloats('twinsize', requiredSize=size(prm%N_tw)) prm%r = config%getFloats('r_twin', requiredSize=size(prm%N_tw)) - + prm%xc_twin = config%getFloat('xc_twin') prm%L_tw = config%getFloat('l0_twin') prm%i_tw = config%getFloat('cmfptwin') - + prm%gamma_char= lattice_characteristicShear_Twin(prm%N_tw,config%getString('lattice_structure'),& config%getFloat('c/a',defaultVal=0.0_pReal)) - + prm%C66_tw = lattice_C66_twin(prm%N_tw,prm%C66,config%getString('lattice_structure'),& config%getFloat('c/a',defaultVal=0.0_pReal)) - + if (.not. prm%fccTwinTransNucleation) then - prm%dot_N_0_tw = config%getFloats('ndot0_twin') + prm%dot_N_0_tw = config%getFloats('ndot0_twin') prm%dot_N_0_tw = math_expand(prm%dot_N_0_tw,prm%N_tw) endif - + ! expand: family => system prm%b_tw = math_expand(prm%b_tw,prm%N_tw) prm%t_tw = math_expand(prm%t_tw,prm%N_tw) prm%r = math_expand(prm%r,prm%N_tw) - + else allocate(prm%gamma_char(0)) allocate(prm%t_tw (0)) @@ -329,7 +304,7 @@ module subroutine plastic_dislotwin_init allocate(prm%r (0)) allocate(prm%h_tw_tw (0,0)) endif - + !-------------------------------------------------------------------------------------------------- ! transformation related parameters prm%N_tr = config%getInts('ntrans', defaultVal=emptyIntArray) @@ -337,29 +312,29 @@ module subroutine plastic_dislotwin_init if (prm%sum_N_tr > 0) then prm%b_tr = config%getFloats('transburgers') prm%b_tr = math_expand(prm%b_tr,prm%N_tr) - + prm%h = config%getFloat('transstackheight', defaultVal=0.0_pReal) ! ToDo: How to handle that??? prm%i_tr = config%getFloat('cmfptrans', defaultVal=0.0_pReal) ! ToDo: How to handle that??? prm%gamma_fcc_hex = config%getFloat('deltag') prm%xc_trans = config%getFloat('xc_trans', defaultVal=0.0_pReal) ! ToDo: How to handle that??? prm%L_tr = config%getFloat('l0_trans') - + prm%h_tr_tr = lattice_interaction_TransByTrans(prm%N_tr,& config%getFloats('interaction_transtrans'), & config%getString('lattice_structure')) - + prm%C66_tr = lattice_C66_trans(prm%N_tr,prm%C66, & config%getString('trans_lattice_structure'), & 0.0_pReal, & config%getFloat('a_bcc', defaultVal=0.0_pReal), & config%getFloat('a_fcc', defaultVal=0.0_pReal)) - + prm%P_tr = lattice_SchmidMatrix_trans(prm%N_tr, & config%getString('trans_lattice_structure'), & 0.0_pReal, & config%getFloat('a_bcc', defaultVal=0.0_pReal), & config%getFloat('a_fcc', defaultVal=0.0_pReal)) - + if (lattice_structure(p) /= LATTICE_fcc_ID) then prm%dot_N_0_tr = config%getFloats('ndot0_trans') prm%dot_N_0_tr = math_expand(prm%dot_N_0_tr,prm%N_tr) @@ -374,59 +349,57 @@ module subroutine plastic_dislotwin_init allocate(prm%s (0)) allocate(prm%h_tr_tr(0,0)) endif - + if (sum(prm%N_tw) > 0 .or. prm%sum_N_tr > 0) then prm%SFE_0K = config%getFloat('sfe_0k') prm%dSFE_dT = config%getFloat('dsfe_dt') prm%V_cs = config%getFloat('vcrossslip') endif - + if (prm%sum_N_sl > 0 .and. prm%sum_N_tw > 0) then prm%h_sl_tw = lattice_interaction_SlipByTwin(prm%N_sl,prm%N_tw,& config%getFloats('interaction_sliptwin'), & config%getString('lattice_structure')) if (prm%fccTwinTransNucleation .and. prm%sum_N_tw > 12) write(6,*) 'mist' ! ToDo: implement better test. The model will fail also if N_tw is [6,6] - endif - - if (prm%sum_N_sl > 0 .and. prm%sum_N_tr > 0) then + endif + + if (prm%sum_N_sl > 0 .and. prm%sum_N_tr > 0) then prm%h_sl_tr = lattice_interaction_SlipByTrans(prm%N_sl,prm%N_tr,& config%getFloats('interaction_sliptrans'), & config%getString('lattice_structure')) if (prm%fccTwinTransNucleation .and. prm%sum_N_tr > 12) write(6,*) 'mist' ! ToDo: implement better test. The model will fail also if N_tr is [6,6] - endif - + endif + !-------------------------------------------------------------------------------------------------- ! shearband related parameters prm%sbVelocity = config%getFloat('shearbandvelocity',defaultVal=0.0_pReal) - if (prm%sbVelocity > 0.0_pReal) then + if (prm%sbVelocity > 0.0_pReal) then prm%sbResistance = config%getFloat('shearbandresistance') prm%sbQedge = config%getFloat('qedgepersbsystem') prm%p_sb = config%getFloat('p_shearband') prm%q_sb = config%getFloat('q_shearband') - + ! sanity checks if (prm%sbResistance < 0.0_pReal) extmsg = trim(extmsg)//' shearbandresistance' if (prm%sbQedge < 0.0_pReal) extmsg = trim(extmsg)//' qedgepersbsystem' if (prm%p_sb <= 0.0_pReal) extmsg = trim(extmsg)//' p_shearband' if (prm%q_sb <= 0.0_pReal) extmsg = trim(extmsg)//' q_shearband' endif - - - + prm%D = config%getFloat('grainsize') - + if (config%keyExists('dipoleformationfactor')) call IO_error(1,ext_msg='use /nodipoleformation/') prm%dipoleformation = .not. config%keyExists('/nodipoleformation/') - - + + !if (Ndot0PerTwinFamily(f,p) < 0.0_pReal) & ! call IO_error(211,el=p,ext_msg='dot_N_0_tw ('//PLASTICITY_DISLOTWIN_label//')') - + if (any(prm%atomicVolume <= 0.0_pReal)) & call IO_error(211,el=p,ext_msg='cAtomicVolume ('//PLASTICITY_DISLOTWIN_label//')') if (prm%sum_N_tw > 0) then if (prm%aTol_rho <= 0.0_pReal) & - call IO_error(211,el=p,ext_msg='aTol_rho ('//PLASTICITY_DISLOTWIN_label//')') + call IO_error(211,el=p,ext_msg='aTol_rho ('//PLASTICITY_DISLOTWIN_label//')') if (prm%aTol_f_tw <= 0.0_pReal) & call IO_error(211,el=p,ext_msg='aTol_f_tw ('//PLASTICITY_DISLOTWIN_label//')') endif @@ -434,50 +407,9 @@ module subroutine plastic_dislotwin_init if (prm%aTol_f_tr <= 0.0_pReal) & call IO_error(211,el=p,ext_msg='aTol_f_tr ('//PLASTICITY_DISLOTWIN_label//')') endif - - outputs = config%getStrings('(output)', defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - do i= 1, size(outputs) - outputID = undefined_ID - select case(outputs(i)) - case ('rho_mob') - outputID = merge(rho_mob_ID,undefined_ID,prm%sum_N_sl > 0) - outputSize = prm%sum_N_sl - case ('rho_dip') - outputID = merge(rho_dip_ID,undefined_ID,prm%sum_N_sl > 0) - outputSize = prm%sum_N_sl - case ('gamma_sl') - outputID = merge(gamma_sl_ID,undefined_ID,prm%sum_N_sl > 0) - outputSize = prm%sum_N_sl - case ('lambda_sl') - outputID = merge(Lambda_sl_ID,undefined_ID,prm%sum_N_sl > 0) - outputSize = prm%sum_N_sl - case ('tau_pass') - outputID= merge(tau_pass_ID,undefined_ID,prm%sum_N_sl > 0) - outputSize = prm%sum_N_sl - - case ('f_tw') - outputID = merge(f_tw_ID,undefined_ID,prm%sum_N_tw >0) - outputSize = prm%sum_N_tw - case ('lambda_tw') - outputID = merge(Lambda_tw_ID,undefined_ID,prm%sum_N_tw >0) - outputSize = prm%sum_N_tw - case ('tau_hat_tw') - outputID = merge(tau_hat_tw_ID,undefined_ID,prm%sum_N_tw >0) - outputSize = prm%sum_N_tw - - case ('f_tr') - outputID = f_tr_ID - outputSize = prm%sum_N_tr - - end select - - if (outputID /= undefined_ID) then - prm%outputID = [prm%outputID, outputID] - endif - - enddo - + + prm%output = config%getStrings('(output)', defaultVal=emptyStringArray) + !-------------------------------------------------------------------------------------------------- ! allocate state arrays NipcMyPhase = count(material_phaseAt == p) * discretization_nIP @@ -485,9 +417,9 @@ module subroutine plastic_dislotwin_init + size(['f_tw']) * prm%sum_N_tw & + size(['f_tr']) * prm%sum_N_tr sizeState = sizeDotState - + call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) - + !-------------------------------------------------------------------------------------------------- ! locally defined state aliases and initialization of state0 and aTolState startIndex = 1 @@ -496,14 +428,14 @@ module subroutine plastic_dislotwin_init stt%rho_mob= spread(prm%rho_mob_0,2,NipcMyPhase) dot%rho_mob=>plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%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,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTol_rho - + startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl stt%gamma_sl=>plasticState(p)%state(startIndex:endIndex,:) @@ -511,66 +443,66 @@ module subroutine plastic_dislotwin_init plasticState(p)%aTolState(startIndex:endIndex) = 1.0e6_pReal !ToDo: better make optional parameter ! global alias plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) - + startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_tw stt%f_tw=>plasticState(p)%state(startIndex:endIndex,:) dot%f_tw=>plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTol_f_tw - + startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_tr stt%f_tr=>plasticState(p)%state(startIndex:endIndex,:) dot%f_tr=>plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTol_f_tr - + allocate(dst%Lambda_sl (prm%sum_N_sl,NipcMyPhase),source=0.0_pReal) allocate(dst%tau_pass (prm%sum_N_sl,NipcMyPhase),source=0.0_pReal) - + allocate(dst%Lambda_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal) allocate(dst%tau_hat_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal) allocate(dst%tau_r_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal) allocate(dst%V_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal) - + allocate(dst%Lambda_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal) allocate(dst%tau_hat_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal) allocate(dst%tau_r_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal) allocate(dst%V_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal) - - + + plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally - + end associate - + enddo end subroutine plastic_dislotwin_init !-------------------------------------------------------------------------------------------------- -!> @brief returns the homogenized elasticity matrix +!> @brief Return the homogenized elasticity matrix. !-------------------------------------------------------------------------------------------------- module function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC) - + real(pReal), dimension(6,6) :: & homogenizedC integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element - + integer :: i, & of real(pReal) :: f_unrotated - + of = material_phasememberAt(ipc,ip,el) associate(prm => param(phase_plasticityInstance(material_phaseAt(ipc,el))),& stt => state(phase_plasticityInstance(material_phaseAT(ipc,el)))) - + f_unrotated = 1.0_pReal & - sum(stt%f_tw(1:prm%sum_N_tw,of)) & - sum(stt%f_tr(1:prm%sum_N_tr,of)) - + homogenizedC = f_unrotated * prm%C66 do i=1,prm%sum_N_tw homogenizedC = homogenizedC & @@ -580,23 +512,23 @@ module function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC) homogenizedC = homogenizedC & + stt%f_tr(i,of)*prm%C66_tr(1:6,1:6,i) enddo - + end associate - + end function plastic_dislotwin_homogenizedC !-------------------------------------------------------------------------------------------------- -!> @brief calculates plastic velocity gradient and its tangent +!> @brief Calculate plastic velocity gradient and its tangent. !-------------------------------------------------------------------------------------------------- module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) - + real(pReal), dimension(3,3), intent(out) :: Lp real(pReal), dimension(3,3,3,3), intent(out) :: dLp_dMp real(pReal), dimension(3,3), intent(in) :: Mp integer, intent(in) :: instance,of real(pReal), intent(in) :: T - + integer :: i,k,l,m,n real(pReal) :: & f_unrotated,StressRatio_p,& @@ -632,16 +564,16 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) 0, 1,-1, & 0, 1, 1 & ],pReal),[ 3,6]) - + associate(prm => param(instance), stt => state(instance)) - + f_unrotated = 1.0_pReal & - sum(stt%f_tw(1:prm%sum_N_tw,of)) & - sum(stt%f_tr(1:prm%sum_N_tr,of)) - + Lp = 0.0_pReal - dLp_dMp = 0.0_pReal - + dLp_dMp = 0.0_pReal + call kinetics_slip(Mp,T,instance,of,dot_gamma_sl,ddot_gamma_dtau_slip) slipContribution: do i = 1, prm%sum_N_sl Lp = Lp + dot_gamma_sl(i)*prm%P_sl(1:3,1:3,i) @@ -649,37 +581,37 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + ddot_gamma_dtau_slip(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i) enddo slipContribution - + !ToDo: Why do this before shear banding? Lp = Lp * f_unrotated dLp_dMp = dLp_dMp * f_unrotated - + shearBandingContribution: if(dNeq0(prm%sbVelocity)) then - + BoltzmannRatio = prm%sbQedge/(kB*T) call math_eigenValuesVectorsSym(Mp,eigValues,eigVectors,error) - + do i = 1,6 P_sb = 0.5_pReal * math_outer(matmul(eigVectors,sb_sComposition(1:3,i)),& matmul(eigVectors,sb_mComposition(1:3,i))) tau = math_mul33xx33(Mp,P_sb) - + significantShearBandStress: if (abs(tau) > tol_math_check) then StressRatio_p = (abs(tau)/prm%sbResistance)**prm%p_sb dot_gamma_sb = sign(prm%sbVelocity*exp(-BoltzmannRatio*(1-StressRatio_p)**prm%q_sb), tau) ddot_gamma_dtau = abs(dot_gamma_sb)*BoltzmannRatio* prm%p_sb*prm%q_sb/ prm%sbResistance & * (abs(tau)/prm%sbResistance)**(prm%p_sb-1.0_pReal) & * (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal) - + Lp = Lp + dot_gamma_sb * P_sb forall (k=1:3,l=1:3,m=1:3,n=1:3) & dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + ddot_gamma_dtau * P_sb(k,l) * P_sb(m,n) endif significantShearBandStress enddo - + endif shearBandingContribution - + call kinetics_twin(Mp,T,dot_gamma_sl,instance,of,dot_gamma_twin,ddot_gamma_dtau_twin) twinContibution: do i = 1, prm%sum_N_tw Lp = Lp + dot_gamma_twin(i)*prm%P_tw(1:3,1:3,i) * f_unrotated @@ -687,7 +619,7 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + ddot_gamma_dtau_twin(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i) * f_unrotated enddo twinContibution - + call kinetics_trans(Mp,T,dot_gamma_sl,instance,of,dot_gamma_tr,ddot_gamma_dtau_trans) transContibution: do i = 1, prm%sum_N_tr Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i) * f_unrotated @@ -695,15 +627,15 @@ module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) & + ddot_gamma_dtau_trans(i)* prm%P_tr(k,l,i)*prm%P_tr(m,n,i) * f_unrotated enddo transContibution - - + + end associate - + end subroutine plastic_dislotwin_LpAndItsTangent !-------------------------------------------------------------------------------------------------- -!> @brief calculates the rate of change of microstructure +!> @brief Calculate the rate of change of microstructure. !-------------------------------------------------------------------------------------------------- module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) @@ -720,10 +652,10 @@ module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) f_unrotated, & VacancyDiffusion, & rho_dip_distance, & - v_cl, & !< climb velocity + v_cl, & !< climb velocity Gamma, & !< stacking fault energy tau, & - sigma_cl, & !< climb stress + sigma_cl, & !< climb stress b_d !< ratio of burgers vector to stacking fault width real(pReal), dimension(param(instance)%sum_N_sl) :: & dot_rho_dip_formation, & @@ -745,9 +677,9 @@ module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) call kinetics_slip(Mp,T,instance,of,dot_gamma_sl) dot%gamma_sl(:,of) = abs(dot_gamma_sl) - + rho_dip_distance_min = prm%CEdgeDipMinDistance*prm%b_sl - + slipState: do i = 1, prm%sum_N_sl tau = math_mul33xx33(Mp,prm%P_sl(1:3,1:3,i)) @@ -771,15 +703,15 @@ module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) else !@details: Refer: Argon & Moffat, Acta Metallurgica, Vol. 29, pg 293 to 299, 1981 sigma_cl = dot_product(prm%n0_sl(1:3,i),matmul(Mp,prm%n0_sl(1:3,i))) - if (prm%ExtendedDislocations) then + if (prm%ExtendedDislocations) then Gamma = prm%SFE_0K + prm%dSFE_dT * T b_d = 24.0_pReal*PI*(1.0_pReal - prm%nu)/(2.0_pReal + prm%nu)* Gamma/(prm%mu*prm%b_sl(i)) else - b_d = 1.0_pReal + b_d = 1.0_pReal endif v_cl = 2.0_pReal*prm%omega*b_d**2.0_pReal*exp(-prm%Qsd/(kB*T)) & * (exp(abs(sigma_cl)*prm%b_sl(i)**3.0_pReal/(kB*T)) - 1.0_pReal) - + dot_rho_dip_climb(i) = 4.0_pReal*v_cl*stt%rho_dip(i,of) & / (rho_dip_distance-rho_dip_distance_min(i)) endif @@ -801,12 +733,12 @@ module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) dot%f_tr(:,of) = f_unrotated*dot_gamma_tr end associate - + end subroutine plastic_dislotwin_dotState !-------------------------------------------------------------------------------------------------- -!> @brief calculates derived quantities from state +!> @brief Calculate derived quantities from state. !-------------------------------------------------------------------------------------------------- module subroutine plastic_dislotwin_dependentState(T,instance,of) @@ -815,7 +747,7 @@ module subroutine plastic_dislotwin_dependentState(T,instance,of) of real(pReal), intent(in) :: & T - + real(pReal) :: & sumf_twin,Gamma,sumf_trans real(pReal), dimension(param(instance)%sum_N_sl) :: & @@ -830,35 +762,35 @@ module subroutine plastic_dislotwin_dependentState(T,instance,of) f_over_t_tr real(pReal), dimension(:), allocatable :: & x0 - - + + associate(prm => param(instance),& stt => state(instance),& dst => dependentState(instance)) - + sumf_twin = sum(stt%f_tw(1:prm%sum_N_tw,of)) sumf_trans = sum(stt%f_tr(1:prm%sum_N_tr,of)) - + Gamma = prm%SFE_0K + prm%dSFE_dT * T - + !* rescaled volume fraction for topology f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,of)/prm%t_tw ! this is per system ... f_over_t_tr = sumf_trans/prm%t_tr ! but this not ! ToDo ...Physically correct, but naming could be adjusted - + inv_lambda_sl_sl = sqrt(matmul(prm%forestProjection, & stt%rho_mob(:,of)+stt%rho_dip(:,of)))/prm%CLambdaSlip - + if (prm%sum_N_tw > 0 .and. prm%sum_N_sl > 0) & inv_lambda_sl_tw = matmul(prm%h_sl_tw,f_over_t_tw)/(1.0_pReal-sumf_twin) - + inv_lambda_tw_tw = matmul(prm%h_tw_tw,f_over_t_tw)/(1.0_pReal-sumf_twin) - + if (prm%sum_N_tr > 0 .and. prm%sum_N_sl > 0) & inv_lambda_sl_tr = matmul(prm%h_sl_tr,f_over_t_tr)/(1.0_pReal-sumf_trans) - + inv_lambda_tr_tr = matmul(prm%h_tr_tr,f_over_t_tr)/(1.0_pReal-sumf_trans) - + if ((prm%sum_N_tw > 0) .or. (prm%sum_N_tr > 0)) then ! ToDo: better logic needed here dst%Lambda_sl(:,of) = prm%D & / (1.0_pReal+prm%D*(inv_lambda_sl_sl + inv_lambda_sl_tw + inv_lambda_sl_tr)) @@ -866,13 +798,13 @@ module subroutine plastic_dislotwin_dependentState(T,instance,of) dst%Lambda_sl(:,of) = prm%D & / (1.0_pReal+prm%D*inv_lambda_sl_sl) !!!!!! correct? endif - + dst%Lambda_tw(:,of) = prm%i_tw*prm%D/(1.0_pReal+prm%D*inv_lambda_tw_tw) dst%Lambda_tr(:,of) = prm%i_tr*prm%D/(1.0_pReal+prm%D*inv_lambda_tr_tr) - + !* threshold stress for dislocation motion dst%tau_pass(:,of) = prm%mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,of)+stt%rho_dip(:,of))) - + !* threshold stress for growing twin/martensite if(prm%sum_N_tw == prm%sum_N_sl) & dst%tau_hat_tw(:,of) = Gamma/(3.0_pReal*prm%b_tw) & @@ -880,66 +812,67 @@ module subroutine plastic_dislotwin_dependentState(T,instance,of) if(prm%sum_N_tr == prm%sum_N_sl) & dst%tau_hat_tr(:,of) = Gamma/(3.0_pReal*prm%b_tr) & + 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_sl) & ! slip burgers here correct? - + prm%h*prm%gamma_fcc_hex/ (3.0_pReal*prm%b_tr) - + + prm%h*prm%gamma_fcc_hex/ (3.0_pReal*prm%b_tr) + dst%V_tw(:,of) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,of)**2.0_pReal dst%V_tr(:,of) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,of)**2.0_pReal - - + + x0 = prm%mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the burgers vector for slip and is the same for twin and trans dst%tau_r_tw(:,of) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%xc_twin)+cos(pi/3.0_pReal)/x0) - + x0 = prm%mu*prm%b_tr**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu) ! ToDo: In the paper, this is the burgers vector for slip dst%tau_r_tr(:,of) = prm%mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%xc_trans)+cos(pi/3.0_pReal)/x0) - + end associate end subroutine plastic_dislotwin_dependentState !-------------------------------------------------------------------------------------------------- -!> @brief writes results to HDF5 output file +!> @brief Write results to HDF5 output file. !-------------------------------------------------------------------------------------------------- module subroutine plastic_dislotwin_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group + integer :: o - associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) + associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) + outputsLoop: do o = 1,size(prm%output) + select case(trim(prm%output(o))) - case (rho_mob_ID) - call results_writeDataset(group,stt%rho_mob,'rho_mob',& - 'mobile dislocation density','1/m²') - case (rho_dip_ID) - call results_writeDataset(group,stt%rho_dip,'rho_dip',& - 'dislocation dipole density''1/m²') - case (gamma_sl_ID) - call results_writeDataset(group,stt%gamma_sl,'gamma_sl',& - 'plastic shear','1') - case (Lambda_sl_ID) - call results_writeDataset(group,dst%Lambda_sl,'Lambda_sl',& - 'mean free path for slip','m') - case (tau_pass_ID) - call results_writeDataset(group,dst%tau_pass,'tau_pass',& - 'passing stress for slip','Pa') + case('rho_mob') + if(prm%sum_N_sl>0) call results_writeDataset(group,stt%rho_mob,'rho_mob',& + 'mobile dislocation density','1/m²') + case('rho_dip') + if(prm%sum_N_sl>0) call results_writeDataset(group,stt%rho_dip,'rho_dip',& + 'dislocation dipole density''1/m²') + case('gamma_sl') + if(prm%sum_N_sl>0) call results_writeDataset(group,stt%gamma_sl,'gamma_sl',& + 'plastic shear','1') + case('lambda_sl') + if(prm%sum_N_sl>0) call results_writeDataset(group,dst%Lambda_sl,'Lambda_sl',& + 'mean free path for slip','m') + case('tau_pass') + if(prm%sum_N_sl>0) call results_writeDataset(group,dst%tau_pass,'tau_pass',& + 'passing stress for slip','Pa') - case (f_tw_ID) - call results_writeDataset(group,stt%f_tw,'f_tw',& - 'twinned volume fraction','m³/m³') - case (Lambda_tw_ID) - call results_writeDataset(group,dst%Lambda_tw,'Lambda_tw',& - 'mean free path for twinning','m') - case (tau_hat_tw_ID) - call results_writeDataset(group,dst%tau_hat_tw,'tau_hat_tw',& - 'threshold stress for twinning','Pa') + case('f_tw') + if(prm%sum_N_tw>0) call results_writeDataset(group,stt%f_tw,'f_tw',& + 'twinned volume fraction','m³/m³') + case('lambda_tw') + if(prm%sum_N_tw>0) call results_writeDataset(group,dst%Lambda_tw,'Lambda_tw',& + 'mean free path for twinning','m') + case('tau_hat_tw') + if(prm%sum_N_tw>0) call results_writeDataset(group,dst%tau_hat_tw,'tau_hat_tw',& + 'threshold stress for twinning','Pa') + + case('f_tr') + if(prm%sum_N_tr>0) call results_writeDataset(group,stt%f_tr,'f_tr',& + 'martensite volume fraction','m³/m³') - case (f_tr_ID) - call results_writeDataset(group,stt%f_tr,'f_tr',& - 'martensite volume fraction','m³/m³') - end select enddo outputsLoop end associate @@ -948,8 +881,8 @@ end subroutine plastic_dislotwin_results !-------------------------------------------------------------------------------------------------- -!> @brief Shear rates on slip systems, their derivatives with respect to resolved stress and the -! resolved stresss +!> @brief Calculate shear rates on slip systems, their derivatives with respect to resolved +! stress, and the resolved stress. !> @details Derivatives and resolved stress are calculated only optionally. ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to ! have the optional arguments at the end @@ -964,7 +897,7 @@ pure subroutine kinetics_slip(Mp,T,instance,of, & integer, intent(in) :: & instance, & of - + real(pReal), dimension(param(instance)%sum_N_sl), intent(out) :: & dot_gamma_sl real(pReal), dimension(param(instance)%sum_N_sl), optional, intent(out) :: & @@ -972,7 +905,7 @@ pure subroutine kinetics_slip(Mp,T,instance,of, & tau_slip real(pReal), dimension(param(instance)%sum_N_sl) :: & ddot_gamma_dtau - + real(pReal), dimension(param(instance)%sum_N_sl) :: & tau, & stressRatio, & @@ -984,25 +917,25 @@ pure subroutine kinetics_slip(Mp,T,instance,of, & dV_run_inverse_dTau, & dV_dTau, & tau_eff !< effective resolved stress - integer :: i - + integer :: i + associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) - + do i = 1, prm%sum_N_sl tau(i) = math_mul33xx33(Mp,prm%P_sl(1:3,1:3,i)) enddo - + tau_eff = abs(tau)-dst%tau_pass(:,of) - + significantStress: where(tau_eff > tol_math_check) stressRatio = tau_eff/prm%tau_0 StressRatio_p = stressRatio** prm%p BoltzmannRatio = prm%Delta_F/(kB*T) v_wait_inverse = prm%v0**(-1.0_pReal) * exp(BoltzmannRatio*(1.0_pReal-StressRatio_p)** prm%q) v_run_inverse = prm%B/(tau_eff*prm%b_sl) - + dot_gamma_sl = sign(stt%rho_mob(:,of)*prm%b_sl/(v_wait_inverse+v_run_inverse),tau) - + dV_wait_inverse_dTau = -1.0_pReal * v_wait_inverse * prm%p * prm%q * BoltzmannRatio & * (stressRatio**(prm%p-1.0_pReal)) & * (1.0_pReal-StressRatio_p)**(prm%q-1.0_pReal) & @@ -1015,17 +948,21 @@ pure subroutine kinetics_slip(Mp,T,instance,of, & dot_gamma_sl = 0.0_pReal ddot_gamma_dtau = 0.0_pReal end where significantStress - + end associate - + if(present(ddot_gamma_dtau_slip)) ddot_gamma_dtau_slip = ddot_gamma_dtau if(present(tau_slip)) tau_slip = tau - + end subroutine kinetics_slip !-------------------------------------------------------------------------------------------------- -!> @brief calculates shear rates on twin systems +!> @brief Calculate shear rates on twin systems and their derivatives with respect to resolved +! stress. +!> @details Derivatives are calculated only optionally. +! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to +! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,of,& dot_gamma_twin,ddot_gamma_dtau_twin) @@ -1039,22 +976,22 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,of,& of real(pReal), dimension(param(instance)%sum_N_sl), intent(in) :: & dot_gamma_sl - + real(pReal), dimension(param(instance)%sum_N_tw), intent(out) :: & dot_gamma_twin real(pReal), dimension(param(instance)%sum_N_tw), optional, intent(out) :: & ddot_gamma_dtau_twin - + real, dimension(param(instance)%sum_N_tw) :: & tau, & Ndot0, & stressRatio_r, & ddot_gamma_dtau - + integer :: i,s1,s2 - + associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) - + do i = 1, prm%sum_N_tw tau(i) = math_mul33xx33(Mp,prm%P_tw(1:3,1:3,i)) isFCC: if (prm%fccTwinTransNucleation) then @@ -1072,7 +1009,7 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,of,& Ndot0=prm%dot_N_0_tw(i) endif isFCC enddo - + significantStress: where(tau > tol_math_check) StressRatio_r = (dst%tau_hat_tw(:,of)/tau)**prm%r dot_gamma_twin = prm%gamma_char * dst%V_tw(:,of) * Ndot0*exp(-StressRatio_r) @@ -1081,16 +1018,20 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,of,& dot_gamma_twin = 0.0_pReal ddot_gamma_dtau = 0.0_pReal end where significantStress - + end associate - + if(present(ddot_gamma_dtau_twin)) ddot_gamma_dtau_twin = ddot_gamma_dtau end subroutine kinetics_twin !-------------------------------------------------------------------------------------------------- -!> @brief calculates shear rates on twin systems +!> @brief Calculate shear rates on transformation systems and their derivatives with respect to +! resolved stress. +!> @details Derivatives are calculated only optionally. +! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to +! have the optional arguments at the end. !-------------------------------------------------------------------------------------------------- pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,& dot_gamma_tr,ddot_gamma_dtau_trans) @@ -1104,21 +1045,21 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,& of real(pReal), dimension(param(instance)%sum_N_sl), intent(in) :: & dot_gamma_sl - + real(pReal), dimension(param(instance)%sum_N_tr), intent(out) :: & dot_gamma_tr real(pReal), dimension(param(instance)%sum_N_tr), optional, intent(out) :: & ddot_gamma_dtau_trans - + real, dimension(param(instance)%sum_N_tr) :: & tau, & Ndot0, & stressRatio_s, & ddot_gamma_dtau - + integer :: i,s1,s2 associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) - + do i = 1, prm%sum_N_tr tau(i) = math_mul33xx33(Mp,prm%P_tr(1:3,1:3,i)) isFCC: if (prm%fccTwinTransNucleation) then @@ -1136,7 +1077,7 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,& Ndot0=prm%dot_N_0_tr(i) endif isFCC enddo - + significantStress: where(tau > tol_math_check) StressRatio_s = (dst%tau_hat_tr(:,of)/tau)**prm%s dot_gamma_tr = dst%V_tr(:,of) * Ndot0*exp(-StressRatio_s) @@ -1145,9 +1086,9 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,& dot_gamma_tr = 0.0_pReal ddot_gamma_dtau = 0.0_pReal end where significantStress - + end associate - + if(present(ddot_gamma_dtau_trans)) ddot_gamma_dtau_trans = ddot_gamma_dtau end subroutine kinetics_trans diff --git a/src/constitutive_plastic_kinehardening.f90 b/src/constitutive_plastic_kinehardening.f90 index 36cf20af7..65751989c 100644 --- a/src/constitutive_plastic_kinehardening.f90 +++ b/src/constitutive_plastic_kinehardening.f90 @@ -377,24 +377,21 @@ module subroutine plastic_kinehardening_results(instance,group) case('resistance') if(prm%totalNslip>0) call results_writeDataset(group,stt%crss,'xi_sl', & 'resistance against plastic slip','Pa') - case('backstress') ! ToDo: should be 'tau_back' if(prm%totalNslip>0) call results_writeDataset(group,stt%crss_back,'tau_back', & 'back stress against plastic slip','Pa') - case ('sense') - if(prm%totalNslip>0) call results_writeDataset(group,stt%sense,'sense_of_shear','tbd','1') - + if(prm%totalNslip>0) call results_writeDataset(group,stt%sense,'sense_of_shear', & + 'tbd','1') case ('chi0') - if(prm%totalNslip>0) call results_writeDataset(group,stt%chi0,'chi0','tbd','Pa') - + if(prm%totalNslip>0) call results_writeDataset(group,stt%chi0,'chi0', & + 'tbd','Pa') case ('gamma0') - if(prm%totalNslip>0) call results_writeDataset(group,stt%gamma0,'gamma0','tbd','1') - + if(prm%totalNslip>0) call results_writeDataset(group,stt%gamma0,'gamma0', & + 'tbd','1') case ('accumulatedshear') if(prm%totalNslip>0) call results_writeDataset(group,stt%accshear,'gamma_sl', & 'plastic shear','1') - end select enddo outputsLoop end associate