diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index ef4914db4..6f28fd7f3 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -110,6 +110,15 @@ submodule(phase) mechanical dLp_dFi !< derivative of Lp with respect to Fi end subroutine plastic_LpAndItsTangents + module subroutine plastic_KinematicJump(ph, en, Jump_occurr,deltaFp) + integer, intent(in) :: & + ph, & + en + logical , intent(out) :: & + Jump_occurr + real(pReal), dimension(3,3), intent(out) :: & + deltaFp + end subroutine plastic_KinematicJump module subroutine plastic_isotropic_result(ph,group) integer, intent(in) :: ph @@ -376,6 +385,13 @@ end subroutine mechanical_result !-------------------------------------------------------------------------------------------------- !> @brief calculation of stress (P) with time integration based on a residuum in Lp and !> intermediate acceleration of the Newton-Raphson correction +!> @modified by Satya and Achal +!> check for detour i.e. if twinning is possible (we are not going ahead in Lp loop consistency) +!> checking by calling the deltaFp subroutine that should return 4 things +!1) deltaFp= Cij (correspondance matrix) representing twinning shear and reorientation +!2) -(twin volume fraction) for each twin system to make it harder for twinned material point to twin again by any twin system +!3) jump in the last sampled volume fraction +!4) logical true if twinning possible false if not occurring !-------------------------------------------------------------------------------------------------- function integrateStress(F,Fp0,Fi0,Delta_t,ph,en) result(status) @@ -989,6 +1005,7 @@ end subroutine mechanical_forward !-------------------------------------------------------------------------------------------------- !> @brief calculate stress (P) +!> @modified by Satya and Achal !-------------------------------------------------------------------------------------------------- module function phase_mechanical_constitutive(Delta_t,co,ce) result(status) @@ -1001,8 +1018,8 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(status) real(pREAL) :: & formerStep integer :: & - ph, en, sizeDotState - logical :: todo + ph, en, sizeDotState, o, sd + logical :: todo, FpJumped real(pREAL) :: stepFrac,step real(pREAL), dimension(3,3) :: & Fp0, & @@ -1010,7 +1027,8 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(status) Lp0, & Li0, & F0, & - F + F, & + deltaFp real(pREAL), dimension(plasticState(material_ID_phase(co,ce))%sizeState) :: state0 @@ -1031,6 +1049,28 @@ module function phase_mechanical_constitutive(Delta_t,co,ce) result(status) todo = .true. cutbackLooping: do while (todo) + ! achal calling Kinematic DeltaFp here + !** Starting to implement changes for accommodating large shear and reorientation caused by twinning** + !if(.not. FpJumped .and. NiterationStressLp>1) then !Achal: Reason for this if statement? + call plastic_KinematicJump(ph, en, FpJumped,deltaFp) + !if(FpJumped) write(6,*) 'element jumped', en + !if(FpJumped) then + !Fp0 = matmul(deltaFp,phase_mechanical_Fp0(ph)%data(1:3,1:3,en)) + + o = plasticState(ph)%offsetDeltaState + sd = plasticState(ph)%sizeDeltaState + + !update current state by jump + !plasticState(ph)%state(o+1:o+sd,en) = plasticState(ph)%state(o+1:o+sd,en) & + !+ plasticState(ph)%deltaState(o+1:o+sd,en) + + !store jumped state as initial value for next iteration + !plasticState(ph)%state0(o+1:o+sd,en) = plasticState(ph)%state(o+1:o+sd,en) + + !store jumped state as initial value for for substate, partitioned state as well + + !endif + if (status == STATUS_OK) then formerStep = step stepFrac = stepFrac + step diff --git a/src/phase_mechanical_plastic.f90 b/src/phase_mechanical_plastic.f90 index 1f0629e36..b2da13bd2 100644 --- a/src/phase_mechanical_plastic.f90 +++ b/src/phase_mechanical_plastic.f90 @@ -407,7 +407,8 @@ module function plastic_deltaState(ph, en) result(status) status = STATUS_OK select case (mechanical_plasticity_type(ph)) - case (MECHANICAL_PLASTICITY_NONLOCAL,MECHANICAL_PLASTICITY_KINEHARDENING) + case (MECHANICAL_PLASTICITY_PHENOPOWERLAW, MECHANICAL_PLASTICITY_NONLOCAL,& + MECHANICAL_PLASTICITY_KINEHARDENING) Mp = matmul(matmul(transpose(phase_mechanical_Fi(ph)%data(1:3,1:3,en)),& phase_mechanical_Fi(ph)%data(1:3,1:3,en)),& diff --git a/src/phase_mechanical_plastic_phenopowerlaw.f90 b/src/phase_mechanical_plastic_phenopowerlaw.f90 index a80ee9315..c22176f90 100644 --- a/src/phase_mechanical_plastic_phenopowerlaw.f90 +++ b/src/phase_mechanical_plastic_phenopowerlaw.f90 @@ -312,7 +312,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) startIndex = endIndex + 1 endIndex = endIndex + 1 stt%frozen => plasticState(ph)%state(startIndex,:) - stt%frozen = 0.0_pReal-1.0_pReal + !stt%frozen = 0.0_pReal-1.0_pReal dlt%frozen => plasticState(ph)%deltaState(startIndex-o,:) plasticState(ph)%atol(startIndex:endIndex) = pl%get_asReal('atol_gamma',defaultVal=1.0e-6_pReal) @@ -424,6 +424,7 @@ module function phenopowerlaw_dotState(Mp,ph,en) result(dotState) call kinetics_sl(Mp,ph,en, dot_gamma_sl) call kinetics_tw(Mp,ph,en, dot_gamma_tw, fdot_twin) + dot_gamma_sl = abs(dot_gamma_sl) sumF = sum(stt%gamma_tw(:,en)/prm%gamma_char) @@ -444,6 +445,48 @@ module function phenopowerlaw_dotState(Mp,ph,en) result(dotState) end function phenopowerlaw_dotState +!-------------------------------------------------------------------------------------------------- +!> @brief calculates (instantaneous) incremental change of microstructure +!> Satya, Achal +!-------------------------------------------------------------------------------------------------- +module subroutine plastic_phenopowerlaw_deltaState(ph,en) + implicit none + + integer, intent(in)::& + ph, & + en + + logical :: & + twinJump + + integer :: & + twin_var + + real(pREAL), dimension(3,3) :: & + deltaFp + + ! These are updated at every strain increment. What should these initilizations be? + + associate(prm => param(ph), stt => state(ph), dlt => deltastate(ph)) + + twin_var = maxloc(stt%f_twin(:,en),dim=1) + + !write(6,*)'twin_var',twin_var + + call plastic_kinematic_deltaFp(ph,en,twinJump,deltaFp) + !write(6,*)'twinJump',twinJump + if(twinJump) then + !write(6,*)'el',en + dlt%f_twin(:,en) = 0.0_pReal - stt%f_twin(:,en) + dlt%fmc_twin(:,en) = 0.0_pReal - stt%fmc_twin(:,en) + dlt%frozen(en) = 1.0_pReal - stt%frozen(en) + write(6,*)'frozen',en,dlt%frozen(en),stt%frozen(en) + dlt%variant_twin(en) = twin_var !- stt%variant_twin(en) + endif + + end associate + + end subroutine plastic_phenopowerlaw_deltaState !-------------------------------------------------------------------------------------------------- !> @brief calculates instantaneous incremental change of kinematics and associated jump state @@ -487,14 +530,16 @@ module subroutine plastic_kinematic_deltaFp(ph,en,twinJump,deltaFp) Ability_Nucleation: if(stt%f_twin(twin_var,en)>(stt%fmc_twin(twin_var,en)+prm%checkstep(twin_var))) then stt%fmc_twin(twin_var,en) = stt%fmc_twin(twin_var,en)+prm%checkstep(twin_var) - Success_Nucleation: if (random <= stt%f_twin(twin_var,en)) then ! Instead of sum take max + write(6,*)'element twinned',en,'random',random,'variant',twin_var,'volume fraction', stt%f_twin(twin_var,en) twinJump = .true. deltaFp = prm%CorrespondenceMatrix(:,:,twin_var) - dlt%f_twin(:,en) = 0.0_pReal - stt%f_twin(:,en) - dlt%fmc_twin(:,en) = 0.0_pReal - stt%fmc_twin(:,en) - dlt%frozen(en) = 1.0_pReal - stt%frozen(en) - dlt%variant_twin(en) = twin_var !- stt%variant_twin(en) ! Achal LHS is real, RHS integer ! why this equation? + !write(6,*) twinJump + !dlt%f_twin(:,en) = 0.0_pReal - stt%f_twin(:,en) + !dlt%fmc_twin(:,en) = 0.0_pReal - stt%fmc_twin(:,en) + !dlt%frozen(en) = 1.0_pReal - stt%frozen(en) + !write(6,*)'frozen',en,dlt%frozen(en),stt%frozen(en) + !dlt%variant_twin(en) = twin_var !- stt%variant_twin(en) ! Achal LHS is real, RHS integer ! why this equation? endif Success_Nucleation @@ -504,30 +549,6 @@ module subroutine plastic_kinematic_deltaFp(ph,en,twinJump,deltaFp) end subroutine plastic_kinematic_deltaFp -!-------------------------------------------------------------------------------------------------- -!> @brief calculates (instantaneous) incremental change of microstructure -!> Satya, Achal -!-------------------------------------------------------------------------------------------------- -module subroutine plastic_phenopowerlaw_deltaState(ph,en) - implicit none - - integer, intent(in)::& - ph, & - en - - ! These are updated at every strain increment. What should these initilizations be? - - associate(dlt => deltastate(ph)) - - dlt%f_twin(:,en) = 0.0_pReal - dlt%fmc_twin(:,en) = 0.0_pReal - !dlt%variant_twin(en) = 0.0_pReal - !dlt%frozen(en) = 0.0_pReal - - end associate - - end subroutine plastic_phenopowerlaw_deltaState - !-------------------------------------------------------------------------------------------------- !> @brief Write results to HDF5 output file. !-------------------------------------------------------------------------------------------------- @@ -559,6 +580,17 @@ module subroutine plastic_phenopowerlaw_result(ph,group) call result_writeDataset(stt%gamma_tw,group,trim(prm%output(ou)), & 'twinning shear','1',prm%systems_tw) + case('f_twin') + call result_writeDataset(stt%f_twin,group,trim(prm%output(ou)), & + 'volume fraction','1',prm%systems_tw) !Achal + + case('variant_twin') + call result_writeDataset(stt%variant_twin,group,trim(prm%output(ou)), & + 'twin variant','1') !Achal + + case('fbinary_twin') + call result_writeDataset(stt%frozen,group,trim(prm%output(ou)), & + 'binary twin flag','1') !Achal end select end do @@ -643,14 +675,14 @@ pure subroutine kinetics_tw(Mp,ph,en,& integer :: i - associate(prm => param(ph), stt => state(ph)) + associate(prm => param(ph), stt => state(ph), dlt => state(ph)) tau_tw = [(math_tensordot(Mp,prm%P_tw(1:3,1:3,i)),i=1,prm%sum_N_tw)] - where(tau_tw > 0.0_pREAL .and. stt%frozen(en) < 0.9_pReal) ! Achal - dot_gamma_tw = (1.0_pREAL-sum(stt%gamma_tw(:,en)/prm%gamma_char)) & ! only twin in untwinned volume fraction - * prm%dot_gamma_0_tw*(tau_tw/stt%xi_tw(:,en))**prm%n_tw - fdot_twin = (0.05_pReal*(abs(tau_tw)/stt%xi_tw(:,en))**prm%n_tw)/prm%gamma_char ! Achal 0.05 is a constant + where(tau_tw > 0.0_pREAL .and. stt%frozen(en) < 0.9_pReal) ! Achal .and. stt%frozen(en) < 0.9_pReal + dot_gamma_tw = 0.0_pREAL !(1.0_pREAL-sum(stt%gamma_tw(:,en)/prm%gamma_char)) & ! only twin in untwinned volume fraction + !* prm%dot_gamma_0_tw*(tau_tw/stt%xi_tw(:,en))**prm%n_tw + fdot_twin = (0.005_pReal*(abs(tau_tw)/stt%xi_tw(:,en))**prm%n_tw)/prm%gamma_char ! Achal 0.005 is reference slip rate else where dot_gamma_tw = 0.0_pREAL fdot_twin = 0.0_pREAL