2540 lines
147 KiB
Fortran
2540 lines
147 KiB
Fortran
!* $Id$
|
|
!************************************
|
|
!* Module: CONSTITUTIVE_NONLOCAL *
|
|
!************************************
|
|
!* contains: *
|
|
!* - constitutive equations *
|
|
!* - parameters definition *
|
|
!************************************
|
|
|
|
|
|
MODULE constitutive_nonlocal
|
|
|
|
!* Include other modules
|
|
use prec, only: pReal,pInt
|
|
implicit none
|
|
|
|
|
|
!* Definition of parameters
|
|
character (len=*), parameter :: constitutive_nonlocal_label = 'nonlocal'
|
|
character(len=22), dimension(10), parameter :: constitutive_nonlocal_listBasicStates = (/'rhoSglEdgePosMobile ', &
|
|
'rhoSglEdgeNegMobile ', &
|
|
'rhoSglScrewPosMobile ', &
|
|
'rhoSglScrewNegMobile ', &
|
|
'rhoSglEdgePosImmobile ', &
|
|
'rhoSglEdgeNegImmobile ', &
|
|
'rhoSglScrewPosImmobile', &
|
|
'rhoSglScrewNegImmobile', &
|
|
'rhoDipEdge ', &
|
|
'rhoDipScrew ' /) ! list of "basic" microstructural state variables that are independent from other state variables
|
|
character(len=15), dimension(3), parameter :: constitutive_nonlocal_listDependentStates = (/'rhoForest ', &
|
|
'tauThreshold ', &
|
|
'Tdislocation_v ' /) ! list of microstructural state variables that depend on other state variables
|
|
real(pReal), parameter :: kB = 1.38e-23_pReal ! Physical parameter, Boltzmann constant in J/Kelvin
|
|
|
|
!* Definition of global variables
|
|
integer(pInt), dimension(:), allocatable :: constitutive_nonlocal_sizeDotState, & ! number of dotStates
|
|
constitutive_nonlocal_sizeState, & ! total number of microstructural state variables
|
|
constitutive_nonlocal_sizePostResults ! cumulative size of post results
|
|
integer(pInt), dimension(:,:), allocatable, target :: constitutive_nonlocal_sizePostResult ! size of each post result output
|
|
character(len=64), dimension(:,:), allocatable, target :: constitutive_nonlocal_output ! name of each post result output
|
|
|
|
character(len=32), dimension(:), allocatable :: constitutive_nonlocal_structureName ! name of the lattice structure
|
|
integer(pInt), dimension(:), allocatable :: constitutive_nonlocal_structure, & ! number representing the kind of lattice structure
|
|
constitutive_nonlocal_totalNslip ! total number of active slip systems for each instance
|
|
integer(pInt), dimension(:,:), allocatable :: constitutive_nonlocal_Nslip, & ! number of active slip systems for each family and instance
|
|
constitutive_nonlocal_slipFamily, & ! lookup table relating active slip system to slip family for each instance
|
|
constitutive_nonlocal_slipSystemLattice ! lookup table relating active slip system index to lattice slip system index for each instance
|
|
|
|
real(pReal), dimension(:), allocatable :: constitutive_nonlocal_CoverA, & ! c/a ratio for hex type lattice
|
|
constitutive_nonlocal_C11, & ! C11 element in elasticity matrix
|
|
constitutive_nonlocal_C12, & ! C12 element in elasticity matrix
|
|
constitutive_nonlocal_C13, & ! C13 element in elasticity matrix
|
|
constitutive_nonlocal_C33, & ! C33 element in elasticity matrix
|
|
constitutive_nonlocal_C44, & ! C44 element in elasticity matrix
|
|
constitutive_nonlocal_Gmod, & ! shear modulus
|
|
constitutive_nonlocal_nu, & ! poisson's ratio
|
|
constitutive_nonlocal_atomicVolume, & ! atomic volume
|
|
constitutive_nonlocal_Dsd0, & ! prefactor for self-diffusion coefficient
|
|
constitutive_nonlocal_Qsd, & ! activation enthalpy for diffusion
|
|
constitutive_nonlocal_aTolRho, & ! absolute tolerance for dislocation density in state integration
|
|
constitutive_nonlocal_R, & ! cutoff radius for dislocation stress
|
|
constitutive_nonlocal_d0, & ! wall depth as multiple of b
|
|
constitutive_nonlocal_tauObs, & ! obstacle strength in Pa
|
|
constitutive_nonlocal_fattack, & ! attack frequency in Hz
|
|
constitutive_nonlocal_vs, & ! maximum dislocation velocity = velocity of sound
|
|
constitutive_nonlocal_rhoSglScatter, & ! standard deviation of scatter in initial dislocation density
|
|
constitutive_nonlocal_surfaceTransmissivity ! transmissivity at free surface
|
|
real(pReal), dimension(:,:,:), allocatable :: constitutive_nonlocal_Cslip_66 ! elasticity matrix in Mandel notation for each instance
|
|
real(pReal), dimension(:,:,:,:,:), allocatable :: constitutive_nonlocal_Cslip_3333 ! elasticity matrix for each instance
|
|
real(pReal), dimension(:,:), allocatable :: constitutive_nonlocal_rhoSglEdgePos0, & ! initial edge_pos dislocation density per slip system for each family and instance
|
|
constitutive_nonlocal_rhoSglEdgeNeg0, & ! initial edge_neg dislocation density per slip system for each family and instance
|
|
constitutive_nonlocal_rhoSglScrewPos0, & ! initial screw_pos dislocation density per slip system for each family and instance
|
|
constitutive_nonlocal_rhoSglScrewNeg0, & ! initial screw_neg dislocation density per slip system for each family and instance
|
|
constitutive_nonlocal_rhoDipEdge0, & ! initial edge dipole dislocation density per slip system for each family and instance
|
|
constitutive_nonlocal_rhoDipScrew0, & ! initial screw dipole dislocation density per slip system for each family and instance
|
|
constitutive_nonlocal_lambda0PerSlipFamily, & ! mean free path prefactor for each family and instance
|
|
constitutive_nonlocal_lambda0PerSlipSystem, & ! mean free path prefactor for each slip system and instance
|
|
constitutive_nonlocal_burgersPerSlipFamily, & ! absolute length of burgers vector [m] for each family and instance
|
|
constitutive_nonlocal_burgersPerSlipSystem, & ! absolute length of burgers vector [m] for each slip system and instance
|
|
constitutive_nonlocal_dLowerEdgePerSlipFamily, & ! minimum stable edge dipole height for each family and instance
|
|
constitutive_nonlocal_dLowerEdgePerSlipSystem, & ! minimum stable edge dipole height for each slip system and instance
|
|
constitutive_nonlocal_dLowerScrewPerSlipFamily, & ! minimum stable screw dipole height for each family and instance
|
|
constitutive_nonlocal_dLowerScrewPerSlipSystem, & ! minimum stable screw dipole height for each slip system and instance
|
|
constitutive_nonlocal_Qeff0, & ! prefactor for activation enthalpy for dislocation glide in J
|
|
constitutive_nonlocal_interactionSlipSlip ! coefficients for slip-slip interaction for each interaction type and instance
|
|
real(pReal), dimension(:,:,:,:,:), allocatable :: constitutive_nonlocal_v, & ! dislocation velocity
|
|
constitutive_nonlocal_rhoDotFlux ! dislocation convection term
|
|
real(pReal), dimension(:,:,:,:,:,:), allocatable :: constitutive_nonlocal_compatibility ! slip system compatibility between me and my neighbors
|
|
real(pReal), dimension(:,:,:), allocatable :: constitutive_nonlocal_forestProjectionEdge, & ! matrix of forest projections of edge dislocations for each instance
|
|
constitutive_nonlocal_forestProjectionScrew, & ! matrix of forest projections of screw dislocations for each instance
|
|
constitutive_nonlocal_interactionMatrixSlipSlip ! interaction matrix of the different slip systems for each instance
|
|
|
|
CONTAINS
|
|
!****************************************
|
|
!* - constitutive_init
|
|
!* - constitutive_stateInit
|
|
!* - constitutive_homogenizedC
|
|
!* - constitutive_microstructure
|
|
!* - constitutive_LpAndItsTangent
|
|
!* - constitutive_dotState
|
|
!* - constitutive_dotTemperature
|
|
!* - constitutive_postResults
|
|
!****************************************
|
|
|
|
|
|
!**************************************
|
|
!* Module initialization *
|
|
!**************************************
|
|
subroutine constitutive_nonlocal_init(file)
|
|
|
|
use prec, only: pInt, pReal
|
|
use math, only: math_Mandel3333to66, &
|
|
math_Voigt66to3333, &
|
|
math_mul3x3
|
|
use IO, only: IO_lc, &
|
|
IO_getTag, &
|
|
IO_isBlank, &
|
|
IO_stringPos, &
|
|
IO_stringValue, &
|
|
IO_floatValue, &
|
|
IO_intValue, &
|
|
IO_error
|
|
use mesh, only: mesh_NcpElems, &
|
|
mesh_maxNips, &
|
|
FE_maxNipNeighbors
|
|
use material, only: homogenization_maxNgrains, &
|
|
phase_constitution, &
|
|
phase_constitutionInstance, &
|
|
phase_Noutput
|
|
use lattice, only: lattice_maxNslipFamily, &
|
|
lattice_maxNtwinFamily, &
|
|
lattice_maxNslip, &
|
|
lattice_maxNtwin, &
|
|
lattice_maxNinteraction, &
|
|
lattice_NslipSystem, &
|
|
lattice_NtwinSystem, &
|
|
lattice_initializeStructure, &
|
|
lattice_Qtwin, &
|
|
lattice_sd, &
|
|
lattice_sn, &
|
|
lattice_st, &
|
|
lattice_interactionSlipSlip
|
|
|
|
!*** output variables
|
|
|
|
!*** input variables
|
|
integer(pInt), intent(in) :: file
|
|
|
|
!*** local variables
|
|
integer(pInt), parameter :: maxNchunks = 21
|
|
integer(pInt), dimension(1+2*maxNchunks) :: positions
|
|
integer(pInt) section, &
|
|
maxNinstance, &
|
|
maxTotalNslip, &
|
|
myStructure, &
|
|
f, & ! index of my slip family
|
|
i, & ! index of my instance of this constitution
|
|
j, &
|
|
k, &
|
|
l, &
|
|
ns, & ! short notation for total number of active slip systems for the current instance
|
|
o, & ! index of my output
|
|
s, & ! index of my slip system
|
|
s1, & ! index of my slip system
|
|
s2, & ! index of my slip system
|
|
it, & ! index of my interaction type
|
|
output, &
|
|
mySize
|
|
character(len=64) tag
|
|
character(len=1024) line
|
|
|
|
|
|
write(6,*)
|
|
write(6,'(a20,a20,a12)') '<<<+- constitutive_',constitutive_nonlocal_label,' init -+>>>'
|
|
write(6,*) '$Id$'
|
|
write(6,*)
|
|
|
|
maxNinstance = count(phase_constitution == constitutive_nonlocal_label)
|
|
if (maxNinstance == 0) return ! we don't have to do anything if there's no instance for this constitutive law
|
|
|
|
|
|
!*** space allocation for global variables
|
|
|
|
allocate(constitutive_nonlocal_sizeDotState(maxNinstance))
|
|
allocate(constitutive_nonlocal_sizeState(maxNinstance))
|
|
allocate(constitutive_nonlocal_sizePostResults(maxNinstance))
|
|
allocate(constitutive_nonlocal_sizePostResult(maxval(phase_Noutput), maxNinstance))
|
|
allocate(constitutive_nonlocal_output(maxval(phase_Noutput), maxNinstance))
|
|
constitutive_nonlocal_sizeDotState = 0_pInt
|
|
constitutive_nonlocal_sizeState = 0_pInt
|
|
constitutive_nonlocal_sizePostResults = 0_pInt
|
|
constitutive_nonlocal_sizePostResult = 0_pInt
|
|
constitutive_nonlocal_output = ''
|
|
|
|
allocate(constitutive_nonlocal_structureName(maxNinstance))
|
|
allocate(constitutive_nonlocal_structure(maxNinstance))
|
|
allocate(constitutive_nonlocal_Nslip(lattice_maxNslipFamily, maxNinstance))
|
|
allocate(constitutive_nonlocal_slipFamily(lattice_maxNslip, maxNinstance))
|
|
allocate(constitutive_nonlocal_slipSystemLattice(lattice_maxNslip, maxNinstance))
|
|
allocate(constitutive_nonlocal_totalNslip(maxNinstance))
|
|
constitutive_nonlocal_structureName = ''
|
|
constitutive_nonlocal_structure = 0_pInt
|
|
constitutive_nonlocal_Nslip = 0_pInt
|
|
constitutive_nonlocal_slipFamily = 0_pInt
|
|
constitutive_nonlocal_slipSystemLattice = 0_pInt
|
|
constitutive_nonlocal_totalNslip = 0_pInt
|
|
|
|
allocate(constitutive_nonlocal_CoverA(maxNinstance))
|
|
allocate(constitutive_nonlocal_C11(maxNinstance))
|
|
allocate(constitutive_nonlocal_C12(maxNinstance))
|
|
allocate(constitutive_nonlocal_C13(maxNinstance))
|
|
allocate(constitutive_nonlocal_C33(maxNinstance))
|
|
allocate(constitutive_nonlocal_C44(maxNinstance))
|
|
allocate(constitutive_nonlocal_Gmod(maxNinstance))
|
|
allocate(constitutive_nonlocal_nu(maxNinstance))
|
|
allocate(constitutive_nonlocal_atomicVolume(maxNinstance))
|
|
allocate(constitutive_nonlocal_Dsd0(maxNinstance))
|
|
allocate(constitutive_nonlocal_Qsd(maxNinstance))
|
|
allocate(constitutive_nonlocal_aTolRho(maxNinstance))
|
|
allocate(constitutive_nonlocal_Cslip_66(6,6,maxNinstance))
|
|
allocate(constitutive_nonlocal_Cslip_3333(3,3,3,3,maxNinstance))
|
|
allocate(constitutive_nonlocal_R(maxNinstance))
|
|
allocate(constitutive_nonlocal_d0(maxNinstance))
|
|
allocate(constitutive_nonlocal_tauObs(maxNinstance))
|
|
allocate(constitutive_nonlocal_vs(maxNinstance))
|
|
allocate(constitutive_nonlocal_fattack(maxNinstance))
|
|
allocate(constitutive_nonlocal_rhoSglScatter(maxNinstance))
|
|
allocate(constitutive_nonlocal_surfaceTransmissivity(maxNinstance))
|
|
constitutive_nonlocal_CoverA = 0.0_pReal
|
|
constitutive_nonlocal_C11 = 0.0_pReal
|
|
constitutive_nonlocal_C12 = 0.0_pReal
|
|
constitutive_nonlocal_C13 = 0.0_pReal
|
|
constitutive_nonlocal_C33 = 0.0_pReal
|
|
constitutive_nonlocal_C44 = 0.0_pReal
|
|
constitutive_nonlocal_Gmod = 0.0_pReal
|
|
constitutive_nonlocal_atomicVolume = 0.0_pReal
|
|
constitutive_nonlocal_Dsd0 = 0.0_pReal
|
|
constitutive_nonlocal_Qsd = 0.0_pReal
|
|
constitutive_nonlocal_aTolRho = 0.0_pReal
|
|
constitutive_nonlocal_nu = 0.0_pReal
|
|
constitutive_nonlocal_Cslip_66 = 0.0_pReal
|
|
constitutive_nonlocal_Cslip_3333 = 0.0_pReal
|
|
constitutive_nonlocal_R = -1.0_pReal
|
|
constitutive_nonlocal_d0 = 0.0_pReal
|
|
constitutive_nonlocal_tauObs = 0.0_pReal
|
|
constitutive_nonlocal_vs = 0.0_pReal
|
|
constitutive_nonlocal_fattack = 0.0_pReal
|
|
constitutive_nonlocal_rhoSglScatter = 0.0_pReal
|
|
constitutive_nonlocal_surfaceTransmissivity = 1.0_pReal
|
|
|
|
allocate(constitutive_nonlocal_rhoSglEdgePos0(lattice_maxNslipFamily, maxNinstance))
|
|
allocate(constitutive_nonlocal_rhoSglEdgeNeg0(lattice_maxNslipFamily, maxNinstance))
|
|
allocate(constitutive_nonlocal_rhoSglScrewPos0(lattice_maxNslipFamily, maxNinstance))
|
|
allocate(constitutive_nonlocal_rhoSglScrewNeg0(lattice_maxNslipFamily, maxNinstance))
|
|
allocate(constitutive_nonlocal_rhoDipEdge0(lattice_maxNslipFamily, maxNinstance))
|
|
allocate(constitutive_nonlocal_rhoDipScrew0(lattice_maxNslipFamily, maxNinstance))
|
|
allocate(constitutive_nonlocal_burgersPerSlipFamily(lattice_maxNslipFamily, maxNinstance))
|
|
allocate(constitutive_nonlocal_Lambda0PerSlipFamily(lattice_maxNslipFamily, maxNinstance))
|
|
allocate(constitutive_nonlocal_interactionSlipSlip(lattice_maxNinteraction, maxNinstance))
|
|
allocate(constitutive_nonlocal_dLowerEdgePerSlipFamily(lattice_maxNslipFamily, maxNinstance))
|
|
allocate(constitutive_nonlocal_dLowerScrewPerSlipFamily(lattice_maxNslipFamily, maxNinstance))
|
|
constitutive_nonlocal_rhoSglEdgePos0 = -1.0_pReal
|
|
constitutive_nonlocal_rhoSglEdgeNeg0 = -1.0_pReal
|
|
constitutive_nonlocal_rhoSglScrewPos0 = -1.0_pReal
|
|
constitutive_nonlocal_rhoSglScrewNeg0 = -1.0_pReal
|
|
constitutive_nonlocal_rhoDipEdge0 = -1.0_pReal
|
|
constitutive_nonlocal_rhoDipScrew0 = -1.0_pReal
|
|
constitutive_nonlocal_burgersPerSlipFamily = 0.0_pReal
|
|
constitutive_nonlocal_lambda0PerSlipFamily = 0.0_pReal
|
|
constitutive_nonlocal_interactionSlipSlip = 0.0_pReal
|
|
constitutive_nonlocal_dLowerEdgePerSlipFamily = 0.0_pReal
|
|
constitutive_nonlocal_dLowerScrewPerSlipFamily = 0.0_pReal
|
|
|
|
!*** readout data from material.config file
|
|
|
|
rewind(file)
|
|
line = ''
|
|
section = 0
|
|
|
|
do while (IO_lc(IO_getTag(line,'<','>')) /= 'phase') ! wind forward to <phase>
|
|
read(file,'(a1024)',END=100) line
|
|
enddo
|
|
|
|
do ! read thru sections of phase part
|
|
read(file,'(a1024)',END=100) line
|
|
if (IO_isBlank(line)) cycle ! skip empty lines
|
|
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
|
|
if (IO_getTag(line,'[',']') /= '') then ! next section
|
|
section = section + 1
|
|
output = 0 ! reset output counter
|
|
cycle
|
|
endif
|
|
if (section > 0 .and. phase_constitution(section) == constitutive_nonlocal_label) then ! one of my sections
|
|
i = phase_constitutionInstance(section) ! which instance of my constitution is present phase
|
|
positions = IO_stringPos(line,maxNchunks)
|
|
tag = IO_lc(IO_stringValue(line,positions,1)) ! extract key
|
|
select case(tag)
|
|
case('constitution','/nonlocal/')
|
|
cycle
|
|
case ('(output)')
|
|
output = output + 1
|
|
constitutive_nonlocal_output(output,i) = IO_lc(IO_stringValue(line,positions,2))
|
|
case ('lattice_structure')
|
|
constitutive_nonlocal_structureName(i) = IO_lc(IO_stringValue(line,positions,2))
|
|
case ('covera_ratio')
|
|
constitutive_nonlocal_CoverA(i) = IO_floatValue(line,positions,2)
|
|
case ('c11')
|
|
constitutive_nonlocal_C11(i) = IO_floatValue(line,positions,2)
|
|
case ('c12')
|
|
constitutive_nonlocal_C12(i) = IO_floatValue(line,positions,2)
|
|
case ('c13')
|
|
constitutive_nonlocal_C13(i) = IO_floatValue(line,positions,2)
|
|
case ('c33')
|
|
constitutive_nonlocal_C33(i) = IO_floatValue(line,positions,2)
|
|
case ('c44')
|
|
constitutive_nonlocal_C44(i) = IO_floatValue(line,positions,2)
|
|
case ('nslip')
|
|
forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_Nslip(f,i) = IO_intValue(line,positions,1+f)
|
|
case ('rhosgledgepos0')
|
|
forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_rhoSglEdgePos0(f,i) = IO_floatValue(line,positions,1+f)
|
|
case ('rhosgledgeneg0')
|
|
forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_rhoSglEdgeNeg0(f,i) = IO_floatValue(line,positions,1+f)
|
|
case ('rhosglscrewpos0')
|
|
forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_rhoSglScrewPos0(f,i) = IO_floatValue(line,positions,1+f)
|
|
case ('rhosglscrewneg0')
|
|
forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_rhoSglScrewNeg0(f,i) = IO_floatValue(line,positions,1+f)
|
|
case ('rhodipedge0')
|
|
forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_rhoDipEdge0(f,i) = IO_floatValue(line,positions,1+f)
|
|
case ('rhodipscrew0')
|
|
forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_rhoDipScrew0(f,i) = IO_floatValue(line,positions,1+f)
|
|
case ('lambda0')
|
|
forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_lambda0PerSlipFamily(f,i) = IO_floatValue(line,positions,1+f)
|
|
case ('burgers')
|
|
forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_burgersPerSlipFamily(f,i) = IO_floatValue(line,positions,1+f)
|
|
case('r')
|
|
constitutive_nonlocal_R(i) = IO_floatValue(line,positions,2)
|
|
case('ddipminedge')
|
|
forall (f = 1:lattice_maxNslipFamily) &
|
|
constitutive_nonlocal_dLowerEdgePerSlipFamily(f,i) = IO_floatValue(line,positions,1+f)
|
|
case('ddipminscrew')
|
|
forall (f = 1:lattice_maxNslipFamily) &
|
|
constitutive_nonlocal_dLowerScrewPerSlipFamily(f,i) = IO_floatValue(line,positions,1+f)
|
|
case('atomicvolume')
|
|
constitutive_nonlocal_atomicVolume(i) = IO_floatValue(line,positions,2)
|
|
case('dsd0')
|
|
constitutive_nonlocal_Dsd0(i) = IO_floatValue(line,positions,2)
|
|
case('qsd')
|
|
constitutive_nonlocal_Qsd(i) = IO_floatValue(line,positions,2)
|
|
case('atol_rho')
|
|
constitutive_nonlocal_aTolRho(i) = IO_floatValue(line,positions,2)
|
|
case ('interaction_slipslip')
|
|
forall (it = 1:lattice_maxNinteraction) constitutive_nonlocal_interactionSlipSlip(it,i) = IO_floatValue(line,positions,1+it)
|
|
case('d0')
|
|
constitutive_nonlocal_d0(i) = IO_floatValue(line,positions,2)
|
|
case('tauobs')
|
|
constitutive_nonlocal_tauObs(i) = IO_floatValue(line,positions,2)
|
|
case('vs')
|
|
constitutive_nonlocal_vs(i) = IO_floatValue(line,positions,2)
|
|
case('fattack')
|
|
constitutive_nonlocal_fattack(i) = IO_floatValue(line,positions,2)
|
|
case('rhosglscatter')
|
|
constitutive_nonlocal_rhoSglScatter(i) = IO_floatValue(line,positions,2)
|
|
case('surfacetransmissivity')
|
|
constitutive_nonlocal_surfaceTransmissivity(i) = IO_floatValue(line,positions,2)
|
|
case default
|
|
call IO_error(236,ext_msg=tag)
|
|
end select
|
|
endif
|
|
enddo
|
|
|
|
|
|
100 do i = 1,maxNinstance
|
|
|
|
constitutive_nonlocal_structure(i) = &
|
|
lattice_initializeStructure(constitutive_nonlocal_structureName(i), constitutive_nonlocal_CoverA(i)) ! our lattice structure is defined in the material.config file by the structureName (and the c/a ratio)
|
|
myStructure = constitutive_nonlocal_structure(i)
|
|
|
|
!*** sanity checks
|
|
|
|
if (myStructure < 1 .or. myStructure > 3) call IO_error(205)
|
|
if (sum(constitutive_nonlocal_Nslip(:,i)) <= 0_pInt) call IO_error(235,ext_msg='Nslip')
|
|
do o = 1,maxval(phase_Noutput)
|
|
if(len(constitutive_nonlocal_output(o,i)) > 64) call IO_error(666)
|
|
enddo
|
|
do f = 1,lattice_maxNslipFamily
|
|
if (constitutive_nonlocal_Nslip(f,i) > 0_pInt) then
|
|
if (constitutive_nonlocal_rhoSglEdgePos0(f,i) < 0.0_pReal) call IO_error(235,ext_msg='rhoSglEdgePos0')
|
|
if (constitutive_nonlocal_rhoSglEdgeNeg0(f,i) < 0.0_pReal) call IO_error(235,ext_msg='rhoSglEdgeNeg0')
|
|
if (constitutive_nonlocal_rhoSglScrewPos0(f,i) < 0.0_pReal) call IO_error(235,ext_msg='rhoSglScrewPos0')
|
|
if (constitutive_nonlocal_rhoSglScrewNeg0(f,i) < 0.0_pReal) call IO_error(235,ext_msg='rhoSglScrewNeg0')
|
|
if (constitutive_nonlocal_rhoDipEdge0(f,i) < 0.0_pReal) call IO_error(235,ext_msg='rhoDipEdge0')
|
|
if (constitutive_nonlocal_rhoDipScrew0(f,i) < 0.0_pReal) call IO_error(235,ext_msg='rhoDipScrew0')
|
|
if (constitutive_nonlocal_burgersPerSlipFamily(f,i) <= 0.0_pReal) call IO_error(235,ext_msg='burgers')
|
|
if (constitutive_nonlocal_lambda0PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(235,ext_msg='lambda0')
|
|
if (constitutive_nonlocal_dLowerEdgePerSlipFamily(f,i) <= 0.0_pReal) call IO_error(235,ext_msg='dDipMinEdge')
|
|
if (constitutive_nonlocal_dLowerScrewPerSlipFamily(f,i) <= 0.0_pReal) call IO_error(235,ext_msg='dDipMinScrew')
|
|
endif
|
|
enddo
|
|
if (any(constitutive_nonlocal_interactionSlipSlip(1:maxval(lattice_interactionSlipSlip(:,:,myStructure)),i) < 0.0_pReal)) &
|
|
call IO_error(235,ext_msg='interaction_SlipSlip')
|
|
if (constitutive_nonlocal_R(i) < 0.0_pReal) call IO_error(235,ext_msg='r')
|
|
if (constitutive_nonlocal_atomicVolume(i) <= 0.0_pReal) call IO_error(235,ext_msg='atomicVolume')
|
|
if (constitutive_nonlocal_Dsd0(i) <= 0.0_pReal) call IO_error(235,ext_msg='Dsd0')
|
|
if (constitutive_nonlocal_Qsd(i) <= 0.0_pReal) call IO_error(235,ext_msg='Qsd')
|
|
if (constitutive_nonlocal_aTolRho(i) <= 0.0_pReal) call IO_error(235,ext_msg='aTol_rho')
|
|
if (constitutive_nonlocal_d0(i) <= 0.0_pReal) call IO_error(235,ext_msg='d0')
|
|
if (constitutive_nonlocal_tauObs(i) <= 0.0_pReal) call IO_error(235,ext_msg='tauObs')
|
|
if (constitutive_nonlocal_vs(i) <= 0.0_pReal) call IO_error(235,ext_msg='vs')
|
|
if (constitutive_nonlocal_fattack(i) <= 0.0_pReal) call IO_error(235,ext_msg='fAttack')
|
|
if (constitutive_nonlocal_rhoSglScatter(i) < 0.0_pReal) call IO_error(235,ext_msg='rhoSglScatter')
|
|
if (constitutive_nonlocal_surfaceTransmissivity(i) < 0.0_pReal &
|
|
.or. constitutive_nonlocal_surfaceTransmissivity(i) > 1.0_pReal) call IO_error(235,ext_msg='surfaceTransmissivity')
|
|
|
|
|
|
!*** determine total number of active slip systems
|
|
|
|
constitutive_nonlocal_Nslip(1:lattice_maxNslipFamily,i) = min( lattice_NslipSystem(1:lattice_maxNslipFamily, myStructure), &
|
|
constitutive_nonlocal_Nslip(1:lattice_maxNslipFamily,i) ) ! we can't use more slip systems per family than specified in lattice
|
|
constitutive_nonlocal_totalNslip(i) = sum(constitutive_nonlocal_Nslip(1:lattice_maxNslipFamily,i))
|
|
|
|
enddo
|
|
|
|
|
|
!*** allocation of variables whose size depends on the total number of active slip systems
|
|
|
|
maxTotalNslip = maxval(constitutive_nonlocal_totalNslip)
|
|
|
|
allocate(constitutive_nonlocal_burgersPerSlipSystem(maxTotalNslip, maxNinstance))
|
|
constitutive_nonlocal_burgersPerSlipSystem = 0.0_pReal
|
|
|
|
allocate(constitutive_nonlocal_lambda0PerSlipSystem(maxTotalNslip, maxNinstance))
|
|
constitutive_nonlocal_lambda0PerSlipSystem = 0.0_pReal
|
|
|
|
allocate(constitutive_nonlocal_dLowerEdgePerSlipSystem(maxTotalNslip, maxNinstance))
|
|
constitutive_nonlocal_dLowerEdgePerSlipSystem = 0.0_pReal
|
|
|
|
allocate(constitutive_nonlocal_dLowerScrewPerSlipSystem(maxTotalNslip, maxNinstance))
|
|
constitutive_nonlocal_dLowerScrewPerSlipSystem = 0.0_pReal
|
|
|
|
allocate(constitutive_nonlocal_Qeff0(maxTotalNslip, maxNinstance))
|
|
constitutive_nonlocal_Qeff0 = 0.0_pReal
|
|
|
|
allocate(constitutive_nonlocal_forestProjectionEdge(maxTotalNslip, maxTotalNslip, maxNinstance))
|
|
constitutive_nonlocal_forestProjectionEdge = 0.0_pReal
|
|
|
|
allocate(constitutive_nonlocal_forestProjectionScrew(maxTotalNslip, maxTotalNslip, maxNinstance))
|
|
constitutive_nonlocal_forestProjectionScrew = 0.0_pReal
|
|
|
|
allocate(constitutive_nonlocal_interactionMatrixSlipSlip(maxTotalNslip, maxTotalNslip, maxNinstance))
|
|
constitutive_nonlocal_interactionMatrixSlipSlip = 0.0_pReal
|
|
|
|
allocate(constitutive_nonlocal_v(maxTotalNslip, 4, homogenization_maxNgrains, mesh_maxNips, mesh_NcpElems))
|
|
constitutive_nonlocal_v = 0.0_pReal
|
|
|
|
allocate(constitutive_nonlocal_rhoDotFlux(maxTotalNslip, 10, homogenization_maxNgrains, mesh_maxNips, mesh_NcpElems))
|
|
constitutive_nonlocal_rhoDotFlux = 0.0_pReal
|
|
|
|
allocate(constitutive_nonlocal_compatibility(2,maxTotalNslip, maxTotalNslip, FE_maxNipNeighbors, mesh_maxNips, mesh_NcpElems))
|
|
constitutive_nonlocal_compatibility = 0.0_pReal
|
|
|
|
do i = 1,maxNinstance
|
|
|
|
myStructure = constitutive_nonlocal_structure(i) ! lattice structure of this instance
|
|
|
|
!*** Inverse lookup of my slip system family and the slip system in lattice
|
|
|
|
l = 0_pInt
|
|
do f = 1,lattice_maxNslipFamily
|
|
do s = 1,constitutive_nonlocal_Nslip(f,i)
|
|
l = l + 1
|
|
constitutive_nonlocal_slipFamily(l,i) = f
|
|
constitutive_nonlocal_slipSystemLattice(l,i) = sum(lattice_NslipSystem(1:f-1, myStructure)) + s
|
|
enddo; enddo
|
|
|
|
|
|
!*** determine size of state array
|
|
|
|
ns = constitutive_nonlocal_totalNslip(i)
|
|
constitutive_nonlocal_sizeState(i) = size(constitutive_nonlocal_listBasicStates) * ns &
|
|
+ (size(constitutive_nonlocal_listDependentStates) - 1_pInt) * ns + 6_pInt
|
|
constitutive_nonlocal_sizeDotState(i) = size(constitutive_nonlocal_listBasicStates) * ns
|
|
|
|
|
|
!*** determine size of postResults array
|
|
|
|
do o = 1,maxval(phase_Noutput)
|
|
select case(constitutive_nonlocal_output(o,i))
|
|
case( 'rho', &
|
|
'delta', &
|
|
'rho_edge', &
|
|
'rho_screw', &
|
|
'rho_sgl', &
|
|
'delta_sgl', &
|
|
'rho_sgl_edge', &
|
|
'rho_sgl_edge_pos', &
|
|
'rho_sgl_edge_neg', &
|
|
'rho_sgl_screw', &
|
|
'rho_sgl_screw_pos', &
|
|
'rho_sgl_screw_neg', &
|
|
'rho_sgl_mobile', &
|
|
'rho_sgl_edge_mobile', &
|
|
'rho_sgl_edge_pos_mobile', &
|
|
'rho_sgl_edge_neg_mobile', &
|
|
'rho_sgl_screw_mobile', &
|
|
'rho_sgl_screw_pos_mobile', &
|
|
'rho_sgl_screw_neg_mobile', &
|
|
'rho_sgl_immobile', &
|
|
'rho_sgl_edge_immobile', &
|
|
'rho_sgl_edge_pos_immobile', &
|
|
'rho_sgl_edge_neg_immobile', &
|
|
'rho_sgl_screw_immobile', &
|
|
'rho_sgl_screw_pos_immobile', &
|
|
'rho_sgl_screw_neg_immobile', &
|
|
'rho_dip', &
|
|
'delta_dip', &
|
|
'rho_dip_edge', &
|
|
'rho_dip_screw', &
|
|
'excess_rho', &
|
|
'excess_rho_edge', &
|
|
'excess_rho_screw', &
|
|
'rho_forest', &
|
|
'shearrate', &
|
|
'resolvedstress', &
|
|
'resolvedstress_internal', &
|
|
'resolvedstress_external', &
|
|
'resistance', &
|
|
'rho_dot', &
|
|
'rho_dot_sgl', &
|
|
'rho_dot_dip', &
|
|
'rho_dot_gen', &
|
|
'rho_dot_gen_edge', &
|
|
'rho_dot_gen_screw', &
|
|
'rho_dot_sgl2dip', &
|
|
'rho_dot_dip2sgl', &
|
|
'rho_dot_ann_ath', &
|
|
'rho_dot_ann_the', &
|
|
'rho_dot_flux', &
|
|
'rho_dot_flux_edge', &
|
|
'rho_dot_flux_screw', &
|
|
'dislocationvelocity', &
|
|
'fluxdensity_edge_pos_x', &
|
|
'fluxdensity_edge_pos_y', &
|
|
'fluxdensity_edge_pos_z', &
|
|
'fluxdensity_edge_neg_x', &
|
|
'fluxdensity_edge_neg_y', &
|
|
'fluxdensity_edge_neg_z', &
|
|
'fluxdensity_screw_pos_x', &
|
|
'fluxdensity_screw_pos_y', &
|
|
'fluxdensity_screw_pos_z', &
|
|
'fluxdensity_screw_neg_x', &
|
|
'fluxdensity_screw_neg_y', &
|
|
'fluxdensity_screw_neg_z', &
|
|
'd_upper_edge', &
|
|
'd_upper_screw', &
|
|
'd_upper_dot_edge', &
|
|
'd_upper_dot_screw' )
|
|
mySize = constitutive_nonlocal_totalNslip(i)
|
|
case default
|
|
mySize = 0_pInt
|
|
end select
|
|
|
|
if (mySize > 0_pInt) then ! any meaningful output found
|
|
constitutive_nonlocal_sizePostResult(o,i) = mySize
|
|
constitutive_nonlocal_sizePostResults(i) = constitutive_nonlocal_sizePostResults(i) + mySize
|
|
endif
|
|
enddo
|
|
|
|
|
|
!*** elasticity matrix and shear modulus according to material.config
|
|
|
|
select case (myStructure)
|
|
case(1:2) ! cubic(s)
|
|
forall(k=1:3)
|
|
forall(j=1:3) constitutive_nonlocal_Cslip_66(k,j,i) = constitutive_nonlocal_C12(i)
|
|
constitutive_nonlocal_Cslip_66(k,k,i) = constitutive_nonlocal_C11(i)
|
|
constitutive_nonlocal_Cslip_66(k+3,k+3,i) = constitutive_nonlocal_C44(i)
|
|
end forall
|
|
case(3:) ! all hex
|
|
constitutive_nonlocal_Cslip_66(1,1,i) = constitutive_nonlocal_C11(i)
|
|
constitutive_nonlocal_Cslip_66(2,2,i) = constitutive_nonlocal_C11(i)
|
|
constitutive_nonlocal_Cslip_66(3,3,i) = constitutive_nonlocal_C33(i)
|
|
constitutive_nonlocal_Cslip_66(1,2,i) = constitutive_nonlocal_C12(i)
|
|
constitutive_nonlocal_Cslip_66(2,1,i) = constitutive_nonlocal_C12(i)
|
|
constitutive_nonlocal_Cslip_66(1,3,i) = constitutive_nonlocal_C13(i)
|
|
constitutive_nonlocal_Cslip_66(3,1,i) = constitutive_nonlocal_C13(i)
|
|
constitutive_nonlocal_Cslip_66(2,3,i) = constitutive_nonlocal_C13(i)
|
|
constitutive_nonlocal_Cslip_66(3,2,i) = constitutive_nonlocal_C13(i)
|
|
constitutive_nonlocal_Cslip_66(4,4,i) = constitutive_nonlocal_C44(i)
|
|
constitutive_nonlocal_Cslip_66(5,5,i) = constitutive_nonlocal_C44(i)
|
|
constitutive_nonlocal_Cslip_66(6,6,i) = 0.5_pReal*(constitutive_nonlocal_C11(i)- constitutive_nonlocal_C12(i))
|
|
end select
|
|
constitutive_nonlocal_Cslip_66(1:6,1:6,i) = math_Mandel3333to66(math_Voigt66to3333(constitutive_nonlocal_Cslip_66(1:6,1:6,i)))
|
|
constitutive_nonlocal_Cslip_3333(1:3,1:3,1:3,1:3,i) = math_Voigt66to3333(constitutive_nonlocal_Cslip_66(1:6,1:6,i))
|
|
|
|
constitutive_nonlocal_Gmod(i) = 0.2_pReal * ( constitutive_nonlocal_C11(i) - constitutive_nonlocal_C12(i) &
|
|
+ 3.0_pReal*constitutive_nonlocal_C44(i) ) ! (C11iso-C12iso)/2 with C11iso=(3*C11+2*C12+4*C44)/5 and C12iso=(C11+4*C12-2*C44)/5
|
|
constitutive_nonlocal_nu(i) = ( constitutive_nonlocal_C11(i) + 4.0_pReal*constitutive_nonlocal_C12(i) &
|
|
- 2.0_pReal*constitutive_nonlocal_C44(i) ) &
|
|
/ ( 4.0_pReal*constitutive_nonlocal_C11(i) + 6.0_pReal*constitutive_nonlocal_C12(i) &
|
|
+ 2.0_pReal*constitutive_nonlocal_C44(i) ) ! C12iso/(C11iso+C12iso) with C11iso=(3*C11+2*C12+4*C44)/5 and C12iso=(C11+4*C12-2*C44)/5
|
|
|
|
|
|
do s1 = 1,ns
|
|
|
|
f = constitutive_nonlocal_slipFamily(s1,i)
|
|
|
|
!*** burgers vector, mean free path prefactor and minimum dipole distance for each slip system
|
|
|
|
constitutive_nonlocal_burgersPerSlipSystem(s1,i) = constitutive_nonlocal_burgersPerSlipFamily(f,i)
|
|
constitutive_nonlocal_lambda0PerSlipSystem(s1,i) = constitutive_nonlocal_lambda0PerSlipFamily(f,i)
|
|
constitutive_nonlocal_dLowerEdgePerSlipSystem(s1,i) = constitutive_nonlocal_dLowerEdgePerSlipFamily(f,i)
|
|
constitutive_nonlocal_dLowerScrewPerSlipSystem(s1,i) = constitutive_nonlocal_dLowerScrewPerSlipFamily(f,i)
|
|
|
|
do s2 = 1,ns
|
|
|
|
!*** calculation of forest projections for edge and screw dislocations. s2 acts as forest to s1
|
|
|
|
constitutive_nonlocal_forestProjectionEdge(s1, s2, i) &
|
|
= abs(math_mul3x3(lattice_sn(1:3, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure), &
|
|
lattice_st(1:3, constitutive_nonlocal_slipSystemLattice(s2,i), myStructure))) ! forest projection of edge dislocations is the projection of (t = b x n) onto the slip normal of the respective slip plane
|
|
|
|
constitutive_nonlocal_forestProjectionScrew(s1, s2, i) &
|
|
= abs(math_mul3x3(lattice_sn(1:3, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure), &
|
|
lattice_sd(1:3, constitutive_nonlocal_slipSystemLattice(s2,i), myStructure))) ! forest projection of screw dislocations is the projection of b onto the slip normal of the respective splip plane
|
|
|
|
!*** calculation of interaction matrices
|
|
|
|
constitutive_nonlocal_interactionMatrixSlipSlip(s1, s2, i) &
|
|
= constitutive_nonlocal_interactionSlipSlip( lattice_interactionSlipSlip(constitutive_nonlocal_slipSystemLattice(s1,i), &
|
|
constitutive_nonlocal_slipSystemLattice(s2,i), &
|
|
myStructure), &
|
|
i )
|
|
|
|
enddo; enddo
|
|
|
|
!*** calculation of prefactor for activation enthalpy for dislocation glide
|
|
|
|
constitutive_nonlocal_Qeff0(1:ns,i) = constitutive_nonlocal_burgersPerSlipSystem(1:ns,i) ** 3.0_pReal &
|
|
* dsqrt(0.5_pReal * constitutive_nonlocal_d0(i) ** 3.0_pReal &
|
|
* constitutive_nonlocal_Gmod(i) * constitutive_nonlocal_tauObs(i))
|
|
|
|
enddo
|
|
|
|
endsubroutine
|
|
|
|
|
|
|
|
!*********************************************************************
|
|
!* initial microstructural state (just the "basic" states) *
|
|
!*********************************************************************
|
|
function constitutive_nonlocal_stateInit(myInstance)
|
|
|
|
use prec, only: pReal, &
|
|
pInt
|
|
use lattice, only: lattice_maxNslipFamily
|
|
use math, only: math_sampleGaussVar
|
|
|
|
implicit none
|
|
|
|
!*** input variables
|
|
integer(pInt), intent(in) :: myInstance ! number specifying the current instance of the constitution
|
|
|
|
!*** output variables
|
|
real(pReal), dimension(constitutive_nonlocal_sizeState(myInstance)) :: &
|
|
constitutive_nonlocal_stateInit
|
|
|
|
!*** local variables
|
|
real(pReal), dimension(constitutive_nonlocal_totalNslip(myInstance)) :: &
|
|
rhoSglEdgePos, & ! positive edge dislocation density
|
|
rhoSglEdgeNeg, & ! negative edge dislocation density
|
|
rhoSglScrewPos, & ! positive screw dislocation density
|
|
rhoSglScrewNeg, & ! negative screw dislocation density
|
|
rhoSglEdgePosUsed, & ! used positive edge dislocation density
|
|
rhoSglEdgeNegUsed, & ! used negative edge dislocation density
|
|
rhoSglScrewPosUsed, & ! used positive screw dislocation density
|
|
rhoSglScrewNegUsed, & ! used negative screw dislocation density
|
|
rhoDipEdge, & ! edge dipole dislocation density
|
|
rhoDipScrew, & ! screw dipole dislocation density
|
|
rhoForest, & ! forest dislocation density
|
|
tauSlipThreshold ! threshold shear stress for slip
|
|
integer(pInt) ns, & ! short notation for total number of active slip systems
|
|
f, & ! index of lattice family
|
|
from, &
|
|
upto, &
|
|
s, & ! index of slip system
|
|
i
|
|
real(pReal), dimension(2) :: noise
|
|
|
|
constitutive_nonlocal_stateInit = 0.0_pReal
|
|
ns = constitutive_nonlocal_totalNslip(myInstance)
|
|
|
|
!*** set the basic state variables
|
|
|
|
do f = 1,lattice_maxNslipFamily
|
|
from = 1 + sum(constitutive_nonlocal_Nslip(1:f-1,myInstance))
|
|
upto = sum(constitutive_nonlocal_Nslip(1:f,myInstance))
|
|
do s = from,upto
|
|
do i = 1,2
|
|
noise(i) = math_sampleGaussVar(0.0_pReal, constitutive_nonlocal_rhoSglScatter(myInstance))
|
|
enddo
|
|
rhoSglEdgePos(s) = constitutive_nonlocal_rhoSglEdgePos0(f, myInstance) + noise(1)
|
|
rhoSglEdgeNeg(s) = constitutive_nonlocal_rhoSglEdgeNeg0(f, myInstance) + noise(1)
|
|
rhoSglScrewPos(s) = constitutive_nonlocal_rhoSglScrewPos0(f, myInstance) + noise(2)
|
|
rhoSglScrewNeg(s) = constitutive_nonlocal_rhoSglScrewNeg0(f, myInstance) + noise(2)
|
|
enddo
|
|
rhoSglEdgePosUsed(from:upto) = 0.0_pReal
|
|
rhoSglEdgeNegUsed(from:upto) = 0.0_pReal
|
|
rhoSglScrewPosUsed(from:upto) = 0.0_pReal
|
|
rhoSglScrewNegUsed(from:upto) = 0.0_pReal
|
|
rhoDipEdge(from:upto) = constitutive_nonlocal_rhoDipEdge0(f, myInstance)
|
|
rhoDipScrew(from:upto) = constitutive_nonlocal_rhoDipScrew0(f, myInstance)
|
|
enddo
|
|
|
|
|
|
!*** put everything together and in right order
|
|
|
|
constitutive_nonlocal_stateInit( 1: ns) = rhoSglEdgePos
|
|
constitutive_nonlocal_stateInit( ns+1: 2*ns) = rhoSglEdgeNeg
|
|
constitutive_nonlocal_stateInit( 2*ns+1: 3*ns) = rhoSglScrewPos
|
|
constitutive_nonlocal_stateInit( 3*ns+1: 4*ns) = rhoSglScrewNeg
|
|
constitutive_nonlocal_stateInit( 4*ns+1: 5*ns) = rhoSglEdgePosUsed
|
|
constitutive_nonlocal_stateInit( 5*ns+1: 6*ns) = rhoSglEdgeNegUsed
|
|
constitutive_nonlocal_stateInit( 6*ns+1: 7*ns) = rhoSglScrewPosUsed
|
|
constitutive_nonlocal_stateInit( 7*ns+1: 8*ns) = rhoSglScrewNegUsed
|
|
constitutive_nonlocal_stateInit( 8*ns+1: 9*ns) = rhoDipEdge
|
|
constitutive_nonlocal_stateInit( 9*ns+1:10*ns) = rhoDipScrew
|
|
|
|
endfunction
|
|
|
|
|
|
|
|
!*********************************************************************
|
|
!* absolute state tolerance *
|
|
!*********************************************************************
|
|
pure function constitutive_nonlocal_aTolState(myInstance)
|
|
|
|
use prec, only: pReal, &
|
|
pInt
|
|
implicit none
|
|
|
|
!*** input variables
|
|
integer(pInt), intent(in) :: myInstance ! number specifying the current instance of the constitution
|
|
|
|
!*** output variables
|
|
real(pReal), dimension(constitutive_nonlocal_sizeState(myInstance)) :: &
|
|
constitutive_nonlocal_aTolState ! absolute state tolerance for the current instance of this constitution
|
|
|
|
!*** local variables
|
|
|
|
constitutive_nonlocal_aTolState = constitutive_nonlocal_aTolRho(myInstance)
|
|
|
|
endfunction
|
|
|
|
|
|
|
|
!*********************************************************************
|
|
!* calculates homogenized elacticity matrix *
|
|
!*********************************************************************
|
|
pure function constitutive_nonlocal_homogenizedC(state,g,ip,el)
|
|
|
|
use prec, only: pReal, &
|
|
pInt, &
|
|
p_vec
|
|
use mesh, only: mesh_NcpElems, &
|
|
mesh_maxNips
|
|
use material, only: homogenization_maxNgrains, &
|
|
material_phase, &
|
|
phase_constitutionInstance
|
|
implicit none
|
|
|
|
!*** input variables
|
|
integer(pInt), intent(in) :: g, & ! current grain ID
|
|
ip, & ! current integration point
|
|
el ! current element
|
|
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state ! microstructural state
|
|
|
|
!*** output variables
|
|
real(pReal), dimension(6,6) :: constitutive_nonlocal_homogenizedC ! homogenized elasticity matrix
|
|
|
|
!*** local variables
|
|
integer(pInt) myInstance ! current instance of this constitution
|
|
|
|
myInstance = phase_constitutionInstance(material_phase(g,ip,el))
|
|
|
|
constitutive_nonlocal_homogenizedC = constitutive_nonlocal_Cslip_66(1:6,1:6,myInstance)
|
|
|
|
endfunction
|
|
|
|
|
|
|
|
!*********************************************************************
|
|
!* calculates quantities characterizing the microstructure *
|
|
!*********************************************************************
|
|
subroutine constitutive_nonlocal_microstructure(state, Temperature, Tstar_v, Fe, Fp, g, ip, el)
|
|
|
|
use prec, only: pReal, &
|
|
pInt, &
|
|
p_vec
|
|
use math, only: math_Plain3333to99, &
|
|
math_Mandel33to6, &
|
|
math_Mandel6to33, &
|
|
math_mul33x33, &
|
|
math_mul3x3, &
|
|
math_mul33x3, &
|
|
math_inv3x3, &
|
|
math_det3x3, &
|
|
math_transpose3x3, &
|
|
pi
|
|
use debug, only: debugger, &
|
|
verboseDebugger
|
|
use mesh, only: mesh_NcpElems, &
|
|
mesh_maxNips, &
|
|
mesh_maxNipNeighbors, &
|
|
mesh_element, &
|
|
FE_NipNeighbors, &
|
|
mesh_ipNeighborhood, &
|
|
mesh_ipVolume, &
|
|
mesh_ipCenterOfGravity
|
|
use material, only: homogenization_maxNgrains, &
|
|
material_phase, &
|
|
phase_localConstitution, &
|
|
phase_constitutionInstance
|
|
use lattice, only: lattice_Sslip, &
|
|
lattice_Sslip_v, &
|
|
lattice_maxNslipFamily, &
|
|
lattice_NslipSystem, &
|
|
lattice_maxNslip, &
|
|
lattice_sd, &
|
|
lattice_sn, &
|
|
lattice_st
|
|
|
|
implicit none
|
|
|
|
!*** input variables
|
|
integer(pInt), intent(in) :: g, & ! current grain ID
|
|
ip, & ! current integration point
|
|
el ! current element
|
|
real(pReal), intent(in) :: Temperature ! temperature
|
|
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
|
|
Fe, & ! elastic deformation gradient
|
|
Fp ! plastic deformation gradient
|
|
real(pReal), dimension(6), intent(in) :: &
|
|
Tstar_v ! 2nd Piola-Kirchhoff stress in Mandel notation
|
|
|
|
|
|
!*** input/output variables
|
|
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(inout) :: &
|
|
state ! microstructural state
|
|
|
|
!*** output variables
|
|
|
|
!*** local variables
|
|
integer(pInt) myInstance, & ! current instance of this constitution
|
|
myStructure, & ! current lattice structure
|
|
myPhase, &
|
|
ns, & ! short notation for the total number of active slip systems
|
|
neighboring_el, & ! element number of my neighbor
|
|
neighboring_ip, & ! integration point of my neighbor
|
|
c, & ! index of dilsocation character (edge, screw)
|
|
n, & ! index of my current neighbor
|
|
s, & ! index of my current slip system
|
|
t, & ! index of dilsocation type (e+, e-, s+, s-, used e+, used e-, used s+, used s-)
|
|
sLattice, & ! index of my current slip system according to lattice order
|
|
i, &
|
|
j
|
|
real(pReal) nu ! poisson's ratio
|
|
real(pReal), dimension(3,2) :: rhoExcessDifference, & ! finite differences of excess density (in 3 directions for edge and screw)
|
|
disloGradients ! spatial gradient in excess dislocation density (in 3 directions for edge and screw)
|
|
real(pReal), dimension(3,3) :: sigma, & ! dislocation stress for one slip system in its slip system frame
|
|
lattice2slip, & ! orthogonal transformation matrix from lattice coordinate system to slip coordinate system with e1=bxn, e2=b, e3=n (passive rotation!!!)
|
|
F, & ! total deformation gradient
|
|
neighboring_F, & ! total deformation gradient of neighbor
|
|
invFe, & ! inverse elastic deformation gradient
|
|
Q ! inverse transpose of 3x3 matrix with finite differences of opposing position vectors
|
|
real(pReal), dimension(6) :: Tdislocation_v ! dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress in Mandel notation
|
|
real(pReal), dimension(2,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: &
|
|
rhoExcess ! central excess density
|
|
real(pReal), dimension(6,2,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: &
|
|
neighboring_rhoExcess ! excess density for each neighbor, dislo character and slip system
|
|
real(pReal), dimension(3,6) :: neighboring_position ! position vector of each neighbor when seen from the centreal material point's lattice frame
|
|
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),8) :: &
|
|
rhoSgl ! single dislocation density (edge+, edge-, screw+, screw-, used edge+, used edge-, used screw+, used screw-)
|
|
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),2) :: &
|
|
rhoDip ! dipole dislocation density (edge, screw)
|
|
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: &
|
|
transmissivity, & ! transmissivity
|
|
rhoForest, & ! forest dislocation density
|
|
tauThreshold, & ! threshold shear stress
|
|
tau ! resolved shear stress
|
|
|
|
myPhase = material_phase(g,ip,el)
|
|
myInstance = phase_constitutionInstance(myPhase)
|
|
myStructure = constitutive_nonlocal_structure(myInstance)
|
|
ns = constitutive_nonlocal_totalNslip(myInstance)
|
|
|
|
|
|
!**********************************************************************
|
|
!*** get basic states
|
|
|
|
forall (t = 1:4) rhoSgl(1:ns,t) = max(state(g,ip,el)%p((t-1)*ns+1:t*ns), 0.0_pReal) ! ensure positive single mobile densities
|
|
forall (t = 5:8) rhoSgl(1:ns,t) = state(g,ip,el)%p((t-1)*ns+1:t*ns)
|
|
forall (c = 1:2) rhoDip(1:ns,c) = max(state(g,ip,el)%p((c+7)*ns+1:(c+8)*ns), 0.0_pReal) ! ensure positive dipole densities
|
|
where(rhoSgl(1:ns,1:4) < min(0.1, 0.01*constitutive_nonlocal_aTolRho(myInstance))) &
|
|
rhoSgl(1:ns,1:4) = 0.0_pReal ! delete non-significant single density
|
|
|
|
|
|
!**********************************************************************
|
|
!*** calculate dependent states
|
|
|
|
!*** calculate the forest dislocation density
|
|
|
|
forall (s = 1:ns) &
|
|
rhoForest(s) = dot_product( ( sum(abs(rhoSgl(1:ns,(/1,2,5,6/))),2) + rhoDip(1:ns,1) ), &
|
|
constitutive_nonlocal_forestProjectionEdge(s, 1:ns, myInstance) ) &
|
|
+ dot_product( ( sum(abs(rhoSgl(1:ns,(/3,4,7,8/))),2) + rhoDip(1:ns,2) ), &
|
|
constitutive_nonlocal_forestProjectionScrew(s, 1:ns, myInstance) ) ! calculation of forest dislocation density as projection of screw and edge dislocations
|
|
! if (debugger) write(6,'(a30,3(i3,x),/,12(e10.3,x),/)') 'forest dislocation density at ',g,ip,el, rhoForest
|
|
|
|
|
|
!*** calculate the threshold shear stress for dislocation slip
|
|
|
|
forall (s = 1:ns) &
|
|
tauThreshold(s) = constitutive_nonlocal_Gmod(myInstance) &
|
|
* constitutive_nonlocal_burgersPerSlipSystem(s, myInstance) &
|
|
* sqrt( dot_product( ( sum(abs(rhoSgl),2) + sum(abs(rhoDip),2) ), &
|
|
constitutive_nonlocal_interactionMatrixSlipSlip(s, 1:ns, myInstance) ) )
|
|
! if (debugger) write(6,'(a22,3(i3,x),/,12(f10.5,x),/)') 'tauThreshold / MPa at ',g,ip,el, tauThreshold/1e6
|
|
|
|
|
|
!*** calculate the dislocation stress of the neighboring excess dislocation densities
|
|
|
|
Tdislocation_v = 0.0_pReal
|
|
|
|
if (.not. phase_localConstitution(myPhase)) then ! only calculate dislocation stress for nonlocal material
|
|
|
|
F = math_mul33x33(Fe(1:3,1:3,g,ip,el), Fp(1:3,1:3,g,ip,el))
|
|
invFe = math_inv3x3(Fe(1:3,1:3,g,ip,el))
|
|
nu = constitutive_nonlocal_nu(myInstance)
|
|
|
|
forall (s = 1:ns, c = 1:2) &
|
|
rhoExcess(c,s) = state(g,ip,el)%p((2*c-2)*ns+s) + abs(state(g,ip,el)%p((2*c+2)*ns+s)) &
|
|
- state(g,ip,el)%p((2*c-1)*ns+s) - abs(state(g,ip,el)%p((2*c+3)*ns+s))
|
|
do n = 1,6
|
|
neighboring_el = mesh_ipNeighborhood(1,n,ip,el)
|
|
neighboring_ip = mesh_ipNeighborhood(2,n,ip,el)
|
|
if ( neighboring_ip == 0 .or. neighboring_el == 0 ) then ! at free surfaces ...
|
|
neighboring_el = el ! ... use central values instead of neighboring values
|
|
neighboring_ip = ip
|
|
neighboring_position(1:3,n) = 0.0_pReal
|
|
neighboring_rhoExcess(n,1:2,1:ns) = rhoExcess
|
|
elseif (phase_localConstitution(material_phase(1,neighboring_ip,neighboring_el))) then ! for neighbors with local constitution
|
|
neighboring_el = el ! ... use central values instead of neighboring values
|
|
neighboring_ip = ip
|
|
neighboring_position(1:3,n) = 0.0_pReal
|
|
neighboring_rhoExcess(n,1:2,1:ns) = rhoExcess
|
|
elseif (myPhase /= material_phase(1,neighboring_ip,neighboring_el)) then ! for neighbors with different phase
|
|
neighboring_el = el ! ... use central values instead of neighboring values
|
|
neighboring_ip = ip
|
|
neighboring_position(1:3,n) = 0.0_pReal
|
|
neighboring_rhoExcess(n,1:2,1:ns) = rhoExcess
|
|
else
|
|
transmissivity = sum(constitutive_nonlocal_compatibility(2,1:ns,1:ns,n,ip,el)**2.0_pReal, 1)
|
|
if ( any(transmissivity < 0.9_pReal) ) then ! at grain boundary (=significantly decreased transmissivity) ...
|
|
neighboring_el = el ! ... use central values instead of neighboring values
|
|
neighboring_ip = ip
|
|
neighboring_position(1:3,n) = 0.0_pReal
|
|
neighboring_rhoExcess(n,1:2,1:ns) = rhoExcess
|
|
else
|
|
neighboring_F = math_mul33x33(Fe(1:3,1:3,g,neighboring_ip,neighboring_el), Fp(1:3,1:3,g,neighboring_ip,neighboring_el))
|
|
neighboring_position(1:3,n) = &
|
|
0.5_pReal * math_mul33x3(math_mul33x33(invFe,neighboring_F) + Fp(1:3,1:3,g,ip,el), &
|
|
mesh_ipCenterOfGravity(1:3,neighboring_ip,neighboring_el) - mesh_ipCenterOfGravity(1:3,ip,el))
|
|
forall (s = 1:ns, c = 1:2) &
|
|
neighboring_rhoExcess(n,c,s) = state(g,neighboring_ip,neighboring_el)%p((2*c-2)*ns+s) &
|
|
+ abs(state(g,neighboring_ip,neighboring_el)%p((2*c+2)*ns+s)) &
|
|
- state(g,neighboring_ip,neighboring_el)%p((2*c-1)*ns+s) &
|
|
- abs(state(g,neighboring_ip,neighboring_el)%p((2*c+3)*ns+s))
|
|
endif
|
|
endif
|
|
enddo
|
|
|
|
Q = math_inv3x3(math_transpose3x3(neighboring_position(1:3,(/1,3,5/)) - neighboring_position(1:3,(/2,4,6/))))
|
|
|
|
do s = 1,ns
|
|
|
|
lattice2slip = math_transpose3x3(reshape((/ &
|
|
lattice_sd(1:3, constitutive_nonlocal_slipSystemLattice(s,myInstance), myStructure), &
|
|
-lattice_st(1:3, constitutive_nonlocal_slipSystemLattice(s,myInstance), myStructure), &
|
|
lattice_sn(1:3, constitutive_nonlocal_slipSystemLattice(s,myInstance), myStructure)/), (/3,3/)))
|
|
|
|
rhoExcessDifference = neighboring_rhoExcess((/1,3,5/),1:2,s) - neighboring_rhoExcess((/2,4,6/),1:2,s)
|
|
forall (c = 1:2) &
|
|
disloGradients(1:3,c) = math_mul33x3(lattice2slip, math_mul33x3(Q, rhoExcessDifference(1:3,c)))
|
|
|
|
sigma = 0.0_pReal
|
|
sigma(1,1) = + 0.375_pReal / (1.0_pReal-nu) * disloGradients(3,1)
|
|
sigma(2,2) = + 0.5_pReal * nu / (1.0_pReal-nu) * disloGradients(3,1)
|
|
sigma(3,3) = + 0.125_pReal / (1.0_pReal-nu) * disloGradients(3,1)
|
|
sigma(1,2) = + 0.25_pReal * disloGradients(3,2)
|
|
sigma(1,3) = - 0.125_pReal / (1.0_pReal-nu) * disloGradients(1,1) - 0.25_pReal * disloGradients(2,2)
|
|
sigma(2,1) = sigma(1,2)
|
|
sigma(3,1) = sigma(1,3)
|
|
|
|
forall (i=1:3, j=1:3) &
|
|
sigma(i,j) = sigma(i,j) * constitutive_nonlocal_Gmod(myInstance) * constitutive_nonlocal_burgersPerSlipSystem(s,myInstance) &
|
|
* constitutive_nonlocal_R(myInstance)**2.0_pReal
|
|
|
|
Tdislocation_v = Tdislocation_v + math_Mandel33to6(math_mul33x33(math_transpose3x3(lattice2slip), &
|
|
math_mul33x33(sigma, lattice2slip)))
|
|
|
|
enddo
|
|
|
|
endif
|
|
|
|
|
|
!**********************************************************************
|
|
!*** set states
|
|
|
|
state(g,ip,el)%p(1:8*ns) = reshape(rhoSgl,(/8*ns/)) ! ensure positive single mobile densities
|
|
state(g,ip,el)%p(8*ns+1:10*ns) = reshape(rhoDip,(/2*ns/)) ! ensure positive dipole densities
|
|
state(g,ip,el)%p(10*ns+1:11*ns) = rhoForest
|
|
state(g,ip,el)%p(11*ns+1:12*ns) = tauThreshold
|
|
state(g,ip,el)%p(12*ns+1:12*ns+6) = Tdislocation_v
|
|
|
|
endsubroutine
|
|
|
|
|
|
|
|
!*********************************************************************
|
|
!* calculates kinetics *
|
|
!*********************************************************************
|
|
subroutine constitutive_nonlocal_kinetics(Tstar_v, Temperature, state, g, ip, el, dv_dtau)
|
|
|
|
use prec, only: pReal, &
|
|
pInt, &
|
|
p_vec
|
|
use math, only: math_mul6x6, &
|
|
math_Mandel6to33
|
|
use debug, only: debugger, &
|
|
verboseDebugger, &
|
|
debug_g, &
|
|
debug_i, &
|
|
debug_e
|
|
use mesh, only: mesh_NcpElems, &
|
|
mesh_maxNips
|
|
use material, only: homogenization_maxNgrains, &
|
|
material_phase, &
|
|
phase_constitutionInstance
|
|
use lattice, only: lattice_Sslip, &
|
|
lattice_Sslip_v
|
|
|
|
implicit none
|
|
|
|
!*** input variables
|
|
integer(pInt), intent(in) :: g, & ! current grain number
|
|
ip, & ! current integration point
|
|
el ! current element number
|
|
real(pReal), intent(in) :: Temperature ! temperature
|
|
type(p_vec), intent(in) :: state ! microstructural state
|
|
real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola-Kirchhoff stress in Mandel notation
|
|
|
|
!*** output variables
|
|
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))), &
|
|
intent(out), optional :: dv_dtau ! velocity derivative with respect to resolved shear stress
|
|
|
|
!*** local variables
|
|
integer(pInt) myInstance, & ! current instance of this constitution
|
|
myStructure, & ! current lattice structure
|
|
ns, & ! short notation for the total number of active slip systems
|
|
t, & ! dislocation type
|
|
s ! index of my current slip system
|
|
real(pReal), dimension(6) :: Tdislocation_v ! dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress
|
|
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: &
|
|
tauThreshold, & ! threshold shear stress
|
|
tau, & ! resolved shear stress
|
|
rhoForest ! forest dislocation density
|
|
real(pReal) boltzmannProbability, &
|
|
tauRel, & ! relative thermally active resolved shear stress
|
|
wallFunc, & ! functions reflecting the shape of the obstacle wall (see PhD thesis Mohles p.53)
|
|
timeRatio ! ratio of travel to dwell time
|
|
|
|
|
|
myInstance = phase_constitutionInstance(material_phase(g,ip,el))
|
|
myStructure = constitutive_nonlocal_structure(myInstance)
|
|
ns = constitutive_nonlocal_totalNslip(myInstance)
|
|
|
|
rhoForest = state%p(10*ns+1:11*ns)
|
|
tauThreshold = state%p(11*ns+1:12*ns)
|
|
Tdislocation_v = state%p(12*ns+1:12*ns+6)
|
|
|
|
tau = 0.0_pReal
|
|
constitutive_nonlocal_v(1:ns,1:4,g,ip,el) = 0.0_pReal
|
|
if (present(dv_dtau)) dv_dtau = 0.0_pReal
|
|
|
|
|
|
if (Temperature > 0.0_pReal) then
|
|
do s = 1,ns
|
|
|
|
tau(s) = math_mul6x6(Tstar_v + Tdislocation_v, &
|
|
lattice_Sslip_v(:,constitutive_nonlocal_slipSystemLattice(s,myInstance),myStructure))
|
|
|
|
!*** Only if the resolved shear stress exceeds the threshold stress, dislocations are able to cut the dislocation forest.
|
|
!*** In contrast to small atomic obstacles the forest can't be overcome by thermal activation.
|
|
!***
|
|
!*** mean travel distance
|
|
!*** The mean dislocation velocity is calculated as: --------------------------
|
|
!*** dwell time + travel time
|
|
!***
|
|
!*** with : mean travel distance = inverse of the root of forest density
|
|
!*** dwell time = inverse of attack frequency times probability of success
|
|
!*** travel time = mean travel distance over velocity of sound
|
|
|
|
tauRel = (abs(tau(s)) - tauThreshold(s)) / constitutive_nonlocal_tauObs(myInstance)
|
|
if (tauRel > 0.0_pReal .and. tauRel < 1.0_pReal) then
|
|
wallFunc = 4.0_pReal * dsqrt(2.0_pReal) / 3.0_pReal * dsqrt(1.0_pReal - tauRel) / tauRel
|
|
boltzmannProbability = dexp(- constitutive_nonlocal_Qeff0(s,myInstance) * wallFunc / (kB * Temperature))
|
|
timeRatio = boltzmannProbability * constitutive_nonlocal_fattack(myInstance) &
|
|
/ (constitutive_nonlocal_vs(myInstance) * dsqrt(rhoForest(s)))
|
|
|
|
constitutive_nonlocal_v(s,:,g,ip,el) = sign(constitutive_nonlocal_vs(myInstance),tau(s)) * timeRatio / (1.0_pReal + timeRatio)
|
|
|
|
if (present(dv_dtau)) then
|
|
dv_dtau(s) = abs(constitutive_nonlocal_v(s,1,g,ip,el)) * constitutive_nonlocal_Qeff0(s,myInstance) &
|
|
/ (kB * Temperature * (1.0_pReal + timeRatio)) &
|
|
* 0.5_pReal * wallFunc * (2.0_pReal - tauRel) &
|
|
/ ((1.0_pReal - tauRel) * (abs(tau(s)) - tauThreshold(s)))
|
|
endif
|
|
|
|
!*** If resolved stress exceeds threshold plus obstacle stress, the probability for thermal activation is 1.
|
|
!*** The tangent is zero, since no dependency of tau.
|
|
|
|
elseif (tauRel >= 1.0_pReal) then
|
|
constitutive_nonlocal_v(s,1:4,g,ip,el) = sign(constitutive_nonlocal_vs(myInstance), tau(s)) &
|
|
* constitutive_nonlocal_fattack(myInstance) &
|
|
/ (constitutive_nonlocal_vs(myInstance) * dsqrt(rhoForest(s)) &
|
|
+ constitutive_nonlocal_fattack(myInstance))
|
|
endif
|
|
enddo
|
|
endif
|
|
|
|
!if (verboseDebugger .and. s) then
|
|
! !$OMP CRITICAL (write2out)
|
|
! write(6,*) '::: kinetics',g,ip,el
|
|
! write(6,*)
|
|
! write(6,'(a,/,3(3(f12.3,x)/))') 'Tdislocation / MPa', math_Mandel6to33(Tdislocation_v/1e6)
|
|
! write(6,'(a,/,3(3(f12.3,x)/))') 'Tstar / MPa', math_Mandel6to33(Tstar_v/1e6)
|
|
! write(6,'(a,/,12(f12.5,x),/)') 'tau / MPa', tau/1e6_pReal
|
|
! write(6,'(a,/,12(e12.5,x),/)') 'rhoForest / 1/m**2', rhoForest
|
|
! write(6,'(a,/,4(12(f12.5,x),/))') 'v / 1e-3m/s', constitutive_nonlocal_v(:,:,g,ip,el)*1e3
|
|
! !$OMP END CRITICAL (write2out)
|
|
!endif
|
|
|
|
endsubroutine
|
|
|
|
|
|
|
|
!*********************************************************************
|
|
!* calculates plastic velocity gradient and its tangent *
|
|
!*********************************************************************
|
|
subroutine constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar99, Tstar_v, Temperature, state, g, ip, el)
|
|
|
|
use prec, only: pReal, &
|
|
pInt, &
|
|
p_vec
|
|
use math, only: math_Plain3333to99, &
|
|
math_mul6x6, &
|
|
math_Mandel6to33
|
|
use debug, only: debugger, &
|
|
verboseDebugger, &
|
|
debug_g, &
|
|
debug_i, &
|
|
debug_e
|
|
use mesh, only: mesh_NcpElems, &
|
|
mesh_maxNips
|
|
use material, only: homogenization_maxNgrains, &
|
|
material_phase, &
|
|
phase_constitutionInstance
|
|
use lattice, only: lattice_Sslip, &
|
|
lattice_Sslip_v
|
|
|
|
implicit none
|
|
|
|
!*** input variables
|
|
integer(pInt), intent(in) :: g, & ! current grain number
|
|
ip, & ! current integration point
|
|
el ! current element number
|
|
real(pReal), intent(in) :: Temperature ! temperature
|
|
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
|
|
state ! microstructural state
|
|
real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola-Kirchhoff stress in Mandel notation
|
|
|
|
!*** output variables
|
|
real(pReal), dimension(3,3), intent(out) :: Lp ! plastic velocity gradient
|
|
real(pReal), dimension(9,9), intent(out) :: dLp_dTstar99 ! derivative of Lp with respect to Tstar (9x9 matrix)
|
|
|
|
!*** local variables
|
|
integer(pInt) myInstance, & ! current instance of this constitution
|
|
myStructure, & ! current lattice structure
|
|
ns, & ! short notation for the total number of active slip systems
|
|
i, &
|
|
j, &
|
|
k, &
|
|
l, &
|
|
t, & ! dislocation type
|
|
s, & ! index of my current slip system
|
|
sLattice ! index of my current slip system according to lattice order
|
|
real(pReal), dimension(3,3,3,3) :: dLp_dTstar3333 ! derivative of Lp with respect to Tstar (3x3x3x3 matrix)
|
|
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),8) :: &
|
|
rhoSgl ! single dislocation densities (including used)
|
|
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: &
|
|
gdot ! shear rate per dislocation type
|
|
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: &
|
|
tauThreshold, & ! threshold shear stress
|
|
gdotTotal, & ! shear rate
|
|
dv_dtau, & ! velocity derivative with respect to the shear stress
|
|
dgdotTotal_dtau, & ! derivative of the shear rate with respect to the shear stress
|
|
rhoForest ! forest dislocation density
|
|
|
|
|
|
!*** initialize local variables
|
|
|
|
gdot = 0.0_pReal
|
|
Lp = 0.0_pReal
|
|
dLp_dTstar3333 = 0.0_pReal
|
|
|
|
myInstance = phase_constitutionInstance(material_phase(g,ip,el))
|
|
myStructure = constitutive_nonlocal_structure(myInstance)
|
|
ns = constitutive_nonlocal_totalNslip(myInstance)
|
|
|
|
!*** shortcut to state variables
|
|
|
|
forall (t = 1:8) &
|
|
rhoSgl(1:ns,t) = state(g,ip,el)%p((t-1)*ns+1:t*ns)
|
|
forall (s = 1:ns, t = 5:8, rhoSgl(s,t) * constitutive_nonlocal_v(s,t-4,g,ip,el) < 0.0_pReal) & ! contribution of used rho for changing sign of v
|
|
rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t))
|
|
|
|
rhoForest = state(g,ip,el)%p(10*ns+1:11*ns)
|
|
tauThreshold = state(g,ip,el)%p(11*ns+1:12*ns)
|
|
|
|
call constitutive_nonlocal_kinetics(Tstar_v, Temperature, state(g,ip,el), g, ip, el, dv_dtau) ! update dislocation velocity
|
|
|
|
!*** Calculation of gdot and its tangent
|
|
|
|
forall (t = 1:4) &
|
|
gdot(1:ns,t) = rhoSgl(1:ns,t) * constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) &
|
|
* constitutive_nonlocal_v(1:ns,t,g,ip,el)
|
|
gdotTotal = sum(gdot,2)
|
|
|
|
dgdotTotal_dtau = sum(rhoSgl,2) * constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) * dv_dtau
|
|
|
|
!*** Calculation of Lp and its tangent
|
|
|
|
do s = 1,ns
|
|
sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance)
|
|
|
|
Lp = Lp + gdotTotal(s) * lattice_Sslip(1:3,1:3,sLattice,myStructure)
|
|
|
|
forall (i=1:3,j=1:3,k=1:3,l=1:3) &
|
|
dLp_dTstar3333(i,j,k,l) = dLp_dTstar3333(i,j,k,l) + dgdotTotal_dtau(s) * lattice_Sslip(i,j, sLattice,myStructure) &
|
|
* lattice_Sslip(k,l, sLattice,myStructure)
|
|
enddo
|
|
|
|
dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333)
|
|
|
|
!if (verboseDebugger .and. (debug_g==g .and. debug_i==i .and. debug_e==e)) then
|
|
! !$OMP CRITICAL (write2out)
|
|
! write(6,*) '::: LpandItsTangent',g,ip,el
|
|
! write(6,*)
|
|
! write(6,'(a,/,12(f12.5,x),/)') 'v / 1e-3m/s', constitutive_nonlocal_v(:,:,g,ip,el)*1e3
|
|
! write(6,'(a,/,12(f12.5,x),/)') 'gdot / 1e-3',gdot*1e3_pReal
|
|
! write(6,'(a,/,12(f12.5,x),/)') 'gdot total / 1e-3',gdotTotal*1e3_pReal
|
|
! write(6,'(a,/,3(3(f12.7,x)/))') 'Lp',Lp
|
|
! ! call flush(6)
|
|
! !$OMP END CRITICAL (write2out)
|
|
!endif
|
|
|
|
endsubroutine
|
|
|
|
|
|
|
|
!*********************************************************************
|
|
!* rate of change of microstructure *
|
|
!*********************************************************************
|
|
subroutine constitutive_nonlocal_dotState(dotState, Tstar_v, previousTstar_v, Fe, Fp, Temperature, dt_previous, &
|
|
state, previousState, aTolState, timestep, orientation, g,ip,el)
|
|
|
|
use prec, only: pReal, &
|
|
pInt, &
|
|
p_vec
|
|
use IO, only: IO_error
|
|
use debug, only: debugger, &
|
|
debug_g, &
|
|
debug_i, &
|
|
debug_e, &
|
|
verboseDebugger
|
|
use math, only: math_norm3, &
|
|
math_mul6x6, &
|
|
math_mul3x3, &
|
|
math_mul33x3, &
|
|
math_mul33x33, &
|
|
math_inv3x3, &
|
|
math_det3x3, &
|
|
math_Mandel6to33, &
|
|
math_QuaternionDisorientation, &
|
|
math_qRot, &
|
|
pi, &
|
|
NaN
|
|
use mesh, only: mesh_NcpElems, &
|
|
mesh_maxNips, &
|
|
mesh_maxNipNeighbors, &
|
|
mesh_element, &
|
|
FE_NipNeighbors, &
|
|
mesh_ipNeighborhood, &
|
|
mesh_ipVolume, &
|
|
mesh_ipArea, &
|
|
mesh_ipAreaNormal, &
|
|
mesh_ipCenterOfGravity
|
|
use material, only: homogenization_maxNgrains, &
|
|
material_phase, &
|
|
phase_constitutionInstance, &
|
|
phase_localConstitution
|
|
use lattice, only: lattice_Sslip, &
|
|
lattice_Sslip_v, &
|
|
lattice_sd, &
|
|
lattice_sn, &
|
|
lattice_st, &
|
|
lattice_maxNslipFamily, &
|
|
lattice_NslipSystem
|
|
use FEsolving, only:theInc, &
|
|
FEsolving_execElem, &
|
|
FEsolving_execIP
|
|
|
|
implicit none
|
|
|
|
!*** input variables
|
|
integer(pInt), intent(in) :: g, & ! current grain number
|
|
ip, & ! current integration point
|
|
el ! current element number
|
|
real(pReal), intent(in) :: Temperature, & ! temperature
|
|
timestep, & ! substepped crystallite time increment
|
|
dt_previous ! time increment between previous and current state
|
|
real(pReal), dimension(6), intent(in) :: Tstar_v, & ! current 2nd Piola-Kirchhoff stress in Mandel notation
|
|
previousTstar_v ! previous 2nd Piola-Kirchhoff stress in Mandel notation
|
|
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
|
|
Fe, & ! elastic deformation gradient
|
|
Fp ! plastic deformation gradient
|
|
real(pReal), dimension(4,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
|
|
orientation ! crystal lattice orientation
|
|
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
|
|
state, & ! current microstructural state
|
|
previousState, & ! previous microstructural state
|
|
aTolState ! absolute state tolerance
|
|
|
|
!*** input/output variables
|
|
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(inout) :: &
|
|
dotState ! evolution of state variables / microstructure
|
|
|
|
!*** output variables
|
|
|
|
!*** local variables
|
|
integer(pInt) myInstance, & ! current instance of this constitution
|
|
myStructure, & ! current lattice structure
|
|
ns, & ! short notation for the total number of active slip systems
|
|
c, & ! character of dislocation
|
|
n, & ! index of my current neighbor
|
|
neighboring_el, & ! element number of my neighbor
|
|
neighboring_ip, & ! integration point of my neighbor
|
|
neighboring_n, & ! neighbor index pointing to me when looking from my neighbor
|
|
opposite_n, & ! index of my opposite neighbor
|
|
opposite_ip, & ! ip of my opposite neighbor
|
|
opposite_el, & ! element index of my opposite neighbor
|
|
t, & ! type of dislocation
|
|
topp, & ! type of dislocation with opposite sign to t
|
|
s, & ! index of my current slip system
|
|
sLattice, & ! index of my current slip system according to lattice order
|
|
i
|
|
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),10) :: &
|
|
rhoDot, & ! density evolution
|
|
rhoDotRemobilization, & ! density evolution by remobilization
|
|
rhoDotMultiplication, & ! density evolution by multiplication
|
|
rhoDotFlux, & ! density evolution by flux
|
|
neighboring_rhoDotFlux, & ! density evolution by flux at neighbor
|
|
rhoDotSingle2DipoleGlide, & ! density evolution by dipole formation (by glide)
|
|
rhoDotAthermalAnnihilation, & ! density evolution by athermal annihilation
|
|
rhoDotThermalAnnihilation, & ! density evolution by thermal annihilation
|
|
rhoDotDipole2SingleStressChange, & ! density evolution by dipole dissociation (by stress increase)
|
|
rhoDotSingle2DipoleStressChange ! density evolution by dipole formation (by stress decrease)
|
|
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),8) :: &
|
|
rhoSgl, & ! current single dislocation densities (positive/negative screw and edge without dipoles)
|
|
previousRhoSgl ! previous single dislocation densities (positive/negative screw and edge without dipoles)
|
|
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: &
|
|
fluxdensity, & ! flux density at central material point
|
|
neighboring_fluxdensity, & ! flux density at neighboring material point
|
|
gdot ! shear rates
|
|
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: &
|
|
rhoForest, & ! forest dislocation density
|
|
tauThreshold, & ! threshold shear stress
|
|
tau, & ! current resolved shear stress
|
|
previousTau, & ! previous resolved shear stress
|
|
invLambda, & ! inverse of mean free path for dislocations
|
|
vClimb ! climb velocity of edge dipoles
|
|
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),2) :: &
|
|
rhoDip, & ! current dipole dislocation densities (screw and edge dipoles)
|
|
previousRhoDip, & ! previous 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
|
|
previousDUpper, & ! previous maximum stable dipole distance for edges and screws
|
|
dUpperDot ! rate of change of the maximum stable dipole distance for edges and screws
|
|
real(pReal), dimension(3,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: &
|
|
m ! direction of dislocation motion
|
|
real(pReal), dimension(3,3) :: my_F, & ! my total deformation gradient
|
|
neighboring_F, & ! total deformation gradient of my neighbor
|
|
my_Fe, & ! my elastic deformation gradient
|
|
neighboring_Fe, & ! elastic deformation gradient of my neighbor
|
|
Favg ! average total deformation gradient of me and my neighbor
|
|
real(pReal), dimension(6) :: Tdislocation_v, & ! current dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress
|
|
previousTdislocation_v ! previous dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress
|
|
real(pReal), dimension(3) :: normal_neighbor2me, & ! interface normal pointing from my neighbor to me in neighbor's lattice configuration
|
|
normal_neighbor2me_defConf, & ! interface normal pointing from my neighbor to me in shared deformed configuration
|
|
normal_me2neighbor, & ! interface normal pointing from me to my neighbor in my lattice configuration
|
|
normal_me2neighbor_defConf ! interface normal pointing from me to my neighbor in shared deformed configuration
|
|
real(pReal) area, & ! area of the current interface
|
|
transmissivity, & ! overall transmissivity of dislocation flux to neighboring material point
|
|
lineLength, & ! dislocation line length leaving the current interface
|
|
D, & ! self diffusion
|
|
correction
|
|
logical considerEnteringFlux, &
|
|
considerLeavingFlux
|
|
|
|
if (verboseDebugger .and. (debug_g==g .and. debug_i==ip .and. debug_e==el)) then
|
|
!$OMP CRITICAL (write2out)
|
|
write(6,*) '::: constitutive_nonlocal_dotState at ',g,ip,el
|
|
write(6,*)
|
|
!$OMP END CRITICAL (write2out)
|
|
endif
|
|
|
|
select case(mesh_element(2,el))
|
|
case (1,6,7,8,9)
|
|
! all fine
|
|
case default
|
|
call IO_error(-1,el,ip,g,'element type not supported for nonlocal constitution')
|
|
|
|
|
|
end select
|
|
|
|
myInstance = phase_constitutionInstance(material_phase(g,ip,el))
|
|
myStructure = constitutive_nonlocal_structure(myInstance)
|
|
ns = constitutive_nonlocal_totalNslip(myInstance)
|
|
|
|
tau = 0.0_pReal
|
|
previousTau = 0.0_pReal
|
|
gdot = 0.0_pReal
|
|
dLower = 0.0_pReal
|
|
dUpper = 0.0_pReal
|
|
previousDUpper = 0.0_pReal
|
|
dUpperDot = 0.0_pReal
|
|
|
|
!*** shortcut to state variables
|
|
|
|
forall (t = 1:8) rhoSgl(1:ns,t) = state(g,ip,el)%p((t-1)*ns+1:t*ns)
|
|
forall (t = 1:8) previousRhoSgl(1:ns,t) = previousState(g,ip,el)%p((t-1)*ns+1:t*ns)
|
|
forall (c = 1:2) rhoDip(1:ns,c) = state(g,ip,el)%p((7+c)*ns+1:(8+c)*ns)
|
|
forall (c = 1:2) previousRhoDip(1:ns,c) = previousState(g,ip,el)%p((7+c)*ns+1:(8+c)*ns)
|
|
rhoForest = state(g,ip,el)%p(10*ns+1:11*ns)
|
|
tauThreshold = state(g,ip,el)%p(11*ns+1:12*ns)
|
|
Tdislocation_v = state(g,ip,el)%p(12*ns+1:12*ns+6)
|
|
previousTdislocation_v = previousState(g,ip,el)%p(12*ns+1:12*ns+6)
|
|
|
|
!*** sanity check for timestep
|
|
|
|
if (timestep <= 0.0_pReal) then ! if illegal timestep...
|
|
dotState(g,ip,el)%p = 0.0_pReal ! ...return without doing anything (-> zero dotState)
|
|
return
|
|
endif
|
|
|
|
|
|
!****************************************************************************
|
|
!*** Calculate shear rate
|
|
|
|
call constitutive_nonlocal_kinetics(Tstar_v, Temperature, state(g,ip,el), g, ip, el) ! get velocities
|
|
|
|
forall (t = 1:4) &
|
|
gdot(1:ns,t) = rhoSgl(1:ns,t) * constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) &
|
|
* constitutive_nonlocal_v(1:ns,t,g,ip,el)
|
|
forall (s = 1:ns, t = 1:4, rhoSgl(s,t+4) * constitutive_nonlocal_v(s,t,g,ip,el) < 0.0_pReal) & ! contribution of used rho for changing sign of v
|
|
gdot(s,t) = gdot(s,t) + abs(rhoSgl(s,t+4)) * constitutive_nonlocal_burgersPerSlipSystem(s,myInstance) &
|
|
* constitutive_nonlocal_v(s,t,g,ip,el)
|
|
|
|
if (verboseDebugger .and. (debug_g==g .and. debug_i==ip .and. debug_e==el)) then
|
|
!$OMP CRITICAL (write2out)
|
|
write(6,'(a,/,10(12(e12.5,x),/))') 'rho / 1/m^2', rhoSgl, rhoDip
|
|
write(6,'(a,/,4(12(e12.5,x),/))') 'v / m/s', constitutive_nonlocal_v(:,:,g,ip,el)
|
|
write(6,'(a,/,4(12(e12.5,x),/))') 'gdot / 1/s',gdot
|
|
!$OMP END CRITICAL (write2out)
|
|
endif
|
|
|
|
|
|
!****************************************************************************
|
|
!*** calculate limits for stable dipole height and its rate of change
|
|
|
|
do s = 1,ns ! loop over slip systems
|
|
sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance)
|
|
|
|
tau(s) = math_mul6x6( Tstar_v + Tdislocation_v, lattice_Sslip_v(1:6,sLattice,myStructure) )
|
|
previousTau(s) = math_mul6x6( previousTstar_v + previousTdislocation_v, lattice_Sslip_v(1:6,sLattice,myStructure) )
|
|
|
|
enddo
|
|
|
|
dLower(1:ns,1) = constitutive_nonlocal_dLowerEdgePerSlipSystem(1:ns,myInstance)
|
|
dLower(1:ns,2) = constitutive_nonlocal_dLowerScrewPerSlipSystem(1:ns,myInstance)
|
|
dUpper(1:ns,2) = min( 1.0_pReal / sqrt( sum(abs(rhoSgl),2)+sum(rhoDip,2) ), &
|
|
constitutive_nonlocal_Gmod(myInstance) * constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) &
|
|
/ ( 8.0_pReal * pi * abs(tau) ) )
|
|
dUpper(1:ns,1) = dUpper(1:ns,2) / ( 1.0_pReal - constitutive_nonlocal_nu(myInstance) )
|
|
previousDUpper(1:ns,2) = min( 1.0_pReal / sqrt( sum(abs(previousRhoSgl),2) + sum(previousRhoDip,2) ), &
|
|
constitutive_nonlocal_Gmod(myInstance) * constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) &
|
|
/ ( 8.0_pReal * pi * abs(previousTau) ) )
|
|
previousDUpper(1:ns,1) = previousDUpper(1:ns,2) / ( 1.0_pReal - constitutive_nonlocal_nu(myInstance) )
|
|
|
|
if (dt_previous > 0.0_pReal) dUpperDot = (dUpper - previousDUpper) / dt_previous
|
|
|
|
|
|
!****************************************************************************
|
|
!*** dislocation remobilization (bauschinger effect)
|
|
|
|
rhoDotRemobilization = 0.0_pReal
|
|
if (timestep > 0.0_pReal) then
|
|
do t = 1,4
|
|
do s = 1,ns
|
|
if (rhoSgl(s,t+4) * constitutive_nonlocal_v(s,t,g,ip,el) < 0.0_pReal) then
|
|
rhoDotRemobilization(s,t) = abs(rhoSgl(s,t+4)) / timestep
|
|
rhoSgl(s,t) = rhoSgl(s,t) + abs(rhoSgl(s,t+4))
|
|
rhoDotRemobilization(s,t+4) = - rhoSgl(s,t+4) / timestep
|
|
rhoSgl(s,t+4) = 0.0_pReal
|
|
endif
|
|
enddo
|
|
enddo
|
|
endif
|
|
|
|
|
|
!****************************************************************************
|
|
!*** calculate dislocation multiplication
|
|
|
|
rhoDotMultiplication = 0.0_pReal
|
|
where (rhoSgl(1:ns,3:4) > 0.0_pReal) &
|
|
rhoDotMultiplication(1:ns,1:2) = spread(0.5_pReal * sum(abs(gdot(1:ns,3:4)),2) * sqrt(rhoForest) &
|
|
/ constitutive_nonlocal_lambda0PerSlipSystem(1:ns,myInstance) &
|
|
/ constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance), 2, 2)
|
|
where (rhoSgl(1:ns,1:2) > 0.0_pReal) &
|
|
rhoDotMultiplication(1:ns,3:4) = spread(0.5_pReal * sum(abs(gdot(1:ns,1:2)),2) * sqrt(rhoForest) &
|
|
/ constitutive_nonlocal_lambda0PerSlipSystem(1:ns,myInstance) &
|
|
/ constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance), 2, 2)
|
|
|
|
|
|
!****************************************************************************
|
|
!*** calculate dislocation fluxes (only for nonlocal constitution)
|
|
|
|
constitutive_nonlocal_rhoDotFlux(:,:,g,ip,el) = 0.0_pReal
|
|
rhoDotFlux = 0.0_pReal
|
|
|
|
if (.not. phase_localConstitution(material_phase(g,ip,el))) then ! only for nonlocal constitution
|
|
|
|
!*** take care of the definition of lattice_st = lattice_sd x lattice_sn !!!
|
|
!*** opposite sign to our p vector in the (s,p,n) triplet !!!
|
|
|
|
m(1:3,1:ns,1) = lattice_sd(1:3, constitutive_nonlocal_slipSystemLattice(1:ns,myInstance), myStructure)
|
|
m(1:3,1:ns,2) = -lattice_sd(1:3, constitutive_nonlocal_slipSystemLattice(1:ns,myInstance), myStructure)
|
|
m(1:3,1:ns,3) = -lattice_st(1:3, constitutive_nonlocal_slipSystemLattice(1:ns,myInstance), myStructure)
|
|
m(1:3,1:ns,4) = lattice_st(1:3, constitutive_nonlocal_slipSystemLattice(1:ns,myInstance), myStructure)
|
|
|
|
my_Fe = Fe(1:3,1:3,g,ip,el)
|
|
my_F = math_mul33x33(my_Fe, Fp(1:3,1:3,g,ip,el))
|
|
|
|
fluxdensity = rhoSgl(1:ns,1:4) * constitutive_nonlocal_v(1:ns,1:4,g,ip,el)
|
|
|
|
do n = 1,FE_NipNeighbors(mesh_element(2,el)) ! loop through my neighbors
|
|
|
|
neighboring_el = mesh_ipNeighborhood(1,n,ip,el)
|
|
neighboring_ip = mesh_ipNeighborhood(2,n,ip,el)
|
|
if (neighboring_el > 0_pInt .and. neighboring_ip > 0_pInt) then ! if neighbor exists ...
|
|
do neighboring_n = 1,FE_NipNeighbors(mesh_element(2,neighboring_el)) ! find neighboring index that points from my neighbor to myself
|
|
if ( el == mesh_ipNeighborhood(1,neighboring_n,neighboring_ip,neighboring_el) &
|
|
.and. ip == mesh_ipNeighborhood(2,neighboring_n,neighboring_ip,neighboring_el) ) &
|
|
exit
|
|
enddo
|
|
endif
|
|
|
|
opposite_n = n + mod(n,2) - mod(n+1,2)
|
|
opposite_el = mesh_ipNeighborhood(1,opposite_n,ip,el)
|
|
opposite_ip = mesh_ipNeighborhood(2,opposite_n,ip,el)
|
|
|
|
if (neighboring_el > 0_pInt .and. neighboring_ip > 0_pInt) then ! if neighbor exists, average deformation gradient
|
|
neighboring_Fe = Fe(1:3,1:3,g,neighboring_ip,neighboring_el)
|
|
neighboring_F = math_mul33x33(neighboring_Fe, Fp(1:3,1:3,g,neighboring_ip,neighboring_el))
|
|
Favg = 0.5_pReal * (my_F + neighboring_F)
|
|
else ! if no neighbor, take my value as average
|
|
Favg = my_F
|
|
endif
|
|
|
|
|
|
!* FLUX FROM ME TO MY NEIGHBOR
|
|
!* This is not considered, if my opposite neighbor has a local constitution.
|
|
!* Then, we assume, that the opposite(!) neighbor sends an equal amount of dislocations to me.
|
|
!* So the net flux in the direction of my neighbor is equal to zero:
|
|
!* leaving flux to neighbor == entering flux from opposite neighbor
|
|
!* In case of reduced transmissivity, part of the leaving flux is stored as dead dislocation density.
|
|
!* That means for an interface of zero transmissivity the leaving flux is fully converted to dead dislocations.
|
|
|
|
considerLeavingFlux = .true.
|
|
if (opposite_el > 0 .and. opposite_ip > 0) then
|
|
if (phase_localConstitution(material_phase(1,opposite_ip,opposite_el))) &
|
|
considerLeavingFlux = .false.
|
|
endif
|
|
|
|
if (considerLeavingFlux) then
|
|
normal_me2neighbor_defConf = math_det3x3(Favg) * math_mul33x3(math_inv3x3(transpose(Favg)), mesh_ipAreaNormal(1:3,n,ip,el)) ! calculate the normal of the interface in (average) deformed configuration (pointing from me to my neighbor!!!)
|
|
normal_me2neighbor = math_mul33x3(transpose(my_Fe), normal_me2neighbor_defConf) / math_det3x3(my_Fe) ! interface normal in my lattice configuration
|
|
area = mesh_ipArea(n,ip,el) * math_norm3(normal_me2neighbor)
|
|
normal_me2neighbor = normal_me2neighbor / math_norm3(normal_me2neighbor) ! normalize the surface normal to unit length
|
|
do s = 1,ns
|
|
do t = 1,4
|
|
c = (t + 1) / 2
|
|
if (fluxdensity(s,t) * math_mul3x3(m(1:3,s,t), normal_me2neighbor) > 0.0_pReal ) then ! flux from me to my neighbor == leaving flux for me (might also be a pure flux from my mobile density to dead density if interface not at all transmissive)
|
|
lineLength = fluxdensity(s,t) * math_mul3x3(m(1:3,s,t), normal_me2neighbor) * area ! positive line length that wants to leave through this interface
|
|
transmissivity = sum(constitutive_nonlocal_compatibility(c,1:ns,s,n,ip,el)**2.0_pReal) ! overall transmissivity from this slip system to my neighbor
|
|
rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / mesh_ipVolume(ip,el) ! subtract dislocation flux from current mobile type
|
|
rhoDotFlux(s,t+4) = rhoDotFlux(s,t+4) + lineLength / mesh_ipVolume(ip,el) * (1.0_pReal - transmissivity) &
|
|
* sign(1.0_pReal, fluxdensity(s,t)) ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point
|
|
endif
|
|
enddo
|
|
enddo
|
|
endif
|
|
|
|
|
|
!* FLUX FROM MY NEIGHBOR TO ME
|
|
!* This is only considered, if I have a neighbor of nonlocal constitution 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.
|
|
if (neighboring_el > 0_pInt .or. neighboring_ip > 0_pInt) then
|
|
if (.not. phase_localConstitution(material_phase(1,neighboring_ip,neighboring_el)) &
|
|
.and. any(constitutive_nonlocal_compatibility(:,:,:,n,ip,el) > 0.0_pReal)) &
|
|
considerEnteringFlux = .true.
|
|
endif
|
|
|
|
if (considerEnteringFlux) then
|
|
forall (t = 1:4) &
|
|
neighboring_fluxdensity(1:ns,t) = state(g,neighboring_ip,neighboring_el)%p((t-1)*ns+1:t*ns) &
|
|
* constitutive_nonlocal_v(1:ns,t,g,neighboring_ip,neighboring_el)
|
|
normal_neighbor2me_defConf = math_det3x3(Favg) &
|
|
* math_mul33x3(math_inv3x3(transpose(Favg)), mesh_ipAreaNormal(1:3,neighboring_n,neighboring_ip,neighboring_el)) ! calculate the normal of the interface in (average) deformed configuration (now pointing from my neighbor to me!!!)
|
|
normal_neighbor2me = math_mul33x3(transpose(neighboring_Fe), normal_neighbor2me_defConf) / math_det3x3(neighboring_Fe) ! interface normal in the lattice configuration of my neighbor
|
|
area = mesh_ipArea(neighboring_n,neighboring_ip,neighboring_el) * math_norm3(normal_neighbor2me)
|
|
normal_neighbor2me = normal_neighbor2me / math_norm3(normal_neighbor2me) ! normalize the surface normal to unit length
|
|
do s = 1,ns
|
|
do t = 1,4
|
|
c = (t + 1) / 2
|
|
topp = t + mod(t,2) - mod(t+1,2)
|
|
if (neighboring_fluxdensity(s,t) * math_mul3x3(m(1:3,s,t), normal_neighbor2me) > 0.0_pReal) then ! flux from my neighbor to me == entering flux for me
|
|
lineLength = neighboring_fluxdensity(s,t) * math_mul3x3(m(1:3,s,t), normal_neighbor2me) * area ! positive line length that wants to enter through this interface
|
|
where (constitutive_nonlocal_compatibility(c,1:ns,s,n,ip,el) > 0.0_pReal) & ! positive compatibility...
|
|
rhoDotFlux(1:ns,t) = rhoDotFlux(1:ns,t) + lineLength / mesh_ipVolume(ip,el) & ! ... transferring to equally signed dislocation type
|
|
* constitutive_nonlocal_compatibility(c,1:ns,s,n,ip,el) ** 2.0_pReal
|
|
where (constitutive_nonlocal_compatibility(c,1:ns,s,n,ip,el) < 0.0_pReal) & ! ..negative compatibility...
|
|
rhoDotFlux(1:ns,topp) = rhoDotFlux(1:ns,topp) + lineLength / mesh_ipVolume(ip,el) & ! ... transferring to opposite signed dislocation type
|
|
* constitutive_nonlocal_compatibility(c,1:ns,s,n,ip,el) ** 2.0_pReal
|
|
endif
|
|
enddo
|
|
enddo
|
|
endif
|
|
|
|
enddo ! neighbor loop
|
|
|
|
constitutive_nonlocal_rhoDotFlux(1:ns,1:10,g,ip,el) = constitutive_nonlocal_rhoDotFlux(1:ns,1:10,g,ip,el) + rhoDotFlux
|
|
|
|
endif
|
|
|
|
|
|
!****************************************************************************
|
|
!*** calculate dipole formation and annihilation
|
|
|
|
!*** formation by glide
|
|
|
|
do c = 1,2
|
|
|
|
rhoDotSingle2DipoleGlide(1:ns,2*c-1) = -2.0_pReal * dUpper(1:ns,c) / constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) &
|
|
* (rhoSgl(1:ns,2*c-1) * abs(gdot(1:ns,2*c)) & ! negative mobile --> positive mobile
|
|
+ rhoSgl(1:ns,2*c) * abs(gdot(1:ns,2*c-1)) & ! positive mobile --> negative mobile
|
|
+ abs(rhoSgl(1:ns,2*c+4)) * abs(gdot(1:ns,2*c-1))) ! positive mobile --> negative immobile
|
|
|
|
rhoDotSingle2DipoleGlide(1:ns,2*c) = -2.0_pReal * dUpper(1:ns,c) / constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) &
|
|
* (rhoSgl(1:ns,2*c-1) * abs(gdot(1:ns,2*c)) & ! negative mobile --> positive mobile
|
|
+ rhoSgl(1:ns,2*c) * abs(gdot(1:ns,2*c-1)) & ! positive mobile --> negative mobile
|
|
+ abs(rhoSgl(1:ns,2*c+3)) * abs(gdot(1:ns,2*c))) ! negative mobile --> positive immobile
|
|
|
|
rhoDotSingle2DipoleGlide(1:ns,2*c+3) = -2.0_pReal * dUpper(1:ns,c) / constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) &
|
|
* rhoSgl(1:ns,2*c+3) * abs(gdot(1:ns,2*c)) ! negative mobile --> positive immobile
|
|
|
|
rhoDotSingle2DipoleGlide(1:ns,2*c+4) = -2.0_pReal * dUpper(1:ns,c) / constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) &
|
|
* rhoSgl(1:ns,2*c+4) * abs(gdot(1:ns,2*c-1)) ! positive mobile --> negative immobile
|
|
|
|
rhoDotSingle2DipoleGlide(1:ns,c+8) = - rhoDotSingle2DipoleGlide(1:ns,2*c-1) - rhoDotSingle2DipoleGlide(1:ns,2*c) &
|
|
+ abs(rhoDotSingle2DipoleGlide(1:ns,2*c+3)) + abs(rhoDotSingle2DipoleGlide(1:ns,2*c+4))
|
|
enddo
|
|
|
|
|
|
!*** athermal annihilation
|
|
|
|
rhoDotAthermalAnnihilation = 0.0_pReal
|
|
|
|
forall (c=1:2) &
|
|
rhoDotAthermalAnnihilation(1:ns,c+8) = -2.0_pReal * dLower(1:ns,c) / constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) &
|
|
* ( 2.0_pReal * (rhoSgl(1:ns,2*c-1) * abs(gdot(1:ns,2*c)) + rhoSgl(1:ns,2*c) * abs(gdot(1:ns,2*c-1))) & ! was single hitting single
|
|
+ 2.0_pReal * (abs(rhoSgl(1:ns,2*c+3)) * abs(gdot(1:ns,2*c)) + abs(rhoSgl(1:ns,2*c+4)) * abs(gdot(1:ns,2*c-1))) & ! was single hitting immobile single or was immobile single hit by single
|
|
+ rhoDip(1:ns,c) * (abs(gdot(1:ns,2*c-1)) + abs(gdot(1:ns,2*c)))) ! single knocks dipole constituent
|
|
|
|
|
|
!*** thermally activated annihilation of dipoles
|
|
|
|
rhoDotThermalAnnihilation = 0.0_pReal
|
|
|
|
D = constitutive_nonlocal_Dsd0(myInstance) * exp(-constitutive_nonlocal_Qsd(myInstance) / (kB * Temperature))
|
|
|
|
vClimb = constitutive_nonlocal_atomicVolume(myInstance) * D / ( kB * Temperature ) &
|
|
* constitutive_nonlocal_Gmod(myInstance) / ( 2.0_pReal * pi * (1.0_pReal-constitutive_nonlocal_nu(myInstance)) ) &
|
|
* 2.0_pReal / ( dUpper(1:ns,1) + dLower(1:ns,1) )
|
|
|
|
rhoDotThermalAnnihilation(1:ns,9) = - 4.0_pReal * rhoDip(1:ns,1) * vClimb / (dUpper(1:ns,1) - dLower(1:ns,1)) ! edge climb
|
|
rhoDotThermalAnnihilation(1:ns,10) = 0.0_pReal !!! cross slipping still has to be implemented !!!
|
|
|
|
|
|
!*** formation/dissociation by stress change = alteration in dUpper
|
|
|
|
rhoDotDipole2SingleStressChange = 0.0_pReal
|
|
forall (c=1:2, s=1:ns, dUpperDot(s,c) < 0.0_pReal) & ! increased stress => dipole dissociation
|
|
rhoDotDipole2SingleStressChange(s,8+c) = rhoDip(s,c) * dUpperDot(s,c) / (previousDUpper(s,c) - dLower(s,c))
|
|
|
|
forall (t=1:4) &
|
|
rhoDotDipole2SingleStressChange(1:ns,t) = -0.5_pReal * rhoDotDipole2SingleStressChange(1:ns,(t-1)/2+9)
|
|
|
|
|
|
rhoDotSingle2DipoleStressChange = 0.0_pReal
|
|
do c = 1,2
|
|
do s = 1,ns
|
|
if (dUpperDot(s,c) > 0.0_pReal) then ! stress decrease => dipole formation
|
|
rhoDotSingle2DipoleStressChange(s,2*(c-1)+1) = -4.0_pReal * dUpperDot(s,c) * previousDUpper(s,c) * rhoSgl(s,2*(c-1)+1) &
|
|
* ( rhoSgl(s,2*(c-1)+2) + abs(rhoSgl(s,2*(c-1)+6)) )
|
|
rhoDotSingle2DipoleStressChange(s,2*(c-1)+2) = -4.0_pReal * dUpperDot(s,c) * previousDUpper(s,c) * rhoSgl(s,2*(c-1)+2) &
|
|
* ( rhoSgl(s,2*(c-1)+1) + abs(rhoSgl(s,2*(c-1)+5)) )
|
|
rhoDotSingle2DipoleStressChange(s,2*(c-1)+5) = -4.0_pReal * dUpperDot(s,c) * previousDUpper(s,c) * rhoSgl(s,2*(c-1)+5) &
|
|
* ( rhoSgl(s,2*(c-1)+2) + abs(rhoSgl(s,2*(c-1)+6)) )
|
|
rhoDotSingle2DipoleStressChange(s,2*(c-1)+6) = -4.0_pReal * dUpperDot(s,c) * previousDUpper(s,c) * rhoSgl(s,2*(c-1)+6) &
|
|
* ( rhoSgl(s,2*(c-1)+1) + abs(rhoSgl(s,2*(c-1)+5)) )
|
|
endif
|
|
enddo
|
|
enddo
|
|
|
|
forall (c = 1:2) &
|
|
rhoDotSingle2DipoleStressChange(1:ns,8+c) = abs(rhoDotSingle2DipoleStressChange(1:ns,2*(c-1)+1)) &
|
|
+ abs(rhoDotSingle2DipoleStressChange(1:ns,2*(c-1)+2)) &
|
|
+ abs(rhoDotSingle2DipoleStressChange(1:ns,2*(c-1)+5)) &
|
|
+ abs(rhoDotSingle2DipoleStressChange(1:ns,2*(c-1)+6))
|
|
|
|
|
|
!****************************************************************************
|
|
!*** assign the rates of dislocation densities to my dotState
|
|
|
|
rhoDot = 0.0_pReal
|
|
forall (t = 1:10) &
|
|
rhoDot(1:ns,t) = rhoDotFlux(1:ns,t) &
|
|
+ rhoDotMultiplication(1:ns,t) &
|
|
+ rhoDotRemobilization(1:ns,t) &
|
|
+ rhoDotSingle2DipoleGlide(1:ns,t) &
|
|
+ rhoDotAthermalAnnihilation(1:ns,t) &
|
|
+ rhoDotThermalAnnihilation(1:ns,t)
|
|
! + rhoDotDipole2SingleStressChange(1:ns,t)
|
|
! + rhoDotSingle2DipoleStressChange(1:ns,t)
|
|
|
|
if (verboseDebugger .and. (debug_g==g .and. debug_i==ip .and. debug_e==el)) then
|
|
!$OMP CRITICAL (write2out)
|
|
write(6,'(a,/,8(12(e12.5,x),/))') 'dislocation remobilization', rhoDotRemobilization(1:ns,1:8) * timestep
|
|
write(6,'(a,/,4(12(e12.5,x),/))') 'dislocation multiplication', rhoDotMultiplication(1:ns,1:4) * timestep
|
|
write(6,'(a,/,8(12(e12.5,x),/))') 'dislocation flux', rhoDotFlux(1:ns,1:8) * timestep
|
|
write(6,'(a,/,10(12(e12.5,x),/))') 'dipole formation by glide', rhoDotSingle2DipoleGlide * timestep
|
|
write(6,'(a,/,2(12(e12.5,x),/))') 'athermal dipole annihilation', rhoDotAthermalAnnihilation(1:ns,1:2) * timestep
|
|
write(6,'(a,/,2(12(e12.5,x),/))') 'thermally activated dipole annihilation', rhoDotThermalAnnihilation(1:ns,9:10) * timestep
|
|
! write(6,'(a,/,10(12(e12.5,x),/))') 'dipole dissociation by stress increase', rhoDotDipole2SingleStressChange * timestep
|
|
! write(6,'(a,/,10(12(e12.5,x),/))') 'dipole formation by stress decrease', rhoDotSingle2DipoleStressChange * timestep
|
|
write(6,'(a,/,10(12(e12.5,x),/))') 'total density change', rhoDot * timestep
|
|
write(6,'(a,/,10(12(f12.7,x),/))') 'relative density change', rhoDot(1:ns,1:8) * timestep / (abs(rhoSgl)+1.0e-10), &
|
|
rhoDot(1:ns,9:10) * timestep / (rhoDip+1.0e-10)
|
|
write(6,*)
|
|
!$OMP END CRITICAL (write2out)
|
|
endif
|
|
|
|
!$OMP CRITICAL (copy2dotState)
|
|
dotState(g,ip,el)%p(1:10*ns) = dotState(g,ip,el)%p(1:10*ns) + reshape(rhoDot,(/10*ns/))
|
|
!$OMP END CRITICAL (copy2dotState)
|
|
|
|
endsubroutine
|
|
|
|
|
|
|
|
!*********************************************************************
|
|
!* COMPATIBILITY UPDATE *
|
|
!* Compatibility is defined as normalized product of signed cosine *
|
|
!* of the angle between the slip plane normals and signed cosine of *
|
|
!* the angle between the slip directions. Only the largest values *
|
|
!* that sum up to a total of 1 are considered, all others are set to *
|
|
!* zero. *
|
|
!*********************************************************************
|
|
subroutine constitutive_nonlocal_updateCompatibility(orientation,i,e)
|
|
|
|
use prec, only: pReal, &
|
|
pInt
|
|
use math, only: math_QuaternionDisorientation, &
|
|
math_mul3x3, &
|
|
math_qRot
|
|
use material, only: material_phase, &
|
|
phase_constitution, &
|
|
phase_localConstitution, &
|
|
phase_constitutionInstance, &
|
|
homogenization_maxNgrains
|
|
use mesh, only: mesh_element, &
|
|
mesh_ipNeighborhood, &
|
|
FE_NipNeighbors, &
|
|
FE_maxNipNeighbors, &
|
|
mesh_maxNips, &
|
|
mesh_NcpElems
|
|
use lattice, only: lattice_sn, &
|
|
lattice_sd, &
|
|
lattice_st
|
|
use debug, only: debugger, &
|
|
debug_e, debug_i, debug_g, &
|
|
verboseDebugger
|
|
|
|
implicit none
|
|
|
|
!* input variables
|
|
integer(pInt), intent(in) :: i, & ! ip index
|
|
e ! element index
|
|
real(pReal), dimension(4,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
|
|
orientation ! crystal orientation in quaternions
|
|
|
|
!* output variables
|
|
|
|
!* local variables
|
|
integer(pInt) Nneighbors, & ! number of neighbors
|
|
n, & ! neighbor index
|
|
neighboring_e, & ! element index of my neighbor
|
|
neighboring_i, & ! integration point index of my neighbor
|
|
my_phase, &
|
|
neighboring_phase, &
|
|
my_structure, & ! lattice structure
|
|
my_instance, & ! instance of constitution
|
|
ns, & ! number of active slip systems
|
|
s1, & ! slip system index (me)
|
|
s2 ! slip system index (my neighbor)
|
|
real(pReal), dimension(4) :: absoluteMisorientation ! absolute misorientation (without symmetry) between me and my neighbor
|
|
real(pReal), dimension(2,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(1,i,e)))) :: &
|
|
compatibility ! compatibility of one specific slip system to all neighbors slip systems's for edges and screws
|
|
real(pReal), dimension(3,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(1,i,e)))) :: &
|
|
slipNormal, &
|
|
slipDirection
|
|
real(pReal) compatibilitySum, &
|
|
thresholdValue, &
|
|
nThresholdValues
|
|
logical, dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(1,i,e)))) :: &
|
|
belowThreshold
|
|
|
|
|
|
Nneighbors = FE_NipNeighbors(mesh_element(2,e))
|
|
my_phase = material_phase(1,i,e)
|
|
my_instance = phase_constitutionInstance(my_phase)
|
|
my_structure = constitutive_nonlocal_structure(my_instance)
|
|
ns = constitutive_nonlocal_totalNslip(my_instance)
|
|
slipNormal(1:3,1:ns) = lattice_sn(1:3, constitutive_nonlocal_slipSystemLattice(1:ns,my_instance), my_structure)
|
|
slipDirection(1:3,1:ns) = lattice_sd(1:3, constitutive_nonlocal_slipSystemLattice(1:ns,my_instance), my_structure)
|
|
|
|
|
|
!*** start out fully compatible
|
|
|
|
constitutive_nonlocal_compatibility(:,:,:,:,i,e) = 0.0_pReal
|
|
forall(s1 = 1:maxval(constitutive_nonlocal_totalNslip)) &
|
|
constitutive_nonlocal_compatibility(1:2,s1,s1,1:Nneighbors,i,e) = 1.0_pReal
|
|
|
|
|
|
!*** Loop thrugh neighbors and check whether there is any compatibility.
|
|
!*** This is only the case for
|
|
|
|
do n = 1,Nneighbors
|
|
neighboring_e = mesh_ipNeighborhood(1,n,i,e)
|
|
neighboring_i = mesh_ipNeighborhood(2,n,i,e)
|
|
|
|
|
|
!* FREE SURFACE
|
|
!* Set surface transmissivity to the value specified in the material.config
|
|
|
|
if (neighboring_e <= 0 .or. neighboring_i <= 0) then
|
|
forall(s1 = 1:ns) &
|
|
constitutive_nonlocal_compatibility(1:2,s1,s1,n,i,e) = sqrt(constitutive_nonlocal_surfaceTransmissivity(my_instance))
|
|
cycle
|
|
endif
|
|
|
|
|
|
!* PHASE BOUNDARY
|
|
!* If we encounter a different nonlocal "cpfem" phase at the neighbor,
|
|
!* we consider this to be a real "physical" phase boundary, so fully incompatible.
|
|
!* If the neighboring "cpfem" phase has a local constitution,
|
|
!* we do not consider this to be a phase boundary, so fully compatible.
|
|
|
|
neighboring_phase = material_phase(1,neighboring_i,neighboring_e)
|
|
if (neighboring_phase /= my_phase) then
|
|
if (.not. phase_localConstitution(neighboring_phase)) then
|
|
constitutive_nonlocal_compatibility(:,:,:,n,i,e) = 0.0_pReal
|
|
endif
|
|
cycle
|
|
endif
|
|
|
|
|
|
!* GRAIN BOUNDARY ?
|
|
!* The compatibility value is defined as the product of the slip normal projection and the slip direction projection.
|
|
!* Its sign is always positive for screws, for edges it has the same sign as the slip normal projection.
|
|
!* Since the sum for each slip system can easily exceed one (which would result in a transmissivity larger than one),
|
|
!* only values above or equal to a certain threshold value are considered. This threshold value is chosen, such that
|
|
!* the number of compatible slip systems is minimized with the sum of the original compatibility values exceeding one.
|
|
!* Finally the smallest compatibility value is decreased until the sum is exactly equal to one.
|
|
!* All values below the threshold are set to zero.
|
|
|
|
absoluteMisorientation = math_QuaternionDisorientation(orientation(1:4,1,i,e), &
|
|
orientation(1:4,1,neighboring_i,neighboring_e), &
|
|
0_pInt) ! no symmetry
|
|
|
|
do s1 = 1,ns ! my slip systems
|
|
|
|
do s2 = 1,ns ! my neighbor's slip systems
|
|
compatibility(1,s2) = math_mul3x3(slipNormal(1:3,s1), math_qRot(absoluteMisorientation, slipNormal(1:3,s2))) &
|
|
* abs(math_mul3x3(slipDirection(1:3,s1), math_qRot(absoluteMisorientation, slipDirection(1:3,s2))))
|
|
compatibility(2,s2) = abs(math_mul3x3(slipNormal(1:3,s1), math_qRot(absoluteMisorientation, slipNormal(1:3,s2)))) &
|
|
* abs(math_mul3x3(slipDirection(1:3,s1), math_qRot(absoluteMisorientation, slipDirection(1:3,s2))))
|
|
enddo
|
|
|
|
compatibilitySum = 0.0_pReal
|
|
belowThreshold = .true.
|
|
do while (compatibilitySum < 1.0_pReal .and. any(belowThreshold(1:ns)))
|
|
thresholdValue = maxval(compatibility(2,1:ns), belowThreshold(1:ns)) ! screws always positive
|
|
nThresholdValues = dble(count(compatibility(2,1:ns) == thresholdValue))
|
|
where (compatibility(2,1:ns) >= thresholdValue) &
|
|
belowThreshold(1:ns) = .false.
|
|
if (compatibilitySum + thresholdValue * nThresholdValues > 1.0_pReal) &
|
|
where (abs(compatibility(1:2,1:ns)) == thresholdValue) &
|
|
compatibility(1:2,1:ns) = sign((1.0_pReal - compatibilitySum) / nThresholdValues, compatibility(1:2,1:ns))
|
|
compatibilitySum = compatibilitySum + nThresholdValues * thresholdValue
|
|
enddo
|
|
where (belowThreshold(1:ns)) compatibility(1,1:ns) = 0.0_pReal
|
|
where (belowThreshold(1:ns)) compatibility(2,1:ns) = 0.0_pReal
|
|
|
|
constitutive_nonlocal_compatibility(1:2,1:ns,s1,n,i,e) = compatibility(1:2,1:ns)
|
|
|
|
enddo ! my slip systems cycle
|
|
|
|
enddo ! neighbor cycle
|
|
|
|
endsubroutine
|
|
|
|
|
|
|
|
!*********************************************************************
|
|
!* rate of change of temperature *
|
|
!*********************************************************************
|
|
pure function constitutive_nonlocal_dotTemperature(Tstar_v,Temperature,state,g,ip,el)
|
|
|
|
use prec, only: pReal, &
|
|
pInt, &
|
|
p_vec
|
|
use mesh, only: mesh_NcpElems, &
|
|
mesh_maxNips
|
|
use material, only: homogenization_maxNgrains
|
|
implicit none
|
|
|
|
!* input variables
|
|
integer(pInt), intent(in) :: g, & ! current grain ID
|
|
ip, & ! current integration point
|
|
el ! current element
|
|
real(pReal), intent(in) :: Temperature ! temperature
|
|
real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola-Kirchhoff stress in Mandel notation
|
|
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
|
|
state ! microstructural state
|
|
|
|
!* output variables
|
|
real(pReal) constitutive_nonlocal_dotTemperature ! evolution of Temperature
|
|
|
|
!* local variables
|
|
|
|
constitutive_nonlocal_dotTemperature = 0.0_pReal
|
|
|
|
endfunction
|
|
|
|
|
|
|
|
!*********************************************************************
|
|
!* return array of constitutive results *
|
|
!*********************************************************************
|
|
function constitutive_nonlocal_postResults(Tstar_v, previousTstar_v, Fe, Fp, Temperature, disorientation, dt, dt_previous, &
|
|
state, previousState, dotState, g,ip,el)
|
|
|
|
use prec, only: pReal, &
|
|
pInt, &
|
|
p_vec
|
|
use math, only: math_norm3, &
|
|
math_mul6x6, &
|
|
math_mul3x3, &
|
|
math_mul33x3, &
|
|
math_mul33x33, &
|
|
math_inv3x3, &
|
|
math_det3x3, &
|
|
math_Mandel6to33, &
|
|
pi
|
|
use mesh, only: mesh_NcpElems, &
|
|
mesh_maxNips, &
|
|
mesh_maxNipNeighbors, &
|
|
mesh_element, &
|
|
FE_NipNeighbors, &
|
|
mesh_ipNeighborhood, &
|
|
mesh_ipVolume, &
|
|
mesh_ipArea, &
|
|
mesh_ipAreaNormal
|
|
use material, only: homogenization_maxNgrains, &
|
|
material_phase, &
|
|
phase_constitutionInstance, &
|
|
phase_Noutput
|
|
use lattice, only: lattice_Sslip, &
|
|
lattice_Sslip_v, &
|
|
lattice_sd, &
|
|
lattice_sn, &
|
|
lattice_st, &
|
|
lattice_maxNslipFamily, &
|
|
lattice_NslipSystem
|
|
implicit none
|
|
|
|
!*** input variables
|
|
integer(pInt), intent(in) :: g, & ! current grain number
|
|
ip, & ! current integration point
|
|
el ! current element number
|
|
real(pReal), intent(in) :: Temperature, & ! temperature
|
|
dt, & ! time increment
|
|
dt_previous ! time increment between previous and current state
|
|
real(pReal), dimension(6), intent(in) :: Tstar_v, & ! current 2nd Piola-Kirchhoff stress in Mandel notation
|
|
previousTstar_v ! previous 2nd Piola-Kirchhoff stress in Mandel notation
|
|
real(pReal), dimension(4,mesh_maxNipNeighbors), intent(in) :: &
|
|
disorientation ! crystal disorientation between me and my neighbor (axis, angle pair)
|
|
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
|
|
previousState, & ! previous microstructural state
|
|
dotState ! evolution rate of microstructural state
|
|
|
|
!*** output variables
|
|
real(pReal), dimension(constitutive_nonlocal_sizePostResults(phase_constitutionInstance(material_phase(g,ip,el)))) :: &
|
|
constitutive_nonlocal_postResults
|
|
|
|
!*** local variables
|
|
integer(pInt) myInstance, & ! current instance of this constitution
|
|
myStructure, & ! current lattice structure
|
|
ns, & ! short notation for the total number of active slip systems
|
|
neighboring_el, & ! element number of my neighbor
|
|
neighboring_ip, & ! integration point of my neighbor
|
|
c, & ! character of dislocation
|
|
cs, & ! constitutive result index
|
|
n, & ! index of my current neighbor
|
|
o, & ! index of current output
|
|
t, & ! type of dislocation
|
|
s, & ! index of my current slip system
|
|
sLattice ! index of my current slip system according to lattice order
|
|
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),6,4) :: &
|
|
fluxes ! outgoing fluxes per slipsystem, neighbor and dislocation type
|
|
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),8) :: &
|
|
rhoSgl, & ! current single dislocation densities (positive/negative screw and edge without dipoles)
|
|
previousRhoSgl, & ! previous single dislocation densities (positive/negative screw and edge without dipoles)
|
|
rhoDotSgl ! evolution rate of single dislocation densities (positive/negative screw and edge without dipoles)
|
|
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: &
|
|
gdot, & ! shear rates
|
|
lineLength ! dislocation line length leaving the current interface
|
|
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: &
|
|
rhoForest, & ! forest dislocation density
|
|
tauThreshold, & ! threshold shear stress
|
|
tau, & ! current resolved shear stress
|
|
previousTau, & ! previous resolved shear stress
|
|
invLambda, & ! inverse of mean free path for dislocations
|
|
vClimb ! climb velocity of edge dipoles
|
|
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),2) :: &
|
|
rhoDip, & ! current dipole dislocation densities (screw and edge dipoles)
|
|
previousRhoDip, & ! previous 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
|
|
previousDUpper, & ! previous maximum stable dipole distance for edges and screws
|
|
dUpperDot ! rate of change of the maximum stable dipole distance for edges and screws
|
|
real(pReal), dimension(3,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),2) :: &
|
|
m, & ! direction of dislocation motion for edge and screw (unit vector)
|
|
m_currentconf ! direction of dislocation motion for edge and screw (unit vector) in current configuration
|
|
real(pReal), dimension(3,3) :: F, & ! total deformation gradient
|
|
neighboring_F, & ! total deformation gradient of my neighbor
|
|
Favg ! average total deformation gradient of me and my neighbor
|
|
real(pReal), dimension(6) :: Tdislocation_v, & ! current dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress
|
|
previousTdislocation_v ! previous dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress
|
|
real(pReal), dimension(3) :: surfaceNormal, & ! surface normal in lattice configuration
|
|
surfaceNormal_currentconf ! surface normal in current configuration
|
|
real(pReal) area, & ! area of the current interface
|
|
detFe, & ! determinant of elastic defornmation gradient
|
|
D ! self diffusion
|
|
|
|
|
|
myInstance = phase_constitutionInstance(material_phase(g,ip,el))
|
|
myStructure = constitutive_nonlocal_structure(myInstance)
|
|
ns = constitutive_nonlocal_totalNslip(myInstance)
|
|
|
|
cs = 0_pInt
|
|
constitutive_nonlocal_postResults = 0.0_pReal
|
|
|
|
|
|
!* short hand notations for state variables
|
|
|
|
forall (t = 1:8) rhoSgl(1:ns,t) = state(g,ip,el)%p((t-1)*ns+1:t*ns)
|
|
forall (t = 1:8) previousRhoSgl(1:ns,t) = previousState(g,ip,el)%p((t-1)*ns+1:t*ns)
|
|
forall (c = 1:2) rhoDip(1:ns,c) = state(g,ip,el)%p((7+c)*ns+1:(8+c)*ns)
|
|
forall (c = 1:2) previousRhoDip(1:ns,c) = previousState(g,ip,el)%p((7+c)*ns+1:(8+c)*ns)
|
|
rhoForest = state(g,ip,el)%p(10*ns+1:11*ns)
|
|
tauThreshold = state(g,ip,el)%p(11*ns+1:12*ns)
|
|
Tdislocation_v = state(g,ip,el)%p(12*ns+1:12*ns+6)
|
|
previousTdislocation_v = previousState(g,ip,el)%p(12*ns+1:12*ns+6)
|
|
forall (t = 1:8) rhoDotSgl(1:ns,t) = dotState(g,ip,el)%p((t-1)*ns+1:t*ns)
|
|
forall (c = 1:2) rhoDotDip(1:ns,c) = dotState(g,ip,el)%p((7+c)*ns+1:(8+c)*ns)
|
|
|
|
|
|
!* Calculate shear rate
|
|
|
|
call constitutive_nonlocal_kinetics(Tstar_v, Temperature, state(g,ip,el), g, ip, el) ! need to calculate dislocation velocity again, because it was overwritten during stiffness calculation
|
|
|
|
do t = 1,4
|
|
do s = 1,ns
|
|
if (rhoSgl(s,t+4) * constitutive_nonlocal_v(s,t,g,ip,el) < 0.0_pReal) then
|
|
rhoSgl(s,t) = rhoSgl(s,t) + abs(rhoSgl(s,t+4)) ! remobilization of immobile singles for changing sign of v (bauschinger effect)
|
|
rhoSgl(s,t+4) = 0.0_pReal ! remobilization of immobile singles for changing sign of v (bauschinger effect)
|
|
endif
|
|
enddo
|
|
enddo
|
|
|
|
forall (t = 1:4) &
|
|
gdot(1:ns,t) = rhoSgl(1:ns,t) * constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) &
|
|
* constitutive_nonlocal_v(1:ns,t,g,ip,el)
|
|
|
|
!* calculate limits for stable dipole height and its rate of change
|
|
|
|
do s = 1,ns
|
|
sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance)
|
|
tau(s) = math_mul6x6( Tstar_v + Tdislocation_v, lattice_Sslip_v(1:6,sLattice,myStructure) )
|
|
previousTau(s) = math_mul6x6( previousTstar_v + previousTdislocation_v, lattice_Sslip_v(1:6,sLattice,myStructure) )
|
|
enddo
|
|
|
|
dLower(1:ns,1) = constitutive_nonlocal_dLowerEdgePerSlipSystem(1:ns,myInstance)
|
|
dLower(1:ns,2) = constitutive_nonlocal_dLowerScrewPerSlipSystem(1:ns,myInstance)
|
|
dUpper(1:ns,2) = min( constitutive_nonlocal_Gmod(myInstance) * constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) &
|
|
/ (8.0_pReal * pi * abs(tau)), &
|
|
1.0_pReal / sqrt(sum(abs(rhoSgl),2)+sum(rhoDip,2)) )
|
|
dUpper(1:ns,1) = dUpper(1:ns,2) / (1.0_pReal - constitutive_nonlocal_nu(myInstance))
|
|
previousDUpper(1:ns,2) = min( constitutive_nonlocal_Gmod(myInstance) * constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) &
|
|
/ (8.0_pReal * pi * abs(previousTau)), &
|
|
1.0_pReal / sqrt(sum(abs(previousRhoSgl),2) + sum(previousRhoDip,2)) )
|
|
previousDUpper(1:ns,1) = previousDUpper(1:ns,2) / (1.0_pReal - constitutive_nonlocal_nu(myInstance))
|
|
|
|
if (dt_previous > 0.0_pReal) then
|
|
dUpperDot = (dUpper - previousDUpper) / dt_previous
|
|
else
|
|
dUpperDot = 0.0_pReal
|
|
endif
|
|
|
|
|
|
!*** dislocation motion
|
|
|
|
m(1:3,1:ns,1) = lattice_sd(1:3,constitutive_nonlocal_slipSystemLattice(1:ns,myInstance),myStructure)
|
|
m(1:3,1:ns,2) = -lattice_st(1:3,constitutive_nonlocal_slipSystemLattice(1:ns,myInstance),myStructure)
|
|
forall (c = 1:2, s = 1:ns) &
|
|
m_currentconf(1:3,s,c) = math_mul33x3(Fe(1:3,1:3,g,ip,el), m(1:3,s,c))
|
|
|
|
|
|
do o = 1,phase_Noutput(material_phase(g,ip,el))
|
|
select case(constitutive_nonlocal_output(o,myInstance))
|
|
|
|
case ('rho')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(abs(rhoSgl),2) + sum(rhoDip,2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(abs(rhoSgl),2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_mobile')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(abs(rhoSgl(1:ns,1:4)),2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_immobile')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(abs(rhoSgl(1:ns,5:8)),2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_dip')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(rhoDip,2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_edge')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(abs(rhoSgl(1:ns,(/1,2,5,6/))),2) + rhoDip(1:ns,1)
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_edge')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(abs(rhoSgl(1:ns,(/1,2,5,6/))),2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_edge_mobile')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(rhoSgl(1:ns,1:2),2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_edge_immobile')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(abs(rhoSgl(1:ns,5:6)),2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_edge_pos')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,1) + abs(rhoSgl(1:ns,5))
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_edge_pos_mobile')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,1)
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_edge_pos_immobile')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = abs(rhoSgl(1:ns,5))
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_edge_neg')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,2) + abs(rhoSgl(1:ns,6))
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_edge_neg_mobile')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_edge_neg_immobile')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = abs(rhoSgl(1:ns,6))
|
|
cs = cs + ns
|
|
|
|
case ('rho_dip_edge')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoDip(1:ns,1)
|
|
cs = cs + ns
|
|
|
|
case ('rho_screw')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(abs(rhoSgl(1:ns,(/3,4,7,8/))),2) + rhoDip(1:ns,2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_screw')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(abs(rhoSgl(1:ns,(/3,4,7,8/))),2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_screw_mobile')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(rhoSgl(1:ns,3:4),2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_screw_immobile')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(abs(rhoSgl(1:ns,7:8)),2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_screw_pos')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,3) + abs(rhoSgl(1:ns,7))
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_screw_pos_mobile')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,3)
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_screw_pos_immobile')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = abs(rhoSgl(1:ns,7))
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_screw_neg')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,4) + abs(rhoSgl(1:ns,8))
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_screw_neg_mobile')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,4)
|
|
cs = cs + ns
|
|
|
|
case ('rho_sgl_screw_neg_immobile')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = abs(rhoSgl(1:ns,8))
|
|
cs = cs + ns
|
|
|
|
case ('rho_dip_screw')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoDip(1:ns,2)
|
|
cs = cs + ns
|
|
|
|
case ('excess_rho')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = (rhoSgl(1:ns,1) + abs(rhoSgl(1:ns,5))) &
|
|
- (rhoSgl(1:ns,2) + abs(rhoSgl(1:ns,6))) &
|
|
+ (rhoSgl(1:ns,3) + abs(rhoSgl(1:ns,7))) &
|
|
- (rhoSgl(1:ns,4) + abs(rhoSgl(1:ns,8)))
|
|
cs = cs + ns
|
|
|
|
case ('excess_rho_edge')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = (rhoSgl(1:ns,1) + abs(rhoSgl(1:ns,5))) &
|
|
- (rhoSgl(1:ns,2) + abs(rhoSgl(1:ns,6)))
|
|
cs = cs + ns
|
|
|
|
case ('excess_rho_screw')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = (rhoSgl(1:ns,3) + abs(rhoSgl(1:ns,7))) &
|
|
- (rhoSgl(1:ns,4) + abs(rhoSgl(1:ns,8)))
|
|
cs = cs + ns
|
|
|
|
case ('rho_forest')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoForest
|
|
cs = cs + ns
|
|
|
|
case ('delta')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = 1.0_pReal / sqrt(sum(abs(rhoSgl),2) + sum(rhoDip,2))
|
|
cs = cs + ns
|
|
|
|
case ('delta_sgl')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = 1.0_pReal / sqrt(sum(abs(rhoSgl),2))
|
|
cs = cs + ns
|
|
|
|
case ('delta_dip')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = 1.0_pReal / sqrt(sum(rhoDip,2))
|
|
cs = cs + ns
|
|
|
|
case ('shearrate')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(gdot,2)
|
|
cs = cs + ns
|
|
|
|
case ('resolvedstress')
|
|
do s = 1,ns
|
|
sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance)
|
|
constitutive_nonlocal_postResults(cs+s) = math_mul6x6(Tstar_v + Tdislocation_v, lattice_Sslip_v(1:6,sLattice,myStructure))
|
|
enddo
|
|
cs = cs + ns
|
|
|
|
case ('resolvedstress_internal')
|
|
do s = 1,ns
|
|
sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance)
|
|
constitutive_nonlocal_postResults(cs+s) = math_mul6x6(Tdislocation_v, lattice_Sslip_v(1:6,sLattice,myStructure))
|
|
enddo
|
|
cs = cs + ns
|
|
|
|
case ('resolvedstress_external')
|
|
do s = 1,ns
|
|
sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance)
|
|
constitutive_nonlocal_postResults(cs+s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,sLattice,myStructure))
|
|
enddo
|
|
cs = cs + ns
|
|
|
|
case ('resistance')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = tauThreshold
|
|
cs = cs + ns
|
|
|
|
case ('rho_dot')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(rhoDotSgl,2) + sum(rhoDotDip,2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_dot_sgl')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(rhoDotSgl,2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_dot_dip')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(rhoDotDip,2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_dot_gen')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(abs(gdot),2) * sqrt(rhoForest) &
|
|
/ constitutive_nonlocal_lambda0PerSlipSystem(1:ns,myInstance) &
|
|
/ constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance)
|
|
cs = cs + ns
|
|
|
|
case ('rho_dot_gen_edge')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(abs(gdot(1:ns,3:4)),2) * sqrt(rhoForest) &
|
|
/ constitutive_nonlocal_lambda0PerSlipSystem(1:ns,myInstance) &
|
|
/ constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance)
|
|
cs = cs + ns
|
|
|
|
case ('rho_dot_gen_screw')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(abs(gdot(1:ns,1:2)),2) * sqrt(rhoForest) &
|
|
/ constitutive_nonlocal_lambda0PerSlipSystem(1:ns,myInstance) &
|
|
/ constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance)
|
|
cs = cs + ns
|
|
|
|
case ('rho_dot_sgl2dip')
|
|
do c=1,2 ! dipole formation by glide
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = constitutive_nonlocal_postResults(cs+1:cs+ns) + &
|
|
2.0_pReal * dUpper(1:ns,c) / constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) &
|
|
* ( 2.0_pReal * ( rhoSgl(1:ns,2*c-1) * abs(gdot(1:ns,2*c)) &
|
|
+ rhoSgl(1:ns,2*c) * abs(gdot(1:ns,2*c-1))) & ! was single hitting single
|
|
+ 2.0_pReal * ( abs(rhoSgl(1:ns,2*c+3)) * abs(gdot(1:ns,2*c)) &
|
|
+ abs(rhoSgl(1:ns,2*c+4)) * abs(gdot(1:ns,2*c-1)))) ! was single hitting immobile/used single
|
|
enddo
|
|
! do c=1,2
|
|
! forall (s=1:ns, dUpperDot(s,c) > 0.0_pReal) & ! dipole formation by stress decrease
|
|
! constitutive_nonlocal_postResults(cs+s) = constitutive_nonlocal_postResults(cs+s) + &
|
|
! 8.0_pReal * rhoSgl(s,2*c-1) * rhoSgl(s,2*c) * previousDUpper(s,c) * dUpperDot(s,c)
|
|
! enddo
|
|
cs = cs + ns
|
|
|
|
case ('rho_dot_dip2sgl')
|
|
do c=1,2
|
|
forall (s=1:ns, dUpperDot(s,c) < 0.0_pReal) &
|
|
constitutive_nonlocal_postResults(cs+s) = constitutive_nonlocal_postResults(cs+s) - &
|
|
rhoDip(s,c) * dUpperDot(s,c) / (previousDUpper(s,c) - dLower(s,c))
|
|
enddo
|
|
cs = cs + ns
|
|
|
|
case ('rho_dot_ann_ath')
|
|
do c=1,2
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = constitutive_nonlocal_postResults(cs+1:cs+ns) + &
|
|
2.0_pReal * dLower(1:ns,c) / constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) &
|
|
* ( 2.0_pReal * ( rhoSgl(1:ns,2*c-1) * abs(gdot(1:ns,2*c)) &
|
|
+ rhoSgl(1:ns,2*c) * abs(gdot(1:ns,2*c-1))) & ! was single hitting single
|
|
+ 2.0_pReal * ( abs(rhoSgl(1:ns,2*c+3)) * abs(gdot(1:ns,2*c)) &
|
|
+ abs(rhoSgl(1:ns,2*c+4)) * abs(gdot(1:ns,2*c-1))) & ! was single hitting immobile/used single
|
|
+ rhoDip(1:ns,c) * (abs(gdot(1:ns,2*c-1)) + abs(gdot(1:ns,2*c)))) ! single knocks dipole constituent
|
|
enddo
|
|
cs = cs + ns
|
|
|
|
case ('rho_dot_ann_the')
|
|
D = constitutive_nonlocal_Dsd0(myInstance) * exp(-constitutive_nonlocal_Qsd(myInstance) / (kB * Temperature))
|
|
|
|
vClimb = constitutive_nonlocal_atomicVolume(myInstance) * D / (kB * Temperature) &
|
|
* constitutive_nonlocal_Gmod(myInstance) / (2.0_pReal * pi * (1.0_pReal-constitutive_nonlocal_nu(myInstance))) &
|
|
* 2.0_pReal / (dUpper(1:ns,1) + dLower(1:ns,1))
|
|
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = 4.0_pReal * rhoDip(1:ns,1) * vClimb / (dUpper(1:ns,1) - dLower(1:ns,1))
|
|
! !!! cross-slip of screws missing !!!
|
|
cs = cs + ns
|
|
|
|
case ('rho_dot_flux')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(constitutive_nonlocal_rhoDotFlux(1:ns,1:4,g,ip,el),2) &
|
|
+ sum(abs(constitutive_nonlocal_rhoDotFlux(1:ns,5:8,g,ip,el)),2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_dot_flux_edge')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(constitutive_nonlocal_rhoDotFlux(1:ns,1:2,g,ip,el),2) &
|
|
+ sum(abs(constitutive_nonlocal_rhoDotFlux(1:ns,5:6,g,ip,el)),2)
|
|
cs = cs + ns
|
|
|
|
case ('rho_dot_flux_screw')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(constitutive_nonlocal_rhoDotFlux(1:ns,3:4,g,ip,el),2) &
|
|
+ sum(abs(constitutive_nonlocal_rhoDotFlux(1:ns,7:8,g,ip,el)),2)
|
|
cs = cs + ns
|
|
|
|
case ('dislocationvelocity')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = constitutive_nonlocal_v(1:ns,1,g,ip,el)
|
|
cs = cs + ns
|
|
|
|
case ('fluxdensity_edge_pos_x')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,1) * constitutive_nonlocal_v(1:ns,1,g,ip,el) &
|
|
* m_currentconf(1,1:ns,1)
|
|
cs = cs + ns
|
|
|
|
case ('fluxdensity_edge_pos_y')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,1) * constitutive_nonlocal_v(1:ns,1,g,ip,el) &
|
|
* m_currentconf(2,1:ns,1)
|
|
cs = cs + ns
|
|
|
|
case ('fluxdensity_edge_pos_z')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,1) * constitutive_nonlocal_v(1:ns,1,g,ip,el) &
|
|
* m_currentconf(3,1:ns,1)
|
|
cs = cs + ns
|
|
|
|
case ('fluxdensity_edge_neg_x')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = - rhoSgl(1:ns,2) * constitutive_nonlocal_v(1:ns,2,g,ip,el) &
|
|
* m_currentconf(1,1:ns,1)
|
|
cs = cs + ns
|
|
|
|
case ('fluxdensity_edge_neg_y')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = - rhoSgl(1:ns,2) * constitutive_nonlocal_v(1:ns,2,g,ip,el) &
|
|
* m_currentconf(2,1:ns,1)
|
|
cs = cs + ns
|
|
|
|
case ('fluxdensity_edge_neg_z')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = - rhoSgl(1:ns,2) * constitutive_nonlocal_v(1:ns,2,g,ip,el) &
|
|
* m_currentconf(3,1:ns,1)
|
|
cs = cs + ns
|
|
|
|
case ('fluxdensity_screw_pos_x')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,3) * constitutive_nonlocal_v(1:ns,3,g,ip,el) &
|
|
* m_currentconf(1,1:ns,2)
|
|
cs = cs + ns
|
|
|
|
case ('fluxdensity_screw_pos_y')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,3) * constitutive_nonlocal_v(1:ns,3,g,ip,el) &
|
|
* m_currentconf(2,1:ns,2)
|
|
cs = cs + ns
|
|
|
|
case ('fluxdensity_screw_pos_z')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoSgl(1:ns,3) * constitutive_nonlocal_v(1:ns,3,g,ip,el) &
|
|
* m_currentconf(3,1:ns,2)
|
|
cs = cs + ns
|
|
|
|
case ('fluxdensity_screw_neg_x')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = - rhoSgl(1:ns,4) * constitutive_nonlocal_v(1:ns,4,g,ip,el) &
|
|
* m_currentconf(1,1:ns,2)
|
|
cs = cs + ns
|
|
|
|
case ('fluxdensity_screw_neg_y')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = - rhoSgl(1:ns,4) * constitutive_nonlocal_v(1:ns,4,g,ip,el) &
|
|
* m_currentconf(2,1:ns,2)
|
|
cs = cs + ns
|
|
|
|
case ('fluxdensity_screw_neg_z')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = - rhoSgl(1:ns,4) * constitutive_nonlocal_v(1:ns,4,g,ip,el) &
|
|
* m_currentconf(3,1:ns,2)
|
|
cs = cs + ns
|
|
|
|
case ('d_upper_edge')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = dUpper(1:ns,1)
|
|
cs = cs + ns
|
|
|
|
case ('d_upper_screw')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = dUpper(1:ns,2)
|
|
cs = cs + ns
|
|
|
|
case ('d_upper_dot_edge')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = dUpperDot(1:ns,1)
|
|
cs = cs + ns
|
|
|
|
case ('d_upper_dot_screw')
|
|
constitutive_nonlocal_postResults(cs+1:cs+ns) = dUpperDot(1:ns,2)
|
|
cs = cs + ns
|
|
|
|
end select
|
|
enddo
|
|
|
|
endfunction
|
|
END MODULE
|