From ced00954fe04ca74def0650ee0111df65d8e6405 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 15 Jan 2015 10:56:15 +0000 Subject: [PATCH] added new model by david (LLNL/UCLA) --- code/Makefile | 34 +- code/commercialFEM_fileList.f90 | 3 +- code/config/Phase_DisloUCLA_Tungsten.config | 48 + code/constitutive.f90 | 44 + code/material.f90 | 5 + ...stic_dislokmc.f90 => plastic_disloKMC.f90} | 0 code/plastic_disloUCLA.f90 | 1919 +++++++++++++++++ 7 files changed, 2038 insertions(+), 15 deletions(-) create mode 100644 code/config/Phase_DisloUCLA_Tungsten.config rename code/{plastic_dislokmc.f90 => plastic_disloKMC.f90} (100%) create mode 100644 code/plastic_disloUCLA.f90 diff --git a/code/Makefile b/code/Makefile index 9aad0ad8d..839e29cef 100644 --- a/code/Makefile +++ b/code/Makefile @@ -85,7 +85,7 @@ endif ifneq "$(FASTBUILD)" "YES" STANDARD_CHECK_ifort ?=-stand f08 -standard-semantics -STANDARD_CHECK_gfortran ?=-std=f2008 -pedantic-errors +STANDARD_CHECK_gfortran ?=-std=f2008ts -pedantic-errors endif #-std=f2008ts: for newer gfortran #-pedantic: more strict on standard, enables some warnings @@ -316,9 +316,10 @@ THERMAL_FILES = \ VACANCY_FILES = \ vacancy_constant.o vacancy_generation.o -CONSTITUTIVE_FILES = \ - plastic_dislotwin.o plastic_dislokmc.o plastic_j2.o plastic_phenopowerlaw.o \ +PLASTIC_FILES = \ + plastic_dislotwin.o plastic_disloKMC.o plastic_disloUCLA.o plastic_j2.o plastic_phenopowerlaw.o \ plastic_titanmod.o plastic_nonlocal.o plastic_none.o constitutive.o + HOMOGENIZATION_FILES = \ homogenization_RGC.o homogenization_isostrain.o homogenization_none.o homogenization.o @@ -335,7 +336,7 @@ DAMASK_spectral.exe: INTERFACENAME := DAMASK_spectral_interface.f90 SPECTRAL_FILES = prec.o DAMASK_interface.o IO.o libs.o numerics.o debug.o math.o \ FEsolving.o mesh.o material.o lattice.o \ - $(DAMAGE_FILES) $(THERMAL_FILES) $(VACANCY_FILES) $(CONSTITUTIVE_FILES) \ + $(DAMAGE_FILES) $(THERMAL_FILES) $(VACANCY_FILES) $(PLASTIC_FILES) \ crystallite.o $(HOMOGENIZATION_FILES) CPFEM.o \ DAMASK_spectral_utilities.o DAMASK_spectral_solverBasic.o \ DAMASK_spectral_solverAL.o DAMASK_spectral_solverBasicPETSc.o DAMASK_spectral_solverPolarisation.o @@ -344,7 +345,8 @@ DAMASK_spectral.exe: DAMASK_spectral_driver.o $(PREFIX) $(LINKERNAME) $(OPENMP_FLAG_$(F90)) $(LINK_OPTIONS_$(F90)) $(STANDARD_CHECK_$(F90)) $(OPTIMIZATION_$(MAXOPTI)_$(F90)) \ -o DAMASK_spectral.exe DAMASK_spectral_driver.o \ $(SPECTRAL_FILES) $(LIBRARIES) $(LIB_DIRS) $(RUN_PATH) $(SUFFIX) - + + DAMASK_spectral_driver.o: DAMASK_spectral_driver.f90 DAMASK_spectral_solverBasic.o \ DAMASK_spectral_solverAL.o \ DAMASK_spectral_solverBasicPETSc.o \ @@ -378,7 +380,7 @@ DAMASK_FEM.exe: INCLUDE_DIRS += -I./ FEM_FILES = prec.o DAMASK_interface.o FEZoo.o IO.o libs.o numerics.o debug.o math.o \ FEsolving.o mesh.o material.o lattice.o \ - $(DAMAGE_FILES) $(THERMAL_FILES) $(VACANCY_FILES) $(CONSTITUTIVE_FILES) \ + $(DAMAGE_FILES) $(THERMAL_FILES) $(VACANCY_FILES) $(PLASTIC_FILES) \ crystallite.o $(HOMOGENIZATION_FILES) CPFEM.o \ FEM_utilities.o FEM_mech.o FEM_thermal.o FEM_damage.o FEM_vacancyDiffusion.o @@ -434,7 +436,8 @@ constitutive.o: constitutive.f90 \ plastic_nonlocal.o \ plastic_titanmod.o \ plastic_dislotwin.o \ - plastic_dislokmc.o \ + plastic_disloKMC.o \ + plastic_disloUCLA.o \ plastic_phenopowerlaw.o \ plastic_j2.o \ plastic_none.o \ @@ -442,25 +445,28 @@ constitutive.o: constitutive.f90 \ $(THERMAL_FILES) \ $(VACANCY_FILES) -plastic_nonlocal.o: plastic_nonlocal.f90 \ +plastic_nonlocal.o: plastic_nonlocal.f90 \ lattice.o -plastic_titanmod.o: plastic_titanmod.f90 \ +plastic_titanmod.o: plastic_titanmod.f90 \ lattice.o -plastic_dislokmc.o: plastic_dislokmc.f90 \ +plastic_disloKMC.o: plastic_disloKMC.f90 \ lattice.o -plastic_dislotwin.o: plastic_dislotwin.f90 \ +plastic_disloUCLA.o: plastic_disloUCLA.f90 \ lattice.o -plastic_phenopowerlaw.o: plastic_phenopowerlaw.f90 \ +plastic_dislotwin.o: plastic_dislotwin.f90 \ lattice.o -plastic_j2.o: plastic_j2.f90 \ +plastic_phenopowerlaw.o: plastic_phenopowerlaw.f90 \ lattice.o -plastic_none.o: plastic_none.f90 \ +plastic_j2.o: plastic_j2.f90 \ + lattice.o + +plastic_none.o: plastic_none.f90 \ lattice.o damage_none.o: damage_none.f90 \ diff --git a/code/commercialFEM_fileList.f90 b/code/commercialFEM_fileList.f90 index c53de713b..608096309 100644 --- a/code/commercialFEM_fileList.f90 +++ b/code/commercialFEM_fileList.f90 @@ -30,7 +30,8 @@ #include "plastic_phenopowerlaw.f90" #include "plastic_titanmod.f90" #include "plastic_dislotwin.f90" -#include "plastic_dislokmc.f90" +#include "plastic_disloKMC.f90" +#include "plastic_disloUCLA.f90" #include "plastic_nonlocal.f90" #include "constitutive.f90" #include "crystallite.f90" diff --git a/code/config/Phase_DisloUCLA_Tungsten.config b/code/config/Phase_DisloUCLA_Tungsten.config new file mode 100644 index 000000000..57497eead --- /dev/null +++ b/code/config/Phase_DisloUCLA_Tungsten.config @@ -0,0 +1,48 @@ +[Tungsten] +elasticity hooke +#plasticity dislotwin +#plasticity dislokmc +plasticity disloucla + +### Material parameters ### +lattice_structure bcc +C11 523.0e9 # From Marinica et al. Journal of Physics: Condensed Matter(2013) +C12 202.0e9 +C44 161.0e9 + +grainsize 2.7e-5 # Average grain size [m] 2.0e-5 +SolidSolutionStrength 0.0 # Strength due to elements in solid solution + +### Dislocation glide parameters ### +#per family +Nslip 12 0 +slipburgers 2.72e-10 2.72e-10 # Burgers vector of slip system [m] +rhoedge0 1.0e12 1.0e12 # Initial edge dislocation density [m/m**3] 1.0e12 +rhoedgedip0 1.0 1.0 # Initial edged dipole dislocation density [m/m**3] +Qedge 2.725e-19 2.725e-19 # Activation energy for dislocation glide [J] +v0 3693.4 755.59 # Initial glide velocity [m/s] 1.0e-4 (kmC 3693.4 ; 755.59) +p_slip 0.86 0.22 # p-exponent in glide velocity +q_slip 1.69 1.01 # q-exponent in glide velocity +u_slip 2.47 0.38 # u-exponent of stress pre-factor (kmC) +s_slip 0.00 0.50 # self hardening in glide velocity (kmC) +tau_peierls 2.03e9 2.03e9 # peierls stress [Pa] +#nonschmid_a1 1.00 1.00 # a1 coeff. non schmid law +#nonschmid_a2 0.00 0.00 # a2 coeff. non schmid law + +#mobility law +kink_height 2.567e-10 2.567e-10 # kink height sqrt(6)/3*lattice_parameter [m] +omega 9.1e11 9.1e11 # attemp frequency (from kMC paper) [s^(-1)] +kink_width 29.95e-10 29.95e-10 # kink pair width ~ 11 b (kMC paper) [m] +dislolength 78.0e-10 78.0e-10 # dislocation length (ideally lambda) [m] initial value 25b +#friction_coeff 8.3e-5 8.3e-5 # friction coeff. B (kMC) [Pa*s] +# + +#hardening +dipoleformationfactor 0 0 # to have hardening due to dipole formation off +CLambdaSlip 10.0 10.0 # Adj. parameter controlling dislocation mean free path +D0 4.0e-5 4.0e-5 # Vacancy diffusion prefactor [m**2/s] +Qsd 4.5e-19 4.5e-19 # Activation energy for climb [J] +Catomicvolume 1.0 1.0 # Adj. parameter controlling the atomic volume [in b] +Cedgedipmindistance 1.0 1.0 # Adj. parameter controlling the minimum dipole distance [in b] +interaction_slipslip 0.0625 0.0625 0.72 0.05 0.09 0.06 # Queyreau +nonschmid_coefficients 0 0.56 0.75 0 0 0 #?????????????? diff --git a/code/constitutive.f90 b/code/constitutive.f90 index beb2c4605..e1d5232f5 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -117,6 +117,7 @@ subroutine constitutive_init(temperature_init) PLASTICITY_phenopowerlaw_ID, & PLASTICITY_dislotwin_ID, & PLASTICITY_dislokmc_ID, & + PLASTICITY_disloucla_ID, & PLASTICITY_titanmod_ID, & PLASTICITY_nonlocal_ID ,& ELASTICITY_HOOKE_label, & @@ -125,6 +126,7 @@ subroutine constitutive_init(temperature_init) PLASTICITY_PHENOPOWERLAW_label, & PLASTICITY_DISLOTWIN_label, & PLASTICITY_DISLOKMC_label, & + PLASTICITY_DISLOUCLA_label, & PLASTICITY_TITANMOD_label, & PLASTICITY_NONLOCAL_label, & LOCAL_DAMAGE_none_ID, & @@ -161,6 +163,7 @@ subroutine constitutive_init(temperature_init) use plastic_phenopowerlaw use plastic_dislotwin use plastic_dislokmc + use plastic_disloucla use plastic_titanmod use plastic_nonlocal use damage_none @@ -199,6 +202,7 @@ subroutine constitutive_init(temperature_init) if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_init(FILEUNIT) if (any(phase_plasticity == PLASTICITY_DISLOTWIN_ID)) call plastic_dislotwin_init(FILEUNIT) if (any(phase_plasticity == PLASTICITY_DISLOKMC_ID)) call plastic_dislokmc_init(FILEUNIT) + if (any(phase_plasticity == PLASTICITY_DISLOUCLA_ID)) call plastic_disloucla_init(FILEUNIT) if (any(phase_plasticity == PLASTICITY_TITANMOD_ID)) call plastic_titanmod_init(FILEUNIT) if (any(phase_plasticity == PLASTICITY_NONLOCAL_ID)) then call plastic_nonlocal_init(FILEUNIT) @@ -274,6 +278,11 @@ subroutine constitutive_init(temperature_init) thisNoutput => plastic_dislokmc_Noutput thisOutput => plastic_dislokmc_output thisSize => plastic_dislokmc_sizePostResult + case (PLASTICITY_DISLOUCLA_ID) + outputName = PLASTICITY_DISLOUCLA_label + thisNoutput => plastic_disloucla_Noutput + thisOutput => plastic_disloucla_output + thisSize => plastic_disloucla_sizePostResult case (PLASTICITY_TITANMOD_ID) outputName = PLASTICITY_TITANMOD_label thisNoutput => plastic_titanmod_Noutput @@ -477,6 +486,7 @@ function constitutive_homogenizedC(ipc,ip,el) PLASTICITY_TITANMOD_ID, & PLASTICITY_DISLOTWIN_ID, & PLASTICITY_DISLOKMC_ID, & + PLASTICITY_DISLOUCLA_ID, & plasticState,& mappingConstitutive @@ -486,6 +496,8 @@ function constitutive_homogenizedC(ipc,ip,el) plastic_dislotwin_homogenizedC use plastic_dislokmc, only: & plastic_dislokmc_homogenizedC + use plastic_disloucla, only: & + plastic_disloucla_homogenizedC use lattice, only: & lattice_C66 @@ -502,6 +514,8 @@ function constitutive_homogenizedC(ipc,ip,el) constitutive_homogenizedC = plastic_dislotwin_homogenizedC(ipc,ip,el) case (PLASTICITY_DISLOKMC_ID) constitutive_homogenizedC = plastic_dislokmc_homogenizedC(ipc,ip,el) + case (PLASTICITY_DISLOUCLA_ID) + constitutive_homogenizedC = plastic_disloucla_homogenizedC(ipc,ip,el) case (PLASTICITY_TITANMOD_ID) constitutive_homogenizedC = plastic_titanmod_homogenizedC (ipc,ip,el) case default @@ -569,6 +583,7 @@ subroutine constitutive_microstructure(Tstar_v, Fe, Fp, subdt, ipc, ip, el) material_phase, & PLASTICITY_dislotwin_ID, & PLASTICITY_dislokmc_ID, & + PLASTICITY_disloucla_ID, & PLASTICITY_titanmod_ID, & PLASTICITY_nonlocal_ID, & LOCAL_DAMAGE_isoBrittle_ID, & @@ -586,6 +601,8 @@ subroutine constitutive_microstructure(Tstar_v, Fe, Fp, subdt, ipc, ip, el) plastic_dislotwin_microstructure use plastic_dislokmc, only: & plastic_dislokmc_microstructure + use plastic_disloucla, only: & + plastic_disloucla_microstructure use damage_isoBrittle, only: & damage_isoBrittle_microstructure use damage_isoDuctile, only: & @@ -620,6 +637,8 @@ subroutine constitutive_microstructure(Tstar_v, Fe, Fp, subdt, ipc, ip, el) call plastic_dislotwin_microstructure(constitutive_getTemperature(ipc,ip,el),ipc,ip,el) case (PLASTICITY_DISLOKMC_ID) call plastic_dislokmc_microstructure(constitutive_getTemperature(ipc,ip,el),ipc,ip,el) + case (PLASTICITY_DISLOUCLA_ID) + call plastic_disloucla_microstructure(constitutive_getTemperature(ipc,ip,el),ipc,ip,el) case (PLASTICITY_TITANMOD_ID) call plastic_titanmod_microstructure (constitutive_getTemperature(ipc,ip,el),ipc,ip,el) case (PLASTICITY_NONLOCAL_ID) @@ -676,6 +695,7 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, ipc, ip, el) PLASTICITY_PHENOPOWERLAW_ID, & PLASTICITY_DISLOTWIN_ID, & PLASTICITY_DISLOKMC_ID, & + PLASTICITY_DISLOUCLA_ID, & PLASTICITY_TITANMOD_ID, & PLASTICITY_NONLOCAL_ID use plastic_j2, only: & @@ -689,6 +709,9 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, ipc, ip, el) use plastic_dislokmc, only: & plastic_dislokmc_LpAndItsTangent, & plastic_dislokmc_totalNslip + use plastic_disloucla, only: & + plastic_disloucla_LpAndItsTangent, & + plastic_disloucla_totalNslip use plastic_titanmod, only: & plastic_titanmod_LpAndItsTangent, & plastic_titanmod_totalNslip @@ -743,6 +766,12 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, ipc, ip, el) constitutive_getTemperature(ipc,ip,el), & nSlip,constitutive_getSlipDamage(nSlip,Tstar_v,ipc,ip,el), & ipc,ip,el) + case (PLASTICITY_DISLOUCLA_ID) + nSlip = plastic_disloucla_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))) + call plastic_disloucla_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, & + constitutive_getTemperature(ipc,ip,el), & + nSlip,constitutive_getSlipDamage(nSlip,Tstar_v,ipc,ip,el), & + ipc,ip,el) case (PLASTICITY_TITANMOD_ID) nSlip = plastic_titanmod_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))) call plastic_titanmod_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v, & @@ -1130,6 +1159,7 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su PLASTICITY_phenopowerlaw_ID, & PLASTICITY_dislotwin_ID, & PLASTICITY_dislokmc_ID, & + PLASTICITY_disloucla_ID, & PLASTICITY_titanmod_ID, & PLASTICITY_nonlocal_ID, & LOCAL_DAMAGE_gurson_ID @@ -1144,6 +1174,9 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su use plastic_dislokmc, only: & plastic_dislokmc_dotState, & plastic_dislokmc_totalNslip + use plastic_disloucla, only: & + plastic_disloucla_dotState, & + plastic_disloucla_totalNslip use plastic_titanmod, only: & plastic_titanmod_dotState, & plastic_titanmod_totalNslip @@ -1194,6 +1227,10 @@ subroutine constitutive_collectDotState(Tstar_v, Lp, FeArray, FpArray, subdt, su nSlip = plastic_dislokmc_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))) call plastic_dislokmc_dotState (Tstar_v,constitutive_getTemperature(ipc,ip,el), & nSlip,constitutive_getSlipDamage(nSlip,Tstar_v,ipc,ip,el),ipc,ip,el) + case (PLASTICITY_DISLOUCLA_ID) + nSlip = plastic_disloucla_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))) + call plastic_disloucla_dotState (Tstar_v,constitutive_getTemperature(ipc,ip,el), & + nSlip,constitutive_getSlipDamage(nSlip,Tstar_v,ipc,ip,el),ipc,ip,el) case (PLASTICITY_TITANMOD_ID) call plastic_titanmod_dotState (Tstar_v,constitutive_getTemperature(ipc,ip,el), & ipc,ip,el) @@ -1827,6 +1864,7 @@ function constitutive_getVacancyMobility33(ipc, ip, el) el !< element number real(pReal), dimension(3,3) :: & constitutive_getVacancyMobility33 + integer(pInt) :: & nSlip @@ -1895,6 +1933,7 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el) PLASTICITY_PHENOPOWERLAW_ID, & PLASTICITY_DISLOTWIN_ID, & PLASTICITY_DISLOKMC_ID, & + PLASTICITY_DISLOUCLA_ID, & PLASTICITY_TITANMOD_ID, & PLASTICITY_NONLOCAL_ID, & LOCAL_DAMAGE_isoBrittle_ID, & @@ -1916,6 +1955,8 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el) plastic_dislotwin_postResults use plastic_dislokmc, only: & plastic_dislokmc_postResults + use plastic_disloucla, only: & + plastic_disloucla_postResults use plastic_titanmod, only: & plastic_titanmod_postResults use plastic_nonlocal, only: & @@ -1979,6 +2020,9 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el) case (PLASTICITY_DISLOKMC_ID) constitutive_postResults(startPos:endPos) = & plastic_dislokmc_postResults(Tstar_v,constitutive_getTemperature(ipc,ip,el),ipc,ip,el) + case (PLASTICITY_DISLOUCLA_ID) + constitutive_postResults(startPos:endPos) = & + plastic_disloucla_postResults(Tstar_v,constitutive_getTemperature(ipc,ip,el),ipc,ip,el) case (PLASTICITY_NONLOCAL_ID) constitutive_postResults(startPos:endPos) = & plastic_nonlocal_postResults (Tstar_v,FeArray,ip,el) diff --git a/code/material.f90 b/code/material.f90 index 167f26dff..4ee6f9205 100644 --- a/code/material.f90 +++ b/code/material.f90 @@ -26,6 +26,7 @@ module material PLASTICITY_phenopowerlaw_label = 'phenopowerlaw', & PLASTICITY_dislotwin_label = 'dislotwin', & PLASTICITY_dislokmc_label = 'dislokmc', & + PLASTICITY_disloucla_label = 'disloucla', & PLASTICITY_titanmod_label = 'titanmod', & PLASTICITY_nonlocal_label = 'nonlocal', & LOCAL_DAMAGE_none_LABEL = 'none', & @@ -62,6 +63,7 @@ module material PLASTICITY_phenopowerlaw_ID, & PLASTICITY_dislotwin_ID, & PLASTICITY_dislokmc_ID, & + PLASTICITY_disloucla_ID, & PLASTICITY_titanmod_ID, & PLASTICITY_nonlocal_ID end enum @@ -240,6 +242,7 @@ module material PLASTICITY_phenopowerlaw_ID, & PLASTICITY_dislotwin_ID, & PLASTICITY_dislokmc_ID, & + PLASTICITY_disloucla_ID, & PLASTICITY_titanmod_ID, & PLASTICITY_nonlocal_ID, & LOCAL_DAMAGE_none_ID, & @@ -843,6 +846,8 @@ subroutine material_parsePhase(fileUnit,myPart) phase_plasticity(section) = PLASTICITY_DISLOTWIN_ID case (PLASTICITY_DISLOKMC_label) phase_plasticity(section) = PLASTICITY_DISLOKMC_ID + case (PLASTICITY_DISLOUCLA_label) + phase_plasticity(section) = PLASTICITY_DISLOUCLA_ID case (PLASTICITY_TITANMOD_label) phase_plasticity(section) = PLASTICITY_TITANMOD_ID case (PLASTICITY_NONLOCAL_label) diff --git a/code/plastic_dislokmc.f90 b/code/plastic_disloKMC.f90 similarity index 100% rename from code/plastic_dislokmc.f90 rename to code/plastic_disloKMC.f90 diff --git a/code/plastic_disloUCLA.f90 b/code/plastic_disloUCLA.f90 new file mode 100644 index 000000000..9acef87a6 --- /dev/null +++ b/code/plastic_disloUCLA.f90 @@ -0,0 +1,1919 @@ + !-------------------------------------------------------------------------------------------------- +! $Id$ +!-------------------------------------------------------------------------------------------------- +!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH +!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH +!> @author David Cereceda, Lawrence Livermore National Laboratory +!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH +!> @brief material subroutine incoprorating dislocation and twinning physics +!> @details to be done +!-------------------------------------------------------------------------------------------------- +module plastic_disloucla + use prec, only: & + pReal, & + pInt + + implicit none + private + integer(pInt), dimension(:), allocatable, public, protected :: & + plastic_disloucla_sizePostResults !< cumulative size of post results + + integer(pInt), dimension(:,:), allocatable, target, public :: & + plastic_disloucla_sizePostResult !< size of each post result output + + character(len=64), dimension(:,:), allocatable, target, public :: & + plastic_disloucla_output !< name of each post result output + + character(len=12), dimension(3), parameter, private :: & + plastic_disloucla_listBasicSlipStates = & + ['rhoEdge ', 'rhoEdgeDip ', 'accshearslip'] + + character(len=12), dimension(2), parameter, private :: & + plastic_disloucla_listBasicTwinStates = & + ['twinFraction', 'accsheartwin'] + + character(len=17), dimension(4), parameter, private :: & + plastic_disloucla_listDependentSlipStates = & + ['invLambdaSlip ', 'invLambdaSlipTwin', 'meanFreePathSlip ', 'tauSlipThreshold '] + + character(len=16), dimension(4), parameter, private :: & + plastic_disloucla_listDependentTwinStates = & + ['invLambdaTwin ', 'meanFreePathTwin', 'tauTwinThreshold', 'twinVolume '] + + real(pReal), parameter, private :: & + kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin + + integer(pInt), dimension(:), allocatable, target, public :: & + plastic_disloucla_Noutput !< number of outputs per instance of this plasticity + + integer(pInt), dimension(:), allocatable, public, protected :: & + plastic_disloucla_totalNslip, & !< total number of active slip systems for each instance + plastic_disloucla_totalNtwin !< total number of active twin systems for each instance + + integer(pInt), dimension(:,:), allocatable, private :: & + plastic_disloucla_Nslip, & !< number of active slip systems for each family and instance + plastic_disloucla_Ntwin !< number of active twin systems for each family and instance + + real(pReal), dimension(:), allocatable, private :: & + plastic_disloucla_CAtomicVolume, & !< atomic volume in Bugers vector unit + plastic_disloucla_D0, & !< prefactor for self-diffusion coefficient + plastic_disloucla_Qsd, & !< activation energy for dislocation climb + plastic_disloucla_GrainSize, & !< grain size + plastic_disloucla_MaxTwinFraction, & !< maximum allowed total twin volume fraction + plastic_disloucla_CEdgeDipMinDistance, & !< + plastic_disloucla_Cmfptwin, & !< + plastic_disloucla_Cthresholdtwin, & !< + plastic_disloucla_SolidSolutionStrength, & !< Strength due to elements in solid solution + plastic_disloucla_L0, & !< Length of twin nuclei in Burgers vectors + plastic_disloucla_xc, & !< critical distance for formation of twin nucleus + plastic_disloucla_VcrossSlip, & !< cross slip volume + plastic_disloucla_SFE_0K, & !< stacking fault energy at zero K + plastic_disloucla_dSFE_dT, & !< temperature dependance of stacking fault energy + plastic_disloucla_dipoleFormationFactor, & !< scaling factor for dipole formation: 0: off, 1: on. other values not useful + plastic_disloucla_aTolRho, & !< absolute tolerance for integration of dislocation density + plastic_disloucla_aTolTwinFrac !< absolute tolerance for integration of twin volume fraction + + real(pReal), dimension(:,:,:,:), allocatable, private :: & + plastic_disloucla_Ctwin66 !< twin elasticity matrix in Mandel notation for each instance + real(pReal), dimension(:,:,:,:,:,:), allocatable, private :: & + plastic_disloucla_Ctwin3333 !< twin elasticity matrix for each instance + real(pReal), dimension(:,:), allocatable, private :: & + plastic_disloucla_rhoEdge0, & !< initial edge dislocation density per slip system for each family and instance + plastic_disloucla_rhoEdgeDip0, & !< initial edge dipole density per slip system for each family and instance + plastic_disloucla_burgersPerSlipFamily, & !< absolute length of burgers vector [m] for each slip family and instance + plastic_disloucla_burgersPerSlipSystem, & !< absolute length of burgers vector [m] for each slip system and instance + plastic_disloucla_burgersPerTwinFamily, & !< absolute length of burgers vector [m] for each twin family and instance + plastic_disloucla_burgersPerTwinSystem, & !< absolute length of burgers vector [m] for each twin system and instance + plastic_disloucla_QedgePerSlipFamily, & !< activation energy for glide [J] for each slip family and instance + plastic_disloucla_QedgePerSlipSystem, & !< activation energy for glide [J] for each slip system and instance + plastic_disloucla_v0PerSlipFamily, & !< dislocation velocity prefactor [m/s] for each family and instance + plastic_disloucla_v0PerSlipSystem, & !< dislocation velocity prefactor [m/s] for each slip system and instance + plastic_disloucla_tau_peierlsPerSlipFamily, & !< Peierls stress [Pa] for each family and instance + plastic_disloucla_Ndot0PerTwinFamily, & !< twin nucleation rate [1/m³s] for each twin family and instance + plastic_disloucla_Ndot0PerTwinSystem, & !< twin nucleation rate [1/m³s] for each twin system and instance + plastic_disloucla_tau_r, & !< stress to bring partial close together for each twin system and instance + plastic_disloucla_twinsizePerTwinFamily, & !< twin thickness [m] for each twin family and instance + plastic_disloucla_twinsizePerTwinSystem, & !< twin thickness [m] for each twin system and instance + plastic_disloucla_CLambdaSlipPerSlipFamily, & !< Adj. parameter for distance between 2 forest dislocations for each slip family and instance + plastic_disloucla_CLambdaSlipPerSlipSystem, & !< Adj. parameter for distance between 2 forest dislocations for each slip system and instance + plastic_disloucla_interaction_SlipSlip, & !< coefficients for slip-slip interaction for each interaction type and instance + plastic_disloucla_interaction_SlipTwin, & !< coefficients for slip-twin interaction for each interaction type and instance + plastic_disloucla_interaction_TwinSlip, & !< coefficients for twin-slip interaction for each interaction type and instance + plastic_disloucla_interaction_TwinTwin, & !< coefficients for twin-twin interaction for each interaction type and instance + plastic_disloucla_pPerSlipFamily, & !< p-exponent in glide velocity + plastic_disloucla_qPerSlipFamily, & !< q-exponent in glide velocity + plastic_disloucla_uPerSlipFamily, & !< u-exponent in glide velocity + plastic_disloucla_sPerSlipFamily, & !< self-hardening in glide velocity + !* mobility law parameters + plastic_disloucla_kinkheight, & !< height of the kink pair + plastic_disloucla_omega, & !< attempt frequency for kink pair nucleation + plastic_disloucla_kinkwidth, & !< width of the kink pair + plastic_disloucla_dislolength, & !< dislocation length (lamda) + !* + plastic_disloucla_rPerTwinFamily, & !< r-exponent in twin nucleation rate + plastic_disloucla_nonSchmidCoeff !< non-Schmid coefficients (bcc) + real(pReal), dimension(:,:,:), allocatable, private :: & + plastic_disloucla_interactionMatrix_SlipSlip, & !< interaction matrix of the different slip systems for each instance + plastic_disloucla_interactionMatrix_SlipTwin, & !< interaction matrix of slip systems with twin systems for each instance + plastic_disloucla_interactionMatrix_TwinSlip, & !< interaction matrix of twin systems with slip systems for each instance + plastic_disloucla_interactionMatrix_TwinTwin, & !< interaction matrix of the different twin systems for each instance + plastic_disloucla_forestProjectionEdge !< matrix of forest projections of edge dislocations for each instance + + enum, bind(c) + enumerator :: undefined_ID, & + edge_density_ID, & + dipole_density_ID, & + shear_rate_slip_ID, & + accumulated_shear_slip_ID, & + mfp_slip_ID, & + resolved_stress_slip_ID, & + threshold_stress_slip_ID, & + edge_dipole_distance_ID, & + stress_exponent_ID, & + twin_fraction_ID, & + shear_rate_twin_ID, & + accumulated_shear_twin_ID, & + mfp_twin_ID, & + resolved_stress_twin_ID, & + threshold_stress_twin_ID + end enum + integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: & + plastic_disloucla_outputID !< ID of each post result output + + + public :: & + plastic_disloucla_init, & + plastic_disloucla_homogenizedC, & + plastic_disloucla_microstructure, & + plastic_disloucla_LpAndItsTangent, & + plastic_disloucla_dotState, & + plastic_disloucla_postResults + private :: & + plastic_disloucla_stateInit, & + plastic_disloucla_aTolState + +contains + + +!-------------------------------------------------------------------------------------------------- +!> @brief module initialization +!> @details reads in material parameters, allocates arrays, and does sanity checks +!-------------------------------------------------------------------------------------------------- +subroutine plastic_disloucla_init(fileUnit) + use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) + use debug, only: & + debug_level,& + debug_constitutive,& + debug_levelBasic + use math, only: & + math_Mandel3333to66, & + math_Voigt66to3333, & + math_mul3x3 + use mesh, only: & + mesh_NcpElems + use IO, only: & + IO_read, & + IO_lc, & + IO_getTag, & + IO_isBlank, & + IO_stringPos, & + IO_stringValue, & + IO_floatValue, & + IO_intValue, & + IO_warning, & + IO_error, & + IO_timeStamp, & + IO_EOF + use material, only: & + phase_plasticity, & + phase_plasticityInstance, & + phase_Noutput, & + PLASTICITY_DISLOUCLA_label, & + PLASTICITY_DISLOUCLA_ID, & + material_phase, & + plasticState, & + MATERIAL_partPhase + use lattice + use numerics,only: & + worldrank, & + numerics_integrator + + implicit none + integer(pInt), intent(in) :: fileUnit + + integer(pInt), parameter :: MAXNCHUNKS = LATTICE_maxNinteraction + 1_pInt + integer(pInt), dimension(1+2*MAXNCHUNKS) :: positions + integer(pInt) :: maxNinstance,mySize=0_pInt,phase,maxTotalNslip,maxTotalNtwin,& + f,instance,j,k,l,m,n,o,p,q,r,s,ns,nt, & + Nchunks_SlipSlip, Nchunks_SlipTwin, Nchunks_TwinSlip, Nchunks_TwinTwin, & + Nchunks_SlipFamilies, Nchunks_TwinFamilies, Nchunks_nonSchmid, & + offset_slip, index_myFamily, index_otherFamily + integer(pInt) :: sizeState, sizeDotState + integer(pInt) :: NofMyPhase + character(len=65536) :: & + tag = '', & + line = '' + real(pReal), dimension(:), allocatable :: tempPerSlip, tempPerTwin + + mainProcess: if (worldrank == 0) then + write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_DISLOUCLA_label//' init -+>>>' + write(6,'(a)') ' $Id$' + write(6,'(a15,a)') ' Current time: ',IO_timeStamp() +#include "compilation_info.f90" + endif mainProcess + + maxNinstance = int(count(phase_plasticity == PLASTICITY_DISLOUCLA_ID),pInt) + if (maxNinstance == 0_pInt) return + + if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & + write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance + + allocate(plastic_disloucla_sizePostResults(maxNinstance), source=0_pInt) + allocate(plastic_disloucla_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt) + allocate(plastic_disloucla_output(maxval(phase_Noutput),maxNinstance)) + plastic_disloucla_output = '' + allocate(plastic_disloucla_outputID(maxval(phase_Noutput),maxNinstance), source=undefined_ID) + allocate(plastic_disloucla_Noutput(maxNinstance), source=0_pInt) + allocate(plastic_disloucla_Nslip(lattice_maxNslipFamily,maxNinstance), source=0_pInt) + allocate(plastic_disloucla_Ntwin(lattice_maxNtwinFamily,maxNinstance), source=0_pInt) + allocate(plastic_disloucla_totalNslip(maxNinstance), source=0_pInt) + allocate(plastic_disloucla_totalNtwin(maxNinstance), source=0_pInt) + allocate(plastic_disloucla_CAtomicVolume(maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_D0(maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_Qsd(maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_GrainSize(maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_MaxTwinFraction(maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_CEdgeDipMinDistance(maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_Cmfptwin(maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_Cthresholdtwin(maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_SolidSolutionStrength(maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_L0(maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_xc(maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_VcrossSlip(maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_aTolRho(maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_aTolTwinFrac(maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_SFE_0K(maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_dSFE_dT(maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_dipoleFormationFactor(maxNinstance), source=1.0_pReal) !should be on by default + allocate(plastic_disloucla_rhoEdge0(lattice_maxNslipFamily,maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_rhoEdgeDip0(lattice_maxNslipFamily,maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_burgersPerSlipFamily(lattice_maxNslipFamily,maxNinstance),source=0.0_pReal) + allocate(plastic_disloucla_kinkheight(lattice_maxNslipFamily,maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_omega(lattice_maxNslipFamily,maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_kinkwidth(lattice_maxNslipFamily,maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_dislolength(lattice_maxNslipFamily,maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_burgersPerTwinFamily(lattice_maxNtwinFamily,maxNinstance),source=0.0_pReal) + allocate(plastic_disloucla_QedgePerSlipFamily(lattice_maxNslipFamily,maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_v0PerSlipFamily(lattice_maxNslipFamily,maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_tau_peierlsPerSlipFamily(lattice_maxNslipFamily,maxNinstance), & + source=0.0_pReal) + allocate(plastic_disloucla_pPerSlipFamily(lattice_maxNslipFamily,maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_qPerSlipFamily(lattice_maxNslipFamily,maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_uPerSlipFamily(lattice_maxNslipFamily,maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_sPerSlipFamily(lattice_maxNslipFamily,maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_Ndot0PerTwinFamily(lattice_maxNtwinFamily,maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_twinsizePerTwinFamily(lattice_maxNtwinFamily,maxNinstance),source=0.0_pReal) + allocate(plastic_disloucla_CLambdaSlipPerSlipFamily(lattice_maxNslipFamily,maxNinstance), & + source=0.0_pReal) + allocate(plastic_disloucla_rPerTwinFamily(lattice_maxNtwinFamily,maxNinstance),source=0.0_pReal) + allocate(plastic_disloucla_interaction_SlipSlip(lattice_maxNinteraction,maxNinstance),source=0.0_pReal) + allocate(plastic_disloucla_interaction_SlipTwin(lattice_maxNinteraction,maxNinstance),source=0.0_pReal) + allocate(plastic_disloucla_interaction_TwinSlip(lattice_maxNinteraction,maxNinstance),source=0.0_pReal) + allocate(plastic_disloucla_interaction_TwinTwin(lattice_maxNinteraction,maxNinstance),source=0.0_pReal) + allocate(plastic_disloucla_nonSchmidCoeff(lattice_maxNnonSchmid,maxNinstance), source=0.0_pReal) + + + rewind(fileUnit) + phase = 0_pInt + do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase) ! wind forward to + line = IO_read(fileUnit) + enddo + + parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part + line = IO_read(fileUnit) + if (IO_isBlank(line)) cycle ! skip empty lines + if (IO_getTag(line,'<','>') /= '') then ! stop at next part + line = IO_read(fileUnit, .true.) ! reset IO_read + exit + endif + if (IO_getTag(line,'[',']') /= '') then ! next phase section + phase = phase + 1_pInt ! advance phase section counter + if (phase_plasticity(phase) == PLASTICITY_DISLOUCLA_ID) then + Nchunks_SlipFamilies = count(lattice_NslipSystem(:,phase) > 0_pInt) + Nchunks_TwinFamilies = count(lattice_NtwinSystem(:,phase) > 0_pInt) + Nchunks_SlipSlip = maxval(lattice_interactionSlipSlip(:,:,phase)) + Nchunks_SlipTwin = maxval(lattice_interactionSlipTwin(:,:,phase)) + Nchunks_TwinSlip = maxval(lattice_interactionTwinSlip(:,:,phase)) + Nchunks_TwinTwin = maxval(lattice_interactionTwinTwin(:,:,phase)) + Nchunks_nonSchmid = lattice_NnonSchmid(phase) + if(allocated(tempPerSlip)) deallocate(tempPerSlip) + if(allocated(tempPerTwin)) deallocate(tempPerTwin) + allocate(tempPerSlip(Nchunks_SlipFamilies)) + allocate(tempPerTwin(Nchunks_TwinFamilies)) + endif + cycle ! skip to next line + endif + if (phase > 0_pInt ) then; if (phase_plasticity(phase) == PLASTICITY_DISLOUCLA_ID) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran + instance = phase_plasticityInstance(phase) ! which instance of my plasticity is present phase + positions = IO_stringPos(line,MAXNCHUNKS) + tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key + select case(tag) + case ('(output)') + select case(IO_lc(IO_stringValue(line,positions,2_pInt))) + case ('edge_density') + plastic_disloucla_Noutput(instance) = plastic_disloucla_Noutput(instance) + 1_pInt + plastic_disloucla_outputID(plastic_disloucla_Noutput(instance),instance) = edge_density_ID + plastic_disloucla_output(plastic_disloucla_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,positions,2_pInt)) + case ('dipole_density') + plastic_disloucla_Noutput(instance) = plastic_disloucla_Noutput(instance) + 1_pInt + plastic_disloucla_outputID(plastic_disloucla_Noutput(instance),instance) = dipole_density_ID + plastic_disloucla_output(plastic_disloucla_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,positions,2_pInt)) + case ('shear_rate_slip') + plastic_disloucla_Noutput(instance) = plastic_disloucla_Noutput(instance) + 1_pInt + plastic_disloucla_outputID(plastic_disloucla_Noutput(instance),instance) = shear_rate_slip_ID + plastic_disloucla_output(plastic_disloucla_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,positions,2_pInt)) + case ('accumulated_shear_slip') + plastic_disloucla_Noutput(instance) = plastic_disloucla_Noutput(instance) + 1_pInt + plastic_disloucla_outputID(plastic_disloucla_Noutput(instance),instance) = accumulated_shear_slip_ID + plastic_disloucla_output(plastic_disloucla_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,positions,2_pInt)) + case ('mfp_slip') + plastic_disloucla_Noutput(instance) = plastic_disloucla_Noutput(instance) + 1_pInt + plastic_disloucla_outputID(plastic_disloucla_Noutput(instance),instance) = mfp_slip_ID + plastic_disloucla_output(plastic_disloucla_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,positions,2_pInt)) + case ('resolved_stress_slip') + plastic_disloucla_Noutput(instance) = plastic_disloucla_Noutput(instance) + 1_pInt + plastic_disloucla_outputID(plastic_disloucla_Noutput(instance),instance) = resolved_stress_slip_ID + plastic_disloucla_output(plastic_disloucla_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,positions,2_pInt)) + case ('edge_dipole_distance') + plastic_disloucla_Noutput(instance) = plastic_disloucla_Noutput(instance) + 1_pInt + plastic_disloucla_outputID(plastic_disloucla_Noutput(instance),instance) = edge_dipole_distance_ID + plastic_disloucla_output(plastic_disloucla_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,positions,2_pInt)) + case ('stress_exponent') + plastic_disloucla_Noutput(instance) = plastic_disloucla_Noutput(instance) + 1_pInt + plastic_disloucla_outputID(plastic_disloucla_Noutput(instance),instance) = stress_exponent_ID + plastic_disloucla_output(plastic_disloucla_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,positions,2_pInt)) + case ('twin_fraction') + plastic_disloucla_Noutput(instance) = plastic_disloucla_Noutput(instance) + 1_pInt + plastic_disloucla_outputID(plastic_disloucla_Noutput(instance),instance) = twin_fraction_ID + plastic_disloucla_output(plastic_disloucla_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,positions,2_pInt)) + case ('shear_rate_twin') + plastic_disloucla_Noutput(instance) = plastic_disloucla_Noutput(instance) + 1_pInt + plastic_disloucla_outputID(plastic_disloucla_Noutput(instance),instance) = shear_rate_twin_ID + plastic_disloucla_output(plastic_disloucla_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,positions,2_pInt)) + case ('accumulated_shear_twin') + plastic_disloucla_Noutput(instance) = plastic_disloucla_Noutput(instance) + 1_pInt + plastic_disloucla_outputID(plastic_disloucla_Noutput(instance),instance) = accumulated_shear_twin_ID + plastic_disloucla_output(plastic_disloucla_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,positions,2_pInt)) + case ('mfp_twin') + plastic_disloucla_Noutput(instance) = plastic_disloucla_Noutput(instance) + 1_pInt + plastic_disloucla_outputID(plastic_disloucla_Noutput(instance),instance) = mfp_twin_ID + plastic_disloucla_output(plastic_disloucla_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,positions,2_pInt)) + case ('resolved_stress_twin') + plastic_disloucla_Noutput(instance) = plastic_disloucla_Noutput(instance) + 1_pInt + plastic_disloucla_outputID(plastic_disloucla_Noutput(instance),instance) = resolved_stress_twin_ID + plastic_disloucla_output(plastic_disloucla_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,positions,2_pInt)) + case ('threshold_stress_twin') + plastic_disloucla_Noutput(instance) = plastic_disloucla_Noutput(instance) + 1_pInt + plastic_disloucla_outputID(plastic_disloucla_Noutput(instance),instance) = threshold_stress_twin_ID + plastic_disloucla_output(plastic_disloucla_Noutput(instance),instance) = & + IO_lc(IO_stringValue(line,positions,2_pInt)) + end select +!-------------------------------------------------------------------------------------------------- +! parameters depending on number of slip system families + case ('nslip') + if (positions(1) < Nchunks_SlipFamilies + 1_pInt) & + call IO_warning(50_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_DISLOUCLA_label//')') + if (positions(1) > Nchunks_SlipFamilies + 1_pInt) & + call IO_error(150_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_DISLOUCLA_label//')') + Nchunks_SlipFamilies = positions(1) - 1_pInt + do j = 1_pInt, Nchunks_SlipFamilies + plastic_disloucla_Nslip(j,instance) = IO_intValue(line,positions,1_pInt+j) + enddo + case ('rhoedge0','rhoedgedip0','slipburgers','qedge','v0','clambdaslip','tau_peierls','p_slip','q_slip',& + 'u_slip','s_slip','kink_height','omega','kink_width','dislolength') + do j = 1_pInt, Nchunks_SlipFamilies + tempPerSlip(j) = IO_floatValue(line,positions,1_pInt+j) + enddo + select case(tag) + case ('rhoedge0') + plastic_disloucla_rhoEdge0(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies) + case ('rhoedgedip0') + plastic_disloucla_rhoEdgeDip0(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies) + case ('slipburgers') + plastic_disloucla_burgersPerSlipFamily(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies) + case ('qedge') + plastic_disloucla_QedgePerSlipFamily(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies) + case ('v0') + plastic_disloucla_v0PerSlipFamily(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies) + case ('clambdaslip') + plastic_disloucla_CLambdaSlipPerSlipFamily(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies) + case ('tau_peierls') + if (lattice_structure(phase) /= LATTICE_bcc_ID) & + call IO_warning(42_pInt,ext_msg=trim(tag)//' for non-bcc ('//PLASTICITY_DISLOUCLA_label//')') + plastic_disloucla_tau_peierlsPerSlipFamily(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies) + case ('p_slip') + plastic_disloucla_pPerSlipFamily(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies) + case ('q_slip') + plastic_disloucla_qPerSlipFamily(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies) + case ('u_slip') + plastic_disloucla_uPerSlipFamily(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies) + case ('s_slip') + plastic_disloucla_sPerSlipFamily(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies) + case ('kink_height') + plastic_disloucla_kinkheight(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies) + case ('omega') + plastic_disloucla_omega(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies) + case ('kink_width') + plastic_disloucla_kinkwidth(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies) + case ('dislolength') + plastic_disloucla_dislolength(1:Nchunks_SlipFamilies,instance) = tempPerSlip(1:Nchunks_SlipFamilies) + end select +!-------------------------------------------------------------------------------------------------- +! parameters depending on slip number of twin families + case ('ntwin') + if (positions(1) < Nchunks_TwinFamilies + 1_pInt) & + call IO_warning(51_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_DISLOUCLA_label//')') + if (positions(1) > Nchunks_TwinFamilies + 1_pInt) & + call IO_error(150_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_DISLOUCLA_label//')') + Nchunks_TwinFamilies = positions(1) - 1_pInt + do j = 1_pInt, Nchunks_TwinFamilies + plastic_disloucla_Ntwin(j,instance) = IO_intValue(line,positions,1_pInt+j) + enddo + case ('ndot0','twinsize','twinburgers','r_twin') + do j = 1_pInt, Nchunks_TwinFamilies + tempPerTwin(j) = IO_floatValue(line,positions,1_pInt+j) + enddo + select case(tag) + case ('ndot0') + if (lattice_structure(phase) == LATTICE_fcc_ID) & + call IO_warning(42_pInt,ext_msg=trim(tag)//' for fcc ('//PLASTICITY_DISLOUCLA_label//')') + plastic_disloucla_Ndot0PerTwinFamily(1:Nchunks_TwinFamilies,instance) = tempPerTwin(1:Nchunks_TwinFamilies) + case ('twinsize') + plastic_disloucla_twinsizePerTwinFamily(1:Nchunks_TwinFamilies,instance) = tempPerTwin(1:Nchunks_TwinFamilies) + case ('twinburgers') + plastic_disloucla_burgersPerTwinFamily(1:Nchunks_TwinFamilies,instance) = tempPerTwin(1:Nchunks_TwinFamilies) + case ('r_twin') + plastic_disloucla_rPerTwinFamily(1:Nchunks_TwinFamilies,instance) = tempPerTwin(1:Nchunks_TwinFamilies) + end select +!-------------------------------------------------------------------------------------------------- +! parameters depending on number of interactions + case ('interaction_slipslip','interactionslipslip') + if (positions(1) < 1_pInt + Nchunks_SlipSlip) & + call IO_warning(52_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_DISLOUCLA_label//')') + do j = 1_pInt, Nchunks_SlipSlip + plastic_disloucla_interaction_SlipSlip(j,instance) = IO_floatValue(line,positions,1_pInt+j) + enddo + case ('interaction_sliptwin','interactionsliptwin') + if (positions(1) < 1_pInt + Nchunks_SlipTwin) & + call IO_warning(52_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_DISLOUCLA_label//')') + do j = 1_pInt, Nchunks_SlipTwin + plastic_disloucla_interaction_SlipTwin(j,instance) = IO_floatValue(line,positions,1_pInt+j) + enddo + case ('interaction_twinslip','interactiontwinslip') + if (positions(1) < 1_pInt + Nchunks_TwinSlip) & + call IO_warning(52_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_DISLOUCLA_label//')') + do j = 1_pInt, Nchunks_TwinSlip + plastic_disloucla_interaction_TwinSlip(j,instance) = IO_floatValue(line,positions,1_pInt+j) + enddo + case ('interaction_twintwin','interactiontwintwin') + if (positions(1) < 1_pInt + Nchunks_TwinTwin) & + call IO_warning(52_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_DISLOUCLA_label//')') + do j = 1_pInt, Nchunks_TwinTwin + plastic_disloucla_interaction_TwinTwin(j,instance) = IO_floatValue(line,positions,1_pInt+j) + enddo + case ('nonschmid_coefficients') + if (positions(1) < 1_pInt + Nchunks_nonSchmid) & + call IO_warning(52_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_DISLOUCLA_label//')') + do j = 1_pInt,Nchunks_nonSchmid + plastic_disloucla_nonSchmidCoeff(j,instance) = IO_floatValue(line,positions,1_pInt+j) + enddo +!-------------------------------------------------------------------------------------------------- +! parameters independent of number of slip/twin systems + case ('grainsize') + plastic_disloucla_GrainSize(instance) = IO_floatValue(line,positions,2_pInt) + case ('maxtwinfraction') + plastic_disloucla_MaxTwinFraction(instance) = IO_floatValue(line,positions,2_pInt) + case ('d0') + plastic_disloucla_D0(instance) = IO_floatValue(line,positions,2_pInt) + case ('qsd') + plastic_disloucla_Qsd(instance) = IO_floatValue(line,positions,2_pInt) + case ('atol_rho') + plastic_disloucla_aTolRho(instance) = IO_floatValue(line,positions,2_pInt) + case ('atol_twinfrac') + plastic_disloucla_aTolTwinFrac(instance) = IO_floatValue(line,positions,2_pInt) + case ('cmfptwin') + plastic_disloucla_Cmfptwin(instance) = IO_floatValue(line,positions,2_pInt) + case ('cthresholdtwin') + plastic_disloucla_Cthresholdtwin(instance) = IO_floatValue(line,positions,2_pInt) + case ('solidsolutionstrength') + plastic_disloucla_SolidSolutionStrength(instance) = IO_floatValue(line,positions,2_pInt) + case ('l0') + plastic_disloucla_L0(instance) = IO_floatValue(line,positions,2_pInt) + case ('xc') + plastic_disloucla_xc(instance) = IO_floatValue(line,positions,2_pInt) + case ('vcrossslip') + plastic_disloucla_VcrossSlip(instance) = IO_floatValue(line,positions,2_pInt) + case ('cedgedipmindistance') + plastic_disloucla_CEdgeDipMinDistance(instance) = IO_floatValue(line,positions,2_pInt) + case ('catomicvolume') + plastic_disloucla_CAtomicVolume(instance) = IO_floatValue(line,positions,2_pInt) + case ('sfe_0k') + plastic_disloucla_SFE_0K(instance) = IO_floatValue(line,positions,2_pInt) + case ('dsfe_dt') + plastic_disloucla_dSFE_dT(instance) = IO_floatValue(line,positions,2_pInt) + case ('dipoleformationfactor') + plastic_disloucla_dipoleFormationFactor(instance) = IO_floatValue(line,positions,2_pInt) + end select + endif; endif + enddo parsingFile + + sanityChecks: do phase = 1_pInt, size(phase_plasticity) + myPhase: if (phase_plasticity(phase) == PLASTICITY_disloucla_ID) then + instance = phase_plasticityInstance(phase) + if (sum(plastic_disloucla_Nslip(:,instance)) < 0_pInt) & + call IO_error(211_pInt,el=instance,ext_msg='Nslip ('//PLASTICITY_DISLOUCLA_label//')') + if (sum(plastic_disloucla_Ntwin(:,instance)) < 0_pInt) & + call IO_error(211_pInt,el=instance,ext_msg='Ntwin ('//PLASTICITY_DISLOUCLA_label//')') + do f = 1_pInt,lattice_maxNslipFamily + if (plastic_disloucla_Nslip(f,instance) > 0_pInt) then + if (plastic_disloucla_rhoEdge0(f,instance) < 0.0_pReal) & + call IO_error(211_pInt,el=instance,ext_msg='rhoEdge0 ('//PLASTICITY_DISLOUCLA_label//')') + if (plastic_disloucla_rhoEdgeDip0(f,instance) < 0.0_pReal) & + call IO_error(211_pInt,el=instance,ext_msg='rhoEdgeDip0 ('//PLASTICITY_DISLOUCLA_label//')') + if (plastic_disloucla_burgersPerSlipFamily(f,instance) <= 0.0_pReal) & + call IO_error(211_pInt,el=instance,ext_msg='slipBurgers ('//PLASTICITY_DISLOUCLA_label//')') + if (plastic_disloucla_v0PerSlipFamily(f,instance) <= 0.0_pReal) & + call IO_error(211_pInt,el=instance,ext_msg='v0 ('//PLASTICITY_DISLOUCLA_label//')') + if (plastic_disloucla_tau_peierlsPerSlipFamily(f,instance) < 0.0_pReal) & + call IO_error(211_pInt,el=instance,ext_msg='tau_peierls ('//PLASTICITY_DISLOUCLA_label//')') + endif + enddo + do f = 1_pInt,lattice_maxNtwinFamily + if (plastic_disloucla_Ntwin(f,instance) > 0_pInt) then + if (plastic_disloucla_burgersPerTwinFamily(f,instance) <= 0.0_pReal) & + call IO_error(211_pInt,el=instance,ext_msg='twinburgers ('//PLASTICITY_DISLOUCLA_label//')') + if (plastic_disloucla_Ndot0PerTwinFamily(f,instance) < 0.0_pReal) & + call IO_error(211_pInt,el=instance,ext_msg='ndot0 ('//PLASTICITY_DISLOUCLA_label//')') + endif + enddo + if (plastic_disloucla_CAtomicVolume(instance) <= 0.0_pReal) & + call IO_error(211_pInt,el=instance,ext_msg='cAtomicVolume ('//PLASTICITY_DISLOUCLA_label//')') + if (plastic_disloucla_D0(instance) <= 0.0_pReal) & + call IO_error(211_pInt,el=instance,ext_msg='D0 ('//PLASTICITY_DISLOUCLA_label//')') + if (plastic_disloucla_Qsd(instance) <= 0.0_pReal) & + call IO_error(211_pInt,el=instance,ext_msg='Qsd ('//PLASTICITY_DISLOUCLA_label//')') + if (sum(plastic_disloucla_Ntwin(:,instance)) > 0_pInt) then + if (plastic_disloucla_SFE_0K(instance) == 0.0_pReal .and. & + plastic_disloucla_dSFE_dT(instance) == 0.0_pReal .and. & + lattice_structure(phase) == LATTICE_fcc_ID) & + call IO_error(211_pInt,el=instance,ext_msg='SFE0K ('//PLASTICITY_DISLOUCLA_label//')') + if (plastic_disloucla_aTolRho(instance) <= 0.0_pReal) & + call IO_error(211_pInt,el=instance,ext_msg='aTolRho ('//PLASTICITY_DISLOUCLA_label//')') + if (plastic_disloucla_aTolTwinFrac(instance) <= 0.0_pReal) & + call IO_error(211_pInt,el=instance,ext_msg='aTolTwinFrac ('//PLASTICITY_DISLOUCLA_label//')') + endif + if (plastic_disloucla_dipoleFormationFactor(instance) /= 0.0_pReal .and. & + plastic_disloucla_dipoleFormationFactor(instance) /= 1.0_pReal) & + call IO_error(211_pInt,el=instance,ext_msg='dipoleFormationFactor ('//PLASTICITY_DISLOUCLA_label//')') + +!-------------------------------------------------------------------------------------------------- +! Determine total number of active slip or twin systems + plastic_disloucla_Nslip(:,instance) = min(lattice_NslipSystem(:,phase),plastic_disloucla_Nslip(:,instance)) + plastic_disloucla_Ntwin(:,instance) = min(lattice_NtwinSystem(:,phase),plastic_disloucla_Ntwin(:,instance)) + plastic_disloucla_totalNslip(instance) = sum(plastic_disloucla_Nslip(:,instance)) + plastic_disloucla_totalNtwin(instance) = sum(plastic_disloucla_Ntwin(:,instance)) + endif myPhase + enddo sanityChecks + +!-------------------------------------------------------------------------------------------------- +! allocation of variables whose size depends on the total number of active slip systems + maxTotalNslip = maxval(plastic_disloucla_totalNslip) + maxTotalNtwin = maxval(plastic_disloucla_totalNtwin) + + allocate(plastic_disloucla_burgersPerSlipSystem(maxTotalNslip, maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_burgersPerTwinSystem(maxTotalNtwin, maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_QedgePerSlipSystem(maxTotalNslip, maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_v0PerSlipSystem(maxTotalNslip, maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_Ndot0PerTwinSystem(maxTotalNtwin, maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_tau_r(maxTotalNtwin, maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_twinsizePerTwinSystem(maxTotalNtwin, maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_CLambdaSlipPerSlipSystem(maxTotalNslip, maxNinstance),source=0.0_pReal) + + allocate(plastic_disloucla_interactionMatrix_SlipSlip(maxval(plastic_disloucla_totalNslip),& ! slip resistance from slip activity + maxval(plastic_disloucla_totalNslip),& + maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_interactionMatrix_SlipTwin(maxval(plastic_disloucla_totalNslip),& ! slip resistance from twin activity + maxval(plastic_disloucla_totalNtwin),& + maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_interactionMatrix_TwinSlip(maxval(plastic_disloucla_totalNtwin),& ! twin resistance from slip activity + maxval(plastic_disloucla_totalNslip),& + maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_interactionMatrix_TwinTwin(maxval(plastic_disloucla_totalNtwin),& ! twin resistance from twin activity + maxval(plastic_disloucla_totalNtwin),& + maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_forestProjectionEdge(maxTotalNslip,maxTotalNslip,maxNinstance), & + source=0.0_pReal) + allocate(plastic_disloucla_Ctwin66(6,6,maxTotalNtwin,maxNinstance), source=0.0_pReal) + allocate(plastic_disloucla_Ctwin3333(3,3,3,3,maxTotalNtwin,maxNinstance), source=0.0_pReal) + + initializeInstances: do phase = 1_pInt, size(phase_plasticity) + myPhase2: if (phase_plasticity(phase) == PLASTICITY_disloucla_ID) then + NofMyPhase=count(material_phase==phase) + instance = phase_plasticityInstance(phase) + + ns = plastic_disloucla_totalNslip(instance) + nt = plastic_disloucla_totalNtwin(instance) + +!-------------------------------------------------------------------------------------------------- +! Determine size of postResults array + outputs: do o = 1_pInt,plastic_disloucla_Noutput(instance) + select case(plastic_disloucla_outputID(o,instance)) + case(edge_density_ID, & + dipole_density_ID, & + shear_rate_slip_ID, & + accumulated_shear_slip_ID, & + mfp_slip_ID, & + resolved_stress_slip_ID, & + threshold_stress_slip_ID, & + edge_dipole_distance_ID, & + stress_exponent_ID & + ) + mySize = ns + case(twin_fraction_ID, & + shear_rate_twin_ID, & + accumulated_shear_twin_ID, & + mfp_twin_ID, & + resolved_stress_twin_ID, & + threshold_stress_twin_ID & + ) + mySize = nt + end select + + if (mySize > 0_pInt) then ! any meaningful output found + plastic_disloucla_sizePostResult(o,instance) = mySize + plastic_disloucla_sizePostResults(instance) = plastic_disloucla_sizePostResults(instance) + mySize + endif + enddo outputs + +!-------------------------------------------------------------------------------------------------- +! allocate state arrays + sizeDotState = int(size(plastic_disloucla_listBasicSlipStates),pInt) * ns & + + int(size(plastic_disloucla_listBasicTwinStates),pInt) * nt + sizeState = sizeDotState & + + int(size(plastic_disloucla_listDependentSlipStates),pInt) * ns & + + int(size(plastic_disloucla_listDependentTwinStates),pInt) * nt + + plasticState(phase)%sizeState = sizeState + plasticState(phase)%sizeDotState = sizeDotState + plasticState(phase)%nSlip = plastic_disloucla_totalNslip(instance) + plasticState(phase)%nTwin = 0_pInt + plasticState(phase)%nTrans= 0_pInt + allocate(plasticState(phase)%aTolState (sizeState), source=0.0_pReal) + allocate(plasticState(phase)%state0 (sizeState,NofMyPhase), source=0.0_pReal) + allocate(plasticState(phase)%partionedState0 (sizeState,NofMyPhase), source=0.0_pReal) + allocate(plasticState(phase)%subState0 (sizeState,NofMyPhase), source=0.0_pReal) + allocate(plasticState(phase)%state (sizeState,NofMyPhase), source=0.0_pReal) + allocate(plasticState(phase)%state_backup (sizeState,NofMyPhase), source=0.0_pReal) + + allocate(plasticState(phase)%dotState (sizeDotState,NofMyPhase), source=0.0_pReal) + allocate(plasticState(phase)%deltaState (sizeDotState,NofMyPhase), source=0.0_pReal) + allocate(plasticState(phase)%dotState_backup (sizeDotState,NofMyPhase), source=0.0_pReal) + if (any(numerics_integrator == 1_pInt)) then + allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase), source=0.0_pReal) + allocate(plasticState(phase)%previousDotState2 (sizeDotState,NofMyPhase), source=0.0_pReal) + endif + if (any(numerics_integrator == 4_pInt)) & + allocate(plasticState(phase)%RK4dotState (sizeDotState,NofMyPhase), source=0.0_pReal) + if (any(numerics_integrator == 5_pInt)) & + allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal) + offset_slip = 2_pInt*plasticState(phase)%nSlip + plasticState(phase)%slipRate => & + plasticState(phase)%dotState(offset_slip+1:offset_slip+plasticState(phase)%nSlip,1:NofMyPhase) + plasticState(phase)%accumulatedSlip => & + plasticState(phase)%state (offset_slip+1:offset_slip+plasticState(phase)%nSlip,1:NofMyPhase) + !* Process slip related parameters ------------------------------------------------ + + mySlipFamilies: do f = 1_pInt,lattice_maxNslipFamily + index_myFamily = sum(plastic_disloucla_Nslip(1:f-1_pInt,instance)) ! index in truncated slip system list + mySlipSystems: do j = 1_pInt,plastic_disloucla_Nslip(f,instance) + + !* Burgers vector, + ! dislocation velocity prefactor, + ! mean free path prefactor, + ! and minimum dipole distance + + plastic_disloucla_burgersPerSlipSystem(index_myFamily+j,instance) = & + plastic_disloucla_burgersPerSlipFamily(f,instance) + + plastic_disloucla_QedgePerSlipSystem(index_myFamily+j,instance) = & + plastic_disloucla_QedgePerSlipFamily(f,instance) + + plastic_disloucla_v0PerSlipSystem(index_myFamily+j,instance) = & + plastic_disloucla_v0PerSlipFamily(f,instance) + + plastic_disloucla_CLambdaSlipPerSlipSystem(index_myFamily+j,instance) = & + plastic_disloucla_CLambdaSlipPerSlipFamily(f,instance) + + !* Calculation of forest projections for edge dislocations + !* Interaction matrices + + otherSlipFamilies: do o = 1_pInt,lattice_maxNslipFamily + index_otherFamily = sum(plastic_disloucla_Nslip(1:o-1_pInt,instance)) + otherSlipSystems: do k = 1_pInt,plastic_disloucla_Nslip(o,instance) + plastic_disloucla_forestProjectionEdge(index_myFamily+j,index_otherFamily+k,instance) = & + abs(math_mul3x3(lattice_sn(:,sum(lattice_NslipSystem(1:f-1,phase))+j,phase), & + lattice_st(:,sum(lattice_NslipSystem(1:o-1,phase))+k,phase))) + plastic_disloucla_interactionMatrix_SlipSlip(index_myFamily+j,index_otherFamily+k,instance) = & + plastic_disloucla_interaction_SlipSlip(lattice_interactionSlipSlip( & + sum(lattice_NslipSystem(1:f-1,phase))+j, & + sum(lattice_NslipSystem(1:o-1,phase))+k, & + phase), instance ) + enddo otherSlipSystems; enddo otherSlipFamilies + + otherTwinFamilies: do o = 1_pInt,lattice_maxNtwinFamily + index_otherFamily = sum(plastic_disloucla_Ntwin(1:o-1_pInt,instance)) + otherTwinSystems: do k = 1_pInt,plastic_disloucla_Ntwin(o,instance) + plastic_disloucla_interactionMatrix_SlipTwin(index_myFamily+j,index_otherFamily+k,instance) = & + plastic_disloucla_interaction_SlipTwin(lattice_interactionSlipTwin( & + sum(lattice_NslipSystem(1:f-1_pInt,phase))+j, & + sum(lattice_NtwinSystem(1:o-1_pInt,phase))+k, & + phase), instance ) + enddo otherTwinSystems; enddo otherTwinFamilies + + enddo mySlipSystems + enddo mySlipFamilies + + !* Process twin related parameters ------------------------------------------------ + + myTwinFamilies: do f = 1_pInt,lattice_maxNtwinFamily + index_myFamily = sum(plastic_disloucla_Ntwin(1:f-1_pInt,instance)) ! index in truncated twin system list + myTwinSystems: do j = 1_pInt,plastic_disloucla_Ntwin(f,instance) + + !* Burgers vector, + ! nucleation rate prefactor, + ! and twin size + + plastic_disloucla_burgersPerTwinSystem(index_myFamily+j,instance) = & + plastic_disloucla_burgersPerTwinFamily(f,instance) + + plastic_disloucla_Ndot0PerTwinSystem(index_myFamily+j,instance) = & + plastic_disloucla_Ndot0PerTwinFamily(f,instance) + + plastic_disloucla_twinsizePerTwinSystem(index_myFamily+j,instance) = & + plastic_disloucla_twinsizePerTwinFamily(f,instance) + + !* Rotate twin elasticity matrices + + index_otherFamily = sum(lattice_NtwinSystem(1:f-1_pInt,phase)) ! index in full lattice twin list + do l = 1_pInt,3_pInt; do m = 1_pInt,3_pInt; do n = 1_pInt,3_pInt; do o = 1_pInt,3_pInt + do p = 1_pInt,3_pInt; do q = 1_pInt,3_pInt; do r = 1_pInt,3_pInt; do s = 1_pInt,3_pInt + plastic_disloucla_Ctwin3333(l,m,n,o,index_myFamily+j,instance) = & + plastic_disloucla_Ctwin3333(l,m,n,o,index_myFamily+j,instance) + & + lattice_C3333(p,q,r,s,instance) * & + lattice_Qtwin(l,p,index_otherFamily+j,phase) * & + lattice_Qtwin(m,q,index_otherFamily+j,phase) * & + lattice_Qtwin(n,r,index_otherFamily+j,phase) * & + lattice_Qtwin(o,s,index_otherFamily+j,phase) + enddo; enddo; enddo; enddo + enddo; enddo; enddo; enddo + plastic_disloucla_Ctwin66(1:6,1:6,index_myFamily+j,instance) = & + math_Mandel3333to66(plastic_disloucla_Ctwin3333(1:3,1:3,1:3,1:3,index_myFamily+j,instance)) + + !* Interaction matrices + otherSlipFamilies2: do o = 1_pInt,lattice_maxNslipFamily + index_otherFamily = sum(plastic_disloucla_Nslip(1:o-1_pInt,instance)) + otherSlipSystems2: do k = 1_pInt,plastic_disloucla_Nslip(o,instance) + plastic_disloucla_interactionMatrix_TwinSlip(index_myFamily+j,index_otherFamily+k,instance) = & + plastic_disloucla_interaction_TwinSlip(lattice_interactionTwinSlip( & + sum(lattice_NtwinSystem(1:f-1_pInt,phase))+j, & + sum(lattice_NslipSystem(1:o-1_pInt,phase))+k, & + phase), instance ) + enddo otherSlipSystems2; enddo otherSlipFamilies2 + + otherTwinFamilies2: do o = 1_pInt,lattice_maxNtwinFamily + index_otherFamily = sum(plastic_disloucla_Ntwin(1:o-1_pInt,instance)) + otherTwinSystems2: do k = 1_pInt,plastic_disloucla_Ntwin(o,instance) + plastic_disloucla_interactionMatrix_TwinTwin(index_myFamily+j,index_otherFamily+k,instance) = & + plastic_disloucla_interaction_TwinTwin(lattice_interactionTwinTwin( & + sum(lattice_NtwinSystem(1:f-1_pInt,phase))+j, & + sum(lattice_NtwinSystem(1:o-1_pInt,phase))+k, & + phase), instance ) + enddo otherTwinSystems2; enddo otherTwinFamilies2 + + enddo myTwinSystems + enddo myTwinFamilies + call plastic_disloucla_stateInit(phase,instance) + call plastic_disloucla_aTolState(phase,instance) + endif myPhase2 + + enddo initializeInstances + +end subroutine plastic_disloucla_init + +!-------------------------------------------------------------------------------------------------- +!> @brief sets the relevant state values for a given instance of this plasticity +!-------------------------------------------------------------------------------------------------- +subroutine plastic_disloucla_stateInit(ph,instance) + use math, only: & + pi + use lattice, only: & + lattice_maxNslipFamily, & + lattice_mu + use material, only: & + plasticState + + implicit none + integer(pInt), intent(in) :: & + instance, & !< number specifying the instance of the plasticity + ph + + real(pReal), dimension(plasticState(ph)%sizeState) :: tempState + + integer(pInt) :: i,j,f,ns,nt, index_myFamily + real(pReal), dimension(plastic_disloucla_totalNslip(instance)) :: & + rhoEdge0, & + rhoEdgeDip0, & + invLambdaSlip0, & + MeanFreePathSlip0, & + tauSlipThreshold0 + real(pReal), dimension(plastic_disloucla_totalNtwin(instance)) :: & + MeanFreePathTwin0,TwinVolume0 + tempState = 0.0_pReal + ns = plastic_disloucla_totalNslip(instance) + nt = plastic_disloucla_totalNtwin(instance) + +!-------------------------------------------------------------------------------------------------- +! initialize basic slip state variables + do f = 1_pInt,lattice_maxNslipFamily + index_myFamily = sum(plastic_disloucla_Nslip(1:f-1_pInt,instance)) ! index in truncated slip system list + rhoEdge0(index_myFamily+1_pInt: & + index_myFamily+plastic_disloucla_Nslip(f,instance)) = & + plastic_disloucla_rhoEdge0(f,instance) + rhoEdgeDip0(index_myFamily+1_pInt: & + index_myFamily+plastic_disloucla_Nslip(f,instance)) = & + plastic_disloucla_rhoEdgeDip0(f,instance) + enddo + + tempState(1_pInt:ns) = rhoEdge0 + tempState(ns+1_pInt:2_pInt*ns) = rhoEdgeDip0 + +!-------------------------------------------------------------------------------------------------- +! initialize dependent slip microstructural variables + forall (i = 1_pInt:ns) & + invLambdaSlip0(i) = sqrt(dot_product((rhoEdge0+rhoEdgeDip0),plastic_disloucla_forestProjectionEdge(1:ns,i,instance)))/ & + plastic_disloucla_CLambdaSlipPerSlipSystem(i,instance) + tempState(3_pInt*ns+2_pInt*nt+1:4_pInt*ns+2_pInt*nt) = invLambdaSlip0 + + forall (i = 1_pInt:ns) & + MeanFreePathSlip0(i) = & + plastic_disloucla_GrainSize(instance)/(1.0_pReal+invLambdaSlip0(i)*plastic_disloucla_GrainSize(instance)) + tempState(5_pInt*ns+3_pInt*nt+1:6_pInt*ns+3_pInt*nt) = MeanFreePathSlip0 + + forall (i = 1_pInt:ns) & + tauSlipThreshold0(i) = & + lattice_mu(ph)*plastic_disloucla_burgersPerSlipSystem(i,instance) * & + sqrt(dot_product((rhoEdge0+rhoEdgeDip0),plastic_disloucla_interactionMatrix_SlipSlip(i,1:ns,instance))) + + tempState(6_pInt*ns+4_pInt*nt+1:7_pInt*ns+4_pInt*nt) = tauSlipThreshold0 + + + +!-------------------------------------------------------------------------------------------------- +! initialize dependent twin microstructural variables + forall (j = 1_pInt:nt) & + MeanFreePathTwin0(j) = plastic_disloucla_GrainSize(instance) + tempState(6_pInt*ns+3_pInt*nt+1_pInt:6_pInt*ns+4_pInt*nt) = MeanFreePathTwin0 + + forall (j = 1_pInt:nt) & + TwinVolume0(j) = & + (pi/4.0_pReal)*plastic_disloucla_twinsizePerTwinSystem(j,instance)*MeanFreePathTwin0(j)**(2.0_pReal) + tempState(7_pInt*ns+5_pInt*nt+1_pInt:7_pInt*ns+6_pInt*nt) = TwinVolume0 + +plasticState(ph)%state0 = spread(tempState,2,size(plasticState(ph)%state(1,:))) + +end subroutine plastic_disloucla_stateInit + +!-------------------------------------------------------------------------------------------------- +!> @brief sets the relevant state values for a given instance of this plasticity +!-------------------------------------------------------------------------------------------------- +subroutine plastic_disloucla_aTolState(ph,instance) + use material, only: & + plasticState + + implicit none + integer(pInt), intent(in) :: & + ph, & + instance ! number specifying the current instance of the plasticity + + ! Tolerance state for dislocation densities + plasticState(ph)%aTolState(1_pInt:2_pInt*plastic_disloucla_totalNslip(instance)) = & + plastic_disloucla_aTolRho(instance) + + ! Tolerance state for accumulated shear due to slip + plasticState(ph)%aTolState(2_pInt*plastic_disloucla_totalNslip(instance)+1_pInt: & + 3_pInt*plastic_disloucla_totalNslip(instance))=1e6_pReal + + + ! Tolerance state for twin volume fraction + plasticState(ph)%aTolState(3_pInt*plastic_disloucla_totalNslip(instance)+1_pInt: & + 3_pInt*plastic_disloucla_totalNslip(instance)+& + plastic_disloucla_totalNtwin(instance)) = & + plastic_disloucla_aTolTwinFrac(instance) + +! Tolerance state for accumulated shear due to twin + plasticState(ph)%aTolState(3_pInt*plastic_disloucla_totalNslip(instance)+ & + plastic_disloucla_totalNtwin(instance)+1_pInt: & + 3_pInt*plastic_disloucla_totalNslip(instance)+ & + 2_pInt*plastic_disloucla_totalNtwin(instance)) = 1e6_pReal + +end subroutine plastic_disloucla_aTolState + + +!-------------------------------------------------------------------------------------------------- +!> @brief returns the homogenized elasticity matrix +!-------------------------------------------------------------------------------------------------- +function plastic_disloucla_homogenizedC(ipc,ip,el) + use material, only: & + homogenization_maxNgrains, & + phase_plasticityInstance, & + plasticState, & + mappingConstitutive + use lattice, only: & + lattice_C66 + + implicit none + real(pReal), dimension(6,6) :: & + plastic_disloucla_homogenizedC + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + + integer(pInt) :: instance,ns,nt,i, & + ph, & + of + real(pReal) :: sumf + + !* Shortened notation + of = mappingConstitutive(1,ipc,ip,el) + ph = mappingConstitutive(2,ipc,ip,el) + instance = phase_plasticityInstance(ph) + ns = plastic_disloucla_totalNslip(instance) + nt = plastic_disloucla_totalNtwin(instance) + + !* Total twin volume fraction + sumf = sum(plasticState(ph)%state((3_pInt*ns+1_pInt):(3_pInt*ns+nt),of)) ! safe for nt == 0 + !* Homogenized elasticity matrix + plastic_disloucla_homogenizedC = (1.0_pReal-sumf)*lattice_C66(1:6,1:6,ph) + do i=1_pInt,nt + plastic_disloucla_homogenizedC = plastic_disloucla_homogenizedC & + + plasticState(ph)%state(3_pInt*ns+i, of)*plastic_disloucla_Ctwin66(1:6,1:6,i,instance) + enddo + +end function plastic_disloucla_homogenizedC + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates derived quantities from state +!-------------------------------------------------------------------------------------------------- +subroutine plastic_disloucla_microstructure(temperature,ipc,ip,el) + use math, only: & + pi + use material, only: & + material_phase, & + phase_plasticityInstance, & + plasticState, & + mappingConstitutive + use lattice, only: & + lattice_mu, & + lattice_nu + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in) :: & + temperature !< temperature at IP + + integer(pInt) :: & + instance, & + ns,nt,s,t, & + ph, & + of + real(pReal) :: & + sumf,sfe,x0 + real(pReal), dimension(plastic_disloucla_totalNtwin(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: fOverStacksize + + !* Shortened notation + of = mappingConstitutive(1,ipc,ip,el) + ph = mappingConstitutive(2,ipc,ip,el) + instance = phase_plasticityInstance(ph) + ns = plastic_disloucla_totalNslip(instance) + nt = plastic_disloucla_totalNtwin(instance) + !* State: 1 : ns rho_edge + !* State: ns+1 : 2*ns rho_dipole + !* State: 2*ns+1 : 3*ns accumulated shear due to slip + !* State: 3*ns+1 : 3*ns+nt f + !* State: 3*ns+nt+1 : 3*ns+2*nt accumulated shear due to twin + !* State: 3*ns+2*nt+1 : 4*ns+2*nt 1/lambda_slip + !* State: 4*ns+2*nt+1 : 5*ns+2*nt 1/lambda_sliptwin + !* State: 5*ns+2*nt+1 : 5*ns+3*nt 1/lambda_twin + !* State: 5*ns+3*nt+1 : 6*ns+3*nt mfp_slip + !* State: 6*ns+3*nt+1 : 6*ns+4*nt mfp_twin + !* State: 6*ns+4*nt+1 : 7*ns+4*nt threshold_stress_slip + !* State: 7*ns+4*nt+1 : 7*ns+5*nt threshold_stress_twin + !* State: 7*ns+5*nt+1 : 7*ns+6*nt twin volume + + !* Total twin volume fraction + sumf = sum(plasticState(ph)%state((3*ns+1):(3*ns+nt), of)) ! safe for nt == 0 + + !* Stacking fault energy + sfe = plastic_disloucla_SFE_0K(instance) + & + plastic_disloucla_dSFE_dT(instance) * Temperature + + !* rescaled twin volume fraction for topology + forall (t = 1_pInt:nt) & + fOverStacksize(t) = & + plasticState(ph)%state(3_pInt*ns+t, of)/plastic_disloucla_twinsizePerTwinSystem(t,instance) + + !* 1/mean free distance between 2 forest dislocations seen by a moving dislocation + forall (s = 1_pInt:ns) & + plasticState(ph)%state(3_pInt*ns+2_pInt*nt+s, of) = & + sqrt(dot_product((plasticState(ph)%state(1:ns,of)+plasticState(ph)%state(ns+1_pInt:2_pInt*ns,of)),& + plastic_disloucla_forestProjectionEdge(1:ns,s,instance)))/ & + plastic_disloucla_CLambdaSlipPerSlipSystem(s,instance) + !* 1/mean free distance between 2 twin stacks from different systems seen by a moving dislocation + !$OMP CRITICAL (evilmatmul) + plasticState(ph)%state((4_pInt*ns+2_pInt*nt+1_pInt):(5_pInt*ns+2_pInt*nt), of) = 0.0_pReal + if (nt > 0_pInt .and. ns > 0_pInt) & + plasticState(ph)%state((4_pInt*ns+2_pInt*nt+1):(5_pInt*ns+2_pInt*nt), of) = & + matmul(plastic_disloucla_interactionMatrix_SlipTwin(1:ns,1:nt,instance),fOverStacksize(1:nt))/(1.0_pReal-sumf) + !$OMP END CRITICAL (evilmatmul) + + !* 1/mean free distance between 2 twin stacks from different systems seen by a growing twin + !$OMP CRITICAL (evilmatmul) + if (nt > 0_pInt) & + plasticState(ph)%state((5_pInt*ns+2_pInt*nt+1_pInt):(5_pInt*ns+3_pInt*nt), of) = & + matmul(plastic_disloucla_interactionMatrix_TwinTwin(1:nt,1:nt,instance),fOverStacksize(1:nt))/(1.0_pReal-sumf) + !$OMP END CRITICAL (evilmatmul) + + !* mean free path between 2 obstacles seen by a moving dislocation + do s = 1_pInt,ns + if (nt > 0_pInt) then + plasticState(ph)%state(5_pInt*ns+3_pInt*nt+s, of) = & + plastic_disloucla_GrainSize(instance)/(1.0_pReal+plastic_disloucla_GrainSize(instance)*& + (plasticState(ph)%state(3_pInt*ns+2_pInt*nt+s, of)+plasticState(ph)%state(4_pInt*ns+2_pInt*nt+s, of))) + else + plasticState(ph)%state(5_pInt*ns+s, of) = & + plastic_disloucla_GrainSize(instance)/& + (1.0_pReal+plastic_disloucla_GrainSize(instance)*(plasticState(ph)%state(3_pInt*ns+s, of))) + endif + enddo + + !* mean free path between 2 obstacles seen by a growing twin + forall (t = 1_pInt:nt) & + plasticState(ph)%state(6_pInt*ns+3_pInt*nt+t, of) = & + (plastic_disloucla_Cmfptwin(instance)*plastic_disloucla_GrainSize(instance))/& + (1.0_pReal+plastic_disloucla_GrainSize(instance)*plasticState(ph)%state(5_pInt*ns+2_pInt*nt+t, of)) + + !* threshold stress for dislocation motion + forall (s = 1_pInt:ns) & + plasticState(ph)%state(6_pInt*ns+4_pInt*nt+s, of) = & + lattice_mu(ph)*plastic_disloucla_burgersPerSlipSystem(s,instance)*& + sqrt(dot_product((plasticState(ph)%state(1:ns, of)+plasticState(ph)%state(ns+1_pInt:2_pInt*ns, of)),& + plastic_disloucla_interactionMatrix_SlipSlip(s,1:ns,instance))) + + !* threshold stress for growing twin + forall (t = 1_pInt:nt) & + plasticState(ph)%state(7_pInt*ns+4_pInt*nt+t, of) = & + plastic_disloucla_Cthresholdtwin(instance)*& + (sfe/(3.0_pReal*plastic_disloucla_burgersPerTwinSystem(t,instance))+& + 3.0_pReal*plastic_disloucla_burgersPerTwinSystem(t,instance)*lattice_mu(ph)/& + (plastic_disloucla_L0(instance)*plastic_disloucla_burgersPerSlipSystem(t,instance))) + + !* final twin volume after growth + forall (t = 1_pInt:nt) & + plasticState(ph)%state(7_pInt*ns+5_pInt*nt+t, of) = & + (pi/4.0_pReal)*plastic_disloucla_twinsizePerTwinSystem(t,instance)*plasticState(ph)%state(6*ns+3*nt+t, of)**(2.0_pReal) + + !* equilibrium seperation of partial dislocations + do t = 1_pInt,nt + x0 = lattice_mu(ph)*plastic_disloucla_burgersPerTwinSystem(t,instance)**(2.0_pReal)/& + (sfe*8.0_pReal*pi)*(2.0_pReal+lattice_nu(ph))/(1.0_pReal-lattice_nu(ph)) + plastic_disloucla_tau_r(t,instance)= & + lattice_mu(ph)*plastic_disloucla_burgersPerTwinSystem(t,instance)/(2.0_pReal*pi)*& + (1/(x0+plastic_disloucla_xc(instance))+cos(pi/3.0_pReal)/x0) !!! used where?? + enddo + +end subroutine plastic_disloucla_microstructure + + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates plastic velocity gradient and its tangent +!-------------------------------------------------------------------------------------------------- +subroutine plastic_disloucla_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature,nSlipDamage,slipDamage, & + ipc,ip,el) + use prec, only: & + tol_math_check + use math, only: & + math_Plain3333to99, & + math_Mandel6to33, & + math_Mandel33to6, & + math_spectralDecompositionSym33, & + math_tensorproduct, & + math_symmetric33, & + math_mul33x3 + use material, only: & + material_phase, & + phase_plasticityInstance, & + plasticState, & + mappingConstitutive + use lattice, only: & + lattice_Sslip, & + lattice_Sslip_v, & + lattice_Stwin, & + lattice_Stwin_v, & + lattice_maxNslipFamily,& + lattice_maxNtwinFamily, & + lattice_NslipSystem, & + lattice_NtwinSystem, & + lattice_NnonSchmid, & + lattice_shearTwin, & + lattice_structure, & + lattice_fcc_twinNucleationSlipPair, & + LATTICE_fcc_ID + + implicit none + integer(pInt), intent(in) :: nSlipDamage,ipc,ip,el + real(pReal), intent(in) :: Temperature + real(pReal), dimension(6), intent(in) :: Tstar_v + real(pReal), dimension(nSlipDamage), intent(in) :: slipDamage + real(pReal), dimension(3,3), intent(out) :: Lp + real(pReal), dimension(9,9), intent(out) :: dLp_dTstar99 + + integer(pInt) :: instance,ph,of,ns,nt,f,i,j,k,l,m,n,index_myFamily,s1,s2 + real(pReal) :: sumf,StressRatio_p,StressRatio_pminus1,StressRatio_r,BoltzmannRatio,DotGamma0,Ndot0, & + StressRatio_u,StressRatio_uminus1,tau_slip_pos,tau_slip_neg,vel_slip,dvel_slip,& + dgdot_dtauslip_pos,dgdot_dtauslip_neg,dgdot_dtautwin,tau_twin,gdot_twin,stressRatio + real(pReal), dimension(3,3,2) :: & + nonSchmid_tensor + real(pReal), dimension(3,3,3,3) :: & + dLp_dTstar3333 + real(pReal), dimension(plastic_disloucla_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & + gdot_slip_pos,gdot_slip_neg + + !* Shortened notation + of = mappingConstitutive(1,ipc,ip,el) + ph = mappingConstitutive(2,ipc,ip,el) + instance = phase_plasticityInstance(ph) + ns = plastic_disloucla_totalNslip(instance) + nt = plastic_disloucla_totalNtwin(instance) + + Lp = 0.0_pReal + dLp_dTstar3333 = 0.0_pReal + +!-------------------------------------------------------------------------------------------------- +! Dislocation glide part + gdot_slip_pos = 0.0_pReal + gdot_slip_neg = 0.0_pReal + dgdot_dtauslip_pos = 0.0_pReal + dgdot_dtauslip_neg = 0.0_pReal + + j = 0_pInt + slipFamilies: do f = 1_pInt,lattice_maxNslipFamily + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family + slipSystems: do i = 1_pInt,plastic_disloucla_Nslip(f,instance) + j = j+1_pInt + !* Boltzmann ratio + BoltzmannRatio = plastic_disloucla_QedgePerSlipSystem(j,instance)/(kB*Temperature) + !* Initial shear rates + DotGamma0 = & + plasticState(ph)%state(j, of)*plastic_disloucla_burgersPerSlipSystem(j,instance)*& + plastic_disloucla_v0PerSlipSystem(j,instance) + !* Resolved shear stress on slip system + tau_slip_pos = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph)) + tau_slip_neg = tau_slip_pos + nonSchmid_tensor(1:3,1:3,1) = lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph) + nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,1) + nonSchmidSystems: do k = 1,lattice_NnonSchmid(ph) + tau_slip_pos = tau_slip_pos + plastic_disloucla_nonSchmidCoeff(k,instance)* & + dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k,index_myFamily+i,ph)) + tau_slip_neg = tau_slip_neg + plastic_disloucla_nonSchmidCoeff(k,instance)* & + dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,ph)) + nonSchmid_tensor(1:3,1:3,1) = nonSchmid_tensor(1:3,1:3,1) + plastic_disloucla_nonSchmidCoeff(k,instance)*& + lattice_Sslip(1:3,1:3,2*k,index_myFamily+i,ph) + nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,2) + plastic_disloucla_nonSchmidCoeff(k,instance)*& + lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+i,ph) + enddo nonSchmidSystems + !* Applying damage to slip system + tau_slip_pos = tau_slip_pos/slipDamage(j) + tau_slip_neg = tau_slip_neg/slipDamage(j) + + significantPostitiveStress: if((abs(tau_slip_pos)-plasticState(ph)%state(6*ns+4*nt+j, of)) > tol_math_check) then + !* Stress ratios + stressRatio = ((abs(tau_slip_pos)-plasticState(ph)%state(6*ns+4*nt+j, of))/& + (plastic_disloucla_SolidSolutionStrength(instance)+& + plastic_disloucla_tau_peierlsPerSlipFamily(f,instance))) + stressRatio_p = stressRatio** plastic_disloucla_pPerSlipFamily(f,instance) + stressRatio_pminus1 = stressRatio**(plastic_disloucla_pPerSlipFamily(f,instance)-1.0_pReal) + stressRatio_u = stressRatio** plastic_disloucla_uPerSlipFamily(f,instance) + stressRatio_uminus1 = stressRatio**(plastic_disloucla_uPerSlipFamily(f,instance)-1.0_pReal) + !* Shear rates due to slip + vel_slip = plastic_disloucla_kinkheight(f,instance) * plastic_disloucla_omega(f,instance) & + * ( plasticState(ph)%state(5_pInt*ns+3_pInt*nt+j,of) - plastic_disloucla_kinkwidth(f,instance) ) & +! ^ ^ ^ +! | | position in state for all points of the same instance +! | offset in state (positin of mfp of j in state) +! state of the phase "ph" (in the case of single phase material ph=1) + / plastic_disloucla_burgersPerSlipFamily(f,instance) & + * exp(-BoltzmannRatio*(1-StressRatio_p) ** plastic_disloucla_qPerSlipFamily(f,instance)) + + gdot_slip_pos(j) = DotGamma0 & + * vel_slip & + * sign(1.0_pReal,tau_slip_pos) + + !* Derivatives of shear rates + dvel_slip = & + (abs(exp(-BoltzmannRatio*(1-StressRatio_p) ** plastic_disloucla_qPerSlipFamily(f,instance)))& + *BoltzmannRatio*plastic_disloucla_pPerSlipFamily(f,instance)& + *plastic_disloucla_qPerSlipFamily(f,instance)/& + (plastic_disloucla_SolidSolutionStrength(instance)+plastic_disloucla_tau_peierlsPerSlipFamily(f,instance))*& + StressRatio_pminus1*(1-StressRatio_p)**(plastic_disloucla_qPerSlipFamily(f,instance)-1.0_pReal) )& + * plastic_disloucla_kinkheight(f,instance) * plastic_disloucla_omega(f,instance) & + * ( plastic_disloucla_dislolength(f,instance) - plastic_disloucla_kinkwidth(f,instance) ) & + / plastic_disloucla_burgersPerSlipFamily(f,instance) + + + dgdot_dtauslip_pos = DotGamma0 * dvel_slip + + + + endif significantPostitiveStress + significantNegativeStress: if((abs(tau_slip_neg)-plasticState(ph)%state(6*ns+4*nt+j, of)) > tol_math_check) then + !* Stress ratios + stressRatio = ((abs(tau_slip_neg)-plasticState(ph)%state(6*ns+4*nt+j, of))/& + (plastic_disloucla_SolidSolutionStrength(instance)+& + plastic_disloucla_tau_peierlsPerSlipFamily(f,instance))) + stressRatio_p = stressRatio** plastic_disloucla_pPerSlipFamily(f,instance) + stressRatio_pminus1 = stressRatio**(plastic_disloucla_pPerSlipFamily(f,instance)-1.0_pReal) + stressRatio_u = stressRatio** plastic_disloucla_uPerSlipFamily(f,instance) + stressRatio_uminus1 = stressRatio**(plastic_disloucla_uPerSlipFamily(f,instance)-1.0_pReal) + !* Shear rates due to slip + vel_slip = plastic_disloucla_kinkheight(f,instance) * plastic_disloucla_omega(f,instance) & + * ( plastic_disloucla_dislolength(f,instance) - plastic_disloucla_kinkwidth(f,instance) ) & + / plastic_disloucla_burgersPerSlipFamily(f,instance) & + * exp(-BoltzmannRatio*(1-StressRatio_p) ** plastic_disloucla_qPerSlipFamily(f,instance)) + + gdot_slip_neg(j) = DotGamma0 & + * vel_slip & + * sign(1.0_pReal,tau_slip_neg) + + !* Derivatives of shear rates + dvel_slip = & + (abs(exp(-BoltzmannRatio*(1-StressRatio_p) ** plastic_disloucla_qPerSlipFamily(f,instance)))& + *BoltzmannRatio*plastic_disloucla_pPerSlipFamily(f,instance)& + *plastic_disloucla_qPerSlipFamily(f,instance)/& + (plastic_disloucla_SolidSolutionStrength(instance)+plastic_disloucla_tau_peierlsPerSlipFamily(f,instance))*& + StressRatio_pminus1*(1-StressRatio_p)**(plastic_disloucla_qPerSlipFamily(f,instance)-1.0_pReal) )& + * plastic_disloucla_kinkheight(f,instance) * plastic_disloucla_omega(f,instance) & + * ( plastic_disloucla_dislolength(f,instance) - plastic_disloucla_kinkwidth(f,instance) ) & + / plastic_disloucla_burgersPerSlipFamily(f,instance) + + + dgdot_dtauslip_neg = DotGamma0 * dvel_slip + endif significantNegativeStress + !* Plastic velocity gradient for dislocation glide + Lp = Lp + (gdot_slip_pos(j)+gdot_slip_neg(j))*0.5_pReal*lattice_Sslip(1:3,1:3,1,index_myFamily+i,ph) + !* Calculation of the tangent of Lp + forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & + dLp_dTstar3333(k,l,m,n) = & + dLp_dTstar3333(k,l,m,n) + (dgdot_dtauslip_pos*nonSchmid_tensor(m,n,1)+& + dgdot_dtauslip_neg*nonSchmid_tensor(m,n,2))*0.5_pReal*& + lattice_Sslip(k,l,1,index_myFamily+i,ph) + enddo slipSystems + enddo slipFamilies + +!-------------------------------------------------------------------------------------------------- +! correct Lp and dLp_dTstar3333 for twinned fraction + !* Total twin volume fraction + sumf = sum(plasticState(ph)%state((3_pInt*ns+1_pInt):(3_pInt*ns+nt), of)) ! safe for nt == 0 + Lp = Lp * (1.0_pReal - sumf) + dLp_dTstar3333 = dLp_dTstar3333 * (1.0_pReal - sumf) + +!-------------------------------------------------------------------------------------------------- +! Mechanical twinning part + gdot_twin = 0.0_pReal + dgdot_dtautwin = 0.0_pReal + j = 0_pInt + twinFamilies: do f = 1_pInt,lattice_maxNtwinFamily + index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,ph)) ! at which index starts my family + twinSystems: do i = 1_pInt,plastic_disloucla_Ntwin(f,instance) + j = j+1_pInt + !* Resolved shear stress on twin system + tau_twin = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,ph)) + + !* Stress ratios + if (tau_twin > tol_math_check) then + StressRatio_r = (plasticState(ph)%state(7*ns+4*nt+j, of)/tau_twin)**plastic_disloucla_rPerTwinFamily(f,instance) + !* Shear rates and their derivatives due to twin + select case(lattice_structure(ph)) + case (LATTICE_fcc_ID) + s1=lattice_fcc_twinNucleationSlipPair(1,index_myFamily+i) + s2=lattice_fcc_twinNucleationSlipPair(2,index_myFamily+i) + if (tau_twin < plastic_disloucla_tau_r(j,instance)) then + Ndot0=(abs(gdot_slip_pos(s1))*(plasticState(ph)%state(s2,of)+plasticState(ph)%state(ns+s2, of))+& !no non-Schmid behavior for fcc, just take the not influenced positive gdot_slip_pos (= gdot_slip_neg) + abs(gdot_slip_pos(s2))*(plasticState(ph)%state(s1,of)+plasticState(ph)%state(ns+s1, of)))/& + (plastic_disloucla_L0(instance)*plastic_disloucla_burgersPerSlipSystem(j,instance))*& + (1.0_pReal-exp(-plastic_disloucla_VcrossSlip(instance)/(kB*Temperature)*& + (plastic_disloucla_tau_r(j,instance)-tau_twin))) + else + Ndot0=0.0_pReal + end if + case default + Ndot0=plastic_disloucla_Ndot0PerTwinSystem(j,instance) + end select + gdot_twin = & + (plastic_disloucla_MaxTwinFraction(instance)-sumf)*lattice_shearTwin(index_myFamily+i,ph)*& + plasticState(ph)%state(7*ns+5*nt+j, of)*Ndot0*exp(-StressRatio_r) + dgdot_dtautwin = ((gdot_twin*plastic_disloucla_rPerTwinFamily(f,instance))/tau_twin)*StressRatio_r + endif + + !* Plastic velocity gradient for mechanical twinning + Lp = Lp + gdot_twin*lattice_Stwin(1:3,1:3,index_myFamily+i,ph) + + !* Calculation of the tangent of Lp + forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & + dLp_dTstar3333(k,l,m,n) = & + dLp_dTstar3333(k,l,m,n) + dgdot_dtautwin*& + lattice_Stwin(k,l,index_myFamily+i,ph)*& + lattice_Stwin(m,n,index_myFamily+i,ph) + enddo twinSystems + enddo twinFamilies + + dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333) + +end subroutine plastic_disloucla_LpAndItsTangent + + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates the rate of change of microstructure +!-------------------------------------------------------------------------------------------------- +subroutine plastic_disloucla_dotState(Tstar_v,Temperature,nSlipDamage,slipDamage,ipc,ip,el) + use prec, only: & + tol_math_check + use math, only: & + pi + use material, only: & + material_phase, & + phase_plasticityInstance, & + plasticState, & + mappingConstitutive + use lattice, only: & + lattice_Sslip_v, & + lattice_Stwin_v, & + lattice_Sslip, & + lattice_maxNslipFamily, & + lattice_maxNtwinFamily, & + lattice_NslipSystem, & + lattice_NtwinSystem, & + lattice_NnonSchmid, & + lattice_sheartwin, & + lattice_mu, & + lattice_structure, & + lattice_fcc_twinNucleationSlipPair, & + LATTICE_fcc_ID + + implicit none + real(pReal), dimension(6), intent(in):: & + Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation + real(pReal), intent(in) :: & + temperature !< temperature at integration point + integer(pInt), intent(in) :: & + nSlipDamage, & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), dimension(nSlipDamage), intent(in) :: & + slipDamage + + integer(pInt) :: instance,ns,nt,f,i,j,k,index_myFamily,s1,s2, & + ph, & + of + real(pReal) :: & + sumf, & + stressRatio_p,& + BoltzmannRatio,& + DotGamma0,& + stressRatio_u,& + stressRatio, & + EdgeDipMinDistance,& + AtomicVolume,& + VacancyDiffusion,& + StressRatio_r,& + Ndot0,& + tau_slip_pos,& + tau_slip_neg,& + DotRhoMultiplication,& + EdgeDipDistance, & + DotRhoEdgeDipAnnihilation, & + DotRhoEdgeEdgeAnnihilation, & + ClimbVelocity, & + DotRhoEdgeDipClimb, & + DotRhoDipFormation, & + tau_twin, & + vel_slip, & + gdot_slip + real(pReal), dimension(plastic_disloucla_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & + gdot_slip_pos, gdot_slip_neg + + !* Shortened notation + of = mappingConstitutive(1,ipc,ip,el) + ph = mappingConstitutive(2,ipc,ip,el) + instance = phase_plasticityInstance(ph) + ns = plastic_disloucla_totalNslip(instance) + nt = plastic_disloucla_totalNtwin(instance) + + !* Total twin volume fraction + sumf = sum(plasticState(ph)%state((3_pInt*ns+1_pInt):(3_pInt*ns+nt), of)) ! safe for nt == 0 + plasticState(ph)%dotState(:,of) = 0.0_pReal + + !* Dislocation density evolution + gdot_slip_pos = 0.0_pReal + gdot_slip_neg = 0.0_pReal + j = 0_pInt + slipFamilies: do f = 1_pInt,lattice_maxNslipFamily + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family + slipSystems: do i = 1_pInt,plastic_disloucla_Nslip(f,instance) + j = j+1_pInt + !* Boltzmann ratio + BoltzmannRatio = plastic_disloucla_QedgePerSlipSystem(j,instance)/(kB*Temperature) + !* Initial shear rates + DotGamma0 = & + plasticState(ph)%state(j, of)*plastic_disloucla_burgersPerSlipSystem(j,instance)*& + plastic_disloucla_v0PerSlipSystem(j,instance) + !* Resolved shear stress on slip system + tau_slip_pos = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,ph)) + tau_slip_neg = tau_slip_pos + + nonSchmidSystems: do k = 1,lattice_NnonSchmid(ph) + tau_slip_pos = tau_slip_pos + plastic_disloucla_nonSchmidCoeff(k,instance)* & + dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k, index_myFamily+i,ph)) + tau_slip_neg = tau_slip_neg + plastic_disloucla_nonSchmidCoeff(k,instance)* & + dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,ph)) + enddo nonSchmidSystems + tau_slip_pos = tau_slip_pos/slipDamage(j) + tau_slip_neg = tau_slip_pos/slipDamage(j) + + significantPositiveStress: if((abs(tau_slip_pos)-plasticState(ph)%state(6*ns+4*nt+j, of)) > tol_math_check) then + !* Stress ratios + stressRatio = ((abs(tau_slip_pos)-plasticState(ph)%state(6*ns+4*nt+j, of))/& + (plastic_disloucla_SolidSolutionStrength(instance)+& + plastic_disloucla_tau_peierlsPerSlipFamily(f,instance))) + stressRatio_p = stressRatio** plastic_disloucla_pPerSlipFamily(f,instance) + stressRatio_u = stressRatio** plastic_disloucla_uPerSlipFamily(f,instance) + !* Shear rates due to slip + + vel_slip = plastic_disloucla_kinkheight(f,instance) * plastic_disloucla_omega(f,instance) & + * ( plastic_disloucla_dislolength(f,instance) - plastic_disloucla_kinkwidth(f,instance) ) & + / plastic_disloucla_burgersPerSlipFamily(f,instance) & + * exp(-BoltzmannRatio*(1.0_pReal-StressRatio_p) ** plastic_disloucla_qPerSlipFamily(f,instance)) + + gdot_slip_pos(j) = DotGamma0 & + * vel_slip & + * sign(1.0_pReal,tau_slip_pos) + endif significantPositiveStress + significantNegativeStress: if((abs(tau_slip_neg)-plasticState(ph)%state(6*ns+4*nt+j, of)) > tol_math_check) then + !* Stress ratios + stressRatio = ((abs(tau_slip_neg)-plasticState(ph)%state(6*ns+4*nt+j, of))/& + (plastic_disloucla_SolidSolutionStrength(instance)+& + plastic_disloucla_tau_peierlsPerSlipFamily(f,instance))) + stressRatio_p = stressRatio** plastic_disloucla_pPerSlipFamily(f,instance) + stressRatio_u = stressRatio** plastic_disloucla_uPerSlipFamily(f,instance) + !* Shear rates due to slip + + vel_slip = plastic_disloucla_kinkheight(f,instance) * plastic_disloucla_omega(f,instance) & + * ( plastic_disloucla_dislolength(f,instance) - plastic_disloucla_kinkwidth(f,instance) ) & + / plastic_disloucla_burgersPerSlipFamily(f,instance) & + * exp(-BoltzmannRatio*(1.0_pReal-StressRatio_p) ** plastic_disloucla_qPerSlipFamily(f,instance)) + + gdot_slip_neg(j) = DotGamma0 & + * vel_slip & + * sign(1.0_pReal,tau_slip_neg) + endif significantNegativeStress + gdot_slip = (gdot_slip_pos(j)+gdot_slip_neg(j))*0.5_pReal + !* Multiplication + DotRhoMultiplication = abs(gdot_slip)/& + (plastic_disloucla_burgersPerSlipSystem(j,instance)* & + plasticState(ph)%state(5*ns+3*nt+j, of)) + + !* Dipole formation + EdgeDipMinDistance = & + plastic_disloucla_CEdgeDipMinDistance(instance)*plastic_disloucla_burgersPerSlipSystem(j,instance) + if (tau_slip_pos == 0.0_pReal) then + DotRhoDipFormation = 0.0_pReal + else + EdgeDipDistance = & + (3.0_pReal*lattice_mu(ph)*plastic_disloucla_burgersPerSlipSystem(j,instance))/& + (16.0_pReal*pi*abs(tau_slip_pos)) + if (EdgeDipDistance>plasticState(ph)%state(5*ns+3*nt+j, of)) EdgeDipDistance=plasticState(ph)%state(5*ns+3*nt+j, of) + if (EdgeDipDistance tol_math_check) then + StressRatio_r = (plasticState(ph)%state(7*ns+4*nt+j, of)/tau_twin)**plastic_disloucla_rPerTwinFamily(f,instance) + !* Shear rates and their derivatives due to twin + + select case(lattice_structure(ph)) + case (LATTICE_fcc_ID) + s1=lattice_fcc_twinNucleationSlipPair(1,index_myFamily+i) + s2=lattice_fcc_twinNucleationSlipPair(2,index_myFamily+i) + if (tau_twin < plastic_disloucla_tau_r(j,instance)) then + Ndot0=(abs(gdot_slip_pos(s1))*(plasticState(ph)%state(s2, of)+plasticState(ph)%state(ns+s2, of))+& !no non-Schmid behavior for fcc, just take the not influenced positive slip (gdot_slip_pos = gdot_slip_neg) + abs(gdot_slip_pos(s2))*(plasticState(ph)%state(s1, of)+plasticState(ph)%state(ns+s1, of)))/& + (plastic_disloucla_L0(instance)*plastic_disloucla_burgersPerSlipSystem(j,instance))*& + (1.0_pReal-exp(-plastic_disloucla_VcrossSlip(instance)/(kB*Temperature)*& + (plastic_disloucla_tau_r(j,instance)-tau_twin))) + else + Ndot0=0.0_pReal + end if + case default + Ndot0=plastic_disloucla_Ndot0PerTwinSystem(j,instance) + end select + + plasticState(ph)%dotState(3_pInt*ns+j, of) = & + (plastic_disloucla_MaxTwinFraction(instance)-sumf)*& + plasticState(ph)%state(7_pInt*ns+5_pInt*nt+j, of)*Ndot0*exp(-StressRatio_r) + !* Dotstate for accumulated shear due to twin + plasticState(ph)%dotState(3_pInt*ns+nt+j, of) = plasticState(ph)%dotState(3_pInt*ns+j, of) * & + lattice_sheartwin(index_myfamily+i,ph) + endif + enddo twinSystems + enddo twinFamilies + +end subroutine plastic_disloucla_dotState + + +!-------------------------------------------------------------------------------------------------- +!> @brief return array of constitutive results +!-------------------------------------------------------------------------------------------------- +function plastic_disloucla_postResults(Tstar_v,Temperature,ipc,ip,el) + use prec, only: & + tol_math_check + use math, only: & + pi + use material, only: & + material_phase, & + phase_plasticityInstance,& + plasticState, & + mappingConstitutive + use lattice, only: & + lattice_Sslip_v, & + lattice_Stwin_v, & + lattice_Sslip, & + lattice_maxNslipFamily, & + lattice_maxNtwinFamily, & + lattice_NslipSystem, & + lattice_NtwinSystem, & + lattice_NnonSchmid, & + lattice_shearTwin, & + lattice_mu, & + lattice_structure, & + lattice_fcc_twinNucleationSlipPair, & + LATTICE_fcc_ID + + implicit none + real(pReal), dimension(6), intent(in) :: & + Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation + real(pReal), intent(in) :: & + temperature !< temperature at integration point + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + + real(pReal), dimension(plastic_disloucla_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & + plastic_disloucla_postResults + + integer(pInt) :: & + instance,& + ns,nt,& + f,o,i,c,j,k,index_myFamily,& + s1,s2, & + ph, & + of + real(pReal) :: sumf,tau_twin,StressRatio_p,StressRatio_pminus1,& + BoltzmannRatio,DotGamma0,StressRatio_r,Ndot0,stressRatio + real(pReal) :: dvel_slip, vel_slip + real(pReal) :: StressRatio_u,StressRatio_uminus1 + real(pReal), dimension(plastic_disloucla_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & + gdot_slip_pos,dgdot_dtauslip_pos,tau_slip_pos,gdot_slip_neg,dgdot_dtauslip_neg,tau_slip_neg + + !* Shortened notation + of = mappingConstitutive(1,ipc,ip,el) + ph = mappingConstitutive(2,ipc,ip,el) + instance = phase_plasticityInstance(ph) + ns = plastic_disloucla_totalNslip(instance) + nt = plastic_disloucla_totalNtwin(instance) + + !* Total twin volume fraction + sumf = sum(plasticState(ph)%state((3_pInt*ns+1_pInt):(3_pInt*ns+nt), of)) ! safe for nt == 0 + + !* Required output + c = 0_pInt + plastic_disloucla_postResults = 0.0_pReal + + do o = 1_pInt,plastic_disloucla_Noutput(instance) + select case(plastic_disloucla_outputID(o,instance)) + + case (edge_density_ID) + plastic_disloucla_postResults(c+1_pInt:c+ns) = plasticState(ph)%state(1_pInt:ns, of) + c = c + ns + case (dipole_density_ID) + plastic_disloucla_postResults(c+1_pInt:c+ns) = plasticState(ph)%state(ns+1_pInt:2_pInt*ns, of) + c = c + ns + case (shear_rate_slip_ID,shear_rate_twin_ID,stress_exponent_ID) + gdot_slip_pos = 0.0_pReal + gdot_slip_neg = 0.0_pReal + dgdot_dtauslip_pos = 0.0_pReal + dgdot_dtauslip_neg = 0.0_pReal + j = 0_pInt + slipFamilies: do f = 1_pInt,lattice_maxNslipFamily + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family + slipSystems: do i = 1_pInt,plastic_disloucla_Nslip(f,instance) + j = j + 1_pInt + !* Boltzmann ratio + BoltzmannRatio = plastic_disloucla_QedgePerSlipSystem(j,instance)/(kB*Temperature) + !* Initial shear rates + DotGamma0 = & + plasticState(ph)%state(j, of)*plastic_disloucla_burgersPerSlipSystem(j,instance)*& + plastic_disloucla_v0PerSlipSystem(j,instance) + !* Resolved shear stress on slip system + tau_slip_pos(j) = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,ph)) + tau_slip_neg(j) = tau_slip_pos(j) + + nonSchmidSystems: do k = 1,lattice_NnonSchmid(ph) + tau_slip_pos = tau_slip_pos + plastic_disloucla_nonSchmidCoeff(k,instance)* & + dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k,index_myFamily+i,ph)) + tau_slip_neg = tau_slip_neg + plastic_disloucla_nonSchmidCoeff(k,instance)* & + dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,ph)) + enddo nonSchmidSystems + + significantPostitiveStress: if((abs(tau_slip_pos(j))-plasticState(ph)%state(6*ns+4*nt+j, of)) > tol_math_check) then + !* Stress ratios + stressRatio = ((abs(tau_slip_pos(j))-plasticState(ph)%state(6*ns+4*nt+j, of))/& + (plastic_disloucla_SolidSolutionStrength(instance)+& + plastic_disloucla_tau_peierlsPerSlipFamily(f,instance))) + stressRatio_p = stressRatio** plastic_disloucla_pPerSlipFamily(f,instance) + stressRatio_pminus1 = stressRatio**(plastic_disloucla_pPerSlipFamily(f,instance)-1.0_pReal) + stressRatio_u = stressRatio** plastic_disloucla_uPerSlipFamily(f,instance) + stressRatio_uminus1 = stressRatio**(plastic_disloucla_uPerSlipFamily(f,instance)-1.0_pReal) + !* Shear rates due to slip + vel_slip = plastic_disloucla_kinkheight(f,instance) * plastic_disloucla_omega(f,instance) & + * ( plastic_disloucla_dislolength(f,instance) - plastic_disloucla_kinkwidth(f,instance) ) & + / plastic_disloucla_burgersPerSlipFamily(f,instance) & + * exp(-BoltzmannRatio*(1.0_pReal-StressRatio_p) ** plastic_disloucla_qPerSlipFamily(f,instance)) + + gdot_slip_pos(j) = DotGamma0 & + * vel_slip & + * sign(1.0_pReal,tau_slip_pos(j)) + !* Derivatives of shear rates + + dvel_slip = & + (abs(exp(-BoltzmannRatio*(1.0_pReal-StressRatio_p) ** plastic_disloucla_qPerSlipFamily(f,instance)))& + *BoltzmannRatio*plastic_disloucla_pPerSlipFamily(f,instance)& + *plastic_disloucla_qPerSlipFamily(f,instance)/& + (plastic_disloucla_SolidSolutionStrength(instance)+plastic_disloucla_tau_peierlsPerSlipFamily(f,instance))*& + StressRatio_pminus1*(1.0_pReal-StressRatio_p)**(plastic_disloucla_qPerSlipFamily(f,instance)-1.0_pReal) )& + * plastic_disloucla_kinkheight(f,instance) * plastic_disloucla_omega(f,instance) & + * ( plastic_disloucla_dislolength(f,instance) - plastic_disloucla_kinkwidth(f,instance) ) & + / plastic_disloucla_burgersPerSlipFamily(f,instance) + + + dgdot_dtauslip_pos(j) = DotGamma0 * dvel_slip + + endif significantPostitiveStress + significantNegativeStress: if((abs(tau_slip_neg(j))-plasticState(ph)%state(6*ns+4*nt+j, of)) > tol_math_check) then + !* Stress ratios + stressRatio = ((abs(tau_slip_neg(j))-plasticState(ph)%state(6*ns+4*nt+j, of))/& + (plastic_disloucla_SolidSolutionStrength(instance)+& + plastic_disloucla_tau_peierlsPerSlipFamily(f,instance))) + stressRatio_p = stressRatio** plastic_disloucla_pPerSlipFamily(f,instance) + stressRatio_pminus1 = stressRatio**(plastic_disloucla_pPerSlipFamily(f,instance)-1.0_pReal) + stressRatio_u = stressRatio** plastic_disloucla_uPerSlipFamily(f,instance) + stressRatio_uminus1 = stressRatio**(plastic_disloucla_uPerSlipFamily(f,instance)-1.0_pReal) + !* Shear rates due to slip + vel_slip = plastic_disloucla_kinkheight(f,instance) * plastic_disloucla_omega(f,instance) & + * ( plastic_disloucla_dislolength(f,instance) - plastic_disloucla_kinkwidth(f,instance) ) & + / plastic_disloucla_burgersPerSlipFamily(f,instance) & + * exp(-BoltzmannRatio*(1.0_pReal-StressRatio_p) ** plastic_disloucla_qPerSlipFamily(f,instance)) + + gdot_slip_neg(j) = DotGamma0 & + * vel_slip & + * sign(1.0_pReal,tau_slip_neg(j)) + !* Derivatives of shear rates + dvel_slip = & + (abs(exp(-BoltzmannRatio*(1-StressRatio_p) ** plastic_disloucla_qPerSlipFamily(f,instance)))& + *BoltzmannRatio*plastic_disloucla_pPerSlipFamily(f,instance)& + *plastic_disloucla_qPerSlipFamily(f,instance)/& + (plastic_disloucla_SolidSolutionStrength(instance)+plastic_disloucla_tau_peierlsPerSlipFamily(f,instance))*& + StressRatio_pminus1*(1-StressRatio_p)**(plastic_disloucla_qPerSlipFamily(f,instance)-1.0_pReal) )& + * plastic_disloucla_kinkheight(f,instance) * plastic_disloucla_omega(f,instance) & + * ( plastic_disloucla_dislolength(f,instance) - plastic_disloucla_kinkwidth(f,instance) ) & + / plastic_disloucla_burgersPerSlipFamily(f,instance) + + + dgdot_dtauslip_neg(j) = DotGamma0 * dvel_slip + + endif significantNegativeStress + enddo slipSystems + enddo slipFamilies + + if (plastic_disloucla_outputID(o,instance) == shear_rate_slip_ID) then + plastic_disloucla_postResults(c+1:c+ns) = (gdot_slip_pos + gdot_slip_neg)*0.5_pReal + c = c + ns + elseif (plastic_disloucla_outputID(o,instance) == shear_rate_twin_ID) then + if (nt > 0_pInt) then + j = 0_pInt + twinFamilies1: do f = 1_pInt,lattice_maxNtwinFamily + index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,ph)) ! at which index starts my family + twinSystems1: do i = 1,plastic_disloucla_Ntwin(f,instance) + j = j + 1_pInt + + !* Resolved shear stress on twin system + tau_twin = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,ph)) + !* Stress ratios + StressRatio_r = (plasticState(ph)%state(7_pInt*ns+4_pInt*nt+j, of)/ & + tau_twin)**plastic_disloucla_rPerTwinFamily(f,instance) + + !* Shear rates due to twin + if ( tau_twin > 0.0_pReal ) then + select case(lattice_structure(ph)) + case (LATTICE_fcc_ID) + s1=lattice_fcc_twinNucleationSlipPair(1,index_myFamily+i) + s2=lattice_fcc_twinNucleationSlipPair(2,index_myFamily+i) + if (tau_twin < plastic_disloucla_tau_r(j,instance)) then + Ndot0=(abs(gdot_slip_pos(s1))*(plasticState(ph)%state(s2, of)+plasticState(ph)%state(ns+s2, of))+& !no non-Schmid behavior for fcc, just take the not influenced positive slip (gdot_slip_pos = gdot_slip_neg) + abs(gdot_slip_pos(s2))*(plasticState(ph)%state(s1, of)+plasticState(ph)%state(ns+s1, of)))/& + (plastic_disloucla_L0(instance)*& + plastic_disloucla_burgersPerSlipSystem(j,instance))*& + (1.0_pReal-exp(-plastic_disloucla_VcrossSlip(instance)/(kB*Temperature)*& + (plastic_disloucla_tau_r(j,instance)-tau_twin))) + else + Ndot0=0.0_pReal + end if + + case default + Ndot0=plastic_disloucla_Ndot0PerTwinSystem(j,instance) + end select + plastic_disloucla_postResults(c+j) = & + (plastic_disloucla_MaxTwinFraction(instance)-sumf)*lattice_shearTwin(index_myFamily+i,ph)*& + plasticState(ph)%state(7_pInt*ns+5_pInt*nt+j, of)*Ndot0*exp(-StressRatio_r) + endif + enddo twinSystems1 + enddo twinFamilies1 + endif + c = c + nt + elseif(plastic_disloucla_outputID(o,instance) == stress_exponent_ID) then + do j = 1_pInt, ns + if ((gdot_slip_pos(j)+gdot_slip_neg(j))*0.5_pReal==0.0_pReal) then + plastic_disloucla_postResults(c+j) = 0.0_pReal + else + plastic_disloucla_postResults(c+j) = (tau_slip_pos(j)+tau_slip_neg(j))/& + (gdot_slip_pos(j)+gdot_slip_neg(j))*& + (dgdot_dtauslip_pos(j)+dgdot_dtauslip_neg(j))* 0.5_pReal + endif + enddo + c = c + ns + endif + + case (accumulated_shear_slip_ID) + plastic_disloucla_postResults(c+1_pInt:c+ns) = & + plasticState(ph)%state((2_pInt*ns+1_pInt):(3_pInt*ns), of) + c = c + ns + case (mfp_slip_ID) + plastic_disloucla_postResults(c+1_pInt:c+ns) =& + plasticState(ph)%state((5_pInt*ns+3_pInt*nt+1_pInt):(6_pInt*ns+3_pInt*nt), of) + c = c + ns + case (resolved_stress_slip_ID) + j = 0_pInt + slipFamilies1: do f = 1_pInt,lattice_maxNslipFamily + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family + slipSystems1: do i = 1_pInt,plastic_disloucla_Nslip(f,instance) + j = j + 1_pInt + plastic_disloucla_postResults(c+j) =& + dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,ph)) + enddo slipSystems1; enddo slipFamilies1 + c = c + ns + case (threshold_stress_slip_ID) + plastic_disloucla_postResults(c+1_pInt:c+ns) = & + plasticState(ph)%state((6_pInt*ns+4_pInt*nt+1_pInt):(7_pInt*ns+4_pInt*nt), of) + c = c + ns + case (edge_dipole_distance_ID) + j = 0_pInt + slipFamilies2: do f = 1_pInt,lattice_maxNslipFamily + index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) ! at which index starts my family + slipSystems2: do i = 1_pInt,plastic_disloucla_Nslip(f,instance) + j = j + 1_pInt + plastic_disloucla_postResults(c+j) = & + (3.0_pReal*lattice_mu(ph)*plastic_disloucla_burgersPerSlipSystem(j,instance))/& + (16.0_pReal*pi*abs(dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,ph)))) + plastic_disloucla_postResults(c+j)=min(plastic_disloucla_postResults(c+j),& + plasticState(ph)%state(5*ns+3*nt+j, of)) + enddo slipSystems2; enddo slipFamilies2 + c = c + ns + case (twin_fraction_ID) + plastic_disloucla_postResults(c+1_pInt:c+nt) = plasticState(ph)%state((3_pInt*ns+1_pInt):(3_pInt*ns+nt), of) + c = c + nt + + case (accumulated_shear_twin_ID) + plastic_disloucla_postResults(c+1_pInt:c+nt) = plasticState(ph)% & + state((3_pInt*ns+nt+1_pInt) :(3_pInt*ns+2_pInt*nt), of) + c = c + nt + case (mfp_twin_ID) + plastic_disloucla_postResults(c+1_pInt:c+nt) = plasticState(ph)% & + state((6_pInt*ns+3_pInt*nt+1_pInt):(6_pInt*ns+4_pInt*nt), of) + c = c + nt + case (resolved_stress_twin_ID) + if (nt > 0_pInt) then + j = 0_pInt + twinFamilies2: do f = 1_pInt,lattice_maxNtwinFamily + index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,ph)) ! at which index starts my family + twinSystems2: do i = 1_pInt,plastic_disloucla_Ntwin(f,instance) + j = j + 1_pInt + plastic_disloucla_postResults(c+j) = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,ph)) + enddo twinSystems2; enddo twinFamilies2 + endif + c = c + nt + case (threshold_stress_twin_ID) + plastic_disloucla_postResults(c+1_pInt:c+nt) = plasticState(ph)% & + state((7_pInt*ns+4_pInt*nt+1_pInt):(7_pInt*ns+5_pInt*nt), of) + c = c + nt + end select + enddo +end function plastic_disloucla_postResults + +end module plastic_disloucla