223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity
DAMASK_EICMD/code/constitutive_nonlocal.f90

2650 lines
147 KiB
Fortran
Raw Blame History

! Copyright 2011 Max-Planck-Institut für Eisenforschung GmbH
!
! This file is part of DAMASK,
! the Düsseldorf Advanced MAterial Simulation Kit.
!
! DAMASK is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
!
! DAMASK is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License
! along with DAMASK. If not, see <http://www.gnu.org/licenses/>.
!
!##############################################################
!* $Id$
!************************************
!* Module: CONSTITUTIVE_NONLOCAL *
!************************************
!* contains: *
!* - constitutive equations *
!* - parameters definition *
!************************************
MODULE constitutive_nonlocal
!* Include other modules
use prec, only: pReal,pInt
implicit none
!* Definition of parameters
character (len=*), parameter :: constitutive_nonlocal_label = 'nonlocal'
character(len=22), dimension(10), parameter :: constitutive_nonlocal_listBasicStates = (/'rhoSglEdgePosMobile ', &
'rhoSglEdgeNegMobile ', &
'rhoSglScrewPosMobile ', &
'rhoSglScrewNegMobile ', &
'rhoSglEdgePosImmobile ', &
'rhoSglEdgeNegImmobile ', &
'rhoSglScrewPosImmobile', &
'rhoSglScrewNegImmobile', &
'rhoDipEdge ', &
'rhoDipScrew ' /) ! list of "basic" microstructural state variables that are independent from other state variables
character(len=15), dimension(3), parameter :: constitutive_nonlocal_listDependentStates = (/'rhoForest ', &
'tauThreshold ', &
'Tdislocation_v ' /) ! list of microstructural state variables that depend on other state variables
character(len=15), dimension(2), parameter :: constitutive_nonlocal_listOtherStates = (/'velocityEdge ', &
'velocityScrew ' /) ! list of other dependent state variables that are not updated by microstructure
real(pReal), parameter :: kB = 1.38e-23_pReal ! Physical parameter, Boltzmann constant in J/Kelvin
!* Definition of global variables
integer(pInt), dimension(:), allocatable :: constitutive_nonlocal_sizeDotState, & ! number of dotStates = number of basic state variables
constitutive_nonlocal_sizeDependentState, & ! number of dependent state variables
constitutive_nonlocal_sizeState, & ! total number of state variables
constitutive_nonlocal_sizePostResults ! cumulative size of post results
integer(pInt), dimension(:,:), allocatable, target :: constitutive_nonlocal_sizePostResult ! size of each post result output
character(len=64), dimension(:,:), allocatable, target :: constitutive_nonlocal_output ! name of each post result output
character(len=32), dimension(:), allocatable :: constitutive_nonlocal_structureName ! name of the lattice structure
integer(pInt), dimension(:), allocatable :: constitutive_nonlocal_structure, & ! number representing the kind of lattice structure
constitutive_nonlocal_totalNslip ! total number of active slip systems for each instance
integer(pInt), dimension(:,:), allocatable :: constitutive_nonlocal_Nslip, & ! number of active slip systems for each family and instance
constitutive_nonlocal_slipFamily, & ! lookup table relating active slip system to slip family for each instance
constitutive_nonlocal_slipSystemLattice ! lookup table relating active slip system index to lattice slip system index for each instance
real(pReal), dimension(:), allocatable :: constitutive_nonlocal_CoverA, & ! c/a ratio for hex type lattice
constitutive_nonlocal_C11, & ! C11 element in elasticity matrix
constitutive_nonlocal_C12, & ! C12 element in elasticity matrix
constitutive_nonlocal_C13, & ! C13 element in elasticity matrix
constitutive_nonlocal_C33, & ! C33 element in elasticity matrix
constitutive_nonlocal_C44, & ! C44 element in elasticity matrix
constitutive_nonlocal_Gmod, & ! shear modulus
constitutive_nonlocal_nu, & ! poisson's ratio
constitutive_nonlocal_atomicVolume, & ! atomic volume
constitutive_nonlocal_Dsd0, & ! prefactor for self-diffusion coefficient
constitutive_nonlocal_Qsd, & ! activation enthalpy for diffusion
constitutive_nonlocal_aTolRho, & ! absolute tolerance for dislocation density in state integration
constitutive_nonlocal_R, & ! cutoff radius for dislocation stress
constitutive_nonlocal_d0, & ! wall depth as multiple of b
constitutive_nonlocal_tauObs, & ! obstacle strength in Pa
constitutive_nonlocal_fattack, & ! attack frequency in Hz
constitutive_nonlocal_vs, & ! maximum dislocation velocity = velocity of sound
constitutive_nonlocal_rhoSglScatter, & ! standard deviation of scatter in initial dislocation density
constitutive_nonlocal_surfaceTransmissivity ! transmissivity at free surface
real(pReal), dimension(:,:,:), allocatable :: constitutive_nonlocal_Cslip_66 ! elasticity matrix in Mandel notation for each instance
real(pReal), dimension(:,:,:,:,:), allocatable :: constitutive_nonlocal_Cslip_3333 ! elasticity matrix for each instance
real(pReal), dimension(:,:), allocatable :: constitutive_nonlocal_rhoSglEdgePos0, & ! initial edge_pos dislocation density per slip system for each family and instance
constitutive_nonlocal_rhoSglEdgeNeg0, & ! initial edge_neg dislocation density per slip system for each family and instance
constitutive_nonlocal_rhoSglScrewPos0, & ! initial screw_pos dislocation density per slip system for each family and instance
constitutive_nonlocal_rhoSglScrewNeg0, & ! initial screw_neg dislocation density per slip system for each family and instance
constitutive_nonlocal_rhoDipEdge0, & ! initial edge dipole dislocation density per slip system for each family and instance
constitutive_nonlocal_rhoDipScrew0, & ! initial screw dipole dislocation density per slip system for each family and instance
constitutive_nonlocal_lambda0PerSlipFamily, & ! mean free path prefactor for each family and instance
constitutive_nonlocal_lambda0PerSlipSystem, & ! mean free path prefactor for each slip system and instance
constitutive_nonlocal_burgersPerSlipFamily, & ! absolute length of burgers vector [m] for each family and instance
constitutive_nonlocal_burgersPerSlipSystem, & ! absolute length of burgers vector [m] for each slip system and instance
constitutive_nonlocal_dLowerEdgePerSlipFamily, & ! minimum stable edge dipole height for each family and instance
constitutive_nonlocal_dLowerEdgePerSlipSystem, & ! minimum stable edge dipole height for each slip system and instance
constitutive_nonlocal_dLowerScrewPerSlipFamily, & ! minimum stable screw dipole height for each family and instance
constitutive_nonlocal_dLowerScrewPerSlipSystem, & ! minimum stable screw dipole height for each slip system and instance
constitutive_nonlocal_Qeff0, & ! prefactor for activation enthalpy for dislocation glide in J
constitutive_nonlocal_interactionSlipSlip ! coefficients for slip-slip interaction for each interaction type and instance
real(pReal), dimension(:,:,:,:,:), allocatable :: constitutive_nonlocal_rhoDotFlux ! dislocation convection term
real(pReal), dimension(:,:,:,:,:,:), allocatable :: constitutive_nonlocal_compatibility ! slip system compatibility between me and my neighbors
real(pReal), dimension(:,:,:), allocatable :: constitutive_nonlocal_forestProjectionEdge, & ! matrix of forest projections of edge dislocations for each instance
constitutive_nonlocal_forestProjectionScrew, & ! matrix of forest projections of screw dislocations for each instance
constitutive_nonlocal_interactionMatrixSlipSlip ! interaction matrix of the different slip systems for each instance
real(pReal), dimension(:,:,:,:), allocatable :: constitutive_nonlocal_lattice2slip ! orthogonal transformation matrix from lattice coordinate system to slip coordinate system (passive rotation !!!)
CONTAINS
!****************************************
!* - constitutive_init
!* - constitutive_stateInit
!* - constitutive_homogenizedC
!* - constitutive_microstructure
!* - constitutive_LpAndItsTangent
!* - constitutive_dotState
!* - constitutive_dotTemperature
!* - constitutive_postResults
!****************************************
!**************************************
!* Module initialization *
!**************************************
subroutine constitutive_nonlocal_init(file)
use prec, only: pInt, pReal
use math, only: math_Mandel3333to66, &
math_Voigt66to3333, &
math_mul3x3, &
math_transpose3x3
use IO, only: IO_lc, &
IO_getTag, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_intValue, &
IO_error
use debug, only: debug_verbosity
use mesh, only: mesh_NcpElems, &
mesh_maxNips, &
FE_maxNipNeighbors
use material, only: homogenization_maxNgrains, &
phase_constitution, &
phase_constitutionInstance, &
phase_Noutput
use lattice, only: lattice_maxNslipFamily, &
lattice_maxNtwinFamily, &
lattice_maxNslip, &
lattice_maxNtwin, &
lattice_maxNinteraction, &
lattice_NslipSystem, &
lattice_NtwinSystem, &
lattice_initializeStructure, &
lattice_Qtwin, &
lattice_sd, &
lattice_sn, &
lattice_st, &
lattice_interactionSlipSlip
!*** output variables
!*** input variables
integer(pInt), intent(in) :: file
!*** local variables
integer(pInt), parameter :: maxNchunks = 21
integer(pInt), dimension(1+2*maxNchunks) :: positions
integer(pInt) section, &
maxNinstance, &
maxTotalNslip, &
myStructure, &
f, & ! index of my slip family
i, & ! index of my instance of this constitution
j, &
k, &
l, &
ns, & ! short notation for total number of active slip systems for the current instance
o, & ! index of my output
s, & ! index of my slip system
s1, & ! index of my slip system
s2, & ! index of my slip system
it, & ! index of my interaction type
output, &
mySize
character(len=64) tag
character(len=1024) line
!$OMP CRITICAL (write2out)
write(6,*)
write(6,*) '<<<+- constitutive_',trim(constitutive_nonlocal_label),' init -+>>>'
write(6,*) '$Id$'
write(6,*)
!$OMP END CRITICAL (write2out)
maxNinstance = count(phase_constitution == constitutive_nonlocal_label)
if (maxNinstance == 0) return ! we don't have to do anything if there's no instance for this constitutive law
if (debug_verbosity > 0) then
!$OMP CRITICAL (write2out)
write(6,'(a16,x,i5)') '# instances:',maxNinstance
!$OMP END CRITICAL (write2out)
endif
!*** space allocation for global variables
allocate(constitutive_nonlocal_sizeDotState(maxNinstance))
allocate(constitutive_nonlocal_sizeDependentState(maxNinstance))
allocate(constitutive_nonlocal_sizeState(maxNinstance))
allocate(constitutive_nonlocal_sizePostResults(maxNinstance))
allocate(constitutive_nonlocal_sizePostResult(maxval(phase_Noutput), maxNinstance))
allocate(constitutive_nonlocal_output(maxval(phase_Noutput), maxNinstance))
constitutive_nonlocal_sizeDotState = 0_pInt
constitutive_nonlocal_sizeDependentState = 0_pInt
constitutive_nonlocal_sizeState = 0_pInt
constitutive_nonlocal_sizePostResults = 0_pInt
constitutive_nonlocal_sizePostResult = 0_pInt
constitutive_nonlocal_output = ''
allocate(constitutive_nonlocal_structureName(maxNinstance))
allocate(constitutive_nonlocal_structure(maxNinstance))
allocate(constitutive_nonlocal_Nslip(lattice_maxNslipFamily, maxNinstance))
allocate(constitutive_nonlocal_slipFamily(lattice_maxNslip, maxNinstance))
allocate(constitutive_nonlocal_slipSystemLattice(lattice_maxNslip, maxNinstance))
allocate(constitutive_nonlocal_totalNslip(maxNinstance))
constitutive_nonlocal_structureName = ''
constitutive_nonlocal_structure = 0_pInt
constitutive_nonlocal_Nslip = 0_pInt
constitutive_nonlocal_slipFamily = 0_pInt
constitutive_nonlocal_slipSystemLattice = 0_pInt
constitutive_nonlocal_totalNslip = 0_pInt
allocate(constitutive_nonlocal_CoverA(maxNinstance))
allocate(constitutive_nonlocal_C11(maxNinstance))
allocate(constitutive_nonlocal_C12(maxNinstance))
allocate(constitutive_nonlocal_C13(maxNinstance))
allocate(constitutive_nonlocal_C33(maxNinstance))
allocate(constitutive_nonlocal_C44(maxNinstance))
allocate(constitutive_nonlocal_Gmod(maxNinstance))
allocate(constitutive_nonlocal_nu(maxNinstance))
allocate(constitutive_nonlocal_atomicVolume(maxNinstance))
allocate(constitutive_nonlocal_Dsd0(maxNinstance))
allocate(constitutive_nonlocal_Qsd(maxNinstance))
allocate(constitutive_nonlocal_aTolRho(maxNinstance))
allocate(constitutive_nonlocal_Cslip_66(6,6,maxNinstance))
allocate(constitutive_nonlocal_Cslip_3333(3,3,3,3,maxNinstance))
allocate(constitutive_nonlocal_R(maxNinstance))
allocate(constitutive_nonlocal_d0(maxNinstance))
allocate(constitutive_nonlocal_tauObs(maxNinstance))
allocate(constitutive_nonlocal_vs(maxNinstance))
allocate(constitutive_nonlocal_fattack(maxNinstance))
allocate(constitutive_nonlocal_rhoSglScatter(maxNinstance))
allocate(constitutive_nonlocal_surfaceTransmissivity(maxNinstance))
constitutive_nonlocal_CoverA = 0.0_pReal
constitutive_nonlocal_C11 = 0.0_pReal
constitutive_nonlocal_C12 = 0.0_pReal
constitutive_nonlocal_C13 = 0.0_pReal
constitutive_nonlocal_C33 = 0.0_pReal
constitutive_nonlocal_C44 = 0.0_pReal
constitutive_nonlocal_Gmod = 0.0_pReal
constitutive_nonlocal_atomicVolume = 0.0_pReal
constitutive_nonlocal_Dsd0 = 0.0_pReal
constitutive_nonlocal_Qsd = 0.0_pReal
constitutive_nonlocal_aTolRho = 0.0_pReal
constitutive_nonlocal_nu = 0.0_pReal
constitutive_nonlocal_Cslip_66 = 0.0_pReal
constitutive_nonlocal_Cslip_3333 = 0.0_pReal
constitutive_nonlocal_R = -1.0_pReal
constitutive_nonlocal_d0 = 0.0_pReal
constitutive_nonlocal_tauObs = 0.0_pReal
constitutive_nonlocal_vs = 0.0_pReal
constitutive_nonlocal_fattack = 0.0_pReal
constitutive_nonlocal_rhoSglScatter = 0.0_pReal
constitutive_nonlocal_surfaceTransmissivity = 1.0_pReal
allocate(constitutive_nonlocal_rhoSglEdgePos0(lattice_maxNslipFamily, maxNinstance))
allocate(constitutive_nonlocal_rhoSglEdgeNeg0(lattice_maxNslipFamily, maxNinstance))
allocate(constitutive_nonlocal_rhoSglScrewPos0(lattice_maxNslipFamily, maxNinstance))
allocate(constitutive_nonlocal_rhoSglScrewNeg0(lattice_maxNslipFamily, maxNinstance))
allocate(constitutive_nonlocal_rhoDipEdge0(lattice_maxNslipFamily, maxNinstance))
allocate(constitutive_nonlocal_rhoDipScrew0(lattice_maxNslipFamily, maxNinstance))
allocate(constitutive_nonlocal_burgersPerSlipFamily(lattice_maxNslipFamily, maxNinstance))
allocate(constitutive_nonlocal_Lambda0PerSlipFamily(lattice_maxNslipFamily, maxNinstance))
allocate(constitutive_nonlocal_interactionSlipSlip(lattice_maxNinteraction, maxNinstance))
allocate(constitutive_nonlocal_dLowerEdgePerSlipFamily(lattice_maxNslipFamily, maxNinstance))
allocate(constitutive_nonlocal_dLowerScrewPerSlipFamily(lattice_maxNslipFamily, maxNinstance))
constitutive_nonlocal_rhoSglEdgePos0 = -1.0_pReal
constitutive_nonlocal_rhoSglEdgeNeg0 = -1.0_pReal
constitutive_nonlocal_rhoSglScrewPos0 = -1.0_pReal
constitutive_nonlocal_rhoSglScrewNeg0 = -1.0_pReal
constitutive_nonlocal_rhoDipEdge0 = -1.0_pReal
constitutive_nonlocal_rhoDipScrew0 = -1.0_pReal
constitutive_nonlocal_burgersPerSlipFamily = 0.0_pReal
constitutive_nonlocal_lambda0PerSlipFamily = 0.0_pReal
constitutive_nonlocal_interactionSlipSlip = 0.0_pReal
constitutive_nonlocal_dLowerEdgePerSlipFamily = 0.0_pReal
constitutive_nonlocal_dLowerScrewPerSlipFamily = 0.0_pReal
!*** readout data from material.config file
rewind(file)
line = ''
section = 0
do while (IO_lc(IO_getTag(line,'<','>')) /= 'phase') ! wind forward to <phase>
read(file,'(a1024)',END=100) line
enddo
do ! read thru sections of phase part
read(file,'(a1024)',END=100) line
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1
output = 0 ! reset output counter
cycle
endif
if (section > 0 .and. phase_constitution(section) == constitutive_nonlocal_label) then ! one of my sections
i = phase_constitutionInstance(section) ! which instance of my constitution is present phase
positions = IO_stringPos(line,maxNchunks)
tag = IO_lc(IO_stringValue(line,positions,1)) ! extract key
select case(tag)
case('constitution','/nonlocal/')
cycle
case ('(output)')
output = output + 1
constitutive_nonlocal_output(output,i) = IO_lc(IO_stringValue(line,positions,2))
case ('lattice_structure')
constitutive_nonlocal_structureName(i) = IO_lc(IO_stringValue(line,positions,2))
case ('covera_ratio')
constitutive_nonlocal_CoverA(i) = IO_floatValue(line,positions,2)
case ('c11')
constitutive_nonlocal_C11(i) = IO_floatValue(line,positions,2)
case ('c12')
constitutive_nonlocal_C12(i) = IO_floatValue(line,positions,2)
case ('c13')
constitutive_nonlocal_C13(i) = IO_floatValue(line,positions,2)
case ('c33')
constitutive_nonlocal_C33(i) = IO_floatValue(line,positions,2)
case ('c44')
constitutive_nonlocal_C44(i) = IO_floatValue(line,positions,2)
case ('nslip')
forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_Nslip(f,i) = IO_intValue(line,positions,1+f)
case ('rhosgledgepos0')
forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_rhoSglEdgePos0(f,i) = IO_floatValue(line,positions,1+f)
case ('rhosgledgeneg0')
forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_rhoSglEdgeNeg0(f,i) = IO_floatValue(line,positions,1+f)
case ('rhosglscrewpos0')
forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_rhoSglScrewPos0(f,i) = IO_floatValue(line,positions,1+f)
case ('rhosglscrewneg0')
forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_rhoSglScrewNeg0(f,i) = IO_floatValue(line,positions,1+f)
case ('rhodipedge0')
forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_rhoDipEdge0(f,i) = IO_floatValue(line,positions,1+f)
case ('rhodipscrew0')
forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_rhoDipScrew0(f,i) = IO_floatValue(line,positions,1+f)
case ('lambda0')
forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_lambda0PerSlipFamily(f,i) = IO_floatValue(line,positions,1+f)
case ('burgers')
forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_burgersPerSlipFamily(f,i) = IO_floatValue(line,positions,1+f)
case('r')
constitutive_nonlocal_R(i) = IO_floatValue(line,positions,2)
case('ddipminedge')
forall (f = 1:lattice_maxNslipFamily) &
constitutive_nonlocal_dLowerEdgePerSlipFamily(f,i) = IO_floatValue(line,positions,1+f)
case('ddipminscrew')
forall (f = 1:lattice_maxNslipFamily) &
constitutive_nonlocal_dLowerScrewPerSlipFamily(f,i) = IO_floatValue(line,positions,1+f)
case('atomicvolume')
constitutive_nonlocal_atomicVolume(i) = IO_floatValue(line,positions,2)
case('dsd0')
constitutive_nonlocal_Dsd0(i) = IO_floatValue(line,positions,2)
case('qsd')
constitutive_nonlocal_Qsd(i) = IO_floatValue(line,positions,2)
case('atol_rho')
constitutive_nonlocal_aTolRho(i) = IO_floatValue(line,positions,2)
case ('interaction_slipslip')
forall (it = 1:lattice_maxNinteraction) constitutive_nonlocal_interactionSlipSlip(it,i) = IO_floatValue(line,positions,1+it)
case('d0')
constitutive_nonlocal_d0(i) = IO_floatValue(line,positions,2)
case('tauobs')
constitutive_nonlocal_tauObs(i) = IO_floatValue(line,positions,2)
case('vs')
constitutive_nonlocal_vs(i) = IO_floatValue(line,positions,2)
case('fattack')
constitutive_nonlocal_fattack(i) = IO_floatValue(line,positions,2)
case('rhosglscatter')
constitutive_nonlocal_rhoSglScatter(i) = IO_floatValue(line,positions,2)
case('surfacetransmissivity')
constitutive_nonlocal_surfaceTransmissivity(i) = IO_floatValue(line,positions,2)
case default
call IO_error(236,ext_msg=tag)
end select
endif
enddo
100 do i = 1,maxNinstance
constitutive_nonlocal_structure(i) = &
lattice_initializeStructure(constitutive_nonlocal_structureName(i), constitutive_nonlocal_CoverA(i)) ! our lattice structure is defined in the material.config file by the structureName (and the c/a ratio)
myStructure = constitutive_nonlocal_structure(i)
!*** sanity checks
if (myStructure < 1 .or. myStructure > 3) call IO_error(205)
if (sum(constitutive_nonlocal_Nslip(:,i)) <= 0_pInt) call IO_error(235,ext_msg='Nslip')
do o = 1,maxval(phase_Noutput)
if(len(constitutive_nonlocal_output(o,i)) > 64) call IO_error(666)
enddo
do f = 1,lattice_maxNslipFamily
if (constitutive_nonlocal_Nslip(f,i) > 0_pInt) then
if (constitutive_nonlocal_rhoSglEdgePos0(f,i) < 0.0_pReal) call IO_error(235,ext_msg='rhoSglEdgePos0')
if (constitutive_nonlocal_rhoSglEdgeNeg0(f,i) < 0.0_pReal) call IO_error(235,ext_msg='rhoSglEdgeNeg0')
if (constitutive_nonlocal_rhoSglScrewPos0(f,i) < 0.0_pReal) call IO_error(235,ext_msg='rhoSglScrewPos0')
if (constitutive_nonlocal_rhoSglScrewNeg0(f,i) < 0.0_pReal) call IO_error(235,ext_msg='rhoSglScrewNeg0')
if (constitutive_nonlocal_rhoDipEdge0(f,i) < 0.0_pReal) call IO_error(235,ext_msg='rhoDipEdge0')
if (constitutive_nonlocal_rhoDipScrew0(f,i) < 0.0_pReal) call IO_error(235,ext_msg='rhoDipScrew0')
if (constitutive_nonlocal_burgersPerSlipFamily(f,i) <= 0.0_pReal) call IO_error(235,ext_msg='burgers')
if (constitutive_nonlocal_lambda0PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(235,ext_msg='lambda0')
if (constitutive_nonlocal_dLowerEdgePerSlipFamily(f,i) <= 0.0_pReal) call IO_error(235,ext_msg='dDipMinEdge')
if (constitutive_nonlocal_dLowerScrewPerSlipFamily(f,i) <= 0.0_pReal) call IO_error(235,ext_msg='dDipMinScrew')
endif
enddo
if (any(constitutive_nonlocal_interactionSlipSlip(1:maxval(lattice_interactionSlipSlip(:,:,myStructure)),i) < 0.0_pReal)) &
call IO_error(235,ext_msg='interaction_SlipSlip')
if (constitutive_nonlocal_R(i) < 0.0_pReal) call IO_error(235,ext_msg='r')
if (constitutive_nonlocal_atomicVolume(i) <= 0.0_pReal) call IO_error(235,ext_msg='atomicVolume')
if (constitutive_nonlocal_Dsd0(i) <= 0.0_pReal) call IO_error(235,ext_msg='Dsd0')
if (constitutive_nonlocal_Qsd(i) <= 0.0_pReal) call IO_error(235,ext_msg='Qsd')
if (constitutive_nonlocal_aTolRho(i) <= 0.0_pReal) call IO_error(235,ext_msg='aTol_rho')
if (constitutive_nonlocal_d0(i) <= 0.0_pReal) call IO_error(235,ext_msg='d0')
if (constitutive_nonlocal_tauObs(i) <= 0.0_pReal) call IO_error(235,ext_msg='tauObs')
if (constitutive_nonlocal_vs(i) <= 0.0_pReal) call IO_error(235,ext_msg='vs')
if (constitutive_nonlocal_fattack(i) <= 0.0_pReal) call IO_error(235,ext_msg='fAttack')
if (constitutive_nonlocal_rhoSglScatter(i) < 0.0_pReal) call IO_error(235,ext_msg='rhoSglScatter')
if (constitutive_nonlocal_surfaceTransmissivity(i) < 0.0_pReal &
.or. constitutive_nonlocal_surfaceTransmissivity(i) > 1.0_pReal) call IO_error(235,ext_msg='surfaceTransmissivity')
!*** determine total number of active slip systems
constitutive_nonlocal_Nslip(1:lattice_maxNslipFamily,i) = min( lattice_NslipSystem(1:lattice_maxNslipFamily, myStructure), &
constitutive_nonlocal_Nslip(1:lattice_maxNslipFamily,i) ) ! we can't use more slip systems per family than specified in lattice
constitutive_nonlocal_totalNslip(i) = sum(constitutive_nonlocal_Nslip(1:lattice_maxNslipFamily,i))
enddo
!*** allocation of variables whose size depends on the total number of active slip systems
maxTotalNslip = maxval(constitutive_nonlocal_totalNslip)
allocate(constitutive_nonlocal_burgersPerSlipSystem(maxTotalNslip, maxNinstance))
constitutive_nonlocal_burgersPerSlipSystem = 0.0_pReal
allocate(constitutive_nonlocal_lambda0PerSlipSystem(maxTotalNslip, maxNinstance))
constitutive_nonlocal_lambda0PerSlipSystem = 0.0_pReal
allocate(constitutive_nonlocal_dLowerEdgePerSlipSystem(maxTotalNslip, maxNinstance))
constitutive_nonlocal_dLowerEdgePerSlipSystem = 0.0_pReal
allocate(constitutive_nonlocal_dLowerScrewPerSlipSystem(maxTotalNslip, maxNinstance))
constitutive_nonlocal_dLowerScrewPerSlipSystem = 0.0_pReal
allocate(constitutive_nonlocal_Qeff0(maxTotalNslip, maxNinstance))
constitutive_nonlocal_Qeff0 = 0.0_pReal
allocate(constitutive_nonlocal_forestProjectionEdge(maxTotalNslip, maxTotalNslip, maxNinstance))
constitutive_nonlocal_forestProjectionEdge = 0.0_pReal
allocate(constitutive_nonlocal_forestProjectionScrew(maxTotalNslip, maxTotalNslip, maxNinstance))
constitutive_nonlocal_forestProjectionScrew = 0.0_pReal
allocate(constitutive_nonlocal_interactionMatrixSlipSlip(maxTotalNslip, maxTotalNslip, maxNinstance))
constitutive_nonlocal_interactionMatrixSlipSlip = 0.0_pReal
allocate(constitutive_nonlocal_lattice2slip(1:3, 1:3, maxTotalNslip, maxNinstance))
constitutive_nonlocal_lattice2slip = 0.0_pReal
allocate(constitutive_nonlocal_rhoDotFlux(maxTotalNslip, 10, homogenization_maxNgrains, mesh_maxNips, mesh_NcpElems))
constitutive_nonlocal_rhoDotFlux = 0.0_pReal
allocate(constitutive_nonlocal_compatibility(2,maxTotalNslip, maxTotalNslip, FE_maxNipNeighbors, mesh_maxNips, mesh_NcpElems))
constitutive_nonlocal_compatibility = 0.0_pReal
do i = 1,maxNinstance
myStructure = constitutive_nonlocal_structure(i) ! lattice structure of this instance
!*** Inverse lookup of my slip system family and the slip system in lattice
l = 0_pInt
do f = 1,lattice_maxNslipFamily
do s = 1,constitutive_nonlocal_Nslip(f,i)
l = l + 1
constitutive_nonlocal_slipFamily(l,i) = f
constitutive_nonlocal_slipSystemLattice(l,i) = sum(lattice_NslipSystem(1:f-1, myStructure)) + s
enddo; enddo
!*** determine size of state array
ns = constitutive_nonlocal_totalNslip(i)
constitutive_nonlocal_sizeDotState(i) = size(constitutive_nonlocal_listBasicStates) * ns
constitutive_nonlocal_sizeDependentState(i) = (size(constitutive_nonlocal_listDependentStates) - 1_pInt) * ns + 6_pInt
constitutive_nonlocal_sizeState(i) = constitutive_nonlocal_sizeDotState(i) &
+ constitutive_nonlocal_sizeDependentState(i) &
+ size(constitutive_nonlocal_listOtherStates) * ns
!*** determine size of postResults array
do o = 1,maxval(phase_Noutput)
select case(constitutive_nonlocal_output(o,i))
case( 'rho', &
'delta', &
'rho_edge', &
'rho_screw', &
'rho_sgl', &
'delta_sgl', &
'rho_sgl_edge', &
'rho_sgl_edge_pos', &
'rho_sgl_edge_neg', &
'rho_sgl_screw', &
'rho_sgl_screw_pos', &
'rho_sgl_screw_neg', &
'rho_sgl_mobile', &
'rho_sgl_edge_mobile', &
'rho_sgl_edge_pos_mobile', &
'rho_sgl_edge_neg_mobile', &
'rho_sgl_screw_mobile', &
'rho_sgl_screw_pos_mobile', &
'rho_sgl_screw_neg_mobile', &
'rho_sgl_immobile', &
'rho_sgl_edge_immobile', &
'rho_sgl_edge_pos_immobile', &
'rho_sgl_edge_neg_immobile', &
'rho_sgl_screw_immobile', &
'rho_sgl_screw_pos_immobile', &
'rho_sgl_screw_neg_immobile', &
'rho_dip', &
'delta_dip', &
'rho_dip_edge', &
'rho_dip_screw', &
'excess_rho', &
'excess_rho_edge', &
'excess_rho_screw', &
'rho_forest', &
'shearrate', &
'resolvedstress', &
'resolvedstress_internal', &
'resolvedstress_external', &
'resistance', &
'rho_dot', &
'rho_dot_sgl', &
'rho_dot_dip', &
'rho_dot_gen', &
'rho_dot_gen_edge', &
'rho_dot_gen_screw', &
'rho_dot_sgl2dip', &
'rho_dot_ann_ath', &
'rho_dot_ann_the', &
'rho_dot_flux', &
'rho_dot_flux_edge', &
'rho_dot_flux_screw', &
'dislocationvelocity', &
'fluxdensity_edge_pos_x', &
'fluxdensity_edge_pos_y', &
'fluxdensity_edge_pos_z', &
'fluxdensity_edge_neg_x', &
'fluxdensity_edge_neg_y', &
'fluxdensity_edge_neg_z', &
'fluxdensity_screw_pos_x', &
'fluxdensity_screw_pos_y', &
'fluxdensity_screw_pos_z', &
'fluxdensity_screw_neg_x', &
'fluxdensity_screw_neg_y', &
'fluxdensity_screw_neg_z', &
'd_upper_edge', &
'd_upper_screw' )
mySize = constitutive_nonlocal_totalNslip(i)
case('internalstress')
mySize = 6_pInt
case default
mySize = 0_pInt
end select
if (mySize > 0_pInt) then ! any meaningful output found
constitutive_nonlocal_sizePostResult(o,i) = mySize
constitutive_nonlocal_sizePostResults(i) = constitutive_nonlocal_sizePostResults(i) + mySize
endif
enddo
!*** elasticity matrix and shear modulus according to material.config
select case (myStructure)
case(1:2) ! cubic(s)
forall(k=1:3)
forall(j=1:3) constitutive_nonlocal_Cslip_66(k,j,i) = constitutive_nonlocal_C12(i)
constitutive_nonlocal_Cslip_66(k,k,i) = constitutive_nonlocal_C11(i)
constitutive_nonlocal_Cslip_66(k+3,k+3,i) = constitutive_nonlocal_C44(i)
end forall
case(3:) ! all hex
constitutive_nonlocal_Cslip_66(1,1,i) = constitutive_nonlocal_C11(i)
constitutive_nonlocal_Cslip_66(2,2,i) = constitutive_nonlocal_C11(i)
constitutive_nonlocal_Cslip_66(3,3,i) = constitutive_nonlocal_C33(i)
constitutive_nonlocal_Cslip_66(1,2,i) = constitutive_nonlocal_C12(i)
constitutive_nonlocal_Cslip_66(2,1,i) = constitutive_nonlocal_C12(i)
constitutive_nonlocal_Cslip_66(1,3,i) = constitutive_nonlocal_C13(i)
constitutive_nonlocal_Cslip_66(3,1,i) = constitutive_nonlocal_C13(i)
constitutive_nonlocal_Cslip_66(2,3,i) = constitutive_nonlocal_C13(i)
constitutive_nonlocal_Cslip_66(3,2,i) = constitutive_nonlocal_C13(i)
constitutive_nonlocal_Cslip_66(4,4,i) = constitutive_nonlocal_C44(i)
constitutive_nonlocal_Cslip_66(5,5,i) = constitutive_nonlocal_C44(i)
constitutive_nonlocal_Cslip_66(6,6,i) = 0.5_pReal*(constitutive_nonlocal_C11(i)- constitutive_nonlocal_C12(i))
end select
constitutive_nonlocal_Cslip_66(1:6,1:6,i) = math_Mandel3333to66(math_Voigt66to3333(constitutive_nonlocal_Cslip_66(1:6,1:6,i)))
constitutive_nonlocal_Cslip_3333(1:3,1:3,1:3,1:3,i) = math_Voigt66to3333(constitutive_nonlocal_Cslip_66(1:6,1:6,i))
constitutive_nonlocal_Gmod(i) = 0.2_pReal * ( constitutive_nonlocal_C11(i) - constitutive_nonlocal_C12(i) &
+ 3.0_pReal*constitutive_nonlocal_C44(i) ) ! (C11iso-C12iso)/2 with C11iso=(3*C11+2*C12+4*C44)/5 and C12iso=(C11+4*C12-2*C44)/5
constitutive_nonlocal_nu(i) = ( constitutive_nonlocal_C11(i) + 4.0_pReal*constitutive_nonlocal_C12(i) &
- 2.0_pReal*constitutive_nonlocal_C44(i) ) &
/ ( 4.0_pReal*constitutive_nonlocal_C11(i) + 6.0_pReal*constitutive_nonlocal_C12(i) &
+ 2.0_pReal*constitutive_nonlocal_C44(i) ) ! C12iso/(C11iso+C12iso) with C11iso=(3*C11+2*C12+4*C44)/5 and C12iso=(C11+4*C12-2*C44)/5
do s1 = 1,ns
f = constitutive_nonlocal_slipFamily(s1,i)
!*** burgers vector, mean free path prefactor and minimum dipole distance for each slip system
constitutive_nonlocal_burgersPerSlipSystem(s1,i) = constitutive_nonlocal_burgersPerSlipFamily(f,i)
constitutive_nonlocal_lambda0PerSlipSystem(s1,i) = constitutive_nonlocal_lambda0PerSlipFamily(f,i)
constitutive_nonlocal_dLowerEdgePerSlipSystem(s1,i) = constitutive_nonlocal_dLowerEdgePerSlipFamily(f,i)
constitutive_nonlocal_dLowerScrewPerSlipSystem(s1,i) = constitutive_nonlocal_dLowerScrewPerSlipFamily(f,i)
do s2 = 1,ns
!*** calculation of forest projections for edge and screw dislocations. s2 acts as forest for s1
constitutive_nonlocal_forestProjectionEdge(s1,s2,i) &
= abs(math_mul3x3(lattice_sn(1:3,constitutive_nonlocal_slipSystemLattice(s1,i),myStructure), &
lattice_st(1:3,constitutive_nonlocal_slipSystemLattice(s2,i),myStructure))) ! forest projection of edge dislocations is the projection of (t = b x n) onto the slip normal of the respective slip plane
constitutive_nonlocal_forestProjectionScrew(s1,s2,i) &
= abs(math_mul3x3(lattice_sn(1:3,constitutive_nonlocal_slipSystemLattice(s1,i),myStructure), &
lattice_sd(1:3,constitutive_nonlocal_slipSystemLattice(s2,i),myStructure))) ! forest projection of screw dislocations is the projection of b onto the slip normal of the respective splip plane
!*** calculation of interaction matrices
constitutive_nonlocal_interactionMatrixSlipSlip(s1,s2,i) &
= constitutive_nonlocal_interactionSlipSlip(lattice_interactionSlipSlip(constitutive_nonlocal_slipSystemLattice(s1,i), &
constitutive_nonlocal_slipSystemLattice(s2,i), &
myStructure), &
i)
enddo
!*** rotation matrix from lattice configuration to slip system
constitutive_nonlocal_lattice2slip(1:3,1:3,s1,i) &
= math_transpose3x3( reshape((/ lattice_sd(1:3, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure), &
-lattice_st(1:3, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure), &
lattice_sn(1:3, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure)/), (/3,3/)))
enddo
!*** calculation of prefactor for activation enthalpy for dislocation glide
constitutive_nonlocal_Qeff0(1:ns,i) = constitutive_nonlocal_burgersPerSlipSystem(1:ns,i) ** 3.0_pReal &
* sqrt(0.5_pReal * constitutive_nonlocal_d0(i) ** 3.0_pReal &
* constitutive_nonlocal_Gmod(i) * constitutive_nonlocal_tauObs(i))
enddo
endsubroutine
!*********************************************************************
!* initial microstructural state (just the "basic" states) *
!*********************************************************************
function constitutive_nonlocal_stateInit(myInstance)
use prec, only: pReal, &
pInt
use lattice, only: lattice_maxNslipFamily
use math, only: math_sampleGaussVar
implicit none
!*** input variables
integer(pInt), intent(in) :: myInstance ! number specifying the current instance of the constitution
!*** output variables
real(pReal), dimension(constitutive_nonlocal_sizeState(myInstance)) :: &
constitutive_nonlocal_stateInit
!*** local variables
real(pReal), dimension(constitutive_nonlocal_totalNslip(myInstance)) :: &
rhoSglEdgePos, & ! positive edge dislocation density
rhoSglEdgeNeg, & ! negative edge dislocation density
rhoSglScrewPos, & ! positive screw dislocation density
rhoSglScrewNeg, & ! negative screw dislocation density
rhoSglEdgePosUsed, & ! used positive edge dislocation density
rhoSglEdgeNegUsed, & ! used negative edge dislocation density
rhoSglScrewPosUsed, & ! used positive screw dislocation density
rhoSglScrewNegUsed, & ! used negative screw dislocation density
rhoDipEdge, & ! edge dipole dislocation density
rhoDipScrew ! screw dipole dislocation density
integer(pInt) ns, & ! short notation for total number of active slip systems
f, & ! index of lattice family
from, &
upto, &
s, & ! index of slip system
i
real(pReal), dimension(2) :: noise
constitutive_nonlocal_stateInit = 0.0_pReal
ns = constitutive_nonlocal_totalNslip(myInstance)
!*** set the basic state variables
do f = 1,lattice_maxNslipFamily
from = 1 + sum(constitutive_nonlocal_Nslip(1:f-1,myInstance))
upto = sum(constitutive_nonlocal_Nslip(1:f,myInstance))
do s = from,upto
do i = 1,2
noise(i) = math_sampleGaussVar(0.0_pReal, constitutive_nonlocal_rhoSglScatter(myInstance))
enddo
rhoSglEdgePos(s) = constitutive_nonlocal_rhoSglEdgePos0(f, myInstance) + noise(1)
rhoSglEdgeNeg(s) = constitutive_nonlocal_rhoSglEdgeNeg0(f, myInstance) + noise(1)
rhoSglScrewPos(s) = constitutive_nonlocal_rhoSglScrewPos0(f, myInstance) + noise(2)
rhoSglScrewNeg(s) = constitutive_nonlocal_rhoSglScrewNeg0(f, myInstance) + noise(2)
enddo
rhoSglEdgePosUsed(from:upto) = 0.0_pReal
rhoSglEdgeNegUsed(from:upto) = 0.0_pReal
rhoSglScrewPosUsed(from:upto) = 0.0_pReal
rhoSglScrewNegUsed(from:upto) = 0.0_pReal
rhoDipEdge(from:upto) = constitutive_nonlocal_rhoDipEdge0(f, myInstance)
rhoDipScrew(from:upto) = constitutive_nonlocal_rhoDipScrew0(f, myInstance)
enddo
!*** put everything together and in right order
constitutive_nonlocal_stateInit( 1: ns) = rhoSglEdgePos
constitutive_nonlocal_stateInit( ns+1: 2*ns) = rhoSglEdgeNeg
constitutive_nonlocal_stateInit( 2*ns+1: 3*ns) = rhoSglScrewPos
constitutive_nonlocal_stateInit( 3*ns+1: 4*ns) = rhoSglScrewNeg
constitutive_nonlocal_stateInit( 4*ns+1: 5*ns) = rhoSglEdgePosUsed
constitutive_nonlocal_stateInit( 5*ns+1: 6*ns) = rhoSglEdgeNegUsed
constitutive_nonlocal_stateInit( 6*ns+1: 7*ns) = rhoSglScrewPosUsed
constitutive_nonlocal_stateInit( 7*ns+1: 8*ns) = rhoSglScrewNegUsed
constitutive_nonlocal_stateInit( 8*ns+1: 9*ns) = rhoDipEdge
constitutive_nonlocal_stateInit( 9*ns+1:10*ns) = rhoDipScrew
endfunction
!*********************************************************************
!* absolute state tolerance *
!*********************************************************************
pure function constitutive_nonlocal_aTolState(myInstance)
use prec, only: pReal, &
pInt
implicit none
!*** input variables
integer(pInt), intent(in) :: myInstance ! number specifying the current instance of the constitution
!*** output variables
real(pReal), dimension(constitutive_nonlocal_sizeState(myInstance)) :: &
constitutive_nonlocal_aTolState ! absolute state tolerance for the current instance of this constitution
!*** local variables
constitutive_nonlocal_aTolState = constitutive_nonlocal_aTolRho(myInstance)
endfunction
!*********************************************************************
!* calculates homogenized elacticity matrix *
!*********************************************************************
pure function constitutive_nonlocal_homogenizedC(state,g,ip,el)
use prec, only: pReal, &
pInt, &
p_vec
use mesh, only: mesh_NcpElems, &
mesh_maxNips
use material, only: homogenization_maxNgrains, &
material_phase, &
phase_constitutionInstance
implicit none
!*** input variables
integer(pInt), intent(in) :: g, & ! current grain ID
ip, & ! current integration point
el ! current element
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state ! microstructural state
!*** output variables
real(pReal), dimension(6,6) :: constitutive_nonlocal_homogenizedC ! homogenized elasticity matrix
!*** local variables
integer(pInt) myInstance ! current instance of this constitution
myInstance = phase_constitutionInstance(material_phase(g,ip,el))
constitutive_nonlocal_homogenizedC = constitutive_nonlocal_Cslip_66(1:6,1:6,myInstance)
endfunction
!*********************************************************************
!* calculates quantities characterizing the microstructure *
!*********************************************************************
subroutine constitutive_nonlocal_microstructure(state, Temperature, Fe, g, ip, el)
use prec, only: pReal, &
pInt, &
p_vec
use IO, only: IO_error
use math, only: math_Mandel33to6, &
math_mul33x33, &
math_mul33x3, &
math_invert3x3, &
math_transpose3x3, &
pi
use debug, only: debug_verbosity, &
debug_selectiveDebugger, &
debug_g, &
debug_i, &
debug_e
use mesh, only: mesh_NcpElems, &
mesh_maxNips, &
mesh_element, &
mesh_node0, &
FE_Nips, &
mesh_ipCenterOfGravity, &
mesh_ipVolume, &
mesh_periodicSurface
use material, only: homogenization_maxNgrains, &
material_phase, &
phase_localConstitution, &
phase_constitutionInstance
use lattice, only: lattice_sd, &
lattice_sn, &
lattice_st
implicit none
!*** input variables
integer(pInt), intent(in) :: g, & ! current grain ID
ip, & ! current integration point
el ! current element
real(pReal), intent(in) :: Temperature ! temperature
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
Fe ! elastic deformation gradient
!*** input/output variables
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(inout) :: &
state ! microstructural state
!*** output variables
!*** local variables
integer(pInt) neighboring_el, & ! element number of neighboring material point
neighboring_ip, & ! integration point of neighboring material point
instance, & ! my instance of this constitution
neighboring_instance, & ! instance of this constitution of neighboring material point
latticeStruct, & ! my lattice structure
neighboring_latticeStruct, & ! lattice structure of neighboring material point
phase, &
neighboring_phase, &
ns, & ! total number of active slip systems at my material point
neighboring_ns, & ! total number of active slip systems at neighboring material point
c, & ! index of dilsocation character (edge, screw)
s, & ! slip system index
t, & ! index of dilsocation type (e+, e-, s+, s-, used e+, used e-, used s+, used s-)
dir, &
deltaX, deltaY, deltaZ, &
side, &
j
integer(pInt), dimension(2,3) :: periodicImages
real(pReal) nu, & ! poisson's ratio
x, y, z, & ! coordinates of connection vector in neighboring lattice frame
xsquare, ysquare, zsquare, & ! squares of respective coordinates
distance, & ! length of connection vector
segmentLength, & ! segment length of dislocations
lambda, &
R, Rsquare, Rcube, &
denominator, &
flipSign, &
neighboring_ipVolumeSideLength, &
detFe
real(pReal), dimension(3) :: connection, & ! connection vector between me and my neighbor in the deformed configuration
connection_neighboringLattice, & ! connection vector between me and my neighbor in the lattice configuration of my neighbor
connection_neighboringSlip, & ! connection vector between me and my neighbor in the slip system frame of my neighbor
maxCoord, minCoord, &
meshSize, &
ipCoords, &
neighboring_ipCoords
real(pReal), dimension(3,3) :: sigma, & ! dislocation stress for one slip system in neighboring material point's slip system frame
Tdislo_neighboringLattice, & ! dislocation stress as 2nd Piola-Kirchhoff stress at neighboring material point
Tdislo, & ! dislocation stress as 2nd Piola-Kirchhoff stress at my material point
invFe, & ! inverse of my elastic deformation gradient
neighboring_invFe, &
neighboringLattice2myLattice ! mapping from neighboring MPs lattice configuration to my lattice configuration
real(pReal), dimension(2,2,maxval(constitutive_nonlocal_totalNslip)) :: &
neighboring_rhoExcess ! excess density at neighboring material point (edge/screw,mobile/dead,slipsystem)
real(pReal), dimension(2,maxval(constitutive_nonlocal_totalNslip)) :: &
rhoExcessDead
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),8) :: &
rhoSgl ! single dislocation density (edge+, edge-, screw+, screw-, used edge+, used edge-, used screw+, used screw-)
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),2) :: &
rhoDip ! dipole dislocation density (edge, screw)
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: &
rhoForest, & ! forest dislocation density
tauThreshold ! threshold shear stress
logical inversionError
phase = material_phase(g,ip,el)
instance = phase_constitutionInstance(phase)
latticeStruct = constitutive_nonlocal_structure(instance)
ns = constitutive_nonlocal_totalNslip(instance)
!*** get basic states
forall (s = 1:ns, t = 1:4) &
rhoSgl(s,t) = max(state(g,ip,el)%p((t-1)*ns+s), 0.0_pReal) ! ensure positive single mobile densities
forall (t = 5:8) &
rhoSgl(1:ns,t) = state(g,ip,el)%p((t-1)*ns+1:t*ns)
forall (s = 1:ns, c = 1:2) &
rhoDip(s,c) = max(state(g,ip,el)%p((7+c)*ns+s), 0.0_pReal) ! ensure positive dipole densities
!*** calculate the forest dislocation density
!*** (= projection of screw and edge dislocations)
forall (s = 1:ns) &
rhoForest(s) = dot_product((sum(abs(rhoSgl(1:ns,(/1,2,5,6/))),2) + rhoDip(1:ns,1)), &
constitutive_nonlocal_forestProjectionEdge(s,1:ns,instance)) &
+ dot_product((sum(abs(rhoSgl(1:ns,(/3,4,7,8/))),2) + rhoDip(1:ns,2)), &
constitutive_nonlocal_forestProjectionScrew(s,1:ns,instance))
!*** calculate the threshold shear stress for dislocation slip
forall (s = 1:ns) &
tauThreshold(s) = constitutive_nonlocal_Gmod(instance) &
* constitutive_nonlocal_burgersPerSlipSystem(s,instance) &
* sqrt(dot_product((sum(abs(rhoSgl),2) + sum(abs(rhoDip),2)), &
constitutive_nonlocal_interactionMatrixSlipSlip(s,1:ns,instance)))
!*** calculate the dislocation stress of the neighboring excess dislocation densities
!*** zero for material points of local constitution
Tdislo = 0.0_pReal
if (.not. phase_localConstitution(phase)) then
call math_invert3x3(Fe(1:3,1:3,1,ip,el), invFe, detFe, inversionError)
! if (inversionError) then
! return
! endif
!* in case of periodic surfaces we have to find out how many periodic images in each direction we need
do dir = 1,3
maxCoord(dir) = maxval(mesh_node0(dir,:))
minCoord(dir) = minval(mesh_node0(dir,:))
enddo
meshSize = maxCoord - minCoord
ipCoords = mesh_ipCenterOfGravity(1:3,ip,el)
periodicImages = 0_pInt
do dir = 1,3
if (mesh_periodicSurface(dir)) then
periodicImages(1,dir) = floor((ipCoords(dir) - constitutive_nonlocal_R(instance) - minCoord(dir)) / meshSize(dir), pInt)
periodicImages(2,dir) = ceiling((ipCoords(dir) + constitutive_nonlocal_R(instance) - maxCoord(dir)) / meshSize(dir), pInt)
endif
enddo
!* loop through all material points (also through their periodic images if present),
!* but only consider nonlocal neighbors within a certain cutoff radius R
do neighboring_el = 1,mesh_NcpElems
do neighboring_ip = 1,FE_Nips(mesh_element(2,neighboring_el))
neighboring_phase = material_phase(g,neighboring_ip,neighboring_el)
if (phase_localConstitution(neighboring_phase)) then
cycle
endif
neighboring_instance = phase_constitutionInstance(neighboring_phase)
neighboring_latticeStruct = constitutive_nonlocal_structure(neighboring_instance)
neighboring_ns = constitutive_nonlocal_totalNslip(neighboring_instance)
call math_invert3x3(Fe(1:3,1:3,1,neighboring_ip,neighboring_el), neighboring_invFe, detFe, inversionError)
! if (inversionError) then
! return
! endif
neighboring_ipVolumeSideLength = mesh_ipVolume(neighboring_ip,neighboring_el) ** (1.0_pReal/3.0_pReal) ! reference volume used here
forall (s = 1:neighboring_ns, c = 1:2) &
neighboring_rhoExcess(c,1,s) = state(g,neighboring_ip,neighboring_el)%p((2*c-2)*neighboring_ns+s) & ! positive mobiles
- state(g,neighboring_ip,neighboring_el)%p((2*c-1)*neighboring_ns+s) ! negative mobiles
forall (s = 1:neighboring_ns, c = 1:2) &
neighboring_rhoExcess(c,2,s) = abs(state(g,neighboring_ip,neighboring_el)%p((2*c+2)*neighboring_ns+s)) & ! positive deads
- abs(state(g,neighboring_ip,neighboring_el)%p((2*c+3)*neighboring_ns+s)) ! negative deads
nu = constitutive_nonlocal_nu(neighboring_instance)
Tdislo_neighboringLattice = 0.0_pReal
do deltaX = periodicImages(1,1),periodicImages(2,1)
do deltaY = periodicImages(1,2),periodicImages(2,2)
do deltaZ = periodicImages(1,3),periodicImages(2,3)
!* regular case
if (neighboring_el /= el .or. neighboring_ip /= ip &
.or. deltaX /= 0_pInt .or. deltaY /= 0_pInt .or. deltaZ /= 0_pInt) then
neighboring_ipCoords = mesh_ipCenterOfGravity(1:3,neighboring_ip,neighboring_el) &
+ (/real(deltaX,pReal), real(deltaY,pReal), real(deltaZ,pReal)/) * meshSize
connection = neighboring_ipCoords - ipCoords
distance = sqrt(sum(connection * connection))
if (distance > constitutive_nonlocal_R(instance)) then
cycle
endif
!* the segment length is the minimum of the third root of the control volume and the ip distance
!* this ensures, that the central MP never sits on a neighboring dislocation segment
connection_neighboringLattice = math_mul33x3(neighboring_invFe, connection)
segmentLength = min(neighboring_ipVolumeSideLength, distance)
!* loop through all slip systems of the neighboring material point
!* and add up the stress contributions from egde and screw excess on these slip systems (if significant)
do s = 1,neighboring_ns
if (all(abs(neighboring_rhoExcess(:,:,s)) < constitutive_nonlocal_aTolRho(instance))) then
cycle ! not significant
endif
!* map the connection vector from the lattice into the slip system frame
connection_neighboringSlip = math_mul33x3(constitutive_nonlocal_lattice2slip(1:3,1:3,s,neighboring_instance), &
connection_neighboringLattice)
!* edge contribution to stress
sigma = 0.0_pReal
x = connection_neighboringSlip(1)
y = connection_neighboringSlip(2)
z = connection_neighboringSlip(3)
xsquare = x * x
ysquare = y * y
zsquare = z * z
do j = 1,2
if (abs(neighboring_rhoExcess(1,j,s)) < constitutive_nonlocal_aTolRho(instance)) then
cycle
elseif (j > 1) then
x = connection_neighboringSlip(1) + sign(0.5_pReal * segmentLength, &
state(g,neighboring_ip,neighboring_el)%p(4*neighboring_ns+s) &
- state(g,neighboring_ip,neighboring_el)%p(5*neighboring_ns+s))
xsquare = x * x
endif
flipSign = sign(1.0_pReal, -y)
do side = 1,-1,-2
lambda = real(side,pReal) * 0.5_pReal * segmentLength - y
R = sqrt(xsquare + zsquare + lambda * lambda)
Rsquare = R * R
Rcube = Rsquare * R
denominator = R * (R + flipSign * lambda)
if (denominator == 0.0_pReal) then
call IO_error(237,el,ip,g)
endif
sigma(1,1) = sigma(1,1) - real(side,pReal) * flipSign * z / denominator &
* (1.0_pReal + xsquare / Rsquare + xsquare / denominator) &
* neighboring_rhoExcess(1,j,s)
sigma(2,2) = sigma(2,2) - real(side,pReal) * (flipSign * 2.0_pReal * nu * z / denominator + z * lambda / Rcube)&
* neighboring_rhoExcess(1,j,s)
sigma(3,3) = sigma(3,3) + real(side,pReal) * flipSign * z / denominator &
* (1.0_pReal - zsquare / Rsquare - zsquare / denominator) &
* neighboring_rhoExcess(1,j,s)
sigma(1,2) = sigma(1,2) + real(side,pReal) * x * z / Rcube * neighboring_rhoExcess(1,j,s)
sigma(1,3) = sigma(1,3) + real(side,pReal) * flipSign * x / denominator &
* (1.0_pReal - zsquare / Rsquare - zsquare / denominator) &
* neighboring_rhoExcess(1,j,s)
sigma(2,3) = sigma(2,3) - real(side,pReal) * (nu / R - zsquare / Rcube) * neighboring_rhoExcess(1,j,s)
enddo
enddo
!* screw contribution to stress
x = connection_neighboringSlip(1) ! have to restore this value, because position might have been adapted for edge deads before
do j = 1,2
if (abs(neighboring_rhoExcess(2,j,s)) < constitutive_nonlocal_aTolRho(instance)) then
cycle
elseif (j > 1) then
y = connection_neighboringSlip(2) + sign(0.5_pReal * segmentLength, &
state(g,neighboring_ip,neighboring_el)%p(6*neighboring_ns+s) &
- state(g,neighboring_ip,neighboring_el)%p(7*neighboring_ns+s))
ysquare = y * y
endif
flipSign = sign(1.0_pReal, x)
do side = 1,-1,-2
lambda = x + real(side,pReal) * 0.5_pReal * segmentLength
R = sqrt(ysquare + zsquare + lambda * lambda)
Rsquare = R * R
Rcube = Rsquare * R
denominator = R * (R + flipSign * lambda)
if (denominator == 0.0_pReal) then
call IO_error(237,el,ip,g)
endif
sigma(1,2) = sigma(1,2) - real(side,pReal) * flipSign * z * (1.0_pReal - nu) / denominator &
* neighboring_rhoExcess(2,j,s)
sigma(1,3) = sigma(1,3) + real(side,pReal) * flipSign * y * (1.0_pReal - nu) / denominator &
* neighboring_rhoExcess(2,j,s)
enddo
enddo
if (all(abs(sigma) < 1.0e-10_pReal)) then ! SIGMA IS NOT A REAL STRESS, THATS WHY WE NEED A REALLY SMALL VALUE HERE
cycle
endif
!* copy symmetric parts
sigma(2,1) = sigma(1,2)
sigma(3,1) = sigma(1,3)
sigma(3,2) = sigma(2,3)
!* scale stresses and map them into the neighboring material point's lattice configuration
sigma = sigma * constitutive_nonlocal_Gmod(neighboring_instance) &
* constitutive_nonlocal_burgersPerSlipSystem(s,neighboring_instance) &
/ (4.0_pReal * pi * (1.0_pReal - nu)) &
* mesh_ipVolume(neighboring_ip,neighboring_el) / segmentLength ! reference volume is used here (according to the segment length calculation)
Tdislo_neighboringLattice = Tdislo_neighboringLattice &
+ math_mul33x33(math_transpose3x3(constitutive_nonlocal_lattice2slip(1:3,1:3,s,neighboring_instance)), &
math_mul33x33(sigma, constitutive_nonlocal_lattice2slip(1:3,1:3,s,neighboring_instance)))
enddo ! slip system loop
!* special case of central ip volume
!* only consider dead dislocations
!* we assume that they all sit at a distance equal to half the third root of V
!* in direction of the according slip direction
else
forall (s = 1:ns, c = 1:2) &
rhoExcessDead(c,s) = state(g,ip,el)%p((2*c+2)*ns+s) & ! positive deads (here we use symmetry: if this has negative sign it is treated as negative density at positive position instead of positive density at negative position)
+ state(g,ip,el)%p((2*c+3)*ns+s) ! negative deads (here we use symmetry: if this has negative sign it is treated as positive density at positive position instead of negative density at negative position)
do s = 1,ns
if (all(abs(rhoExcessDead(:,s)) < constitutive_nonlocal_aTolRho(instance))) then
cycle ! not significant
endif
sigma = 0.0_pReal ! all components except for sigma13 are zero
sigma(1,3) = - (rhoExcessDead(1,s) + rhoExcessDead(2,s) * (1.0_pReal - nu)) * neighboring_ipVolumeSideLength &
* constitutive_nonlocal_Gmod(instance) * constitutive_nonlocal_burgersPerSlipSystem(s,instance) &
/ (sqrt(2.0_pReal) * pi * (1.0_pReal - nu))
sigma(3,1) = sigma(1,3)
Tdislo_neighboringLattice = Tdislo_neighboringLattice &
+ math_mul33x33(math_transpose3x3(constitutive_nonlocal_lattice2slip(1:3,1:3,s,instance)), &
math_mul33x33(sigma, constitutive_nonlocal_lattice2slip(1:3,1:3,s,instance)))
enddo ! slip system loop
endif
enddo ! deltaZ loop
enddo ! deltaY loop
enddo ! deltaX loop
!* map the stress from the neighboring MP's lattice configuration into the deformed configuration
!* and back into my lattice configuration
neighboringLattice2myLattice = math_mul33x33(invFe, Fe(1:3,1:3,1,neighboring_ip,neighboring_el))
Tdislo = Tdislo + math_mul33x33(neighboringLattice2myLattice, &
math_mul33x33(Tdislo_neighboringLattice, math_transpose3x3(neighboringLattice2myLattice)))
enddo ! ip loop
enddo ! element loop
endif
!*** set dependent states
state(g,ip,el)%p(10*ns+1:11*ns) = rhoForest
state(g,ip,el)%p(11*ns+1:12*ns) = tauThreshold
state(g,ip,el)%p(12*ns+1:12*ns+6) = math_Mandel33to6(Tdislo)
#ifndef _OPENMP
if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then
write(6,*)
write(6,'(a,i8,x,i2,x,i1)') '<< CONST >> nonlocal_microstructure at el ip g',el,ip,g
write(6,*)
write(6,'(a,/,12(x),12(e10.3,x))') '<< CONST >> rhoForest', rhoForest
write(6,'(a,/,12(x),12(f10.5,x))') '<< CONST >> tauThreshold / MPa', tauThreshold/1e6
write(6,'(a,/,3(12(x),3(f10.5,x),/))') '<< CONST >> Tdislocation / MPa', Tdislo/1e6
endif
#endif
endsubroutine
!*********************************************************************
!* calculates kinetics *
!*********************************************************************
subroutine constitutive_nonlocal_kinetics(v, Tstar_v, Temperature, state, g, ip, el, dv_dtau)
use prec, only: pReal, &
pInt, &
p_vec
use math, only: math_mul6x6, &
math_Mandel6to33
use debug, only: debug_verbosity, &
debug_selectiveDebugger, &
debug_g, &
debug_i, &
debug_e
use mesh, only: mesh_NcpElems, &
mesh_maxNips
use material, only: homogenization_maxNgrains, &
material_phase, &
phase_constitutionInstance
use lattice, only: lattice_Sslip, &
lattice_Sslip_v
implicit none
!*** input variables
integer(pInt), intent(in) :: g, & ! current grain number
ip, & ! current integration point
el ! current element number
real(pReal), intent(in) :: Temperature ! temperature
real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola-Kirchhoff stress in Mandel notation
!*** input/output variables
type(p_vec), intent(inout) :: state ! microstructural state
!*** output variables
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4), &
intent(out) :: v ! velocity
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4), &
intent(out), optional :: dv_dtau ! velocity derivative with respect to resolved shear stress
!*** local variables
integer(pInt) myInstance, & ! current instance of this constitution
myStructure, & ! current lattice structure
ns, & ! short notation for the total number of active slip systems
s, & ! index of my current slip system
t ! index of dislocation character
real(pReal), dimension(6) :: Tdislocation_v ! dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: &
tauThreshold, & ! threshold shear stress
tau, & ! resolved shear stress
rhoForest ! forest dislocation density
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: &
rhoSgl
real(pReal) boltzmannProbability, &
tauRel, & ! relative thermally active resolved shear stress
wallFunc, & ! functions reflecting the shape of the obstacle wall (see PhD thesis Mohles p.53)
timeRatio ! ratio of travel to dwell time
myInstance = phase_constitutionInstance(material_phase(g,ip,el))
myStructure = constitutive_nonlocal_structure(myInstance)
ns = constitutive_nonlocal_totalNslip(myInstance)
forall (s = 1:ns, t = 1:4) &
rhoSgl(s,t) = max(state%p((t-1)*ns+s), 0.0_pReal)
rhoForest = state%p(10*ns+1:11*ns)
tauThreshold = state%p(11*ns+1:12*ns)
Tdislocation_v = state%p(12*ns+1:12*ns+6)
tau = 0.0_pReal
v = 0.0_pReal
if (present(dv_dtau)) dv_dtau = 0.0_pReal
if (Temperature > 0.0_pReal) then
do s = 1,ns
tau(s) = math_mul6x6(Tstar_v + Tdislocation_v, &
lattice_Sslip_v(:,constitutive_nonlocal_slipSystemLattice(s,myInstance),myStructure))
!*** Only if the resolved shear stress exceeds the threshold stress, dislocations are able to cut the dislocation forest.
!*** In contrast to small atomic obstacles the forest can't be overcome by thermal activation.
!***
!*** mean travel distance
!*** The mean dislocation velocity is calculated as: --------------------------
!*** dwell time + travel time
!***
!*** with : mean travel distance = inverse of the root of forest density
!*** dwell time = inverse of attack frequency times probability of success
!*** travel time = mean travel distance over velocity of sound
tauRel = (abs(tau(s)) - tauThreshold(s)) / constitutive_nonlocal_tauObs(myInstance)
if (tauRel > 0.0_pReal .and. tauRel < 1.0_pReal) then
wallFunc = 4.0_pReal * sqrt(2.0_pReal) / 3.0_pReal * sqrt(1.0_pReal - tauRel) / tauRel
boltzmannProbability = exp(- constitutive_nonlocal_Qeff0(s,myInstance) * wallFunc / (kB * Temperature))
timeRatio = boltzmannProbability * constitutive_nonlocal_fattack(myInstance) &
/ (constitutive_nonlocal_vs(myInstance) * sqrt(rhoForest(s)))
v(s,1:4) = sign(constitutive_nonlocal_vs(myInstance),tau(s)) * timeRatio / (1.0_pReal + timeRatio)
if (present(dv_dtau)) then
dv_dtau(s,1:4) = abs(v(s,1:4)) * constitutive_nonlocal_Qeff0(s,myInstance) / (kB * Temperature * (1.0_pReal + timeRatio)) &
* 0.5_pReal * wallFunc * (2.0_pReal - tauRel) / ((1.0_pReal - tauRel) * (abs(tau(s)) - tauThreshold(s)))
endif
!*** If resolved stress exceeds threshold plus obstacle stress, the probability for thermal activation is 1.
!*** The tangent is zero, since no dependency of tau.
elseif (tauRel >= 1.0_pReal) then
v(s,1:4) = sign(constitutive_nonlocal_vs(myInstance), tau(s)) * constitutive_nonlocal_fattack(myInstance) &
/ (constitutive_nonlocal_vs(myInstance) * sqrt(rhoForest(s)) + constitutive_nonlocal_fattack(myInstance))
endif
enddo
endif
state%p(12*ns+6+1:13*ns+6) = v(1:ns,1)
state%p(13*ns+6+1:14*ns+6) = v(1:ns,3)
#ifndef _OPENMP
if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then
write(6,*)
write(6,'(a,i8,x,i2,x,i1)') '<< CONST >> nonlocal_kinetics at el ip g',el,ip,g
write(6,*)
write(6,'(a,/,12(x),12(f12.5,x))') '<< CONST >> tau / MPa', tau / 1e6_pReal
write(6,'(a,/,4(12(x),12(f12.5,x),/))') '<< CONST >> v / 1e-3m/s', v * 1e3
endif
#endif
endsubroutine
!*********************************************************************
!* calculates plastic velocity gradient and its tangent *
!*********************************************************************
subroutine constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar99, Tstar_v, Temperature, state, g, ip, el)
use prec, only: pReal, &
pInt, &
p_vec
use math, only: math_Plain3333to99, &
math_mul6x6, &
math_Mandel6to33
use debug, only: debug_verbosity, &
debug_selectiveDebugger, &
debug_g, &
debug_i, &
debug_e
use mesh, only: mesh_NcpElems, &
mesh_maxNips
use material, only: homogenization_maxNgrains, &
material_phase, &
phase_constitutionInstance
use lattice, only: lattice_Sslip, &
lattice_Sslip_v
implicit none
!*** input variables
integer(pInt), intent(in) :: g, & ! current grain number
ip, & ! current integration point
el ! current element number
real(pReal), intent(in) :: Temperature ! temperature
real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola-Kirchhoff stress in Mandel notation
!*** input/output variables
type(p_vec), intent(inout) :: state ! microstructural state
!*** output variables
real(pReal), dimension(3,3), intent(out) :: Lp ! plastic velocity gradient
real(pReal), dimension(9,9), intent(out) :: dLp_dTstar99 ! derivative of Lp with respect to Tstar (9x9 matrix)
!*** local variables
integer(pInt) myInstance, & ! current instance of this constitution
myStructure, & ! current lattice structure
ns, & ! short notation for the total number of active slip systems
i, &
j, &
k, &
l, &
t, & ! dislocation type
s, & ! index of my current slip system
sLattice ! index of my current slip system according to lattice order
real(pReal), dimension(3,3,3,3) :: dLp_dTstar3333 ! derivative of Lp with respect to Tstar (3x3x3x3 matrix)
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: &
rhoSgl, & ! single dislocation densities (including used)
v, & ! velocity
dv_dtau ! velocity derivative with respect to the shear stress
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: &
gdotTotal, & ! shear rate
dgdotTotal_dtau ! derivative of the shear rate with respect to the shear stress
!*** initialize local variables
Lp = 0.0_pReal
dLp_dTstar3333 = 0.0_pReal
myInstance = phase_constitutionInstance(material_phase(g,ip,el))
myStructure = constitutive_nonlocal_structure(myInstance)
ns = constitutive_nonlocal_totalNslip(myInstance)
!*** update dislocation velocity
call constitutive_nonlocal_kinetics(v, Tstar_v, Temperature, state, g, ip, el, dv_dtau)
!*** shortcut to state variables
forall (s = 1:ns, t = 1:4) &
rhoSgl(s,t) = max(state%p((t-1)*ns+s), 0.0_pReal)
forall (s = 1:ns, t = 5:8, state%p((t-1)*ns+s) * v(s,t-4) < 0.0_pReal) & ! dead dislocations cntribute instantaneously to slip when sign of velocity changes
rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(state%p((t-1)*ns+s))
!*** Calculation of gdot and its tangent
gdotTotal = sum(rhoSgl * v(1:ns,1:4), 2) * constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance)
dgdotTotal_dtau = sum(rhoSgl * dv_dtau, 2) * constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance)
!*** Calculation of Lp and its tangent
do s = 1,ns
sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance)
Lp = Lp + gdotTotal(s) * lattice_Sslip(1:3,1:3,sLattice,myStructure)
forall (i=1:3,j=1:3,k=1:3,l=1:3) &
dLp_dTstar3333(i,j,k,l) = dLp_dTstar3333(i,j,k,l) + dgdotTotal_dtau(s) * lattice_Sslip(i,j, sLattice,myStructure) &
* lattice_Sslip(k,l, sLattice,myStructure)
enddo
dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333)
#ifndef _OPENMP
if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then
write(6,*)
write(6,'(a,i8,x,i2,x,i1)') '<< CONST >> nonlocal_LpandItsTangent at el ip g ',el,ip,g
write(6,*)
write(6,'(a,/,12(x),12(f12.5,x))') '<< CONST >> gdot total / 1e-3',gdotTotal*1e3_pReal
write(6,'(a,/,3(12(x),3(f12.7,x),/))') '<< CONST >> Lp',Lp
endif
#endif
endsubroutine
!*********************************************************************
!* rate of change of microstructure *
!*********************************************************************
subroutine constitutive_nonlocal_dotState(dotState, Tstar_v, Fe, Fp, Temperature, state, aTolState, timestep, orientation, g,ip,el)
use prec, only: pReal, &
pInt, &
p_vec, &
DAMASK_NaN
use numerics, only: numerics_integrationMode
use IO, only: IO_error
use debug, only: debug_verbosity, &
debug_selectiveDebugger, &
debug_g, &
debug_i, &
debug_e
use math, only: math_norm3, &
math_mul6x6, &
math_mul3x3, &
math_mul33x3, &
math_mul33x33, &
math_inv3x3, &
math_det3x3, &
math_Mandel6to33, &
math_QuaternionDisorientation, &
math_qRot, &
pi
use mesh, only: mesh_NcpElems, &
mesh_maxNips, &
mesh_maxNipNeighbors, &
mesh_element, &
FE_NipNeighbors, &
mesh_ipNeighborhood, &
mesh_ipVolume, &
mesh_ipArea, &
mesh_ipAreaNormal, &
mesh_ipCenterOfGravity
use material, only: homogenization_maxNgrains, &
material_phase, &
phase_constitutionInstance, &
phase_localConstitution, &
phase_constitution
use lattice, only: lattice_Sslip, &
lattice_Sslip_v, &
lattice_sd, &
lattice_sn, &
lattice_st, &
lattice_maxNslipFamily, &
lattice_NslipSystem
use FEsolving, only:theInc, &
FEsolving_execElem, &
FEsolving_execIP
implicit none
!*** input variables
integer(pInt), intent(in) :: g, & ! current grain number
ip, & ! current integration point
el ! current element number
real(pReal), intent(in) :: Temperature, & ! temperature
timestep ! substepped crystallite time increment
real(pReal), dimension(6), intent(in) :: Tstar_v ! current 2nd Piola-Kirchhoff stress in Mandel notation
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
Fe, & ! elastic deformation gradient
Fp ! plastic deformation gradient
real(pReal), dimension(4,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
orientation ! crystal lattice orientation
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
state, & ! current microstructural state
aTolState ! absolute state tolerance
!*** input/output variables
type(p_vec), intent(inout) :: dotState ! evolution of state variables / microstructure
!*** output variables
!*** local variables
integer(pInt) myInstance, & ! current instance of this constitution
myStructure, & ! current lattice structure
ns, & ! short notation for the total number of active slip systems
c, & ! character of dislocation
n, & ! index of my current neighbor
neighboring_el, & ! element number of my neighbor
neighboring_ip, & ! integration point of my neighbor
neighboring_n, & ! neighbor index pointing to me when looking from my neighbor
opposite_n, & ! index of my opposite neighbor
opposite_ip, & ! ip of my opposite neighbor
opposite_el, & ! element index of my opposite neighbor
t, & ! type of dislocation
topp, & ! type of dislocation with opposite sign to t
s, & ! index of my current slip system
sLattice ! index of my current slip system according to lattice order
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),10) :: &
rhoDot, & ! density evolution
rhoDotRemobilization, & ! density evolution by remobilization
rhoDotMultiplication, & ! density evolution by multiplication
rhoDotFlux, & ! density evolution by flux
rhoDotSingle2DipoleGlide, & ! density evolution by dipole formation (by glide)
rhoDotAthermalAnnihilation, & ! density evolution by athermal annihilation
rhoDotThermalAnnihilation ! density evolution by thermal annihilation
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),8) :: &
rhoSgl ! current single dislocation densities (positive/negative screw and edge without dipoles)
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: &
v, & ! dislocation glide velocity
fluxdensity, & ! flux density at central material point
neighboring_fluxdensity, & ! flux density at neighboring material point
gdot ! shear rates
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: &
rhoForest, & ! forest dislocation density
tauThreshold, & ! threshold shear stress
tau, & ! current resolved shear stress
vClimb ! climb velocity of edge dipoles
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),2) :: &
rhoDip, & ! current dipole dislocation densities (screw and edge dipoles)
dLower, & ! minimum stable dipole distance for edges and screws
dUpper ! current maximum stable dipole distance for edges and screws
real(pReal), dimension(3,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: &
m ! direction of dislocation motion
real(pReal), dimension(3,3) :: my_F, & ! my total deformation gradient
neighboring_F, & ! total deformation gradient of my neighbor
my_Fe, & ! my elastic deformation gradient
neighboring_Fe, & ! elastic deformation gradient of my neighbor
Favg ! average total deformation gradient of me and my neighbor
real(pReal), dimension(6) :: Tdislocation_v ! current dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress
real(pReal), dimension(3) :: normal_neighbor2me, & ! interface normal pointing from my neighbor to me in neighbor's lattice configuration
normal_neighbor2me_defConf, & ! interface normal pointing from my neighbor to me in shared deformed configuration
normal_me2neighbor, & ! interface normal pointing from me to my neighbor in my lattice configuration
normal_me2neighbor_defConf ! interface normal pointing from me to my neighbor in shared deformed configuration
real(pReal) area, & ! area of the current interface
transmissivity, & ! overall transmissivity of dislocation flux to neighboring material point
lineLength, & ! dislocation line length leaving the current interface
D ! self diffusion
logical considerEnteringFlux, &
considerLeavingFlux
#ifndef _OPENMP
if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then
write(6,*)
write(6,'(a,i8,x,i2,x,i1)') '<< CONST >> nonlocal_dotState at el ip g ',el,ip,g
write(6,*)
endif
#endif
select case(mesh_element(2,el))
case (1,6,7,8,9)
! all fine
case default
call IO_error(-1,el,ip,g,'element type not supported for nonlocal constitution')
end select
myInstance = phase_constitutionInstance(material_phase(g,ip,el))
myStructure = constitutive_nonlocal_structure(myInstance)
ns = constitutive_nonlocal_totalNslip(myInstance)
tau = 0.0_pReal
gdot = 0.0_pReal
dLower = 0.0_pReal
dUpper = 0.0_pReal
!*** shortcut to state variables
forall (s = 1:ns, t = 1:4) &
rhoSgl(s,t) = max(state(g,ip,el)%p((t-1)*ns+s), 0.0_pReal)
forall (s = 1:ns, t = 5:8) &
rhoSgl(s,t) = state(g,ip,el)%p((t-1)*ns+s)
forall (s = 1:ns, c = 1:2) &
rhoDip(s,c) = max(state(g,ip,el)%p((7+c)*ns+s), 0.0_pReal)
rhoForest = state(g,ip,el)%p(10*ns+1:11*ns)
tauThreshold = state(g,ip,el)%p(11*ns+1:12*ns)
Tdislocation_v = state(g,ip,el)%p(12*ns+1:12*ns+6)
v(1:ns,1) = state(g,ip,el)%p(12*ns+6+1:13*ns+6)
v(1:ns,2) = state(g,ip,el)%p(12*ns+6+1:13*ns+6)
v(1:ns,3) = state(g,ip,el)%p(13*ns+6+1:14*ns+6)
v(1:ns,4) = state(g,ip,el)%p(13*ns+6+1:14*ns+6)
!*** sanity check for timestep
if (timestep <= 0.0_pReal) then ! if illegal timestep...
dotState%p = 0.0_pReal ! ...return without doing anything (-> zero dotState)
return
endif
!****************************************************************************
!*** Calculate shear rate
forall (t = 1:4) &
gdot(1:ns,t) = rhoSgl(1:ns,t) * constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) * v(1:ns,t)
forall (s = 1:ns, t = 1:4, rhoSgl(s,t+4) * v(s,t) < 0.0_pReal) & ! contribution of used rho for changing sign of v
gdot(s,t) = gdot(s,t) + abs(rhoSgl(s,t+4)) * constitutive_nonlocal_burgersPerSlipSystem(s,myInstance) * v(s,t)
#ifndef _OPENMP
if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then
write(6,'(a,/,10(12(x),12(e12.5,x),/))') '<< CONST >> rho / 1/m^2', rhoSgl, rhoDip
write(6,'(a,/,4(12(x),12(e12.5,x),/))') '<< CONST >> gdot / 1/s',gdot
endif
#endif
!****************************************************************************
!*** check CFL condition for flux
if (any(1.2_pReal * abs(v) * timestep > mesh_ipVolume(ip,el) / maxval(mesh_ipArea(:,ip,el)))) then ! safety factor 1.2 (we use the reference volume and are for simplicity here)
#ifndef _OPENMP
if (debug_verbosity > 6) then
write(6,*) '<< CONST >> CFL condition not fullfilled'
endif
#endif
dotState%p = DAMASK_NaN
return
endif
!****************************************************************************
!*** calculate limits for stable dipole height
do s = 1,ns ! loop over slip systems
sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance)
tau(s) = math_mul6x6( Tstar_v + Tdislocation_v, lattice_Sslip_v(1:6,sLattice,myStructure) )
enddo
dLower(1:ns,1) = constitutive_nonlocal_dLowerEdgePerSlipSystem(1:ns,myInstance)
dLower(1:ns,2) = constitutive_nonlocal_dLowerScrewPerSlipSystem(1:ns,myInstance)
dUpper(1:ns,2) = min( 1.0_pReal / sqrt( sum(abs(rhoSgl),2)+sum(rhoDip,2) ), &
constitutive_nonlocal_Gmod(myInstance) * constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) &
/ ( 8.0_pReal * pi * abs(tau) ) )
dUpper(1:ns,1) = dUpper(1:ns,2) / ( 1.0_pReal - constitutive_nonlocal_nu(myInstance) )
!****************************************************************************
!*** dislocation remobilization (bauschinger effect)
rhoDotRemobilization = 0.0_pReal
if (timestep > 0.0_pReal) then
do t = 1,4
do s = 1,ns
if (rhoSgl(s,t+4) * v(s,t) < 0.0_pReal) then
rhoDotRemobilization(s,t) = abs(rhoSgl(s,t+4)) / timestep
rhoSgl(s,t) = rhoSgl(s,t) + abs(rhoSgl(s,t+4))
rhoDotRemobilization(s,t+4) = - rhoSgl(s,t+4) / timestep
rhoSgl(s,t+4) = 0.0_pReal
endif
enddo
enddo
endif
!****************************************************************************
!*** calculate dislocation multiplication
rhoDotMultiplication = 0.0_pReal
where (rhoSgl(1:ns,3:4) > 0.0_pReal) &
rhoDotMultiplication(1:ns,1:2) = spread(0.5_pReal * sum(abs(gdot(1:ns,3:4)),2) * sqrt(rhoForest) &
/ constitutive_nonlocal_lambda0PerSlipSystem(1:ns,myInstance) &
/ constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance), 2, 2)
where (rhoSgl(1:ns,1:2) > 0.0_pReal) &
rhoDotMultiplication(1:ns,3:4) = spread(0.5_pReal * sum(abs(gdot(1:ns,1:2)),2) * sqrt(rhoForest) &
/ constitutive_nonlocal_lambda0PerSlipSystem(1:ns,myInstance) &
/ constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance), 2, 2)
!****************************************************************************
!*** calculate dislocation fluxes (only for nonlocal constitution)
rhoDotFlux = 0.0_pReal
if (.not. phase_localConstitution(material_phase(g,ip,el))) then ! only for nonlocal constitution
!*** take care of the definition of lattice_st = lattice_sd x lattice_sn !!!
!*** opposite sign to our p vector in the (s,p,n) triplet !!!
m(1:3,1:ns,1) = lattice_sd(1:3, constitutive_nonlocal_slipSystemLattice(1:ns,myInstance), myStructure)
m(1:3,1:ns,2) = -lattice_sd(1:3, constitutive_nonlocal_slipSystemLattice(1:ns,myInstance), myStructure)
m(1:3,1:ns,3) = -lattice_st(1:3, constitutive_nonlocal_slipSystemLattice(1:ns,myInstance), myStructure)
m(1:3,1:ns,4) = lattice_st(1:3, constitutive_nonlocal_slipSystemLattice(1:ns,myInstance), myStructure)
my_Fe = Fe(1:3,1:3,g,ip,el)
my_F = math_mul33x33(my_Fe, Fp(1:3,1:3,g,ip,el))
fluxdensity = rhoSgl(1:ns,1:4) * v
do n = 1,FE_NipNeighbors(mesh_element(2,el)) ! loop through my neighbors
neighboring_el = mesh_ipNeighborhood(1,n,ip,el)
neighboring_ip = mesh_ipNeighborhood(2,n,ip,el)
if (neighboring_el > 0_pInt .and. neighboring_ip > 0_pInt) then ! if neighbor exists ...
do neighboring_n = 1,FE_NipNeighbors(mesh_element(2,neighboring_el)) ! find neighboring index that points from my neighbor to myself
if ( el == mesh_ipNeighborhood(1,neighboring_n,neighboring_ip,neighboring_el) &
.and. ip == mesh_ipNeighborhood(2,neighboring_n,neighboring_ip,neighboring_el)) then ! possible candidate
if (math_mul3x3(mesh_ipAreaNormal(1:3,n,ip,el),&
mesh_ipAreaNormal(1:3,neighboring_n,neighboring_ip,neighboring_el)) < 0.0_pReal) then ! area normals have opposite orientation (we have to check that because of special case for single element with two ips and periodicity. In this case the neighbor is identical in two different directions.)
exit
endif
endif
enddo
endif
opposite_n = n + mod(n,2) - mod(n+1,2)
opposite_el = mesh_ipNeighborhood(1,opposite_n,ip,el)
opposite_ip = mesh_ipNeighborhood(2,opposite_n,ip,el)
if (neighboring_el > 0_pInt .and. neighboring_ip > 0_pInt) then ! if neighbor exists, average deformation gradient
neighboring_Fe = Fe(1:3,1:3,g,neighboring_ip,neighboring_el)
neighboring_F = math_mul33x33(neighboring_Fe, Fp(1:3,1:3,g,neighboring_ip,neighboring_el))
Favg = 0.5_pReal * (my_F + neighboring_F)
else ! if no neighbor, take my value as average
Favg = my_F
endif
!* FLUX FROM MY NEIGHBOR TO ME
!* This is only considered, if I have a neighbor of nonlocal constitution (also nonlocal constitutive law with local properties) that is at least a little bit compatible.
!* If it's not at all compatible, no flux is arriving, because everything is dammed in front of my neighbor's interface.
!* The entering flux from my neighbor will be distributed on my slip systems according to the compatibility
considerEnteringFlux = .false.
neighboring_fluxdensity = 0.0_pReal ! needed for check of sign change in flux density below
if (neighboring_el > 0_pInt .or. neighboring_ip > 0_pInt) then
if (phase_constitution(material_phase(1,neighboring_ip,neighboring_el)) == constitutive_nonlocal_label &
.and. any(constitutive_nonlocal_compatibility(:,:,:,n,ip,el) > 0.0_pReal)) &
considerEnteringFlux = .true.
endif
if (considerEnteringFlux) then
neighboring_fluxdensity(1:ns,1) = state(g,neighboring_ip,neighboring_el)%p(1:ns) &
* state(g,neighboring_ip,neighboring_el)%p(12*ns+6+1:13*ns+6)
neighboring_fluxdensity(1:ns,2) = state(g,neighboring_ip,neighboring_el)%p(ns+1:2*ns) &
* state(g,neighboring_ip,neighboring_el)%p(12*ns+6+1:13*ns+6)
neighboring_fluxdensity(1:ns,3) = state(g,neighboring_ip,neighboring_el)%p(2*ns+1:3*ns) &
* state(g,neighboring_ip,neighboring_el)%p(13*ns+6+1:14*ns+6)
neighboring_fluxdensity(1:ns,4) = state(g,neighboring_ip,neighboring_el)%p(3*ns+1:4*ns) &
* state(g,neighboring_ip,neighboring_el)%p(13*ns+6+1:14*ns+6)
normal_neighbor2me_defConf = math_det3x3(Favg) &
* math_mul33x3(math_inv3x3(transpose(Favg)), mesh_ipAreaNormal(1:3,neighboring_n,neighboring_ip,neighboring_el)) ! calculate the normal of the interface in (average) deformed configuration (now pointing from my neighbor to me!!!)
normal_neighbor2me = math_mul33x3(transpose(neighboring_Fe), normal_neighbor2me_defConf) / math_det3x3(neighboring_Fe) ! interface normal in the lattice configuration of my neighbor
area = mesh_ipArea(neighboring_n,neighboring_ip,neighboring_el) * math_norm3(normal_neighbor2me)
normal_neighbor2me = normal_neighbor2me / math_norm3(normal_neighbor2me) ! normalize the surface normal to unit length
do s = 1,ns
do t = 1,4
c = (t + 1) / 2
topp = t + mod(t,2) - mod(t+1,2)
if (neighboring_fluxdensity(s,t) * math_mul3x3(m(1:3,s,t), normal_neighbor2me) > 0.0_pReal & ! flux from my neighbor to me == entering flux for me
.and. fluxdensity(s,t) * neighboring_fluxdensity(s,t) >= 0.0_pReal ) then ! ... only if no sign change in flux density
lineLength = neighboring_fluxdensity(s,t) * math_mul3x3(m(1:3,s,t), normal_neighbor2me) * area ! positive line length that wants to enter through this interface
where (constitutive_nonlocal_compatibility(c,1:ns,s,n,ip,el) > 0.0_pReal) & ! positive compatibility...
rhoDotFlux(1:ns,t) = rhoDotFlux(1:ns,t) + lineLength / mesh_ipVolume(ip,el) & ! ... transferring to equally signed dislocation type
* constitutive_nonlocal_compatibility(c,1:ns,s,n,ip,el) ** 2.0_pReal
where (constitutive_nonlocal_compatibility(c,1:ns,s,n,ip,el) < 0.0_pReal) & ! ..negative compatibility...
rhoDotFlux(1:ns,topp) = rhoDotFlux(1:ns,topp) + lineLength / mesh_ipVolume(ip,el) & ! ... transferring to opposite signed dislocation type
* constitutive_nonlocal_compatibility(c,1:ns,s,n,ip,el) ** 2.0_pReal
endif
enddo
enddo
endif
!* FLUX FROM ME TO MY NEIGHBOR
!* This is not considered, if my opposite neighbor has a different constitutive law than nonlocal (still considered for nonlocal law with lcal properties).
!* Then, we assume, that the opposite(!) neighbor sends an equal amount of dislocations to me.
!* So the net flux in the direction of my neighbor is equal to zero:
!* leaving flux to neighbor == entering flux from opposite neighbor
!* In case of reduced transmissivity, part of the leaving flux is stored as dead dislocation density.
!* That means for an interface of zero transmissivity the leaving flux is fully converted to dead dislocations.
considerLeavingFlux = .true.
if (opposite_el > 0 .and. opposite_ip > 0) then
if (phase_constitution(material_phase(1,opposite_ip,opposite_el)) /= constitutive_nonlocal_label) &
considerLeavingFlux = .false.
endif
if (considerLeavingFlux) then
normal_me2neighbor_defConf = math_det3x3(Favg) * math_mul33x3(math_inv3x3(transpose(Favg)), mesh_ipAreaNormal(1:3,n,ip,el)) ! calculate the normal of the interface in (average) deformed configuration (pointing from me to my neighbor!!!)
normal_me2neighbor = math_mul33x3(transpose(my_Fe), normal_me2neighbor_defConf) / math_det3x3(my_Fe) ! interface normal in my lattice configuration
area = mesh_ipArea(n,ip,el) * math_norm3(normal_me2neighbor)
normal_me2neighbor = normal_me2neighbor / math_norm3(normal_me2neighbor) ! normalize the surface normal to unit length
do s = 1,ns
do t = 1,4
c = (t + 1) / 2
if (fluxdensity(s,t) * math_mul3x3(m(1:3,s,t), normal_me2neighbor) > 0.0_pReal ) then ! flux from me to my neighbor == leaving flux for me (might also be a pure flux from my mobile density to dead density if interface not at all transmissive)
lineLength = fluxdensity(s,t) * math_mul3x3(m(1:3,s,t), normal_me2neighbor) * area ! positive line length that wants to leave through this interface
if (fluxdensity(s,t) * neighboring_fluxdensity(s,t) >= 0.0_pReal) then ! no sign change in flux density
transmissivity = sum(constitutive_nonlocal_compatibility(c,1:ns,s,n,ip,el)**2.0_pReal) ! overall transmissivity from this slip system to my neighbor
else ! sign change in flux density means sign change in stress which does not allow for dislocations to arive at the neighbor
transmissivity = 0.0_pReal
endif
rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / mesh_ipVolume(ip,el) ! subtract dislocation flux from current mobile type
rhoDotFlux(s,t+4) = rhoDotFlux(s,t+4) + lineLength / mesh_ipVolume(ip,el) * (1.0_pReal - transmissivity) &
* sign(1.0_pReal, fluxdensity(s,t)) ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point
endif
enddo
enddo
endif
enddo ! neighbor loop
endif
if (numerics_integrationMode == 1_pInt) then
constitutive_nonlocal_rhoDotFlux(1:ns,1:10,g,ip,el) = rhoDotFlux(1:ns,1:10) ! save flux calculation for output (if in central integration mode)
endif
!****************************************************************************
!*** calculate dipole formation and annihilation
!*** formation by glide
do c = 1,2
rhoDotSingle2DipoleGlide(1:ns,2*c-1) = -2.0_pReal * dUpper(1:ns,c) / constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) &
* (rhoSgl(1:ns,2*c-1) * abs(gdot(1:ns,2*c)) & ! negative mobile --> positive mobile
+ rhoSgl(1:ns,2*c) * abs(gdot(1:ns,2*c-1)) & ! positive mobile --> negative mobile
+ abs(rhoSgl(1:ns,2*c+4)) * abs(gdot(1:ns,2*c-1))) ! positive mobile --> negative immobile
rhoDotSingle2DipoleGlide(1:ns,2*c) = -2.0_pReal * dUpper(1:ns,c) / constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) &
* (rhoSgl(1:ns,2*c-1) * abs(gdot(1:ns,2*c)) & ! negative mobile --> positive mobile
+ rhoSgl(1:ns,2*c) * abs(gdot(1:ns,2*c-1)) & ! positive mobile --> negative mobile
+ abs(rhoSgl(1:ns,2*c+3)) * abs(gdot(1:ns,2*c))) ! negative mobile --> positive immobile
rhoDotSingle2DipoleGlide(1:ns,2*c+3) = -2.0_pReal * dUpper(1:ns,c) / constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) &
* rhoSgl(1:ns,2*c+3) * abs(gdot(1:ns,2*c)) ! negative mobile --> positive immobile
rhoDotSingle2DipoleGlide(1:ns,2*c+4) = -2.0_pReal * dUpper(1:ns,c) / constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) &
* rhoSgl(1:ns,2*c+4) * abs(gdot(1:ns,2*c-1)) ! positive mobile --> negative immobile
rhoDotSingle2DipoleGlide(1:ns,c+8) = - rhoDotSingle2DipoleGlide(1:ns,2*c-1) - rhoDotSingle2DipoleGlide(1:ns,2*c) &
+ abs(rhoDotSingle2DipoleGlide(1:ns,2*c+3)) + abs(rhoDotSingle2DipoleGlide(1:ns,2*c+4))
enddo
!*** athermal annihilation
rhoDotAthermalAnnihilation = 0.0_pReal
forall (c=1:2) &
rhoDotAthermalAnnihilation(1:ns,c+8) = -2.0_pReal * dLower(1:ns,c) / constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) &
* ( 2.0_pReal * (rhoSgl(1:ns,2*c-1) * abs(gdot(1:ns,2*c)) + rhoSgl(1:ns,2*c) * abs(gdot(1:ns,2*c-1))) & ! was single hitting single
+ 2.0_pReal * (abs(rhoSgl(1:ns,2*c+3)) * abs(gdot(1:ns,2*c)) + abs(rhoSgl(1:ns,2*c+4)) * abs(gdot(1:ns,2*c-1))) & ! was single hitting immobile single or was immobile single hit by single
+ rhoDip(1:ns,c) * (abs(gdot(1:ns,2*c-1)) + abs(gdot(1:ns,2*c)))) ! single knocks dipole constituent
!*** thermally activated annihilation of dipoles
rhoDotThermalAnnihilation = 0.0_pReal
D = constitutive_nonlocal_Dsd0(myInstance) * exp(-constitutive_nonlocal_Qsd(myInstance) / (kB * Temperature))
vClimb = constitutive_nonlocal_atomicVolume(myInstance) * D / ( kB * Temperature ) &
* constitutive_nonlocal_Gmod(myInstance) / ( 2.0_pReal * pi * (1.0_pReal-constitutive_nonlocal_nu(myInstance)) ) &
* 2.0_pReal / ( dUpper(1:ns,1) + dLower(1:ns,1) )
rhoDotThermalAnnihilation(1:ns,9) = - 4.0_pReal * rhoDip(1:ns,1) * vClimb / (dUpper(1:ns,1) - dLower(1:ns,1)) ! edge climb
rhoDotThermalAnnihilation(1:ns,10) = 0.0_pReal !!! cross slipping still has to be implemented !!!
!****************************************************************************
!*** assign the rates of dislocation densities to my dotState
rhoDot = 0.0_pReal
rhoDot = rhoDotFlux &
+ rhoDotMultiplication &
+ rhoDotRemobilization &
+ rhoDotSingle2DipoleGlide &
+ rhoDotAthermalAnnihilation &
+ rhoDotThermalAnnihilation
dotState%p(1:10*ns) = dotState%p(1:10*ns) + reshape(rhoDot,(/10*ns/))
#ifndef _OPENMP
if (debug_verbosity > 6 .and. ((debug_e == el .and. debug_i == ip .and. debug_g == g) .or. .not. debug_selectiveDebugger)) then
write(6,'(a,/,8(12(x),12(e12.5,x),/))') '<< CONST >> dislocation remobilization', rhoDotRemobilization(1:ns,1:8) * timestep
write(6,'(a,/,4(12(x),12(e12.5,x),/))') '<< CONST >> dislocation multiplication', rhoDotMultiplication(1:ns,1:4) * timestep
write(6,'(a,/,8(12(x),12(e12.5,x),/))') '<< CONST >> dislocation flux', rhoDotFlux(1:ns,1:8) * timestep
write(6,'(a,/,10(12(x),12(e12.5,x),/))') '<< CONST >> dipole formation by glide', rhoDotSingle2DipoleGlide * timestep
write(6,'(a,/,2(12(x),12(e12.5,x),/))') '<< CONST >> athermal dipole annihilation', &
rhoDotAthermalAnnihilation(1:ns,1:2) * timestep
write(6,'(a,/,2(12(x),12(e12.5,x),/))') '<< CONST >> thermally activated dipole annihilation', &
rhoDotThermalAnnihilation(1:ns,9:10) * timestep
write(6,'(a,/,10(12(x),12(e12.5,x),/))') '<< CONST >> total density change', rhoDot * timestep
write(6,'(a,/,10(12(x),12(f12.7,x),/))') '<< CONST >> relative density change', &
rhoDot(1:ns,1:8) * timestep / (abs(rhoSgl)+1.0e-10), &
rhoDot(1:ns,9:10) * timestep / (rhoDip+1.0e-10)
endif
#endif
endsubroutine
!*********************************************************************
!* COMPATIBILITY UPDATE *
!* Compatibility is defined as normalized product of signed cosine *
!* of the angle between the slip plane normals and signed cosine of *
!* the angle between the slip directions. Only the largest values *
!* that sum up to a total of 1 are considered, all others are set to *
!* zero. *
!*********************************************************************
subroutine constitutive_nonlocal_updateCompatibility(orientation,i,e)
use prec, only: pReal, &
pInt
use math, only: math_QuaternionDisorientation, &
math_mul3x3, &
math_qRot
use material, only: material_phase, &
phase_constitution, &
phase_localConstitution, &
phase_constitutionInstance, &
homogenization_maxNgrains
use mesh, only: mesh_element, &
mesh_ipNeighborhood, &
FE_NipNeighbors, &
FE_maxNipNeighbors, &
mesh_maxNips, &
mesh_NcpElems
use lattice, only: lattice_sn, &
lattice_sd, &
lattice_st
implicit none
!* input variables
integer(pInt), intent(in) :: i, & ! ip index
e ! element index
real(pReal), dimension(4,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
orientation ! crystal orientation in quaternions
!* output variables
!* local variables
integer(pInt) Nneighbors, & ! number of neighbors
n, & ! neighbor index
neighboring_e, & ! element index of my neighbor
neighboring_i, & ! integration point index of my neighbor
my_phase, &
neighboring_phase, &
my_structure, & ! lattice structure
my_instance, & ! instance of constitution
ns, & ! number of active slip systems
s1, & ! slip system index (me)
s2 ! slip system index (my neighbor)
real(pReal), dimension(4) :: absoluteMisorientation ! absolute misorientation (without symmetry) between me and my neighbor
real(pReal), dimension(2,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(1,i,e))),&
constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(1,i,e))),&
FE_NipNeighbors(mesh_element(2,e))) :: &
compatibility ! compatibility for current element and ip
real(pReal), dimension(3,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(1,i,e)))) :: &
slipNormal, &
slipDirection
real(pReal) compatibilitySum, &
thresholdValue, &
nThresholdValues
logical, dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(1,i,e)))) :: &
belowThreshold
Nneighbors = FE_NipNeighbors(mesh_element(2,e))
my_phase = material_phase(1,i,e)
my_instance = phase_constitutionInstance(my_phase)
my_structure = constitutive_nonlocal_structure(my_instance)
ns = constitutive_nonlocal_totalNslip(my_instance)
slipNormal(1:3,1:ns) = lattice_sn(1:3, constitutive_nonlocal_slipSystemLattice(1:ns,my_instance), my_structure)
slipDirection(1:3,1:ns) = lattice_sd(1:3, constitutive_nonlocal_slipSystemLattice(1:ns,my_instance), my_structure)
!*** start out fully compatible
compatibility = 0.0_pReal
forall(s1 = 1:ns) &
compatibility(1:2,s1,s1,1:Nneighbors) = 1.0_pReal
!*** Loop thrugh neighbors and check whether there is any compatibility.
do n = 1,Nneighbors
neighboring_e = mesh_ipNeighborhood(1,n,i,e)
neighboring_i = mesh_ipNeighborhood(2,n,i,e)
!* FREE SURFACE
!* Set surface transmissivity to the value specified in the material.config
if (neighboring_e <= 0 .or. neighboring_i <= 0) then
forall(s1 = 1:ns) &
compatibility(1:2,s1,s1,n) = sqrt(constitutive_nonlocal_surfaceTransmissivity(my_instance))
cycle
endif
!* PHASE BOUNDARY
!* If we encounter a different nonlocal "cpfem" phase at the neighbor,
!* we consider this to be a real "physical" phase boundary, so completely incompatible.
!* If the neighboring "cpfem" phase has a local constitution,
!* we do not consider this to be a phase boundary, so completely compatible.
neighboring_phase = material_phase(1,neighboring_i,neighboring_e)
if (neighboring_phase /= my_phase) then
if (.not. phase_localConstitution(neighboring_phase)) then
forall(s1 = 1:ns) &
compatibility(1:2,s1,s1,n) = 0.0_pReal ! = sqrt(0.0)
endif
cycle
endif
!* GRAIN BOUNDARY ?
!* The compatibility value is defined as the product of the slip normal projection and the slip direction projection.
!* Its sign is always positive for screws, for edges it has the same sign as the slip normal projection.
!* Since the sum for each slip system can easily exceed one (which would result in a transmissivity larger than one),
!* only values above or equal to a certain threshold value are considered. This threshold value is chosen, such that
!* the number of compatible slip systems is minimized with the sum of the original compatibility values exceeding one.
!* Finally the smallest compatibility value is decreased until the sum is exactly equal to one.
!* All values below the threshold are set to zero.
absoluteMisorientation = math_QuaternionDisorientation(orientation(1:4,1,i,e), &
orientation(1:4,1,neighboring_i,neighboring_e), &
0_pInt) ! no symmetry
do s1 = 1,ns ! my slip systems
do s2 = 1,ns ! my neighbor's slip systems
compatibility(1,s2,s1,n) = math_mul3x3(slipNormal(1:3,s1), math_qRot(absoluteMisorientation, slipNormal(1:3,s2))) &
* abs(math_mul3x3(slipDirection(1:3,s1), math_qRot(absoluteMisorientation, slipDirection(1:3,s2))))
compatibility(2,s2,s1,n) = abs(math_mul3x3(slipNormal(1:3,s1), math_qRot(absoluteMisorientation, slipNormal(1:3,s2)))) &
* abs(math_mul3x3(slipDirection(1:3,s1), math_qRot(absoluteMisorientation, slipDirection(1:3,s2))))
enddo
compatibilitySum = 0.0_pReal
belowThreshold = .true.
do while (compatibilitySum < 1.0_pReal .and. any(belowThreshold(1:ns)))
thresholdValue = maxval(compatibility(2,1:ns,s1,n), belowThreshold(1:ns)) ! screws always positive
nThresholdValues = real(count(compatibility(2,1:ns,s1,n) == thresholdValue),pReal)
where (compatibility(2,1:ns,s1,n) >= thresholdValue) &
belowThreshold(1:ns) = .false.
if (compatibilitySum + thresholdValue * nThresholdValues > 1.0_pReal) &
where (abs(compatibility(1:2,1:ns,s1,n)) == thresholdValue) &
compatibility(1:2,1:ns,s1,n) = sign((1.0_pReal - compatibilitySum) / nThresholdValues, compatibility(1:2,1:ns,s1,n))
compatibilitySum = compatibilitySum + nThresholdValues * thresholdValue
enddo
where (belowThreshold(1:ns)) compatibility(1,1:ns,s1,n) = 0.0_pReal
where (belowThreshold(1:ns)) compatibility(2,1:ns,s1,n) = 0.0_pReal
enddo ! my slip systems cycle
enddo ! neighbor cycle
constitutive_nonlocal_compatibility(1:2,1:ns,1:ns,1:Nneighbors,i,e) = compatibility
endsubroutine
!*********************************************************************
!* rate of change of temperature *
!*********************************************************************
pure function constitutive_nonlocal_dotTemperature(Tstar_v,Temperature,state,g,ip,el)
use prec, only: pReal, &
pInt, &
p_vec
use mesh, only: mesh_NcpElems, &
mesh_maxNips
use material, only: homogenization_maxNgrains
implicit none
!* input variables
integer(pInt), intent(in) :: g, & ! current grain ID
ip, & ! current integration point
el ! current element
real(pReal), intent(in) :: Temperature ! temperature
real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola-Kirchhoff stress in Mandel notation
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
state ! microstructural state
!* output variables
real(pReal) constitutive_nonlocal_dotTemperature ! evolution of Temperature
!* local variables
constitutive_nonlocal_dotTemperature = 0.0_pReal
endfunction
!*********************************************************************
!* return array of constitutive results *
!*********************************************************************
function constitutive_nonlocal_postResults(Tstar_v, Fe, Temperature, dt, state, dotState, g,ip,el)
use prec, only: pReal, &
pInt, &
p_vec
use math, only: math_norm3, &
math_mul6x6, &
math_mul3x3, &
math_mul33x3, &
math_mul33x33, &
math_inv3x3, &
math_det3x3, &
math_Mandel6to33, &
math_transpose3x3, &
pi
use mesh, only: mesh_NcpElems, &
mesh_maxNips, &
mesh_element
use material, only: homogenization_maxNgrains, &
material_phase, &
phase_constitutionInstance, &
phase_Noutput
use lattice, only: lattice_Sslip, &
lattice_Sslip_v, &
lattice_sd, &
lattice_st, &
lattice_maxNslipFamily, &
lattice_NslipSystem
implicit none
!*** input variables
integer(pInt), intent(in) :: g, & ! current grain number
ip, & ! current integration point
el ! current element number
real(pReal), intent(in) :: Temperature, & ! temperature
dt ! time increment
real(pReal), dimension(6), intent(in) :: Tstar_v ! current 2nd Piola-Kirchhoff stress in Mandel notation
real(pReal), dimension(3,3), intent(in) :: Fe ! elastic deformation gradient
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
state, & ! current microstructural state
dotState ! evolution rate of microstructural state
!*** output variables
real(pReal), dimension(constitutive_nonlocal_sizePostResults(phase_constitutionInstance(material_phase(g,ip,el)))) :: &
constitutive_nonlocal_postResults
!*** local variables
integer(pInt) myInstance, & ! current instance of this constitution
myStructure, & ! current lattice structure
ns, & ! short notation for the total number of active slip systems
c, & ! character of dislocation
cs, & ! constitutive result index
o, & ! index of current output
t, & ! type of dislocation
s, & ! index of my current slip system
sLattice ! index of my current slip system according to lattice order
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),8) :: &
rhoSgl, & ! current single dislocation densities (positive/negative screw and edge without dipoles)
rhoDotSgl ! evolution rate of single dislocation densities (positive/negative screw and edge without dipoles)
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: &
gdot, & ! shear rates
v ! velocities
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: &
rhoForest, & ! forest dislocation density
tauThreshold, & ! threshold shear stress
tau, & ! current resolved shear stress
vClimb ! climb velocity of edge dipoles
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),2) :: &
rhoDip, & ! current dipole dislocation densities (screw and edge dipoles)
rhoDotDip, & ! evolution rate of dipole dislocation densities (screw and edge dipoles)
dLower, & ! minimum stable dipole distance for edges and screws
dUpper ! current maximum stable dipole distance for edges and screws
real(pReal), dimension(3,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),2) :: &
m, & ! direction of dislocation motion for edge and screw (unit vector)
m_currentconf ! direction of dislocation motion for edge and screw (unit vector) in current configuration
real(pReal), dimension(6) :: Tdislocation_v ! current dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress
real(pReal) D ! self diffusion
real(pReal), dimension(3,3) :: sigma
myInstance = phase_constitutionInstance(material_phase(g,ip,el))
myStructure = constitutive_nonlocal_structure(myInstance)
ns = constitutive_nonlocal_totalNslip(myInstance)
cs = 0_pInt
constitutive_nonlocal_postResults = 0.0_pReal
!* short hand notations for state variables
forall (t = 1:8) rhoSgl(1:ns,t) = state(g,ip,el)%p((t-1)*ns+1:t*ns)
forall (c = 1:2) rhoDip(1:ns,c) = state(g,ip,el)%p((7+c)*ns+1:(8+c)*ns)
rhoForest = state(g,ip,el)%p(10*ns+1:11*ns)
tauThreshold = state(g,ip,el)%p(11*ns+1:12*ns)
Tdislocation_v = state(g,ip,el)%p(12*ns+1:12*ns+6)
forall (t = 1:8) rhoDotSgl(1:ns,t) = dotState(g,ip,el)%p((t-1)*ns+1:t*ns)
forall (c = 1:2) rhoDotDip(1:ns,c) = dotState(g,ip,el)%p((7+c)*ns+1:(8+c)*ns)
forall (t = 1:4) v(1:ns,t) = state(g,ip,el)%p((12+(t-1)/2)*ns+6+1:(13+(t-1)/2)*ns+6)
!* Calculate shear rate
do t = 1,4
do s = 1,ns
if (rhoSgl(s,t+4) * v(s,t) < 0.0_pReal) then
rhoSgl(s,t) = rhoSgl(s,t) + abs(rhoSgl(s,t+4)) ! remobilization of immobile singles for changing sign of v (bauschinger effect)
rhoSgl(s,t+4) = 0.0_pReal ! remobilization of immobile singles for changing sign of v (bauschinger effect)
endif
enddo
enddo
forall (t = 1:4) &
gdot(1:ns,t) = rhoSgl(1:ns,t) * constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) * v(1:ns,t)
!* calculate limits for stable dipole height
do s = 1,ns
sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance)
tau(s) = math_mul6x6(Tstar_v + Tdislocation_v, lattice_Sslip_v(1:6,sLattice,myStructure))
enddo
dLower(1:ns,1) = constitutive_nonlocal_dLowerEdgePerSlipSystem(1:ns,myInstance)
dLower(1:ns,2) = constitutive_nonlocal_dLowerScrewPerSlipSystem(1:ns,myInstance)
dUpper(1:ns,2) = min( constitutive_nonlocal_Gmod(myInstance) * constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) &
/ (8.0_pReal * pi * abs(tau)), &
1.0_pReal / sqrt(sum(abs(rhoSgl),2)+sum(rhoDip,2)) )
dUpper(1:ns,1) = dUpper(1:ns,2) / (1.0_pReal - constitutive_nonlocal_nu(myInstance))
!*** dislocation motion
m(1:3,1:ns,1) = lattice_sd(1:3,constitutive_nonlocal_slipSystemLattice(1:ns,myInstance),myStructure)
m(1:3,1:ns,2) = -lattice_st(1:3,constitutive_nonlocal_slipSystemLattice(1:ns,myInstance),myStructure)
forall (c = 1:2, s = 1:ns) &
m_currentconf(1:3,s,c) = math_mul33x3(Fe, m(1:3,s,c))
do o = 1,phase_Noutput(material_phase(g,ip,el))
select case(constitutive_nonlocal_output(o,myInstance))
case ('rho')
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(abs(rhoSgl),2) + sum(rhoDip,2)
cs = cs + ns
case ('rho_sgl')
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(abs(rhoSgl),2)
cs = cs + ns
case ('rho_sgl_mobile')
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(rhoSgl(1:ns,1:4),2)
cs = cs + ns
case ('rho_sgl_immobile')
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(rhoSgl(1:ns,5:8),2)
cs = cs + ns
case ('rho_dip')
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(rhoDip,2)
cs = cs + ns
case ('rho_edge')
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(abs(rhoSgl(1:ns,(/1,2,5,6/))),2) + rhoDip(1:ns,1)
cs = cs + ns
case ('rho_sgl_edge')
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(abs(rhoSgl(1:ns,(/1,2,5,6/))),2)
cs = cs + ns
case ('rho_sgl_edge_mobile')
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(rhoSgl(1:ns,1:2),2)
cs = cs + ns
case ('rho_sgl_edge_immobile')
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(rhoSgl(1:ns,5:6),2)
cs = cs + ns
case ('rho_sgl_edge_pos')
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,1) + abs(rhoSgl(1:ns,5))
cs = cs + ns
case ('rho_sgl_edge_pos_mobile')
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,1)
cs = cs + ns
case ('rho_sgl_edge_pos_immobile')
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,5)
cs = cs + ns
case ('rho_sgl_edge_neg')
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,2) + abs(rhoSgl(1:ns,6))
cs = cs + ns
case ('rho_sgl_edge_neg_mobile')
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,2)
cs = cs + ns
case ('rho_sgl_edge_neg_immobile')
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,6)
cs = cs + ns
case ('rho_dip_edge')
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoDip(1:ns,1)
cs = cs + ns
case ('rho_screw')
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(abs(rhoSgl(1:ns,(/3,4,7,8/))),2) + rhoDip(1:ns,2)
cs = cs + ns
case ('rho_sgl_screw')
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(abs(rhoSgl(1:ns,(/3,4,7,8/))),2)
cs = cs + ns
case ('rho_sgl_screw_mobile')
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(rhoSgl(1:ns,3:4),2)
cs = cs + ns
case ('rho_sgl_screw_immobile')
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(rhoSgl(1:ns,7:8),2)
cs = cs + ns
case ('rho_sgl_screw_pos')
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,3) + abs(rhoSgl(1:ns,7))
cs = cs + ns
case ('rho_sgl_screw_pos_mobile')
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,3)
cs = cs + ns
case ('rho_sgl_screw_pos_immobile')
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,7)
cs = cs + ns
case ('rho_sgl_screw_neg')
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,4) + abs(rhoSgl(1:ns,8))
cs = cs + ns
case ('rho_sgl_screw_neg_mobile')
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,4)
cs = cs + ns
case ('rho_sgl_screw_neg_immobile')
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,8)
cs = cs + ns
case ('rho_dip_screw')
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoDip(1:ns,2)
cs = cs + ns
case ('excess_rho')
constitutive_nonlocal_postResults(cs+1:cs+ns) = (rhoSgl(1:ns,1) + abs(rhoSgl(1:ns,5))) &
- (rhoSgl(1:ns,2) + abs(rhoSgl(1:ns,6))) &
+ (rhoSgl(1:ns,3) + abs(rhoSgl(1:ns,7))) &
- (rhoSgl(1:ns,4) + abs(rhoSgl(1:ns,8)))
cs = cs + ns
case ('excess_rho_edge')
constitutive_nonlocal_postResults(cs+1:cs+ns) = (rhoSgl(1:ns,1) + abs(rhoSgl(1:ns,5))) &
- (rhoSgl(1:ns,2) + abs(rhoSgl(1:ns,6)))
cs = cs + ns
case ('excess_rho_screw')
constitutive_nonlocal_postResults(cs+1:cs+ns) = (rhoSgl(1:ns,3) + abs(rhoSgl(1:ns,7))) &
- (rhoSgl(1:ns,4) + abs(rhoSgl(1:ns,8)))
cs = cs + ns
case ('rho_forest')
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoForest
cs = cs + ns
case ('delta')
constitutive_nonlocal_postResults(cs+1:cs+ns) = 1.0_pReal / sqrt(sum(abs(rhoSgl),2) + sum(rhoDip,2))
cs = cs + ns
case ('delta_sgl')
constitutive_nonlocal_postResults(cs+1:cs+ns) = 1.0_pReal / sqrt(sum(abs(rhoSgl),2))
cs = cs + ns
case ('delta_dip')
constitutive_nonlocal_postResults(cs+1:cs+ns) = 1.0_pReal / sqrt(sum(rhoDip,2))
cs = cs + ns
case ('shearrate')
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(gdot,2)
cs = cs + ns
case ('resolvedstress')
do s = 1,ns
sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance)
constitutive_nonlocal_postResults(cs+s) = math_mul6x6(Tstar_v + Tdislocation_v, lattice_Sslip_v(1:6,sLattice,myStructure))
enddo
cs = cs + ns
case ('resolvedstress_internal')
do s = 1,ns
sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance)
constitutive_nonlocal_postResults(cs+s) = math_mul6x6(Tdislocation_v, lattice_Sslip_v(1:6,sLattice,myStructure))
enddo
cs = cs + ns
case ('resolvedstress_external')
do s = 1,ns
sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance)
constitutive_nonlocal_postResults(cs+s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,sLattice,myStructure))
enddo
cs = cs + ns
case ('resistance')
constitutive_nonlocal_postResults(cs+1:cs+ns) = tauThreshold
cs = cs + ns
case ('rho_dot')
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(rhoDotSgl,2) + sum(rhoDotDip,2)
cs = cs + ns
case ('rho_dot_sgl')
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(rhoDotSgl,2)
cs = cs + ns
case ('rho_dot_dip')
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(rhoDotDip,2)
cs = cs + ns
case ('rho_dot_gen')
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(abs(gdot),2) * sqrt(rhoForest) &
/ constitutive_nonlocal_lambda0PerSlipSystem(1:ns,myInstance) &
/ constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance)
cs = cs + ns
case ('rho_dot_gen_edge')
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(abs(gdot(1:ns,3:4)),2) * sqrt(rhoForest) &
/ constitutive_nonlocal_lambda0PerSlipSystem(1:ns,myInstance) &
/ constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance)
cs = cs + ns
case ('rho_dot_gen_screw')
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(abs(gdot(1:ns,1:2)),2) * sqrt(rhoForest) &
/ constitutive_nonlocal_lambda0PerSlipSystem(1:ns,myInstance) &
/ constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance)
cs = cs + ns
case ('rho_dot_sgl2dip')
do c=1,2 ! dipole formation by glide
constitutive_nonlocal_postResults(cs+1:cs+ns) = constitutive_nonlocal_postResults(cs+1:cs+ns) + &
2.0_pReal * dUpper(1:ns,c) / constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) &
* ( 2.0_pReal * ( rhoSgl(1:ns,2*c-1) * abs(gdot(1:ns,2*c)) &
+ rhoSgl(1:ns,2*c) * abs(gdot(1:ns,2*c-1))) & ! was single hitting single
+ 2.0_pReal * ( abs(rhoSgl(1:ns,2*c+3)) * abs(gdot(1:ns,2*c)) &
+ abs(rhoSgl(1:ns,2*c+4)) * abs(gdot(1:ns,2*c-1)))) ! was single hitting immobile/used single
enddo
cs = cs + ns
case ('rho_dot_ann_ath')
do c=1,2
constitutive_nonlocal_postResults(cs+1:cs+ns) = constitutive_nonlocal_postResults(cs+1:cs+ns) + &
2.0_pReal * dLower(1:ns,c) / constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) &
* ( 2.0_pReal * ( rhoSgl(1:ns,2*c-1) * abs(gdot(1:ns,2*c)) &
+ rhoSgl(1:ns,2*c) * abs(gdot(1:ns,2*c-1))) & ! was single hitting single
+ 2.0_pReal * ( abs(rhoSgl(1:ns,2*c+3)) * abs(gdot(1:ns,2*c)) &
+ abs(rhoSgl(1:ns,2*c+4)) * abs(gdot(1:ns,2*c-1))) & ! was single hitting immobile/used single
+ rhoDip(1:ns,c) * (abs(gdot(1:ns,2*c-1)) + abs(gdot(1:ns,2*c)))) ! single knocks dipole constituent
enddo
cs = cs + ns
case ('rho_dot_ann_the')
D = constitutive_nonlocal_Dsd0(myInstance) * exp(-constitutive_nonlocal_Qsd(myInstance) / (kB * Temperature))
vClimb = constitutive_nonlocal_atomicVolume(myInstance) * D / (kB * Temperature) &
* constitutive_nonlocal_Gmod(myInstance) / (2.0_pReal * pi * (1.0_pReal-constitutive_nonlocal_nu(myInstance))) &
* 2.0_pReal / (dUpper(1:ns,1) + dLower(1:ns,1))
constitutive_nonlocal_postResults(cs+1:cs+ns) = 4.0_pReal * rhoDip(1:ns,1) * vClimb / (dUpper(1:ns,1) - dLower(1:ns,1))
! !!! cross-slip of screws missing !!!
cs = cs + ns
case ('rho_dot_flux')
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(constitutive_nonlocal_rhoDotFlux(1:ns,1:4,g,ip,el),2) &
+ sum(abs(constitutive_nonlocal_rhoDotFlux(1:ns,5:8,g,ip,el)),2)
cs = cs + ns
case ('rho_dot_flux_edge')
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(constitutive_nonlocal_rhoDotFlux(1:ns,1:2,g,ip,el),2) &
+ sum(abs(constitutive_nonlocal_rhoDotFlux(1:ns,5:6,g,ip,el)),2)
cs = cs + ns
case ('rho_dot_flux_screw')
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(constitutive_nonlocal_rhoDotFlux(1:ns,3:4,g,ip,el),2) &
+ sum(abs(constitutive_nonlocal_rhoDotFlux(1:ns,7:8,g,ip,el)),2)
cs = cs + ns
case ('dislocationvelocity')
constitutive_nonlocal_postResults(cs+1:cs+ns) = v(1:ns,1)
cs = cs + ns
case ('fluxdensity_edge_pos_x')
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,1) * v(1:ns,1) * m_currentconf(1,1:ns,1)
cs = cs + ns
case ('fluxdensity_edge_pos_y')
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,1) * v(1:ns,1) * m_currentconf(2,1:ns,1)
cs = cs + ns
case ('fluxdensity_edge_pos_z')
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,1) * v(1:ns,1) * m_currentconf(3,1:ns,1)
cs = cs + ns
case ('fluxdensity_edge_neg_x')
constitutive_nonlocal_postResults(cs+1:cs+ns) = - rhoSgl(1:ns,2) * v(1:ns,2) * m_currentconf(1,1:ns,1)
cs = cs + ns
case ('fluxdensity_edge_neg_y')
constitutive_nonlocal_postResults(cs+1:cs+ns) = - rhoSgl(1:ns,2) * v(1:ns,2) * m_currentconf(2,1:ns,1)
cs = cs + ns
case ('fluxdensity_edge_neg_z')
constitutive_nonlocal_postResults(cs+1:cs+ns) = - rhoSgl(1:ns,2) * v(1:ns,2) * m_currentconf(3,1:ns,1)
cs = cs + ns
case ('fluxdensity_screw_pos_x')
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,3) * v(1:ns,3) * m_currentconf(1,1:ns,2)
cs = cs + ns
case ('fluxdensity_screw_pos_y')
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,3) * v(1:ns,3) * m_currentconf(2,1:ns,2)
cs = cs + ns
case ('fluxdensity_screw_pos_z')
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,3) * v(1:ns,3) * m_currentconf(3,1:ns,2)
cs = cs + ns
case ('fluxdensity_screw_neg_x')
constitutive_nonlocal_postResults(cs+1:cs+ns) = - rhoSgl(1:ns,4) * v(1:ns,4) * m_currentconf(1,1:ns,2)
cs = cs + ns
case ('fluxdensity_screw_neg_y')
constitutive_nonlocal_postResults(cs+1:cs+ns) = - rhoSgl(1:ns,4) * v(1:ns,4) * m_currentconf(2,1:ns,2)
cs = cs + ns
case ('fluxdensity_screw_neg_z')
constitutive_nonlocal_postResults(cs+1:cs+ns) = - rhoSgl(1:ns,4) * v(1:ns,4) * m_currentconf(3,1:ns,2)
cs = cs + ns
case ('d_upper_edge')
constitutive_nonlocal_postResults(cs+1:cs+ns) = dUpper(1:ns,1)
cs = cs + ns
case ('d_upper_screw')
constitutive_nonlocal_postResults(cs+1:cs+ns) = dUpper(1:ns,2)
cs = cs + ns
case('internalstress')
sigma = math_mul33x33(Fe, math_mul33x33(math_Mandel6to33(Tdislocation_v), math_transpose3x3(Fe))) / math_det3x3(Fe)
constitutive_nonlocal_postResults(cs+1) = sigma(1,1)
constitutive_nonlocal_postResults(cs+2) = sigma(2,2)
constitutive_nonlocal_postResults(cs+3) = sigma(3,3)
constitutive_nonlocal_postResults(cs+4) = sigma(1,2)
constitutive_nonlocal_postResults(cs+5) = sigma(2,3)
constitutive_nonlocal_postResults(cs+6) = sigma(3,1)
cs = cs + 6_pInt
end select
enddo
endfunction
END MODULE
Reference in New Issue View Git Blame Copy Permalink