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

* internal stress calculation now considers dead dislocations always at the interface, not at the center of the IP volume; used to merge them together with "normal" dislocations for stress calculation. 

* dislocation flux is blocked if we encounter a sign change in the resolved shear stress from the central ip to the neighbor
* do not set density to zero if below certain threshold; this creates an artificial sink term
This commit is contained in:
Christoph Kords 2011-08-02 11:17:45 +00:00
parent 0373fa64e4
commit 6f859e99de
1 changed files with 181 additions and 146 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

327
code/constitutive_nonlocal.f90
View File

@ -108,6 +108,8 @@ real(pReal), dimension(:,:,:,:,:,:), allocatable :: constitutive_nonlocal_
real(pReal), dimension(:,:,:), allocatable :: constitutive_nonlocal_forestProjectionEdge, & ! matrix of forest projections of edge dislocations for each instance 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_forestProjectionScrew, & ! matrix of forest projections of screw dislocations for each instance
constitutive_nonlocal_interactionMatrixSlipSlip ! interaction matrix of the different slip systems for each instance constitutive_nonlocal_interactionMatrixSlipSlip ! interaction matrix of the different slip systems for each instance
real(pReal), dimension(:,:,:,:), allocatable :: constitutive_nonlocal_lattice2slip ! orthogonal transformation matrix from lattice coordinate system to slip coordinate system (passive rotation !!!)
CONTAINS CONTAINS
!**************************************** !****************************************
@ -130,7 +132,8 @@ subroutine constitutive_nonlocal_init(file)
use prec, only: pInt, pReal use prec, only: pInt, pReal
use math, only: math_Mandel3333to66, & use math, only: math_Mandel3333to66, &
math_Voigt66to3333, & math_Voigt66to3333, &
math_mul3x3 math_mul3x3, &
math_transpose3x3
use IO, only: IO_lc, & use IO, only: IO_lc, &
IO_getTag, & IO_getTag, &
IO_isBlank, & IO_isBlank, &
@ -479,6 +482,9 @@ constitutive_nonlocal_forestProjectionScrew = 0.0_pReal
allocate(constitutive_nonlocal_interactionMatrixSlipSlip(maxTotalNslip, maxTotalNslip, maxNinstance)) allocate(constitutive_nonlocal_interactionMatrixSlipSlip(maxTotalNslip, maxTotalNslip, maxNinstance))
constitutive_nonlocal_interactionMatrixSlipSlip = 0.0_pReal constitutive_nonlocal_interactionMatrixSlipSlip = 0.0_pReal
allocate(constitutive_nonlocal_lattice2slip(1:3, 1:3, maxTotalNslip, maxNinstance))
constitutive_nonlocal_lattice2slip = 0.0_pReal
allocate(constitutive_nonlocal_v(maxTotalNslip, 4, homogenization_maxNgrains, mesh_maxNips, mesh_NcpElems)) allocate(constitutive_nonlocal_v(maxTotalNslip, 4, homogenization_maxNgrains, mesh_maxNips, mesh_NcpElems))
constitutive_nonlocal_v = 0.0_pReal constitutive_nonlocal_v = 0.0_pReal
@ -658,6 +664,13 @@ do i = 1,maxNinstance
i) i)
enddo enddo
!*** rotation matrix from lattice configuration to slip system
constitutive_nonlocal_lattice2slip(1:3,1:3,s1,i) &
= math_transpose3x3( reshape((/ lattice_sd(1:3, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure), &
-lattice_st(1:3, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure), &
lattice_sn(1:3, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure)/), (/3,3/)))
enddo enddo
@ -879,22 +892,22 @@ integer(pInt) neighboring_el, & ! element number o
neighboring_ns, & ! total number of active slip systems at neighboring material point neighboring_ns, & ! total number of active slip systems at neighboring material point
c, & ! index of dilsocation character (edge, screw) c, & ! index of dilsocation character (edge, screw)
s, & ! slip system index s, & ! slip system index
s2, & ! slip system index according to ordering in "lattice.f90"
t, & ! index of dilsocation type (e+, e-, s+, s-, used e+, used e-, used s+, used s-) t, & ! index of dilsocation type (e+, e-, s+, s-, used e+, used e-, used s+, used s-)
dir, & dir, &
deltaX, deltaY, deltaZ, & deltaX, deltaY, deltaZ, &
side side, &
j
integer(pInt), dimension(2,3) :: periodicImages integer(pInt), dimension(2,3) :: periodicImages
real(pReal) nu, & ! poisson's ratio real(pReal) nu, & ! poisson's ratio
x, y, z, & ! coordinates of connection vector in neighboring lattice frame x, y, z, & ! coordinates of connection vector in neighboring lattice frame
xsquare, ysquare, zsquare, & ! squares of respective coordinates xsquare, ysquare, zsquare, & ! squares of respective coordinates
distance, & ! length of connection vector distance, & ! length of connection vector
segmentLength, & ! segment length of dislocations segmentLength, & ! segment length of dislocations
neighboring_Nexcess, & ! excess number of dislocation segments at neighboring material point for specific slip system and dislocation character
lambda, & lambda, &
R, Rsquare, Rcube, & R, Rsquare, Rcube, &
denominator, & denominator, &
flipSign flipSign, &
neighboring_ipVolumeSideLength
real(pReal), dimension(3) :: connection, & ! connection vector between me and my neighbor in the deformed configuration real(pReal), dimension(3) :: connection, & ! connection vector between me and my neighbor in the deformed configuration
connection_neighboringLattice, & ! connection vector between me and my neighbor in the lattice configuration of my neighbor connection_neighboringLattice, & ! connection vector between me and my neighbor in the lattice configuration of my neighbor
connection_neighboringSlip, & ! connection vector between me and my neighbor in the slip system frame of my neighbor connection_neighboringSlip, & ! connection vector between me and my neighbor in the slip system frame of my neighbor
@ -903,14 +916,15 @@ real(pReal), dimension(3) :: connection, & ! connection vecto
ipCoords, & ipCoords, &
neighboring_ipCoords neighboring_ipCoords
real(pReal), dimension(3,3) :: sigma, & ! dislocation stress for one slip system in neighboring material point's slip system frame real(pReal), dimension(3,3) :: sigma, & ! dislocation stress for one slip system in neighboring material point's slip system frame
neighboringLattice2neighboringSlip, & ! orthogonal transformation matrix from lattice coordinate system to slip coordinate system (passive rotation ! ! !)
lattice2slip , & ! orthogonal transformation matrix from lattice coordinate system to slip coordinate system (passive rotation ! ! !)
Tdislo_neighboringLattice, & ! dislocation stress as 2nd Piola-Kirchhoff stress at neighboring material point Tdislo_neighboringLattice, & ! dislocation stress as 2nd Piola-Kirchhoff stress at neighboring material point
Tdislo, & ! dislocation stress as 2nd Piola-Kirchhoff stress at my material point Tdislo, & ! dislocation stress as 2nd Piola-Kirchhoff stress at my material point
invFe, & ! inverse of my elastic deformation gradient invFe, & ! inverse of my elastic deformation gradient
neighboring_invFe, &
neighboringLattice2myLattice ! mapping from neighboring MPs lattice configuration to my lattice configuration neighboringLattice2myLattice ! mapping from neighboring MPs lattice configuration to my lattice configuration
real(pReal), dimension(2,2,maxval(constitutive_nonlocal_totalNslip)) :: &
neighboring_rhoExcess ! excess density at neighboring material point (edge/screw,mobile/dead,slipsystem)
real(pReal), dimension(2,maxval(constitutive_nonlocal_totalNslip)) :: & real(pReal), dimension(2,maxval(constitutive_nonlocal_totalNslip)) :: &
neighboring_rhoExcess ! excess density at neighboring material point rhoExcessDead
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),8) :: & 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-) 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) :: & real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),2) :: &
@ -931,8 +945,6 @@ ns = constitutive_nonlocal_totalNslip(instance)
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 = 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 (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 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(instance))) &
rhoSgl(1:ns,1:4) = 0.0_pReal ! delete non-significant single density
@ -995,6 +1007,14 @@ if (.not. phase_localConstitution(phase)) then
neighboring_instance = phase_constitutionInstance(neighboring_phase) neighboring_instance = phase_constitutionInstance(neighboring_phase)
neighboring_latticeStruct = constitutive_nonlocal_structure(neighboring_instance) neighboring_latticeStruct = constitutive_nonlocal_structure(neighboring_instance)
neighboring_ns = constitutive_nonlocal_totalNslip(neighboring_instance) neighboring_ns = constitutive_nonlocal_totalNslip(neighboring_instance)
neighboring_invFe = math_inv3x3(Fe(1:3,1:3,1,neighboring_ip,neighboring_el))
neighboring_ipVolumeSideLength = mesh_ipVolume(neighboring_ip,neighboring_el) ** (1.0_pReal/3.0_pReal)
forall (s = 1:neighboring_ns, c = 1:2) &
neighboring_rhoExcess(c,1,s) = state(g,neighboring_ip,neighboring_el)%p((2*c-2)*neighboring_ns+s) & ! positive mobiles
- state(g,neighboring_ip,neighboring_el)%p((2*c-1)*neighboring_ns+s) ! negative mobiles
forall (s = 1:neighboring_ns, c = 1:2) &
neighboring_rhoExcess(c,2,s) = abs(state(g,neighboring_ip,neighboring_el)%p((2*c+2)*neighboring_ns+s)) & ! positive deads
- abs(state(g,neighboring_ip,neighboring_el)%p((2*c+3)*neighboring_ns+s)) ! negative deads
nu = constitutive_nonlocal_nu(neighboring_instance) nu = constitutive_nonlocal_nu(neighboring_instance)
Tdislo_neighboringLattice = 0.0_pReal Tdislo_neighboringLattice = 0.0_pReal
do deltaX = periodicImages(1,1),periodicImages(2,1) do deltaX = periodicImages(1,1),periodicImages(2,1)
@ -1002,35 +1022,30 @@ if (.not. phase_localConstitution(phase)) then
do deltaZ = periodicImages(1,3),periodicImages(2,3) do deltaZ = periodicImages(1,3),periodicImages(2,3)
!* special case of dead dislocations in the central ip volume !* special case of central ip volume
!* only consider dead dislocations
!* we assume that they all sit at a distance equal to half the third root of V !* we assume that they all sit at a distance equal to half the third root of V
!* the direction is determined by the character of dislocation !* in direction of the according slip direction
if (neighboring_el == el .and. neighboring_ip == ip & if (neighboring_el == el .and. neighboring_ip == ip &
.and. deltaX == 0 .and. deltaY == 0 .and. deltaZ == 0) then .and. deltaX == 0 .and. deltaY == 0 .and. deltaZ == 0) then
forall (s = 1:ns, c = 1:2) & forall (s = 1:ns, c = 1:2) &
neighboring_rhoExcess(c,s) = state(g,ip,el)%p((2*c+2)*ns+s) & ! positive deads (here we use symmetry: if this has negative sign it is treated as negative density at positive position instead of positive density at negative position) rhoExcessDead(c,s) = state(g,ip,el)%p((2*c+2)*ns+s) & ! positive deads (here we use symmetry: if this has negative sign it is treated as negative density at positive position instead of positive density at negative position)
- state(g,ip,el)%p((2*c+3)*ns+s) ! negative deads (here we use symmetry: if this has negative sign it is treated as positive density at positive position instead of negative density at negative position) + state(g,ip,el)%p((2*c+3)*ns+s) ! negative deads (here we use symmetry: if this has negative sign it is treated as positive density at positive position instead of negative density at negative position)
segmentLength = mesh_ipVolume(ip,el)**(1.0_pReal/3.0_pReal)
distance = 0.5_pReal * mesh_ipVolume(ip,el)**(1.0_pReal/3.0_pReal)
do s = 1,ns do s = 1,ns
if (all(abs(neighboring_rhoExcess(:,s)) < 1.0_pReal)) then if (all(abs(rhoExcessDead(:,s)) < constitutive_nonlocal_aTolRho(instance))) then
cycle ! not significant cycle ! not significant
endif endif
sigma = 0.0_pReal ! all components except for sigma13 are zero sigma = 0.0_pReal ! all components except for sigma13 are zero
sigma(1,3) = (neighboring_rhoExcess(1,s) + neighboring_rhoExcess(2,s) * (1.0_pReal - nu)) * mesh_ipVolume(ip,el) & sigma(1,3) = - (rhoExcessDead(1,s) + rhoExcessDead(2,s) * (1.0_pReal - nu)) * neighboring_ipVolumeSideLength &
/ (distance * sqrt(distance**2.0_pReal + 0.25_pReal * segmentLength**2.0_pReal)) & * constitutive_nonlocal_Gmod(instance) * constitutive_nonlocal_burgersPerSlipSystem(s,instance) &
* constitutive_nonlocal_Gmod(instance) * constitutive_nonlocal_burgersPerSlipSystem(s,instance) & / (sqrt(2.0_pReal) * pi * (1.0_pReal - nu))
/ (4.0_pReal * pi * (1.0_pReal - nu))
sigma(3,1) = sigma(1,3) sigma(3,1) = sigma(1,3)
s2 = constitutive_nonlocal_slipSystemLattice(s,instance)
lattice2slip = math_transpose3x3( reshape((/ lattice_sd(1:3, s2, latticeStruct), &
-lattice_st(1:3, s2, latticeStruct), &
lattice_sn(1:3, s2, latticeStruct)/), (/3,3/)))
Tdislo_neighboringLattice = Tdislo_neighboringLattice & Tdislo_neighboringLattice = Tdislo_neighboringLattice &
+ math_mul33x33(math_transpose3x3(lattice2slip), math_mul33x33(sigma, lattice2slip)) + math_mul33x33(math_transpose3x3(constitutive_nonlocal_lattice2slip(1:3,1:3,s,instance)), &
math_mul33x33(sigma, constitutive_nonlocal_lattice2slip(1:3,1:3,s,instance)))
enddo ! slip system loop enddo ! slip system loop
@ -1042,99 +1057,114 @@ if (.not. phase_localConstitution(phase)) then
neighboring_ipCoords = mesh_ipCenterOfGravity(1:3,neighboring_ip,neighboring_el) & neighboring_ipCoords = mesh_ipCenterOfGravity(1:3,neighboring_ip,neighboring_el) &
+ (/real(deltaX,pReal), real(deltaY,pReal), real(deltaZ,pReal)/) * meshSize + (/real(deltaX,pReal), real(deltaY,pReal), real(deltaZ,pReal)/) * meshSize
connection = neighboring_ipCoords - ipCoords connection = neighboring_ipCoords - ipCoords
distance = sqrt(sum(connection ** 2.0_pReal)) distance = sqrt(sum(connection * connection))
if (distance > constitutive_nonlocal_R(instance)) then if (distance > constitutive_nonlocal_R(instance)) then
cycle cycle
endif endif
!* determine the effective number of dislocations
!* the segment length is the minimum of the third root of the control volume and the ip distance !* the segment length is the minimum of the third root of the control volume and the ip distance
!* this ensures, that the central MP never sits on a neighboring dislocation segment !* this ensures, that the central MP never sits on a neighboring dislocation segment
connection_neighboringLattice = math_mul33x3(math_inv3x3(Fe(1:3,1:3,1,neighboring_ip,neighboring_el)), connection) connection_neighboringLattice = math_mul33x3(neighboring_invFe, connection)
forall (s = 1:neighboring_ns, c = 1:2) & segmentLength = min(neighboring_ipVolumeSideLength, distance)
neighboring_rhoExcess(c,s) = state(g,neighboring_ip,neighboring_el)%p((2*c-2)*neighboring_ns+s) & ! positive mobiles
+ abs(state(g,neighboring_ip,neighboring_el)%p((2*c+2)*neighboring_ns+s)) & ! positive deads
- state(g,neighboring_ip,neighboring_el)%p((2*c-1)*neighboring_ns+s) & ! negative mobiles
- abs(state(g,neighboring_ip,neighboring_el)%p((2*c+3)*neighboring_ns+s)) ! negative deads
segmentLength = min(mesh_ipVolume(neighboring_ip,neighboring_el)**(1.0_pReal/3.0_pReal), distance)
!* loop through all slip systems of the neighboring material point !* loop through all slip systems of the neighboring material point
!* and add up the stress contributions from egde and screw excess on these slip systems !* and add up the stress contributions from egde and screw excess on these slip systems (if significant)
do s = 1,neighboring_ns do s = 1,neighboring_ns
if (all(abs(neighboring_rhoExcess(:,s)) < 1.0_pReal)) then if (all(abs(neighboring_rhoExcess(:,:,s)) < constitutive_nonlocal_aTolRho(instance))) then
cycle ! not significant cycle ! not significant
endif endif
!* map the connection vector from the lattice into the slip system frame !* map the connection vector from the lattice into the slip system frame
s2 = constitutive_nonlocal_slipSystemLattice(s,neighboring_instance) connection_neighboringSlip = math_mul33x3(constitutive_nonlocal_lattice2slip(1:3,1:3,s,neighboring_instance), &
neighboringLattice2neighboringSlip = math_transpose3x3( & connection_neighboringLattice)
reshape((/ lattice_sd(1:3, s2, neighboring_latticeStruct), &
-lattice_st(1:3, s2, neighboring_latticeStruct), &
lattice_sn(1:3, s2, neighboring_latticeStruct)/), (/3,3/)))
connection_neighboringSlip = math_mul33x3(neighboringLattice2neighboringSlip, connection_neighboringLattice)
x = connection_neighboringSlip(1)
y = connection_neighboringSlip(2)
z = connection_neighboringSlip(3)
xsquare = x ** 2.0_pReal
ysquare = y ** 2.0_pReal
zsquare = z ** 2.0_pReal
!* edge contribution to stress !* edge contribution to stress
sigma = 0.0_pReal sigma = 0.0_pReal
neighboring_Nexcess = neighboring_rhoExcess(1,s) * mesh_ipVolume(neighboring_ip,neighboring_el) / segmentLength
flipSign = sign(1.0_pReal, -y)
do side = 1,-1,-2
lambda = real(side,pReal) * 0.5_pReal * segmentLength - y
R = sqrt(xsquare + zsquare + lambda**2.0_pReal)
Rsquare = R ** 2.0_pReal
Rcube = R**3.0_pReal
denominator = R * (R + flipSign * lambda)
if (denominator == 0.0_pReal) then
call IO_error(237,el,ip,g)
endif
sigma(1,1) = sigma(1,1) - real(side,pReal) * flipSign * z / denominator &
* (1.0_pReal + xsquare / Rsquare + xsquare / denominator) &
* neighboring_Nexcess
sigma(2,2) = sigma(2,2) - real(side,pReal) * (flipSign * 2.0_pReal * nu * z / denominator + z * lambda / Rcube) &
* neighboring_Nexcess
sigma(3,3) = sigma(3,3) + real(side,pReal) * flipSign * z / denominator &
* (1.0_pReal - zsquare / Rsquare - zsquare / denominator) &
* neighboring_Nexcess
sigma(1,2) = sigma(1,2) + real(side,pReal) * x * z / Rcube * neighboring_Nexcess
sigma(1,3) = sigma(1,3) + real(side,pReal) * flipSign * x / denominator &
* (1.0_pReal - zsquare / Rsquare - zsquare / denominator) &
* neighboring_Nexcess
sigma(2,3) = sigma(2,3) - real(side,pReal) * (nu / R - zsquare / Rcube) * neighboring_Nexcess
enddo
x = connection_neighboringSlip(1)
y = connection_neighboringSlip(2)
z = connection_neighboringSlip(3)
xsquare = x * x
ysquare = y * y
zsquare = z * z
do j = 1,2
if (abs(neighboring_rhoExcess(1,j,s)) < constitutive_nonlocal_aTolRho(instance)) then
cycle
elseif (j > 1) then
x = connection_neighboringSlip(1) + sign(0.5_pReal * segmentLength, &
state(g,neighboring_ip,neighboring_el)%p(4*neighboring_ns+s) &
- state(g,neighboring_ip,neighboring_el)%p(5*neighboring_ns+s))
xsquare = x * x
endif
flipSign = sign(1.0_pReal, -y)
do side = 1,-1,-2
lambda = real(side,pReal) * 0.5_pReal * segmentLength - y
R = sqrt(xsquare + zsquare + lambda * lambda)
Rsquare = R * R
Rcube = Rsquare * R
denominator = R * (R + flipSign * lambda)
if (denominator == 0.0_pReal) then
call IO_error(237,el,ip,g)
endif
sigma(1,1) = sigma(1,1) - real(side,pReal) * flipSign * z / denominator &
* (1.0_pReal + xsquare / Rsquare + xsquare / denominator) &
* neighboring_rhoExcess(1,j,s)
sigma(2,2) = sigma(2,2) - real(side,pReal) * (flipSign * 2.0_pReal * nu * z / denominator + z * lambda / Rcube)&
* neighboring_rhoExcess(1,j,s)
sigma(3,3) = sigma(3,3) + real(side,pReal) * flipSign * z / denominator &
* (1.0_pReal - zsquare / Rsquare - zsquare / denominator) &
* neighboring_rhoExcess(1,j,s)
sigma(1,2) = sigma(1,2) + real(side,pReal) * x * z / Rcube * neighboring_rhoExcess(1,j,s)
sigma(1,3) = sigma(1,3) + real(side,pReal) * flipSign * x / denominator &
* (1.0_pReal - zsquare / Rsquare - zsquare / denominator) &
* neighboring_rhoExcess(1,j,s)
sigma(2,3) = sigma(2,3) - real(side,pReal) * (nu / R - zsquare / Rcube) * neighboring_rhoExcess(1,j,s)
enddo
enddo
!* screw contribution to stress !* screw contribution to stress
neighboring_Nexcess = neighboring_rhoExcess(2,s) * mesh_ipVolume(neighboring_ip,neighboring_el) / segmentLength x = connection_neighboringSlip(1) ! have to restore this value, because position might have been adapted for edge deads before
flipSign = sign(1.0_pReal, x) do j = 1,2
do side = 1,-1,-2 if (abs(neighboring_rhoExcess(2,j,s)) < constitutive_nonlocal_aTolRho(instance)) then
lambda = x + real(side,pReal) * 0.5_pReal * segmentLength cycle
R = sqrt(ysquare + zsquare + lambda**2.0_pReal) elseif (j > 1) then
Rsquare = R ** 2.0_pReal y = connection_neighboringSlip(2) + sign(0.5_pReal * segmentLength, &
Rcube = R**3.0_pReal state(g,neighboring_ip,neighboring_el)%p(6*neighboring_ns+s) &
denominator = R * (R + flipSign * lambda) - state(g,neighboring_ip,neighboring_el)%p(7*neighboring_ns+s))
if (denominator == 0.0_pReal) then ysquare = y * y
call IO_error(237,el,ip,g)
endif endif
sigma(1,2) = sigma(1,2) - real(side,pReal) * flipSign * z * (1.0_pReal - nu) / denominator * neighboring_Nexcess flipSign = sign(1.0_pReal, x)
sigma(1,3) = sigma(1,3) + real(side,pReal) * flipSign * y * (1.0_pReal - nu) / denominator * neighboring_Nexcess do side = 1,-1,-2
lambda = x + real(side,pReal) * 0.5_pReal * segmentLength
R = sqrt(ysquare + zsquare + lambda * lambda)
Rsquare = R * R
Rcube = Rsquare * R
denominator = R * (R + flipSign * lambda)
if (denominator == 0.0_pReal) then
call IO_error(237,el,ip,g)
endif
sigma(1,2) = sigma(1,2) - real(side,pReal) * flipSign * z * (1.0_pReal - nu) / denominator &
* neighboring_rhoExcess(2,j,s)
sigma(1,3) = sigma(1,3) + real(side,pReal) * flipSign * y * (1.0_pReal - nu) / denominator &
* neighboring_rhoExcess(2,j,s)
enddo
enddo enddo
if (all(abs(sigma) < 1.0e-10_pReal)) then ! SIGMA IS NOT A REAL STRESS, THATS WHY WE NEED A REALLY SMALL VALUE HERE
cycle
endif
!* copy symmetric parts !* copy symmetric parts
@ -1147,14 +1177,15 @@ if (.not. phase_localConstitution(phase)) then
sigma = sigma * constitutive_nonlocal_Gmod(neighboring_instance) & sigma = sigma * constitutive_nonlocal_Gmod(neighboring_instance) &
* constitutive_nonlocal_burgersPerSlipSystem(s,neighboring_instance) & * constitutive_nonlocal_burgersPerSlipSystem(s,neighboring_instance) &
/ (4.0_pReal * pi * (1.0_pReal - nu)) / (4.0_pReal * pi * (1.0_pReal - nu)) &
* mesh_ipVolume(neighboring_ip,neighboring_el) / segmentLength
Tdislo_neighboringLattice = Tdislo_neighboringLattice & Tdislo_neighboringLattice = Tdislo_neighboringLattice &
+ math_mul33x33(math_transpose3x3(neighboringLattice2neighboringSlip), & + math_mul33x33(math_transpose3x3(constitutive_nonlocal_lattice2slip(1:3,1:3,s,neighboring_instance)), &
math_mul33x33(sigma, neighboringLattice2neighboringSlip)) math_mul33x33(sigma, constitutive_nonlocal_lattice2slip(1:3,1:3,s,neighboring_instance)))
enddo ! slip system loop enddo ! slip system loop
endif
endif
enddo ! deltaZ loop enddo ! deltaZ loop
enddo ! deltaY loop enddo ! deltaY loop
@ -1174,10 +1205,8 @@ if (.not. phase_localConstitution(phase)) then
endif endif
!*** set states !*** set dependent 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(10*ns+1:11*ns) = rhoForest
state(g,ip,el)%p(11*ns+1:12*ns) = tauThreshold state(g,ip,el)%p(11*ns+1:12*ns) = tauThreshold
state(g,ip,el)%p(12*ns+1:12*ns+6) = math_Mandel33to6(Tdislo) state(g,ip,el)%p(12*ns+1:12*ns+6) = math_Mandel33to6(Tdislo)
@ -1419,7 +1448,7 @@ forall (t = 1:4) &
gdot(1:ns,t) = rhoSgl(1:ns,t) * constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) & gdot(1:ns,t) = rhoSgl(1:ns,t) * constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) &
* constitutive_nonlocal_v(1:ns,t,g,ip,el) * constitutive_nonlocal_v(1:ns,t,g,ip,el)
gdotTotal = sum(gdot,2) gdotTotal = sum(gdot,2)
dgdotTotal_dtau = sum(rhoSgl,2) * constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) * dv_dtau dgdotTotal_dtau = sum(rhoSgl(1:ns,1:4),2) * constitutive_nonlocal_burgersPerSlipSystem(1:ns,myInstance) * dv_dtau
!*** Calculation of Lp and its tangent !*** Calculation of Lp and its tangent
@ -1755,6 +1784,47 @@ if (.not. phase_localConstitution(material_phase(g,ip,el))) then
endif 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.
neighboring_fluxdensity = 0.0_pReal ! needed for check of sign change in flux density below
if (neighboring_el > 0_pInt .or. neighboring_ip > 0_pInt) then
if (.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 & ! flux from my neighbor to me == entering flux for me
.and. fluxdensity(s,t) * neighboring_fluxdensity(s,t) >= 0.0_pReal ) then ! ... only if no sign change in flux density
lineLength = neighboring_fluxdensity(s,t) * math_mul3x3(m(1:3,s,t), normal_neighbor2me) * area ! positive line length that wants to enter through this interface
where (constitutive_nonlocal_compatibility(c,1:ns,s,n,ip,el) > 0.0_pReal) & ! positive compatibility...
rhoDotFlux(1:ns,t) = rhoDotFlux(1:ns,t) + lineLength / mesh_ipVolume(ip,el) & ! ... transferring to equally signed dislocation type
* constitutive_nonlocal_compatibility(c,1:ns,s,n,ip,el) ** 2.0_pReal
where (constitutive_nonlocal_compatibility(c,1:ns,s,n,ip,el) < 0.0_pReal) & ! ..negative compatibility...
rhoDotFlux(1:ns,topp) = rhoDotFlux(1:ns,topp) + lineLength / mesh_ipVolume(ip,el) & ! ... transferring to opposite signed dislocation type
* constitutive_nonlocal_compatibility(c,1:ns,s,n,ip,el) ** 2.0_pReal
endif
enddo
enddo
endif
!* FLUX FROM ME TO MY NEIGHBOR !* FLUX FROM ME TO MY NEIGHBOR
!* This is not considered, if my opposite neighbor has a local constitution. !* 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. !* Then, we assume, that the opposite(!) neighbor sends an equal amount of dislocations to me.
@ -1779,7 +1849,11 @@ if (.not. phase_localConstitution(material_phase(g,ip,el))) then
c = (t + 1) / 2 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) 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 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 if (fluxdensity(s,t) * neighboring_fluxdensity(s,t) >= 0.0_pReal) then ! no sign change in flux density
transmissivity = sum(constitutive_nonlocal_compatibility(c,1:ns,s,n,ip,el)**2.0_pReal) ! overall transmissivity from this slip system to my neighbor
else ! sign change in flux density means sign change in stress which does not allow for dislocations to arive at the neighbor
transmissivity = 0.0_pReal
endif
rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / mesh_ipVolume(ip,el) ! subtract dislocation flux from current mobile type rhoDotFlux(s,t) = 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) & 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 * 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
@ -1788,45 +1862,6 @@ if (.not. phase_localConstitution(material_phase(g,ip,el))) then
enddo enddo
endif 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 enddo ! neighbor loop
endif endif
@ -1893,13 +1928,12 @@ rhoDotThermalAnnihilation(1:ns,10) = 0.0_pReal
!*** assign the rates of dislocation densities to my dotState !*** assign the rates of dislocation densities to my dotState
rhoDot = 0.0_pReal rhoDot = 0.0_pReal
forall (t = 1:10) & rhoDot = rhoDotFlux &
rhoDot(1:ns,t) = rhoDotFlux(1:ns,t) & + rhoDotMultiplication &
+ rhoDotMultiplication(1:ns,t) & + rhoDotRemobilization &
+ rhoDotRemobilization(1:ns,t) & + rhoDotSingle2DipoleGlide &
+ rhoDotSingle2DipoleGlide(1:ns,t) & + rhoDotAthermalAnnihilation &
+ rhoDotAthermalAnnihilation(1:ns,t) & + rhoDotThermalAnnihilation
+ rhoDotThermalAnnihilation(1:ns,t)
dotState%p(1:10*ns) = dotState%p(1:10*ns) + reshape(rhoDot,(/10*ns/)) dotState%p(1:10*ns) = dotState%p(1:10*ns) + reshape(rhoDot,(/10*ns/))
@ -2033,7 +2067,8 @@ do n = 1,Nneighbors
neighboring_phase = material_phase(1,neighboring_i,neighboring_e) neighboring_phase = material_phase(1,neighboring_i,neighboring_e)
if (neighboring_phase /= my_phase) then if (neighboring_phase /= my_phase) then
if (.not. phase_localConstitution(neighboring_phase)) then if (.not. phase_localConstitution(neighboring_phase)) then
compatibility(1:2,1:ns,1:ns,n) = 0.0_pReal forall(s1 = 1:ns) &
compatibility(1:2,s1,s1,n) = 0.0_pReal ! = sqrt(0.0)
endif endif
cycle cycle
endif endif