3752 lines
189 KiB
Fortran
3752 lines
189 KiB
Fortran
! Copyright 2011-13 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, &
|
|
p_vec
|
|
|
|
implicit none
|
|
private
|
|
|
|
|
|
!* Definition of parameters
|
|
|
|
character (len=*), parameter, public :: &
|
|
CONSTITUTIVE_NONLOCAL_LABEL = 'nonlocal'
|
|
|
|
character(len=22), dimension(11), parameter, private :: &
|
|
BASICSTATES = (/'rhoSglEdgePosMobile ', &
|
|
'rhoSglEdgeNegMobile ', &
|
|
'rhoSglScrewPosMobile ', &
|
|
'rhoSglScrewNegMobile ', &
|
|
'rhoSglEdgePosImmobile ', &
|
|
'rhoSglEdgeNegImmobile ', &
|
|
'rhoSglScrewPosImmobile', &
|
|
'rhoSglScrewNegImmobile', &
|
|
'rhoDipEdge ', &
|
|
'rhoDipScrew ', &
|
|
'accumulatedshear ' /) !< list of "basic" microstructural state variables that are independent from other state variables
|
|
|
|
character(len=16), dimension(3), parameter, private :: &
|
|
DEPENDENTSTATES = (/'rhoForest ', &
|
|
'tauThreshold ', &
|
|
'tauBack ' /) !< list of microstructural state variables that depend on other state variables
|
|
|
|
character(len=20), dimension(6), parameter, private :: &
|
|
OTHERSTATES = (/'velocityEdgePos ', &
|
|
'velocityEdgeNeg ', &
|
|
'velocityScrewPos ', &
|
|
'velocityScrewNeg ', &
|
|
'maxDipoleHeightEdge ', &
|
|
'maxDipoleHeightScrew' /) !< list of other dependent state variables that are not updated by microstructure
|
|
|
|
real(pReal), parameter, private :: &
|
|
KB = 1.38e-23_pReal !< Physical parameter, Boltzmann constant in J/Kelvin
|
|
|
|
|
|
!* Definition of global variables
|
|
|
|
integer(pInt), dimension(:), allocatable, public, protected :: &
|
|
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, public :: &
|
|
constitutive_nonlocal_sizePostResult !< size of each post result output
|
|
|
|
character(len=64), dimension(:,:), allocatable, target, public :: &
|
|
constitutive_nonlocal_output !< name of each post result output
|
|
|
|
integer(pInt), dimension(:), allocatable, private :: &
|
|
Noutput !< number of outputs per instance of this plasticity
|
|
|
|
integer(pInt), dimension(:,:), allocatable, private :: &
|
|
iGamma, & !< state indices for accumulated shear
|
|
iRhoF, & !< state indices for forest density
|
|
iTauF, & !< state indices for critical resolved shear stress
|
|
iTauB !< state indices for backstress
|
|
integer(pInt), dimension(:,:,:), allocatable, private :: &
|
|
iRhoU, & !< state indices for unblocked density
|
|
iRhoB, & !< state indices for blocked density
|
|
iRhoD, & !< state indices for dipole density
|
|
iV, & !< state indices for dislcation velocities
|
|
iD !< state indices for stable dipole height
|
|
|
|
|
|
character(len=32), dimension(:), allocatable, public :: &
|
|
constitutive_nonlocal_structureName !< name of the lattice structure
|
|
|
|
integer(pInt), dimension(:), allocatable, public :: &
|
|
constitutive_nonlocal_structure !< number representing the kind of lattice structure
|
|
|
|
integer(pInt), dimension(:), allocatable, private :: &
|
|
totalNslip !< total number of active slip systems for each instance
|
|
|
|
integer(pInt), dimension(:,:), allocatable, private :: &
|
|
Nslip, & !< number of active slip systems for each family and instance
|
|
slipFamily, & !< lookup table relating active slip system to slip family for each instance
|
|
slipSystemLattice, & !< lookup table relating active slip system index to lattice slip system index for each instance
|
|
colinearSystem !< colinear system to the active slip system (only valid for fcc!)
|
|
|
|
real(pReal), dimension(:), allocatable, private :: &
|
|
CoverA, & !< c/a ratio for hex type lattice
|
|
mu, & !< shear modulus
|
|
nu, & !< poisson's ratio
|
|
atomicVolume, & !< atomic volume
|
|
Dsd0, & !< prefactor for self-diffusion coefficient
|
|
selfDiffusionEnergy, & !< activation enthalpy for diffusion
|
|
aTolRho, & !< absolute tolerance for dislocation density in state integration
|
|
aTolShear, & !< absolute tolerance for accumulated shear in state integration
|
|
significantRho, & !< density considered significant
|
|
significantN, & !< number of dislocations considered significant
|
|
cutoffRadius, & !< cutoff radius for dislocation stress
|
|
doublekinkwidth, & !< width of a doubkle kink in multiples of the burgers vector length b
|
|
solidSolutionEnergy, & !< activation energy for solid solution in J
|
|
solidSolutionSize, & !< solid solution obstacle size in multiples of the burgers vector length
|
|
solidSolutionConcentration, & !< concentration of solid solution in atomic parts
|
|
pParam, & !< parameter for kinetic law (Kocks,Argon,Ashby)
|
|
qParam, & !< parameter for kinetic law (Kocks,Argon,Ashby)
|
|
viscosity, & !< viscosity for dislocation glide in Pa s
|
|
fattack, & !< attack frequency in Hz
|
|
rhoSglScatter, & !< standard deviation of scatter in initial dislocation density
|
|
surfaceTransmissivity, & !< transmissivity at free surface
|
|
grainboundaryTransmissivity, & !< transmissivity at grain boundary (identified by different texture)
|
|
CFLfactor, & !< safety factor for CFL flux condition
|
|
fEdgeMultiplication, & !< factor that determines how much edge dislocations contribute to multiplication (0...1)
|
|
rhoSglRandom, &
|
|
rhoSglRandomBinning, &
|
|
linetensionEffect, &
|
|
edgeJogFactor
|
|
|
|
real(pReal), dimension(:,:), allocatable, private :: &
|
|
rhoSglEdgePos0, & !< initial edge_pos dislocation density per slip system for each family and instance
|
|
rhoSglEdgeNeg0, & !< initial edge_neg dislocation density per slip system for each family and instance
|
|
rhoSglScrewPos0, & !< initial screw_pos dislocation density per slip system for each family and instance
|
|
rhoSglScrewNeg0, & !< initial screw_neg dislocation density per slip system for each family and instance
|
|
rhoDipEdge0, & !< initial edge dipole dislocation density per slip system for each family and instance
|
|
rhoDipScrew0, & !< initial screw dipole dislocation density per slip system for each family and instance
|
|
lambda0PerSlipFamily, & !< mean free path prefactor for each family and instance
|
|
lambda0, & !< mean free path prefactor for each slip system and instance
|
|
burgersPerSlipFamily, & !< absolute length of burgers vector [m] for each family and instance
|
|
burgers, & !< absolute length of burgers vector [m] for each slip system and instance
|
|
interactionSlipSlip !< coefficients for slip-slip interaction for each interaction type and instance
|
|
|
|
real(pReal), dimension(:,:,:), allocatable, private :: &
|
|
Cslip66, & !< elasticity matrix in Mandel notation for each instance
|
|
minDipoleHeightPerSlipFamily, & !< minimum stable edge/screw dipole height for each family and instance
|
|
minDipoleHeight, & !< minimum stable edge/screw dipole height for each slip system and instance
|
|
peierlsStressPerSlipFamily, & !< Peierls stress (edge and screw)
|
|
peierlsStress, & !< Peierls stress (edge and screw)
|
|
forestProjectionEdge, & !< matrix of forest projections of edge dislocations for each instance
|
|
forestProjectionScrew, & !< matrix of forest projections of screw dislocations for each instance
|
|
interactionMatrixSlipSlip !< interaction matrix of the different slip systems for each instance
|
|
|
|
real(pReal), dimension(:,:,:,:), allocatable, private :: &
|
|
lattice2slip, & !< orthogonal transformation matrix from lattice coordinate system to slip coordinate system (passive rotation !!!)
|
|
rhoDotEdgeJogsOutput, &
|
|
sourceProbability
|
|
|
|
real(pReal), dimension(:,:,:,:,:), allocatable, private :: &
|
|
Cslip3333, & !< elasticity matrix for each instance
|
|
rhoDotFluxOutput, &
|
|
rhoDotMultiplicationOutput, &
|
|
rhoDotSingle2DipoleGlideOutput, &
|
|
rhoDotAthermalAnnihilationOutput, &
|
|
rhoDotThermalAnnihilationOutput, &
|
|
nonSchmidProjection !< combined projection of Schmid and non-Schmid contributions to the resolved shear stress (only for screws)
|
|
|
|
real(pReal), dimension(:,:,:,:,:,:), allocatable, private :: &
|
|
compatibility !< slip system compatibility between me and my neighbors
|
|
|
|
real(pReal), dimension(:,:), allocatable, private :: &
|
|
nonSchmidCoeff
|
|
|
|
logical, dimension(:), allocatable, private :: &
|
|
shortRangeStressCorrection, & !< flag indicating the use of the short range stress correction by a excess density gradient term
|
|
probabilisticMultiplication
|
|
|
|
public :: &
|
|
constitutive_nonlocal_init, &
|
|
constitutive_nonlocal_stateInit, &
|
|
constitutive_nonlocal_aTolState, &
|
|
constitutive_nonlocal_homogenizedC, &
|
|
constitutive_nonlocal_microstructure, &
|
|
constitutive_nonlocal_LpAndItsTangent, &
|
|
constitutive_nonlocal_dotState, &
|
|
constitutive_nonlocal_deltaState, &
|
|
constitutive_nonlocal_updateCompatibility, &
|
|
constitutive_nonlocal_postResults
|
|
|
|
private :: &
|
|
constitutive_nonlocal_kinetics, &
|
|
constitutive_nonlocal_dislocationstress
|
|
|
|
|
|
CONTAINS
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
!> @brief module initialization
|
|
!> @details reads in material parameters, allocates arrays, and does sanity checks
|
|
!--------------------------------------------------------------------------------------------------
|
|
subroutine constitutive_nonlocal_init(myFile)
|
|
|
|
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
|
|
use math, only: math_Mandel3333to66, &
|
|
math_Voigt66to3333, &
|
|
math_mul3x3, &
|
|
math_transpose33
|
|
use IO, only: IO_read, &
|
|
IO_lc, &
|
|
IO_getTag, &
|
|
IO_isBlank, &
|
|
IO_stringPos, &
|
|
IO_stringValue, &
|
|
IO_floatValue, &
|
|
IO_intValue, &
|
|
IO_error, &
|
|
IO_warning, &
|
|
IO_timeStamp
|
|
use debug, only: debug_level, &
|
|
debug_constitutive, &
|
|
debug_levelBasic
|
|
use mesh, only: mesh_NcpElems, &
|
|
mesh_maxNips, &
|
|
mesh_maxNipNeighbors
|
|
use material, only: homogenization_maxNgrains, &
|
|
phase_plasticity, &
|
|
phase_plasticityInstance, &
|
|
phase_Noutput
|
|
use lattice
|
|
|
|
!*** output variables
|
|
|
|
!*** input variables
|
|
integer(pInt), intent(in) :: myFile
|
|
|
|
!*** local variables
|
|
integer(pInt), parameter :: MAXNCHUNKS = LATTICE_maxNinteraction + 1_pInt
|
|
integer(pInt), &
|
|
dimension(1_pInt+2_pInt*MAXNCHUNKS) :: positions
|
|
integer(pInt), dimension(7) :: configNchunks
|
|
integer(pInt) :: section = 0_pInt, &
|
|
maxNinstance, &
|
|
maxTotalNslip, &
|
|
structID, &
|
|
f, & ! index of my slip family
|
|
i, & ! index of my instance of this plasticity
|
|
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
|
|
t, & ! index of dislocation type
|
|
c, & ! index of dislocation character
|
|
Nchunks_SlipSlip = 0_pInt, &
|
|
Nchunks_SlipFamilies = 0_pInt, &
|
|
Nchunks_nonSchmid = 0_pInt, &
|
|
mySize = 0_pInt ! to suppress warnings, safe as init is called only once
|
|
character(len=65536) tag
|
|
character(len=65536) :: line = '' ! to start initialized
|
|
|
|
write(6,'(/,a)') ' <<<+- constitutive_'//CONSTITUTIVE_NONLOCAL_label//' init -+>>>'
|
|
write(6,'(a)') ' $Id$'
|
|
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
|
|
#include "compilation_info.f90"
|
|
|
|
maxNinstance = int(count(phase_plasticity == CONSTITUTIVE_NONLOCAL_LABEL),pInt)
|
|
if (maxNinstance == 0) return ! we don't have to do anything if there's no instance for this constitutive law
|
|
|
|
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
|
|
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
|
|
|
|
!*** memory 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))
|
|
allocate(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 = ''
|
|
Noutput = 0_pInt
|
|
|
|
allocate(constitutive_nonlocal_structureName(maxNinstance))
|
|
allocate(constitutive_nonlocal_structure(maxNinstance))
|
|
allocate(Nslip(lattice_maxNslipFamily, maxNinstance))
|
|
allocate(slipFamily(lattice_maxNslip, maxNinstance))
|
|
allocate(slipSystemLattice(lattice_maxNslip, maxNinstance))
|
|
allocate(totalNslip(maxNinstance))
|
|
constitutive_nonlocal_structureName = ''
|
|
constitutive_nonlocal_structure = 0_pInt
|
|
Nslip = 0_pInt
|
|
slipFamily = 0_pInt
|
|
slipSystemLattice = 0_pInt
|
|
totalNslip = 0_pInt
|
|
|
|
allocate(CoverA(maxNinstance))
|
|
allocate(mu(maxNinstance))
|
|
allocate(nu(maxNinstance))
|
|
allocate(atomicVolume(maxNinstance))
|
|
allocate(Dsd0(maxNinstance))
|
|
allocate(selfDiffusionEnergy(maxNinstance))
|
|
allocate(aTolRho(maxNinstance))
|
|
allocate(aTolShear(maxNinstance))
|
|
allocate(significantRho(maxNinstance))
|
|
allocate(significantN(maxNinstance))
|
|
allocate(Cslip66(6,6,maxNinstance))
|
|
allocate(Cslip3333(3,3,3,3,maxNinstance))
|
|
allocate(cutoffRadius(maxNinstance))
|
|
allocate(doublekinkwidth(maxNinstance))
|
|
allocate(solidSolutionEnergy(maxNinstance))
|
|
allocate(solidSolutionSize(maxNinstance))
|
|
allocate(solidSolutionConcentration(maxNinstance))
|
|
allocate(pParam(maxNinstance))
|
|
allocate(qParam(maxNinstance))
|
|
allocate(viscosity(maxNinstance))
|
|
allocate(fattack(maxNinstance))
|
|
allocate(rhoSglScatter(maxNinstance))
|
|
allocate(rhoSglRandom(maxNinstance))
|
|
allocate(rhoSglRandomBinning(maxNinstance))
|
|
allocate(surfaceTransmissivity(maxNinstance))
|
|
allocate(grainboundaryTransmissivity(maxNinstance))
|
|
allocate(shortRangeStressCorrection(maxNinstance))
|
|
allocate(probabilisticMultiplication(maxNinstance))
|
|
allocate(CFLfactor(maxNinstance))
|
|
allocate(fEdgeMultiplication(maxNinstance))
|
|
allocate(linetensionEffect(maxNinstance))
|
|
allocate(edgeJogFactor(maxNinstance))
|
|
CoverA = 0.0_pReal
|
|
mu = 0.0_pReal
|
|
atomicVolume = 0.0_pReal
|
|
Dsd0 = -1.0_pReal
|
|
selfDiffusionEnergy = 0.0_pReal
|
|
aTolRho = 0.0_pReal
|
|
aTolShear = 0.0_pReal
|
|
significantRho = 0.0_pReal
|
|
significantN = 0.0_pReal
|
|
nu = 0.0_pReal
|
|
Cslip66 = 0.0_pReal
|
|
Cslip3333 = 0.0_pReal
|
|
cutoffRadius = -1.0_pReal
|
|
doublekinkwidth = 0.0_pReal
|
|
solidSolutionEnergy = 0.0_pReal
|
|
solidSolutionSize = 0.0_pReal
|
|
solidSolutionConcentration = 0.0_pReal
|
|
pParam = 1.0_pReal
|
|
qParam = 1.0_pReal
|
|
viscosity = 0.0_pReal
|
|
fattack = 0.0_pReal
|
|
rhoSglScatter = 0.0_pReal
|
|
rhoSglRandom = 0.0_pReal
|
|
rhoSglRandomBinning = 1.0_pReal
|
|
surfaceTransmissivity = 1.0_pReal
|
|
grainboundaryTransmissivity = -1.0_pReal
|
|
CFLfactor = 2.0_pReal
|
|
fEdgeMultiplication = 0.0_pReal
|
|
linetensionEffect = 0.0_pReal
|
|
edgeJogFactor = 1.0_pReal
|
|
shortRangeStressCorrection = .false.
|
|
probabilisticMultiplication = .false.
|
|
|
|
allocate(rhoSglEdgePos0(lattice_maxNslipFamily,maxNinstance))
|
|
allocate(rhoSglEdgeNeg0(lattice_maxNslipFamily,maxNinstance))
|
|
allocate(rhoSglScrewPos0(lattice_maxNslipFamily,maxNinstance))
|
|
allocate(rhoSglScrewNeg0(lattice_maxNslipFamily,maxNinstance))
|
|
allocate(rhoDipEdge0(lattice_maxNslipFamily,maxNinstance))
|
|
allocate(rhoDipScrew0(lattice_maxNslipFamily,maxNinstance))
|
|
allocate(burgersPerSlipFamily(lattice_maxNslipFamily,maxNinstance))
|
|
allocate(lambda0PerSlipFamily(lattice_maxNslipFamily,maxNinstance))
|
|
allocate(interactionSlipSlip(lattice_maxNinteraction,maxNinstance))
|
|
rhoSglEdgePos0 = -1.0_pReal
|
|
rhoSglEdgeNeg0 = -1.0_pReal
|
|
rhoSglScrewPos0 = -1.0_pReal
|
|
rhoSglScrewNeg0 = -1.0_pReal
|
|
rhoDipEdge0 = -1.0_pReal
|
|
rhoDipScrew0 = -1.0_pReal
|
|
burgersPerSlipFamily = 0.0_pReal
|
|
lambda0PerSlipFamily = 0.0_pReal
|
|
interactionSlipSlip = 0.0_pReal
|
|
|
|
allocate(minDipoleHeightPerSlipFamily(lattice_maxNslipFamily,2,maxNinstance))
|
|
allocate(peierlsStressPerSlipFamily(lattice_maxNslipFamily,2,maxNinstance))
|
|
minDipoleHeightPerSlipFamily = -1.0_pReal
|
|
peierlsStressPerSlipFamily = 0.0_pReal
|
|
|
|
allocate(nonSchmidCoeff(lattice_maxNnonSchmid,maxNinstance))
|
|
nonSchmidCoeff = 0.0_pReal
|
|
|
|
!*** readout data from material.config file
|
|
|
|
rewind(myFile)
|
|
do while (trim(line) /= '#EOF#' .and. IO_lc(IO_getTag(line,'<','>')) /= 'phase') ! wind forward to <phase>
|
|
line = IO_read(myFile)
|
|
enddo
|
|
|
|
do while (trim(line) /= '#EOF#') ! read thru sections of phase part
|
|
line = IO_read(myFile)
|
|
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_pInt ! advance section counter
|
|
cycle
|
|
endif
|
|
if (section > 0_pInt ) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
|
|
if (trim(phase_plasticity(section)) == CONSTITUTIVE_NONLOCAL_LABEL) then ! one of my sections
|
|
i = phase_plasticityInstance(section) ! which instance of my plasticity is present phase
|
|
positions = IO_stringPos(line,MAXNCHUNKS)
|
|
tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key
|
|
select case(tag)
|
|
case('plasticity','elasticity','/nonlocal/')
|
|
cycle
|
|
case ('(output)')
|
|
Noutput(i) = Noutput(i) + 1_pInt
|
|
constitutive_nonlocal_output(Noutput(i),i) = IO_lc(IO_stringValue(line,positions,2_pInt))
|
|
case ('lattice_structure')
|
|
constitutive_nonlocal_structureName(i) = IO_lc(IO_stringValue(line,positions,2_pInt))
|
|
configNchunks = lattice_configNchunks(constitutive_nonlocal_structureName(i))
|
|
Nchunks_SlipFamilies = configNchunks(1)
|
|
Nchunks_SlipSlip = configNchunks(3)
|
|
Nchunks_nonSchmid = configNchunks(7)
|
|
case ('c/a_ratio','covera_ratio')
|
|
CoverA(i) = IO_floatValue(line,positions,2_pInt)
|
|
case ('c11')
|
|
Cslip66(1,1,i) = IO_floatValue(line,positions,2_pInt)
|
|
case ('c12')
|
|
Cslip66(1,2,i) = IO_floatValue(line,positions,2_pInt)
|
|
case ('c13')
|
|
Cslip66(1,3,i) = IO_floatValue(line,positions,2_pInt)
|
|
case ('c22')
|
|
Cslip66(2,2,i) = IO_floatValue(line,positions,2_pInt)
|
|
case ('c23')
|
|
Cslip66(2,3,i) = IO_floatValue(line,positions,2_pInt)
|
|
case ('c33')
|
|
Cslip66(3,3,i) = IO_floatValue(line,positions,2_pInt)
|
|
case ('c44')
|
|
Cslip66(4,4,i) = IO_floatValue(line,positions,2_pInt)
|
|
case ('c55')
|
|
Cslip66(5,5,i) = IO_floatValue(line,positions,2_pInt)
|
|
case ('c66')
|
|
Cslip66(6,6,i) = IO_floatValue(line,positions,2_pInt)
|
|
case ('nslip')
|
|
if (positions(1) < 1_pInt + Nchunks_SlipFamilies) &
|
|
call IO_warning(50_pInt,ext_msg=trim(tag)//' ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
Nchunks_SlipFamilies = positions(1) - 1_pInt
|
|
do f = 1_pInt, Nchunks_SlipFamilies
|
|
Nslip(f,i) = IO_intValue(line,positions,1_pInt+f)
|
|
enddo
|
|
case ('rhosgledgepos0')
|
|
do f = 1_pInt, Nchunks_SlipFamilies
|
|
rhoSglEdgePos0(f,i) = IO_floatValue(line,positions,1_pInt+f)
|
|
enddo
|
|
case ('rhosgledgeneg0')
|
|
do f = 1_pInt, Nchunks_SlipFamilies
|
|
rhoSglEdgeNeg0(f,i) = IO_floatValue(line,positions,1_pInt+f)
|
|
enddo
|
|
case ('rhosglscrewpos0')
|
|
do f = 1_pInt, Nchunks_SlipFamilies
|
|
rhoSglScrewPos0(f,i) = IO_floatValue(line,positions,1_pInt+f)
|
|
enddo
|
|
case ('rhosglscrewneg0')
|
|
do f = 1_pInt, Nchunks_SlipFamilies
|
|
rhoSglScrewNeg0(f,i) = IO_floatValue(line,positions,1_pInt+f)
|
|
enddo
|
|
case ('rhodipedge0')
|
|
do f = 1_pInt, Nchunks_SlipFamilies
|
|
rhoDipEdge0(f,i) = IO_floatValue(line,positions,1_pInt+f)
|
|
enddo
|
|
case ('rhodipscrew0')
|
|
do f = 1_pInt, Nchunks_SlipFamilies
|
|
rhoDipScrew0(f,i) = IO_floatValue(line,positions,1_pInt+f)
|
|
enddo
|
|
case ('lambda0')
|
|
do f = 1_pInt, Nchunks_SlipFamilies
|
|
lambda0PerSlipFamily(f,i) = IO_floatValue(line,positions,1_pInt+f)
|
|
enddo
|
|
case ('burgers')
|
|
do f = 1_pInt, Nchunks_SlipFamilies
|
|
burgersPerSlipFamily(f,i) = IO_floatValue(line,positions,1_pInt+f)
|
|
enddo
|
|
case('cutoffradius','r')
|
|
cutoffRadius(i) = IO_floatValue(line,positions,2_pInt)
|
|
case('minimumdipoleheightedge','ddipminedge')
|
|
do f = 1_pInt, Nchunks_SlipFamilies
|
|
minDipoleHeightPerSlipFamily(f,1_pInt,i) = IO_floatValue(line,positions,1_pInt+f)
|
|
enddo
|
|
case('minimumdipoleheightscrew','ddipminscrew')
|
|
do f = 1_pInt, Nchunks_SlipFamilies
|
|
minDipoleHeightPerSlipFamily(f,2_pInt,i) = IO_floatValue(line,positions,1_pInt+f)
|
|
enddo
|
|
case('atomicvolume')
|
|
atomicVolume(i) = IO_floatValue(line,positions,2_pInt)
|
|
case('selfdiffusionprefactor','dsd0')
|
|
Dsd0(i) = IO_floatValue(line,positions,2_pInt)
|
|
case('selfdiffusionenergy','qsd')
|
|
selfDiffusionEnergy(i) = IO_floatValue(line,positions,2_pInt)
|
|
case('atol_rho','atol_density','absolutetolerancedensity','absolutetolerance_density')
|
|
aTolRho(i) = IO_floatValue(line,positions,2_pInt)
|
|
case('atol_shear','atol_plasticshear','atol_accumulatedshear','absolutetoleranceshear','absolutetolerance_shear')
|
|
aTolShear(i) = IO_floatValue(line,positions,2_pInt)
|
|
case('significantrho','significant_rho','significantdensity','significant_density')
|
|
significantRho(i) = IO_floatValue(line,positions,2_pInt)
|
|
case('significantn','significant_n','significantdislocations','significant_dislcations')
|
|
significantN(i) = IO_floatValue(line,positions,2_pInt)
|
|
case ('interaction_slipslip')
|
|
if (positions(1) < 1_pInt + Nchunks_SlipSlip) &
|
|
call IO_warning(52_pInt,ext_msg=trim(tag)//' ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
do it = 1_pInt,Nchunks_SlipSlip
|
|
interactionSlipSlip(it,i) = IO_floatValue(line,positions,1_pInt+it)
|
|
enddo
|
|
case('linetension','linetensioneffect','linetension_effect')
|
|
linetensionEffect(i) = IO_floatValue(line,positions,2_pInt)
|
|
case('edgejog','edgejogs','edgejogeffect','edgejog_effect')
|
|
edgeJogFactor(i) = IO_floatValue(line,positions,2_pInt)
|
|
case('peierlsstressedge','peierlsstress_edge')
|
|
do f = 1_pInt, Nchunks_SlipFamilies
|
|
peierlsStressPerSlipFamily(f,1_pInt,i) = IO_floatValue(line,positions,1_pInt+f)
|
|
enddo
|
|
case('peierlsstressscrew','peierlsstress_screw')
|
|
do f = 1_pInt, Nchunks_SlipFamilies
|
|
peierlsStressPerSlipFamily(f,2_pInt,i) = IO_floatValue(line,positions,1_pInt+f)
|
|
enddo
|
|
case('doublekinkwidth')
|
|
doublekinkwidth(i) = IO_floatValue(line,positions,2_pInt)
|
|
case('solidsolutionenergy')
|
|
solidSolutionEnergy(i) = IO_floatValue(line,positions,2_pInt)
|
|
case('solidsolutionsize')
|
|
solidSolutionSize(i) = IO_floatValue(line,positions,2_pInt)
|
|
case('solidsolutionconcentration')
|
|
solidSolutionConcentration(i) = IO_floatValue(line,positions,2_pInt)
|
|
case('p')
|
|
pParam(i) = IO_floatValue(line,positions,2_pInt)
|
|
case('q')
|
|
qParam(i) = IO_floatValue(line,positions,2_pInt)
|
|
case('viscosity','glideviscosity')
|
|
viscosity(i) = IO_floatValue(line,positions,2_pInt)
|
|
case('attackfrequency','fattack')
|
|
fattack(i) = IO_floatValue(line,positions,2_pInt)
|
|
case('rhosglscatter')
|
|
rhoSglScatter(i) = IO_floatValue(line,positions,2_pInt)
|
|
case('rhosglrandom')
|
|
rhoSglRandom(i) = IO_floatValue(line,positions,2_pInt)
|
|
case('rhosglrandombinning')
|
|
rhoSglRandomBinning(i) = IO_floatValue(line,positions,2_pInt)
|
|
case('surfacetransmissivity')
|
|
surfaceTransmissivity(i) = IO_floatValue(line,positions,2_pInt)
|
|
case('grainboundarytransmissivity')
|
|
grainboundaryTransmissivity(i) = IO_floatValue(line,positions,2_pInt)
|
|
case('cflfactor')
|
|
CFLfactor(i) = IO_floatValue(line,positions,2_pInt)
|
|
case('fedgemultiplication','edgemultiplicationfactor','edgemultiplication')
|
|
fEdgeMultiplication(i) = IO_floatValue(line,positions,2_pInt)
|
|
case('shortrangestresscorrection')
|
|
shortRangeStressCorrection(i) = IO_floatValue(line,positions,2_pInt) > 0.0_pReal
|
|
case ('nonschmid_coefficients')
|
|
if (positions(1) < 1_pInt + Nchunks_nonSchmid) &
|
|
call IO_warning(52_pInt,ext_msg=trim(tag)//' ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
do f = 1_pInt,Nchunks_nonSchmid
|
|
nonSchmidCoeff(f,i) = IO_floatValue(line,positions,1_pInt+f)
|
|
enddo
|
|
case('probabilisticmultiplication','randomsources','randommultiplication','discretesources')
|
|
probabilisticMultiplication(i) = IO_floatValue(line,positions,2_pInt) > 0.0_pReal
|
|
case default
|
|
call IO_error(210_pInt,ext_msg=trim(tag)//' ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
end select
|
|
endif
|
|
endif
|
|
enddo
|
|
|
|
|
|
do i = 1_pInt,maxNinstance
|
|
|
|
constitutive_nonlocal_structure(i) = &
|
|
lattice_initializeStructure(constitutive_nonlocal_structureName(i), CoverA(i)) ! our lattice structure is defined in the material.config file by the structureName (and the c/a ratio)
|
|
structID = constitutive_nonlocal_structure(i)
|
|
|
|
|
|
!*** sanity checks
|
|
|
|
if (structID < 1_pInt) &
|
|
call IO_error(205_pInt,el=i)
|
|
if (sum(Nslip(:,i)) <= 0_pInt) &
|
|
call IO_error(211_pInt,ext_msg='Nslip ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
do o = 1_pInt,maxval(phase_Noutput)
|
|
if(len(constitutive_nonlocal_output(o,i)) > 64_pInt) &
|
|
call IO_error(666_pInt)
|
|
enddo
|
|
do f = 1_pInt,lattice_maxNslipFamily
|
|
if (Nslip(f,i) > 0_pInt) then
|
|
if (rhoSglEdgePos0(f,i) < 0.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='rhoSglEdgePos0 ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (rhoSglEdgeNeg0(f,i) < 0.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='rhoSglEdgeNeg0 ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (rhoSglScrewPos0(f,i) < 0.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='rhoSglScrewPos0 ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (rhoSglScrewNeg0(f,i) < 0.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='rhoSglScrewNeg0 ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (rhoDipEdge0(f,i) < 0.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='rhoDipEdge0 ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (rhoDipScrew0(f,i) < 0.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='rhoDipScrew0 ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (burgersPerSlipFamily(f,i) <= 0.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='Burgers ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (lambda0PerSlipFamily(f,i) <= 0.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='lambda0 ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (minDipoleHeightPerSlipFamily(f,1,i) < 0.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='minimumDipoleHeightEdge ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (minDipoleHeightPerSlipFamily(f,2,i) < 0.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='minimumDipoleHeightScrew ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (peierlsStressPerSlipFamily(f,1,i) <= 0.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='peierlsStressEdge ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (peierlsStressPerSlipFamily(f,2,i) <= 0.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='peierlsStressScrew ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
endif
|
|
enddo
|
|
if (any(interactionSlipSlip(1:maxval(lattice_interactionSlipSlip(:,:,structID)),i) < 0.0_pReal)) &
|
|
call IO_error(211_pInt,ext_msg='interaction_SlipSlip ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (linetensionEffect(i) < 0.0_pReal .or. linetensionEffect(i) > 1.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='linetension ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (edgeJogFactor(i) < 0.0_pReal .or. edgeJogFactor(i) > 1.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='edgejog ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (cutoffRadius(i) < 0.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='r ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (atomicVolume(i) <= 0.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='atomicVolume ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (Dsd0(i) < 0.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='selfDiffusionPrefactor ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (selfDiffusionEnergy(i) <= 0.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='selfDiffusionEnergy ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (aTolRho(i) <= 0.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='aTol_rho ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (aTolShear(i) <= 0.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='aTol_shear ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (significantRho(i) < 0.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='significantRho ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (significantN(i) < 0.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='significantN ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (doublekinkwidth(i) <= 0.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='doublekinkwidth ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (solidSolutionEnergy(i) <= 0.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='solidSolutionEnergy ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (solidSolutionSize(i) <= 0.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='solidSolutionSize ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (solidSolutionConcentration(i) <= 0.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='solidSolutionConcentration ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (pParam(i) <= 0.0_pReal .or. pParam(i) > 1.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='p ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (qParam(i) < 1.0_pReal .or. qParam(i) > 2.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='q ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (viscosity(i) <= 0.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='viscosity ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (fattack(i) <= 0.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='attackFrequency ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (rhoSglScatter(i) < 0.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='rhoSglScatter ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (rhoSglRandom(i) < 0.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='rhoSglRandom ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (rhoSglRandomBinning(i) <= 0.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='rhoSglRandomBinning ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (surfaceTransmissivity(i) < 0.0_pReal .or. surfaceTransmissivity(i) > 1.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='surfaceTransmissivity ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (grainboundaryTransmissivity(i) > 1.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='grainboundaryTransmissivity ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (CFLfactor(i) < 0.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='CFLfactor ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
if (fEdgeMultiplication(i) < 0.0_pReal .or. fEdgeMultiplication(i) > 1.0_pReal) &
|
|
call IO_error(211_pInt,ext_msg='edgemultiplicationfactor ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
|
|
|
|
!*** determine total number of active slip systems
|
|
|
|
Nslip(1:lattice_maxNslipFamily,i) = min(lattice_NslipSystem(1:lattice_maxNslipFamily,structID), &
|
|
Nslip(1:lattice_maxNslipFamily,i) ) ! we can't use more slip systems per family than specified in lattice
|
|
totalNslip(i) = sum(Nslip(1:lattice_maxNslipFamily,i))
|
|
|
|
enddo
|
|
|
|
|
|
!*** allocation of variables whose size depends on the total number of active slip systems
|
|
|
|
maxTotalNslip = maxval(totalNslip)
|
|
|
|
allocate(iRhoU(maxTotalNslip,4,maxNinstance))
|
|
allocate(iRhoB(maxTotalNslip,4,maxNinstance))
|
|
allocate(iRhoD(maxTotalNslip,2,maxNinstance))
|
|
allocate(iV(maxTotalNslip,4,maxNinstance))
|
|
allocate(iD(maxTotalNslip,2,maxNinstance))
|
|
allocate(iGamma(maxTotalNslip,maxNinstance))
|
|
allocate(iRhoF(maxTotalNslip,maxNinstance))
|
|
allocate(iTauF(maxTotalNslip,maxNinstance))
|
|
allocate(iTauB(maxTotalNslip,maxNinstance))
|
|
iRhoU = 0_pInt
|
|
iRhoB = 0_pInt
|
|
iRhoD = 0_pInt
|
|
iV = 0_pInt
|
|
iD = 0_pInt
|
|
iGamma = 0_pInt
|
|
iRhoF = 0_pInt
|
|
iTauF = 0_pInt
|
|
iTauB = 0_pInt
|
|
|
|
allocate(burgers(maxTotalNslip, maxNinstance))
|
|
burgers = 0.0_pReal
|
|
|
|
allocate(lambda0(maxTotalNslip, maxNinstance))
|
|
lambda0 = 0.0_pReal
|
|
|
|
allocate(minDipoleHeight(maxTotalNslip,2,maxNinstance))
|
|
minDipoleHeight = -1.0_pReal
|
|
|
|
allocate(forestProjectionEdge(maxTotalNslip, maxTotalNslip, maxNinstance))
|
|
forestProjectionEdge = 0.0_pReal
|
|
|
|
allocate(forestProjectionScrew(maxTotalNslip, maxTotalNslip, maxNinstance))
|
|
forestProjectionScrew = 0.0_pReal
|
|
|
|
allocate(interactionMatrixSlipSlip(maxTotalNslip, maxTotalNslip, maxNinstance))
|
|
interactionMatrixSlipSlip = 0.0_pReal
|
|
|
|
allocate(lattice2slip(1:3, 1:3, maxTotalNslip, maxNinstance))
|
|
lattice2slip = 0.0_pReal
|
|
|
|
allocate(sourceProbability(maxTotalNslip, homogenization_maxNgrains, mesh_maxNips, mesh_NcpElems))
|
|
sourceProbability = 2.0_pReal
|
|
|
|
allocate(rhoDotFluxOutput(maxTotalNslip, 8, homogenization_maxNgrains, mesh_maxNips, mesh_NcpElems))
|
|
allocate(rhoDotMultiplicationOutput(maxTotalNslip, 2, homogenization_maxNgrains, mesh_maxNips, mesh_NcpElems))
|
|
allocate(rhoDotSingle2DipoleGlideOutput(maxTotalNslip, 2, homogenization_maxNgrains, mesh_maxNips, mesh_NcpElems))
|
|
allocate(rhoDotAthermalAnnihilationOutput(maxTotalNslip, 2, homogenization_maxNgrains, mesh_maxNips, mesh_NcpElems))
|
|
allocate(rhoDotThermalAnnihilationOutput(maxTotalNslip, 2, homogenization_maxNgrains, mesh_maxNips, mesh_NcpElems))
|
|
allocate(rhoDotEdgeJogsOutput(maxTotalNslip, homogenization_maxNgrains, mesh_maxNips, mesh_NcpElems))
|
|
rhoDotFluxOutput = 0.0_pReal
|
|
rhoDotMultiplicationOutput = 0.0_pReal
|
|
rhoDotSingle2DipoleGlideOutput = 0.0_pReal
|
|
rhoDotAthermalAnnihilationOutput = 0.0_pReal
|
|
rhoDotThermalAnnihilationOutput = 0.0_pReal
|
|
rhoDotEdgeJogsOutput = 0.0_pReal
|
|
|
|
allocate(compatibility(2,maxTotalNslip, maxTotalNslip, mesh_maxNipNeighbors, mesh_maxNips, mesh_NcpElems))
|
|
compatibility = 0.0_pReal
|
|
|
|
allocate(peierlsStress(maxTotalNslip,2,maxNinstance))
|
|
peierlsStress = 0.0_pReal
|
|
|
|
allocate(colinearSystem(maxTotalNslip,maxNinstance))
|
|
colinearSystem = 0_pInt
|
|
|
|
allocate(nonSchmidProjection(3,3,4,maxTotalNslip,maxNinstance))
|
|
nonSchmidProjection = 0.0_pReal
|
|
|
|
|
|
do i = 1,maxNinstance
|
|
|
|
structID = 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_pInt,lattice_maxNslipFamily
|
|
do s = 1_pInt,Nslip(f,i)
|
|
l = l + 1_pInt
|
|
slipFamily(l,i) = f
|
|
slipSystemLattice(l,i) = sum(lattice_NslipSystem(1:f-1_pInt, structID)) + s
|
|
enddo; enddo
|
|
|
|
|
|
!*** determine size of state array
|
|
|
|
ns = totalNslip(i)
|
|
constitutive_nonlocal_sizeDotState(i) = int(size(BASICSTATES),pInt) * ns
|
|
constitutive_nonlocal_sizeDependentState(i) = int(size(DEPENDENTSTATES),pInt) * ns
|
|
constitutive_nonlocal_sizeState(i) = constitutive_nonlocal_sizeDotState(i) &
|
|
+ constitutive_nonlocal_sizeDependentState(i) &
|
|
+ int(size(OTHERSTATES),pInt) * ns
|
|
|
|
!*** determine indices to state array
|
|
|
|
l = 0_pInt
|
|
do t = 1_pInt,4_pInt
|
|
do s = 1_pInt,ns
|
|
l = l + 1_pInt
|
|
iRhoU(s,t,i) = l
|
|
enddo
|
|
enddo
|
|
do t = 1_pInt,4_pInt
|
|
do s = 1_pInt,ns
|
|
l = l + 1_pInt
|
|
iRhoB(s,t,i) = l
|
|
enddo
|
|
enddo
|
|
do c = 1_pInt,2_pInt
|
|
do s = 1_pInt,ns
|
|
l = l + 1_pInt
|
|
iRhoD(s,c,i) = l
|
|
enddo
|
|
enddo
|
|
do s = 1_pInt,ns
|
|
l = l + 1_pInt
|
|
iGamma(s,i) = l
|
|
enddo
|
|
do s = 1_pInt,ns
|
|
l = l + 1_pInt
|
|
iRhoF(s,i) = l
|
|
enddo
|
|
do s = 1_pInt,ns
|
|
l = l + 1_pInt
|
|
iTauF(s,i) = l
|
|
enddo
|
|
do s = 1_pInt,ns
|
|
l = l + 1_pInt
|
|
iTauB(s,i) = l
|
|
enddo
|
|
do t = 1_pInt,4_pInt
|
|
do s = 1_pInt,ns
|
|
l = l + 1_pInt
|
|
iV(s,t,i) = l
|
|
enddo
|
|
enddo
|
|
do c = 1_pInt,2_pInt
|
|
do s = 1_pInt,ns
|
|
l = l + 1_pInt
|
|
iD(s,c,i) = l
|
|
enddo
|
|
enddo
|
|
if (iD(ns,2,i) /= constitutive_nonlocal_sizeState(i)) & ! check if last index is equal to size of state
|
|
call IO_error(0_pInt, ext_msg = 'state indices not properly set ('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
|
|
|
|
!*** determine size of postResults array
|
|
|
|
do o = 1_pInt,Noutput(i)
|
|
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_external', &
|
|
'resolvedstress_back', &
|
|
'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_sgl2dip_edge', &
|
|
'rho_dot_sgl2dip_screw', &
|
|
'rho_dot_ann_ath', &
|
|
'rho_dot_ann_the', &
|
|
'rho_dot_ann_the_edge', &
|
|
'rho_dot_ann_the_screw', &
|
|
'rho_dot_edgejogs', &
|
|
'rho_dot_flux', &
|
|
'rho_dot_flux_edge', &
|
|
'rho_dot_flux_screw', &
|
|
'velocity_edge_pos', &
|
|
'velocity_edge_neg', &
|
|
'velocity_screw_pos', &
|
|
'velocity_screw_neg', &
|
|
'slipdirection.x', &
|
|
'slipdirection.y', &
|
|
'slipdirection.z', &
|
|
'slipnormal.x', &
|
|
'slipnormal.y', &
|
|
'slipnormal.z', &
|
|
'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', &
|
|
'maximumdipoleheight_edge', &
|
|
'maximumdipoleheight_screw', &
|
|
'accumulatedshear', &
|
|
'boundarylayer' )
|
|
mySize = totalNslip(i)
|
|
case('dislocationstress')
|
|
mySize = 6_pInt
|
|
case default
|
|
call IO_error(212_pInt,ext_msg=constitutive_nonlocal_output(o,i)//&
|
|
'('//CONSTITUTIVE_NONLOCAL_LABEL//')')
|
|
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
|
|
|
|
Cslip66(:,:,i) = lattice_symmetrizeC66(constitutive_nonlocal_structureName(i), Cslip66(:,:,i))
|
|
mu(i) = 0.2_pReal * ( Cslip66(1,1,i) - Cslip66(1,2,i) + 3.0_pReal*Cslip66(4,4,i)) ! (C11iso-C12iso)/2 with C11iso=(3*C11+2*C12+4*C44)/5 and C12iso=(C11+4*C12-2*C44)/5
|
|
nu(i) = (Cslip66(1,1,i) + 4.0_pReal*Cslip66(1,2,i) - 2.0_pReal*Cslip66(4,4,i)) &
|
|
/ (4.0_pReal*Cslip66(1,1,i) + 6.0_pReal*Cslip66(1,2,i) + 2.0_pReal*Cslip66(4,4,i)) ! C12iso/(C11iso+C12iso) with C11iso=(3*C11+2*C12+4*C44)/5 and C12iso=(C11+4*C12-2*C44)/5
|
|
Cslip66(1:6,1:6,i) = math_Mandel3333to66(math_Voigt66to3333(Cslip66(1:6,1:6,i)))
|
|
Cslip3333(1:3,1:3,1:3,1:3,i) = math_Voigt66to3333(Cslip66(1:6,1:6,i))
|
|
|
|
do s1 = 1_pInt,ns
|
|
f = slipFamily(s1,i)
|
|
|
|
!*** burgers vector, mean free path prefactor and minimum dipole distance for each slip system
|
|
|
|
burgers(s1,i) = burgersPerSlipFamily(f,i)
|
|
lambda0(s1,i) = lambda0PerSlipFamily(f,i)
|
|
minDipoleHeight(s1,1:2,i) = minDipoleHeightPerSlipFamily(f,1:2,i)
|
|
peierlsStress(s1,1:2,i) = peierlsStressPerSlipFamily(f,1:2,i)
|
|
|
|
do s2 = 1_pInt,ns
|
|
|
|
!*** calculation of forest projections for edge and screw dislocations. s2 acts as forest for s1
|
|
|
|
forestProjectionEdge(s1,s2,i) &
|
|
= abs(math_mul3x3(lattice_sn(1:3,slipSystemLattice(s1,i),structID), &
|
|
lattice_st(1:3,slipSystemLattice(s2,i),structID))) ! forest projection of edge dislocations is the projection of (t = b x n) onto the slip normal of the respective slip plane
|
|
|
|
forestProjectionScrew(s1,s2,i) &
|
|
= abs(math_mul3x3(lattice_sn(1:3,slipSystemLattice(s1,i),structID), &
|
|
lattice_sd(1:3,slipSystemLattice(s2,i),structID))) ! forest projection of screw dislocations is the projection of b onto the slip normal of the respective splip plane
|
|
|
|
!*** calculation of interaction matrices
|
|
|
|
interactionMatrixSlipSlip(s1,s2,i) &
|
|
= interactionSlipSlip(lattice_interactionSlipSlip(slipSystemLattice(s1,i), &
|
|
slipSystemLattice(s2,i), &
|
|
structID), i)
|
|
|
|
!*** colinear slip system (only makes sense for fcc like it is defined here)
|
|
|
|
if (lattice_interactionSlipSlip(slipSystemLattice(s1,i), &
|
|
slipSystemLattice(s2,i), &
|
|
structID) == 3_pInt) then
|
|
colinearSystem(s1,i) = s2
|
|
endif
|
|
|
|
enddo
|
|
|
|
!*** rotation matrix from lattice configuration to slip system
|
|
|
|
lattice2slip(1:3,1:3,s1,i) &
|
|
= math_transpose33( reshape([ lattice_sd(1:3, slipSystemLattice(s1,i), structID), &
|
|
-lattice_st(1:3, slipSystemLattice(s1,i), structID), &
|
|
lattice_sn(1:3, slipSystemLattice(s1,i), structID)], [3,3]))
|
|
enddo
|
|
|
|
|
|
!*** combined projection of Schmid and non-Schmid contributions to the resolved shear stress (only for screws)
|
|
!* four types t:
|
|
!* 1) positive screw at positive resolved stress
|
|
!* 2) positive screw at negative resolved stress
|
|
!* 3) negative screw at positive resolved stress
|
|
!* 4) negative screw at negative resolved stress
|
|
|
|
do s = 1_pInt,ns
|
|
do l = 1_pInt,lattice_NnonSchmid(structID)
|
|
nonSchmidProjection(1:3,1:3,1,s,i) = nonSchmidProjection(1:3,1:3,1,s,i) &
|
|
+ nonSchmidCoeff(l,i) * lattice_Sslip(1:3,1:3,2*l,slipSystemLattice(s,i),structID)
|
|
nonSchmidProjection(1:3,1:3,2,s,i) = nonSchmidProjection(1:3,1:3,2,s,i) &
|
|
+ nonSchmidCoeff(l,i) * lattice_Sslip(1:3,1:3,2*l+1,slipSystemLattice(s,i),structID)
|
|
enddo
|
|
nonSchmidProjection(1:3,1:3,3,s,i) = -nonSchmidProjection(1:3,1:3,2,s,i)
|
|
nonSchmidProjection(1:3,1:3,4,s,i) = -nonSchmidProjection(1:3,1:3,1,s,i)
|
|
forall (t = 1:4) &
|
|
nonSchmidProjection(1:3,1:3,t,s,i) = nonSchmidProjection(1:3,1:3,t,s,i) &
|
|
+ lattice_Sslip(1:3,1:3,1,slipSystemLattice(s,i),structID)
|
|
enddo
|
|
|
|
enddo
|
|
|
|
end subroutine constitutive_nonlocal_init
|
|
|
|
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
!> @brief sets the initial microstructural state for a given instance of this plasticity
|
|
!--------------------------------------------------------------------------------------------------
|
|
subroutine constitutive_nonlocal_stateInit(state)
|
|
|
|
use IO, only: IO_error
|
|
use lattice, only: lattice_maxNslipFamily
|
|
use math, only: math_sampleGaussVar
|
|
use mesh, only: mesh_ipVolume, &
|
|
mesh_NcpElems, &
|
|
mesh_maxNips, &
|
|
mesh_element, &
|
|
FE_Nips, &
|
|
FE_geomtype
|
|
use material, only: material_phase, &
|
|
phase_plasticityInstance, &
|
|
phase_plasticity, &
|
|
homogenization_Ngrains
|
|
|
|
implicit none
|
|
|
|
!*** input/output variables
|
|
type(p_vec), dimension(1,mesh_maxNips,mesh_NcpElems), intent(inout) :: &
|
|
state ! microstructural state
|
|
|
|
!*** local variables
|
|
integer(pInt) el, &
|
|
ip, &
|
|
e, &
|
|
i, &
|
|
idx, &
|
|
ns, & ! short notation for total number of active slip systems
|
|
f, & ! index of lattice family
|
|
from, &
|
|
upto, &
|
|
s, & ! index of slip system
|
|
t, &
|
|
j, &
|
|
matID, &
|
|
maxNinstance
|
|
real(pReal), dimension(2) :: noise
|
|
real(pReal), dimension(4) :: rnd
|
|
real(pReal) meanDensity, &
|
|
totalVolume, &
|
|
densityBinning, &
|
|
minimumIpVolume
|
|
|
|
|
|
maxNinstance = int(count(phase_plasticity == CONSTITUTIVE_NONLOCAL_LABEL),pInt)
|
|
|
|
|
|
! ititalize all states to zero
|
|
|
|
do e = 1_pInt,mesh_NcpElems
|
|
do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e)))
|
|
if (CONSTITUTIVE_NONLOCAL_LABEL == phase_plasticity(material_phase(1,i,e))) &
|
|
state(1,i,e)%p = 0.0_pReal
|
|
enddo
|
|
enddo
|
|
|
|
|
|
do matID = 1_pInt,maxNinstance
|
|
ns = totalNslip(matID)
|
|
|
|
! randomly distribute dislocation segments on random slip system and of random type in the volume
|
|
if (rhoSglRandom(matID) > 0.0_pReal) then
|
|
|
|
! get the total volume of the instance
|
|
|
|
minimumIpVolume = 1e99_pReal
|
|
totalVolume = 0.0_pReal
|
|
do e = 1_pInt,mesh_NcpElems
|
|
do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e)))
|
|
if (CONSTITUTIVE_NONLOCAL_LABEL == phase_plasticity(material_phase(1,i,e)) &
|
|
.and. matID == phase_plasticityInstance(material_phase(1,i,e))) then
|
|
totalVolume = totalVolume + mesh_ipVolume(i,e)
|
|
minimumIpVolume = min(minimumIpVolume, mesh_ipVolume(i,e))
|
|
endif
|
|
enddo
|
|
enddo
|
|
densityBinning = rhoSglRandomBinning(matID) / minimumIpVolume ** (2.0_pReal / 3.0_pReal)
|
|
|
|
! subsequently fill random ips with dislocation segments until we reach the desired overall density
|
|
|
|
meanDensity = 0.0_pReal
|
|
do while(meanDensity < rhoSglRandom(matID))
|
|
call random_number(rnd)
|
|
el = nint(rnd(1)*real(mesh_NcpElems,pReal)+0.5_pReal,pInt)
|
|
ip = nint(rnd(2)*real(FE_Nips(FE_geomtype(mesh_element(2,el))),pReal)+0.5_pReal,pInt)
|
|
if (CONSTITUTIVE_NONLOCAL_LABEL == phase_plasticity(material_phase(1,ip,el)) &
|
|
.and. matID == phase_plasticityInstance(material_phase(1,ip,el))) then
|
|
s = nint(rnd(3)*real(ns,pReal)+0.5_pReal,pInt)
|
|
t = nint(rnd(4)*4.0_pReal+0.5_pReal,pInt)
|
|
meanDensity = meanDensity + densityBinning * mesh_ipVolume(ip,el) / totalVolume
|
|
state(1,ip,el)%p(iRhoU(s,t,matID)) = state(1,ip,el)%p(iRhoU(s,t,matID)) + densityBinning
|
|
endif
|
|
enddo
|
|
|
|
! homogeneous distribution of density with some noise
|
|
else
|
|
do e = 1_pInt,mesh_NcpElems
|
|
do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e)))
|
|
if (CONSTITUTIVE_NONLOCAL_LABEL == phase_plasticity(material_phase(1,i,e)) &
|
|
.and. matID == phase_plasticityInstance(material_phase(1,i,e))) then
|
|
do f = 1_pInt,lattice_maxNslipFamily
|
|
from = 1_pInt + sum(Nslip(1:f-1_pInt,matID))
|
|
upto = sum(Nslip(1:f,matID))
|
|
do s = from,upto
|
|
do j = 1_pInt,2_pInt
|
|
noise(j) = math_sampleGaussVar(0.0_pReal, rhoSglScatter(matID))
|
|
enddo
|
|
state(1,i,e)%p(iRhoU(s,1,matID)) = rhoSglEdgePos0(f, matID) + noise(1)
|
|
state(1,i,e)%p(iRhoU(s,2,matID)) = rhoSglEdgeNeg0(f, matID) + noise(1)
|
|
state(1,i,e)%p(iRhoU(s,3,matID)) = rhoSglScrewPos0(f, matID) + noise(2)
|
|
state(1,i,e)%p(iRhoU(s,4,matID)) = rhoSglScrewNeg0(f, matID) + noise(2)
|
|
enddo
|
|
state(1,i,e)%p(iRhoD(from:upto,1,matID)) = rhoDipEdge0(f, matID)
|
|
state(1,i,e)%p(iRhoD(from:upto,2,matID)) = rhoDipScrew0(f, matID)
|
|
enddo
|
|
endif
|
|
enddo
|
|
enddo
|
|
endif
|
|
enddo
|
|
|
|
end subroutine constitutive_nonlocal_stateInit
|
|
|
|
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
!> @brief sets the relevant state values for a given instance of this plasticity
|
|
!--------------------------------------------------------------------------------------------------
|
|
pure function constitutive_nonlocal_aTolState(matID)
|
|
|
|
implicit none
|
|
|
|
!*** input variables
|
|
integer(pInt), intent(in) :: matID ! number specifying the current instance of the plasticity
|
|
|
|
!*** output variables
|
|
real(pReal), dimension(constitutive_nonlocal_sizeState(matID)) :: &
|
|
constitutive_nonlocal_aTolState ! absolute state tolerance for the current instance of this plasticity
|
|
|
|
!*** local variables
|
|
integer(pInt) :: ns, t, c
|
|
|
|
ns = totalNslip(matID)
|
|
constitutive_nonlocal_aTolState = 0.0_pReal
|
|
forall (t = 1_pInt:4_pInt)
|
|
constitutive_nonlocal_aTolState(iRhoU(1:ns,t,matID)) = aTolRho(matID)
|
|
constitutive_nonlocal_aTolState(iRhoB(1:ns,t,matID)) = aTolRho(matID)
|
|
endforall
|
|
forall (c = 1_pInt:2_pInt) &
|
|
constitutive_nonlocal_aTolState(iRhoD(1:ns,c,matID)) = aTolRho(matID)
|
|
constitutive_nonlocal_aTolState(iGamma(1:ns,matID)) = aTolShear(matID)
|
|
|
|
end function constitutive_nonlocal_aTolState
|
|
|
|
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
!> @brief returns the homogenized elasticity matrix
|
|
!--------------------------------------------------------------------------------------------------
|
|
pure function constitutive_nonlocal_homogenizedC(state,g,ip,el)
|
|
|
|
use mesh, only: mesh_NcpElems, &
|
|
mesh_maxNips
|
|
use material, only: homogenization_maxNgrains, &
|
|
material_phase, &
|
|
phase_plasticityInstance
|
|
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) matID ! current instance of this plasticity
|
|
|
|
matID = phase_plasticityInstance(material_phase(g,ip,el))
|
|
|
|
constitutive_nonlocal_homogenizedC = Cslip66(1:6,1:6,matID)
|
|
|
|
end function constitutive_nonlocal_homogenizedC
|
|
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
!> @brief calculates quantities characterizing the microstructure
|
|
!--------------------------------------------------------------------------------------------------
|
|
subroutine constitutive_nonlocal_microstructure(state, Temperature, Fe, Fp, gr, ip, el)
|
|
|
|
use IO, only: &
|
|
IO_error
|
|
use math, only: &
|
|
pi, &
|
|
math_mul33x3, &
|
|
math_mul3x3, &
|
|
math_norm3, &
|
|
math_invert33, &
|
|
math_transpose33
|
|
use debug, only: &
|
|
debug_level, &
|
|
debug_constitutive, &
|
|
debug_levelBasic, &
|
|
debug_levelExtensive, &
|
|
debug_levelSelective, &
|
|
debug_g, &
|
|
debug_i, &
|
|
debug_e
|
|
use mesh, only: &
|
|
mesh_NcpElems, &
|
|
mesh_maxNips, &
|
|
mesh_element, &
|
|
mesh_ipNeighborhood, &
|
|
mesh_ipCoordinates, &
|
|
mesh_ipVolume, &
|
|
mesh_ipAreaNormal, &
|
|
mesh_ipArea, &
|
|
FE_NipNeighbors, &
|
|
mesh_maxNipNeighbors, &
|
|
FE_geomtype, &
|
|
FE_celltype
|
|
use material, only: &
|
|
homogenization_maxNgrains, &
|
|
material_phase, &
|
|
phase_localPlasticity, &
|
|
phase_plasticityInstance
|
|
use lattice, only: &
|
|
lattice_sd, &
|
|
lattice_st, &
|
|
lattice_interactionSlipSlip
|
|
|
|
implicit none
|
|
|
|
!*** input variables
|
|
integer(pInt), intent(in) :: gr, & ! current grain ID
|
|
ip, & ! current integration point
|
|
el ! current element
|
|
real(pReal), intent(in) :: Temperature ! temperature
|
|
real(pReal), dimension(3,3), intent(in) :: &
|
|
Fe, & ! elastic deformation gradient
|
|
Fp ! 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 plasticity
|
|
neighboring_instance, & ! instance of this plasticity 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, &
|
|
n, &
|
|
nRealNeighbors ! number of really existing neighbors
|
|
integer(pInt), dimension(2) :: neighbor
|
|
real(pReal) detFe, &
|
|
detFp, &
|
|
FVsize, &
|
|
temp, &
|
|
correction, &
|
|
myRhoForest
|
|
real(pReal), dimension(2) :: rhoExcessGradient, &
|
|
rhoExcessGradient_over_rho, &
|
|
rhoTotal
|
|
real(pReal), dimension(3) :: rhoExcessDifferences, &
|
|
normal_latticeConf
|
|
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(gr,ip,el)))) :: &
|
|
rhoForest, & ! forest dislocation density
|
|
tauBack, & ! back stress from pileup on same slip system
|
|
tauThreshold ! threshold shear stress
|
|
real(pReal), dimension(3,3) :: invFe, & ! inverse of elastic deformation gradient
|
|
invFp, & ! inverse of plastic deformation gradient
|
|
connections, &
|
|
invConnections
|
|
real(pReal), dimension(3,mesh_maxNipNeighbors) :: &
|
|
connection_latticeConf
|
|
real(pReal), dimension(2,totalNslip(phase_plasticityInstance(material_phase(gr,ip,el)))) :: &
|
|
rhoExcess
|
|
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(gr,ip,el))),2) :: &
|
|
rhoDip ! dipole dislocation density (edge, screw)
|
|
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(gr,ip,el))),8) :: &
|
|
rhoSgl ! single dislocation density (edge+, edge-, screw+, screw-, used edge+, used edge-, used screw+, used screw-)
|
|
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(gr,ip,el))), &
|
|
totalNslip(phase_plasticityInstance(material_phase(gr,ip,el)))) :: &
|
|
myInteractionMatrix ! corrected slip interaction matrix
|
|
real(pReal), dimension(2,maxval(totalNslip),mesh_maxNipNeighbors) :: &
|
|
neighboring_rhoExcess, & ! excess density at neighboring material point
|
|
neighboring_rhoTotal ! total density at neighboring material point
|
|
real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(gr,ip,el))),2) :: &
|
|
m ! direction of dislocation motion
|
|
logical inversionError
|
|
|
|
|
|
phase = material_phase(gr,ip,el)
|
|
instance = phase_plasticityInstance(phase)
|
|
latticeStruct = constitutive_nonlocal_structure(instance)
|
|
ns = totalNslip(instance)
|
|
|
|
|
|
!*** get basic states
|
|
|
|
forall (s = 1_pInt:ns, t = 1_pInt:4_pInt)
|
|
rhoSgl(s,t) = max(state(gr,ip,el)%p(iRhoU(s,t,instance)), 0.0_pReal) ! ensure positive single mobile densities
|
|
rhoSgl(s,t+4_pInt) = state(gr,ip,el)%p(iRhoB(s,t,instance))
|
|
endforall
|
|
forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) &
|
|
rhoDip(s,c) = max(state(gr,ip,el)%p(iRhoD(s,c,instance)), 0.0_pReal) ! ensure positive dipole densities
|
|
where (abs(rhoSgl) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance) &
|
|
.or. abs(rhoSgl) < significantRho(instance)) &
|
|
rhoSgl = 0.0_pReal
|
|
where (abs(rhoDip) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance) &
|
|
.or. abs(rhoDip) < significantRho(instance)) &
|
|
rhoDip = 0.0_pReal
|
|
|
|
|
|
!*** calculate the forest dislocation density
|
|
!*** (= projection of screw and edge dislocations)
|
|
|
|
forall (s = 1_pInt:ns) &
|
|
rhoForest(s) = dot_product((sum(abs(rhoSgl(1:ns,[1,2,5,6])),2) + rhoDip(1:ns,1)), &
|
|
forestProjectionEdge(s,1:ns,instance)) &
|
|
+ dot_product((sum(abs(rhoSgl(1:ns,[3,4,7,8])),2) + rhoDip(1:ns,2)), &
|
|
forestProjectionScrew(s,1:ns,instance))
|
|
|
|
|
|
|
|
!*** calculate the threshold shear stress for dislocation slip
|
|
!*** coefficients are corrected for the line tension effect
|
|
!*** (see Kubin,Devincre,Hoc; 2008; Modeling dislocation storage rates and mean free paths in face-centered cubic crystals)
|
|
|
|
myInteractionMatrix = 0.0_pReal
|
|
myInteractionMatrix(1:ns,1:ns) = interactionMatrixSlipSlip(1:ns,1:ns,instance)
|
|
if (latticeStruct < 3_pInt) then ! only fcc and bcc
|
|
do s = 1_pInt,ns
|
|
myRhoForest = max(rhoForest(s),significantRho(instance))
|
|
correction = ( 1.0_pReal - linetensionEffect(instance) &
|
|
+ linetensionEffect(instance) &
|
|
* log(0.35_pReal * burgers(s,instance) * sqrt(myRhoForest)) &
|
|
/ log(0.35_pReal * burgers(s,instance) * 1e6_pReal)) ** 2.0_pReal
|
|
myInteractionMatrix(s,1:ns) = correction * myInteractionMatrix(s,1:ns)
|
|
enddo
|
|
endif
|
|
forall (s = 1_pInt:ns) &
|
|
tauThreshold(s) = mu(instance) * burgers(s,instance) &
|
|
* sqrt(dot_product((sum(abs(rhoSgl),2) + sum(abs(rhoDip),2)), myInteractionMatrix(s,1:ns)))
|
|
|
|
|
|
|
|
!*** calculate the dislocation stress of the neighboring excess dislocation densities
|
|
!*** zero for material points of local plasticity
|
|
|
|
tauBack = 0.0_pReal
|
|
|
|
if (.not. phase_localPlasticity(phase) .and. shortRangeStressCorrection(instance)) then
|
|
call math_invert33(Fe, invFe, detFe, inversionError)
|
|
call math_invert33(Fp, invFp, detFp, inversionError)
|
|
rhoExcess(1,1:ns) = rhoSgl(1:ns,1) - rhoSgl(1:ns,2)
|
|
rhoExcess(2,1:ns) = rhoSgl(1:ns,3) - rhoSgl(1:ns,4)
|
|
FVsize = mesh_ipVolume(ip,el) ** (1.0_pReal/3.0_pReal)
|
|
|
|
!* loop through my neighborhood and get the connection vectors (in lattice frame) and the excess densities
|
|
|
|
nRealNeighbors = 0_pInt
|
|
neighboring_rhoTotal = 0.0_pReal
|
|
do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el))))
|
|
neighboring_el = mesh_ipNeighborhood(1,n,ip,el)
|
|
neighboring_ip = mesh_ipNeighborhood(2,n,ip,el)
|
|
if (neighboring_el > 0 .and. neighboring_ip > 0) then
|
|
neighboring_phase = material_phase(gr,neighboring_ip,neighboring_el)
|
|
neighboring_instance = phase_plasticityInstance(neighboring_phase)
|
|
neighboring_latticeStruct = constitutive_nonlocal_structure(neighboring_instance)
|
|
neighboring_ns = totalNslip(neighboring_instance)
|
|
if (.not. phase_localPlasticity(neighboring_phase) &
|
|
.and. neighboring_latticeStruct == latticeStruct &
|
|
.and. neighboring_instance == instance) then
|
|
if (neighboring_ns == ns) then
|
|
nRealNeighbors = nRealNeighbors + 1_pInt
|
|
forall (s = 1_pInt:ns, c = 1_pInt:2_pInt)
|
|
neighboring_rhoExcess(c,s,n) = &
|
|
max(state(gr,neighboring_ip,neighboring_el)%p(iRhoU(s,2*c-1,neighboring_instance)), 0.0_pReal) &! positive mobiles
|
|
- max(state(gr,neighboring_ip,neighboring_el)%p(iRhoU(s,2*c,neighboring_instance)), 0.0_pReal) ! negative mobiles
|
|
neighboring_rhoTotal(c,s,n) = &
|
|
max(state(gr,neighboring_ip,neighboring_el)%p(iRhoU(s,2*c-1,neighboring_instance)), 0.0_pReal) &! positive mobiles
|
|
+ max(state(gr,neighboring_ip,neighboring_el)%p(iRhoU(s,2*c,neighboring_instance)), 0.0_pReal) & ! negative mobiles
|
|
+ abs(state(gr,neighboring_ip,neighboring_el)%p(iRhoB(s,2*c-1,neighboring_instance))) & ! positive deads
|
|
+ abs(state(gr,neighboring_ip,neighboring_el)%p(iRhoB(s,2*c,neighboring_instance))) & ! negative deads
|
|
+ max(state(gr,neighboring_ip,neighboring_el)%p(iRhoD(s,c,neighboring_instance)), 0.0_pReal) ! dipoles
|
|
endforall
|
|
connection_latticeConf(1:3,n) = &
|
|
math_mul33x3(invFe, mesh_ipCoordinates(1:3,neighboring_ip,neighboring_el) &
|
|
- mesh_ipCoordinates(1:3,ip,el))
|
|
normal_latticeConf = math_mul33x3(math_transpose33(invFp), mesh_ipAreaNormal(1:3,n,ip,el))
|
|
if (math_mul3x3(normal_latticeConf,connection_latticeConf(1:3,n)) < 0.0_pReal) then ! neighbor connection points in opposite direction to face normal: must be periodic image
|
|
connection_latticeConf(1:3,n) = normal_latticeConf * mesh_ipVolume(ip,el) &
|
|
/ mesh_ipArea(n,ip,el) ! instead take the surface normal scaled with the diameter of the cell
|
|
endif
|
|
else
|
|
! different number of active slip systems
|
|
call IO_error(-1_pInt,ext_msg='different number of active slip systems in neighboring IPs of same crystal structure')
|
|
endif
|
|
else
|
|
! local neighbor or different lattice structure or different constitution instance -> use central values instead
|
|
connection_latticeConf(1:3,n) = 0.0_pReal
|
|
neighboring_rhoExcess(1:2,1:ns,n) = rhoExcess
|
|
endif
|
|
else
|
|
! free surface -> use central values instead
|
|
connection_latticeConf(1:3,n) = 0.0_pReal
|
|
neighboring_rhoExcess(1:2,1:ns,n) = rhoExcess
|
|
endif
|
|
enddo
|
|
|
|
|
|
!* loop through the slip systems and calculate the dislocation gradient by
|
|
!* 1. interpolation of the excess density in the neighorhood
|
|
!* 2. interpolation of the dead dislocation density in the central volume
|
|
|
|
m(1:3,1:ns,1) = lattice_sd(1:3,slipSystemLattice(1:ns,instance),latticeStruct)
|
|
m(1:3,1:ns,2) = -lattice_st(1:3,slipSystemLattice(1:ns,instance),latticeStruct)
|
|
|
|
do s = 1_pInt,ns
|
|
|
|
!* gradient from interpolation of neighboring excess density
|
|
|
|
do c = 1_pInt,2_pInt
|
|
do dir = 1_pInt,3_pInt
|
|
neighbor(1) = 2_pInt * dir - 1_pInt
|
|
neighbor(2) = 2_pInt * dir
|
|
connections(dir,1:3) = connection_latticeConf(1:3,neighbor(1)) &
|
|
- connection_latticeConf(1:3,neighbor(2))
|
|
rhoExcessDifferences(dir) = neighboring_rhoExcess(c,s,neighbor(1)) &
|
|
- neighboring_rhoExcess(c,s,neighbor(2))
|
|
enddo
|
|
call math_invert33(connections,invConnections,temp,inversionError)
|
|
if (inversionError) then
|
|
call IO_error(-1_pInt,ext_msg='back stress calculation: inversion error')
|
|
endif
|
|
rhoExcessGradient(c) = math_mul3x3(m(1:3,s,c), &
|
|
math_mul33x3(invConnections,rhoExcessDifferences))
|
|
enddo
|
|
|
|
!* plus gradient from deads
|
|
|
|
do t = 1_pInt,4_pInt
|
|
c = (t - 1_pInt) / 2_pInt + 1_pInt
|
|
rhoExcessGradient(c) = rhoExcessGradient(c) + rhoSgl(s,t+4_pInt) / FVsize
|
|
enddo
|
|
|
|
!* normalized with the total density
|
|
|
|
rhoExcessGradient_over_rho = 0.0_pReal
|
|
forall (c = 1_pInt:2_pInt) &
|
|
rhoTotal(c) = (sum(abs(rhoSgl(s,[2*c-1,2*c,2*c+3,2*c+4]))) + rhoDip(s,c) + sum(neighboring_rhoTotal(c,s,:))) &
|
|
/ real(1_pInt + nRealNeighbors,pReal)
|
|
forall (c = 1_pInt:2_pInt, rhoTotal(c) > 0.0_pReal) &
|
|
rhoExcessGradient_over_rho(c) = rhoExcessGradient(c) / rhoTotal(c)
|
|
|
|
!* gives the local stress correction when multiplied with a factor
|
|
|
|
tauBack(s) = - mu(instance) * burgers(s,instance) / (2.0_pReal * pi) &
|
|
* (rhoExcessGradient_over_rho(1) / (1.0_pReal - nu(instance)) + rhoExcessGradient_over_rho(2))
|
|
|
|
enddo
|
|
endif
|
|
|
|
|
|
!*** set dependent states
|
|
|
|
state(gr,ip,el)%p(iRhoF(1:ns,instance)) = rhoForest
|
|
state(gr,ip,el)%p(iTauF(1:ns,instance)) = tauThreshold
|
|
state(gr,ip,el)%p(iTauB(1:ns,instance)) = tauBack
|
|
|
|
|
|
#ifndef _OPENMP
|
|
if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt &
|
|
.and. ((debug_e == el .and. debug_i == ip .and. debug_g == gr)&
|
|
.or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then
|
|
write(6,*)
|
|
write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_microstructure at el ip g',el,ip,gr
|
|
write(6,*)
|
|
write(6,'(a,/,12x,12(e10.3,1x))') '<< CONST >> rhoForest', rhoForest
|
|
write(6,'(a,/,12x,12(f10.5,1x))') '<< CONST >> tauThreshold / MPa', tauThreshold/1e6
|
|
write(6,'(a,/,12x,12(f10.5,1x))') '<< CONST >> tauBack / MPa', tauBack/1e6
|
|
write(6,*)
|
|
endif
|
|
#endif
|
|
|
|
end subroutine constitutive_nonlocal_microstructure
|
|
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
!> @brief calculates kinetics
|
|
!--------------------------------------------------------------------------------------------------
|
|
subroutine constitutive_nonlocal_kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, &
|
|
tauThreshold, c, Temperature, g, ip, el)
|
|
|
|
use debug, only: debug_level, &
|
|
debug_constitutive, &
|
|
debug_levelBasic, &
|
|
debug_levelExtensive, &
|
|
debug_levelSelective, &
|
|
debug_g, &
|
|
debug_i, &
|
|
debug_e
|
|
use material, only: material_phase, &
|
|
phase_plasticityInstance
|
|
|
|
implicit none
|
|
|
|
!*** input variables
|
|
integer(pInt), intent(in) :: g, & !< current grain number
|
|
ip, & !< current integration point
|
|
el, & !< current element number
|
|
c !< dislocation character (1:edge, 2:screw)
|
|
real(pReal), intent(in) :: Temperature !< temperature
|
|
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el)))), &
|
|
intent(in) :: tau, & !< resolved external shear stress (without non Schmid effects)
|
|
tauNS, & !< resolved external shear stress (including non Schmid effects)
|
|
tauThreshold !< threshold shear stress
|
|
|
|
!*** output variables
|
|
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el)))), &
|
|
intent(out) :: v, & !< velocity
|
|
dv_dtau, & !< velocity derivative with respect to resolved shear stress (without non Schmid contributions)
|
|
dv_dtauNS !< velocity derivative with respect to resolved shear stress (including non Schmid contributions)
|
|
|
|
!*** local variables
|
|
integer(pInt) :: instance, & !< current instance of this plasticity
|
|
ns, & !< short notation for the total number of active slip systems
|
|
s !< index of my current slip system
|
|
real(pReal) tauRel_P, &
|
|
tauRel_S, &
|
|
tauEff, & !< effective shear stress
|
|
tPeierls, & !< waiting time in front of a peierls barriers
|
|
tSolidSolution, & !< waiting time in front of a solid solution obstacle
|
|
vViscous, & !< viscous glide velocity
|
|
dtPeierls_dtau, & !< derivative with respect to resolved shear stress
|
|
dtSolidSolution_dtau, & !< derivative with respect to resolved shear stress
|
|
meanfreepath_S, & !< mean free travel distance for dislocations between two solid solution obstacles
|
|
meanfreepath_P, & !< mean free travel distance for dislocations between two Peierls barriers
|
|
jumpWidth_P, & !< depth of activated area
|
|
jumpWidth_S, & !< depth of activated area
|
|
activationLength_P, & !< length of activated dislocation line
|
|
activationLength_S, & !< length of activated dislocation line
|
|
activationVolume_P, & !< volume that needs to be activated to overcome barrier
|
|
activationVolume_S, & !< volume that needs to be activated to overcome barrier
|
|
activationEnergy_P, & !< energy that is needed to overcome barrier
|
|
activationEnergy_S, & !< energy that is needed to overcome barrier
|
|
criticalStress_P, & !< maximum obstacle strength
|
|
criticalStress_S, & !< maximum obstacle strength
|
|
mobility !< dislocation mobility
|
|
|
|
|
|
instance = phase_plasticityInstance(material_phase(g,ip,el))
|
|
ns = totalNslip(instance)
|
|
|
|
v = 0.0_pReal
|
|
dv_dtau = 0.0_pReal
|
|
dv_dtauNS = 0.0_pReal
|
|
|
|
|
|
if (Temperature > 0.0_pReal) then
|
|
do s = 1_pInt,ns
|
|
if (abs(tau(s)) > tauThreshold(s)) then
|
|
|
|
!* Peierls contribution
|
|
!* Effective stress includes non Schmid constributions
|
|
!* The derivative only gives absolute values; the correct sign is taken care of in the formula for the derivative of the velocity
|
|
|
|
tauEff = max(0.0_pReal, abs(tauNS(s)) - tauThreshold(s)) ! ensure that the effective stress is positive
|
|
meanfreepath_P = burgers(s,instance)
|
|
jumpWidth_P = burgers(s,instance)
|
|
activationLength_P = doublekinkwidth(instance) * burgers(s,instance)
|
|
activationVolume_P = activationLength_P * jumpWidth_P * burgers(s,instance)
|
|
criticalStress_P = peierlsStress(s,c,instance)
|
|
activationEnergy_P = criticalStress_P * activationVolume_P
|
|
tauRel_P = min(1.0_pReal, tauEff / criticalStress_P) ! ensure that the activation probability cannot become greater than one
|
|
tPeierls = 1.0_pReal / fattack(instance) &
|
|
* exp(activationEnergy_P / (KB * Temperature) &
|
|
* (1.0_pReal - tauRel_P**pParam(instance))**qParam(instance))
|
|
if (tauEff < criticalStress_P) then
|
|
dtPeierls_dtau = tPeierls * pParam(instance) * qParam(instance) * activationVolume_P / (KB * Temperature) &
|
|
* (1.0_pReal - tauRel_P**pParam(instance))**(qParam(instance)-1.0_pReal) &
|
|
* tauRel_P**(pParam(instance)-1.0_pReal)
|
|
else
|
|
dtPeierls_dtau = 0.0_pReal
|
|
endif
|
|
|
|
|
|
!* Contribution from solid solution strengthening
|
|
!* The derivative only gives absolute values; the correct sign is taken care of in the formula for the derivative of the velocity
|
|
|
|
tauEff = abs(tau(s)) - tauThreshold(s)
|
|
meanfreepath_S = burgers(s,instance) / sqrt(solidSolutionConcentration(instance))
|
|
jumpWidth_S = solidSolutionSize(instance) * burgers(s,instance)
|
|
activationLength_S = burgers(s,instance) / sqrt(solidSolutionConcentration(instance))
|
|
activationVolume_S = activationLength_S * jumpWidth_S * burgers(s,instance)
|
|
activationEnergy_S = solidSolutionEnergy(instance)
|
|
criticalStress_S = activationEnergy_S / activationVolume_S
|
|
tauRel_S = min(1.0_pReal, tauEff / criticalStress_S) ! ensure that the activation probability cannot become greater than one
|
|
tSolidSolution = 1.0_pReal / fattack(instance) &
|
|
* exp(activationEnergy_S / (KB * Temperature) &
|
|
* (1.0_pReal - tauRel_S**pParam(instance))**qParam(instance))
|
|
if (tauEff < criticalStress_S) then
|
|
dtSolidSolution_dtau = tSolidSolution * pParam(instance) * qParam(instance) &
|
|
* activationVolume_S / (KB * Temperature) &
|
|
* (1.0_pReal - tauRel_S**pParam(instance))**(qParam(instance)-1.0_pReal) &
|
|
* tauRel_S**(pParam(instance)-1.0_pReal)
|
|
else
|
|
dtSolidSolution_dtau = 0.0_pReal
|
|
endif
|
|
|
|
|
|
!* viscous glide velocity
|
|
|
|
tauEff = abs(tau(s)) - tauThreshold(s)
|
|
mobility = burgers(s,instance) / viscosity(instance)
|
|
vViscous = mobility * tauEff
|
|
|
|
|
|
!* Mean velocity results from waiting time at peierls barriers and solid solution obstacles with respective meanfreepath of
|
|
!* free flight at glide velocity in between.
|
|
!* adopt sign from resolved stress
|
|
|
|
v(s) = sign(1.0_pReal,tau(s)) &
|
|
/ (tPeierls / meanfreepath_P + tSolidSolution / meanfreepath_S + 1.0_pReal / vViscous)
|
|
dv_dtau(s) = v(s) * v(s) * (dtSolidSolution_dtau / meanfreepath_S &
|
|
+ mobility / (vViscous * vViscous))
|
|
dv_dtauNS(s) = v(s) * v(s) * dtPeierls_dtau / meanfreepath_P
|
|
endif
|
|
enddo
|
|
endif
|
|
|
|
|
|
#ifndef _OPENMP
|
|
if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt &
|
|
.and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)&
|
|
.or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then
|
|
write(6,*)
|
|
write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_kinetics at el ip g',el,ip,g
|
|
write(6,*)
|
|
write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> tauThreshold / MPa', tauThreshold / 1e6_pReal
|
|
write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> tau / MPa', tau / 1e6_pReal
|
|
write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> tauNS / MPa', tauNS / 1e6_pReal
|
|
write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> v / 1e-3m/s', v * 1e3
|
|
write(6,'(a,/,12x,12(e12.5,1x))') '<< CONST >> dv_dtau', dv_dtau
|
|
write(6,'(a,/,12x,12(e12.5,1x))') '<< CONST >> dv_dtauNS', dv_dtauNS
|
|
endif
|
|
#endif
|
|
|
|
end subroutine constitutive_nonlocal_kinetics
|
|
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
!> @brief calculates plastic velocity gradient and its tangent
|
|
!--------------------------------------------------------------------------------------------------
|
|
subroutine constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar99, Tstar_v, Temperature, state, g, ip, el)
|
|
|
|
use math, only: math_Plain3333to99, &
|
|
math_mul6x6, &
|
|
math_mul33xx33, &
|
|
math_Mandel6to33
|
|
use debug, only: debug_level, &
|
|
debug_constitutive, &
|
|
debug_levelBasic, &
|
|
debug_levelExtensive, &
|
|
debug_levelSelective, &
|
|
debug_g, &
|
|
debug_i, &
|
|
debug_e
|
|
use material, only: homogenization_maxNgrains, &
|
|
material_phase, &
|
|
phase_plasticityInstance
|
|
use lattice, only: lattice_Sslip, &
|
|
lattice_Sslip_v, &
|
|
lattice_NnonSchmid
|
|
use mesh, only: mesh_ipVolume
|
|
|
|
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) matID, & !< current instance of this plasticity
|
|
structID, & !< 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(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),8) :: &
|
|
rhoSgl !< single dislocation densities (including blocked)
|
|
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),4) :: &
|
|
v, & !< velocity
|
|
tauNS, & !< resolved shear stress including non Schmid and backstress terms
|
|
dv_dtau, & !< velocity derivative with respect to the shear stress
|
|
dv_dtauNS !< velocity derivative with respect to the shear stress
|
|
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el)))) :: &
|
|
tau, & !< resolved shear stress including backstress terms
|
|
gdotTotal, & !< shear rate
|
|
tauBack, & !< back stress from dislocation gradients on same slip system
|
|
tauThreshold !< threshold shear stress
|
|
|
|
|
|
!*** initialize local variables
|
|
|
|
Lp = 0.0_pReal
|
|
dLp_dTstar3333 = 0.0_pReal
|
|
|
|
matID = phase_plasticityInstance(material_phase(g,ip,el))
|
|
structID = constitutive_nonlocal_structure(matID)
|
|
ns = totalNslip(matID)
|
|
|
|
|
|
!*** shortcut to state variables
|
|
|
|
|
|
forall (s = 1_pInt:ns, t = 1_pInt:4_pInt)
|
|
rhoSgl(s,t) = max(state%p(iRhoU(s,t,matID)), 0.0_pReal) ! ensure positive single mobile densities
|
|
rhoSgl(s,t+4_pInt) = state%p(iRhoB(s,t,matID))
|
|
endforall
|
|
where (abs(rhoSgl) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(matID) &
|
|
.or. abs(rhoSgl) < significantRho(matID)) &
|
|
rhoSgl = 0.0_pReal
|
|
|
|
tauBack = state%p(iTauB(1:ns,matID))
|
|
tauThreshold = state%p(iTauF(1:ns,matID))
|
|
|
|
|
|
!*** get resolved shear stress
|
|
!*** for screws possible non-schmid contributions are also taken into account
|
|
|
|
do s = 1_pInt,ns
|
|
sLattice = slipSystemLattice(s,matID)
|
|
tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,structID))
|
|
tauNS(s,1) = tau(s)
|
|
tauNS(s,2) = tau(s)
|
|
if (tau(s) > 0.0_pReal) then
|
|
tauNS(s,3) = math_mul33xx33(math_Mandel6to33(Tstar_v), nonSchmidProjection(1:3,1:3,1,s,matID))
|
|
tauNS(s,4) = math_mul33xx33(math_Mandel6to33(Tstar_v), nonSchmidProjection(1:3,1:3,3,s,matID))
|
|
else
|
|
tauNS(s,3) = math_mul33xx33(math_Mandel6to33(Tstar_v), nonSchmidProjection(1:3,1:3,2,s,matID))
|
|
tauNS(s,4) = math_mul33xx33(math_Mandel6to33(Tstar_v), nonSchmidProjection(1:3,1:3,4,s,matID))
|
|
endif
|
|
enddo
|
|
forall (t = 1_pInt:4_pInt) &
|
|
tauNS(1:ns,t) = tauNS(1:ns,t) + tauBack ! add backstress
|
|
tau = tau + tauBack ! add backstress
|
|
|
|
|
|
!*** get dislocation velocity and its tangent and store the velocity in the state array
|
|
|
|
! edges
|
|
call constitutive_nonlocal_kinetics(v(1:ns,1), dv_dtau(1:ns,1), dv_dtauNS(1:ns,1), &
|
|
tau(1:ns), tauNS(1:ns,1), tauThreshold(1:ns), &
|
|
1_pInt, Temperature, g, ip, el)
|
|
v(1:ns,2) = v(1:ns,1)
|
|
dv_dtau(1:ns,2) = dv_dtau(1:ns,1)
|
|
dv_dtauNS(1:ns,2) = dv_dtauNS(1:ns,1)
|
|
|
|
!screws
|
|
if (lattice_NnonSchmid(structID) == 0_pInt) then ! no non-Schmid contributions
|
|
forall(t = 3_pInt:4_pInt)
|
|
v(1:ns,t) = v(1:ns,1)
|
|
dv_dtau(1:ns,t) = dv_dtau(1:ns,1)
|
|
dv_dtauNS(1:ns,t) = dv_dtauNS(1:ns,1)
|
|
endforall
|
|
else ! take non-Schmid contributions into account
|
|
do t = 3_pInt,4_pInt
|
|
call constitutive_nonlocal_kinetics(v(1:ns,t), dv_dtau(1:ns,t), dv_dtauNS(1:ns,t), &
|
|
tau(1:ns), tauNS(1:ns,t), tauThreshold(1:ns), &
|
|
2_pInt , Temperature, g, ip, el)
|
|
enddo
|
|
endif
|
|
|
|
|
|
!*** store velocity in state
|
|
|
|
forall (t = 1_pInt:4_pInt) &
|
|
state%p(iV(1:ns,t,matID)) = v(1:ns,t)
|
|
|
|
|
|
!*** Bauschinger effect
|
|
|
|
forall (s = 1_pInt:ns, t = 5_pInt:8_pInt, rhoSgl(s,t) * v(s,t-4_pInt) < 0.0_pReal) &
|
|
rhoSgl(s,t-4_pInt) = rhoSgl(s,t-4_pInt) + abs(rhoSgl(s,t))
|
|
|
|
|
|
!*** Calculation of Lp and its tangent
|
|
|
|
gdotTotal = sum(rhoSgl(1:ns,1:4) * v, 2) * burgers(1:ns,matID)
|
|
|
|
do s = 1_pInt,ns
|
|
sLattice = slipSystemLattice(s,matID)
|
|
Lp = Lp + gdotTotal(s) * lattice_Sslip(1:3,1:3,1,sLattice,structID)
|
|
|
|
! Schmid contributions to tangent
|
|
forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt,k=1_pInt:3_pInt,l=1_pInt:3_pInt) &
|
|
dLp_dTstar3333(i,j,k,l) = dLp_dTstar3333(i,j,k,l) &
|
|
+ lattice_Sslip(i,j,1,sLattice,structID) * lattice_Sslip(k,l,1,sLattice,structID) &
|
|
* sum(rhoSgl(s,1:4) * dv_dtau(s,1:4)) * burgers(s,matID)
|
|
|
|
! non Schmid contributions to tangent
|
|
if (tau(s) > 0.0_pReal) then
|
|
forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt,k=1_pInt:3_pInt,l=1_pInt:3_pInt) &
|
|
dLp_dTstar3333(i,j,k,l) = dLp_dTstar3333(i,j,k,l) &
|
|
+ lattice_Sslip(i,j,1,sLattice,structID) &
|
|
* ( nonSchmidProjection(k,l,1,s,matID) * rhoSgl(s,3) * dv_dtauNS(s,3) &
|
|
+ nonSchmidProjection(k,l,3,s,matID) * rhoSgl(s,4) * dv_dtauNS(s,4) ) &
|
|
* burgers(s,matID)
|
|
else
|
|
forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt,k=1_pInt:3_pInt,l=1_pInt:3_pInt) &
|
|
dLp_dTstar3333(i,j,k,l) = dLp_dTstar3333(i,j,k,l) &
|
|
+ lattice_Sslip(i,j,1,sLattice,structID) &
|
|
* ( nonSchmidProjection(k,l,2,s,matID) * rhoSgl(s,3) * dv_dtauNS(s,3) &
|
|
+ nonSchmidProjection(k,l,4,s,matID) * rhoSgl(s,4) * dv_dtauNS(s,4) ) &
|
|
* burgers(s,matID)
|
|
endif
|
|
enddo
|
|
dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333)
|
|
|
|
|
|
#ifndef _OPENMP
|
|
if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt &
|
|
.and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)&
|
|
.or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then
|
|
write(6,*)
|
|
write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_LpandItsTangent at el ip g ',el,ip,g
|
|
write(6,*)
|
|
write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> gdot total / 1e-3',gdotTotal*1e3_pReal
|
|
write(6,'(a,/,3(12x,3(f12.7,1x),/))') '<< CONST >> Lp',transpose(Lp)
|
|
endif
|
|
#endif
|
|
|
|
end subroutine constitutive_nonlocal_LpAndItsTangent
|
|
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
!> @brief (instantaneous) incremental change of microstructure
|
|
!--------------------------------------------------------------------------------------------------
|
|
subroutine constitutive_nonlocal_deltaState(deltaState, state, Tstar_v, Temperature, g,ip,el)
|
|
|
|
use debug, only: debug_level, &
|
|
debug_constitutive, &
|
|
debug_levelBasic, &
|
|
debug_levelExtensive, &
|
|
debug_levelSelective, &
|
|
debug_g, &
|
|
debug_i, &
|
|
debug_e
|
|
use math, only: pi, &
|
|
math_mul6x6
|
|
use lattice, only: lattice_Sslip_v
|
|
use mesh, only: mesh_NcpElems, &
|
|
mesh_maxNips, &
|
|
mesh_ipVolume
|
|
use material, only: homogenization_maxNgrains, &
|
|
material_phase, &
|
|
phase_plasticityInstance
|
|
|
|
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 ! current 2nd Piola-Kirchhoff stress in Mandel notation
|
|
|
|
!*** input/output variables
|
|
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(inout) :: &
|
|
state ! current microstructural state
|
|
|
|
!*** output variables
|
|
type(p_vec), intent(out) :: deltaState ! change of state variables / microstructure
|
|
|
|
!*** local variables
|
|
integer(pInt) matID, & ! current instance of this plasticity
|
|
structID, & ! current lattice structure
|
|
ns, & ! short notation for the total number of active slip systems
|
|
c, & ! character of dislocation
|
|
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(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),10) :: &
|
|
deltaRho, & ! density increment
|
|
deltaRhoRemobilization, & ! density increment by remobilization
|
|
deltaRhoDipole2SingleStress ! density increment by dipole dissociation (by stress change)
|
|
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),8) :: &
|
|
rhoSgl ! current single dislocation densities (positive/negative screw and edge without dipoles)
|
|
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),4) :: &
|
|
v ! dislocation glide velocity
|
|
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el)))) :: &
|
|
tau, & ! current resolved shear stress
|
|
tauBack ! current back stress from pileups on same slip system
|
|
real(pReal), dimension(totalNslip(phase_plasticityInstance(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
|
|
dUpperOld, & ! old maximum stable dipole distance for edges and screws
|
|
deltaDUpper ! change in maximum stable dipole distance for edges and screws
|
|
|
|
|
|
#ifndef _OPENMP
|
|
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt &
|
|
.and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)&
|
|
.or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then
|
|
write(6,*)
|
|
write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_deltaState at el ip g ',el,ip,g
|
|
write(6,*)
|
|
endif
|
|
#endif
|
|
|
|
matID = phase_plasticityInstance(material_phase(g,ip,el))
|
|
structID = constitutive_nonlocal_structure(matID)
|
|
ns = totalNslip(matID)
|
|
|
|
|
|
!*** shortcut to state variables
|
|
|
|
|
|
forall (s = 1_pInt:ns, t = 1_pInt:4_pInt)
|
|
rhoSgl(s,t) = max(state(g,ip,el)%p(iRhoU(s,t,matID)), 0.0_pReal) ! ensure positive single mobile densities
|
|
rhoSgl(s,t+4_pInt) = state(g,ip,el)%p(iRhoB(s,t,matID))
|
|
v(s,t) = state(g,ip,el)%p(iV(s,t,matID))
|
|
endforall
|
|
forall (s = 1_pInt:ns, c = 1_pInt:2_pInt)
|
|
rhoDip(s,c) = max(state(g,ip,el)%p(iRhoD(s,c,matID)), 0.0_pReal) ! ensure positive dipole densities
|
|
dUpperOld(s,c) = state(g,ip,el)%p(iD(s,c,matID))
|
|
endforall
|
|
tauBack = state(g,ip,el)%p(iTauB(1:ns,matID))
|
|
|
|
where (abs(rhoSgl) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(matID) &
|
|
.or. abs(rhoSgl) < significantRho(matID)) &
|
|
rhoSgl = 0.0_pReal
|
|
where (abs(rhoDip) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(matID) &
|
|
.or. abs(rhoDip) < significantRho(matID)) &
|
|
rhoDip = 0.0_pReal
|
|
|
|
|
|
|
|
|
|
!****************************************************************************
|
|
!*** dislocation remobilization (bauschinger effect)
|
|
|
|
deltaRhoRemobilization = 0.0_pReal
|
|
do t = 1_pInt,4_pInt
|
|
do s = 1_pInt,ns
|
|
if (rhoSgl(s,t+4_pInt) * v(s,t) < 0.0_pReal) then
|
|
deltaRhoRemobilization(s,t) = abs(rhoSgl(s,t+4_pInt))
|
|
rhoSgl(s,t) = rhoSgl(s,t) + abs(rhoSgl(s,t+4_pInt))
|
|
deltaRhoRemobilization(s,t+4_pInt) = - rhoSgl(s,t+4_pInt)
|
|
rhoSgl(s,t+4_pInt) = 0.0_pReal
|
|
endif
|
|
enddo
|
|
enddo
|
|
|
|
|
|
|
|
!****************************************************************************
|
|
!*** calculate dipole formation and dissociation by stress change
|
|
|
|
!*** calculate limits for stable dipole height
|
|
|
|
do s = 1_pInt,ns
|
|
sLattice = slipSystemLattice(s,matID)
|
|
tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,structID)) + tauBack(s)
|
|
if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal
|
|
enddo
|
|
dLower = minDipoleHeight(1:ns,1:2,matID)
|
|
dUpper(1:ns,1) = mu(matID) * burgers(1:ns,matID) &
|
|
/ (8.0_pReal * pi * (1.0_pReal - nu(matID)) * abs(tau))
|
|
dUpper(1:ns,2) = mu(matID) * burgers(1:ns,matID) / (4.0_pReal * pi * abs(tau))
|
|
forall (c = 1_pInt:2_pInt) &
|
|
dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) &
|
|
+ abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), &
|
|
dUpper(1:ns,c))
|
|
dUpper = max(dUpper,dLower)
|
|
deltaDUpper = dUpper - dUpperOld
|
|
|
|
|
|
!*** dissociation by stress increase
|
|
|
|
deltaRhoDipole2SingleStress = 0.0_pReal
|
|
forall (c=1_pInt:2_pInt, s=1_pInt:ns, deltaDUpper(s,c) < 0.0_pReal) &
|
|
deltaRhoDipole2SingleStress(s,8_pInt+c) = rhoDip(s,c) * deltaDUpper(s,c) / (dUpperOld(s,c) - dLower(s,c))
|
|
|
|
forall (t=1_pInt:4_pInt) &
|
|
deltaRhoDipole2SingleStress(1_pInt:ns,t) = -0.5_pReal * deltaRhoDipole2SingleStress(1_pInt:ns,(t-1_pInt)/2_pInt+9_pInt)
|
|
|
|
|
|
|
|
!*** store new maximum dipole height in state
|
|
|
|
forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) &
|
|
state(g,ip,el)%p(iD(s,c,matID)) = dUpper(s,c)
|
|
|
|
|
|
|
|
!****************************************************************************
|
|
!*** assign the changes in the dislocation densities to deltaState
|
|
|
|
deltaRho = deltaRhoRemobilization &
|
|
+ deltaRhoDipole2SingleStress
|
|
|
|
deltaState%p = 0.0_pReal
|
|
forall (s = 1:ns, t = 1_pInt:4_pInt)
|
|
deltaState%p(iRhoU(s,t,matID)) = deltaRho(s,t)
|
|
deltaState%p(iRhoB(s,t,matID)) = deltaRho(s,t+4_pInt)
|
|
endforall
|
|
forall (s = 1:ns, c = 1_pInt:2_pInt) &
|
|
deltaState%p(iRhoD(s,c,matID)) = deltaRho(s,c+8_pInt)
|
|
|
|
|
|
#ifndef _OPENMP
|
|
if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt &
|
|
.and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)&
|
|
.or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then
|
|
write(6,'(a,/,8(12x,12(e12.5,1x),/))') '<< CONST >> dislocation remobilization', deltaRhoRemobilization(1:ns,1:8)
|
|
write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> dipole dissociation by stress increase', deltaRhoDipole2SingleStress
|
|
write(6,*)
|
|
endif
|
|
#endif
|
|
|
|
end subroutine constitutive_nonlocal_deltaState
|
|
|
|
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
!> @brief calculates the rate of change of microstructure
|
|
!--------------------------------------------------------------------------------------------------
|
|
function constitutive_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature, state, state0, timestep, subfrac, g,ip,el)
|
|
|
|
use prec, only: DAMASK_NaN
|
|
use numerics, only: numerics_integrationMode, &
|
|
numerics_timeSyncing
|
|
use IO, only: IO_error
|
|
use debug, only: debug_level, &
|
|
debug_constitutive, &
|
|
debug_levelBasic, &
|
|
debug_levelExtensive, &
|
|
debug_levelSelective, &
|
|
debug_g, &
|
|
debug_i, &
|
|
debug_e
|
|
use math, only: math_norm3, &
|
|
math_mul6x6, &
|
|
math_mul3x3, &
|
|
math_mul33x3, &
|
|
math_mul33x33, &
|
|
math_inv33, &
|
|
math_det33, &
|
|
math_transpose33, &
|
|
pi
|
|
use mesh, only: mesh_NcpElems, &
|
|
mesh_maxNips, &
|
|
mesh_element, &
|
|
mesh_ipNeighborhood, &
|
|
mesh_ipVolume, &
|
|
mesh_ipArea, &
|
|
mesh_ipAreaNormal, &
|
|
FE_NipNeighbors, &
|
|
FE_geomtype, &
|
|
FE_celltype
|
|
use material, only: homogenization_maxNgrains, &
|
|
material_phase, &
|
|
phase_plasticityInstance, &
|
|
phase_localPlasticity, &
|
|
phase_plasticity
|
|
use lattice, only: lattice_Sslip_v, &
|
|
lattice_sd, &
|
|
lattice_st
|
|
|
|
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(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
|
|
subfrac !< fraction of timestep at the beginning of the substepped crystallite time increment
|
|
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
|
|
Fe, & !< elastic deformation gradient
|
|
Fp !< plastic deformation gradient
|
|
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
|
|
state, & !< current microstructural state
|
|
state0 !< microstructural state at beginning of crystallite increment
|
|
|
|
!*** input/output variables
|
|
|
|
!*** output variables
|
|
real(pReal), dimension(constitutive_nonlocal_sizeDotState(phase_plasticityInstance(material_phase(g,ip,el)))) :: &
|
|
constitutive_nonlocal_dotState !< evolution of state variables / microstructure
|
|
|
|
!*** local variables
|
|
integer(pInt) matID, & !< current instance of this plasticity
|
|
neighbor_instance, & !< instance of my neighbor's plasticity
|
|
structID, & !< current lattice structure
|
|
ns, & !< short notation for the total number of active slip systems
|
|
c, & !< character of dislocation
|
|
n, & !< index of my current neighbor
|
|
neighbor_el, & !< element number of my neighbor
|
|
neighbor_ip, & !< integration point of my neighbor
|
|
neighbor_n, & !< neighbor index pointing to me when looking from my neighbor
|
|
opposite_neighbor, & !< index of my opposite neighbor
|
|
opposite_ip, & !< ip of my opposite neighbor
|
|
opposite_el, & !< element index of my opposite neighbor
|
|
opposite_n, & !< neighbor index pointing to me when looking from 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
|
|
deads
|
|
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),10) :: &
|
|
rhoDot, & !< density evolution
|
|
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(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),8) :: &
|
|
rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles)
|
|
rhoSglOriginal, &
|
|
neighbor_rhoSgl, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles)
|
|
rhoSgl0, & !< single dislocation densities at start of cryst inc (positive/negative screw and edge without dipoles)
|
|
my_rhoSgl !< single dislocation densities of central ip (positive/negative screw and edge without dipoles)
|
|
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),4) :: &
|
|
v, & !< current dislocation glide velocity
|
|
v0, & !< dislocation glide velocity at start of cryst inc
|
|
my_v, & !< dislocation glide velocity of central ip
|
|
neighbor_v, & !< dislocation glide velocity of enighboring ip
|
|
gdot !< shear rates
|
|
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el)))) :: &
|
|
rhoForest, & !< forest dislocation density
|
|
tauThreshold, & !< threshold shear stress
|
|
tau, & !< current resolved shear stress
|
|
tauBack, & !< current back stress from pileups on same slip system
|
|
vClimb, & !< climb velocity of edge dipoles
|
|
nSources
|
|
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),2) :: &
|
|
rhoDip, & !< current dipole dislocation densities (screw and edge dipoles)
|
|
rhoDipOriginal, &
|
|
dLower, & !< minimum stable dipole distance for edges and screws
|
|
dUpper !< current maximum stable dipole distance for edges and screws
|
|
real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),4) :: &
|
|
m !< direction of dislocation motion
|
|
real(pReal), dimension(3,3) :: my_F, & !< my total deformation gradient
|
|
neighbor_F, & !< total deformation gradient of my neighbor
|
|
my_Fe, & !< my elastic deformation gradient
|
|
neighbor_Fe, & !< elastic deformation gradient of my neighbor
|
|
Favg !< average total deformation gradient of me and my neighbor
|
|
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
|
|
selfDiffusion, & !< self diffusion
|
|
rnd, &
|
|
meshlength
|
|
logical considerEnteringFlux, &
|
|
considerLeavingFlux
|
|
|
|
#ifndef _OPENMP
|
|
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt &
|
|
.and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)&
|
|
.or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then
|
|
write(6,*)
|
|
write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_dotState at el ip g ',el,ip,g
|
|
write(6,*)
|
|
endif
|
|
#endif
|
|
|
|
|
|
matID = phase_plasticityInstance(material_phase(g,ip,el))
|
|
structID = constitutive_nonlocal_structure(matID)
|
|
ns = totalNslip(matID)
|
|
|
|
tau = 0.0_pReal
|
|
gdot = 0.0_pReal
|
|
|
|
|
|
!*** shortcut to state variables
|
|
|
|
|
|
forall (s = 1_pInt:ns, t = 1_pInt:4_pInt)
|
|
rhoSgl(s,t) = max(state(g,ip,el)%p(iRhoU(s,t,matID)), 0.0_pReal) ! ensure positive single mobile densities
|
|
rhoSgl(s,t+4_pInt) = state(g,ip,el)%p(iRhoB(s,t,matID))
|
|
v(s,t) = state(g,ip,el)%p(iV(s,t,matID))
|
|
endforall
|
|
forall (s = 1_pInt:ns, c = 1_pInt:2_pInt)
|
|
rhoDip(s,c) = max(state(g,ip,el)%p(iRhoD(s,c,matID)), 0.0_pReal) ! ensure positive dipole densities
|
|
endforall
|
|
rhoForest = state(g,ip,el)%p(iRhoF(1:ns,matID))
|
|
tauThreshold = state(g,ip,el)%p(iTauF(1:ns,matID))
|
|
tauBack = state(g,ip,el)%p(iTauB(1:ns,matID))
|
|
|
|
rhoSglOriginal = rhoSgl
|
|
rhoDipOriginal = rhoDip
|
|
where (abs(rhoSgl) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(matID) &
|
|
.or. abs(rhoSgl) < significantRho(matID)) &
|
|
rhoSgl = 0.0_pReal
|
|
where (abs(rhoDip) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(matID) &
|
|
.or. abs(rhoDip) < significantRho(matID)) &
|
|
rhoDip = 0.0_pReal
|
|
|
|
if (numerics_timeSyncing) then
|
|
forall (s = 1_pInt:ns, t = 1_pInt:4_pInt)
|
|
rhoSgl0(s,t) = max(state0(g,ip,el)%p(iRhoU(s,t,matID)), 0.0_pReal)
|
|
rhoSgl0(s,t+4_pInt) = state0(g,ip,el)%p(iRhoB(s,t,matID))
|
|
v0(s,t) = state0(g,ip,el)%p(iV(s,t,matID))
|
|
endforall
|
|
where (abs(rhoSgl0) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(matID) &
|
|
.or. abs(rhoSgl0) < significantRho(matID)) &
|
|
rhoSgl0 = 0.0_pReal
|
|
endif
|
|
|
|
|
|
|
|
!*** sanity check for timestep
|
|
|
|
if (timestep <= 0.0_pReal) then ! if illegal timestep...
|
|
constitutive_nonlocal_dotState = 0.0_pReal ! ...return without doing anything (-> zero dotState)
|
|
return
|
|
endif
|
|
|
|
|
|
|
|
!****************************************************************************
|
|
!*** Calculate shear rate
|
|
|
|
forall (t = 1_pInt:4_pInt) &
|
|
gdot(1_pInt:ns,t) = rhoSgl(1_pInt:ns,t) * burgers(1:ns,matID) * v(1:ns,t)
|
|
|
|
#ifndef _OPENMP
|
|
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt &
|
|
.and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)&
|
|
.or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then
|
|
write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> rho / 1/m^2', rhoSgl, rhoDip
|
|
write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> gdot / 1/s',gdot
|
|
endif
|
|
#endif
|
|
|
|
|
|
|
|
!****************************************************************************
|
|
!*** calculate limits for stable dipole height
|
|
|
|
do s = 1_pInt,ns ! loop over slip systems
|
|
sLattice = slipSystemLattice(s,matID)
|
|
tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,structID)) + tauBack(s)
|
|
if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal
|
|
enddo
|
|
|
|
dLower = minDipoleHeight(1:ns,1:2,matID)
|
|
dUpper(1:ns,1) = mu(matID) * burgers(1:ns,matID) &
|
|
/ (8.0_pReal * pi * (1.0_pReal - nu(matID)) * abs(tau))
|
|
dUpper(1:ns,2) = mu(matID) * burgers(1:ns,matID) &
|
|
/ (4.0_pReal * pi * abs(tau))
|
|
forall (c = 1_pInt:2_pInt) &
|
|
dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) &
|
|
+ abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), &
|
|
dUpper(1:ns,c))
|
|
dUpper = max(dUpper,dLower)
|
|
|
|
|
|
|
|
!****************************************************************************
|
|
!*** calculate dislocation multiplication
|
|
|
|
rhoDotMultiplication = 0.0_pReal
|
|
if (structID == 2_pInt) then ! BCC
|
|
forall (s = 1:ns, sum(abs(v(s,1:4))) > 0.0_pReal)
|
|
rhoDotMultiplication(s,1:2) = sum(abs(gdot(s,3:4))) / burgers(s,matID) & ! assuming double-cross-slip of screws to be decisive for multiplication
|
|
* sqrt(rhoForest(s)) / lambda0(s,matID) ! & ! mean free path
|
|
! * 2.0_pReal * sum(abs(v(s,3:4))) / sum(abs(v(s,1:4))) ! ratio of screw to overall velocity determines edge generation
|
|
rhoDotMultiplication(s,3:4) = sum(abs(gdot(s,3:4))) / burgers(s,matID) & ! assuming double-cross-slip of screws to be decisive for multiplication
|
|
* sqrt(rhoForest(s)) / lambda0(s,matID) ! & ! mean free path
|
|
! * 2.0_pReal * sum(abs(v(s,1:2))) / sum(abs(v(s,1:4))) ! ratio of edge to overall velocity determines screw generation
|
|
endforall
|
|
|
|
else ! ALL OTHER STRUCTURES
|
|
if (probabilisticMultiplication(matID)) then
|
|
meshlength = mesh_ipVolume(ip,el)**0.333_pReal
|
|
where(sum(rhoSgl(1:ns,1:4),2) > 0.0_pReal)
|
|
nSources = (sum(rhoSgl(1:ns,1:2),2) * fEdgeMultiplication(matID) + sum(rhoSgl(1:ns,3:4),2)) &
|
|
/ sum(rhoSgl(1:ns,1:4),2) * meshlength / lambda0(1:ns,matID) * sqrt(rhoForest(1:ns))
|
|
elsewhere
|
|
nSources = meshlength / lambda0(1:ns,matID) * sqrt(rhoForest(1:ns))
|
|
endwhere
|
|
do s = 1_pInt,ns
|
|
if (nSources(s) < 1.0_pReal) then
|
|
if (sourceProbability(s,g,ip,el) > 1.0_pReal) then
|
|
call random_number(rnd)
|
|
sourceProbability(s,g,ip,el) = rnd
|
|
!$OMP FLUSH(sourceProbability)
|
|
endif
|
|
if (sourceProbability(s,g,ip,el) > 1.0_pReal - nSources(s)) then
|
|
rhoDotMultiplication(s,1:4) = sum(rhoSglOriginal(s,1:4) * abs(v(s,1:4))) / meshlength
|
|
endif
|
|
else
|
|
sourceProbability(s,g,ip,el) = 2.0_pReal
|
|
rhoDotMultiplication(s,1:4) = &
|
|
(sum(abs(gdot(s,1:2))) * fEdgeMultiplication(matID) + sum(abs(gdot(s,3:4)))) &
|
|
/ burgers(s,matID) * sqrt(rhoForest(s)) / lambda0(s,matID)
|
|
endif
|
|
enddo
|
|
#ifndef _OPENMP
|
|
if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt &
|
|
.and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)&
|
|
.or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then
|
|
write(6,'(a,/,4(12x,12(f12.5,1x),/))') '<< CONST >> sources', nSources
|
|
write(6,*)
|
|
endif
|
|
#endif
|
|
else
|
|
rhoDotMultiplication(1:ns,1:4) = spread( &
|
|
(sum(abs(gdot(1:ns,1:2)),2) * fEdgeMultiplication(matID) + sum(abs(gdot(1:ns,3:4)),2)) &
|
|
* sqrt(rhoForest(1:ns)) / lambda0(1:ns,matID) / burgers(1:ns,matID), 2, 4)
|
|
endif
|
|
endif
|
|
|
|
|
|
|
|
!****************************************************************************
|
|
!*** calculate dislocation fluxes (only for nonlocal plasticity)
|
|
|
|
rhoDotFlux = 0.0_pReal
|
|
|
|
if (.not. phase_localPlasticity(material_phase(g,ip,el))) then ! only for nonlocal plasticity
|
|
|
|
|
|
!*** check CFL (Courant-Friedrichs-Lewy) condition for flux
|
|
|
|
if (any( abs(gdot) > 0.0_pReal & ! any active slip system ...
|
|
.and. CFLfactor(matID) * abs(v) * timestep &
|
|
> mesh_ipVolume(ip,el) / maxval(mesh_ipArea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here)
|
|
#ifndef _OPENMP
|
|
if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt) then
|
|
write(6,'(a,i5,a,i2)') '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip
|
|
write(6,'(a,e10.3,a,e10.3)') '<< CONST >> velocity is at ', &
|
|
maxval(abs(v), abs(gdot) > 0.0_pReal &
|
|
.and. CFLfactor(matID) * abs(v) * timestep &
|
|
> mesh_ipVolume(ip,el) / maxval(mesh_ipArea(:,ip,el))), &
|
|
' at a timestep of ',timestep
|
|
write(6,'(a)') '<< CONST >> enforcing cutback !!!'
|
|
endif
|
|
#endif
|
|
constitutive_nonlocal_dotState = DAMASK_NaN ! -> return NaN and, hence, enforce cutback
|
|
return
|
|
endif
|
|
|
|
|
|
!*** be aware 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, slipSystemLattice(1:ns,matID), structID)
|
|
m(1:3,1:ns,2) = -lattice_sd(1:3, slipSystemLattice(1:ns,matID), structID)
|
|
m(1:3,1:ns,3) = -lattice_st(1:3, slipSystemLattice(1:ns,matID), structID)
|
|
m(1:3,1:ns,4) = lattice_st(1:3, slipSystemLattice(1:ns,matID), structID)
|
|
|
|
my_Fe = Fe(1:3,1:3,g,ip,el)
|
|
my_F = math_mul33x33(my_Fe, Fp(1:3,1:3,g,ip,el))
|
|
|
|
do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el)))) ! loop through my neighbors
|
|
neighbor_el = mesh_ipNeighborhood(1,n,ip,el)
|
|
neighbor_ip = mesh_ipNeighborhood(2,n,ip,el)
|
|
neighbor_n = mesh_ipNeighborhood(3,n,ip,el)
|
|
|
|
opposite_neighbor = n + mod(n,2_pInt) - mod(n+1_pInt,2_pInt)
|
|
opposite_el = mesh_ipNeighborhood(1,opposite_neighbor,ip,el)
|
|
opposite_ip = mesh_ipNeighborhood(2,opposite_neighbor,ip,el)
|
|
opposite_n = mesh_ipNeighborhood(3,opposite_neighbor,ip,el)
|
|
|
|
if (neighbor_n > 0_pInt) then ! if neighbor exists, average deformation gradient
|
|
neighbor_instance = phase_plasticityInstance(material_phase(g,neighbor_ip,neighbor_el))
|
|
neighbor_Fe = Fe(1:3,1:3,g,neighbor_ip,neighbor_el)
|
|
neighbor_F = math_mul33x33(neighbor_Fe, Fp(1:3,1:3,g,neighbor_ip,neighbor_el))
|
|
Favg = 0.5_pReal * (my_F + neighbor_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 plasticity (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.
|
|
neighbor_v = 0.0_pReal ! needed for check of sign change in flux density below
|
|
neighbor_rhoSgl = 0.0_pReal
|
|
if (neighbor_n > 0_pInt) then
|
|
if (phase_plasticity(material_phase(1,neighbor_ip,neighbor_el)) == CONSTITUTIVE_NONLOCAL_LABEL &
|
|
.and. any(compatibility(:,:,:,n,ip,el) > 0.0_pReal)) &
|
|
considerEnteringFlux = .true.
|
|
endif
|
|
|
|
if (considerEnteringFlux) then
|
|
if(numerics_timeSyncing .and. (subfrac(g,neighbor_ip,neighbor_el) /= subfrac(g,ip,el))) then ! for timesyncing: in case of a timestep at the interface we have to use "state0" to make sure that fluxes n both sides are equal
|
|
forall (s = 1:ns, t = 1_pInt:4_pInt)
|
|
neighbor_v(s,t) = state0(g,neighbor_ip,neighbor_el)%p(iV(s,t,neighbor_instance))
|
|
neighbor_rhoSgl(s,t) = max(state0(g,neighbor_ip,neighbor_el)%p(iRhoU(s,t,neighbor_instance)), 0.0_pReal)
|
|
neighbor_rhoSgl(s,t+4_pInt) = state0(g,neighbor_ip,neighbor_el)%p(iRhoB(s,t,neighbor_instance))
|
|
endforall
|
|
else
|
|
forall (s = 1:ns, t = 1_pInt:4_pInt)
|
|
neighbor_v(s,t) = state(g,neighbor_ip,neighbor_el)%p(iV(s,t,neighbor_instance))
|
|
neighbor_rhoSgl(s,t) = max(state(g,neighbor_ip,neighbor_el)%p(iRhoU(s,t,neighbor_instance)), 0.0_pReal)
|
|
neighbor_rhoSgl(s,t+4_pInt) = state(g,neighbor_ip,neighbor_el)%p(iRhoB(s,t,neighbor_instance))
|
|
endforall
|
|
endif
|
|
where (abs(neighbor_rhoSgl) * mesh_ipVolume(neighbor_ip,neighbor_el) ** 0.667_pReal &
|
|
< significantN(matID) &
|
|
.or. abs(neighbor_rhoSgl) < significantRho(matID)) &
|
|
neighbor_rhoSgl = 0.0_pReal
|
|
normal_neighbor2me_defConf = math_det33(Favg) * math_mul33x3(math_inv33(transpose(Favg)), &
|
|
mesh_ipAreaNormal(1:3,neighbor_n,neighbor_ip,neighbor_el)) ! calculate the normal of the interface in (average) deformed configuration (now pointing from my neighbor to me!!!)
|
|
normal_neighbor2me = math_mul33x3(transpose(neighbor_Fe), normal_neighbor2me_defConf) &
|
|
/ math_det33(neighbor_Fe) ! interface normal in the lattice configuration of my neighbor
|
|
area = mesh_ipArea(neighbor_n,neighbor_ip,neighbor_el) * math_norm3(normal_neighbor2me)
|
|
normal_neighbor2me = normal_neighbor2me / math_norm3(normal_neighbor2me) ! normalize the surface normal to unit length
|
|
do s = 1_pInt,ns
|
|
do t = 1_pInt,4_pInt
|
|
c = (t + 1_pInt) / 2
|
|
topp = t + mod(t,2_pInt) - mod(t+1_pInt,2_pInt)
|
|
if (neighbor_v(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. v(s,t) * neighbor_v(s,t) > 0.0_pReal ) then ! ... only if no sign change in flux density
|
|
do deads = 0_pInt,4_pInt,4_pInt
|
|
lineLength = abs(neighbor_rhoSgl(s,t+deads)) * neighbor_v(s,t) &
|
|
* math_mul3x3(m(1:3,s,t), normal_neighbor2me) * area ! positive line length that wants to enter through this interface
|
|
where (compatibility(c,1_pInt:ns,s,n,ip,el) > 0.0_pReal) & ! positive compatibility...
|
|
rhoDotFlux(1_pInt:ns,t) = rhoDotFlux(1_pInt:ns,t) &
|
|
+ lineLength / mesh_ipVolume(ip,el) & ! ... transferring to equally signed mobile dislocation type
|
|
* compatibility(c,1_pInt:ns,s,n,ip,el) ** 2.0_pReal
|
|
where (compatibility(c,1_pInt:ns,s,n,ip,el) < 0.0_pReal) & ! ..negative compatibility...
|
|
rhoDotFlux(1_pInt:ns,topp) = rhoDotFlux(1_pInt:ns,topp) &
|
|
+ lineLength / mesh_ipVolume(ip,el) & ! ... transferring to opposite signed mobile dislocation type
|
|
* compatibility(c,1_pInt:ns,s,n,ip,el) ** 2.0_pReal
|
|
enddo
|
|
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_n > 0_pInt) then
|
|
if (phase_plasticity(material_phase(1,opposite_ip,opposite_el)) /= CONSTITUTIVE_NONLOCAL_LABEL) &
|
|
considerLeavingFlux = .false.
|
|
endif
|
|
|
|
if (considerLeavingFlux) then
|
|
|
|
!* timeSyncing mode: If the central ip has zero subfraction, always use "state0". This is needed in case of
|
|
!* a synchronization step for the central ip, because then "state" contains the values at the end of the
|
|
!* previously converged full time step. Also, if either me or my neighbor has zero subfraction, we have to
|
|
!* use "state0" to make sure that fluxes on both sides of the (potential) timestep are equal.
|
|
my_rhoSgl = rhoSgl
|
|
my_v = v
|
|
if(numerics_timeSyncing) then
|
|
if (subfrac(g,ip,el) == 0.0_pReal) then
|
|
my_rhoSgl = rhoSgl0
|
|
my_v = v0
|
|
elseif (neighbor_n > 0_pInt) then
|
|
if (subfrac(g,neighbor_ip,neighbor_el) == 0.0_pReal) then
|
|
my_rhoSgl = rhoSgl0
|
|
my_v = v0
|
|
endif
|
|
endif
|
|
endif
|
|
|
|
normal_me2neighbor_defConf = math_det33(Favg) &
|
|
* math_mul33x3(math_inv33(math_transpose33(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(math_transpose33(my_Fe), normal_me2neighbor_defConf) &
|
|
/ math_det33(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_pInt,ns
|
|
do t = 1_pInt,4_pInt
|
|
c = (t + 1_pInt) / 2_pInt
|
|
if (my_v(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)
|
|
if (my_v(s,t) * neighbor_v(s,t) > 0.0_pReal) then ! no sign change in flux density
|
|
transmissivity = sum(compatibility(c,1_pInt: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
|
|
lineLength = my_rhoSgl(s,t) * my_v(s,t) &
|
|
* math_mul3x3(m(1:3,s,t), normal_me2neighbor) * area ! positive line length of mobiles that wants to leave through this interface
|
|
rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / mesh_ipVolume(ip,el) ! subtract dislocation flux from current type
|
|
rhoDotFlux(s,t+4_pInt) = rhoDotFlux(s,t+4_pInt) &
|
|
+ lineLength / mesh_ipVolume(ip,el) * (1.0_pReal - transmissivity) &
|
|
* sign(1.0_pReal, my_v(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
|
|
lineLength = my_rhoSgl(s,t+4_pInt) * my_v(s,t) &
|
|
* math_mul3x3(m(1:3,s,t), normal_me2neighbor) * area ! positive line length of deads that wants to leave through this interface
|
|
rhoDotFlux(s,t+4_pInt) = rhoDotFlux(s,t+4_pInt) &
|
|
- lineLength / mesh_ipVolume(ip,el) * transmissivity ! dead dislocations leaving through this interface
|
|
endif
|
|
enddo
|
|
enddo
|
|
endif
|
|
|
|
enddo ! neighbor loop
|
|
endif
|
|
|
|
|
|
|
|
!****************************************************************************
|
|
!*** calculate dipole formation and annihilation
|
|
|
|
!*** formation by glide
|
|
|
|
do c = 1_pInt,2_pInt
|
|
|
|
rhoDotSingle2DipoleGlide(1:ns,2*c-1) = -2.0_pReal * dUpper(1:ns,c) / burgers(1:ns,matID) &
|
|
* (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) / burgers(1:ns,matID) &
|
|
* (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) / burgers(1:ns,matID) &
|
|
* 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) / burgers(1:ns,matID) &
|
|
* 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_pInt:2_pInt) &
|
|
rhoDotAthermalAnnihilation(1:ns,c+8_pInt) = -2.0_pReal * dLower(1:ns,c) / burgers(1:ns,matID) &
|
|
* ( 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
|
|
! annihilated screw dipoles leave edge jogs behind on the colinear system
|
|
if (structID == 1_pInt) then ! only fcc
|
|
forall (s = 1:ns, colinearSystem(s,matID) > 0_pInt) &
|
|
rhoDotAthermalAnnihilation(colinearSystem(s,matID),1:2) = - rhoDotAthermalAnnihilation(s,10) &
|
|
* 0.25_pReal * sqrt(rhoForest(s)) * (dUpper(s,2) + dLower(s,2)) * edgeJogFactor(matID)
|
|
endif
|
|
|
|
|
|
!*** thermally activated annihilation of edge dipoles by climb
|
|
|
|
rhoDotThermalAnnihilation = 0.0_pReal
|
|
selfDiffusion = Dsd0(matID) * exp(-selfDiffusionEnergy(matID) / (KB * Temperature))
|
|
vClimb = atomicVolume(matID) * selfDiffusion / ( KB * Temperature ) &
|
|
* mu(matID) / ( 2.0_pReal * PI * (1.0_pReal-nu(matID)) ) &
|
|
* 2.0_pReal / ( dUpper(1:ns,1) + dLower(1:ns,1) )
|
|
forall (s = 1_pInt:ns, dUpper(s,1) > dLower(s,1)) &
|
|
rhoDotThermalAnnihilation(s,9) = max(- 4.0_pReal * rhoDip(s,1) * vClimb(s) / (dUpper(s,1) - dLower(s,1)), &
|
|
- rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) - rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have
|
|
|
|
|
|
|
|
!****************************************************************************
|
|
!*** assign the rates of dislocation densities to my dotState
|
|
!*** if evolution rates lead to negative densities, a cutback is enforced
|
|
|
|
rhoDot = 0.0_pReal
|
|
rhoDot = rhoDotFlux &
|
|
+ rhoDotMultiplication &
|
|
+ rhoDotSingle2DipoleGlide &
|
|
+ rhoDotAthermalAnnihilation &
|
|
+ rhoDotThermalAnnihilation
|
|
|
|
if (numerics_integrationMode == 1_pInt) then ! save rates for output if in central integration mode
|
|
rhoDotFluxOutput(1:ns,1:8,g,ip,el) = rhoDotFlux(1:ns,1:8)
|
|
rhoDotMultiplicationOutput(1:ns,1:2,g,ip,el) = rhoDotMultiplication(1:ns,[1,3])
|
|
rhoDotSingle2DipoleGlideOutput(1:ns,1:2,g,ip,el) = rhoDotSingle2DipoleGlide(1:ns,9:10)
|
|
rhoDotAthermalAnnihilationOutput(1:ns,1:2,g,ip,el) = rhoDotAthermalAnnihilation(1:ns,9:10)
|
|
rhoDotThermalAnnihilationOutput(1:ns,1:2,g,ip,el) = rhoDotThermalAnnihilation(1:ns,9:10)
|
|
rhoDotEdgeJogsOutput(1:ns,g,ip,el) = 2.0_pReal * rhoDotThermalAnnihilation(1:ns,1)
|
|
endif
|
|
|
|
|
|
#ifndef _OPENMP
|
|
if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt &
|
|
.and. ((debug_e == el .and. debug_i == ip .and. debug_g == g)&
|
|
.or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then
|
|
write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> dislocation multiplication', rhoDotMultiplication(1:ns,1:4) * timestep
|
|
write(6,'(a,/,8(12x,12(e12.5,1x),/))') '<< CONST >> dislocation flux', rhoDotFlux(1:ns,1:8) * timestep
|
|
write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> dipole formation by glide', rhoDotSingle2DipoleGlide * timestep
|
|
write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> athermal dipole annihilation', &
|
|
rhoDotAthermalAnnihilation * timestep
|
|
write(6,'(a,/,2(12x,12(e12.5,1x),/))') '<< CONST >> thermally activated dipole annihilation', &
|
|
rhoDotThermalAnnihilation(1:ns,9:10) * timestep
|
|
write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> total density change', rhoDot * timestep
|
|
write(6,'(a,/,10(12x,12(f12.5,1x),/))') '<< CONST >> relative density change', &
|
|
rhoDot(1:ns,1:8) * timestep / (abs(rhoSglOriginal)+1.0e-10), &
|
|
rhoDot(1:ns,9:10) * timestep / (rhoDipOriginal+1.0e-10)
|
|
write(6,*)
|
|
endif
|
|
#endif
|
|
|
|
|
|
if ( any(rhoSglOriginal(1:ns,1:4) + rhoDot(1:ns,1:4) * timestep < -aTolRho(matID)) &
|
|
.or. any(rhoDipOriginal(1:ns,1:2) + rhoDot(1:ns,9:10) * timestep < -aTolRho(matID))) then
|
|
#ifndef _OPENMP
|
|
if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt) then
|
|
write(6,'(a,i5,a,i2)') '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip
|
|
write(6,'(a)') '<< CONST >> enforcing cutback !!!'
|
|
endif
|
|
#endif
|
|
constitutive_nonlocal_dotState = DAMASK_NaN
|
|
return
|
|
else
|
|
forall (s = 1:ns, t = 1_pInt:4_pInt)
|
|
constitutive_nonlocal_dotState(iRhoU(s,t,matID)) = rhoDot(s,t)
|
|
constitutive_nonlocal_dotState(iRhoB(s,t,matID)) = rhoDot(s,t+4_pInt)
|
|
endforall
|
|
forall (s = 1:ns, c = 1_pInt:2_pInt) &
|
|
constitutive_nonlocal_dotState(iRhoD(s,c,matID)) = rhoDot(s,c+8_pInt)
|
|
forall (s = 1:ns) &
|
|
constitutive_nonlocal_dotState(iGamma(s,matID)) = sum(gdot(s,1:4))
|
|
endif
|
|
|
|
end function constitutive_nonlocal_dotState
|
|
|
|
|
|
|
|
!*********************************************************************
|
|
!* 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 math, only: math_qDisorientation, &
|
|
math_mul3x3, &
|
|
math_qRot
|
|
use material, only: material_phase, &
|
|
material_texture, &
|
|
phase_localPlasticity, &
|
|
phase_plasticityInstance, &
|
|
homogenization_maxNgrains
|
|
use mesh, only: mesh_element, &
|
|
mesh_ipNeighborhood, &
|
|
mesh_maxNips, &
|
|
mesh_NcpElems, &
|
|
FE_NipNeighbors, &
|
|
FE_geomtype, &
|
|
FE_celltype
|
|
use lattice, only: lattice_sn, &
|
|
lattice_sd
|
|
|
|
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
|
|
neighbor_e, & ! element index of my neighbor
|
|
neighbor_i, & ! integration point index of my neighbor
|
|
my_phase, &
|
|
neighbor_phase, &
|
|
my_texture, &
|
|
neighbor_texture, &
|
|
my_structure, & ! lattice structure
|
|
my_instance, & ! instance of plasticity
|
|
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,totalNslip(phase_plasticityInstance(material_phase(1,i,e))),&
|
|
totalNslip(phase_plasticityInstance(material_phase(1,i,e))),&
|
|
FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,e))))) :: &
|
|
my_compatibility ! my_compatibility for current element and ip
|
|
real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(1,i,e)))) :: &
|
|
slipNormal, &
|
|
slipDirection
|
|
real(pReal) my_compatibilitySum, &
|
|
thresholdValue, &
|
|
nThresholdValues
|
|
logical, dimension(totalNslip(phase_plasticityInstance(material_phase(1,i,e)))) :: &
|
|
belowThreshold
|
|
|
|
|
|
Nneighbors = FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,e))))
|
|
my_phase = material_phase(1,i,e)
|
|
my_texture = material_texture(1,i,e)
|
|
my_instance = phase_plasticityInstance(my_phase)
|
|
my_structure = constitutive_nonlocal_structure(my_instance)
|
|
ns = totalNslip(my_instance)
|
|
slipNormal(1:3,1:ns) = lattice_sn(1:3, slipSystemLattice(1:ns,my_instance), my_structure)
|
|
slipDirection(1:3,1:ns) = lattice_sd(1:3, slipSystemLattice(1:ns,my_instance), my_structure)
|
|
|
|
|
|
!*** start out fully compatible
|
|
|
|
my_compatibility = 0.0_pReal
|
|
forall(s1 = 1_pInt:ns) &
|
|
my_compatibility(1:2,s1,s1,1:Nneighbors) = 1.0_pReal
|
|
|
|
|
|
!*** Loop thrugh neighbors and check whether there is any my_compatibility.
|
|
|
|
do n = 1_pInt,Nneighbors
|
|
neighbor_e = mesh_ipNeighborhood(1,n,i,e)
|
|
neighbor_i = mesh_ipNeighborhood(2,n,i,e)
|
|
|
|
|
|
!* FREE SURFACE
|
|
!* Set surface transmissivity to the value specified in the material.config
|
|
|
|
if (neighbor_e <= 0_pInt .or. neighbor_i <= 0_pInt) then
|
|
forall(s1 = 1_pInt:ns) &
|
|
my_compatibility(1:2,s1,s1,n) = sqrt(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 one of the two "CPFEM" phases has a local plasticity law,
|
|
!* we do not consider this to be a phase boundary, so completely compatible.
|
|
|
|
neighbor_phase = material_phase(1,neighbor_i,neighbor_e)
|
|
if (neighbor_phase /= my_phase) then
|
|
if (.not. phase_localPlasticity(neighbor_phase) .and. .not. phase_localPlasticity(my_phase)) then
|
|
forall(s1 = 1_pInt:ns) &
|
|
my_compatibility(1:2,s1,s1,n) = 0.0_pReal ! = sqrt(0.0)
|
|
endif
|
|
cycle
|
|
endif
|
|
|
|
|
|
!* GRAIN BOUNDARY !
|
|
!* fixed transmissivity for adjacent ips with different texture (only if explicitly given in material.config)
|
|
|
|
if (grainboundaryTransmissivity(my_instance) >= 0.0_pReal) then
|
|
neighbor_texture = material_texture(1,neighbor_i,neighbor_e)
|
|
if (neighbor_texture /= my_texture) then
|
|
if (.not. phase_localPlasticity(neighbor_phase)) then
|
|
forall(s1 = 1_pInt:ns) &
|
|
my_compatibility(1:2,s1,s1,n) = sqrt(grainboundaryTransmissivity(my_instance))
|
|
endif
|
|
cycle
|
|
endif
|
|
|
|
|
|
!* GRAIN BOUNDARY ?
|
|
!* Compatibility defined by relative orientation of slip systems:
|
|
!* The my_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 my_compatibility values exceeding one.
|
|
!* Finally the smallest my_compatibility value is decreased until the sum is exactly equal to one.
|
|
!* All values below the threshold are set to zero.
|
|
else
|
|
absoluteMisorientation = math_qDisorientation(orientation(1:4,1,i,e), &
|
|
orientation(1:4,1,neighbor_i,neighbor_e), &
|
|
0_pInt) ! no symmetry
|
|
do s1 = 1_pInt,ns ! my slip systems
|
|
do s2 = 1_pInt,ns ! my neighbor's slip systems
|
|
my_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))))
|
|
my_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
|
|
|
|
my_compatibilitySum = 0.0_pReal
|
|
belowThreshold = .true.
|
|
do while (my_compatibilitySum < 1.0_pReal .and. any(belowThreshold(1:ns)))
|
|
thresholdValue = maxval(my_compatibility(2,1:ns,s1,n), belowThreshold(1:ns)) ! screws always positive
|
|
nThresholdValues = real(count(my_compatibility(2,1:ns,s1,n) == thresholdValue),pReal)
|
|
where (my_compatibility(2,1:ns,s1,n) >= thresholdValue) &
|
|
belowThreshold(1:ns) = .false.
|
|
if (my_compatibilitySum + thresholdValue * nThresholdValues > 1.0_pReal) &
|
|
where (abs(my_compatibility(1:2,1:ns,s1,n)) == thresholdValue) &
|
|
my_compatibility(1:2,1:ns,s1,n) = sign((1.0_pReal - my_compatibilitySum) &
|
|
/ nThresholdValues, my_compatibility(1:2,1:ns,s1,n))
|
|
my_compatibilitySum = my_compatibilitySum + nThresholdValues * thresholdValue
|
|
enddo
|
|
where (belowThreshold(1:ns)) my_compatibility(1,1:ns,s1,n) = 0.0_pReal
|
|
where (belowThreshold(1:ns)) my_compatibility(2,1:ns,s1,n) = 0.0_pReal
|
|
enddo ! my slip systems cycle
|
|
endif
|
|
|
|
enddo ! neighbor cycle
|
|
|
|
compatibility(1:2,1:ns,1:ns,1:Nneighbors,i,e) = my_compatibility
|
|
|
|
end subroutine constitutive_nonlocal_updateCompatibility
|
|
|
|
|
|
!*********************************************************************
|
|
!* calculates quantities characterizing the microstructure *
|
|
!*********************************************************************
|
|
function constitutive_nonlocal_dislocationstress(state, Fe, g, ip, el)
|
|
|
|
use math, only: math_mul33x33, &
|
|
math_mul33x3, &
|
|
math_invert33, &
|
|
math_transpose33, &
|
|
pi
|
|
use mesh, only: mesh_NcpElems, &
|
|
mesh_maxNips, &
|
|
mesh_element, &
|
|
mesh_node0, &
|
|
mesh_cellCenterCoordinates, &
|
|
mesh_ipVolume, &
|
|
mesh_periodicSurface, &
|
|
FE_Nips, &
|
|
FE_geomtype
|
|
use material, only: homogenization_maxNgrains, &
|
|
material_phase, &
|
|
phase_localPlasticity, &
|
|
phase_plasticityInstance
|
|
|
|
implicit none
|
|
|
|
|
|
!*** input variables
|
|
integer(pInt), intent(in) :: g, & ! current grain ID
|
|
ip, & ! current integration point
|
|
el ! current element
|
|
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
|
|
Fe ! elastic deformation gradient
|
|
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
|
|
state ! microstructural state
|
|
|
|
!*** input/output variables
|
|
|
|
!*** output variables
|
|
real(pReal), dimension(3,3) :: constitutive_nonlocal_dislocationstress
|
|
|
|
!*** local variables
|
|
integer(pInt) neighbor_el, & ! element number of neighbor material point
|
|
neighbor_ip, & ! integration point of neighbor material point
|
|
instance, & ! my instance of this plasticity
|
|
neighbor_instance, & ! instance of this plasticity of neighbor material point
|
|
latticeStruct, & ! my lattice structure
|
|
neighbor_latticeStruct, & ! lattice structure of neighbor material point
|
|
phase, &
|
|
neighbor_phase, &
|
|
ns, & ! total number of active slip systems at my material point
|
|
neighbor_ns, & ! total number of active slip systems at neighbor 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) x, y, z, & ! coordinates of connection vector in neighbor 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, &
|
|
neighbor_ipVolumeSideLength, &
|
|
detFe
|
|
real(pReal), dimension(3) :: connection, & ! connection vector between me and my neighbor in the deformed configuration
|
|
connection_neighborLattice, & ! connection vector between me and my neighbor in the lattice configuration of my neighbor
|
|
connection_neighborSlip, & ! connection vector between me and my neighbor in the slip system frame of my neighbor
|
|
maxCoord, minCoord, &
|
|
meshSize, &
|
|
coords, & ! x,y,z coordinates of cell center of ip volume
|
|
neighbor_coords ! x,y,z coordinates of cell center of neighbor ip volume
|
|
real(pReal), dimension(3,3) :: sigma, & ! dislocation stress for one slip system in neighbor material point's slip system frame
|
|
Tdislo_neighborLattice, & ! dislocation stress as 2nd Piola-Kirchhoff stress at neighbor material point
|
|
invFe, & ! inverse of my elastic deformation gradient
|
|
neighbor_invFe, &
|
|
neighborLattice2myLattice ! mapping from neighbor MPs lattice configuration to my lattice configuration
|
|
real(pReal), dimension(2,2,maxval(totalNslip)) :: &
|
|
neighbor_rhoExcess ! excess density at neighbor material point (edge/screw,mobile/dead,slipsystem)
|
|
real(pReal), dimension(2,maxval(totalNslip)) :: &
|
|
rhoExcessDead
|
|
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),8) :: &
|
|
rhoSgl ! single dislocation density (edge+, edge-, screw+, screw-, used edge+, used edge-, used screw+, used screw-)
|
|
logical inversionError
|
|
|
|
phase = material_phase(g,ip,el)
|
|
instance = phase_plasticityInstance(phase)
|
|
latticeStruct = constitutive_nonlocal_structure(instance)
|
|
ns = totalNslip(instance)
|
|
|
|
|
|
|
|
!*** get basic states
|
|
|
|
forall (s = 1_pInt:ns, t = 1_pInt:4_pInt)
|
|
rhoSgl(s,t) = max(state(g,ip,el)%p(iRhoU(s,t,instance)), 0.0_pReal) ! ensure positive single mobile densities
|
|
rhoSgl(s,t+4_pInt) = state(g,ip,el)%p(iRhoB(s,t,instance))
|
|
endforall
|
|
|
|
|
|
|
|
!*** calculate the dislocation stress of the neighboring excess dislocation densities
|
|
!*** zero for material points of local plasticity
|
|
|
|
constitutive_nonlocal_dislocationstress = 0.0_pReal
|
|
|
|
if (.not. phase_localPlasticity(phase)) then
|
|
call math_invert33(Fe(1:3,1:3,g,ip,el), invFe, detFe, inversionError)
|
|
|
|
!* in case of periodic surfaces we have to find out how many periodic images in each direction we need
|
|
|
|
do dir = 1_pInt,3_pInt
|
|
maxCoord(dir) = maxval(mesh_node0(dir,:))
|
|
minCoord(dir) = minval(mesh_node0(dir,:))
|
|
enddo
|
|
meshSize = maxCoord - minCoord
|
|
coords = mesh_cellCenterCoordinates(ip,el)
|
|
periodicImages = 0_pInt
|
|
do dir = 1_pInt,3_pInt
|
|
if (mesh_periodicSurface(dir)) then
|
|
periodicImages(1,dir) = floor((coords(dir) - cutoffRadius(instance) - minCoord(dir)) / meshSize(dir), pInt)
|
|
periodicImages(2,dir) = ceiling((coords(dir) + cutoffRadius(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 neighbor_el = 1_pInt,mesh_NcpElems
|
|
ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el)))
|
|
neighbor_phase = material_phase(g,neighbor_ip,neighbor_el)
|
|
if (phase_localPlasticity(neighbor_phase)) then
|
|
cycle
|
|
endif
|
|
neighbor_instance = phase_plasticityInstance(neighbor_phase)
|
|
neighbor_latticeStruct = constitutive_nonlocal_structure(neighbor_instance)
|
|
neighbor_ns = totalNslip(neighbor_instance)
|
|
call math_invert33(Fe(1:3,1:3,1,neighbor_ip,neighbor_el), neighbor_invFe, detFe, inversionError)
|
|
neighbor_ipVolumeSideLength = mesh_ipVolume(neighbor_ip,neighbor_el) ** (1.0_pReal/3.0_pReal) ! reference volume used here
|
|
forall (s = 1_pInt:neighbor_ns, c = 1_pInt:2_pInt)
|
|
neighbor_rhoExcess(c,1,s) = state(g,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c-1,neighbor_instance)) & ! positive mobiles
|
|
- state(g,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c,neighbor_instance)) ! negative mobiles
|
|
neighbor_rhoExcess(c,2,s) = abs(state(g,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c-1,neighbor_instance))) & ! positive deads
|
|
- abs(state(g,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c,neighbor_instance))) ! negative deads
|
|
endforall
|
|
Tdislo_neighborLattice = 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 (neighbor_el /= el .or. neighbor_ip /= ip &
|
|
.or. deltaX /= 0_pInt .or. deltaY /= 0_pInt .or. deltaZ /= 0_pInt) then
|
|
|
|
neighbor_coords = mesh_cellCenterCoordinates(neighbor_ip,neighbor_el) &
|
|
+ (/real(deltaX,pReal), real(deltaY,pReal), real(deltaZ,pReal)/) * meshSize
|
|
connection = neighbor_coords - coords
|
|
distance = sqrt(sum(connection * connection))
|
|
if (distance > cutoffRadius(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 neighbor dislocation segment
|
|
|
|
connection_neighborLattice = math_mul33x3(neighbor_invFe, connection)
|
|
segmentLength = min(neighbor_ipVolumeSideLength, distance)
|
|
|
|
|
|
!* loop through all slip systems of the neighbor material point
|
|
!* and add up the stress contributions from egde and screw excess on these slip systems (if significant)
|
|
|
|
do s = 1_pInt,neighbor_ns
|
|
if (all(abs(neighbor_rhoExcess(:,:,s)) < significantRho(instance))) then
|
|
cycle ! not significant
|
|
endif
|
|
|
|
|
|
!* map the connection vector from the lattice into the slip system frame
|
|
|
|
connection_neighborSlip = math_mul33x3(lattice2slip(1:3,1:3,s,neighbor_instance), &
|
|
connection_neighborLattice)
|
|
|
|
|
|
!* edge contribution to stress
|
|
sigma = 0.0_pReal
|
|
|
|
x = connection_neighborSlip(1)
|
|
y = connection_neighborSlip(2)
|
|
z = connection_neighborSlip(3)
|
|
xsquare = x * x
|
|
ysquare = y * y
|
|
zsquare = z * z
|
|
|
|
do j = 1_pInt,2_pInt
|
|
if (abs(neighbor_rhoExcess(1,j,s)) < significantRho(instance)) then
|
|
cycle
|
|
elseif (j > 1_pInt) then
|
|
x = connection_neighborSlip(1) + sign(0.5_pReal * segmentLength, &
|
|
state(g,neighbor_ip,neighbor_el)%p(iRhoB(s,1,neighbor_instance)) &
|
|
- state(g,neighbor_ip,neighbor_el)%p(iRhoB(s,2,neighbor_instance)))
|
|
xsquare = x * x
|
|
endif
|
|
|
|
flipSign = sign(1.0_pReal, -y)
|
|
do side = 1_pInt,-1_pInt,-2_pInt
|
|
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
|
|
exit ipLoop
|
|
endif
|
|
|
|
sigma(1,1) = sigma(1,1) - real(side,pReal) &
|
|
* flipSign * z / denominator &
|
|
* (1.0_pReal + xsquare / Rsquare + xsquare / denominator) &
|
|
* neighbor_rhoExcess(1,j,s)
|
|
sigma(2,2) = sigma(2,2) - real(side,pReal) &
|
|
* (flipSign * 2.0_pReal * nu(instance) * z / denominator + z * lambda / Rcube) &
|
|
* neighbor_rhoExcess(1,j,s)
|
|
sigma(3,3) = sigma(3,3) + real(side,pReal) &
|
|
* flipSign * z / denominator &
|
|
* (1.0_pReal - zsquare / Rsquare - zsquare / denominator) &
|
|
* neighbor_rhoExcess(1,j,s)
|
|
sigma(1,2) = sigma(1,2) + real(side,pReal) &
|
|
* x * z / Rcube * neighbor_rhoExcess(1,j,s)
|
|
sigma(1,3) = sigma(1,3) + real(side,pReal) &
|
|
* flipSign * x / denominator &
|
|
* (1.0_pReal - zsquare / Rsquare - zsquare / denominator) &
|
|
* neighbor_rhoExcess(1,j,s)
|
|
sigma(2,3) = sigma(2,3) - real(side,pReal) &
|
|
* (nu(instance) / R - zsquare / Rcube) * neighbor_rhoExcess(1,j,s)
|
|
enddo
|
|
enddo
|
|
|
|
!* screw contribution to stress
|
|
|
|
x = connection_neighborSlip(1) ! have to restore this value, because position might have been adapted for edge deads before
|
|
do j = 1_pInt,2_pInt
|
|
if (abs(neighbor_rhoExcess(2,j,s)) < significantRho(instance)) then
|
|
cycle
|
|
elseif (j > 1_pInt) then
|
|
y = connection_neighborSlip(2) + sign(0.5_pReal * segmentLength, &
|
|
state(g,neighbor_ip,neighbor_el)%p(iRhoB(s,3,neighbor_instance)) &
|
|
- state(g,neighbor_ip,neighbor_el)%p(iRhoB(s,4,neighbor_instance)))
|
|
ysquare = y * y
|
|
endif
|
|
|
|
flipSign = sign(1.0_pReal, x)
|
|
do side = 1_pInt,-1_pInt,-2_pInt
|
|
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
|
|
exit ipLoop
|
|
endif
|
|
|
|
sigma(1,2) = sigma(1,2) - real(side,pReal) * flipSign * z * (1.0_pReal - nu(instance)) / denominator &
|
|
* neighbor_rhoExcess(2,j,s)
|
|
sigma(1,3) = sigma(1,3) + real(side,pReal) * flipSign * y * (1.0_pReal - nu(instance)) / denominator &
|
|
* neighbor_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 neighbor material point's lattice configuration
|
|
|
|
sigma = sigma * mu(neighbor_instance) * burgers(s,neighbor_instance) &
|
|
/ (4.0_pReal * pi * (1.0_pReal - nu(neighbor_instance))) &
|
|
* mesh_ipVolume(neighbor_ip,neighbor_el) / segmentLength ! reference volume is used here (according to the segment length calculation)
|
|
Tdislo_neighborLattice = Tdislo_neighborLattice &
|
|
+ math_mul33x33(math_transpose33(lattice2slip(1:3,1:3,s,neighbor_instance)), &
|
|
math_mul33x33(sigma, lattice2slip(1:3,1:3,s,neighbor_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_pInt:ns, c = 1_pInt:2_pInt) &
|
|
rhoExcessDead(c,s) = state(g,ip,el)%p(iRhoB(s,2*c-1,instance)) & ! 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(iRhoB(s,2*c,instance)) ! 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_pInt,ns
|
|
if (all(abs(rhoExcessDead(:,s)) < significantRho(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(instance))) &
|
|
* neighbor_ipVolumeSideLength * mu(instance) * burgers(s,instance) &
|
|
/ (sqrt(2.0_pReal) * pi * (1.0_pReal - nu(instance)))
|
|
sigma(3,1) = sigma(1,3)
|
|
|
|
Tdislo_neighborLattice = Tdislo_neighborLattice &
|
|
+ math_mul33x33(math_transpose33(lattice2slip(1:3,1:3,s,instance)), &
|
|
math_mul33x33(sigma, 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 neighbor MP's lattice configuration into the deformed configuration
|
|
!* and back into my lattice configuration
|
|
|
|
neighborLattice2myLattice = math_mul33x33(invFe, Fe(1:3,1:3,1,neighbor_ip,neighbor_el))
|
|
constitutive_nonlocal_dislocationstress = constitutive_nonlocal_dislocationstress &
|
|
+ math_mul33x33(neighborLattice2myLattice, &
|
|
math_mul33x33(Tdislo_neighborLattice, &
|
|
math_transpose33(neighborLattice2myLattice)))
|
|
|
|
enddo ipLoop
|
|
enddo ! element loop
|
|
|
|
endif
|
|
|
|
end function constitutive_nonlocal_dislocationstress
|
|
|
|
|
|
!*********************************************************************
|
|
!* return array of constitutive results *
|
|
!*********************************************************************
|
|
function constitutive_nonlocal_postResults(Tstar_v, Fe, Temperature, dt, state, dotState, g,ip,el)
|
|
|
|
use math, only: math_mul6x6, &
|
|
math_mul33x3, &
|
|
math_mul33x33, &
|
|
pi
|
|
use mesh, only: mesh_NcpElems, &
|
|
mesh_maxNips, &
|
|
mesh_ipVolume
|
|
use material, only: homogenization_maxNgrains, &
|
|
material_phase, &
|
|
phase_plasticityInstance, &
|
|
phase_Noutput
|
|
use lattice, only: lattice_Sslip_v, &
|
|
lattice_sd, &
|
|
lattice_st, &
|
|
lattice_sn
|
|
|
|
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,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
|
|
Fe ! elastic deformation gradient
|
|
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
|
|
state ! current microstructural state
|
|
type(p_vec), intent(in) :: dotState ! evolution rate of microstructural state
|
|
|
|
!*** output variables
|
|
real(pReal), dimension(constitutive_nonlocal_sizePostResults(phase_plasticityInstance(material_phase(g,ip,el)))) :: &
|
|
constitutive_nonlocal_postResults
|
|
|
|
!*** local variables
|
|
integer(pInt) matID, & ! current instance of this plasticity
|
|
structID, & ! 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(totalNslip(phase_plasticityInstance(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(totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),4) :: &
|
|
gdot, & ! shear rates
|
|
v ! velocities
|
|
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(g,ip,el)))) :: &
|
|
rhoForest, & ! forest dislocation density
|
|
tauThreshold, & ! threshold shear stress
|
|
tau, & ! current resolved shear stress
|
|
tauBack ! back stress from pileups on same slip system
|
|
real(pReal), dimension(totalNslip(phase_plasticityInstance(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,totalNslip(phase_plasticityInstance(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(3,totalNslip(phase_plasticityInstance(material_phase(g,ip,el)))) :: &
|
|
n_currentconf ! slip system normal (unit vector) in current configuration
|
|
real(pReal), dimension(3,3) :: sigma
|
|
|
|
matID = phase_plasticityInstance(material_phase(g,ip,el))
|
|
structID = constitutive_nonlocal_structure(matID)
|
|
ns = totalNslip(matID)
|
|
|
|
cs = 0_pInt
|
|
constitutive_nonlocal_postResults = 0.0_pReal
|
|
|
|
|
|
!* short hand notations for state variables
|
|
|
|
forall (s = 1_pInt:ns, t = 1_pInt:4_pInt)
|
|
rhoSgl(s,t) = state(g,ip,el)%p(iRhoU(s,t,matID))
|
|
rhoSgl(s,t+4_pInt) = state(g,ip,el)%p(iRhoB(s,t,matID))
|
|
v(s,t) = state(g,ip,el)%p(iV(s,t,matID))
|
|
rhoDotSgl(s,t) = dotState%p(iRhoU(s,t,matID))
|
|
rhoDotSgl(s,t+4_pInt) = dotState%p(iRhoB(s,t,matID))
|
|
endforall
|
|
forall (s = 1_pInt:ns, c = 1_pInt:2_pInt)
|
|
rhoDip(s,c) = state(g,ip,el)%p(iRhoD(s,c,matID))
|
|
rhoDotDip(s,c) = dotState%p(iRhoD(s,c,matID))
|
|
endforall
|
|
rhoForest = state(g,ip,el)%p(iRhoF(1:ns,matID))
|
|
tauThreshold = state(g,ip,el)%p(iTauF(1:ns,matID))
|
|
tauBack = state(g,ip,el)%p(iTauB(1:ns,matID))
|
|
|
|
|
|
|
|
!* Calculate shear rate
|
|
|
|
forall (t = 1_pInt:4_pInt) &
|
|
gdot(1:ns,t) = rhoSgl(1:ns,t) * burgers(1:ns,matID) * v(1:ns,t)
|
|
|
|
|
|
!* calculate limits for stable dipole height
|
|
|
|
do s = 1_pInt,ns
|
|
sLattice = slipSystemLattice(s,matID)
|
|
tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,structID)) + tauBack(s)
|
|
if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal
|
|
enddo
|
|
|
|
dLower = minDipoleHeight(1:ns,1:2,matID)
|
|
dUpper(1:ns,1) = mu(matID) * burgers(1:ns,matID) &
|
|
/ (8.0_pReal * pi * (1.0_pReal - nu(matID)) * abs(tau))
|
|
dUpper(1:ns,2) = mu(matID) * burgers(1:ns,matID) &
|
|
/ (4.0_pReal * pi * abs(tau))
|
|
forall (c = 1_pInt:2_pInt) &
|
|
dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) &
|
|
+ abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), &
|
|
dUpper(1:ns,c))
|
|
dUpper = max(dUpper,dLower)
|
|
|
|
|
|
!*** dislocation motion
|
|
|
|
m(1:3,1:ns,1) = lattice_sd(1:3,slipSystemLattice(1:ns,matID),structID)
|
|
m(1:3,1:ns,2) = -lattice_st(1:3,slipSystemLattice(1:ns,matID),structID)
|
|
forall (c = 1_pInt:2_pInt, s = 1_pInt:ns) &
|
|
m_currentconf(1:3,s,c) = math_mul33x3(Fe(1:3,1:3,g,ip,el), m(1:3,s,c))
|
|
forall (s = 1_pInt:ns) &
|
|
n_currentconf(1:3,s) = math_mul33x3(Fe(1:3,1:3,g,ip,el), &
|
|
lattice_sn(1:3,slipSystemLattice(s,matID),structID))
|
|
|
|
|
|
do o = 1_pInt,phase_Noutput(material_phase(g,ip,el))
|
|
select case(constitutive_nonlocal_output(o,matID))
|
|
|
|
case ('rho')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(abs(rhoSgl),2) + sum(rhoDip,2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(abs(rhoSgl),2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_mobile')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(abs(rhoSgl(1:ns,1:4)),2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_immobile')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoSgl(1:ns,5:8),2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_dip')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDip,2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_edge')
|
|
constitutive_nonlocal_postResults(cs+1_pInt: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_pInt: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_pInt:cs+ns) = sum(rhoSgl(1:ns,1:2),2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_edge_immobile')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoSgl(1:ns,5:6),2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_edge_pos')
|
|
constitutive_nonlocal_postResults(cs+1_pInt: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_pInt:cs+ns) = rhoSgl(1:ns,1)
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_edge_pos_immobile')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,5)
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_edge_neg')
|
|
constitutive_nonlocal_postResults(cs+1_pInt: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_pInt:cs+ns) = rhoSgl(1:ns,2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_edge_neg_immobile')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,6)
|
|
cs = cs + ns
|
|
|
|
case ('rho_dip_edge')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDip(1:ns,1)
|
|
cs = cs + ns
|
|
|
|
case ('rho_screw')
|
|
constitutive_nonlocal_postResults(cs+1_pInt: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_pInt: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_pInt:cs+ns) = sum(rhoSgl(1:ns,3:4),2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_screw_immobile')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoSgl(1:ns,7:8),2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_screw_pos')
|
|
constitutive_nonlocal_postResults(cs+1_pInt: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_pInt:cs+ns) = rhoSgl(1:ns,3)
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_screw_pos_immobile')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,7)
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_screw_neg')
|
|
constitutive_nonlocal_postResults(cs+1_pInt: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_pInt:cs+ns) = rhoSgl(1:ns,4)
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_screw_neg_immobile')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoSgl(1:ns,8)
|
|
cs = cs + ns
|
|
|
|
case ('rho_dip_screw')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDip(1:ns,2)
|
|
cs = cs + ns
|
|
|
|
case ('excess_rho')
|
|
constitutive_nonlocal_postResults(cs+1_pInt: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_pInt: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_pInt: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_pInt:cs+ns) = rhoForest
|
|
cs = cs + ns
|
|
|
|
case ('delta')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = 1.0_pReal / sqrt(sum(abs(rhoSgl),2) + sum(rhoDip,2))
|
|
cs = cs + ns
|
|
|
|
case ('delta_sgl')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = 1.0_pReal / sqrt(sum(abs(rhoSgl),2))
|
|
cs = cs + ns
|
|
|
|
case ('delta_dip')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = 1.0_pReal / sqrt(sum(rhoDip,2))
|
|
cs = cs + ns
|
|
|
|
case ('shearrate')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(gdot,2)
|
|
cs = cs + ns
|
|
|
|
case ('resolvedstress')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = tau
|
|
cs = cs + ns
|
|
|
|
case ('resolvedstress_back')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = tauBack
|
|
cs = cs + ns
|
|
|
|
case ('resolvedstress_external')
|
|
do s = 1_pInt,ns
|
|
sLattice = slipSystemLattice(s,matID)
|
|
constitutive_nonlocal_postResults(cs+s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,structID))
|
|
enddo
|
|
cs = cs + ns
|
|
|
|
case ('resistance')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = tauThreshold
|
|
cs = cs + ns
|
|
|
|
case ('rho_dot')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotSgl,2) + sum(rhoDotDip,2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_dot_sgl')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotSgl,2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_dot_dip')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotDip,2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_dot_gen')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotMultiplicationOutput(1:ns,1,g,ip,el) &
|
|
+ rhoDotMultiplicationOutput(1:ns,2,g,ip,el)
|
|
cs = cs + ns
|
|
|
|
case ('rho_dot_gen_edge')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotMultiplicationOutput(1:ns,1,g,ip,el)
|
|
cs = cs + ns
|
|
|
|
case ('rho_dot_gen_screw')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotMultiplicationOutput(1:ns,2,g,ip,el)
|
|
cs = cs + ns
|
|
|
|
case ('rho_dot_sgl2dip')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotSingle2DipoleGlideOutput(1:ns,1,g,ip,el) &
|
|
+ rhoDotSingle2DipoleGlideOutput(1:ns,2,g,ip,el)
|
|
cs = cs + ns
|
|
|
|
case ('rho_dot_sgl2dip_edge')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotSingle2DipoleGlideOutput(1:ns,1,g,ip,el)
|
|
cs = cs + ns
|
|
|
|
case ('rho_dot_sgl2dip_screw')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotSingle2DipoleGlideOutput(1:ns,2,g,ip,el)
|
|
cs = cs + ns
|
|
|
|
case ('rho_dot_ann_ath')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotAthermalAnnihilationOutput(1:ns,1,g,ip,el) &
|
|
+ rhoDotAthermalAnnihilationOutput(1:ns,2,g,ip,el)
|
|
cs = cs + ns
|
|
|
|
case ('rho_dot_ann_the')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotThermalAnnihilationOutput(1:ns,1,g,ip,el) &
|
|
+ rhoDotThermalAnnihilationOutput(1:ns,2,g,ip,el)
|
|
cs = cs + ns
|
|
|
|
case ('rho_dot_ann_the_edge')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotThermalAnnihilationOutput(1:ns,1,g,ip,el)
|
|
cs = cs + ns
|
|
|
|
case ('rho_dot_ann_the_screw')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotThermalAnnihilationOutput(1:ns,2,g,ip,el)
|
|
cs = cs + ns
|
|
|
|
case ('rho_dot_edgejogs')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotEdgeJogsOutput(1:ns,g,ip,el)
|
|
cs = cs + ns
|
|
|
|
case ('rho_dot_flux')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotFluxOutput(1:ns,1:4,g,ip,el),2) &
|
|
+ sum(abs(rhoDotFluxOutput(1:ns,5:8,g,ip,el)),2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_dot_flux_edge')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotFluxOutput(1:ns,1:2,g,ip,el),2) &
|
|
+ sum(abs(rhoDotFluxOutput(1:ns,5:6,g,ip,el)),2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_dot_flux_screw')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotFluxOutput(1:ns,3:4,g,ip,el),2) &
|
|
+ sum(abs(rhoDotFluxOutput(1:ns,7:8,g,ip,el)),2)
|
|
cs = cs + ns
|
|
|
|
case ('velocity_edge_pos')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = v(1:ns,1)
|
|
cs = cs + ns
|
|
|
|
case ('velocity_edge_neg')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = v(1:ns,2)
|
|
cs = cs + ns
|
|
|
|
case ('velocity_screw_pos')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = v(1:ns,3)
|
|
cs = cs + ns
|
|
|
|
case ('velocity_screw_neg')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = v(1:ns,4)
|
|
cs = cs + ns
|
|
|
|
case ('slipdirection.x')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = m_currentconf(1,1:ns,1)
|
|
cs = cs + ns
|
|
|
|
case ('slipdirection.y')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = m_currentconf(2,1:ns,1)
|
|
cs = cs + ns
|
|
|
|
case ('slipdirection.z')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = m_currentconf(3,1:ns,1)
|
|
cs = cs + ns
|
|
|
|
case ('slipnormal.x')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = n_currentconf(1,1:ns)
|
|
cs = cs + ns
|
|
|
|
case ('slipnormal.y')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = n_currentconf(2,1:ns)
|
|
cs = cs + ns
|
|
|
|
case ('slipnormal.z')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = n_currentconf(3,1:ns)
|
|
cs = cs + ns
|
|
|
|
case ('fluxdensity_edge_pos.x')
|
|
constitutive_nonlocal_postResults(cs+1_pInt: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_pInt: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_pInt: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_pInt: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_pInt: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_pInt: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_pInt: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_pInt: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_pInt: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_pInt: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_pInt: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_pInt:cs+ns) = - rhoSgl(1:ns,4) * v(1:ns,4) * m_currentconf(3,1:ns,2)
|
|
cs = cs + ns
|
|
|
|
case ('maximumdipoleheight_edge')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = dUpper(1:ns,1)
|
|
cs = cs + ns
|
|
|
|
case ('maximumdipoleheight_screw')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = dUpper(1:ns,2)
|
|
cs = cs + ns
|
|
|
|
case('dislocationstress')
|
|
sigma = constitutive_nonlocal_dislocationstress(state, Fe, g, ip, el)
|
|
constitutive_nonlocal_postResults(cs+1_pInt) = sigma(1,1)
|
|
constitutive_nonlocal_postResults(cs+2_pInt) = sigma(2,2)
|
|
constitutive_nonlocal_postResults(cs+3_pInt) = sigma(3,3)
|
|
constitutive_nonlocal_postResults(cs+4_pInt) = sigma(1,2)
|
|
constitutive_nonlocal_postResults(cs+5_pInt) = sigma(2,3)
|
|
constitutive_nonlocal_postResults(cs+6_pInt) = sigma(3,1)
|
|
cs = cs + 6_pInt
|
|
|
|
case('accumulatedshear')
|
|
constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = state(g,ip,el)%p(iGamma(1:ns,matID))
|
|
cs = cs + ns
|
|
|
|
case('boundarylayer')
|
|
do s = 1_pInt,ns
|
|
if (sum(abs(rhoSgl(s,1:8))) > 0.0_pReal) then
|
|
constitutive_nonlocal_postResults(cs+s) = maxval(abs(rhoSgl(s,5:8))/(rhoSgl(s,1:4)+abs(rhoSgl(s,5:8))))
|
|
else
|
|
constitutive_nonlocal_postResults(cs+s) = 0.0_pReal
|
|
endif
|
|
enddo
|
|
cs = cs + ns
|
|
|
|
end select
|
|
enddo
|
|
|
|
endfunction
|
|
|
|
END MODULE
|