! 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 . ! !############################################################## !* $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, & maxNmatIDs, & 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" maxNmatIDs = int(count(phase_plasticity == CONSTITUTIVE_NONLOCAL_LABEL),pInt) if (maxNmatIDs == 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:',maxNmatIDs !*** memory allocation for global variables allocate(constitutive_nonlocal_sizeDotState(maxNmatIDs)) allocate(constitutive_nonlocal_sizeDependentState(maxNmatIDs)) allocate(constitutive_nonlocal_sizeState(maxNmatIDs)) allocate(constitutive_nonlocal_sizePostResults(maxNmatIDs)) allocate(constitutive_nonlocal_sizePostResult(maxval(phase_Noutput), maxNmatIDs)) allocate(constitutive_nonlocal_output(maxval(phase_Noutput), maxNmatIDs)) allocate(Noutput(maxNmatIDs)) 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(maxNmatIDs)) allocate(constitutive_nonlocal_structure(maxNmatIDs)) allocate(Nslip(lattice_maxNslipFamily, maxNmatIDs)) allocate(slipFamily(lattice_maxNslip, maxNmatIDs)) allocate(slipSystemLattice(lattice_maxNslip, maxNmatIDs)) allocate(totalNslip(maxNmatIDs)) constitutive_nonlocal_structureName = '' constitutive_nonlocal_structure = 0_pInt Nslip = 0_pInt slipFamily = 0_pInt slipSystemLattice = 0_pInt totalNslip = 0_pInt allocate(CoverA(maxNmatIDs)) allocate(mu(maxNmatIDs)) allocate(nu(maxNmatIDs)) allocate(atomicVolume(maxNmatIDs)) allocate(Dsd0(maxNmatIDs)) allocate(selfDiffusionEnergy(maxNmatIDs)) allocate(aTolRho(maxNmatIDs)) allocate(aTolShear(maxNmatIDs)) allocate(significantRho(maxNmatIDs)) allocate(significantN(maxNmatIDs)) allocate(Cslip66(6,6,maxNmatIDs)) allocate(Cslip3333(3,3,3,3,maxNmatIDs)) allocate(cutoffRadius(maxNmatIDs)) allocate(doublekinkwidth(maxNmatIDs)) allocate(solidSolutionEnergy(maxNmatIDs)) allocate(solidSolutionSize(maxNmatIDs)) allocate(solidSolutionConcentration(maxNmatIDs)) allocate(pParam(maxNmatIDs)) allocate(qParam(maxNmatIDs)) allocate(viscosity(maxNmatIDs)) allocate(fattack(maxNmatIDs)) allocate(rhoSglScatter(maxNmatIDs)) allocate(rhoSglRandom(maxNmatIDs)) allocate(rhoSglRandomBinning(maxNmatIDs)) allocate(surfaceTransmissivity(maxNmatIDs)) allocate(grainboundaryTransmissivity(maxNmatIDs)) allocate(shortRangeStressCorrection(maxNmatIDs)) allocate(probabilisticMultiplication(maxNmatIDs)) allocate(CFLfactor(maxNmatIDs)) allocate(fEdgeMultiplication(maxNmatIDs)) allocate(linetensionEffect(maxNmatIDs)) allocate(edgeJogFactor(maxNmatIDs)) 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,maxNmatIDs)) allocate(rhoSglEdgeNeg0(lattice_maxNslipFamily,maxNmatIDs)) allocate(rhoSglScrewPos0(lattice_maxNslipFamily,maxNmatIDs)) allocate(rhoSglScrewNeg0(lattice_maxNslipFamily,maxNmatIDs)) allocate(rhoDipEdge0(lattice_maxNslipFamily,maxNmatIDs)) allocate(rhoDipScrew0(lattice_maxNslipFamily,maxNmatIDs)) allocate(burgersPerSlipFamily(lattice_maxNslipFamily,maxNmatIDs)) allocate(lambda0PerSlipFamily(lattice_maxNslipFamily,maxNmatIDs)) allocate(interactionSlipSlip(lattice_maxNinteraction,maxNmatIDs)) 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,maxNmatIDs)) allocate(peierlsStressPerSlipFamily(lattice_maxNslipFamily,2,maxNmatIDs)) minDipoleHeightPerSlipFamily = -1.0_pReal peierlsStressPerSlipFamily = 0.0_pReal allocate(nonSchmidCoeff(lattice_maxNnonSchmid,maxNmatIDs)) 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 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,maxNmatIDs 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,maxNmatIDs)) allocate(iRhoB(maxTotalNslip,4,maxNmatIDs)) allocate(iRhoD(maxTotalNslip,2,maxNmatIDs)) allocate(iV(maxTotalNslip,4,maxNmatIDs)) allocate(iD(maxTotalNslip,2,maxNmatIDs)) allocate(iGamma(maxTotalNslip,maxNmatIDs)) allocate(iRhoF(maxTotalNslip,maxNmatIDs)) allocate(iTauF(maxTotalNslip,maxNmatIDs)) allocate(iTauB(maxTotalNslip,maxNmatIDs)) 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,maxNmatIDs)) burgers = 0.0_pReal allocate(lambda0(maxTotalNslip,maxNmatIDs)) lambda0 = 0.0_pReal allocate(minDipoleHeight(maxTotalNslip,2,maxNmatIDs)) minDipoleHeight = -1.0_pReal allocate(forestProjectionEdge(maxTotalNslip,maxTotalNslip,maxNmatIDs)) forestProjectionEdge = 0.0_pReal allocate(forestProjectionScrew(maxTotalNslip,maxTotalNslip,maxNmatIDs)) forestProjectionScrew = 0.0_pReal allocate(interactionMatrixSlipSlip(maxTotalNslip,maxTotalNslip,maxNmatIDs)) interactionMatrixSlipSlip = 0.0_pReal allocate(lattice2slip(1:3, 1:3, maxTotalNslip, maxNmatIDs)) 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,maxNmatIDs)) peierlsStress = 0.0_pReal allocate(colinearSystem(maxTotalNslip,maxNmatIDs)) colinearSystem = 0_pInt allocate(nonSchmidProjection(3,3,4,maxTotalNslip,maxNmatIDs)) nonSchmidProjection = 0.0_pReal do i = 1,maxNmatIDs 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 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, & 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, & maxNmatIDs real(pReal), dimension(2) :: noise real(pReal), dimension(4) :: rnd real(pReal) meanDensity, & totalVolume, & densityBinning, & minimumIpVolume maxNmatIDs = 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,maxNmatIDs 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(ipc,ip,el) use mesh, only: & mesh_NcpElems, & mesh_maxNips use material, only: & homogenization_maxNgrains, & material_phase, & phase_plasticityInstance implicit none integer(pInt), intent(in) :: & ipc, & ! current grain ID ip, & ! current integration point el ! current element real(pReal), dimension(6,6) :: & constitutive_nonlocal_homogenizedC constitutive_nonlocal_homogenizedC = & Cslip66(1:6,1:6,phase_plasticityInstance(material_phase(ipc,ip,el))) end function constitutive_nonlocal_homogenizedC !-------------------------------------------------------------------------------------------------- !> @brief calculates quantities characterizing the microstructure !-------------------------------------------------------------------------------------------------- subroutine constitutive_nonlocal_microstructure(state, 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 implicit none !*** input variables integer(pInt), intent(in) :: gr, & ! current grain ID ip, & ! current integration point el ! current element 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) neighbor_el, & ! element number of neighboring material point neighbor_ip, & ! integration point of neighboring material point matID, & ! my instance of this plasticity neighbor_matID, & ! instance of this plasticity of neighboring material point structID, & ! my lattice structure neighbor_structID, & ! lattice structure of neighboring material point phase, & neighbor_phaseID, & ns, & ! total number of active slip systems at my material point neighbor_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) :: neighbors 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) :: & neighbor_rhoExcess, & ! excess density at neighboring material point neighbor_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) matID = phase_plasticityInstance(phase) structID = constitutive_nonlocal_structure(matID) ns = totalNslip(matID) !*** 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,matID)), 0.0_pReal) ! ensure positive single mobile densities rhoSgl(s,t+4_pInt) = state(gr,ip,el)%p(iRhoB(s,t,matID)) endforall forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) & rhoDip(s,c) = max(state(gr,ip,el)%p(iRhoD(s,c,matID)), 0.0_pReal) ! ensure positive dipole densities 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 !*** 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,matID)) & + dot_product((sum(abs(rhoSgl(1:ns,[3,4,7,8])),2) + rhoDip(1:ns,2)), & forestProjectionScrew(s,1:ns,matID)) !*** 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,matID) if (structID < 3_pInt) then ! only fcc and bcc do s = 1_pInt,ns myRhoForest = max(rhoForest(s),significantRho(matID)) correction = ( 1.0_pReal - linetensionEffect(matID) & + linetensionEffect(matID) & * log(0.35_pReal * burgers(s,matID) * sqrt(myRhoForest)) & / log(0.35_pReal * burgers(s,matID) * 1e6_pReal)) ** 2.0_pReal myInteractionMatrix(s,1:ns) = correction * myInteractionMatrix(s,1:ns) enddo endif forall (s = 1_pInt:ns) & tauThreshold(s) = mu(matID) * burgers(s,matID) & * 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(matID)) 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 neighbor_rhoTotal = 0.0_pReal do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el)))) neighbor_el = mesh_ipNeighborhood(1,n,ip,el) neighbor_ip = mesh_ipNeighborhood(2,n,ip,el) if (neighbor_el > 0 .and. neighbor_ip > 0) then neighbor_phaseID = material_phase(gr,neighbor_ip,neighbor_el) neighbor_matID = phase_plasticityInstance(neighbor_phaseID) neighbor_structID = constitutive_nonlocal_structure(neighbor_matID) neighbor_ns = totalNslip(neighbor_matID) if (.not. phase_localPlasticity(neighbor_phaseID) & .and. neighbor_structID == structID & .and. neighbor_matID == matID) then if (neighbor_ns == ns) then nRealNeighbors = nRealNeighbors + 1_pInt forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) neighbor_rhoExcess(c,s,n) = & max(state(gr,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c-1,neighbor_matID)), 0.0_pReal) &! positive mobiles - max(state(gr,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c,neighbor_matID)), 0.0_pReal) ! negative mobiles neighbor_rhoTotal(c,s,n) = & max(state(gr,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c-1,neighbor_matID)), 0.0_pReal) &! positive mobiles + max(state(gr,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c,neighbor_matID)), 0.0_pReal) & ! negative mobiles + abs(state(gr,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c-1,neighbor_matID))) & ! positive deads + abs(state(gr,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c,neighbor_matID))) & ! negative deads + max(state(gr,neighbor_ip,neighbor_el)%p(iRhoD(s,c,neighbor_matID)), 0.0_pReal) ! dipoles endforall connection_latticeConf(1:3,n) = & math_mul33x3(invFe, mesh_ipCoordinates(1:3,neighbor_ip,neighbor_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 ! neighboring 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 neighbor_rhoExcess(1:2,1:ns,n) = rhoExcess endif else ! free surface -> use central values instead connection_latticeConf(1:3,n) = 0.0_pReal neighbor_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,matID),structID) m(1:3,1:ns,2) = -lattice_st(1:3,slipSystemLattice(1:ns,matID),structID) do s = 1_pInt,ns !* gradient from interpolation of neighboring excess density do c = 1_pInt,2_pInt do dir = 1_pInt,3_pInt neighbors(1) = 2_pInt * dir - 1_pInt neighbors(2) = 2_pInt * dir connections(dir,1:3) = connection_latticeConf(1:3,neighbors(1)) & - connection_latticeConf(1:3,neighbors(2)) rhoExcessDifferences(dir) = neighbor_rhoExcess(c,s,neighbors(1)) & - neighbor_rhoExcess(c,s,neighbors(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(neighbor_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(matID) * burgers(s,matID) / (2.0_pReal * pi) & * (rhoExcessGradient_over_rho(1) / (1.0_pReal - nu(matID)) + rhoExcessGradient_over_rho(2)) enddo endif !*** set dependent states state(gr,ip,el)%p(iRhoF(1:ns,matID)) = rhoForest state(gr,ip,el)%p(iTauF(1:ns,matID)) = tauThreshold state(gr,ip,el)%p(iTauB(1:ns,matID)) = 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, ipc, 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) :: ipc, & !< 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(ipc,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(ipc,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) :: matID, & !< 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 matID = phase_plasticityInstance(material_phase(ipc,ip,el)) ns = totalNslip(matID) 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,matID) jumpWidth_P = burgers(s,matID) activationLength_P = doublekinkwidth(matID) * burgers(s,matID) activationVolume_P = activationLength_P * jumpWidth_P * burgers(s,matID) criticalStress_P = peierlsStress(s,c,matID) 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(matID) & * exp(activationEnergy_P / (KB * Temperature) & * (1.0_pReal - tauRel_P**pParam(matID))**qParam(matID)) if (tauEff < criticalStress_P) then dtPeierls_dtau = tPeierls * pParam(matID) * qParam(matID) * activationVolume_P / (KB * Temperature) & * (1.0_pReal - tauRel_P**pParam(matID))**(qParam(matID)-1.0_pReal) & * tauRel_P**(pParam(matID)-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,matID) / sqrt(solidSolutionConcentration(matID)) jumpWidth_S = solidSolutionSize(matID) * burgers(s,matID) activationLength_S = burgers(s,matID) / sqrt(solidSolutionConcentration(matID)) activationVolume_S = activationLength_S * jumpWidth_S * burgers(s,matID) activationEnergy_S = solidSolutionEnergy(matID) 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(matID) & * exp(activationEnergy_S / (KB * Temperature) & * (1.0_pReal - tauRel_S**pParam(matID))**qParam(matID)) if (tauEff < criticalStress_S) then dtSolidSolution_dtau = tSolidSolution * pParam(matID) * qParam(matID) & * activationVolume_S / (KB * Temperature) & * (1.0_pReal - tauRel_S**pParam(matID))**(qParam(matID)-1.0_pReal) & * tauRel_S**(pParam(matID)-1.0_pReal) else dtSolidSolution_dtau = 0.0_pReal endif !* viscous glide velocity tauEff = abs(tau(s)) - tauThreshold(s) mobility = burgers(s,matID) / viscosity(matID) 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 == ipc)& .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 ipc',el,ip,ipc 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, ipc, 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: 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) :: ipc, & !< 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(ipc,ip,el))),8) :: & rhoSgl !< single dislocation densities (including blocked) real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,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(ipc,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(ipc,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, ipc, 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, ipc, 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 == ipc)& .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 ipc ',el,ip,ipc 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, ipc,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) :: ipc, & ! current grain number ip, & ! current integration point el ! current element number 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(ipc,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(ipc,ip,el))),8) :: & rhoSgl ! current single dislocation densities (positive/negative screw and edge without dipoles) real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),4) :: & v ! dislocation glide velocity real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,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(ipc,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 == ipc)& .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 ipc ',el,ip,ipc write(6,*) endif #endif matID = phase_plasticityInstance(material_phase(ipc,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(ipc,ip,el)%p(iRhoU(s,t,matID)), 0.0_pReal) ! ensure positive single mobile densities rhoSgl(s,t+4_pInt) = state(ipc,ip,el)%p(iRhoB(s,t,matID)) v(s,t) = state(ipc,ip,el)%p(iV(s,t,matID)) endforall forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) rhoDip(s,c) = max(state(ipc,ip,el)%p(iRhoD(s,c,matID)), 0.0_pReal) ! ensure positive dipole densities dUpperOld(s,c) = state(ipc,ip,el)%p(iD(s,c,matID)) endforall tauBack = state(ipc,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(ipc,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 == ipc)& .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, ipc,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) :: ipc, & !< 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(ipc,ip,el)))) :: & constitutive_nonlocal_dotState !< evolution of state variables / microstructure !*** local variables integer(pInt) matID, & !< current instance of this plasticity neighbor_matID, & !< 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(ipc,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(ipc,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(ipc,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(ipc,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(ipc,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(ipc,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 == ipc)& .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 ipc ',el,ip,ipc write(6,*) endif #endif matID = phase_plasticityInstance(material_phase(ipc,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(ipc,ip,el)%p(iRhoU(s,t,matID)), 0.0_pReal) ! ensure positive single mobile densities rhoSgl(s,t+4_pInt) = state(ipc,ip,el)%p(iRhoB(s,t,matID)) v(s,t) = state(ipc,ip,el)%p(iV(s,t,matID)) endforall forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) rhoDip(s,c) = max(state(ipc,ip,el)%p(iRhoD(s,c,matID)), 0.0_pReal) ! ensure positive dipole densities endforall rhoForest = state(ipc,ip,el)%p(iRhoF(1:ns,matID)) tauThreshold = state(ipc,ip,el)%p(iTauF(1:ns,matID)) tauBack = state(ipc,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(ipc,ip,el)%p(iRhoU(s,t,matID)), 0.0_pReal) rhoSgl0(s,t+4_pInt) = state0(ipc,ip,el)%p(iRhoB(s,t,matID)) v0(s,t) = state0(ipc,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 == ipc)& .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,ipc,ip,el) > 1.0_pReal) then call random_number(rnd) sourceProbability(s,ipc,ip,el) = rnd !$OMP FLUSH(sourceProbability) endif if (sourceProbability(s,ipc,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,ipc,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 == ipc)& .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(ipc,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,ipc,ip,el) my_F = math_mul33x33(my_Fe, Fp(1:3,1:3,ipc,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_matID = phase_plasticityInstance(material_phase(ipc,neighbor_ip,neighbor_el)) neighbor_Fe = Fe(1:3,1:3,ipc,neighbor_ip,neighbor_el) neighbor_F = math_mul33x33(neighbor_Fe, Fp(1:3,1:3,ipc,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(ipc,neighbor_ip,neighbor_el) /= subfrac(ipc,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(ipc,neighbor_ip,neighbor_el)%p(iV(s,t,neighbor_matID)) neighbor_rhoSgl(s,t) = max(state0(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,t,neighbor_matID)), 0.0_pReal) neighbor_rhoSgl(s,t+4_pInt) = state0(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,t,neighbor_matID)) endforall else forall (s = 1:ns, t = 1_pInt:4_pInt) neighbor_v(s,t) = state(ipc,neighbor_ip,neighbor_el)%p(iV(s,t,neighbor_matID)) neighbor_rhoSgl(s,t) = max(state(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,t,neighbor_matID)), 0.0_pReal) neighbor_rhoSgl(s,t+4_pInt) = state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,t,neighbor_matID)) 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(ipc,ip,el) == 0.0_pReal) then my_rhoSgl = rhoSgl0 my_v = v0 elseif (neighbor_n > 0_pInt) then if (subfrac(ipc,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,ipc,ip,el) = rhoDotFlux(1:ns,1:8) rhoDotMultiplicationOutput(1:ns,1:2,ipc,ip,el) = rhoDotMultiplication(1:ns,[1,3]) rhoDotSingle2DipoleGlideOutput(1:ns,1:2,ipc,ip,el) = rhoDotSingle2DipoleGlide(1:ns,9:10) rhoDotAthermalAnnihilationOutput(1:ns,1:2,ipc,ip,el) = rhoDotAthermalAnnihilation(1:ns,9:10) rhoDotThermalAnnihilationOutput(1:ns,1:2,ipc,ip,el) = rhoDotThermalAnnihilation(1:ns,9:10) rhoDotEdgeJogsOutput(1:ns,ipc,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 == ipc)& .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 phaseID, & neighbor_phaseID, & textureID, & neighbor_textureID, & structID, & ! lattice structure matID, & ! 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)))) phaseID = material_phase(1,i,e) textureID = material_texture(1,i,e) matID = phase_plasticityInstance(phaseID) structID = constitutive_nonlocal_structure(matID) ns = totalNslip(matID) slipNormal(1:3,1:ns) = lattice_sn(1:3, slipSystemLattice(1:ns,matID), structID) slipDirection(1:3,1:ns) = lattice_sd(1:3, slipSystemLattice(1:ns,matID), structID) !*** 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(matID)) 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_phaseID = material_phase(1,neighbor_i,neighbor_e) if (neighbor_phaseID /= phaseID) then if (.not. phase_localPlasticity(neighbor_phaseID) .and. .not. phase_localPlasticity(phaseID)) 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(matID) >= 0.0_pReal) then neighbor_textureID = material_texture(1,neighbor_i,neighbor_e) if (neighbor_textureID /= textureID) then if (.not. phase_localPlasticity(neighbor_phaseID)) then forall(s1 = 1_pInt:ns) & my_compatibility(1:2,s1,s1,n) = sqrt(grainboundaryTransmissivity(matID)) 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 * !********************************************************************* pure function constitutive_nonlocal_dislocationstress(state, Fe, ipc, 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) :: ipc, & ! 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 matID, & ! my instance of this plasticity neighbor_matID, & ! instance of this plasticity of neighbor material point structID, & ! my lattice structure neighbor_structID, & ! lattice structure of neighbor material point phase, & neighbor_phaseID, & 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(ipc,ip,el))),8) :: & rhoSgl ! single dislocation density (edge+, edge-, screw+, screw-, used edge+, used edge-, used screw+, used screw-) logical inversionError phase = material_phase(ipc,ip,el) matID = phase_plasticityInstance(phase) structID = constitutive_nonlocal_structure(matID) ns = totalNslip(matID) !*** get basic states forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) rhoSgl(s,t) = max(state(ipc,ip,el)%p(iRhoU(s,t,matID)), 0.0_pReal) ! ensure positive single mobile densities rhoSgl(s,t+4_pInt) = state(ipc,ip,el)%p(iRhoB(s,t,matID)) 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,ipc,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(matID) - minCoord(dir)) / meshSize(dir), pInt) periodicImages(2,dir) = ceiling((coords(dir) + cutoffRadius(matID) - 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_phaseID = material_phase(ipc,neighbor_ip,neighbor_el) if (phase_localPlasticity(neighbor_phaseID)) then cycle endif neighbor_matID = phase_plasticityInstance(neighbor_phaseID) neighbor_structID = constitutive_nonlocal_structure(neighbor_matID) neighbor_ns = totalNslip(neighbor_matID) 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(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c-1,neighbor_matID)) & ! positive mobiles - state(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c,neighbor_matID)) ! negative mobiles neighbor_rhoExcess(c,2,s) = abs(state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c-1,neighbor_matID))) & ! positive deads - abs(state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c,neighbor_matID))) ! 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(matID)) 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(matID))) 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_matID), & 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(matID)) then cycle elseif (j > 1_pInt) then x = connection_neighborSlip(1) + sign(0.5_pReal * segmentLength, & state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,1,neighbor_matID)) & - state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,2,neighbor_matID))) 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(matID) * 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(matID) / 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(matID)) then cycle elseif (j > 1_pInt) then y = connection_neighborSlip(2) + sign(0.5_pReal * segmentLength, & state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,3,neighbor_matID)) & - state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,4,neighbor_matID))) 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(matID)) / denominator & * neighbor_rhoExcess(2,j,s) sigma(1,3) = sigma(1,3) + real(side,pReal) * flipSign * y * (1.0_pReal - nu(matID)) / 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_matID) * burgers(s,neighbor_matID) & / (4.0_pReal * pi * (1.0_pReal - nu(neighbor_matID))) & * 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_matID)), & math_mul33x33(sigma, lattice2slip(1:3,1:3,s,neighbor_matID))) 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(ipc,ip,el)%p(iRhoB(s,2*c-1,matID)) & ! 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(ipc,ip,el)%p(iRhoB(s,2*c,matID)) ! 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(matID))) 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(matID))) & * neighbor_ipVolumeSideLength * mu(matID) * burgers(s,matID) & / (sqrt(2.0_pReal) * pi * (1.0_pReal - nu(matID))) sigma(3,1) = sigma(1,3) Tdislo_neighborLattice = Tdislo_neighborLattice & + math_mul33x33(math_transpose33(lattice2slip(1:3,1:3,s,matID)), & math_mul33x33(sigma, lattice2slip(1:3,1:3,s,matID))) 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 !-------------------------------------------------------------------------------------------------- !> @brief return array of constitutive results !-------------------------------------------------------------------------------------------------- pure function constitutive_nonlocal_postResults(Tstar_v,Fe,state,dotState,ipc,ip,el) use math, only: & math_mul6x6, & math_mul33x3, & math_mul33x33, & pi use mesh, only: & mesh_NcpElems, & mesh_maxNips 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 real(pReal), dimension(6), intent(in) :: & Tstar_v !< 2nd Piola Kirchhoff stress tensor 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 !< microstructure state type(p_vec), intent(in) :: dotState ! evolution rate of microstructural state integer(pInt), intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point el !< element real(pReal), dimension(constitutive_nonlocal_sizePostResults(& phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & constitutive_nonlocal_postResults 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(ipc,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(ipc,ip,el))),4) :: & gdot, & !< shear rates v !< velocities real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,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(ipc,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(ipc,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(ipc,ip,el)))) :: & n_currentconf !< slip system normal (unit vector) in current configuration real(pReal), dimension(3,3) :: & sigma matID = phase_plasticityInstance(material_phase(ipc,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(ipc,ip,el)%p(iRhoU(s,t,matID)) rhoSgl(s,t+4_pInt) = state(ipc,ip,el)%p(iRhoB(s,t,matID)) v(s,t) = state(ipc,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(ipc,ip,el)%p(iRhoD(s,c,matID)) rhoDotDip(s,c) = dotState%p(iRhoD(s,c,matID)) endforall rhoForest = state(ipc,ip,el)%p(iRhoF(1:ns,matID)) tauThreshold = state(ipc,ip,el)%p(iTauF(1:ns,matID)) tauBack = state(ipc,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,ipc,ip,el), m(1:3,s,c)) forall (s = 1_pInt:ns) & n_currentconf(1:3,s) = math_mul33x3(Fe(1:3,1:3,ipc,ip,el), & lattice_sn(1:3,slipSystemLattice(s,matID),structID)) do o = 1_pInt,phase_Noutput(material_phase(ipc,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,ipc,ip,el) & + rhoDotMultiplicationOutput(1:ns,2,ipc,ip,el) cs = cs + ns case ('rho_dot_gen_edge') constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotMultiplicationOutput(1:ns,1,ipc,ip,el) cs = cs + ns case ('rho_dot_gen_screw') constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotMultiplicationOutput(1:ns,2,ipc,ip,el) cs = cs + ns case ('rho_dot_sgl2dip') constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotSingle2DipoleGlideOutput(1:ns,1,ipc,ip,el) & + rhoDotSingle2DipoleGlideOutput(1:ns,2,ipc,ip,el) cs = cs + ns case ('rho_dot_sgl2dip_edge') constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotSingle2DipoleGlideOutput(1:ns,1,ipc,ip,el) cs = cs + ns case ('rho_dot_sgl2dip_screw') constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotSingle2DipoleGlideOutput(1:ns,2,ipc,ip,el) cs = cs + ns case ('rho_dot_ann_ath') constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotAthermalAnnihilationOutput(1:ns,1,ipc,ip,el) & + rhoDotAthermalAnnihilationOutput(1:ns,2,ipc,ip,el) cs = cs + ns case ('rho_dot_ann_the') constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotThermalAnnihilationOutput(1:ns,1,ipc,ip,el) & + rhoDotThermalAnnihilationOutput(1:ns,2,ipc,ip,el) cs = cs + ns case ('rho_dot_ann_the_edge') constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotThermalAnnihilationOutput(1:ns,1,ipc,ip,el) cs = cs + ns case ('rho_dot_ann_the_screw') constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotThermalAnnihilationOutput(1:ns,2,ipc,ip,el) cs = cs + ns case ('rho_dot_edgejogs') constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotEdgeJogsOutput(1:ns,ipc,ip,el) cs = cs + ns case ('rho_dot_flux') constitutive_nonlocal_postResults(cs+1_pInt:cs+ns) = sum(rhoDotFluxOutput(1:ns,1:4,ipc,ip,el),2) & + sum(abs(rhoDotFluxOutput(1:ns,5:8,ipc,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,ipc,ip,el),2) & + sum(abs(rhoDotFluxOutput(1:ns,5:6,ipc,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,ipc,ip,el),2) & + sum(abs(rhoDotFluxOutput(1:ns,7:8,ipc,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, ipc, 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(ipc,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 end function constitutive_nonlocal_postResults end module constitutive_nonlocal