From 8259cb3cdc194a55b0fd6ff8d8373b6a1ae75364 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 21 Jul 2021 15:03:28 +0200 Subject: [PATCH] notation as in paper --- ...phase_mechanical_plastic_kinehardening.f90 | 12 +- ...phase_mechanical_plastic_phenopowerlaw.f90 | 185 ++++++++---------- 2 files changed, 93 insertions(+), 104 deletions(-) diff --git a/src/phase_mechanical_plastic_kinehardening.f90 b/src/phase_mechanical_plastic_kinehardening.f90 index 71799b8f0..44670163c 100644 --- a/src/phase_mechanical_plastic_kinehardening.f90 +++ b/src/phase_mechanical_plastic_kinehardening.f90 @@ -18,6 +18,8 @@ submodule(phase:plastic) kinehardening h_inf_b, & !< asymptotic hardening rate of back stress for each slip xi_inf_f, & xi_inf_b + real(pReal), allocatable, dimension(:,:) :: & + h_sl_sl !< slip resistance from slip activity real(pReal), allocatable, dimension(:,:,:) :: & P, & nonSchmid_pos, & @@ -122,9 +124,9 @@ module function plastic_kinehardening_init() result(myPlasticity) prm%nonSchmid_pos = prm%P prm%nonSchmid_neg = prm%P endif - prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(N_sl, & - pl%get_as1dFloat('h_sl-sl'), & - phase%get_asString('lattice')) + prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl, & + pl%get_as1dFloat('h_sl-sl'), & + phase%get_asString('lattice')) xi_0 = pl%get_as1dFloat('xi_0', requiredSize=size(N_sl)) prm%xi_inf_f = pl%get_as1dFloat('xi_inf_f', requiredSize=size(N_sl)) @@ -158,7 +160,7 @@ module function plastic_kinehardening_init() result(myPlasticity) else slipActive xi_0 = emptyRealArray allocate(prm%xi_inf_f,prm%xi_inf_b,prm%h_0_f,prm%h_inf_f,prm%h_0_b,prm%h_inf_b,source=emptyRealArray) - allocate(prm%interaction_SlipSlip(0,0)) + allocate(prm%h_sl_sl(0,0)) endif slipActive !-------------------------------------------------------------------------------------------------- @@ -288,7 +290,7 @@ module subroutine plastic_kinehardening_dotState(Mp,ph,en) sumGamma = sum(stt%accshear(:,en)) - dot%crss(:,en) = matmul(prm%interaction_SlipSlip,dot%accshear(:,en)) & + dot%crss(:,en) = matmul(prm%h_sl_sl,dot%accshear(:,en)) & * ( prm%h_inf_f & + (prm%h_0_f - prm%h_inf_f + prm%h_0_f*prm%h_inf_f*sumGamma/prm%xi_inf_f) & * exp(-sumGamma*prm%h_0_f/prm%xi_inf_f) & diff --git a/src/phase_mechanical_plastic_phenopowerlaw.f90 b/src/phase_mechanical_plastic_phenopowerlaw.f90 index 45e7a5f14..d82f11908 100644 --- a/src/phase_mechanical_plastic_phenopowerlaw.f90 +++ b/src/phase_mechanical_plastic_phenopowerlaw.f90 @@ -46,10 +46,10 @@ submodule(phase:plastic) phenopowerlaw type :: tPhenopowerlawState real(pReal), pointer, dimension(:,:) :: & - xi_slip, & - xi_twin, & - gamma_slip, & - gamma_twin + xi_sl, & + xi_tw, & + gamma_sl, & + gamma_tw end type tPhenopowerlawState !-------------------------------------------------------------------------------------------------- @@ -235,30 +235,30 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) ! state aliases and initialization startIndex = 1 endIndex = prm%sum_N_sl - stt%xi_slip => plasticState(ph)%state (startIndex:endIndex,:) - stt%xi_slip = spread(xi_0_sl, 2, Nmembers) - dot%xi_slip => plasticState(ph)%dotState(startIndex:endIndex,:) + stt%xi_sl => plasticState(ph)%state (startIndex:endIndex,:) + stt%xi_sl = spread(xi_0_sl, 2, Nmembers) + dot%xi_sl => plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_tw - stt%xi_twin => plasticState(ph)%state (startIndex:endIndex,:) - stt%xi_twin = spread(xi_0_tw, 2, Nmembers) - dot%xi_twin => plasticState(ph)%dotState(startIndex:endIndex,:) + stt%xi_tw => plasticState(ph)%state (startIndex:endIndex,:) + stt%xi_tw = spread(xi_0_tw, 2, Nmembers) + dot%xi_tw => plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_sl - stt%gamma_slip => plasticState(ph)%state (startIndex:endIndex,:) - dot%gamma_slip => plasticState(ph)%dotState(startIndex:endIndex,:) + stt%gamma_sl => plasticState(ph)%state (startIndex:endIndex,:) + dot%gamma_sl => plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma' startIndex = endIndex + 1 endIndex = endIndex + prm%sum_N_tw - stt%gamma_twin => plasticState(ph)%state (startIndex:endIndex,:) - dot%gamma_twin => plasticState(ph)%dotState(startIndex:endIndex,:) + stt%gamma_tw => plasticState(ph)%state (startIndex:endIndex,:) + dot%gamma_tw => plasticState(ph)%dotState(startIndex:endIndex,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal) end associate @@ -293,28 +293,28 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en) integer :: & i,k,l,m,n real(pReal), dimension(param(ph)%sum_N_sl) :: & - gdot_slip_pos,gdot_slip_neg, & + gdot_sl_pos,gdot_sl_neg, & dgdot_dtauslip_pos,dgdot_dtauslip_neg real(pReal), dimension(param(ph)%sum_N_tw) :: & - gdot_twin,dgdot_dtautwin + gdot_tw,dgdot_dtautwin Lp = 0.0_pReal dLp_dMp = 0.0_pReal associate(prm => param(ph)) - call kinetics_slip(Mp,ph,en,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg) + call kinetics_sl(Mp,ph,en,gdot_sl_pos,gdot_sl_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg) slipSystems: do i = 1, prm%sum_N_sl - Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%P_sl(1:3,1:3,i) + Lp = Lp + (gdot_sl_pos(i)+gdot_sl_neg(i))*prm%P_sl(1:3,1:3,i) 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) & + dgdot_dtauslip_pos(i) * prm%P_sl(k,l,i) * prm%nonSchmid_pos(m,n,i) & + dgdot_dtauslip_neg(i) * prm%P_sl(k,l,i) * prm%nonSchmid_neg(m,n,i) enddo slipSystems - call kinetics_twin(Mp,ph,en,gdot_twin,dgdot_dtautwin) + call kinetics_tw(Mp,ph,en,gdot_tw,dgdot_dtautwin) twinSystems: do i = 1, prm%sum_N_tw - Lp = Lp + gdot_twin(i)*prm%P_tw(1:3,1:3,i) + Lp = Lp + gdot_tw(i)*prm%P_tw(1:3,1:3,i) 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) & + dgdot_dtautwin(i)*prm%P_tw(k,l,i)*prm%P_tw(m,n,i) @@ -337,46 +337,33 @@ module subroutine phenopowerlaw_dotState(Mp,ph,en) en real(pReal) :: & - c_SlipSlip,c_TwinSlip,c_TwinTwin, & - xi_slip_sat_offset,& - sumGamma,sumF + xi_sl_sat_offset,& + sumF real(pReal), dimension(param(ph)%sum_N_sl) :: & - left_SlipSlip,right_SlipSlip, & - gdot_slip_pos,gdot_slip_neg + gdot_sl_pos,gdot_sl_neg, & + right_SlipSlip - associate(prm => param(ph), stt => state(ph), & - dot => dotState(ph)) - sumGamma = sum(stt%gamma_slip(:,en)) - sumF = sum(stt%gamma_twin(:,en)/prm%gamma_char) + associate(prm => param(ph), stt => state(ph), dot => dotState(ph)) -!-------------------------------------------------------------------------------------------------- -! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices - c_SlipSlip = prm%h_0_sl_sl * (1.0_pReal + prm%c_1*sumF** prm%c_2) - c_TwinSlip = prm%h_0_tw_sl * sumGamma**prm%c_3 - c_TwinTwin = prm%h_0_tw_tw * sumF**prm%c_4 + call kinetics_sl(Mp,ph,en,gdot_sl_pos,gdot_sl_neg) + dot%gamma_sl(:,en) = abs(gdot_sl_pos+gdot_sl_neg) + call kinetics_tw(Mp,ph,en,dot%gamma_tw(:,en)) -!-------------------------------------------------------------------------------------------------- -! calculate left and right vectors - left_SlipSlip = 1.0_pReal + prm%h_int - xi_slip_sat_offset = prm%f_sat_sl_tw*sqrt(sumF) - right_SlipSlip = sign(abs(1.0_pReal-stt%xi_slip(:,en) / (prm%xi_inf_sl+xi_slip_sat_offset)) **prm%a_sl, & - 1.0_pReal-stt%xi_slip(:,en) / (prm%xi_inf_sl+xi_slip_sat_offset)) -!-------------------------------------------------------------------------------------------------- -! shear rates - call kinetics_slip(Mp,ph,en,gdot_slip_pos,gdot_slip_neg) - dot%gamma_slip(:,en) = abs(gdot_slip_pos+gdot_slip_neg) - call kinetics_twin(Mp,ph,en,dot%gamma_twin(:,en)) + sumF = sum(stt%gamma_tw(:,en)/prm%gamma_char) + xi_sl_sat_offset = prm%f_sat_sl_tw*sqrt(sumF) + right_SlipSlip = sign(abs(1.0_pReal-stt%xi_sl(:,en) / (prm%xi_inf_sl+xi_sl_sat_offset)) **prm%a_sl, & + 1.0_pReal-stt%xi_sl(:,en) / (prm%xi_inf_sl+xi_sl_sat_offset)) -!-------------------------------------------------------------------------------------------------- -! hardening - dot%xi_slip(:,en) = c_SlipSlip * left_SlipSlip * & - matmul(prm%h_sl_sl,dot%gamma_slip(:,en)*right_SlipSlip) & - + matmul(prm%h_sl_tw,dot%gamma_twin(:,en)) + dot%xi_sl(:,en) = prm%h_0_sl_sl * (1.0_pReal + prm%c_1*sumF** prm%c_2) * (1.0_pReal + prm%h_int) & + * matmul(prm%h_sl_sl,dot%gamma_sl(:,en)*right_SlipSlip) & + + matmul(prm%h_sl_tw,dot%gamma_tw(:,en)) + + dot%xi_tw(:,en) = prm%h_0_tw_sl * sum(stt%gamma_sl(:,en))**prm%c_3 & + * matmul(prm%h_tw_sl,dot%gamma_sl(:,en)) & + + prm%h_0_tw_tw * sumF**prm%c_4 * matmul(prm%h_tw_tw,dot%gamma_tw(:,en)) - dot%xi_twin(:,en) = c_TwinSlip * matmul(prm%h_tw_sl,dot%gamma_slip(:,en)) & - + c_TwinTwin * matmul(prm%h_tw_tw,dot%gamma_twin(:,en)) end associate end subroutine phenopowerlaw_dotState @@ -397,17 +384,17 @@ module subroutine plastic_phenopowerlaw_results(ph,group) select case(trim(prm%output(o))) case('xi_sl') - if(prm%sum_N_sl>0) call results_writeDataset(stt%xi_slip,group,trim(prm%output(o)), & + if(prm%sum_N_sl>0) call results_writeDataset(stt%xi_sl,group,trim(prm%output(o)), & 'resistance against plastic slip','Pa') case('gamma_sl') - if(prm%sum_N_sl>0) call results_writeDataset(stt%gamma_slip,group,trim(prm%output(o)), & + if(prm%sum_N_sl>0) call results_writeDataset(stt%gamma_sl,group,trim(prm%output(o)), & 'plastic shear','1') case('xi_tw') - if(prm%sum_N_tw>0) call results_writeDataset(stt%xi_twin,group,trim(prm%output(o)), & + if(prm%sum_N_tw>0) call results_writeDataset(stt%xi_tw,group,trim(prm%output(o)), & 'resistance against twinning','Pa') case('gamma_tw') - if(prm%sum_N_tw>0) call results_writeDataset(stt%gamma_twin,group,trim(prm%output(o)), & + if(prm%sum_N_tw>0) call results_writeDataset(stt%gamma_tw,group,trim(prm%output(o)), & 'twinning shear','1') end select @@ -424,8 +411,8 @@ end subroutine plastic_phenopowerlaw_results ! 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_slip(Mp,ph,en, & - gdot_slip_pos,gdot_slip_neg,dgdot_dtau_slip_pos,dgdot_dtau_slip_neg) +pure subroutine kinetics_sl(Mp,ph,en, & + gdot_sl_pos,gdot_sl_neg,dgdot_dtau_sl_pos,dgdot_dtau_sl_neg) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -434,56 +421,56 @@ pure subroutine kinetics_slip(Mp,ph,en, & en real(pReal), intent(out), dimension(param(ph)%sum_N_sl) :: & - gdot_slip_pos, & - gdot_slip_neg + gdot_sl_pos, & + gdot_sl_neg real(pReal), intent(out), optional, dimension(param(ph)%sum_N_sl) :: & - dgdot_dtau_slip_pos, & - dgdot_dtau_slip_neg + dgdot_dtau_sl_pos, & + dgdot_dtau_sl_neg real(pReal), dimension(param(ph)%sum_N_sl) :: & - tau_slip_pos, & - tau_slip_neg + tau_sl_pos, & + tau_sl_neg integer :: i associate(prm => param(ph), stt => state(ph)) do i = 1, prm%sum_N_sl - tau_slip_pos(i) = math_tensordot(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - tau_slip_neg(i) = merge(math_tensordot(Mp,prm%nonSchmid_neg(1:3,1:3,i)), & - 0.0_pReal, prm%nonSchmidActive) + tau_sl_pos(i) = math_tensordot(Mp,prm%nonSchmid_pos(1:3,1:3,i)) + tau_sl_neg(i) = merge(math_tensordot(Mp,prm%nonSchmid_neg(1:3,1:3,i)), & + 0.0_pReal, prm%nonSchmidActive) enddo - where(dNeq0(tau_slip_pos)) - gdot_slip_pos = prm%dot_gamma_0_sl * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active - * sign(abs(tau_slip_pos/stt%xi_slip(:,en))**prm%n_sl, tau_slip_pos) + where(dNeq0(tau_sl_pos)) + gdot_sl_pos = prm%dot_gamma_0_sl * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) & ! 1/2 if non-Schmid active + * sign(abs(tau_sl_pos/stt%xi_sl(:,en))**prm%n_sl, tau_sl_pos) else where - gdot_slip_pos = 0.0_pReal + gdot_sl_pos = 0.0_pReal end where - where(dNeq0(tau_slip_neg)) - gdot_slip_neg = prm%dot_gamma_0_sl * 0.5_pReal & ! only used if non-Schmid active, always 1/2 - * sign(abs(tau_slip_neg/stt%xi_slip(:,en))**prm%n_sl, tau_slip_neg) + where(dNeq0(tau_sl_neg)) + gdot_sl_neg = prm%dot_gamma_0_sl * 0.5_pReal & ! only used if non-Schmid active, always 1/2 + * sign(abs(tau_sl_neg/stt%xi_sl(:,en))**prm%n_sl, tau_sl_neg) else where - gdot_slip_neg = 0.0_pReal + gdot_sl_neg = 0.0_pReal end where - if (present(dgdot_dtau_slip_pos)) then - where(dNeq0(gdot_slip_pos)) - dgdot_dtau_slip_pos = gdot_slip_pos*prm%n_sl/tau_slip_pos + if (present(dgdot_dtau_sl_pos)) then + where(dNeq0(gdot_sl_pos)) + dgdot_dtau_sl_pos = gdot_sl_pos*prm%n_sl/tau_sl_pos else where - dgdot_dtau_slip_pos = 0.0_pReal + dgdot_dtau_sl_pos = 0.0_pReal end where endif - if (present(dgdot_dtau_slip_neg)) then - where(dNeq0(gdot_slip_neg)) - dgdot_dtau_slip_neg = gdot_slip_neg*prm%n_sl/tau_slip_neg + if (present(dgdot_dtau_sl_neg)) then + where(dNeq0(gdot_sl_neg)) + dgdot_dtau_sl_neg = gdot_sl_neg*prm%n_sl/tau_sl_neg else where - dgdot_dtau_slip_neg = 0.0_pReal + dgdot_dtau_sl_neg = 0.0_pReal end where endif end associate -end subroutine kinetics_slip +end subroutine kinetics_sl !-------------------------------------------------------------------------------------------------- @@ -493,8 +480,8 @@ end subroutine kinetics_slip ! 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,ph,en,& - gdot_twin,dgdot_dtau_twin) +pure subroutine kinetics_tw(Mp,ph,en,& + gdot_tw,dgdot_dtau_tw) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -503,37 +490,37 @@ pure subroutine kinetics_twin(Mp,ph,en,& en real(pReal), dimension(param(ph)%sum_N_tw), intent(out) :: & - gdot_twin + gdot_tw real(pReal), dimension(param(ph)%sum_N_tw), intent(out), optional :: & - dgdot_dtau_twin + dgdot_dtau_tw real(pReal), dimension(param(ph)%sum_N_tw) :: & - tau_twin + tau_tw integer :: i associate(prm => param(ph), stt => state(ph)) do i = 1, prm%sum_N_tw - tau_twin(i) = math_tensordot(Mp,prm%P_tw(1:3,1:3,i)) + tau_tw(i) = math_tensordot(Mp,prm%P_tw(1:3,1:3,i)) enddo - where(tau_twin > 0.0_pReal) - gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,en)/prm%gamma_char)) & ! only twin in untwinned volume fraction - * prm%dot_gamma_0_tw*(abs(tau_twin)/stt%xi_twin(:,en))**prm%n_tw + where(tau_tw > 0.0_pReal) + gdot_tw = (1.0_pReal-sum(stt%gamma_tw(:,en)/prm%gamma_char)) & ! only twin in untwinned volume fraction + * prm%dot_gamma_0_tw*(abs(tau_tw)/stt%xi_tw(:,en))**prm%n_tw else where - gdot_twin = 0.0_pReal + gdot_tw = 0.0_pReal end where - if (present(dgdot_dtau_twin)) then - where(dNeq0(gdot_twin)) - dgdot_dtau_twin = gdot_twin*prm%n_tw/tau_twin + if (present(dgdot_dtau_tw)) then + where(dNeq0(gdot_tw)) + dgdot_dtau_tw = gdot_tw*prm%n_tw/tau_tw else where - dgdot_dtau_twin = 0.0_pReal + dgdot_dtau_tw = 0.0_pReal end where endif end associate -end subroutine kinetics_twin +end subroutine kinetics_tw end submodule phenopowerlaw