diff --git a/VERSION b/VERSION index 8d5912448..6efd0b994 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.2-1277-g53bc24cc +v2.0.2-1291-g19df6f8a diff --git a/src/config.f90 b/src/config.f90 index c7fd95b43..dcc14e015 100644 --- a/src/config.f90 +++ b/src/config.f90 @@ -318,7 +318,7 @@ subroutine show(this) do while (associated(item%next)) write(6,'(a)') ' '//trim(item%string%val) item => item%next - end do + enddo end subroutine show @@ -391,7 +391,7 @@ logical function keyExists(this,key) do while (associated(item%next) .and. .not. keyExists) keyExists = trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key) item => item%next - end do + enddo end function keyExists @@ -417,7 +417,7 @@ integer(pInt) function countKeys(this,key) if (trim(IO_stringValue(item%string%val,item%string%pos,1)) == trim(key)) & countKeys = countKeys + 1_pInt item => item%next - end do + enddo end function countKeys @@ -451,7 +451,7 @@ real(pReal) function getFloat(this,key,defaultVal) getFloat = IO_FloatValue(item%string%val,item%string%pos,2) endif item => item%next - end do + enddo if (.not. found) call IO_error(140_pInt,ext_msg=key) @@ -487,7 +487,7 @@ integer(pInt) function getInt(this,key,defaultVal) getInt = IO_IntValue(item%string%val,item%string%pos,2) endif item => item%next - end do + enddo if (.not. found) call IO_error(140_pInt,ext_msg=key) @@ -538,7 +538,7 @@ character(len=65536) function getString(this,key,defaultVal,raw) endif endif item => item%next - end do + enddo if (.not. found) call IO_error(140_pInt,ext_msg=key) @@ -584,7 +584,7 @@ function getFloats(this,key,defaultVal,requiredShape,requiredSize) enddo endif item => item%next - end do + enddo if (.not. found) then if (present(defaultVal)) then; getFloats = defaultVal; else; call IO_error(140_pInt,ext_msg=key); endif @@ -635,7 +635,7 @@ function getInts(this,key,defaultVal,requiredShape,requiredSize) enddo endif item => item%next - end do + enddo if (.not. found) then if (present(defaultVal)) then; getInts = defaultVal; else; call IO_error(140_pInt,ext_msg=key); endif @@ -712,7 +712,7 @@ function getStrings(this,key,defaultVal,requiredShape,raw) endif notAllocated endif item => item%next - end do + enddo if (.not. found) then if (present(defaultVal)) then; getStrings = defaultVal; else; call IO_error(140_pInt,ext_msg=key); endif diff --git a/src/plastic_disloUCLA.f90 b/src/plastic_disloUCLA.f90 index 9006b092d..67adb083b 100644 --- a/src/plastic_disloUCLA.f90 +++ b/src/plastic_disloUCLA.f90 @@ -3,8 +3,7 @@ !> @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 +!> @brief crystal plasticity model for bcc metals, especially Tungsten !-------------------------------------------------------------------------------------------------- module plastic_disloUCLA use prec, only: & @@ -15,7 +14,6 @@ module plastic_disloUCLA private 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 @@ -23,46 +21,46 @@ module plastic_disloUCLA kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin enum, bind(c) - enumerator :: undefined_ID, & + enumerator :: & + undefined_ID, & rho_ID, & rhoDip_ID, & shearrate_ID, & accumulatedshear_ID, & mfp_ID, & - resolvedstress_ID, & thresholdstress_ID, & - dipoledistance_ID, & - stressexponent_ID + dipoledistance_ID end enum type, private :: tParameters real(pReal) :: & aTolRho, & grainSize, & - SolidSolutionStrength, & !< Strength due to elements in solid solution + SolidSolutionStrength, & !< Strength due to elements in solid solution mu, & - D0, & !< prefactor for self-diffusion coefficient - Qsd !< activation energy for dislocation climb + D0, & !< prefactor for self-diffusion coefficient + Qsd !< activation energy for dislocation climb real(pReal), allocatable, dimension(:) :: & - rho0, & !< initial edge dislocation density per slip system for each family and instance - rhoDip0, & !< initial edge dipole density per slip system for each family and instance - burgers, & !< absolute length of burgers vector [m] for each slip system and instance + rho0, & !< initial edge dislocation density + rhoDip0, & !< initial edge dipole density + burgers, & !< absolute length of burgers vector [m] nonSchmidCoeff, & minDipDistance, & - CLambda, & !< Adj. parameter for distance between 2 forest dislocations for each slip system and instance + CLambda, & !< Adj. parameter for distance between 2 forest dislocations atomicVolume, & + tau_Peierls, & + tau0, & !* mobility law parameters - H0kp, & !< activation energy for glide [J] for each slip system and instance - v0, & !< dislocation velocity prefactor [m/s] for each family and instance - p, & !< p-exponent in glide velocity - q, & !< q-exponent in glide velocity - B, & !< friction coeff. B (kMC) - kink_height, & !< height of the kink pair - kink_width, & !< width of the kink pair - omega, & !< attempt frequency for kink pair nucleation - tau_Peierls + H0kp, & !< activation energy for glide [J] + v0, & !< dislocation velocity prefactor [m/s] + p, & !< p-exponent in glide velocity + q, & !< q-exponent in glide velocity + B, & !< friction coefficient + kink_height, & !< height of the kink pair + w, & !< width of the kink pair + omega !< attempt frequency for kink pair nucleation real(pReal), allocatable, dimension(:,:) :: & - interaction_SlipSlip, & !< slip resistance from slip activity + interaction_SlipSlip, & !< slip resistance from slip activity forestProjectionEdge real(pReal), allocatable, dimension(:,:,:) :: & Schmid_slip, & @@ -85,20 +83,19 @@ module plastic_disloUCLA real(pReal), pointer, dimension(:,:) :: & rhoEdge, & rhoEdgeDip, & - accshear_slip, & - whole + accshear_slip end type type, private :: tDisloUCLAdependentState - real(pReal), allocatable, dimension(:,:) :: & + real(pReal), allocatable, dimension(:,:) :: & mfp, & + dislocationSpacing, & threshold_stress end type tDisloUCLAdependentState type(tDisloUCLAState ), allocatable, dimension(:), private :: & - state, & - dotState - + dotState, & + state type(tDisloUCLAdependentState), allocatable, dimension(:), private :: & dependentState @@ -111,10 +108,8 @@ module plastic_disloUCLA private :: & kinetics - contains - !-------------------------------------------------------------------------------------------------- !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks @@ -141,42 +136,46 @@ subroutine plastic_disloUCLA_init() phase_plasticity, & phase_plasticityInstance, & phase_Noutput, & + material_allocatePlasticState, & PLASTICITY_DISLOUCLA_label, & PLASTICITY_DISLOUCLA_ID, & material_phase, & - plasticState, & - material_allocatePlasticState + plasticState use config, only: & MATERIAL_partPhase, & config_phase use lattice implicit none - integer(pInt) :: Ninstance,& - f,j,k,o, i, & - outputSize, & - offset_slip, index_myFamily, index_otherFamily, & - startIndex, endIndex, p, & - sizeState, sizeDotState, & - NipcMyPhase - character(len=pStringLen) :: & - structure = '',& - extmsg = '' - character(len=65536), dimension(:), allocatable :: outputs - integer(kind(undefined_ID)) :: outputID + integer(pInt) :: & + index_myFamily, index_otherFamily, & + f,j,k,o, & + Ninstance, & + p, i, & + NipcMyPhase, outputSize, & + sizeState, sizeDotState, & + startIndex, endIndex + integer(pInt), dimension(0), parameter :: emptyIntArray = [integer(pInt)::] real(pReal), dimension(0), parameter :: emptyRealArray = [real(pReal)::] character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] - write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_DISLOUCLA_label//' init -+>>>' + integer(kind(undefined_ID)) :: & + outputID + + character(len=pStringLen) :: & + structure = '',& + extmsg = '' + character(len=65536), dimension(:), allocatable :: & + outputs + + write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_DISLOUCLA_label//' init -+>>>' write(6,'(/,a)') ' Cereceda et al., International Journal of Plasticity 78, 2016, 242-256' write(6,'(/,a)') ' http://dx.doi.org/10.1016/j.ijplas.2015.09.002' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" Ninstance = int(count(phase_plasticity == PLASTICITY_DISLOUCLA_ID),pInt) - if (Ninstance == 0_pInt) return - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & write(6,'(a16,1x,i5,/)') '# instances:',Ninstance @@ -190,26 +189,34 @@ subroutine plastic_disloUCLA_init() allocate(dependentState(Ninstance)) - do p = 1_pInt, size(phase_plasticityInstance) + do p = 1_pInt, size(phase_plasticity) if (phase_plasticity(p) /= PLASTICITY_DISLOUCLA_ID) cycle associate(prm => param(phase_plasticityInstance(p)), & dot => dotState(phase_plasticityInstance(p)), & stt => state(phase_plasticityInstance(p)), & - dst => dependentState(phase_plasticityInstance(p))) + dst => dependentState(phase_plasticityInstance(p)), & + config => config_phase(p)) - structure = config_phase(p)%getString('lattice_structure') + structure = config%getString('lattice_structure') + +!-------------------------------------------------------------------------------------------------- +! optional parameters that need to be defined prm%mu = lattice_mu(p) - prm%aTolRho = config_phase(p)%getFloat('atol_rho') + prm%aTolRho = config%getFloat('atol_rho') + + ! sanity checks + if (prm%aTolRho <= 0.0_pReal) extmsg = trim(extmsg)//' atol_rho' + !-------------------------------------------------------------------------------------------------- ! slip related parameters - prm%Nslip = config_phase(p)%getInts('nslip',defaultVal=emptyIntArray) + prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) prm%totalNslip = sum(prm%Nslip) slipActive: if (prm%totalNslip > 0_pInt) then prm%Schmid_slip = lattice_SchmidMatrix_slip(prm%Nslip,structure(1:3),& - config_phase(p)%getFloat('c/a',defaultVal=0.0_pReal)) + config%getFloat('c/a',defaultVal=0.0_pReal)) if(structure=='bcc') then - prm%nonSchmidCoeff = config_phase(p)%getFloats('nonschmid_coefficients',& + prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',& defaultVal = emptyRealArray) prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt) prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt) @@ -218,32 +225,32 @@ subroutine plastic_disloUCLA_init() prm%nonSchmid_neg = prm%Schmid_slip endif prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, & - config_phase(p)%getFloats('interaction_slipslip'), & + config%getFloats('interaction_slipslip'), & structure(1:3)) - prm%rho0 = config_phase(p)%getFloats('rhoedge0', requiredShape=shape(prm%Nslip)) - prm%rhoDip0 = config_phase(p)%getFloats('rhoedgedip0', requiredShape=shape(prm%Nslip)) - prm%v0 = config_phase(p)%getFloats('v0', requiredShape=shape(prm%Nslip)) - prm%burgers = config_phase(p)%getFloats('slipburgers', requiredShape=shape(prm%Nslip)) - prm%H0kp = config_phase(p)%getFloats('qedge', requiredShape=shape(prm%Nslip)) + prm%rho0 = config%getFloats('rhoedge0', requiredShape=shape(prm%Nslip)) + prm%rhoDip0 = config%getFloats('rhoedgedip0', requiredShape=shape(prm%Nslip)) + prm%v0 = config%getFloats('v0', requiredShape=shape(prm%Nslip)) + prm%burgers = config%getFloats('slipburgers', requiredShape=shape(prm%Nslip)) + prm%H0kp = config%getFloats('qedge', requiredShape=shape(prm%Nslip)) - prm%clambda = config_phase(p)%getFloats('clambdaslip', requiredShape=shape(prm%Nslip)) - prm%tau_Peierls = config_phase(p)%getFloats('tau_peierls', requiredShape=shape(prm%Nslip)) - prm%p = config_phase(p)%getFloats('p_slip', requiredShape=shape(prm%Nslip), & - defaultVal=[(1.0_pReal,i=1_pInt,size(prm%Nslip))]) - prm%q = config_phase(p)%getFloats('q_slip', requiredShape=shape(prm%Nslip), & - defaultVal=[(1.0_pReal,i=1_pInt,size(prm%Nslip))]) - prm%kink_height = config_phase(p)%getFloats('kink_height', requiredShape=shape(prm%Nslip)) - prm%kink_width = config_phase(p)%getFloats('kink_width', requiredShape=shape(prm%Nslip)) - prm%omega = config_phase(p)%getFloats('omega', requiredShape=shape(prm%Nslip)) - prm%B = config_phase(p)%getFloats('friction_coeff', requiredShape=shape(prm%Nslip)) + prm%clambda = config%getFloats('clambdaslip', requiredShape=shape(prm%Nslip)) + prm%tau_Peierls = config%getFloats('tau_peierls', requiredShape=shape(prm%Nslip)) ! ToDo: Deprecated + prm%p = config%getFloats('p_slip', requiredShape=shape(prm%Nslip), & + defaultVal=[(1.0_pReal,i=1_pInt,size(prm%Nslip))]) + prm%q = config%getFloats('q_slip', requiredShape=shape(prm%Nslip), & + defaultVal=[(1.0_pReal,i=1_pInt,size(prm%Nslip))]) + prm%kink_height = config%getFloats('kink_height', requiredShape=shape(prm%Nslip)) + prm%w = config%getFloats('kink_width', requiredShape=shape(prm%Nslip)) + prm%omega = config%getFloats('omega', requiredShape=shape(prm%Nslip)) + prm%B = config%getFloats('friction_coeff', requiredShape=shape(prm%Nslip)) - prm%SolidSolutionStrength = config_phase(p)%getFloat('solidsolutionstrength') - prm%grainSize = config_phase(p)%getFloat('grainsize') - prm%D0 = config_phase(p)%getFloat('d0') - prm%Qsd = config_phase(p)%getFloat('qsd') - prm%atomicVolume = config_phase(p)%getFloat('catomicvolume') * prm%burgers**3.0_pReal - prm%minDipDistance = config_phase(p)%getFloat('cedgedipmindistance') * prm%burgers - prm%dipoleformation = config_phase(p)%getFloat('dipoleformationfactor') > 0.0_pReal !should be on by default + prm%SolidSolutionStrength = config%getFloat('solidsolutionstrength') ! ToDo: Deprecated + prm%grainSize = config%getFloat('grainsize') + prm%D0 = config%getFloat('d0') + prm%Qsd = config%getFloat('qsd') + prm%atomicVolume = config%getFloat('catomicvolume') * prm%burgers**3.0_pReal + prm%minDipDistance = config%getFloat('cedgedipmindistance') * prm%burgers + prm%dipoleformation = config%getFloat('dipoleformationfactor') > 0.0_pReal !should be on by default, ToDo: change to /key/-key ! expand: family => system prm%rho0 = math_expand(prm%rho0, prm%Nslip) @@ -253,7 +260,7 @@ subroutine plastic_disloUCLA_init() prm%H0kp = math_expand(prm%H0kp, prm%Nslip) prm%burgers = math_expand(prm%burgers, prm%Nslip) prm%kink_height = math_expand(prm%kink_height, prm%Nslip) - prm%kink_width = math_expand(prm%kink_width, prm%Nslip) + prm%w = math_expand(prm%w, prm%Nslip) prm%omega = math_expand(prm%omega, prm%Nslip) prm%tau_Peierls = math_expand(prm%tau_Peierls, prm%Nslip) prm%v0 = math_expand(prm%v0, prm%Nslip) @@ -261,43 +268,40 @@ subroutine plastic_disloUCLA_init() prm%clambda = math_expand(prm%clambda, prm%Nslip) prm%atomicVolume = math_expand(prm%atomicVolume, prm%Nslip) prm%minDipDistance = math_expand(prm%minDipDistance, prm%Nslip) + + prm%tau0 = prm%tau_peierls + prm%SolidSolutionStrength ! sanity checks - if (any(prm%rho0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoedge0' - if (any(prm%rhoDip0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoedgedip0' - if (any(prm%v0 < 0.0_pReal)) extmsg = trim(extmsg)//' v0' - if (any(prm%burgers <= 0.0_pReal)) extmsg = trim(extmsg)//' slipburgers' - if (any(prm%H0kp <= 0.0_pReal)) extmsg = trim(extmsg)//' qedge' - if (any(prm%tau_peierls < 0.0_pReal)) extmsg = trim(extmsg)//' tau_peierls' - - if ( prm%D0 <= 0.0_pReal) extmsg = trim(extmsg)//' d0' - if ( prm%Qsd <= 0.0_pReal) extmsg = trim(extmsg)//' qsd' + if (any(prm%rho0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoedge0' + if (any(prm%rhoDip0 < 0.0_pReal)) extmsg = trim(extmsg)//' rhoedgedip0' + if (any(prm%v0 < 0.0_pReal)) extmsg = trim(extmsg)//' v0' + if (any(prm%burgers <= 0.0_pReal)) extmsg = trim(extmsg)//' slipburgers' + if (any(prm%H0kp <= 0.0_pReal)) extmsg = trim(extmsg)//' qedge' + if (any(prm%tau_peierls < 0.0_pReal)) extmsg = trim(extmsg)//' tau_peierls' + if ( prm%D0 <= 0.0_pReal) extmsg = trim(extmsg)//' d0' + if ( prm%Qsd <= 0.0_pReal) extmsg = trim(extmsg)//' qsd' !if (plastic_disloUCLA_CAtomicVolume(instance) <= 0.0_pReal) & ! call IO_error(211_pInt,el=instance,ext_msg='cAtomicVolume ('//PLASTICITY_DISLOUCLA_label//')') - - ! if (plastic_disloUCLA_aTolRho(instance) <= 0.0_pReal) & - ! call IO_error(211_pInt,el=instance,ext_msg='aTolRho ('//PLASTICITY_DISLOUCLA_label//')') - else slipActive allocate(prm%rho0(0)) allocate(prm%rhoDip0(0)) endif slipActive +!-------------------------------------------------------------------------------------------------- +! exit if any parameter is out of range + if (extmsg /= '') & + call IO_error(211_pInt,ext_msg=trim(extmsg)//'('//PLASTICITY_DISLOUCLA_label//')') -#if defined(__GFORTRAN__) - outputs = ['GfortranBug86277'] - outputs = config_phase(p)%getStrings('(output)',defaultVal=outputs) - if (outputs(1) == 'GfortranBug86277') outputs = emptyStringArray -#else - outputs = config_phase(p)%getStrings('(output)',defaultVal=emptyStringArray) -#endif +!-------------------------------------------------------------------------------------------------- +! output pararameters + outputs = config%getStrings('(output)',defaultVal=emptyStringArray) allocate(prm%outputID(0)) - - do i = 1_pInt, size(outputs) + do i=1_pInt, size(outputs) outputID = undefined_ID outputSize = prm%totalNslip select case(trim(outputs(i))) + case ('edge_density') outputID = merge(rho_ID,undefined_ID,prm%totalNslip>0_pInt) case ('dipole_density') @@ -308,14 +312,11 @@ subroutine plastic_disloUCLA_init() outputID = merge(accumulatedshear_ID,undefined_ID,prm%totalNslip>0_pInt) case ('mfp','mfp_slip') outputID = merge(mfp_ID,undefined_ID,prm%totalNslip>0_pInt) - case ('resolved_stress','resolved_stress_slip') - outputID = merge(resolvedstress_ID,undefined_ID,prm%totalNslip>0_pInt) case ('threshold_stress','threshold_stress_slip') outputID = merge(thresholdstress_ID,undefined_ID,prm%totalNslip>0_pInt) case ('edge_dipole_distance') outputID = merge(dipoleDistance_ID,undefined_ID,prm%totalNslip>0_pInt) - case ('stress_exponent') - outputID = merge(stressexponent_ID,undefined_ID,prm%totalNslip>0_pInt) + end select if (outputID /= undefined_ID) then @@ -326,17 +327,14 @@ subroutine plastic_disloUCLA_init() enddo - NipcMyPhase=count(material_phase==p) - !-------------------------------------------------------------------------------------------------- ! allocate state arrays - + NipcMyPhase = count(material_phase==p) sizeDotState = int(size(['rhoEdge ','rhoEdgeDip ','accshearslip']),pInt) * prm%totalNslip - sizeState = sizeDotState + sizeState = sizeDotState call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0_pInt, & prm%totalNslip,0_pInt,0_pInt) - plasticState(p)%sizePostResults = sum(plastic_disloUCLA_sizePostResult(:,phase_plasticityInstance(p))) allocate(prm%forestProjectionEdge(prm%totalNslip,prm%totalNslip),source = 0.0_pReal) @@ -355,43 +353,41 @@ subroutine plastic_disloUCLA_init() lattice_st(:,sum(lattice_NslipSystem(1:o-1,p))+k,p))) enddo; enddo enddo slipSystemsLoop - enddo mySlipFamilies + enddo mySlipFamilies - offset_slip = 2_pInt*plasticState(p)%nSlip - plasticState(p)%slipRate => & - plasticState(p)%dotState(offset_slip+1:offset_slip+plasticState(p)%nSlip,1:NipcMyPhase) - plasticState(p)%accumulatedSlip => & - plasticState(p)%state (offset_slip+1:offset_slip+plasticState(p)%nSlip,1:NipcMyPhase) - - startIndex=1_pInt - endIndex=prm%totalNslip +!-------------------------------------------------------------------------------------------------- +! locally defined state aliases and initialization of state0 and aTolState + startIndex = 1_pInt + endIndex = prm%totalNslip stt%rhoEdge=>plasticState(p)%state(startIndex:endIndex,:) stt%rhoEdge= spread(prm%rho0,2,NipcMyPhase) dot%rhoEdge=>plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolRho - startIndex=endIndex+1_pInt - endIndex=endIndex+prm%totalNslip + startIndex = endIndex + 1_pInt + endIndex = endIndex + prm%totalNslip stt%rhoEdgeDip=>plasticState(p)%state(startIndex:endIndex,:) stt%rhoEdgeDip= spread(prm%rhoDip0,2,NipcMyPhase) dot%rhoEdgeDip=>plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = prm%aTolRho - startIndex=endIndex+1_pInt - endIndex=endIndex+prm%totalNslip + startIndex = endIndex + 1_pInt + endIndex = endIndex + prm%totalNslip stt%accshear_slip=>plasticState(p)%state(startIndex:endIndex,:) dot%accshear_slip=>plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%aTolState(startIndex:endIndex) = 1e6_pReal + ! global alias + plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) + plasticState(p)%accumulatedSlip => plasticState(p)%state(startIndex:endIndex,:) - dot%whole => plasticState(p)%dotState + allocate(dst%mfp(prm%totalNslip,NipcMyPhase), source=0.0_pReal) + allocate(dst%dislocationSpacing(prm%totalNslip,NipcMyPhase),source=0.0_pReal) + allocate(dst%threshold_stress(prm%totalNslip,NipcMyPhase), source=0.0_pReal) + plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally - allocate(dst%mfp(prm%totalNslip,NipcMyPhase),source=0.0_pReal) - allocate(dst%threshold_stress(prm%totalNslip,NipcMyPhase),source=0.0_pReal) - - - plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally end associate + enddo end subroutine plastic_disloUCLA_init @@ -403,25 +399,24 @@ end subroutine plastic_disloUCLA_init subroutine plastic_disloUCLA_dependentState(instance,of) implicit none - integer(pInt), intent(in) :: instance, of + integer(pInt), intent(in) :: instance, of integer(pInt) :: & i - real(pReal), dimension(param(instance)%totalNslip) :: & - invLambdaSlip ! 1/mean free distance between 2 forest dislocations seen by a moving dislocation associate(prm => param(instance), stt => state(instance),dst => dependentState(instance)) forall (i = 1_pInt:prm%totalNslip) - invLambdaSlip(i) = sqrt(dot_product(stt%rhoEdge(:,of)+stt%rhoEdgeDip(:,of), & - prm%forestProjectionEdge(:,i))) & - / prm%Clambda(i) + dst%dislocationSpacing(i,of) = sqrt(dot_product(stt%rhoEdge(:,of)+stt%rhoEdgeDip(:,of), & + prm%forestProjectionEdge(:,i))) dst%threshold_stress(i,of) = prm%mu*prm%burgers(i) & * sqrt(dot_product(stt%rhoEdge(:,of)+stt%rhoEdgeDip(:,of), & prm%interaction_SlipSlip(i,:))) end forall - dst%mfp(:,of) = prm%grainSize/(1.0_pReal+prm%grainSize*invLambdaSlip) + dst%mfp(:,of) = prm%grainSize/(1.0_pReal+prm%grainSize*dst%dislocationSpacing(:,of)/prm%Clambda) + dst%dislocationSpacing(:,of) = dst%mfp(:,of) ! ToDo: Hack to recover wrong behavior for the moment + end associate @@ -434,24 +429,31 @@ end subroutine plastic_disloUCLA_dependentState pure subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp,Mp,Temperature,instance,of) implicit none - integer(pInt), intent(in) :: instance, of - real(pReal), intent(in) :: Temperature - real(pReal), dimension(3,3), intent(in) :: Mp - real(pReal), dimension(3,3), intent(out) :: Lp - real(pReal), dimension(3,3,3,3), intent(out) :: dLp_dMp + real(pReal), dimension(3,3), intent(out) :: & + Lp !< plastic velocity gradient + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp !< derivative of Lp with respect to the Mandel stress - integer(pInt) :: i,k,l,m,n + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + temperature !< temperature + integer(pInt), intent(in) :: & + instance, & + of + integer(pInt) :: & + i,k,l,m,n real(pReal), dimension(param(instance)%totalNslip) :: & - gdot_slip_pos,gdot_slip_neg,tau_slip_pos,tau_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg - - associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) + dgdot_dtauslip_pos,dgdot_dtauslip_neg, & + gdot_slip_pos,gdot_slip_neg Lp = 0.0_pReal dLp_dMp = 0.0_pReal + + associate(prm => param(instance)) - call kinetics(prm,stt,dst,Mp,Temperature,of, & - gdot_slip_pos,dgdot_dtauslip_pos,tau_slip_pos,gdot_slip_neg,dgdot_dtauslip_neg,tau_slip_neg) + call kinetics(Mp,Temperature,instance,of,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg) slipSystems: do i = 1_pInt, prm%totalNslip Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%Schmid_slip(1:3,1:3,i) forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & @@ -459,11 +461,9 @@ pure subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp,Mp,Temperature,inst + dgdot_dtauslip_pos(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_pos(m,n,i) & + dgdot_dtauslip_neg(i) * prm%Schmid_slip(k,l,i) * prm%nonSchmid_neg(m,n,i) enddo slipSystems + end associate - Lp = 0.5_pReal * Lp - dLp_dMp = 0.5_pReal * dLp_dMp - end subroutine plastic_disloUCLA_LpAndItsTangent @@ -480,9 +480,9 @@ subroutine plastic_disloUCLA_dotState(Mp,Temperature,instance,of) implicit none real(pReal), dimension(3,3), intent(in):: & - Mp !< 2nd Piola Kirchhoff stress tensor in Mandel notation + Mp !< Mandel stress real(pReal), intent(in) :: & - temperature !< temperature at integration point + temperature !< temperature integer(pInt), intent(in) :: & instance, of @@ -492,44 +492,41 @@ subroutine plastic_disloUCLA_dotState(Mp,Temperature,instance,of) gdot_slip_pos, gdot_slip_neg,& tau_slip_pos,& tau_slip_neg, & - dgdot_dtauslip_neg,dgdot_dtauslip_pos,DotRhoDipFormation, ClimbVelocity, EdgeDipDistance, & + DotRhoDipFormation, ClimbVelocity, EdgeDipDistance, & DotRhoEdgeDipClimb associate(prm => param(instance), stt => state(instance),dot => dotState(instance), dst => dependentState(instance)) - call kinetics(prm,stt,dst,Mp,Temperature,of, & - gdot_slip_pos,dgdot_dtauslip_pos,tau_slip_pos,gdot_slip_neg,dgdot_dtauslip_neg,tau_slip_neg) - - dot%whole(:,of) = 0.0_pReal - dot%accshear_slip(:,of) = (gdot_slip_pos+gdot_slip_neg)*0.5_pReal + call kinetics(Mp,Temperature,instance,of,& + gdot_slip_pos,gdot_slip_neg, & + tau_slip_pos1 = tau_slip_pos,tau_slip_neg1 = tau_slip_neg) + dot%accshear_slip(:,of) = (gdot_slip_pos+gdot_slip_neg) ! ToDo: needs to be abs VacancyDiffusion = prm%D0*exp(-prm%Qsd/(kB*Temperature)) - where(dEq0(tau_slip_pos)) - EdgeDipDistance = dst%mfp(:,of) !ToDo MD@FR: correct? was not handled properly before + where(dEq0(tau_slip_pos)) ! ToDo: use avg of pos and neg DotRhoDipFormation = 0.0_pReal DotRhoEdgeDipClimb = 0.0_pReal else where EdgeDipDistance = math_clip((3.0_pReal*prm%mu*prm%burgers)/(16.0_pReal*PI*abs(tau_slip_pos)), & - prm%minDipDistance, & ! lower limit - dst%mfp(:,of)) ! upper limit - DotRhoDipFormation = merge(((2.0_pReal*EdgeDipDistance)/prm%burgers)* stt%rhoEdge(:,of)*abs(dot%accshear_slip(:,of)), & + prm%minDipDistance, & ! lower limit + dst%mfp(:,of)) ! upper limit + DotRhoDipFormation = merge(((2.0_pReal*EdgeDipDistance)/prm%burgers)* stt%rhoEdge(:,of)*abs(dot%accshear_slip(:,of)), & ! ToDo: ignore region of spontaneous annihilation 0.0_pReal, & prm%dipoleformation) ClimbVelocity = (3.0_pReal*prm%mu*VacancyDiffusion*prm%atomicVolume/(2.0_pReal*pi*kB*Temperature)) & * (1.0_pReal/(EdgeDipDistance+prm%minDipDistance)) - DotRhoEdgeDipClimb = (4.0_pReal*ClimbVelocity*stt%rhoEdgeDip(:,of))/(EdgeDipDistance-prm%minDipDistance) + DotRhoEdgeDipClimb = (4.0_pReal*ClimbVelocity*stt%rhoEdgeDip(:,of))/(EdgeDipDistance-prm%minDipDistance) ! ToDo: Discuss with Franz: Stress dependency? end where - dot%rhoEdge(:,of) = abs(dot%accshear_slip(:,of))/(prm%burgers*dst%mfp(:,of)) & ! multiplication + dot%rhoEdge(:,of) = abs(dot%accshear_slip(:,of))/(prm%burgers*dst%mfp(:,of)) & ! multiplication - DotRhoDipFormation & - (2.0_pReal*prm%minDipDistance)/prm%burgers*stt%rhoEdge(:,of)*abs(dot%accshear_slip(:,of)) !* Spontaneous annihilation of 2 single edge dislocations - dot%rhoEdgeDip(:,of) = DotRhoDipFormation & - - (2.0_pReal*prm%minDipDistance)/prm%burgers* stt%rhoEdgeDip(:,of)*abs(dot%accshear_slip(:,of)) & !* Spontaneous annihilation of a single edge dislocation with a dipole constituent + - (2.0_pReal*prm%minDipDistance)/prm%burgers*stt%rhoEdgeDip(:,of)*abs(dot%accshear_slip(:,of)) & !* Spontaneous annihilation of a single edge dislocation with a dipole constituent - DotRhoEdgeDipClimb -end associate + end associate end subroutine plastic_disloUCLA_dotState @@ -541,14 +538,14 @@ function plastic_disloUCLA_postResults(Mp,Temperature,instance,of) result(postRe use prec, only: & dEq, dNeq0 use math, only: & - pi, & + PI, & math_mul33xx33 implicit none real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress real(pReal), intent(in) :: & - Temperature !< Mandel stress + temperature !< temperature integer(pInt), intent(in) :: & instance, & of @@ -559,14 +556,12 @@ function plastic_disloUCLA_postResults(Mp,Temperature,instance,of) result(postRe integer(pInt) :: & o,c,i real(pReal), dimension(param(instance)%totalNslip) :: & - gdot_slip_pos,dgdot_dtauslip_pos,tau_slip_pos, & - gdot_slip_neg,dgdot_dtauslip_neg,tau_slip_neg + gdot_slip_pos,gdot_slip_neg + + c = 0_pInt associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) - postResults = 0.0_pReal - c = 0_pInt - outputsLoop: do o = 1_pInt,size(prm%outputID) select case(prm%outputID(o)) @@ -574,32 +569,16 @@ function plastic_disloUCLA_postResults(Mp,Temperature,instance,of) result(postRe postResults(c+1_pInt:c+prm%totalNslip) = stt%rhoEdge(1_pInt:prm%totalNslip,of) case (rhoDip_ID) postResults(c+1_pInt:c+prm%totalNslip) = stt%rhoEdgeDip(1_pInt:prm%totalNslip,of) - case (shearrate_ID,stressexponent_ID) - call kinetics(prm,stt,dst,Mp,Temperature,of, & - gdot_slip_pos,dgdot_dtauslip_pos,tau_slip_pos,gdot_slip_neg,dgdot_dtauslip_neg,tau_slip_neg) - - if (prm%outputID(o) == shearrate_ID) then - postResults(c+1:c+prm%totalNslip) = (gdot_slip_pos + gdot_slip_neg)*0.5_pReal - elseif(prm%outputID(o) == stressexponent_ID) then - where (dNeq0(gdot_slip_pos+gdot_slip_neg)) - postResults(c+1_pInt:c + prm%totalNslip) = (tau_slip_pos+tau_slip_neg) * 0.5_pReal & - / (gdot_slip_pos+gdot_slip_neg) & - * (dgdot_dtauslip_pos+dgdot_dtauslip_neg) - else where - postResults(c+1_pInt:c + prm%totalNslip) = 0.0_pReal - end where - endif + case (shearrate_ID) + call kinetics(Mp,Temperature,instance,of,gdot_slip_pos,gdot_slip_neg) + postResults(c+1:c+prm%totalNslip) = gdot_slip_pos + gdot_slip_neg case (accumulatedshear_ID) postResults(c+1_pInt:c+prm%totalNslip) = stt%accshear_slip(1_pInt:prm%totalNslip, of) case (mfp_ID) postResults(c+1_pInt:c+prm%totalNslip) = dst%mfp(1_pInt:prm%totalNslip, of) - case (resolvedstress_ID) - do i = 1_pInt, prm%totalNslip - postResults(c+i) =math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - enddo case (thresholdstress_ID) postResults(c+1_pInt:c+prm%totalNslip) = dst%threshold_stress(1_pInt:prm%totalNslip,of) - case (dipoleDistance_ID) + case (dipoleDistance_ID) ! ToDo: Discuss required changes with Franz do i = 1_pInt, prm%totalNslip if (dNeq0(abs(math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,i))))) then postResults(c+i) = (3.0_pReal*prm%mu*prm%burgers(i)) & @@ -609,183 +588,161 @@ function plastic_disloUCLA_postResults(Mp,Temperature,instance,of) result(postRe endif postResults(c+i)=min(postResults(c+i),dst%mfp(i,of)) enddo + end select c = c + prm%totalNslip + enddo outputsLoop + end associate end function plastic_disloUCLA_postResults +!-------------------------------------------------------------------------------------------------- +!> @brief Shear rates on slip systems, their derivatives with respect to resolved stress and the +! resolved stresss +!> @details Derivatives and resolved stress are calculated only optionally. +! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to +! have the optional arguments at the end !-------------------------------------------------------------------------------------------------- -!> @brief return array of constitutive results -!-------------------------------------------------------------------------------------------------- -pure subroutine kinetics(prm,stt,dst,Mp,Temperature,of, & - gdot_slip_pos,dgdot_dtauslip_pos,tau_slip_pos,gdot_slip_neg,dgdot_dtauslip_neg,tau_slip_neg) +pure subroutine kinetics(Mp,Temperature,instance,of, & + gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg,tau_slip_pos1,tau_slip_neg1) use prec, only: & tol_math_check, & dEq, dNeq0 use math, only: & - pi, & -math_mul33xx33 + PI, & + math_mul33xx33 implicit none - type(tParameters), intent(in) :: & - prm - type(tDisloUCLAState), intent(in) :: & - stt - type(tDisloUCLAdependentState), intent(in) :: & - dst real(pReal), dimension(3,3), intent(in) :: & - Mp !< 2nd Piola Kirchhoff stress tensor in Mandel notation - real(pReal), intent(in) :: & - temperature !< temperature at integration point - integer(pInt), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + temperature !< temperature + integer(pInt), intent(in) :: & + instance, & of - integer(pInt) :: & - j - real(pReal), intent(out), dimension(prm%totalNslip) :: & - gdot_slip_pos,dgdot_dtauslip_pos,tau_slip_pos,gdot_slip_neg,dgdot_dtauslip_neg,tau_slip_neg - real(pReal), dimension(prm%totalNslip) :: & - StressRatio, BoltzmannRatio, & + real(pReal), intent(out), dimension(param(instance)%totalNslip) :: & + gdot_slip_pos, & + gdot_slip_neg + real(pReal), intent(out), optional, dimension(param(instance)%totalNslip) :: & + dgdot_dtauslip_pos, & + dgdot_dtauslip_neg, & + tau_slip_pos1, & + tau_slip_neg1 + real(pReal), dimension(param(instance)%totalNslip) :: & + StressRatio, & StressRatio_p,StressRatio_pminus1, & - DotGamma0, dvel_slip, vel_slip + dvel_slip, vel_slip, & + tau_slip_pos,tau_slip_neg, & + needsGoodName ! ToDo: @Karo: any idea? + integer(pInt) :: j + associate(prm => param(instance), stt => state(instance), dst => dependentState(instance)) + do j = 1_pInt, prm%totalNslip tau_slip_pos(j) = math_mul33xx33(Mp,prm%nonSchmid_pos(1:3,1:3,j)) tau_slip_neg(j) = math_mul33xx33(Mp,prm%nonSchmid_neg(1:3,1:3,j)) enddo + + + if (present(tau_slip_pos1)) tau_slip_pos1 = tau_slip_pos + if (present(tau_slip_neg1)) tau_slip_neg1 = tau_slip_neg - BoltzmannRatio = prm%H0kp/(kB*Temperature) - DotGamma0 = stt%rhoEdge(:,of)*prm%burgers*prm%v0 + associate(BoltzmannRatio => prm%H0kp/(kB*Temperature), & + DotGamma0 => stt%rhoEdge(:,of)*prm%burgers*prm%v0, & + effectiveLength => dst%mfp(:,of) - prm%w) - significantPositiveTau: where(abs(tau_slip_pos)-dst%threshold_stress(:,of) > tol_math_check) - StressRatio = (abs(tau_slip_pos)-dst%threshold_stress(:,of)) & - / (prm%solidSolutionStrength+prm%tau_Peierls) - StressRatio_p = StressRatio** prm%p - StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) + significantPositiveTau: where(abs(tau_slip_pos)-dst%threshold_stress(:,of) > tol_math_check) + StressRatio = (abs(tau_slip_pos)-dst%threshold_stress(:,of))/prm%tau0 + StressRatio_p = StressRatio** prm%p + StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) + needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) - vel_slip = 2.0_pReal*prm%burgers * prm%kink_height * prm%omega & - * ( dst%mfp(:,of) - prm%kink_width ) & - * (tau_slip_pos & - * exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) ) & - / ( & - 2.0_pReal*(prm%burgers**2.0_pReal)*tau_slip_pos & - + prm%omega * prm%B & - *(( dst%mfp(:,of) - prm%kink_width )**2.0_pReal) & - * exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) & + vel_slip = 2.0_pReal*prm%burgers * prm%kink_height * prm%omega & + * effectiveLength * tau_slip_pos * needsGoodName & + / ( 2.0_pReal*(prm%burgers**2.0_pReal)*tau_slip_pos & + + prm%omega * prm%B * effectiveLength**2.0_pReal* needsGoodName & ) - gdot_slip_pos = DotGamma0 * sign(vel_slip,tau_slip_pos) + gdot_slip_pos = DotGamma0 * sign(vel_slip,tau_slip_pos) * 0.5_pReal + else where significantPositiveTau + gdot_slip_pos = 0.0_pReal + end where significantPositiveTau - dvel_slip = 2.0_pReal*prm%burgers * prm%kink_height * prm%omega & - * ( dst%mfp(:,of) - prm%kink_width ) & - * ( & - (exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) & - + tau_slip_pos & - * (abs(exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q))& - *BoltzmannRatio*prm%p& - *prm%q/& - (prm%solidSolutionStrength+prm%tau_Peierls)*& - StressRatio_pminus1*(1-StressRatio_p)**(prm%q-1.0_pReal) ) & - ) & - * (2.0_pReal*(prm%burgers**2.0_pReal)*tau_slip_pos & - + prm%omega * prm%B & - *(( dst%mfp(:,of) - prm%kink_width )**2.0_pReal) & - * exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) & - ) & - - (tau_slip_pos & - * exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) ) & - * (2.0_pReal*(prm%burgers**2.0_pReal) & - + prm%omega * prm%B & - *(( dst%mfp(:,of) - prm%kink_width )**2.0_pReal) & - * (abs(exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q))& - *BoltzmannRatio*prm%p& - *prm%q/& - (prm%solidSolutionStrength+prm%tau_Peierls)*& - StressRatio_pminus1*(1-StressRatio_p)**(prm%q-1.0_pReal) )& - ) & - ) & - / ( & - ( & - 2.0_pReal*(prm%burgers**2.0_pReal)*tau_slip_pos & - + prm%omega * prm%B & - *(( dst%mfp(:,of) - prm%kink_width )**2.0_pReal) & - * exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) & - )**2.0_pReal & - ) + if (present(dgdot_dtauslip_pos)) then + significantPositiveTau2: where(abs(tau_slip_pos)-dst%threshold_stress(:,of) > tol_math_check) + dvel_slip = 2.0_pReal*prm%burgers * prm%kink_height * prm%omega* effectiveLength & + * ( & + (needsGoodName + tau_slip_pos * abs(needsGoodName)*BoltzmannRatio*prm%p & + * prm%q/prm%tau0 & + * StressRatio_pminus1*(1-StressRatio_p)**(prm%q-1.0_pReal) & + ) & + * ( 2.0_pReal*(prm%burgers**2.0_pReal)*tau_slip_pos & + + prm%omega * prm%B* effectiveLength **2.0_pReal* needsGoodName & + ) & + - tau_slip_pos * needsGoodName * (2.0_pReal*prm%burgers**2.0_pReal & + + prm%omega * prm%B *effectiveLength **2.0_pReal& + * (abs(needsGoodName)*BoltzmannRatio*prm%p *prm%q/prm%tau0 & + *StressRatio_pminus1*(1-StressRatio_p)**(prm%q-1.0_pReal) )& + ) & + ) & + /(2.0_pReal*prm%burgers**2.0_pReal*tau_slip_pos & + + prm%omega * prm%B* effectiveLength**2.0_pReal* needsGoodName )**2.0_pReal - dgdot_dtauslip_pos = DotGamma0 * dvel_slip -else where significantPositiveTau - gdot_slip_pos = 0.0_pReal - dgdot_dtauslip_pos = 0.0_pReal -end where significantPositiveTau + dgdot_dtauslip_pos = DotGamma0 * dvel_slip* 0.5_pReal + else where significantPositiveTau2 + dgdot_dtauslip_pos = 0.0_pReal + end where significantPositiveTau2 + endif - significantNegativeTau: where(abs(tau_slip_neg)-dst%threshold_stress(:,of) > tol_math_check) - StressRatio = (abs(tau_slip_neg)-dst%threshold_stress(:,of)) & - / (prm%solidSolutionStrength+prm%tau_Peierls) - StressRatio_p = StressRatio** prm%p - StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) + significantNegativeTau: where(abs(tau_slip_neg)-dst%threshold_stress(:,of) > tol_math_check) + StressRatio = (abs(tau_slip_neg)-dst%threshold_stress(:,of))/prm%tau0 + StressRatio_p = StressRatio** prm%p + StressRatio_pminus1 = StressRatio**(prm%p-1.0_pReal) + needsGoodName = exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) - vel_slip = 2.0_pReal*prm%burgers * prm%kink_height * prm%omega & - * ( dst%mfp(:,of) - prm%kink_width ) & - * (tau_slip_neg & - * exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) ) & - / ( & - 2.0_pReal*(prm%burgers**2.0_pReal)*tau_slip_neg & - + prm%omega * prm%B & - *(( dst%mfp(:,of) - prm%kink_width )**2.0_pReal) & - * exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) & + vel_slip = 2.0_pReal*prm%burgers * prm%kink_height * prm%omega & + * effectiveLength * tau_slip_neg * needsGoodName & + / ( 2.0_pReal*(prm%burgers**2.0_pReal)*tau_slip_neg & + + prm%omega * prm%B * effectiveLength**2.0_pReal* needsGoodName & ) - gdot_slip_neg = DotGamma0 * sign(vel_slip,tau_slip_neg) - - dvel_slip = 2.0_pReal*prm%burgers * prm%kink_height * prm%omega & - * ( dst%mfp(:,of) - prm%kink_width ) & - * ( & - (exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) & - + tau_slip_neg & - * (abs(exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q))& - *BoltzmannRatio*prm%p& - *prm%q/& - (prm%solidSolutionStrength+prm%tau_Peierls)*& - StressRatio_pminus1*(1-StressRatio_p)**(prm%q-1.0_pReal) ) & - ) & - * (2.0_pReal*(prm%burgers**2.0_pReal)*tau_slip_neg & - + prm%omega * prm%B & - *(( dst%mfp(:,of) - prm%kink_width )**2.0_pReal) & - * exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) & - ) & - - (tau_slip_neg & - * exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) ) & - * (2.0_pReal*(prm%burgers**2.0_pReal) & - + prm%omega * prm%B & - *(( dst%mfp(:,of) - prm%kink_width )**2.0_pReal) & - * (abs(exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q))& - *BoltzmannRatio*prm%p& - *prm%q/& - (prm%solidSolutionStrength+prm%tau_Peierls)*& - StressRatio_pminus1*(1-StressRatio_p)**(prm%q-1.0_pReal) )& - ) & - ) & - / ( & - ( & - 2.0_pReal*(prm%burgers**2.0_pReal)*tau_slip_neg & - + prm%omega * prm%B & - *(( dst%mfp(:,of) - prm%kink_width )**2.0_pReal) & - * exp(-BoltzmannRatio*(1-StressRatio_p) ** prm%q) & - )**2.0_pReal & - ) - - - dgdot_dtauslip_neg = DotGamma0 * dvel_slip -else where significantNegativeTau - gdot_slip_neg = 0.0_pReal - dgdot_dtauslip_neg = 0.0_pReal -end where significantNegativeTau + gdot_slip_neg = DotGamma0 * sign(vel_slip,tau_slip_neg) * 0.5_pReal + else where significantNegativeTau + gdot_slip_neg = 0.0_pReal + end where significantNegativeTau + if (present(dgdot_dtauslip_neg)) then + significantNegativeTau2: where(abs(tau_slip_neg)-dst%threshold_stress(:,of) > tol_math_check) + dvel_slip = 2.0_pReal*prm%burgers * prm%kink_height * prm%omega* effectiveLength & + * ( & + (needsGoodName + tau_slip_neg * abs(needsGoodName)*BoltzmannRatio*prm%p & + * prm%q/prm%tau0 & + * StressRatio_pminus1*(1-StressRatio_p)**(prm%q-1.0_pReal) & + ) & + * ( 2.0_pReal*(prm%burgers**2.0_pReal)*tau_slip_neg & + + prm%omega * prm%B* effectiveLength **2.0_pReal* needsGoodName & + ) & + - tau_slip_neg * needsGoodName * (2.0_pReal*prm%burgers**2.0_pReal & + + prm%omega * prm%B *effectiveLength **2.0_pReal& + * (abs(needsGoodName)*BoltzmannRatio*prm%p *prm%q/prm%tau0 & + *StressRatio_pminus1*(1-StressRatio_p)**(prm%q-1.0_pReal) )& + ) & + ) & + /(2.0_pReal*prm%burgers**2.0_pReal*tau_slip_neg & + + prm%omega * prm%B* effectiveLength**2.0_pReal* needsGoodName )**2.0_pReal + + dgdot_dtauslip_neg = DotGamma0 * dvel_slip * 0.5_pReal + else where significantNegativeTau2 + dgdot_dtauslip_neg = 0.0_pReal + end where significantNegativeTau2 + end if + end associate + end associate end subroutine kinetics diff --git a/src/plastic_isotropic.f90 b/src/plastic_isotropic.f90 index a44681676..c7d92651a 100644 --- a/src/plastic_isotropic.f90 +++ b/src/plastic_isotropic.f90 @@ -9,7 +9,7 @@ !-------------------------------------------------------------------------------------------------- module plastic_isotropic use prec, only: & - pReal,& + pReal, & pInt implicit none @@ -71,7 +71,6 @@ module plastic_isotropic contains - !-------------------------------------------------------------------------------------------------- !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks @@ -92,7 +91,7 @@ subroutine plastic_isotropic_init() debug_levelExtensive, & #endif debug_level, & - debug_constitutive,& + debug_constitutive, & debug_levelBasic use IO, only: & IO_error, & @@ -122,7 +121,7 @@ subroutine plastic_isotropic_init() sizeState, sizeDotState character(len=65536), dimension(0), parameter :: emptyStringArray = [character(len=65536)::] - + integer(kind(undefined_ID)) :: & outputID @@ -149,7 +148,7 @@ subroutine plastic_isotropic_init() allocate(state(Ninstance)) allocate(dotState(Ninstance)) - do p = 1_pInt, size(phase_plasticityInstance) + do p = 1_pInt, size(phase_plasticity) if (phase_plasticity(p) /= PLASTICITY_ISOTROPIC_ID) cycle associate(prm => param(phase_plasticityInstance(p)), & dot => dotState(phase_plasticityInstance(p)), & @@ -204,25 +203,27 @@ subroutine plastic_isotropic_init() do i=1_pInt, size(outputs) outputID = undefined_ID select case(outputs(i)) + case ('flowstress') outputID = flowstress_ID case ('strainrate') outputID = strainrate_ID + end select - + if (outputID /= undefined_ID) then - plastic_isotropic_output(i,phase_plasticityInstance(p)) = outputs(i) - plastic_isotropic_sizePostResult(i,phase_plasticityInstance(p)) = 1_pInt - prm%outputID = [prm%outputID , outputID] + plastic_isotropic_output(i,phase_plasticityInstance(p)) = outputs(i) + plastic_isotropic_sizePostResult(i,phase_plasticityInstance(p)) = 1_pInt + prm%outputID = [prm%outputID, outputID] endif - end do + enddo !-------------------------------------------------------------------------------------------------- ! allocate state arrays NipcMyPhase = count(material_phase == p) - sizeState = size(["flowstress ","accumulated_shear"]) - sizeDotState = sizeState + sizeDotState = size(['flowstress ','accumulated_shear']) + sizeState = sizeDotState call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0_pInt, & 1_pInt,0_pInt,0_pInt) @@ -243,8 +244,8 @@ subroutine plastic_isotropic_init() plasticState(p)%accumulatedSlip => plasticState(p)%state (2:2,1:NipcMyPhase) plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally - - end associate + + end associate enddo @@ -319,7 +320,7 @@ subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) end if end associate - + end subroutine plastic_isotropic_LpAndItsTangent @@ -373,9 +374,9 @@ subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar,Tstar,instance,of) Li = 0.0_pReal dLi_dTstar = 0.0_pReal endif - + end associate - + end subroutine plastic_isotropic_LiAndItsTangent @@ -471,6 +472,7 @@ function plastic_isotropic_postResults(Mp,instance,of) result(postResults) outputsLoop: do o = 1_pInt,size(prm%outputID) select case(prm%outputID(o)) + case (flowstress_ID) postResults(c+1_pInt) = stt%flowstress(of) c = c + 1_pInt @@ -478,6 +480,7 @@ function plastic_isotropic_postResults(Mp,instance,of) result(postResults) postResults(c+1_pInt) = prm%gdot0 & * (sqrt(1.5_pReal) * norm_Mp /(prm%fTaylor * stt%flowstress(of)))**prm%n c = c + 1_pInt + end select enddo outputsLoop diff --git a/src/plastic_none.f90 b/src/plastic_none.f90 index 5b5bb49d1..2c6ca6e93 100644 --- a/src/plastic_none.f90 +++ b/src/plastic_none.f90 @@ -1,7 +1,8 @@ !-------------------------------------------------------------------------------------------------- !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine for purely elastic material +!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH +!> @brief Dummy plasticity for purely elastic material !-------------------------------------------------------------------------------------------------- module plastic_none @@ -13,7 +14,6 @@ module plastic_none contains - !-------------------------------------------------------------------------------------------------- !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks @@ -32,52 +32,40 @@ subroutine plastic_none_init debug_levelBasic use IO, only: & IO_timeStamp - use numerics, only: & - numerics_integrator use material, only: & phase_plasticity, & + material_allocatePlasticState, & PLASTICITY_NONE_label, & + PLASTICITY_NONE_ID, & material_phase, & - plasticState, & - PLASTICITY_none_ID + plasticState implicit none - integer(pInt) :: & - maxNinstance, & - phase, & - NofMyPhase - + Ninstance, & + p, & + NipcMyPhase + write(6,'(/,a)') ' <<<+- plastic_'//PLASTICITY_NONE_label//' init -+>>>' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" - - maxNinstance = int(count(phase_plasticity == PLASTICITY_none_ID),pInt) + + Ninstance = int(count(phase_plasticity == PLASTICITY_NONE_ID),pInt) if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & - write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance + write(6,'(a16,1x,i5,/)') '# instances:',Ninstance - initializeInstances: do phase = 1_pInt, size(phase_plasticity) - if (phase_plasticity(phase) == PLASTICITY_none_ID) then - NofMyPhase=count(material_phase==phase) + do p = 1_pInt, size(phase_plasticity) + if (phase_plasticity(p) /= PLASTICITY_NONE_ID) cycle - allocate(plasticState(phase)%aTolState (0_pInt)) - allocate(plasticState(phase)%state0 (0_pInt,NofMyPhase)) - allocate(plasticState(phase)%partionedState0 (0_pInt,NofMyPhase)) - allocate(plasticState(phase)%subState0 (0_pInt,NofMyPhase)) - allocate(plasticState(phase)%state (0_pInt,NofMyPhase)) +!-------------------------------------------------------------------------------------------------- +! allocate state arrays + NipcMyPhase = count(material_phase == p) - allocate(plasticState(phase)%dotState (0_pInt,NofMyPhase)) - allocate(plasticState(phase)%deltaState (0_pInt,NofMyPhase)) - if (any(numerics_integrator == 1_pInt)) then - allocate(plasticState(phase)%previousDotState (0_pInt,NofMyPhase)) - allocate(plasticState(phase)%previousDotState2(0_pInt,NofMyPhase)) - endif - if (any(numerics_integrator == 4_pInt)) & - allocate(plasticState(phase)%RK4dotState (0_pInt,NofMyPhase)) - if (any(numerics_integrator == 5_pInt)) & - allocate(plasticState(phase)%RKCK45dotState (6,0_pInt,NofMyPhase)) - endif - enddo initializeInstances + call material_allocatePlasticState(p,NipcMyPhase,0_pInt,0_pInt,0_pInt, & + 0_pInt,0_pInt,0_pInt) + plasticState(p)%sizePostResults = 0_pInt + + enddo end subroutine plastic_none_init diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index 16aac1ead..e2b56cce6 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -2,12 +2,11 @@ !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!> @brief material subroutine for phenomenological crystal plasticity formulation using a powerlaw -!! fitting +!> @brief phenomenological crystal plasticity formulation using a powerlaw fitting !-------------------------------------------------------------------------------------------------- module plastic_phenopowerlaw use prec, only: & - pReal,& + pReal, & pInt implicit none @@ -101,7 +100,6 @@ module plastic_phenopowerlaw contains - !-------------------------------------------------------------------------------------------------- !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks @@ -174,7 +172,7 @@ subroutine plastic_phenopowerlaw_init allocate(state(Ninstance)) allocate(dotState(Ninstance)) - do p = 1_pInt, size(phase_plasticityInstance) + do p = 1_pInt, size(phase_plasticity) if (phase_plasticity(p) /= PLASTICITY_PHENOPOWERLAW_ID) cycle associate(prm => param(phase_plasticityInstance(p)), & dot => dotState(phase_plasticityInstance(p)), & @@ -195,9 +193,9 @@ subroutine plastic_phenopowerlaw_init prm%aTolTwinfrac = config%getFloat('atol_twinfrac', defaultVal=1.0e-6_pReal) ! sanity checks - if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//'aTolresistance ' - if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//'aTolShear ' - if (prm%aTolTwinfrac <= 0.0_pReal) extmsg = trim(extmsg)//'atoltwinfrac ' + if (prm%aTolResistance <= 0.0_pReal) extmsg = trim(extmsg)//' aTolresistance' + if (prm%aTolShear <= 0.0_pReal) extmsg = trim(extmsg)//' aTolShear' + if (prm%aTolTwinfrac <= 0.0_pReal) extmsg = trim(extmsg)//' atoltwinfrac' !-------------------------------------------------------------------------------------------------- ! slip related parameters @@ -235,11 +233,11 @@ subroutine plastic_phenopowerlaw_init prm%H_int = math_expand(prm%H_int, prm%Nslip) ! sanity checks - if (prm%gdot0_slip <= 0.0_pReal) extmsg = trim(extmsg)//'gdot0_slip ' - if (prm%a_slip <= 0.0_pReal) extmsg = trim(extmsg)//'a_slip ' - if (prm%n_slip <= 0.0_pReal) extmsg = trim(extmsg)//'n_slip ' - if (any(prm%xi_slip_0 <= 0.0_pReal)) extmsg = trim(extmsg)//'xi_slip_0 ' - if (any(prm%xi_slip_sat < prm%xi_slip_0)) extmsg = trim(extmsg)//'xi_slip_sat ' + if (prm%gdot0_slip <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0_slip' + if (prm%a_slip <= 0.0_pReal) extmsg = trim(extmsg)//' a_slip' + if (prm%n_slip <= 0.0_pReal) extmsg = trim(extmsg)//' n_slip' + if (any(prm%xi_slip_0 <= 0.0_pReal)) extmsg = trim(extmsg)//' xi_slip_0' + if (any(prm%xi_slip_sat < prm%xi_slip_0)) extmsg = trim(extmsg)//' xi_slip_sat' else slipActive allocate(prm%interaction_SlipSlip(0,0)) allocate(prm%xi_slip_0(0)) @@ -269,8 +267,8 @@ subroutine plastic_phenopowerlaw_init prm%xi_twin_0 = math_expand(prm%xi_twin_0, prm%Ntwin) ! sanity checks - if (prm%gdot0_twin <= 0.0_pReal) extmsg = trim(extmsg)//'gdot0_twin ' - if (prm%n_twin <= 0.0_pReal) extmsg = trim(extmsg)//'n_twin ' + if (prm%gdot0_twin <= 0.0_pReal) extmsg = trim(extmsg)//' gdot0_twin' + if (prm%n_twin <= 0.0_pReal) extmsg = trim(extmsg)//' n_twin' else twinActive allocate(prm%interaction_TwinTwin(0,0)) allocate(prm%xi_twin_0(0)) @@ -304,48 +302,49 @@ subroutine plastic_phenopowerlaw_init do i=1_pInt, size(outputs) outputID = undefined_ID select case(outputs(i)) - case ('resistance_slip') - outputID = merge(resistance_slip_ID,undefined_ID,prm%totalNslip>0_pInt) - outputSize = prm%totalNslip - case ('accumulatedshear_slip') - outputID = merge(accumulatedshear_slip_ID,undefined_ID,prm%totalNslip>0_pInt) - outputSize = prm%totalNslip - case ('shearrate_slip') - outputID = merge(shearrate_slip_ID,undefined_ID,prm%totalNslip>0_pInt) - outputSize = prm%totalNslip - case ('resolvedstress_slip') - outputID = merge(resolvedstress_slip_ID,undefined_ID,prm%totalNslip>0_pInt) - outputSize = prm%totalNslip - case ('resistance_twin') - outputID = merge(resistance_twin_ID,undefined_ID,prm%totalNtwin>0_pInt) - outputSize = prm%totalNtwin - case ('accumulatedshear_twin') - outputID = merge(accumulatedshear_twin_ID,undefined_ID,prm%totalNtwin>0_pInt) - outputSize = prm%totalNtwin - case ('shearrate_twin') - outputID = merge(shearrate_twin_ID,undefined_ID,prm%totalNtwin>0_pInt) - outputSize = prm%totalNtwin - case ('resolvedstress_twin') - outputID = merge(resolvedstress_twin_ID,undefined_ID,prm%totalNtwin>0_pInt) - outputSize = prm%totalNtwin + case ('resistance_slip') + outputID = merge(resistance_slip_ID,undefined_ID,prm%totalNslip>0_pInt) + outputSize = prm%totalNslip + case ('accumulatedshear_slip') + outputID = merge(accumulatedshear_slip_ID,undefined_ID,prm%totalNslip>0_pInt) + outputSize = prm%totalNslip + case ('shearrate_slip') + outputID = merge(shearrate_slip_ID,undefined_ID,prm%totalNslip>0_pInt) + outputSize = prm%totalNslip + case ('resolvedstress_slip') + outputID = merge(resolvedstress_slip_ID,undefined_ID,prm%totalNslip>0_pInt) + outputSize = prm%totalNslip - end select + case ('resistance_twin') + outputID = merge(resistance_twin_ID,undefined_ID,prm%totalNtwin>0_pInt) + outputSize = prm%totalNtwin + case ('accumulatedshear_twin') + outputID = merge(accumulatedshear_twin_ID,undefined_ID,prm%totalNtwin>0_pInt) + outputSize = prm%totalNtwin + case ('shearrate_twin') + outputID = merge(shearrate_twin_ID,undefined_ID,prm%totalNtwin>0_pInt) + outputSize = prm%totalNtwin + case ('resolvedstress_twin') + outputID = merge(resolvedstress_twin_ID,undefined_ID,prm%totalNtwin>0_pInt) + outputSize = prm%totalNtwin - if (outputID /= undefined_ID) then - plastic_phenopowerlaw_output(i,phase_plasticityInstance(p)) = outputs(i) - plastic_phenopowerlaw_sizePostResult(i,phase_plasticityInstance(p)) = outputSize - prm%outputID = [prm%outputID , outputID] - endif + end select - end do + if (outputID /= undefined_ID) then + plastic_phenopowerlaw_output(i,phase_plasticityInstance(p)) = outputs(i) + plastic_phenopowerlaw_sizePostResult(i,phase_plasticityInstance(p)) = outputSize + prm%outputID = [prm%outputID, outputID] + endif + + enddo !-------------------------------------------------------------------------------------------------- ! allocate state arrays NipcMyPhase = count(material_phase == p) - sizeState = size(['tau_slip ','gamma_slip']) * prm%TotalNslip & - + size(['tau_twin ','gamma_twin']) * prm%TotalNtwin - sizeDotState = sizeState + sizeDotState = size(['tau_slip ','gamma_slip']) * prm%TotalNslip & + + size(['tau_twin ','gamma_twin']) * prm%TotalNtwin + sizeState = sizeDotState call material_allocatePlasticState(p,NipcMyPhase,sizeState,sizeDotState,0_pInt, & prm%totalNslip,prm%totalNtwin,0_pInt) @@ -385,6 +384,7 @@ subroutine plastic_phenopowerlaw_init plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally end associate + enddo end subroutine plastic_phenopowerlaw_init @@ -516,14 +516,14 @@ function plastic_phenopowerlaw_postResults(Mp,instance,of) result(postResults) math_mul33xx33 implicit none - real(pReal), dimension(3,3), intent(in) :: & + real(pReal), dimension(3,3), intent(in) :: & Mp !< Mandel stress - integer(pInt), intent(in) :: & + integer(pInt), intent(in) :: & instance, & of real(pReal), dimension(sum(plastic_phenopowerlaw_sizePostResult(:,instance))) :: & - postResults + postResults integer(pInt) :: & o,c,i @@ -531,8 +531,8 @@ function plastic_phenopowerlaw_postResults(Mp,instance,of) result(postResults) gdot_slip_pos,gdot_slip_neg c = 0_pInt - - associate( prm => param(instance), stt => state(instance)) + + associate(prm => param(instance), stt => state(instance)) outputsLoop: do o = 1_pInt,size(prm%outputID) select case(prm%outputID(o)) @@ -570,7 +570,7 @@ function plastic_phenopowerlaw_postResults(Mp,instance,of) result(postResults) end select enddo outputsLoop - + end associate end function plastic_phenopowerlaw_postResults @@ -585,7 +585,7 @@ end function plastic_phenopowerlaw_postResults pure subroutine kinetics_slip(Mp,instance,of, & gdot_slip_pos,gdot_slip_neg,dgdot_dtau_slip_pos,dgdot_dtau_slip_neg) use prec, only: & - dNeq0 + dNeq0 use math, only: & math_mul33xx33 @@ -596,13 +596,12 @@ pure subroutine kinetics_slip(Mp,instance,of, & instance, & of - real(pReal), dimension(param(instance)%totalNslip), intent(out) :: & + real(pReal), intent(out), dimension(param(instance)%totalNslip) :: & gdot_slip_pos, & gdot_slip_neg - real(pReal), dimension(param(instance)%totalNslip), intent(out), optional :: & + real(pReal), intent(out), optional, dimension(param(instance)%totalNslip) :: & dgdot_dtau_slip_pos, & dgdot_dtau_slip_neg - real(pReal), dimension(param(instance)%totalNslip) :: & tau_slip_pos, & tau_slip_neg