diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index 31bddd141..c22a37a0f 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -36,14 +36,14 @@ #include "kinematics_cleavage_opening.f90" #include "kinematics_slipplane_opening.f90" #include "kinematics_thermal_expansion.f90" -#include "plastic_none.f90" -#include "plastic_isotropic.f90" -#include "plastic_phenopowerlaw.f90" -#include "plastic_kinematichardening.f90" -#include "plastic_dislotwin.f90" -#include "plastic_disloUCLA.f90" -#include "plastic_nonlocal.f90" #include "constitutive.f90" +#include "constitutive_plastic_none.f90" +#include "constitutive_plastic_isotropic.f90" +#include "constitutive_plastic_phenopowerlaw.f90" +#include "constitutive_plastic_kinehardening.f90" +#include "constitutive_plastic_dislotwin.f90" +#include "constitutive_plastic_disloUCLA.f90" +#include "constitutive_plastic_nonlocal.f90" #include "crystallite.f90" #include "thermal_isothermal.f90" #include "thermal_adiabatic.f90" diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 7ddceb098..681691c9e 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -4,7 +4,9 @@ !> @brief elasticity, plasticity, internal microstructure state !-------------------------------------------------------------------------------------------------- module constitutive + use prec use math + use rotations use debug use numerics use IO @@ -13,13 +15,6 @@ module constitutive use results use lattice use discretization - use plastic_none - use plastic_isotropic - use plastic_phenopowerlaw - use plastic_kinehardening - use plastic_dislotwin - use plastic_disloucla - use plastic_nonlocal use geometry_plastic_nonlocal use source_thermal_dissipation use source_thermal_externalheat @@ -38,7 +33,284 @@ module constitutive constitutive_plasticity_maxSizeDotState, & constitutive_source_maxSizeDotState + interface + + module subroutine plastic_none_init + end subroutine plastic_none_init + + module subroutine plastic_isotropic_init + end subroutine plastic_isotropic_init + + module subroutine plastic_phenopowerlaw_init + end subroutine plastic_phenopowerlaw_init + + module subroutine plastic_kinehardening_init + end subroutine plastic_kinehardening_init + + module subroutine plastic_dislotwin_init + end subroutine plastic_dislotwin_init + + module subroutine plastic_disloUCLA_init + end subroutine plastic_disloUCLA_init + + module subroutine plastic_nonlocal_init + end subroutine plastic_nonlocal_init + + + module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) + real(pReal), dimension(3,3), intent(out) :: & + Lp !< plastic velocity gradient + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp !< derivative of Lp with respect to the Mandel stress + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + end subroutine plastic_isotropic_LpAndItsTangent + + pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) + real(pReal), dimension(3,3), intent(out) :: & + Lp !< plastic velocity gradient + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp !< derivative of Lp with respect to the Mandel stress + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + end subroutine plastic_phenopowerlaw_LpAndItsTangent + + pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) + real(pReal), dimension(3,3), intent(out) :: & + Lp !< plastic velocity gradient + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp !< derivative of Lp with respect to the Mandel stress + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + end subroutine plastic_kinehardening_LpAndItsTangent + + module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) + real(pReal), dimension(3,3), intent(out) :: & + Lp !< plastic velocity gradient + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp !< derivative of Lp with respect to the Mandel stress + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T + integer, intent(in) :: & + instance, & + of + end subroutine plastic_dislotwin_LpAndItsTangent + + pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) + real(pReal), dimension(3,3), intent(out) :: & + Lp !< plastic velocity gradient + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp !< derivative of Lp with respect to the Mandel stress + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T + integer, intent(in) :: & + instance, & + of + end subroutine plastic_disloUCLA_LpAndItsTangent + + module subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, & + Mp, Temperature, volume, ip, el) + real(pReal), dimension(3,3), intent(out) :: & + Lp !< plastic velocity gradient + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp !< derivative of Lp with respect to the Mandel stress + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + Temperature, & + volume + integer, intent(in) :: & + ip, & !< current integration point + el !< current element number + end subroutine plastic_nonlocal_LpAndItsTangent + + + module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) + real(pReal), dimension(3,3), intent(out) :: & + Li !< inleastic velocity gradient + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLi_dMi !< derivative of Li with respect to Mandel stress + + real(pReal), dimension(3,3), intent(in) :: & + Mi !< Mandel stress + integer, intent(in) :: & + instance, & + of + end subroutine plastic_isotropic_LiAndItsTangent + + + module subroutine plastic_isotropic_dotState(Mp,instance,of) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + end subroutine plastic_isotropic_dotState + + module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + end subroutine plastic_phenopowerlaw_dotState + + module subroutine plastic_kinehardening_dotState(Mp,instance,of) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + end subroutine plastic_kinehardening_dotState + + module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T + integer, intent(in) :: & + instance, & + of + end subroutine plastic_dislotwin_dotState + + module subroutine plastic_disloUCLA_dotState(Mp,T,instance,of) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T + integer, intent(in) :: & + instance, & + of + end subroutine plastic_disloUCLA_dotState + + module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & + timestep,ip,el) + integer, intent(in) :: & + ip, & !< current integration point + el !< current element number + real(pReal), intent(in) :: & + Temperature, & !< temperature + timestep !< substepped crystallite time increment + real(pReal), dimension(3,3), intent(in) ::& + Mp !< MandelStress + real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: & + Fe, & !< elastic deformation gradient + Fp !< plastic deformation gradient + end subroutine plastic_nonlocal_dotState + + + module subroutine plastic_dislotwin_dependentState(T,instance,of) + integer, intent(in) :: & + instance, & + of + real(pReal), intent(in) :: & + T + end subroutine plastic_dislotwin_dependentState + + module subroutine plastic_disloUCLA_dependentState(instance,of) + integer, intent(in) :: & + instance, & + of + end subroutine plastic_disloUCLA_dependentState + + module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) + integer, intent(in) :: & + ip, & + el + real(pReal), dimension(3,3), intent(in) :: & + Fe, & + Fp + end subroutine plastic_nonlocal_dependentState + + + module subroutine plastic_kinehardening_deltaState(Mp,instance,of) + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + end subroutine plastic_kinehardening_deltaState + + module subroutine plastic_nonlocal_deltaState(Mp,ip,el) + integer, intent(in) :: & + ip, & + el + real(pReal), dimension(3,3), intent(in) :: & + Mp + end subroutine plastic_nonlocal_deltaState + + + 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 + end function plastic_dislotwin_homogenizedC + + module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) + integer, intent(in) :: & + i, & + e + type(rotation), dimension(1,discretization_nIP,discretization_nElem), intent(in) :: & + orientation !< crystal orientation + end subroutine plastic_nonlocal_updateCompatibility + + + module subroutine plastic_isotropic_results(instance,group) + integer, intent(in) :: instance + character(len=*), intent(in) :: group + end subroutine plastic_isotropic_results + + module subroutine plastic_phenopowerlaw_results(instance,group) + integer, intent(in) :: instance + character(len=*), intent(in) :: group + end subroutine plastic_phenopowerlaw_results + + module subroutine plastic_kinehardening_results(instance,group) + integer, intent(in) :: instance + character(len=*), intent(in) :: group + end subroutine plastic_kinehardening_results + + module subroutine plastic_dislotwin_results(instance,group) + integer, intent(in) :: instance + character(len=*), intent(in) :: group + end subroutine plastic_dislotwin_results + + module subroutine plastic_disloUCLA_results(instance,group) + integer, intent(in) :: instance + character(len=*), intent(in) :: group + end subroutine plastic_disloUCLA_results + + module subroutine plastic_nonlocal_results(instance,group) + integer, intent(in) :: instance + character(len=*), intent(in) :: group + end subroutine plastic_nonlocal_results + + end interface + public :: & + plastic_nonlocal_updateCompatibility, & constitutive_init, & constitutive_homogenizedC, & constitutive_microstructure, & @@ -179,8 +451,7 @@ end subroutine constitutive_microstructure ! Mp in, dLp_dMp out !-------------------------------------------------------------------------------------------------- subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & - S, Fi, ipc, ip, el) - + S, Fi, ipc, ip, el) integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point diff --git a/src/plastic_disloUCLA.f90 b/src/constitutive_plastic_disloUCLA.f90 similarity index 96% rename from src/plastic_disloUCLA.f90 rename to src/constitutive_plastic_disloUCLA.f90 index d1291d853..8247d20f5 100644 --- a/src/plastic_disloUCLA.f90 +++ b/src/constitutive_plastic_disloUCLA.f90 @@ -5,21 +5,9 @@ !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief crystal plasticity model for bcc metals, especially Tungsten !-------------------------------------------------------------------------------------------------- -module plastic_disloUCLA - use prec - use debug - use math - use IO - use material - use config - use lattice - use discretization - use results +submodule(constitutive) plastic_disloUCLA - implicit none - private - - real(pReal), parameter, private :: & + real(pReal), parameter :: & kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin enum, bind(c) @@ -33,7 +21,7 @@ module plastic_disloUCLA tau_pass_ID end enum - type, private :: tParameters + type :: tParameters real(pReal) :: & aTol_rho, & D, & !< grain size @@ -75,14 +63,14 @@ module plastic_disloUCLA dipoleFormation !< flag indicating consideration of dipole formation end type !< container type for internal constitutive parameters - type, private :: tDisloUCLAState + type :: tDisloUCLAState real(pReal), dimension(:,:), pointer :: & rho_mob, & rho_dip, & gamma_sl end type tDisloUCLAState - type, private :: tDisloUCLAdependentState + type :: tDisloUCLAdependentState real(pReal), dimension(:,:), allocatable :: & Lambda_sl, & threshold_stress @@ -90,20 +78,11 @@ module plastic_disloUCLA !-------------------------------------------------------------------------------------------------- ! containers for parameters and state - type(tParameters), allocatable, dimension(:), private :: param - type(tDisloUCLAState), allocatable, dimension(:), private :: & + type(tParameters), allocatable, dimension(:) :: param + type(tDisloUCLAState), allocatable, dimension(:) :: & dotState, & state - type(tDisloUCLAdependentState), allocatable, dimension(:), private :: dependentState - - public :: & - plastic_disloUCLA_init, & - plastic_disloUCLA_dependentState, & - plastic_disloUCLA_LpAndItsTangent, & - plastic_disloUCLA_dotState, & - plastic_disloUCLA_results - private :: & - kinetics + type(tDisloUCLAdependentState), allocatable, dimension(:) :: dependentState contains @@ -112,7 +91,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine plastic_disloUCLA_init() +module subroutine plastic_disloUCLA_init integer :: & Ninstance, & @@ -328,7 +307,7 @@ end subroutine plastic_disloUCLA_init !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- -pure subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp, & +pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp, & Mp,T,instance,of) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient @@ -371,7 +350,7 @@ end subroutine plastic_disloUCLA_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine plastic_disloUCLA_dotState(Mp,T,instance,of) +module subroutine plastic_disloUCLA_dotState(Mp,T,instance,of) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -431,7 +410,7 @@ end subroutine plastic_disloUCLA_dotState !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -subroutine plastic_disloUCLA_dependentState(instance,of) +module subroutine plastic_disloUCLA_dependentState(instance,of) integer, intent(in) :: & instance, & @@ -457,7 +436,7 @@ end subroutine plastic_disloUCLA_dependentState !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- -subroutine plastic_disloUCLA_results(instance,group) +module subroutine plastic_disloUCLA_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group @@ -604,4 +583,4 @@ pure subroutine kinetics(Mp,T,instance,of, & end subroutine kinetics -end module plastic_disloUCLA +end submodule plastic_disloUCLA diff --git a/src/constitutive_plastic_dislotwin.f90 b/src/constitutive_plastic_dislotwin.f90 new file mode 100644 index 000000000..068323f3e --- /dev/null +++ b/src/constitutive_plastic_dislotwin.f90 @@ -0,0 +1,1156 @@ +!-------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH +!> @author Su Leen Wong, Max-Planck-Institut für Eisenforschung GmbH +!> @author Nan Jia, Max-Planck-Institut für Eisenforschung GmbH +!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH +!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH +!> @brief material subroutine incoprorating dislocation and twinning physics +!> @details to be done +!-------------------------------------------------------------------------------------------------- +submodule(constitutive) plastic_dislotwin + + 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 + D, & !< grain size + p_sb, & !< p-exponent in shear band velocity + q_sb, & !< q-exponent in shear band velocity + CEdgeDipMinDistance, & !< + i_tw, & !< + tau_0, & !< strength due to elements in solid solution + L_tw, & !< Length of twin nuclei in Burgers vectors + L_tr, & !< Length of trans nuclei in Burgers vectors + xc_twin, & !< critical distance for formation of twin nucleus + xc_trans, & !< critical distance for formation of trans nucleus + V_cs, & !< cross slip volume + sbResistance, & !< value for shearband resistance (might become an internal state variable at some point) + sbVelocity, & !< value for shearband velocity_0 + sbQedge, & !< activation energy for shear bands + SFE_0K, & !< stacking fault energy at zero K + dSFE_dT, & !< temperature dependance of stacking fault energy + aTol_rho, & !< absolute tolerance for integration of dislocation density + aTol_f_tw, & !< absolute tolerance for integration of twin volume fraction + 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 :: & + 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 + b_tw, & !< absolute length of burgers vector [m] for each twin system + b_tr, & !< absolute length of burgers vector [m] for each transformation system + Delta_F,& !< activation energy for glide [J] for each slip system + v0, & !< dislocation velocity prefactor [m/s] for each slip system + dot_N_0_tw, & !< twin nucleation rate [1/m³s] for each twin system + dot_N_0_tr, & !< trans nucleation rate [1/m³s] for each trans system + t_tw, & !< twin thickness [m] for each twin system + CLambdaSlip, & !< Adj. parameter for distance between 2 forest dislocations for each slip system + atomicVolume, & + t_tr, & !< martensite lamellar thickness [m] for each trans system and instance + p, & !< p-exponent in glide velocity + q, & !< q-exponent in glide velocity + r, & !< r-exponent in twin nucleation rate + s, & !< s-exponent in trans nucleation rate + gamma_char, & !< characteristic shear for twins + B !< drag coefficient + real(pReal), dimension(:,:), allocatable :: & + 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 :: & + n0_sl, & !< slip system normal + forestProjection, & + C66 + real(pReal), dimension(:,:,:), allocatable :: & + P_tr, & + P_sl, & + P_tw, & + C66_tw, & + C66_tr + 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 :: & + 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 + 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, & + rho_dip, & + gamma_sl, & + 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 + Lambda_tw, & !< mean free path between 2 obstacles seen by a growing twin + Lambda_tr, & !< mean free path between 2 obstacles seen by a growing martensite + tau_pass, & + tau_hat_tw, & + tau_hat_tr, & + V_tw, & !< volume of a new twin + V_tr, & !< volume of a new martensite disc + 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 + type(tDislotwinState), allocatable, dimension(:) :: & + dotState, & + state + type(tDislotwinMicrostructure), allocatable, dimension(:) :: dependentState + +contains + + +!-------------------------------------------------------------------------------------------------- +!> @brief module initialization +!> @details reads in material parameters, allocates arrays, and does sanity checks +!-------------------------------------------------------------------------------------------------- +module subroutine plastic_dislotwin_init + + integer :: & + Ninstance, & + p, i, & + NipcMyPhase, outputSize, & + 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' + + 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)), & + dot => dotState(phase_plasticityInstance(p)), & + 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) + prm%sum_N_sl = sum(prm%N_sl) + slipActive: if (prm%sum_N_sl > 0) then + prm%P_sl = lattice_SchmidMatrix_slip(prm%N_sl,config%getString('lattice_structure'),& + config%getFloat('c/a',defaultVal=0.0_pReal)) + prm%h_sl_sl = lattice_interaction_SlipBySlip(prm%N_sl, & + config%getFloats('interaction_slipslip'), & + config%getString('lattice_structure')) + 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)) + 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)) + prm%b_sl = config%getFloats('slipburgers',requiredSize=size(prm%N_sl)) + prm%Delta_F = config%getFloats('qedge', requiredSize=size(prm%N_sl)) + prm%CLambdaSlip = config%getFloats('clambdaslip',requiredSize=size(prm%N_sl)) + prm%p = config%getFloats('p_slip', requiredSize=size(prm%N_sl)) + 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') + prm%Qsd = config%getFloat('qsd') + prm%atomicVolume = config%getFloat('catomicvolume') * prm%b_sl**3.0_pReal + prm%ExtendedDislocations = config%keyExists('/extend_dislocations/') + if (prm%ExtendedDislocations) then + 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) + prm%v0 = math_expand(prm%v0, prm%N_sl) + prm%b_sl = math_expand(prm%b_sl, prm%N_sl) + prm%Delta_F = math_expand(prm%Delta_F, prm%N_sl) + prm%CLambdaSlip = math_expand(prm%CLambdaSlip, prm%N_sl) + 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) + + ! sanity checks + if ( prm%D0 <= 0.0_pReal) extmsg = trim(extmsg)//' D0' + if ( prm%Qsd <= 0.0_pReal) extmsg = trim(extmsg)//' Qsd' + if (any(prm%rho_mob_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_mob_0' + if (any(prm%rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_dip_0' + 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)//' Delta_F' + if (any(prm%CLambdaSlip <= 0.0_pReal)) extmsg = trim(extmsg)//' CLambdaSlip' + 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) + prm%sum_N_tw = sum(prm%N_tw) + if (prm%sum_N_tw > 0) then + prm%P_tw = lattice_SchmidMatrix_twin(prm%N_tw,config%getString('lattice_structure'),& + config%getFloat('c/a',defaultVal=0.0_pReal)) + 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 = 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)) + allocate(prm%b_tw (0)) + allocate(prm%r (0)) + allocate(prm%h_tw_tw (0,0)) + endif + +!-------------------------------------------------------------------------------------------------- +! transformation related parameters + prm%N_tr = config%getInts('ntrans', defaultVal=emptyIntArray) + prm%sum_N_tr = sum(prm%N_tr) + 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) + endif + prm%t_tr = config%getFloats('lamellarsize') + prm%t_tr = math_expand(prm%t_tr,prm%N_tr) + prm%s = config%getFloats('s_trans',defaultVal=[0.0_pReal]) + prm%s = math_expand(prm%s,prm%N_tr) + else + allocate(prm%t_tr (0)) + allocate(prm%b_tr (0)) + 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 + 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 + +!-------------------------------------------------------------------------------------------------- +! shearband related parameters + prm%sbVelocity = config%getFloat('shearbandvelocity',defaultVal=0.0_pReal) + 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//')') + if (prm%aTol_f_tw <= 0.0_pReal) & + call IO_error(211,el=p,ext_msg='aTol_f_tw ('//PLASTICITY_DISLOTWIN_label//')') + endif + if (prm%sum_N_tr > 0) then + 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 + +!-------------------------------------------------------------------------------------------------- +! allocate state arrays + NipcMyPhase = count(material_phaseAt == p) * discretization_nIP + sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl & + + 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 + 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,:) + 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,:) + dot%gamma_sl=>plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = 1.0e6_pReal !ToDo: better make optional parameter + ! global alias + plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%accumulatedSlip => plasticState(p)%state(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 +!-------------------------------------------------------------------------------------------------- +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 & + + stt%f_tw(i,of)*prm%C66_tw(1:6,1:6,i) + enddo + do i=1,prm%sum_N_tr + 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 +!-------------------------------------------------------------------------------------------------- +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,& + BoltzmannRatio, & + ddot_gamma_dtau, & + tau + real(pReal), dimension(param(instance)%sum_N_sl) :: & + dot_gamma_sl,ddot_gamma_dtau_slip + real(pReal), dimension(param(instance)%sum_N_tw) :: & + dot_gamma_twin,ddot_gamma_dtau_twin + real(pReal), dimension(param(instance)%sum_N_tr) :: & + dot_gamma_tr,ddot_gamma_dtau_trans + real(pReal):: dot_gamma_sb + real(pReal), dimension(3,3) :: eigVectors, P_sb + real(pReal), dimension(3) :: eigValues + logical :: error + real(pReal), dimension(3,6), parameter :: & + sb_sComposition = & + reshape(real([& + 1, 0, 1, & + 1, 0,-1, & + 1, 1, 0, & + 1,-1, 0, & + 0, 1, 1, & + 0, 1,-1 & + ],pReal),[ 3,6]), & + sb_mComposition = & + reshape(real([& + 1, 0,-1, & + 1, 0,+1, & + 1,-1, 0, & + 1, 1, 0, & + 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 + + 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) + 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_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 + 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_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 + 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_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 +!-------------------------------------------------------------------------------------------------- +module subroutine plastic_dislotwin_dotState(Mp,T,instance,of) + + real(pReal), dimension(3,3), intent(in):: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T !< temperature at integration point + integer, intent(in) :: & + instance, & + of + + integer :: i + real(pReal) :: & + f_unrotated, & + VacancyDiffusion, & + rho_dip_distance, & + v_cl, & !< climb velocity + Gamma, & !< stacking fault energy + tau, & + 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, & + dot_rho_dip_climb, & + rho_dip_distance_min, & + dot_gamma_sl + real(pReal), dimension(param(instance)%sum_N_tw) :: & + dot_gamma_twin + real(pReal), dimension(param(instance)%sum_N_tr) :: & + dot_gamma_tr + + associate(prm => param(instance), stt => state(instance), & + dot => dotState(instance), dst => dependentState(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)) + VacancyDiffusion = prm%D0*exp(-prm%Qsd/(kB*T)) + + 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)) + + significantSlipStress: if (dEq0(tau)) then + dot_rho_dip_formation(i) = 0.0_pReal + dot_rho_dip_climb(i) = 0.0_pReal + else significantSlipStress + rho_dip_distance = 3.0_pReal*prm%mu*prm%b_sl(i)/(16.0_pReal*PI*abs(tau)) + rho_dip_distance = math_clip(rho_dip_distance, right = dst%Lambda_sl(i,of)) + rho_dip_distance = math_clip(rho_dip_distance, left = rho_dip_distance_min(i)) + + if (prm%dipoleFormation) then + dot_rho_dip_formation(i) = 2.0_pReal*(rho_dip_distance-rho_dip_distance_min(i))/prm%b_sl(i) & + * stt%rho_mob(i,of)*abs(dot_gamma_sl(i)) + else + dot_rho_dip_formation(i) = 0.0_pReal + endif + + if (dEq(rho_dip_distance,rho_dip_distance_min(i))) then + dot_rho_dip_climb(i) = 0.0_pReal + 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 + 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 + 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 + endif significantSlipStress + enddo slipState + + dot%rho_mob(:,of) = abs(dot_gamma_sl)/(prm%b_sl*dst%Lambda_sl(:,of)) & + - dot_rho_dip_formation & + - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_mob(:,of)*abs(dot_gamma_sl) + + dot%rho_dip(:,of) = dot_rho_dip_formation & + - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,of)*abs(dot_gamma_sl) & + - dot_rho_dip_climb + + call kinetics_twin(Mp,T,dot_gamma_sl,instance,of,dot_gamma_twin) + dot%f_tw(:,of) = f_unrotated*dot_gamma_twin/prm%gamma_char + + call kinetics_trans(Mp,T,dot_gamma_sl,instance,of,dot_gamma_tr) + dot%f_tr(:,of) = f_unrotated*dot_gamma_tr + + end associate + +end subroutine plastic_dislotwin_dotState + + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates derived quantities from state +!-------------------------------------------------------------------------------------------------- +module subroutine plastic_dislotwin_dependentState(T,instance,of) + + integer, intent(in) :: & + instance, & + of + real(pReal), intent(in) :: & + T + + real(pReal) :: & + sumf_twin,Gamma,sumf_trans + real(pReal), dimension(param(instance)%sum_N_sl) :: & + inv_lambda_sl_sl, & !< 1/mean free distance between 2 forest dislocations seen by a moving dislocation + inv_lambda_sl_tw, & !< 1/mean free distance between 2 twin stacks from different systems seen by a moving dislocation + inv_lambda_sl_tr !< 1/mean free distance between 2 martensite lamellar from different systems seen by a moving dislocation + real(pReal), dimension(param(instance)%sum_N_tw) :: & + inv_lambda_tw_tw, & !< 1/mean free distance between 2 twin stacks from different systems seen by a growing twin + f_over_t_tw + real(pReal), dimension(param(instance)%sum_N_tr) :: & + inv_lambda_tr_tr, & !< 1/mean free distance between 2 martensite stacks from different systems seen by a growing martensite + 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)) + else + 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) & + + 3.0_pReal*prm%b_tw*prm%mu/(prm%L_tw*prm%b_sl) ! slip burgers here correct? + 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) + + 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 +!-------------------------------------------------------------------------------------------------- +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)) + + 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 (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_tr_ID) + call results_writeDataset(group,stt%f_tr,'f_tr',& + 'martensite volume fraction','m³/m³') + + end select + enddo outputsLoop + end associate + +end subroutine plastic_dislotwin_results + + +!-------------------------------------------------------------------------------------------------- +!> @brief Shear rates on slip systems, their derivatives with respect to resolved stress and the +! resolved stresss +!> @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 +!-------------------------------------------------------------------------------------------------- +pure subroutine kinetics_slip(Mp,T,instance,of, & + dot_gamma_sl,ddot_gamma_dtau_slip,tau_slip) + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T !< temperature + 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) :: & + ddot_gamma_dtau_slip, & + tau_slip + real(pReal), dimension(param(instance)%sum_N_sl) :: & + ddot_gamma_dtau + + real(pReal), dimension(param(instance)%sum_N_sl) :: & + tau, & + stressRatio, & + StressRatio_p, & + BoltzmannRatio, & + v_wait_inverse, & !< inverse of the effective velocity of a dislocation waiting at obstacles (unsigned) + v_run_inverse, & !< inverse of the velocity of a free moving dislocation (unsigned) + dV_wait_inverse_dTau, & + dV_run_inverse_dTau, & + dV_dTau, & + tau_eff !< effective resolved stress + 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) & + / prm%tau_0 + dV_run_inverse_dTau = -1.0_pReal * v_run_inverse/tau_eff + dV_dTau = -1.0_pReal * (dV_wait_inverse_dTau+dV_run_inverse_dTau) & + / (v_wait_inverse+v_run_inverse)**2.0_pReal + ddot_gamma_dtau = dV_dTau*stt%rho_mob(:,of)*prm%b_sl + else where significantStress + 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 +!-------------------------------------------------------------------------------------------------- +pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,of,& + dot_gamma_twin,ddot_gamma_dtau_twin) + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T !< temperature + integer, intent(in) :: & + instance, & + 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 + s1=prm%fcc_twinNucleationSlipPair(1,i) + s2=prm%fcc_twinNucleationSlipPair(2,i) + if (tau(i) < dst%tau_r_tw(i,of)) then ! ToDo: correct? + Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,of)+stt%rho_dip(s2,of))+& + abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,of)+stt%rho_dip(s1,of)))/& ! ToDo: MD: it would be more consistent to use shearrates from state + (prm%L_tw*prm%b_sl(i))*& + (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tw(i,of)-tau(i)))) ! P_ncs + else + Ndot0=0.0_pReal + end if + else isFCC + 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) + ddot_gamma_dtau = (dot_gamma_twin*prm%r/tau)*StressRatio_r + else where significantStress + 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 +!-------------------------------------------------------------------------------------------------- +pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,& + dot_gamma_tr,ddot_gamma_dtau_trans) + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T !< temperature + integer, intent(in) :: & + instance, & + 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 + s1=prm%fcc_twinNucleationSlipPair(1,i) + s2=prm%fcc_twinNucleationSlipPair(2,i) + if (tau(i) < dst%tau_r_tr(i,of)) then ! ToDo: correct? + Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,of)+stt%rho_dip(s2,of))+& + abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,of)+stt%rho_dip(s1,of)))/& ! ToDo: MD: it would be more consistent to use shearrates from state + (prm%L_tr*prm%b_sl(i))*& + (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tr(i,of)-tau(i)))) ! P_ncs + else + Ndot0=0.0_pReal + end if + else isFCC + 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) + ddot_gamma_dtau = (dot_gamma_tr*prm%s/tau)*StressRatio_s + else where significantStress + 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 + +end submodule plastic_dislotwin diff --git a/src/plastic_isotropic.f90 b/src/constitutive_plastic_isotropic.f90 similarity index 95% rename from src/plastic_isotropic.f90 rename to src/constitutive_plastic_isotropic.f90 index 96d70be4a..a7a9acfee 100644 --- a/src/plastic_isotropic.f90 +++ b/src/constitutive_plastic_isotropic.f90 @@ -7,18 +7,7 @@ !! resolving the stress on the slip systems. Will give the response of phenopowerlaw for an !! untextured polycrystal !-------------------------------------------------------------------------------------------------- -module plastic_isotropic - use prec - use debug - use math - use IO - use material - use config - use discretization - use results - - implicit none - private +submodule(constitutive) plastic_isotropic enum, bind(c) enumerator :: & @@ -64,20 +53,13 @@ module plastic_isotropic dotState, & state - public :: & - plastic_isotropic_init, & - plastic_isotropic_LpAndItsTangent, & - plastic_isotropic_LiAndItsTangent, & - plastic_isotropic_dotState, & - plastic_isotropic_results - contains !-------------------------------------------------------------------------------------------------- !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine plastic_isotropic_init +module subroutine plastic_isotropic_init integer :: & Ninstance, & @@ -207,7 +189,7 @@ end subroutine plastic_isotropic_init !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- -subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) +module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) real(pReal), dimension(3,3), intent(out) :: & Lp !< plastic velocity gradient @@ -268,7 +250,7 @@ end subroutine plastic_isotropic_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- -subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) +module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) real(pReal), dimension(3,3), intent(out) :: & Li !< inleastic velocity gradient @@ -320,7 +302,7 @@ subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) !-------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine plastic_isotropic_dotState(Mp,instance,of) +module subroutine plastic_isotropic_dotState(Mp,instance,of) real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress @@ -369,7 +351,7 @@ end subroutine plastic_isotropic_dotState !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- -subroutine plastic_isotropic_results(instance,group) +module subroutine plastic_isotropic_results(instance,group) integer, intent(in) :: instance character(len=*), intent(in) :: group @@ -388,4 +370,4 @@ subroutine plastic_isotropic_results(instance,group) end subroutine plastic_isotropic_results -end module plastic_isotropic +end submodule plastic_isotropic diff --git a/src/constitutive_plastic_kinehardening.f90 b/src/constitutive_plastic_kinehardening.f90 new file mode 100644 index 000000000..4949c6533 --- /dev/null +++ b/src/constitutive_plastic_kinehardening.f90 @@ -0,0 +1,522 @@ +!-------------------------------------------------------------------------------------------------- +!> @author Philip Eisenlohr, Michigan State University +!> @author Zhuowen Zhao, Michigan State University +!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH +!> @brief Phenomenological crystal plasticity using a power law formulation for the shear rates +!! and a Voce-type kinematic hardening rule +!-------------------------------------------------------------------------------------------------- +submodule(constitutive) plastic_kinehardening + + enum, bind(c) + enumerator :: & + undefined_ID, & + crss_ID, & !< critical resolved stress + crss_back_ID, & !< critical resolved back stress + sense_ID, & !< sense of acting shear stress (-1 or +1) + chi0_ID, & !< backstress at last switch of stress sense (positive?) + gamma0_ID, & !< accumulated shear at last switch of stress sense (at current switch?) + accshear_ID, & + shearrate_ID, & + resolvedstress_ID + end enum + + type :: tParameters + real(pReal) :: & + gdot0, & !< reference shear strain rate for slip + n, & !< stress exponent for slip + aTolResistance, & + aTolShear + real(pReal), allocatable, dimension(:) :: & + crss0, & !< initial critical shear stress for slip + theta0, & !< initial hardening rate of forward stress for each slip + theta1, & !< asymptotic hardening rate of forward stress for each slip + theta0_b, & !< initial hardening rate of back stress for each slip + theta1_b, & !< asymptotic hardening rate of back stress for each slip + tau1, & + tau1_b, & + nonSchmidCoeff + real(pReal), allocatable, dimension(:,:) :: & + interaction_slipslip !< slip resistance from slip activity + real(pReal), allocatable, dimension(:,:,:) :: & + Schmid, & + nonSchmid_pos, & + nonSchmid_neg + integer :: & + totalNslip, & !< total number of active slip system + of_debug = 0 + integer, allocatable, dimension(:) :: & + Nslip !< number of active slip systems for each family + integer(kind(undefined_ID)), allocatable, dimension(:) :: & + outputID !< ID of each post result output + end type tParameters + + type :: tKinehardeningState + real(pReal), pointer, dimension(:,:) :: & !< vectors along NipcMyInstance + crss, & !< critical resolved stress + crss_back, & !< critical resolved back stress + sense, & !< sense of acting shear stress (-1 or +1) + chi0, & !< backstress at last switch of stress sense + gamma0, & !< accumulated shear at last switch of stress sense + accshear !< accumulated (absolute) shear + end type tKinehardeningState + +!-------------------------------------------------------------------------------------------------- +! containers for parameters and state + type(tParameters), allocatable, dimension(:) :: param + type(tKinehardeningState), allocatable, dimension(:) :: & + dotState, & + deltaState, & + state + +contains + + +!-------------------------------------------------------------------------------------------------- +!> @brief module initialization +!> @details reads in material parameters, allocates arrays, and does sanity checks +!-------------------------------------------------------------------------------------------------- +module subroutine plastic_kinehardening_init + + integer :: & + Ninstance, & + p, i, o, & + NipcMyPhase, & + sizeState, sizeDeltaState, sizeDotState, & + startIndex, endIndex + + integer(kind(undefined_ID)) :: & + outputID + + character(len=pStringLen) :: & + extmsg = '' + character(len=pStringLen), dimension(:), allocatable :: & + outputs + + write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_KINEHARDENING_label//' init -+>>>' + + Ninstance = count(phase_plasticity == PLASTICITY_KINEHARDENING_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(deltaState(Ninstance)) + + do p = 1, size(phase_plasticityInstance) + if (phase_plasticity(p) /= PLASTICITY_KINEHARDENING_ID) cycle + associate(prm => param(phase_plasticityInstance(p)), & + dot => dotState(phase_plasticityInstance(p)), & + dlt => deltaState(phase_plasticityInstance(p)), & + stt => state(phase_plasticityInstance(p)),& + config => config_phase(p)) + +#ifdef DEBUG + if (p==material_phaseAt(debug_g,debug_e)) then + prm%of_debug = material_phasememberAt(debug_g,debug_i,debug_e) + endif +#endif + +!-------------------------------------------------------------------------------------------------- +! optional parameters that need to be defined + prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal) + prm%aTolShear = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal) + + ! sanity checks + if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//' aTolresistance' + if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//' aTolShear' + +!-------------------------------------------------------------------------------------------------- +! slip related parameters + prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) + prm%totalNslip = sum(prm%Nslip) + slipActive: if (prm%totalNslip > 0) then + prm%Schmid = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),& + 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(prm%Nslip,prm%nonSchmidCoeff,+1) + prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1) + else + prm%nonSchmid_pos = prm%Schmid + prm%nonSchmid_neg = prm%Schmid + endif + prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, & + config%getFloats('interaction_slipslip'), & + config%getString('lattice_structure')) + + prm%crss0 = config%getFloats('crss0', requiredSize=size(prm%Nslip)) + prm%tau1 = config%getFloats('tau1', requiredSize=size(prm%Nslip)) + prm%tau1_b = config%getFloats('tau1_b', requiredSize=size(prm%Nslip)) + prm%theta0 = config%getFloats('theta0', requiredSize=size(prm%Nslip)) + prm%theta1 = config%getFloats('theta1', requiredSize=size(prm%Nslip)) + prm%theta0_b = config%getFloats('theta0_b', requiredSize=size(prm%Nslip)) + prm%theta1_b = config%getFloats('theta1_b', requiredSize=size(prm%Nslip)) + + prm%gdot0 = config%getFloat('gdot0') + prm%n = config%getFloat('n_slip') + + ! expand: family => system + prm%crss0 = math_expand(prm%crss0, prm%Nslip) + prm%tau1 = math_expand(prm%tau1, prm%Nslip) + prm%tau1_b = math_expand(prm%tau1_b, prm%Nslip) + prm%theta0 = math_expand(prm%theta0, prm%Nslip) + prm%theta1 = math_expand(prm%theta1, prm%Nslip) + prm%theta0_b = math_expand(prm%theta0_b,prm%Nslip) + prm%theta1_b = math_expand(prm%theta1_b,prm%Nslip) + + + +!-------------------------------------------------------------------------------------------------- +! sanity checks + if ( prm%gdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0' + if ( prm%n <= 0.0_pReal) extmsg = trim(extmsg)//' n_slip' + if (any(prm%crss0 <= 0.0_pReal)) extmsg = trim(extmsg)//' crss0' + if (any(prm%tau1 <= 0.0_pReal)) extmsg = trim(extmsg)//' tau1' + if (any(prm%tau1_b <= 0.0_pReal)) extmsg = trim(extmsg)//' tau1_b' + + !ToDo: Any sensible checks for theta? + + endif slipActive + +!-------------------------------------------------------------------------------------------------- +! exit if any parameter is out of range + if (extmsg /= '') & + call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_KINEHARDENING_label//')') + +!-------------------------------------------------------------------------------------------------- +! output pararameters + outputs = config%getStrings('(output)',defaultVal=emptyStringArray) + allocate(prm%outputID(0)) + do i=1, size(outputs) + outputID = undefined_ID + select case(outputs(i)) + + case ('resistance') + outputID = merge(crss_ID,undefined_ID,prm%totalNslip>0) + case ('accumulatedshear') + outputID = merge(accshear_ID,undefined_ID,prm%totalNslip>0) + case ('shearrate') + outputID = merge(shearrate_ID,undefined_ID,prm%totalNslip>0) + case ('resolvedstress') + outputID = merge(resolvedstress_ID,undefined_ID,prm%totalNslip>0) + case ('backstress') + outputID = merge(crss_back_ID,undefined_ID,prm%totalNslip>0) + case ('sense') + outputID = merge(sense_ID,undefined_ID,prm%totalNslip>0) + case ('chi0') + outputID = merge(chi0_ID,undefined_ID,prm%totalNslip>0) + case ('gamma0') + outputID = merge(gamma0_ID,undefined_ID,prm%totalNslip>0) + + end select + + if (outputID /= undefined_ID) then + prm%outputID = [prm%outputID , outputID] + endif + + enddo + +!-------------------------------------------------------------------------------------------------- +! allocate state arrays + NipcMyPhase = count(material_phaseAt == p) * discretization_nIP + sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%totalNslip + sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%totalNslip + sizeState = sizeDotState + sizeDeltaState + + call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) + +!-------------------------------------------------------------------------------------------------- +! locally defined state aliases and initialization of state0 and aTolState + startIndex = 1 + endIndex = prm%totalNslip + stt%crss => plasticState(p)%state (startIndex:endIndex,:) + stt%crss = spread(prm%crss0, 2, NipcMyPhase) + dot%crss => plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance + + startIndex = endIndex + 1 + endIndex = endIndex + prm%totalNslip + stt%crss_back => plasticState(p)%state (startIndex:endIndex,:) + dot%crss_back => plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance + + startIndex = endIndex + 1 + endIndex = endIndex + prm%totalNslip + stt%accshear => plasticState(p)%state (startIndex:endIndex,:) + dot%accshear => plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear + ! global alias + plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%accumulatedSlip => plasticState(p)%state(startIndex:endIndex,:) + + o = plasticState(p)%offsetDeltaState + startIndex = endIndex + 1 + endIndex = endIndex + prm%totalNslip + stt%sense => plasticState(p)%state (startIndex :endIndex ,:) + dlt%sense => plasticState(p)%deltaState(startIndex-o:endIndex-o,:) + + startIndex = endIndex + 1 + endIndex = endIndex + prm%totalNslip + stt%chi0 => plasticState(p)%state (startIndex :endIndex ,:) + dlt%chi0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:) + + startIndex = endIndex + 1 + endIndex = endIndex + prm%totalNslip + stt%gamma0 => plasticState(p)%state (startIndex :endIndex ,:) + dlt%gamma0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:) + + plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally + + end associate + + enddo + +end subroutine plastic_kinehardening_init + + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates plastic velocity gradient and its tangent +!-------------------------------------------------------------------------------------------------- +pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) + + real(pReal), dimension(3,3), intent(out) :: & + Lp !< plastic velocity gradient + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp !< derivative of Lp with respect to the Mandel stress + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + + integer :: & + i,k,l,m,n + real(pReal), dimension(param(instance)%totalNslip) :: & + gdot_pos,gdot_neg, & + dgdot_dtau_pos,dgdot_dtau_neg + + Lp = 0.0_pReal + dLp_dMp = 0.0_pReal + + associate(prm => param(instance)) + + call kinetics(Mp,instance,of,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) + + do i = 1, prm%totalNslip + Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%Schmid(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_dtau_pos(i) * prm%Schmid(k,l,i) * prm%nonSchmid_pos(m,n,i) & + + dgdot_dtau_neg(i) * prm%Schmid(k,l,i) * prm%nonSchmid_neg(m,n,i) + enddo + + end associate + +end subroutine plastic_kinehardening_LpAndItsTangent + + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates the rate of change of microstructure +!-------------------------------------------------------------------------------------------------- +module subroutine plastic_kinehardening_dotState(Mp,instance,of) + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + + real(pReal) :: & + sumGamma + real(pReal), dimension(param(instance)%totalNslip) :: & + gdot_pos,gdot_neg + + + associate(prm => param(instance), stt => state(instance), dot => dotState(instance)) + + call kinetics(Mp,instance,of,gdot_pos,gdot_neg) + dot%accshear(:,of) = abs(gdot_pos+gdot_neg) + sumGamma = sum(stt%accshear(:,of)) + + + dot%crss(:,of) = matmul(prm%interaction_SlipSlip,dot%accshear(:,of)) & + * ( prm%theta1 & + + (prm%theta0 - prm%theta1 + prm%theta0*prm%theta1*sumGamma/prm%tau1) & + * exp(-sumGamma*prm%theta0/prm%tau1) & + ) + + dot%crss_back(:,of) = stt%sense(:,of)*dot%accshear(:,of) * & + ( prm%theta1_b + & + (prm%theta0_b - prm%theta1_b & + + prm%theta0_b*prm%theta1_b/(prm%tau1_b+stt%chi0(:,of))*(stt%accshear(:,of)-stt%gamma0(:,of))& + ) *exp(-(stt%accshear(:,of)-stt%gamma0(:,of)) *prm%theta0_b/(prm%tau1_b+stt%chi0(:,of))) & + ) + + end associate + +end subroutine plastic_kinehardening_dotState + + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates (instantaneous) incremental change of microstructure +!-------------------------------------------------------------------------------------------------- +module subroutine plastic_kinehardening_deltaState(Mp,instance,of) + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + + real(pReal), dimension(param(instance)%totalNslip) :: & + gdot_pos,gdot_neg, & + sense + + associate(prm => param(instance), stt => state(instance), dlt => deltaState(instance)) + + call kinetics(Mp,instance,of,gdot_pos,gdot_neg) + sense = merge(state(instance)%sense(:,of), & ! keep existing... + sign(1.0_pReal,gdot_pos+gdot_neg), & ! ...or have a defined + dEq0(gdot_pos+gdot_neg,1e-10_pReal)) ! current sense of shear direction + +#ifdef DEBUG + if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0 & + .and. (of == prm%of_debug & + .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0)) then + write(6,'(a)') '======= kinehardening delta state =======' + write(6,*) sense,state(instance)%sense(:,of) + endif +#endif + +!-------------------------------------------------------------------------------------------------- +! switch in sense of shear? + where(dNeq(sense,stt%sense(:,of),0.1_pReal)) + dlt%sense (:,of) = sense - stt%sense(:,of) ! switch sense + dlt%chi0 (:,of) = abs(stt%crss_back(:,of)) - stt%chi0(:,of) ! remember current backstress magnitude + dlt%gamma0(:,of) = stt%accshear(:,of) - stt%gamma0(:,of) ! remember current accumulated shear + else where + dlt%sense (:,of) = 0.0_pReal + dlt%chi0 (:,of) = 0.0_pReal + dlt%gamma0(:,of) = 0.0_pReal + end where + + end associate + +end subroutine plastic_kinehardening_deltaState + + +!-------------------------------------------------------------------------------------------------- +!> @brief writes results to HDF5 output file +!-------------------------------------------------------------------------------------------------- +module subroutine plastic_kinehardening_results(instance,group) + + integer, intent(in) :: instance + character(len=*), intent(in) :: group + integer :: o + + associate(prm => param(instance), stt => state(instance)) + outputsLoop: do o = 1,size(prm%outputID) + select case(prm%outputID(o)) + case (crss_ID) + call results_writeDataset(group,stt%crss,'xi_sl', & + 'resistance against plastic slip','Pa') + + case(crss_back_ID) + call results_writeDataset(group,stt%crss_back,'tau_back', & + 'back stress against plastic slip','Pa') + + case (sense_ID) + call results_writeDataset(group,stt%sense,'sense_of_shear','tbd','1') + + case (chi0_ID) + call results_writeDataset(group,stt%chi0,'chi0','tbd','Pa') + + case (gamma0_ID) + call results_writeDataset(group,stt%gamma0,'gamma0','tbd','1') + + case (accshear_ID) + call results_writeDataset(group,stt%accshear,'gamma_sl', & + 'plastic shear','1') + + end select + enddo outputsLoop + end associate + +end subroutine plastic_kinehardening_results + + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates shear rates on slip systems and derivatives with respect to resolved stress +!> @details: Shear rates 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(Mp,instance,of, & + gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + + real(pReal), intent(out), dimension(param(instance)%totalNslip) :: & + gdot_pos, & + gdot_neg + real(pReal), intent(out), optional, dimension(param(instance)%totalNslip) :: & + dgdot_dtau_pos, & + dgdot_dtau_neg + + real(pReal), dimension(param(instance)%totalNslip) :: & + tau_pos, & + tau_neg + integer :: i + logical :: nonSchmidActive + + associate(prm => param(instance), stt => state(instance)) + + nonSchmidActive = size(prm%nonSchmidCoeff) > 0 + + do i = 1, prm%totalNslip + tau_pos(i) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - stt%crss_back(i,of) + tau_neg(i) = merge(math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)) - stt%crss_back(i,of), & + 0.0_pReal, nonSchmidActive) + enddo + + where(dNeq0(tau_pos)) + gdot_pos = prm%gdot0 * merge(0.5_pReal,1.0_pReal, nonSchmidActive) & ! 1/2 if non-Schmid active + * sign(abs(tau_pos/stt%crss(:,of))**prm%n, tau_pos) + else where + gdot_pos = 0.0_pReal + end where + + where(dNeq0(tau_neg)) + gdot_neg = prm%gdot0 * 0.5_pReal & ! only used if non-Schmid active, always 1/2 + * sign(abs(tau_neg/stt%crss(:,of))**prm%n, tau_neg) + else where + gdot_neg = 0.0_pReal + end where + + if (present(dgdot_dtau_pos)) then + where(dNeq0(gdot_pos)) + dgdot_dtau_pos = gdot_pos*prm%n/tau_pos + else where + dgdot_dtau_pos = 0.0_pReal + end where + endif + if (present(dgdot_dtau_neg)) then + where(dNeq0(gdot_neg)) + dgdot_dtau_neg = gdot_neg*prm%n/tau_neg + else where + dgdot_dtau_neg = 0.0_pReal + end where + endif + end associate + +end subroutine kinetics + +end submodule plastic_kinehardening diff --git a/src/plastic_none.f90 b/src/constitutive_plastic_none.f90 similarity index 53% rename from src/plastic_none.f90 rename to src/constitutive_plastic_none.f90 index f77b19c09..578fb4149 100644 --- a/src/plastic_none.f90 +++ b/src/constitutive_plastic_none.f90 @@ -4,16 +4,7 @@ !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief Dummy plasticity for purely elastic material !-------------------------------------------------------------------------------------------------- -module plastic_none - use material - use discretization - use debug - - implicit none - private - - public :: & - plastic_none_init +submodule(constitutive) plastic_none contains @@ -21,27 +12,27 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine plastic_none_init +module subroutine plastic_none_init - integer :: & - Ninstance, & - p, & - NipcMyPhase - - write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_NONE_label//' init -+>>>' - - Ninstance = count(phase_plasticity == PLASTICITY_NONE_ID) - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & - write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - - do p = 1, size(phase_plasticity) - if (phase_plasticity(p) /= PLASTICITY_NONE_ID) cycle - - NipcMyPhase = count(material_phaseAt == p) * discretization_nIP - call material_allocatePlasticState(p,NipcMyPhase,0,0,0) - - enddo + integer :: & + Ninstance, & + p, & + NipcMyPhase + + write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_NONE_label//' init -+>>>' + + Ninstance = count(phase_plasticity == PLASTICITY_NONE_ID) + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0) & + write(6,'(a16,1x,i5,/)') '# instances:',Ninstance + + do p = 1, size(phase_plasticity) + if (phase_plasticity(p) /= PLASTICITY_NONE_ID) cycle + + NipcMyPhase = count(material_phaseAt == p) * discretization_nIP + call material_allocatePlasticState(p,NipcMyPhase,0,0,0) + + enddo end subroutine plastic_none_init -end module plastic_none +end submodule plastic_none diff --git a/src/plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 similarity index 98% rename from src/plastic_nonlocal.f90 rename to src/constitutive_plastic_nonlocal.f90 index 3977f1fa7..42b8a34f5 100644 --- a/src/plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -4,18 +4,7 @@ !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @brief material subroutine for plasticity including dislocation flux !-------------------------------------------------------------------------------------------------- -module plastic_nonlocal - use prec - use IO - use math - use debug - use material - use lattice - use rotations - use results - use config - use lattice - use discretization +submodule(constitutive) plastic_nonlocal use geometry_plastic_nonlocal, only: & nIPneighbors => geometry_plastic_nonlocal_nIPneighbors, & IPneighborhood => geometry_plastic_nonlocal_IPneighborhood, & @@ -23,11 +12,9 @@ module plastic_nonlocal IParea => geometry_plastic_nonlocal_IParea0, & IPareaNormal => geometry_plastic_nonlocal_IPareaNormal0 - implicit none - private real(pReal), parameter :: & KB = 1.38e-23_pReal !< Physical parameter, Boltzmann constant in J/Kelvin - + ! storage order of dislocation types integer, dimension(8), parameter :: & sgl = [1,2,3,4,5,6,7,8] !< signed (single) @@ -199,25 +186,13 @@ module plastic_nonlocal type(tNonlocalMicrostructure), dimension(:), allocatable :: microstructure - public :: & - plastic_nonlocal_init, & - plastic_nonlocal_dependentState, & - plastic_nonlocal_LpAndItsTangent, & - plastic_nonlocal_dotState, & - plastic_nonlocal_deltaState, & - plastic_nonlocal_updateCompatibility, & - plastic_nonlocal_results - - private :: & - plastic_nonlocal_kinetics - contains !-------------------------------------------------------------------------------------------------- !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine plastic_nonlocal_init +module subroutine plastic_nonlocal_init integer :: & sizeState, sizeDotState,sizeDependentState, sizeDeltaState, & @@ -254,7 +229,6 @@ subroutine plastic_nonlocal_init allocate(dotState(maxNinstances)) allocate(deltaState(maxNinstances)) allocate(microstructure(maxNinstances)) - allocate(totalNslip(maxNinstances), source=0) @@ -738,7 +712,7 @@ end subroutine plastic_nonlocal_init !-------------------------------------------------------------------------------------------------- !> @brief calculates quantities characterizing the microstructure !-------------------------------------------------------------------------------------------------- -subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) +module subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el) integer, intent(in) :: & ip, & @@ -1100,8 +1074,8 @@ end subroutine plastic_nonlocal_kinetics !-------------------------------------------------------------------------------------------------- !> @brief calculates plastic velocity gradient and its tangent !-------------------------------------------------------------------------------------------------- -subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, & - Mp, Temperature, volume, ip, el) +module subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, & + Mp, Temperature, volume, ip, el) integer, intent(in) :: & ip, & !< current integration point @@ -1230,7 +1204,7 @@ end subroutine plastic_nonlocal_LpAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief (instantaneous) incremental change of microstructure !-------------------------------------------------------------------------------------------------- -subroutine plastic_nonlocal_deltaState(Mp,ip,el) +module subroutine plastic_nonlocal_deltaState(Mp,ip,el) integer, intent(in) :: & ip, & @@ -1346,8 +1320,8 @@ end subroutine plastic_nonlocal_deltaState !--------------------------------------------------------------------------------------------------- !> @brief calculates the rate of change of microstructure !--------------------------------------------------------------------------------------------------- -subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & - timestep,ip,el) +module subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, & + timestep,ip,el) integer, intent(in) :: & ip, & !< current integration point @@ -1795,13 +1769,13 @@ end subroutine plastic_nonlocal_dotState ! plane normals and signed cosine of the angle between the slip directions. Only the largest values ! that sum up to a total of 1 are considered, all others are set to zero. !-------------------------------------------------------------------------------------------------- -subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) +module subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) integer, intent(in) :: & i, & e type(rotation), dimension(1,discretization_nIP,discretization_nElem), intent(in) :: & - orientation ! crystal orientation in quaternions + orientation ! crystal orientation integer :: & Nneighbors, & ! number of neighbors @@ -1960,10 +1934,10 @@ end function getRho !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- -subroutine plastic_nonlocal_results(instance,group) +module subroutine plastic_nonlocal_results(instance,group) - integer, intent(in) :: instance - character(len=*) :: group + integer, intent(in) :: instance + character(len=*),intent(in) :: group integer :: o associate(prm => param(instance),dst => microstructure(instance),stt=>state(instance)) @@ -2026,4 +2000,4 @@ subroutine plastic_nonlocal_results(instance,group) end subroutine plastic_nonlocal_results -end module plastic_nonlocal +end submodule plastic_nonlocal diff --git a/src/constitutive_plastic_phenopowerlaw.f90 b/src/constitutive_plastic_phenopowerlaw.f90 new file mode 100644 index 000000000..95fc2b46b --- /dev/null +++ b/src/constitutive_plastic_phenopowerlaw.f90 @@ -0,0 +1,596 @@ +!-------------------------------------------------------------------------------------------------- +!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH +!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH +!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH +!> @brief phenomenological crystal plasticity formulation using a powerlaw fitting +!-------------------------------------------------------------------------------------------------- +submodule(constitutive) plastic_phenopowerlaw + + enum, bind(c) + enumerator :: & + undefined_ID, & + resistance_slip_ID, & + accumulatedshear_slip_ID, & + shearrate_slip_ID, & + resolvedstress_slip_ID, & + resistance_twin_ID, & + accumulatedshear_twin_ID, & + shearrate_twin_ID, & + resolvedstress_twin_ID + end enum + + 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, & + aTolResistance, & !< absolute tolerance for integration of xi + aTolShear, & !< absolute tolerance for integration of gamma + aTolTwinfrac !< absolute tolerance for integration of f + 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(:,:) :: & + interaction_SlipSlip, & !< slip resistance from slip activity + interaction_SlipTwin, & !< slip resistance from twin activity + interaction_TwinSlip, & !< twin resistance from slip activity + interaction_TwinTwin !< twin resistance from twin activity + real(pReal), allocatable, dimension(:,:,:) :: & + Schmid_slip, & + Schmid_twin, & + nonSchmid_pos, & + nonSchmid_neg + integer :: & + totalNslip, & !< total number of active slip system + totalNtwin !< total number of active twin systems + integer, allocatable, dimension(:) :: & + Nslip, & !< number of active slip systems for each family + Ntwin !< number of active twin systems for each family + integer(kind(undefined_ID)), allocatable, dimension(:) :: & + outputID !< ID of each post result output + end type tParameters + + type :: tPhenopowerlawState + real(pReal), pointer, dimension(:,:) :: & + xi_slip, & + xi_twin, & + gamma_slip, & + gamma_twin + end type tPhenopowerlawState + +!-------------------------------------------------------------------------------------------------- +! containers for parameters and state + type(tParameters), allocatable, dimension(:) :: param + type(tPhenopowerlawState), allocatable, dimension(:) :: & + dotState, & + state + +contains + + +!-------------------------------------------------------------------------------------------------- +!> @brief module initialization +!> @details reads in material parameters, allocates arrays, and does sanity checks +!-------------------------------------------------------------------------------------------------- +module subroutine plastic_phenopowerlaw_init + + integer :: & + Ninstance, & + p, i, & + NipcMyPhase, outputSize, & + sizeState, sizeDotState, & + startIndex, endIndex + + integer(kind(undefined_ID)) :: & + outputID + + character(len=pStringLen) :: & + extmsg = '' + character(len=pStringLen), dimension(:), allocatable :: & + outputs + + write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_PHENOPOWERLAW_label//' init -+>>>' + + Ninstance = count(phase_plasticity == PLASTICITY_PHENOPOWERLAW_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)) + + do p = 1, size(phase_plasticity) + if (phase_plasticity(p) /= PLASTICITY_PHENOPOWERLAW_ID) cycle + associate(prm => param(phase_plasticityInstance(p)), & + dot => dotState(phase_plasticityInstance(p)), & + 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) + + prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal) + prm%aTolShear = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal) + prm%aTolTwinfrac = config%getFloat('atol_twinfrac', defaultVal=1.0e-6_pReal) + + ! sanity checks + if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//' aTolresistance' + if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//' aTolShear' + if (prm%aTolTwinfrac <= 0.0_pReal) extmsg = trim(extmsg)//' atoltwinfrac' + +!-------------------------------------------------------------------------------------------------- +! slip related parameters + prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) + prm%totalNslip = sum(prm%Nslip) + slipActive: if (prm%totalNslip > 0) then + prm%Schmid_slip = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),& + 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(prm%Nslip,prm%nonSchmidCoeff,+1) + prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1) + else + allocate(prm%nonSchmidCoeff(0)) + prm%nonSchmid_pos = prm%Schmid_slip + prm%nonSchmid_neg = prm%Schmid_slip + endif + prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, & + config%getFloats('interaction_slipslip'), & + config%getString('lattice_structure')) + + prm%xi_slip_0 = config%getFloats('tau0_slip', requiredSize=size(prm%Nslip)) + prm%xi_slip_sat = config%getFloats('tausat_slip', requiredSize=size(prm%Nslip)) + prm%H_int = config%getFloats('h_int', requiredSize=size(prm%Nslip), & + defaultVal=[(0.0_pReal,i=1,size(prm%Nslip))]) + + prm%gdot0_slip = config%getFloat('gdot0_slip') + prm%n_slip = config%getFloat('n_slip') + prm%a_slip = config%getFloat('a_slip') + prm%h0_SlipSlip = config%getFloat('h0_slipslip') + + ! expand: family => system + prm%xi_slip_0 = math_expand(prm%xi_slip_0, prm%Nslip) + prm%xi_slip_sat = math_expand(prm%xi_slip_sat,prm%Nslip) + prm%H_int = math_expand(prm%H_int, prm%Nslip) + + ! sanity checks + 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(prm%xi_slip_sat <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_slip_sat' + else slipActive + allocate(prm%interaction_SlipSlip(0,0)) + allocate(prm%xi_slip_0(0)) + endif slipActive + +!-------------------------------------------------------------------------------------------------- +! twin related parameters + prm%Ntwin = config%getInts('ntwin', defaultVal=emptyIntArray) + prm%totalNtwin = sum(prm%Ntwin) + twinActive: if (prm%totalNtwin > 0) then + prm%Schmid_twin = lattice_SchmidMatrix_twin(prm%Ntwin,config%getString('lattice_structure'),& + config%getFloat('c/a',defaultVal=0.0_pReal)) + prm%interaction_TwinTwin = lattice_interaction_TwinByTwin(prm%Ntwin,& + config%getFloats('interaction_twintwin'), & + config%getString('lattice_structure')) + prm%gamma_twin_char = lattice_characteristicShear_twin(prm%Ntwin,config%getString('lattice_structure'),& + config%getFloat('c/a')) + + prm%xi_twin_0 = config%getFloats('tau0_twin',requiredSize=size(prm%Ntwin)) + + 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, prm%Ntwin) + + ! 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 + allocate(prm%interaction_TwinTwin(0,0)) + allocate(prm%xi_twin_0(0)) + allocate(prm%gamma_twin_char(0)) + endif twinActive + +!-------------------------------------------------------------------------------------------------- +! slip-twin related parameters + slipAndTwinActive: if (prm%totalNslip > 0 .and. prm%totalNtwin > 0) then + prm%h0_TwinSlip = config%getFloat('h0_twinslip') + prm%interaction_SlipTwin = lattice_interaction_SlipByTwin(prm%Nslip,prm%Ntwin,& + config%getFloats('interaction_sliptwin'), & + config%getString('lattice_structure')) + prm%interaction_TwinSlip = lattice_interaction_TwinBySlip(prm%Ntwin,prm%Nslip,& + config%getFloats('interaction_twinslip'), & + config%getString('lattice_structure')) + else slipAndTwinActive + allocate(prm%interaction_SlipTwin(prm%TotalNslip,prm%TotalNtwin)) ! at least one dimension is 0 + allocate(prm%interaction_TwinSlip(prm%TotalNtwin,prm%TotalNslip)) ! at least one dimension is 0 + prm%h0_TwinSlip = 0.0_pReal + endif slipAndTwinActive + +!-------------------------------------------------------------------------------------------------- +! exit if any parameter is out of range + if (extmsg /= '') & + call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_PHENOPOWERLAW_label//')') + +!-------------------------------------------------------------------------------------------------- +! output pararameters + outputs = config%getStrings('(output)',defaultVal=emptyStringArray) + allocate(prm%outputID(0)) + do i=1, size(outputs) + outputID = undefined_ID + select case(outputs(i)) + + case ('resistance_slip') + outputID = merge(resistance_slip_ID,undefined_ID,prm%totalNslip>0) + outputSize = prm%totalNslip + case ('accumulatedshear_slip') + outputID = merge(accumulatedshear_slip_ID,undefined_ID,prm%totalNslip>0) + outputSize = prm%totalNslip + case ('shearrate_slip') + outputID = merge(shearrate_slip_ID,undefined_ID,prm%totalNslip>0) + outputSize = prm%totalNslip + case ('resolvedstress_slip') + outputID = merge(resolvedstress_slip_ID,undefined_ID,prm%totalNslip>0) + outputSize = prm%totalNslip + + case ('resistance_twin') + outputID = merge(resistance_twin_ID,undefined_ID,prm%totalNtwin>0) + outputSize = prm%totalNtwin + case ('accumulatedshear_twin') + outputID = merge(accumulatedshear_twin_ID,undefined_ID,prm%totalNtwin>0) + outputSize = prm%totalNtwin + case ('shearrate_twin') + outputID = merge(shearrate_twin_ID,undefined_ID,prm%totalNtwin>0) + outputSize = prm%totalNtwin + case ('resolvedstress_twin') + outputID = merge(resolvedstress_twin_ID,undefined_ID,prm%totalNtwin>0) + outputSize = prm%totalNtwin + + end select + + if (outputID /= undefined_ID) then + prm%outputID = [prm%outputID, outputID] + endif + + enddo + +!-------------------------------------------------------------------------------------------------- +! allocate state arrays + NipcMyPhase = count(material_phaseAt == p) * discretization_nIP + sizeDotState = size(['tau_slip ','gamma_slip']) * prm%totalNslip & + + size(['tau_twin ','gamma_twin']) * prm%totalNtwin + sizeState = sizeDotState + + call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) + +!-------------------------------------------------------------------------------------------------- +! locally defined state aliases and initialization of state0 and aTolState + startIndex = 1 + endIndex = prm%totalNslip + stt%xi_slip => plasticState(p)%state (startIndex:endIndex,:) + stt%xi_slip = spread(prm%xi_slip_0, 2, NipcMyPhase) + dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance + + startIndex = endIndex + 1 + endIndex = endIndex + prm%totalNtwin + stt%xi_twin => plasticState(p)%state (startIndex:endIndex,:) + stt%xi_twin = spread(prm%xi_twin_0, 2, NipcMyPhase) + dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance + + startIndex = endIndex + 1 + endIndex = endIndex + prm%totalNslip + stt%gamma_slip => plasticState(p)%state (startIndex:endIndex,:) + dot%gamma_slip => plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear + ! global alias + plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%accumulatedSlip => plasticState(p)%state(startIndex:endIndex,:) + + startIndex = endIndex + 1 + endIndex = endIndex + prm%totalNtwin + stt%gamma_twin => plasticState(p)%state (startIndex:endIndex,:) + dot%gamma_twin => plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear + + plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally + + end associate + + enddo + +end subroutine plastic_phenopowerlaw_init + + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates plastic velocity gradient and its tangent +!> @details asummes that deformation by dislocation glide affects twinned and untwinned volume +! equally (Taylor assumption). Twinning happens only in untwinned volume +!-------------------------------------------------------------------------------------------------- +pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) + + real(pReal), dimension(3,3), intent(out) :: & + Lp !< plastic velocity gradient + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp !< derivative of Lp with respect to the Mandel stress + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + + integer :: & + i,k,l,m,n + real(pReal), dimension(param(instance)%totalNslip) :: & + gdot_slip_pos,gdot_slip_neg, & + dgdot_dtauslip_pos,dgdot_dtauslip_neg + real(pReal), dimension(param(instance)%totalNtwin) :: & + gdot_twin,dgdot_dtautwin + + Lp = 0.0_pReal + dLp_dMp = 0.0_pReal + + associate(prm => param(instance)) + + call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg) + slipSystems: do i = 1, prm%totalNslip + Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%Schmid_slip(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%Schmid_slip(k,l,i) * prm%nonSchmid_pos(m,n,i) & + + dgdot_dtauslip_neg(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_neg(m,n,i) + enddo slipSystems + + call kinetics_twin(Mp,instance,of,gdot_twin,dgdot_dtautwin) + twinSystems: do i = 1, prm%totalNtwin + Lp = Lp + gdot_twin(i)*prm%Schmid_twin(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%Schmid_twin(k,l,i)*prm%Schmid_twin(m,n,i) + enddo twinSystems + + end associate + +end subroutine plastic_phenopowerlaw_LpAndItsTangent + + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates the rate of change of microstructure +!-------------------------------------------------------------------------------------------------- +module subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + + real(pReal) :: & + c_SlipSlip,c_TwinSlip,c_TwinTwin, & + xi_slip_sat_offset,& + sumGamma,sumF + real(pReal), dimension(param(instance)%totalNslip) :: & + left_SlipSlip,right_SlipSlip, & + gdot_slip_pos,gdot_slip_neg + + associate(prm => param(instance), stt => state(instance), dot => dotState(instance)) + + sumGamma = sum(stt%gamma_slip(:,of)) + sumF = sum(stt%gamma_twin(:,of)/prm%gamma_twin_char) + +!-------------------------------------------------------------------------------------------------- +! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices + c_SlipSlip = prm%h0_slipslip * (1.0_pReal + prm%c_1*sumF** prm%c_2) + c_TwinSlip = prm%h0_TwinSlip * sumGamma**prm%c_3 + c_TwinTwin = prm%h0_TwinTwin * sumF**prm%c_4 + +!-------------------------------------------------------------------------------------------------- +! calculate left and right vectors + left_SlipSlip = 1.0_pReal + prm%H_int + xi_slip_sat_offset = prm%spr*sqrt(sumF) + right_SlipSlip = abs(1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset)) **prm%a_slip & + * sign(1.0_pReal,1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset)) + +!-------------------------------------------------------------------------------------------------- +! shear rates + call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg) + dot%gamma_slip(:,of) = abs(gdot_slip_pos+gdot_slip_neg) + call kinetics_twin(Mp,instance,of,dot%gamma_twin(:,of)) + +!-------------------------------------------------------------------------------------------------- +! hardening + dot%xi_slip(:,of) = c_SlipSlip * left_SlipSlip * & + matmul(prm%interaction_SlipSlip,dot%gamma_slip(:,of)*right_SlipSlip) & + + matmul(prm%interaction_SlipTwin,dot%gamma_twin(:,of)) + + dot%xi_twin(:,of) = c_TwinSlip * matmul(prm%interaction_TwinSlip,dot%gamma_slip(:,of)) & + + c_TwinTwin * matmul(prm%interaction_TwinTwin,dot%gamma_twin(:,of)) + end associate + +end subroutine plastic_phenopowerlaw_dotState + + +!-------------------------------------------------------------------------------------------------- +!> @brief writes results to HDF5 output file +!-------------------------------------------------------------------------------------------------- +module subroutine plastic_phenopowerlaw_results(instance,group) + + integer, intent(in) :: instance + character(len=*), intent(in) :: group + + integer :: o + + associate(prm => param(instance), stt => state(instance)) + outputsLoop: do o = 1,size(prm%outputID) + select case(prm%outputID(o)) + + case (resistance_slip_ID) + call results_writeDataset(group,stt%xi_slip, 'xi_sl', & + 'resistance against plastic slip','Pa') + case (accumulatedshear_slip_ID) + call results_writeDataset(group,stt%gamma_slip,'gamma_sl', & + 'plastic shear','1') + + case (resistance_twin_ID) + call results_writeDataset(group,stt%xi_twin, 'xi_tw', & + 'resistance against twinning','Pa') + case (accumulatedshear_twin_ID) + call results_writeDataset(group,stt%gamma_twin,'gamma_tw', & + 'twinning shear','1') + + end select + enddo outputsLoop + end associate + +end subroutine plastic_phenopowerlaw_results + + +!-------------------------------------------------------------------------------------------------- +!> @brief Shear rates on slip 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_slip(Mp,instance,of, & + gdot_slip_pos,gdot_slip_neg,dgdot_dtau_slip_pos,dgdot_dtau_slip_neg) + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + + real(pReal), intent(out), dimension(param(instance)%totalNslip) :: & + gdot_slip_pos, & + gdot_slip_neg + real(pReal), intent(out), optional, dimension(param(instance)%totalNslip) :: & + dgdot_dtau_slip_pos, & + dgdot_dtau_slip_neg + + real(pReal), dimension(param(instance)%totalNslip) :: & + 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%totalNslip + 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) + 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 + * sign(abs(tau_slip_pos/stt%xi_slip(:,of))**prm%n_slip, tau_slip_pos) + else where + gdot_slip_pos = 0.0_pReal + end where + + where(dNeq0(tau_slip_neg)) + gdot_slip_neg = prm%gdot0_slip * 0.5_pReal & ! only used if non-Schmid active, always 1/2 + * sign(abs(tau_slip_neg/stt%xi_slip(:,of))**prm%n_slip, tau_slip_neg) + else where + gdot_slip_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_slip/tau_slip_pos + else where + dgdot_dtau_slip_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_slip/tau_slip_neg + else where + dgdot_dtau_slip_neg = 0.0_pReal + end where + endif + end associate + +end subroutine kinetics_slip + + +!-------------------------------------------------------------------------------------------------- +!> @brief Shear rates on twin systems and their derivatives with respect to resolved stress. +! twinning is assumed to take place only in untwinned volume. +!> @details Derivates 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,instance,of,& + gdot_twin,dgdot_dtau_twin) + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + + real(pReal), dimension(param(instance)%totalNtwin), intent(out) :: & + gdot_twin + real(pReal), dimension(param(instance)%totalNtwin), intent(out), optional :: & + dgdot_dtau_twin + + real(pReal), dimension(param(instance)%totalNtwin) :: & + tau_twin + integer :: i + + associate(prm => param(instance), stt => state(instance)) + + do i = 1, prm%totalNtwin + tau_twin(i) = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,i)) + enddo + + where(tau_twin > 0.0_pReal) + gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,of)/prm%gamma_twin_char)) & ! only twin in untwinned volume fraction + * prm%gdot0_twin*(abs(tau_twin)/stt%xi_twin(:,of))**prm%n_twin + else where + gdot_twin = 0.0_pReal + end where + + if (present(dgdot_dtau_twin)) then + where(dNeq0(gdot_twin)) + dgdot_dtau_twin = gdot_twin*prm%n_twin/tau_twin + else where + dgdot_dtau_twin = 0.0_pReal + end where + endif + + end associate + +end subroutine kinetics_twin + +end submodule plastic_phenopowerlaw diff --git a/src/crystallite.f90 b/src/crystallite.f90 index eb6fc7f94..8356d0a65 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -21,7 +21,6 @@ module crystallite use constitutive use discretization use lattice - use plastic_nonlocal use results implicit none diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 84c3c4be7..08e02290b 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -160,7 +160,7 @@ subroutine homogenization_init allocate(materialpoint_converged(discretization_nIP,discretization_nElem), source=.true.) allocate(materialpoint_doneAndHappy(2,discretization_nIP,discretization_nElem), source=.true.) - write(6,'(/,a)') ' <<<+- homogenization init -+>>>'; flush(6) + write(6,'(/,a)') ' <<<+- homogenization init -+>>>'; flush(6) if (iand(debug_level(debug_homogenization), debug_levelBasic) /= 0) then write(6,'(a32,1x,7(i8,1x))') 'materialpoint_dPdF: ', shape(materialpoint_dPdF) diff --git a/src/plastic_dislotwin.f90 b/src/plastic_dislotwin.f90 deleted file mode 100644 index ae89dcc1c..000000000 --- a/src/plastic_dislotwin.f90 +++ /dev/null @@ -1,1175 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!> @author Su Leen Wong, Max-Planck-Institut für Eisenforschung GmbH -!> @author Nan Jia, Max-Planck-Institut für Eisenforschung GmbH -!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH -!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine incoprorating dislocation and twinning physics -!> @details to be done -!-------------------------------------------------------------------------------------------------- -module plastic_dislotwin - use prec - use debug - use math - use IO - use material - use config - use lattice - use discretization - use results - - implicit none - private - - 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 - D, & !< grain size - p_sb, & !< p-exponent in shear band velocity - q_sb, & !< q-exponent in shear band velocity - CEdgeDipMinDistance, & !< - i_tw, & !< - tau_0, & !< strength due to elements in solid solution - L_tw, & !< Length of twin nuclei in Burgers vectors - L_tr, & !< Length of trans nuclei in Burgers vectors - xc_twin, & !< critical distance for formation of twin nucleus - xc_trans, & !< critical distance for formation of trans nucleus - V_cs, & !< cross slip volume - sbResistance, & !< value for shearband resistance (might become an internal state variable at some point) - sbVelocity, & !< value for shearband velocity_0 - sbQedge, & !< activation energy for shear bands - SFE_0K, & !< stacking fault energy at zero K - dSFE_dT, & !< temperature dependance of stacking fault energy - aTol_rho, & !< absolute tolerance for integration of dislocation density - aTol_f_tw, & !< absolute tolerance for integration of twin volume fraction - 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 :: & - 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 - b_tw, & !< absolute length of burgers vector [m] for each twin system - b_tr, & !< absolute length of burgers vector [m] for each transformation system - Delta_F,& !< activation energy for glide [J] for each slip system - v0, & !< dislocation velocity prefactor [m/s] for each slip system - dot_N_0_tw, & !< twin nucleation rate [1/m³s] for each twin system - dot_N_0_tr, & !< trans nucleation rate [1/m³s] for each trans system - t_tw, & !< twin thickness [m] for each twin system - CLambdaSlip, & !< Adj. parameter for distance between 2 forest dislocations for each slip system - atomicVolume, & - t_tr, & !< martensite lamellar thickness [m] for each trans system and instance - p, & !< p-exponent in glide velocity - q, & !< q-exponent in glide velocity - r, & !< r-exponent in twin nucleation rate - s, & !< s-exponent in trans nucleation rate - gamma_char, & !< characteristic shear for twins - B !< drag coefficient - real(pReal), dimension(:,:), allocatable :: & - 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 :: & - n0_sl, & !< slip system normal - forestProjection, & - C66 - real(pReal), dimension(:,:,:), allocatable :: & - P_tr, & - P_sl, & - P_tw, & - C66_tw, & - C66_tr - 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 :: & - 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 - 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, & - rho_dip, & - gamma_sl, & - 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 - Lambda_tw, & !< mean free path between 2 obstacles seen by a growing twin - Lambda_tr, & !< mean free path between 2 obstacles seen by a growing martensite - tau_pass, & - tau_hat_tw, & - tau_hat_tr, & - V_tw, & !< volume of a new twin - V_tr, & !< volume of a new martensite disc - 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 - type(tDislotwinState), allocatable, dimension(:) :: & - dotState, & - state - type(tDislotwinMicrostructure), allocatable, dimension(:) :: dependentState - - public :: & - plastic_dislotwin_init, & - plastic_dislotwin_homogenizedC, & - plastic_dislotwin_dependentState, & - plastic_dislotwin_LpAndItsTangent, & - plastic_dislotwin_dotState, & - plastic_dislotwin_results - -contains - - -!-------------------------------------------------------------------------------------------------- -!> @brief module initialization -!> @details reads in material parameters, allocates arrays, and does sanity checks -!-------------------------------------------------------------------------------------------------- -subroutine plastic_dislotwin_init - - integer :: & - Ninstance, & - p, i, & - NipcMyPhase, outputSize, & - 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' - - 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)), & - dot => dotState(phase_plasticityInstance(p)), & - 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) - prm%sum_N_sl = sum(prm%N_sl) - slipActive: if (prm%sum_N_sl > 0) then - prm%P_sl = lattice_SchmidMatrix_slip(prm%N_sl,config%getString('lattice_structure'),& - config%getFloat('c/a',defaultVal=0.0_pReal)) - prm%h_sl_sl = lattice_interaction_SlipBySlip(prm%N_sl, & - config%getFloats('interaction_slipslip'), & - config%getString('lattice_structure')) - 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)) - 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)) - prm%b_sl = config%getFloats('slipburgers',requiredSize=size(prm%N_sl)) - prm%Delta_F = config%getFloats('qedge', requiredSize=size(prm%N_sl)) - prm%CLambdaSlip = config%getFloats('clambdaslip',requiredSize=size(prm%N_sl)) - prm%p = config%getFloats('p_slip', requiredSize=size(prm%N_sl)) - 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') - prm%Qsd = config%getFloat('qsd') - prm%atomicVolume = config%getFloat('catomicvolume') * prm%b_sl**3.0_pReal - prm%ExtendedDislocations = config%keyExists('/extend_dislocations/') - if (prm%ExtendedDislocations) then - 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) - prm%v0 = math_expand(prm%v0, prm%N_sl) - prm%b_sl = math_expand(prm%b_sl, prm%N_sl) - prm%Delta_F = math_expand(prm%Delta_F, prm%N_sl) - prm%CLambdaSlip = math_expand(prm%CLambdaSlip, prm%N_sl) - 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) - - ! sanity checks - if ( prm%D0 <= 0.0_pReal) extmsg = trim(extmsg)//' D0' - if ( prm%Qsd <= 0.0_pReal) extmsg = trim(extmsg)//' Qsd' - if (any(prm%rho_mob_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_mob_0' - if (any(prm%rho_dip_0 < 0.0_pReal)) extmsg = trim(extmsg)//' rho_dip_0' - 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)//' Delta_F' - if (any(prm%CLambdaSlip <= 0.0_pReal)) extmsg = trim(extmsg)//' CLambdaSlip' - 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) - prm%sum_N_tw = sum(prm%N_tw) - if (prm%sum_N_tw > 0) then - prm%P_tw = lattice_SchmidMatrix_twin(prm%N_tw,config%getString('lattice_structure'),& - config%getFloat('c/a',defaultVal=0.0_pReal)) - 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 = 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)) - allocate(prm%b_tw (0)) - allocate(prm%r (0)) - allocate(prm%h_tw_tw (0,0)) - endif - -!-------------------------------------------------------------------------------------------------- -! transformation related parameters - prm%N_tr = config%getInts('ntrans', defaultVal=emptyIntArray) - prm%sum_N_tr = sum(prm%N_tr) - 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) - endif - prm%t_tr = config%getFloats('lamellarsize') - prm%t_tr = math_expand(prm%t_tr,prm%N_tr) - prm%s = config%getFloats('s_trans',defaultVal=[0.0_pReal]) - prm%s = math_expand(prm%s,prm%N_tr) - else - allocate(prm%t_tr (0)) - allocate(prm%b_tr (0)) - 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 - 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 - -!-------------------------------------------------------------------------------------------------- -! shearband related parameters - prm%sbVelocity = config%getFloat('shearbandvelocity',defaultVal=0.0_pReal) - 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//')') - if (prm%aTol_f_tw <= 0.0_pReal) & - call IO_error(211,el=p,ext_msg='aTol_f_tw ('//PLASTICITY_DISLOTWIN_label//')') - endif - if (prm%sum_N_tr > 0) then - 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 - -!-------------------------------------------------------------------------------------------------- -! allocate state arrays - NipcMyPhase = count(material_phaseAt == p) * discretization_nIP - sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl & - + 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 - 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,:) - 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,:) - dot%gamma_sl=>plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%aTolState(startIndex:endIndex) = 1.0e6_pReal !ToDo: better make optional parameter - ! global alias - plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%accumulatedSlip => plasticState(p)%state(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 -!-------------------------------------------------------------------------------------------------- -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 & - + stt%f_tw(i,of)*prm%C66_tw(1:6,1:6,i) - enddo - do i=1,prm%sum_N_tr - 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 -!-------------------------------------------------------------------------------------------------- -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,& - BoltzmannRatio, & - ddot_gamma_dtau, & - tau - real(pReal), dimension(param(instance)%sum_N_sl) :: & - dot_gamma_sl,ddot_gamma_dtau_slip - real(pReal), dimension(param(instance)%sum_N_tw) :: & - dot_gamma_twin,ddot_gamma_dtau_twin - real(pReal), dimension(param(instance)%sum_N_tr) :: & - dot_gamma_tr,ddot_gamma_dtau_trans - real(pReal):: dot_gamma_sb - real(pReal), dimension(3,3) :: eigVectors, P_sb - real(pReal), dimension(3) :: eigValues - logical :: error - real(pReal), dimension(3,6), parameter :: & - sb_sComposition = & - reshape(real([& - 1, 0, 1, & - 1, 0,-1, & - 1, 1, 0, & - 1,-1, 0, & - 0, 1, 1, & - 0, 1,-1 & - ],pReal),[ 3,6]), & - sb_mComposition = & - reshape(real([& - 1, 0,-1, & - 1, 0,+1, & - 1,-1, 0, & - 1, 1, 0, & - 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 - - 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) - 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_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 - 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_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 - 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_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 -!-------------------------------------------------------------------------------------------------- -subroutine plastic_dislotwin_dotState(Mp,T,instance,of) - - real(pReal), dimension(3,3), intent(in):: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - T !< temperature at integration point - integer, intent(in) :: & - instance, & - of - - integer :: i - real(pReal) :: & - f_unrotated, & - VacancyDiffusion, & - rho_dip_distance, & - v_cl, & !< climb velocity - Gamma, & !< stacking fault energy - tau, & - 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, & - dot_rho_dip_climb, & - rho_dip_distance_min, & - dot_gamma_sl - real(pReal), dimension(param(instance)%sum_N_tw) :: & - dot_gamma_twin - real(pReal), dimension(param(instance)%sum_N_tr) :: & - dot_gamma_tr - - associate(prm => param(instance), stt => state(instance), & - dot => dotState(instance), dst => dependentState(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)) - VacancyDiffusion = prm%D0*exp(-prm%Qsd/(kB*T)) - - 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)) - - significantSlipStress: if (dEq0(tau)) then - dot_rho_dip_formation(i) = 0.0_pReal - dot_rho_dip_climb(i) = 0.0_pReal - else significantSlipStress - rho_dip_distance = 3.0_pReal*prm%mu*prm%b_sl(i)/(16.0_pReal*PI*abs(tau)) - rho_dip_distance = math_clip(rho_dip_distance, right = dst%Lambda_sl(i,of)) - rho_dip_distance = math_clip(rho_dip_distance, left = rho_dip_distance_min(i)) - - if (prm%dipoleFormation) then - dot_rho_dip_formation(i) = 2.0_pReal*(rho_dip_distance-rho_dip_distance_min(i))/prm%b_sl(i) & - * stt%rho_mob(i,of)*abs(dot_gamma_sl(i)) - else - dot_rho_dip_formation(i) = 0.0_pReal - endif - - if (dEq(rho_dip_distance,rho_dip_distance_min(i))) then - dot_rho_dip_climb(i) = 0.0_pReal - 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 - 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 - 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 - endif significantSlipStress - enddo slipState - - dot%rho_mob(:,of) = abs(dot_gamma_sl)/(prm%b_sl*dst%Lambda_sl(:,of)) & - - dot_rho_dip_formation & - - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_mob(:,of)*abs(dot_gamma_sl) - - dot%rho_dip(:,of) = dot_rho_dip_formation & - - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,of)*abs(dot_gamma_sl) & - - dot_rho_dip_climb - - call kinetics_twin(Mp,T,dot_gamma_sl,instance,of,dot_gamma_twin) - dot%f_tw(:,of) = f_unrotated*dot_gamma_twin/prm%gamma_char - - call kinetics_trans(Mp,T,dot_gamma_sl,instance,of,dot_gamma_tr) - dot%f_tr(:,of) = f_unrotated*dot_gamma_tr - - end associate - -end subroutine plastic_dislotwin_dotState - - -!-------------------------------------------------------------------------------------------------- -!> @brief calculates derived quantities from state -!-------------------------------------------------------------------------------------------------- -subroutine plastic_dislotwin_dependentState(T,instance,of) - - integer, intent(in) :: & - instance, & - of - real(pReal), intent(in) :: & - T - - real(pReal) :: & - sumf_twin,Gamma,sumf_trans - real(pReal), dimension(param(instance)%sum_N_sl) :: & - inv_lambda_sl_sl, & !< 1/mean free distance between 2 forest dislocations seen by a moving dislocation - inv_lambda_sl_tw, & !< 1/mean free distance between 2 twin stacks from different systems seen by a moving dislocation - inv_lambda_sl_tr !< 1/mean free distance between 2 martensite lamellar from different systems seen by a moving dislocation - real(pReal), dimension(param(instance)%sum_N_tw) :: & - inv_lambda_tw_tw, & !< 1/mean free distance between 2 twin stacks from different systems seen by a growing twin - f_over_t_tw - real(pReal), dimension(param(instance)%sum_N_tr) :: & - inv_lambda_tr_tr, & !< 1/mean free distance between 2 martensite stacks from different systems seen by a growing martensite - 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)) - else - 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) & - + 3.0_pReal*prm%b_tw*prm%mu/(prm%L_tw*prm%b_sl) ! slip burgers here correct? - 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) - - 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 -!-------------------------------------------------------------------------------------------------- -subroutine plastic_dislotwin_results(instance,group) - - integer, intent(in) :: instance - character(len=*) :: 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)) - - 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 (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_tr_ID) - call results_writeDataset(group,stt%f_tr,'f_tr',& - 'martensite volume fraction','m³/m³') - - end select - enddo outputsLoop - end associate - -end subroutine plastic_dislotwin_results - - -!-------------------------------------------------------------------------------------------------- -!> @brief Shear rates on slip systems, their derivatives with respect to resolved stress and the -! resolved stresss -!> @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 -!-------------------------------------------------------------------------------------------------- -pure subroutine kinetics_slip(Mp,T,instance,of, & - dot_gamma_sl,ddot_gamma_dtau_slip,tau_slip) - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - T !< temperature - 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) :: & - ddot_gamma_dtau_slip, & - tau_slip - real(pReal), dimension(param(instance)%sum_N_sl) :: & - ddot_gamma_dtau - - real(pReal), dimension(param(instance)%sum_N_sl) :: & - tau, & - stressRatio, & - StressRatio_p, & - BoltzmannRatio, & - v_wait_inverse, & !< inverse of the effective velocity of a dislocation waiting at obstacles (unsigned) - v_run_inverse, & !< inverse of the velocity of a free moving dislocation (unsigned) - dV_wait_inverse_dTau, & - dV_run_inverse_dTau, & - dV_dTau, & - tau_eff !< effective resolved stress - 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) & - / prm%tau_0 - dV_run_inverse_dTau = -1.0_pReal * v_run_inverse/tau_eff - dV_dTau = -1.0_pReal * (dV_wait_inverse_dTau+dV_run_inverse_dTau) & - / (v_wait_inverse+v_run_inverse)**2.0_pReal - ddot_gamma_dtau = dV_dTau*stt%rho_mob(:,of)*prm%b_sl - else where significantStress - 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 -!-------------------------------------------------------------------------------------------------- -pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,instance,of,& - dot_gamma_twin,ddot_gamma_dtau_twin) - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - T !< temperature - integer, intent(in) :: & - instance, & - 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 - s1=prm%fcc_twinNucleationSlipPair(1,i) - s2=prm%fcc_twinNucleationSlipPair(2,i) - if (tau(i) < dst%tau_r_tw(i,of)) then ! ToDo: correct? - Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,of)+stt%rho_dip(s2,of))+& - abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,of)+stt%rho_dip(s1,of)))/& ! ToDo: MD: it would be more consistent to use shearrates from state - (prm%L_tw*prm%b_sl(i))*& - (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tw(i,of)-tau(i)))) ! P_ncs - else - Ndot0=0.0_pReal - end if - else isFCC - 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) - ddot_gamma_dtau = (dot_gamma_twin*prm%r/tau)*StressRatio_r - else where significantStress - 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 -!-------------------------------------------------------------------------------------------------- -pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,instance,of,& - dot_gamma_tr,ddot_gamma_dtau_trans) - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - T !< temperature - integer, intent(in) :: & - instance, & - 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 - s1=prm%fcc_twinNucleationSlipPair(1,i) - s2=prm%fcc_twinNucleationSlipPair(2,i) - if (tau(i) < dst%tau_r_tr(i,of)) then ! ToDo: correct? - Ndot0=(abs(dot_gamma_sl(s1))*(stt%rho_mob(s2,of)+stt%rho_dip(s2,of))+& - abs(dot_gamma_sl(s2))*(stt%rho_mob(s1,of)+stt%rho_dip(s1,of)))/& ! ToDo: MD: it would be more consistent to use shearrates from state - (prm%L_tr*prm%b_sl(i))*& - (1.0_pReal-exp(-prm%V_cs/(kB*T)*(dst%tau_r_tr(i,of)-tau(i)))) ! P_ncs - else - Ndot0=0.0_pReal - end if - else isFCC - 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) - ddot_gamma_dtau = (dot_gamma_tr*prm%s/tau)*StressRatio_s - else where significantStress - 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 - -end module plastic_dislotwin diff --git a/src/plastic_kinematichardening.f90 b/src/plastic_kinematichardening.f90 deleted file mode 100644 index 721073746..000000000 --- a/src/plastic_kinematichardening.f90 +++ /dev/null @@ -1,541 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Philip Eisenlohr, Michigan State University -!> @author Zhuowen Zhao, Michigan State University -!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!> @brief Phenomenological crystal plasticity using a power law formulation for the shear rates -!! and a Voce-type kinematic hardening rule -!-------------------------------------------------------------------------------------------------- -module plastic_kinehardening - use prec - use debug - use math - use IO - use material - use config - use lattice - use discretization - use results - - implicit none - private - - enum, bind(c) - enumerator :: & - undefined_ID, & - crss_ID, & !< critical resolved stress - crss_back_ID, & !< critical resolved back stress - sense_ID, & !< sense of acting shear stress (-1 or +1) - chi0_ID, & !< backstress at last switch of stress sense (positive?) - gamma0_ID, & !< accumulated shear at last switch of stress sense (at current switch?) - accshear_ID, & - shearrate_ID, & - resolvedstress_ID - end enum - - type :: tParameters - real(pReal) :: & - gdot0, & !< reference shear strain rate for slip - n, & !< stress exponent for slip - aTolResistance, & - aTolShear - real(pReal), allocatable, dimension(:) :: & - crss0, & !< initial critical shear stress for slip - theta0, & !< initial hardening rate of forward stress for each slip - theta1, & !< asymptotic hardening rate of forward stress for each slip - theta0_b, & !< initial hardening rate of back stress for each slip - theta1_b, & !< asymptotic hardening rate of back stress for each slip - tau1, & - tau1_b, & - nonSchmidCoeff - real(pReal), allocatable, dimension(:,:) :: & - interaction_slipslip !< slip resistance from slip activity - real(pReal), allocatable, dimension(:,:,:) :: & - Schmid, & - nonSchmid_pos, & - nonSchmid_neg - integer :: & - totalNslip, & !< total number of active slip system - of_debug = 0 - integer, allocatable, dimension(:) :: & - Nslip !< number of active slip systems for each family - integer(kind(undefined_ID)), allocatable, dimension(:) :: & - outputID !< ID of each post result output - end type tParameters - - type :: tKinehardeningState - real(pReal), pointer, dimension(:,:) :: & !< vectors along NipcMyInstance - crss, & !< critical resolved stress - crss_back, & !< critical resolved back stress - sense, & !< sense of acting shear stress (-1 or +1) - chi0, & !< backstress at last switch of stress sense - gamma0, & !< accumulated shear at last switch of stress sense - accshear !< accumulated (absolute) shear - end type tKinehardeningState - -!-------------------------------------------------------------------------------------------------- -! containers for parameters and state - type(tParameters), allocatable, dimension(:) :: param - type(tKinehardeningState), allocatable, dimension(:) :: & - dotState, & - deltaState, & - state - - public :: & - plastic_kinehardening_init, & - plastic_kinehardening_LpAndItsTangent, & - plastic_kinehardening_dotState, & - plastic_kinehardening_deltaState, & - plastic_kinehardening_results - -contains - - -!-------------------------------------------------------------------------------------------------- -!> @brief module initialization -!> @details reads in material parameters, allocates arrays, and does sanity checks -!-------------------------------------------------------------------------------------------------- -subroutine plastic_kinehardening_init - - integer :: & - Ninstance, & - p, i, o, & - NipcMyPhase, & - sizeState, sizeDeltaState, sizeDotState, & - startIndex, endIndex - - integer(kind(undefined_ID)) :: & - outputID - - character(len=pStringLen) :: & - extmsg = '' - character(len=pStringLen), dimension(:), allocatable :: & - outputs - - write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_KINEHARDENING_label//' init -+>>>' - - Ninstance = count(phase_plasticity == PLASTICITY_KINEHARDENING_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(deltaState(Ninstance)) - - do p = 1, size(phase_plasticityInstance) - if (phase_plasticity(p) /= PLASTICITY_KINEHARDENING_ID) cycle - associate(prm => param(phase_plasticityInstance(p)), & - dot => dotState(phase_plasticityInstance(p)), & - dlt => deltaState(phase_plasticityInstance(p)), & - stt => state(phase_plasticityInstance(p)),& - config => config_phase(p)) - -#ifdef DEBUG - if (p==material_phaseAt(debug_g,debug_e)) then - prm%of_debug = material_phasememberAt(debug_g,debug_i,debug_e) - endif -#endif - -!-------------------------------------------------------------------------------------------------- -! optional parameters that need to be defined - prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal) - prm%aTolShear = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal) - - ! sanity checks - if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//' aTolresistance' - if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//' aTolShear' - -!-------------------------------------------------------------------------------------------------- -! slip related parameters - prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) - prm%totalNslip = sum(prm%Nslip) - slipActive: if (prm%totalNslip > 0) then - prm%Schmid = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),& - 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(prm%Nslip,prm%nonSchmidCoeff,+1) - prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1) - else - prm%nonSchmid_pos = prm%Schmid - prm%nonSchmid_neg = prm%Schmid - endif - prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, & - config%getFloats('interaction_slipslip'), & - config%getString('lattice_structure')) - - prm%crss0 = config%getFloats('crss0', requiredSize=size(prm%Nslip)) - prm%tau1 = config%getFloats('tau1', requiredSize=size(prm%Nslip)) - prm%tau1_b = config%getFloats('tau1_b', requiredSize=size(prm%Nslip)) - prm%theta0 = config%getFloats('theta0', requiredSize=size(prm%Nslip)) - prm%theta1 = config%getFloats('theta1', requiredSize=size(prm%Nslip)) - prm%theta0_b = config%getFloats('theta0_b', requiredSize=size(prm%Nslip)) - prm%theta1_b = config%getFloats('theta1_b', requiredSize=size(prm%Nslip)) - - prm%gdot0 = config%getFloat('gdot0') - prm%n = config%getFloat('n_slip') - - ! expand: family => system - prm%crss0 = math_expand(prm%crss0, prm%Nslip) - prm%tau1 = math_expand(prm%tau1, prm%Nslip) - prm%tau1_b = math_expand(prm%tau1_b, prm%Nslip) - prm%theta0 = math_expand(prm%theta0, prm%Nslip) - prm%theta1 = math_expand(prm%theta1, prm%Nslip) - prm%theta0_b = math_expand(prm%theta0_b,prm%Nslip) - prm%theta1_b = math_expand(prm%theta1_b,prm%Nslip) - - - -!-------------------------------------------------------------------------------------------------- -! sanity checks - if ( prm%gdot0 <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0' - if ( prm%n <= 0.0_pReal) extmsg = trim(extmsg)//' n_slip' - if (any(prm%crss0 <= 0.0_pReal)) extmsg = trim(extmsg)//' crss0' - if (any(prm%tau1 <= 0.0_pReal)) extmsg = trim(extmsg)//' tau1' - if (any(prm%tau1_b <= 0.0_pReal)) extmsg = trim(extmsg)//' tau1_b' - - !ToDo: Any sensible checks for theta? - - endif slipActive - -!-------------------------------------------------------------------------------------------------- -! exit if any parameter is out of range - if (extmsg /= '') & - call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_KINEHARDENING_label//')') - -!-------------------------------------------------------------------------------------------------- -! output pararameters - outputs = config%getStrings('(output)',defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - do i=1, size(outputs) - outputID = undefined_ID - select case(outputs(i)) - - case ('resistance') - outputID = merge(crss_ID,undefined_ID,prm%totalNslip>0) - case ('accumulatedshear') - outputID = merge(accshear_ID,undefined_ID,prm%totalNslip>0) - case ('shearrate') - outputID = merge(shearrate_ID,undefined_ID,prm%totalNslip>0) - case ('resolvedstress') - outputID = merge(resolvedstress_ID,undefined_ID,prm%totalNslip>0) - case ('backstress') - outputID = merge(crss_back_ID,undefined_ID,prm%totalNslip>0) - case ('sense') - outputID = merge(sense_ID,undefined_ID,prm%totalNslip>0) - case ('chi0') - outputID = merge(chi0_ID,undefined_ID,prm%totalNslip>0) - case ('gamma0') - outputID = merge(gamma0_ID,undefined_ID,prm%totalNslip>0) - - end select - - if (outputID /= undefined_ID) then - prm%outputID = [prm%outputID , outputID] - endif - - enddo - -!-------------------------------------------------------------------------------------------------- -! allocate state arrays - NipcMyPhase = count(material_phaseAt == p) * discretization_nIP - sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%totalNslip - sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%totalNslip - sizeState = sizeDotState + sizeDeltaState - - call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) - -!-------------------------------------------------------------------------------------------------- -! locally defined state aliases and initialization of state0 and aTolState - startIndex = 1 - endIndex = prm%totalNslip - stt%crss => plasticState(p)%state (startIndex:endIndex,:) - stt%crss = spread(prm%crss0, 2, NipcMyPhase) - dot%crss => plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance - - startIndex = endIndex + 1 - endIndex = endIndex + prm%totalNslip - stt%crss_back => plasticState(p)%state (startIndex:endIndex,:) - dot%crss_back => plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance - - startIndex = endIndex + 1 - endIndex = endIndex + prm%totalNslip - stt%accshear => plasticState(p)%state (startIndex:endIndex,:) - dot%accshear => plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear - ! global alias - plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%accumulatedSlip => plasticState(p)%state(startIndex:endIndex,:) - - o = plasticState(p)%offsetDeltaState - startIndex = endIndex + 1 - endIndex = endIndex + prm%totalNslip - stt%sense => plasticState(p)%state (startIndex :endIndex ,:) - dlt%sense => plasticState(p)%deltaState(startIndex-o:endIndex-o,:) - - startIndex = endIndex + 1 - endIndex = endIndex + prm%totalNslip - stt%chi0 => plasticState(p)%state (startIndex :endIndex ,:) - dlt%chi0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:) - - startIndex = endIndex + 1 - endIndex = endIndex + prm%totalNslip - stt%gamma0 => plasticState(p)%state (startIndex :endIndex ,:) - dlt%gamma0 => plasticState(p)%deltaState(startIndex-o:endIndex-o,:) - - plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally - - end associate - - enddo - -end subroutine plastic_kinehardening_init - - -!-------------------------------------------------------------------------------------------------- -!> @brief calculates plastic velocity gradient and its tangent -!-------------------------------------------------------------------------------------------------- -pure subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) - - real(pReal), dimension(3,3), intent(out) :: & - Lp !< plastic velocity gradient - real(pReal), dimension(3,3,3,3), intent(out) :: & - dLp_dMp !< derivative of Lp with respect to the Mandel stress - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - - integer :: & - i,k,l,m,n - real(pReal), dimension(param(instance)%totalNslip) :: & - gdot_pos,gdot_neg, & - dgdot_dtau_pos,dgdot_dtau_neg - - Lp = 0.0_pReal - dLp_dMp = 0.0_pReal - - associate(prm => param(instance)) - - call kinetics(Mp,instance,of,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) - - do i = 1, prm%totalNslip - Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%Schmid(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_dtau_pos(i) * prm%Schmid(k,l,i) * prm%nonSchmid_pos(m,n,i) & - + dgdot_dtau_neg(i) * prm%Schmid(k,l,i) * prm%nonSchmid_neg(m,n,i) - enddo - - end associate - -end subroutine plastic_kinehardening_LpAndItsTangent - - -!-------------------------------------------------------------------------------------------------- -!> @brief calculates the rate of change of microstructure -!-------------------------------------------------------------------------------------------------- -subroutine plastic_kinehardening_dotState(Mp,instance,of) - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - - real(pReal) :: & - sumGamma - real(pReal), dimension(param(instance)%totalNslip) :: & - gdot_pos,gdot_neg - - - associate(prm => param(instance), stt => state(instance), dot => dotState(instance)) - - call kinetics(Mp,instance,of,gdot_pos,gdot_neg) - dot%accshear(:,of) = abs(gdot_pos+gdot_neg) - sumGamma = sum(stt%accshear(:,of)) - - - dot%crss(:,of) = matmul(prm%interaction_SlipSlip,dot%accshear(:,of)) & - * ( prm%theta1 & - + (prm%theta0 - prm%theta1 + prm%theta0*prm%theta1*sumGamma/prm%tau1) & - * exp(-sumGamma*prm%theta0/prm%tau1) & - ) - - dot%crss_back(:,of) = stt%sense(:,of)*dot%accshear(:,of) * & - ( prm%theta1_b + & - (prm%theta0_b - prm%theta1_b & - + prm%theta0_b*prm%theta1_b/(prm%tau1_b+stt%chi0(:,of))*(stt%accshear(:,of)-stt%gamma0(:,of))& - ) *exp(-(stt%accshear(:,of)-stt%gamma0(:,of)) *prm%theta0_b/(prm%tau1_b+stt%chi0(:,of))) & - ) - - end associate - -end subroutine plastic_kinehardening_dotState - - -!-------------------------------------------------------------------------------------------------- -!> @brief calculates (instantaneous) incremental change of microstructure -!-------------------------------------------------------------------------------------------------- -subroutine plastic_kinehardening_deltaState(Mp,instance,of) - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - - real(pReal), dimension(param(instance)%totalNslip) :: & - gdot_pos,gdot_neg, & - sense - - associate(prm => param(instance), stt => state(instance), dlt => deltaState(instance)) - - call kinetics(Mp,instance,of,gdot_pos,gdot_neg) - sense = merge(state(instance)%sense(:,of), & ! keep existing... - sign(1.0_pReal,gdot_pos+gdot_neg), & ! ...or have a defined - dEq0(gdot_pos+gdot_neg,1e-10_pReal)) ! current sense of shear direction - -#ifdef DEBUG - if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0 & - .and. (of == prm%of_debug & - .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0)) then - write(6,'(a)') '======= kinehardening delta state =======' - write(6,*) sense,state(instance)%sense(:,of) - endif -#endif - -!-------------------------------------------------------------------------------------------------- -! switch in sense of shear? - where(dNeq(sense,stt%sense(:,of),0.1_pReal)) - dlt%sense (:,of) = sense - stt%sense(:,of) ! switch sense - dlt%chi0 (:,of) = abs(stt%crss_back(:,of)) - stt%chi0(:,of) ! remember current backstress magnitude - dlt%gamma0(:,of) = stt%accshear(:,of) - stt%gamma0(:,of) ! remember current accumulated shear - else where - dlt%sense (:,of) = 0.0_pReal - dlt%chi0 (:,of) = 0.0_pReal - dlt%gamma0(:,of) = 0.0_pReal - end where - - end associate - -end subroutine plastic_kinehardening_deltaState - - -!-------------------------------------------------------------------------------------------------- -!> @brief writes results to HDF5 output file -!-------------------------------------------------------------------------------------------------- -subroutine plastic_kinehardening_results(instance,group) - - integer, intent(in) :: instance - character(len=*) :: group - integer :: o - - associate(prm => param(instance), stt => state(instance)) - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) - case (crss_ID) - call results_writeDataset(group,stt%crss,'xi_sl', & - 'resistance against plastic slip','Pa') - - case(crss_back_ID) - call results_writeDataset(group,stt%crss_back,'tau_back', & - 'back stress against plastic slip','Pa') - - case (sense_ID) - call results_writeDataset(group,stt%sense,'sense_of_shear','tbd','1') - - case (chi0_ID) - call results_writeDataset(group,stt%chi0,'chi0','tbd','Pa') - - case (gamma0_ID) - call results_writeDataset(group,stt%gamma0,'gamma0','tbd','1') - - case (accshear_ID) - call results_writeDataset(group,stt%accshear,'gamma_sl', & - 'plastic shear','1') - - end select - enddo outputsLoop - end associate - -end subroutine plastic_kinehardening_results - - -!-------------------------------------------------------------------------------------------------- -!> @brief calculates shear rates on slip systems and derivatives with respect to resolved stress -!> @details: Shear rates 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(Mp,instance,of, & - gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg) - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - - real(pReal), intent(out), dimension(param(instance)%totalNslip) :: & - gdot_pos, & - gdot_neg - real(pReal), intent(out), optional, dimension(param(instance)%totalNslip) :: & - dgdot_dtau_pos, & - dgdot_dtau_neg - - real(pReal), dimension(param(instance)%totalNslip) :: & - tau_pos, & - tau_neg - integer :: i - logical :: nonSchmidActive - - associate(prm => param(instance), stt => state(instance)) - - nonSchmidActive = size(prm%nonSchmidCoeff) > 0 - - do i = 1, prm%totalNslip - tau_pos(i) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - stt%crss_back(i,of) - tau_neg(i) = merge(math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,i)) - stt%crss_back(i,of), & - 0.0_pReal, nonSchmidActive) - enddo - - where(dNeq0(tau_pos)) - gdot_pos = prm%gdot0 * merge(0.5_pReal,1.0_pReal, nonSchmidActive) & ! 1/2 if non-Schmid active - * sign(abs(tau_pos/stt%crss(:,of))**prm%n, tau_pos) - else where - gdot_pos = 0.0_pReal - end where - - where(dNeq0(tau_neg)) - gdot_neg = prm%gdot0 * 0.5_pReal & ! only used if non-Schmid active, always 1/2 - * sign(abs(tau_neg/stt%crss(:,of))**prm%n, tau_neg) - else where - gdot_neg = 0.0_pReal - end where - - if (present(dgdot_dtau_pos)) then - where(dNeq0(gdot_pos)) - dgdot_dtau_pos = gdot_pos*prm%n/tau_pos - else where - dgdot_dtau_pos = 0.0_pReal - end where - endif - if (present(dgdot_dtau_neg)) then - where(dNeq0(gdot_neg)) - dgdot_dtau_neg = gdot_neg*prm%n/tau_neg - else where - dgdot_dtau_neg = 0.0_pReal - end where - endif - end associate - -end subroutine kinetics - -end module plastic_kinehardening diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 deleted file mode 100644 index a8e459f63..000000000 --- a/src/plastic_phenopowerlaw.f90 +++ /dev/null @@ -1,614 +0,0 @@ -!-------------------------------------------------------------------------------------------------- -!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH -!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH -!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!> @brief phenomenological crystal plasticity formulation using a powerlaw fitting -!-------------------------------------------------------------------------------------------------- -module plastic_phenopowerlaw - use prec - use debug - use math - use IO - use material - use config - use lattice - use discretization - use results - - implicit none - private - - enum, bind(c) - enumerator :: & - undefined_ID, & - resistance_slip_ID, & - accumulatedshear_slip_ID, & - shearrate_slip_ID, & - resolvedstress_slip_ID, & - resistance_twin_ID, & - accumulatedshear_twin_ID, & - shearrate_twin_ID, & - resolvedstress_twin_ID - end enum - - 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 - twinB, & - twinC, & - twinD, & - twinE, & - h0_SlipSlip, & !< reference hardening slip - slip - h0_TwinSlip, & !< reference hardening twin - slip - h0_TwinTwin, & !< reference hardening twin - twin - a_slip, & - aTolResistance, & !< absolute tolerance for integration of xi - aTolShear, & !< absolute tolerance for integration of gamma - aTolTwinfrac !< absolute tolerance for integration of f - 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(:,:) :: & - interaction_SlipSlip, & !< slip resistance from slip activity - interaction_SlipTwin, & !< slip resistance from twin activity - interaction_TwinSlip, & !< twin resistance from slip activity - interaction_TwinTwin !< twin resistance from twin activity - real(pReal), allocatable, dimension(:,:,:) :: & - Schmid_slip, & - Schmid_twin, & - nonSchmid_pos, & - nonSchmid_neg - integer :: & - totalNslip, & !< total number of active slip system - totalNtwin !< total number of active twin systems - integer, allocatable, dimension(:) :: & - Nslip, & !< number of active slip systems for each family - Ntwin !< number of active twin systems for each family - integer(kind(undefined_ID)), allocatable, dimension(:) :: & - outputID !< ID of each post result output - end type tParameters - - type :: tPhenopowerlawState - real(pReal), pointer, dimension(:,:) :: & - xi_slip, & - xi_twin, & - gamma_slip, & - gamma_twin - end type tPhenopowerlawState - -!-------------------------------------------------------------------------------------------------- -! containers for parameters and state - type(tParameters), allocatable, dimension(:) :: param - type(tPhenopowerlawState), allocatable, dimension(:) :: & - dotState, & - state - - public :: & - plastic_phenopowerlaw_init, & - plastic_phenopowerlaw_LpAndItsTangent, & - plastic_phenopowerlaw_dotState, & - plastic_phenopowerlaw_results - -contains - - -!-------------------------------------------------------------------------------------------------- -!> @brief module initialization -!> @details reads in material parameters, allocates arrays, and does sanity checks -!-------------------------------------------------------------------------------------------------- -subroutine plastic_phenopowerlaw_init - - integer :: & - Ninstance, & - p, i, & - NipcMyPhase, outputSize, & - sizeState, sizeDotState, & - startIndex, endIndex - - integer(kind(undefined_ID)) :: & - outputID - - character(len=pStringLen) :: & - extmsg = '' - character(len=pStringLen), dimension(:), allocatable :: & - outputs - - write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_PHENOPOWERLAW_label//' init -+>>>' - - Ninstance = count(phase_plasticity == PLASTICITY_PHENOPOWERLAW_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)) - - do p = 1, size(phase_plasticity) - if (phase_plasticity(p) /= PLASTICITY_PHENOPOWERLAW_ID) cycle - associate(prm => param(phase_plasticityInstance(p)), & - dot => dotState(phase_plasticityInstance(p)), & - stt => state(phase_plasticityInstance(p)), & - config => config_phase(p)) - -!-------------------------------------------------------------------------------------------------- -! optional parameters that need to be defined - prm%twinB = config%getFloat('twin_b',defaultVal=1.0_pReal) - prm%twinC = config%getFloat('twin_c',defaultVal=0.0_pReal) - prm%twinD = config%getFloat('twin_d',defaultVal=0.0_pReal) - prm%twinE = config%getFloat('twin_e',defaultVal=0.0_pReal) - - prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal) - prm%aTolShear = config%getFloat('atol_shear', defaultVal=1.0e-6_pReal) - prm%aTolTwinfrac = config%getFloat('atol_twinfrac', defaultVal=1.0e-6_pReal) - - ! sanity checks - if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//' aTolresistance' - if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//' aTolShear' - if (prm%aTolTwinfrac <= 0.0_pReal) extmsg = trim(extmsg)//' atoltwinfrac' - -!-------------------------------------------------------------------------------------------------- -! slip related parameters - prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) - prm%totalNslip = sum(prm%Nslip) - slipActive: if (prm%totalNslip > 0) then - prm%Schmid_slip = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),& - 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(prm%Nslip,prm%nonSchmidCoeff,+1) - prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1) - else - allocate(prm%nonSchmidCoeff(0)) - prm%nonSchmid_pos = prm%Schmid_slip - prm%nonSchmid_neg = prm%Schmid_slip - endif - prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(prm%Nslip, & - config%getFloats('interaction_slipslip'), & - config%getString('lattice_structure')) - - prm%xi_slip_0 = config%getFloats('tau0_slip', requiredSize=size(prm%Nslip)) - prm%xi_slip_sat = config%getFloats('tausat_slip', requiredSize=size(prm%Nslip)) - prm%H_int = config%getFloats('h_int', requiredSize=size(prm%Nslip), & - defaultVal=[(0.0_pReal,i=1,size(prm%Nslip))]) - - prm%gdot0_slip = config%getFloat('gdot0_slip') - prm%n_slip = config%getFloat('n_slip') - prm%a_slip = config%getFloat('a_slip') - prm%h0_SlipSlip = config%getFloat('h0_slipslip') - - ! expand: family => system - prm%xi_slip_0 = math_expand(prm%xi_slip_0, prm%Nslip) - prm%xi_slip_sat = math_expand(prm%xi_slip_sat,prm%Nslip) - prm%H_int = math_expand(prm%H_int, prm%Nslip) - - ! sanity checks - 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(prm%xi_slip_sat <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_slip_sat' - else slipActive - allocate(prm%interaction_SlipSlip(0,0)) - allocate(prm%xi_slip_0(0)) - endif slipActive - -!-------------------------------------------------------------------------------------------------- -! twin related parameters - prm%Ntwin = config%getInts('ntwin', defaultVal=emptyIntArray) - prm%totalNtwin = sum(prm%Ntwin) - twinActive: if (prm%totalNtwin > 0) then - prm%Schmid_twin = lattice_SchmidMatrix_twin(prm%Ntwin,config%getString('lattice_structure'),& - config%getFloat('c/a',defaultVal=0.0_pReal)) - prm%interaction_TwinTwin = lattice_interaction_TwinByTwin(prm%Ntwin,& - config%getFloats('interaction_twintwin'), & - config%getString('lattice_structure')) - prm%gamma_twin_char = lattice_characteristicShear_twin(prm%Ntwin,config%getString('lattice_structure'),& - config%getFloat('c/a')) - - prm%xi_twin_0 = config%getFloats('tau0_twin',requiredSize=size(prm%Ntwin)) - - 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, prm%Ntwin) - - ! 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 - allocate(prm%interaction_TwinTwin(0,0)) - allocate(prm%xi_twin_0(0)) - allocate(prm%gamma_twin_char(0)) - endif twinActive - -!-------------------------------------------------------------------------------------------------- -! slip-twin related parameters - slipAndTwinActive: if (prm%totalNslip > 0 .and. prm%totalNtwin > 0) then - prm%h0_TwinSlip = config%getFloat('h0_twinslip') - prm%interaction_SlipTwin = lattice_interaction_SlipByTwin(prm%Nslip,prm%Ntwin,& - config%getFloats('interaction_sliptwin'), & - config%getString('lattice_structure')) - prm%interaction_TwinSlip = lattice_interaction_TwinBySlip(prm%Ntwin,prm%Nslip,& - config%getFloats('interaction_twinslip'), & - config%getString('lattice_structure')) - else slipAndTwinActive - allocate(prm%interaction_SlipTwin(prm%TotalNslip,prm%TotalNtwin)) ! at least one dimension is 0 - allocate(prm%interaction_TwinSlip(prm%TotalNtwin,prm%TotalNslip)) ! at least one dimension is 0 - prm%h0_TwinSlip = 0.0_pReal - endif slipAndTwinActive - -!-------------------------------------------------------------------------------------------------- -! exit if any parameter is out of range - if (extmsg /= '') & - call IO_error(211,ext_msg=trim(extmsg)//'('//PLASTICITY_PHENOPOWERLAW_label//')') - -!-------------------------------------------------------------------------------------------------- -! output pararameters - outputs = config%getStrings('(output)',defaultVal=emptyStringArray) - allocate(prm%outputID(0)) - do i=1, size(outputs) - outputID = undefined_ID - select case(outputs(i)) - - case ('resistance_slip') - outputID = merge(resistance_slip_ID,undefined_ID,prm%totalNslip>0) - outputSize = prm%totalNslip - case ('accumulatedshear_slip') - outputID = merge(accumulatedshear_slip_ID,undefined_ID,prm%totalNslip>0) - outputSize = prm%totalNslip - case ('shearrate_slip') - outputID = merge(shearrate_slip_ID,undefined_ID,prm%totalNslip>0) - outputSize = prm%totalNslip - case ('resolvedstress_slip') - outputID = merge(resolvedstress_slip_ID,undefined_ID,prm%totalNslip>0) - outputSize = prm%totalNslip - - case ('resistance_twin') - outputID = merge(resistance_twin_ID,undefined_ID,prm%totalNtwin>0) - outputSize = prm%totalNtwin - case ('accumulatedshear_twin') - outputID = merge(accumulatedshear_twin_ID,undefined_ID,prm%totalNtwin>0) - outputSize = prm%totalNtwin - case ('shearrate_twin') - outputID = merge(shearrate_twin_ID,undefined_ID,prm%totalNtwin>0) - outputSize = prm%totalNtwin - case ('resolvedstress_twin') - outputID = merge(resolvedstress_twin_ID,undefined_ID,prm%totalNtwin>0) - outputSize = prm%totalNtwin - - end select - - if (outputID /= undefined_ID) then - prm%outputID = [prm%outputID, outputID] - endif - - enddo - -!-------------------------------------------------------------------------------------------------- -! allocate state arrays - NipcMyPhase = count(material_phaseAt == p) * discretization_nIP - sizeDotState = size(['tau_slip ','gamma_slip']) * prm%totalNslip & - + size(['tau_twin ','gamma_twin']) * prm%totalNtwin - sizeState = sizeDotState - - call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0) - -!-------------------------------------------------------------------------------------------------- -! locally defined state aliases and initialization of state0 and aTolState - startIndex = 1 - endIndex = prm%totalNslip - stt%xi_slip => plasticState(p)%state (startIndex:endIndex,:) - stt%xi_slip = spread(prm%xi_slip_0, 2, NipcMyPhase) - dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance - - startIndex = endIndex + 1 - endIndex = endIndex + prm%totalNtwin - stt%xi_twin => plasticState(p)%state (startIndex:endIndex,:) - stt%xi_twin = spread(prm%xi_twin_0, 2, NipcMyPhase) - dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolResistance - - startIndex = endIndex + 1 - endIndex = endIndex + prm%totalNslip - stt%gamma_slip => plasticState(p)%state (startIndex:endIndex,:) - dot%gamma_slip => plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear - ! global alias - plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%accumulatedSlip => plasticState(p)%state(startIndex:endIndex,:) - - startIndex = endIndex + 1 - endIndex = endIndex + prm%totalNtwin - stt%gamma_twin => plasticState(p)%state (startIndex:endIndex,:) - dot%gamma_twin => plasticState(p)%dotState(startIndex:endIndex,:) - plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolShear - - plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally - - end associate - - enddo - -end subroutine plastic_phenopowerlaw_init - - -!-------------------------------------------------------------------------------------------------- -!> @brief calculates plastic velocity gradient and its tangent -!> @details asummes that deformation by dislocation glide affects twinned and untwinned volume -! equally (Taylor assumption). Twinning happens only in untwinned volume -!-------------------------------------------------------------------------------------------------- -pure subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) - - real(pReal), dimension(3,3), intent(out) :: & - Lp !< plastic velocity gradient - real(pReal), dimension(3,3,3,3), intent(out) :: & - dLp_dMp !< derivative of Lp with respect to the Mandel stress - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - - integer :: & - i,k,l,m,n - real(pReal), dimension(param(instance)%totalNslip) :: & - gdot_slip_pos,gdot_slip_neg, & - dgdot_dtauslip_pos,dgdot_dtauslip_neg - real(pReal), dimension(param(instance)%totalNtwin) :: & - gdot_twin,dgdot_dtautwin - - Lp = 0.0_pReal - dLp_dMp = 0.0_pReal - - associate(prm => param(instance)) - - call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg) - slipSystems: do i = 1, prm%totalNslip - Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%Schmid_slip(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%Schmid_slip(k,l,i) * prm%nonSchmid_pos(m,n,i) & - + dgdot_dtauslip_neg(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_neg(m,n,i) - enddo slipSystems - - call kinetics_twin(Mp,instance,of,gdot_twin,dgdot_dtautwin) - twinSystems: do i = 1, prm%totalNtwin - Lp = Lp + gdot_twin(i)*prm%Schmid_twin(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%Schmid_twin(k,l,i)*prm%Schmid_twin(m,n,i) - enddo twinSystems - - end associate - -end subroutine plastic_phenopowerlaw_LpAndItsTangent - - -!-------------------------------------------------------------------------------------------------- -!> @brief calculates the rate of change of microstructure -!-------------------------------------------------------------------------------------------------- -subroutine plastic_phenopowerlaw_dotState(Mp,instance,of) - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - - real(pReal) :: & - c_SlipSlip,c_TwinSlip,c_TwinTwin, & - xi_slip_sat_offset,& - sumGamma,sumF - real(pReal), dimension(param(instance)%totalNslip) :: & - left_SlipSlip,right_SlipSlip, & - gdot_slip_pos,gdot_slip_neg - - associate(prm => param(instance), stt => state(instance), dot => dotState(instance)) - - sumGamma = sum(stt%gamma_slip(:,of)) - sumF = sum(stt%gamma_twin(:,of)/prm%gamma_twin_char) - -!-------------------------------------------------------------------------------------------------- -! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices - c_SlipSlip = prm%h0_slipslip * (1.0_pReal + prm%twinC*sumF** prm%twinB) - c_TwinSlip = prm%h0_TwinSlip * sumGamma**prm%twinE - c_TwinTwin = prm%h0_TwinTwin * sumF**prm%twinD - -!-------------------------------------------------------------------------------------------------- -! calculate left and right vectors - left_SlipSlip = 1.0_pReal + prm%H_int - xi_slip_sat_offset = prm%spr*sqrt(sumF) - right_SlipSlip = abs(1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset)) **prm%a_slip & - * sign(1.0_pReal,1.0_pReal-stt%xi_slip(:,of) / (prm%xi_slip_sat+xi_slip_sat_offset)) - -!-------------------------------------------------------------------------------------------------- -! shear rates - call kinetics_slip(Mp,instance,of,gdot_slip_pos,gdot_slip_neg) - dot%gamma_slip(:,of) = abs(gdot_slip_pos+gdot_slip_neg) - call kinetics_twin(Mp,instance,of,dot%gamma_twin(:,of)) - -!-------------------------------------------------------------------------------------------------- -! hardening - dot%xi_slip(:,of) = c_SlipSlip * left_SlipSlip * & - matmul(prm%interaction_SlipSlip,dot%gamma_slip(:,of)*right_SlipSlip) & - + matmul(prm%interaction_SlipTwin,dot%gamma_twin(:,of)) - - dot%xi_twin(:,of) = c_TwinSlip * matmul(prm%interaction_TwinSlip,dot%gamma_slip(:,of)) & - + c_TwinTwin * matmul(prm%interaction_TwinTwin,dot%gamma_twin(:,of)) - end associate - -end subroutine plastic_phenopowerlaw_dotState - - -!-------------------------------------------------------------------------------------------------- -!> @brief writes results to HDF5 output file -!-------------------------------------------------------------------------------------------------- -subroutine plastic_phenopowerlaw_results(instance,group) - - integer, intent(in) :: instance - character(len=*), intent(in) :: group - - integer :: o - - associate(prm => param(instance), stt => state(instance)) - outputsLoop: do o = 1,size(prm%outputID) - select case(prm%outputID(o)) - - case (resistance_slip_ID) - call results_writeDataset(group,stt%xi_slip, 'xi_sl', & - 'resistance against plastic slip','Pa') - case (accumulatedshear_slip_ID) - call results_writeDataset(group,stt%gamma_slip,'gamma_sl', & - 'plastic shear','1') - - case (resistance_twin_ID) - call results_writeDataset(group,stt%xi_twin, 'xi_tw', & - 'resistance against twinning','Pa') - case (accumulatedshear_twin_ID) - call results_writeDataset(group,stt%gamma_twin,'gamma_tw', & - 'twinning shear','1') - - end select - enddo outputsLoop - end associate - -end subroutine plastic_phenopowerlaw_results - - -!-------------------------------------------------------------------------------------------------- -!> @brief Shear rates on slip 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_slip(Mp,instance,of, & - gdot_slip_pos,gdot_slip_neg,dgdot_dtau_slip_pos,dgdot_dtau_slip_neg) - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - - real(pReal), intent(out), dimension(param(instance)%totalNslip) :: & - gdot_slip_pos, & - gdot_slip_neg - real(pReal), intent(out), optional, dimension(param(instance)%totalNslip) :: & - dgdot_dtau_slip_pos, & - dgdot_dtau_slip_neg - - real(pReal), dimension(param(instance)%totalNslip) :: & - 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%totalNslip - 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) - 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 - * sign(abs(tau_slip_pos/stt%xi_slip(:,of))**prm%n_slip, tau_slip_pos) - else where - gdot_slip_pos = 0.0_pReal - end where - - where(dNeq0(tau_slip_neg)) - gdot_slip_neg = prm%gdot0_slip * 0.5_pReal & ! only used if non-Schmid active, always 1/2 - * sign(abs(tau_slip_neg/stt%xi_slip(:,of))**prm%n_slip, tau_slip_neg) - else where - gdot_slip_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_slip/tau_slip_pos - else where - dgdot_dtau_slip_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_slip/tau_slip_neg - else where - dgdot_dtau_slip_neg = 0.0_pReal - end where - endif - end associate - -end subroutine kinetics_slip - - -!-------------------------------------------------------------------------------------------------- -!> @brief Shear rates on twin systems and their derivatives with respect to resolved stress. -! twinning is assumed to take place only in untwinned volume. -!> @details Derivates 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,instance,of,& - gdot_twin,dgdot_dtau_twin) - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - - real(pReal), dimension(param(instance)%totalNtwin), intent(out) :: & - gdot_twin - real(pReal), dimension(param(instance)%totalNtwin), intent(out), optional :: & - dgdot_dtau_twin - - real(pReal), dimension(param(instance)%totalNtwin) :: & - tau_twin - integer :: i - - associate(prm => param(instance), stt => state(instance)) - - do i = 1, prm%totalNtwin - tau_twin(i) = math_mul33xx33(Mp,prm%Schmid_twin(1:3,1:3,i)) - enddo - - where(tau_twin > 0.0_pReal) - gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,of)/prm%gamma_twin_char)) & ! only twin in untwinned volume fraction - * prm%gdot0_twin*(abs(tau_twin)/stt%xi_twin(:,of))**prm%n_twin - else where - gdot_twin = 0.0_pReal - end where - - if (present(dgdot_dtau_twin)) then - where(dNeq0(gdot_twin)) - dgdot_dtau_twin = gdot_twin*prm%n_twin/tau_twin - else where - dgdot_dtau_twin = 0.0_pReal - end where - endif - - end associate - -end subroutine kinetics_twin - -end module plastic_phenopowerlaw