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 in nonlocal model: instead of integration of excess gradients (->Bayley) we now sum up contributions from adjacent superdislocations within a certain radius R. When periodicity is used, also periodic images are considered in stress evaluation.

This commit is contained in:
Christoph Kords 2011-03-29 07:34:33 +00:00
parent 3d51dd36fa
commit 314ca3fe7f
3 changed files with 269 additions and 198 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

9
code/constitutive.f90
View File

@ -434,8 +434,7 @@ endfunction
!********************************************************************* !*********************************************************************
!* This function calculates from state needed variables * !* This function calculates from state needed variables *
!********************************************************************* !*********************************************************************
subroutine constitutive_microstructure(Temperature,Tstar_v,Fe,Fp,ipc,ip,el) subroutine constitutive_microstructure(Temperature, Fe, ipc, ip, el)
use prec, only: pReal,pInt use prec, only: pReal,pInt
use material, only: phase_constitution, & use material, only: phase_constitution, &
material_phase, & material_phase, &
@ -460,10 +459,8 @@ integer(pInt), intent(in):: ipc, & ! component-ID of curren
ip, & ! current integration point ip, & ! current integration point
el ! current element el ! current element
real(pReal), intent(in) :: Temperature real(pReal), intent(in) :: Temperature
real(pReal), intent(in), dimension(6) :: Tstar_v ! 2nd Piola-Kirchhoff stress
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
Fe, & ! elastic deformation gradient Fe ! elastic deformation gradient
Fp ! plastic deformation gradient
!*** output variables ***! !*** output variables ***!
@ -485,7 +482,7 @@ select case (phase_constitution(material_phase(ipc,ip,el)))
call constitutive_dislotwin_microstructure(Temperature,constitutive_state,ipc,ip,el) call constitutive_dislotwin_microstructure(Temperature,constitutive_state,ipc,ip,el)
case (constitutive_nonlocal_label) case (constitutive_nonlocal_label)
call constitutive_nonlocal_microstructure(constitutive_state, Temperature, Tstar_v, Fe, Fp, ipc, ip, el) call constitutive_nonlocal_microstructure(constitutive_state, Temperature, Fe, ipc, ip, el)
end select end select

371
code/constitutive_nonlocal.f90
View File

@ -800,19 +800,16 @@ endfunction
!********************************************************************* !*********************************************************************
!* calculates quantities characterizing the microstructure * !* calculates quantities characterizing the microstructure *
!********************************************************************* !*********************************************************************
subroutine constitutive_nonlocal_microstructure(state, Temperature, Tstar_v, Fe, Fp, g, ip, el) subroutine constitutive_nonlocal_microstructure(state, Temperature, Fe, g, ip, el)
use prec, only: pReal, & use prec, only: pReal, &
pInt, & pInt, &
p_vec p_vec
use math, only: math_Plain3333to99, & use IO, only: IO_error
math_Mandel33to6, & use math, only: math_Mandel33to6, &
math_Mandel6to33, &
math_mul33x33, & math_mul33x33, &
math_mul3x3, &
math_mul33x3, & math_mul33x3, &
math_inv3x3, & math_inv3x3, &
math_det3x3, &
math_transpose3x3, & math_transpose3x3, &
pi pi
use debug, only: debug_verbosity, & use debug, only: debug_verbosity, &
@ -822,23 +819,17 @@ use debug, only: debug_verbosity, &
debug_e debug_e
use mesh, only: mesh_NcpElems, & use mesh, only: mesh_NcpElems, &
mesh_maxNips, & mesh_maxNips, &
mesh_maxNipNeighbors, &
mesh_element, & mesh_element, &
FE_NipNeighbors, & mesh_node, &
mesh_ipNeighborhood, & FE_Nips, &
mesh_ipVolume, &
mesh_ipCenterOfGravity, & mesh_ipCenterOfGravity, &
mesh_ipAreaNormal mesh_ipVolume, &
mesh_periodicSurface
use material, only: homogenization_maxNgrains, & use material, only: homogenization_maxNgrains, &
material_phase, & material_phase, &
phase_localConstitution, & phase_localConstitution, &
phase_constitutionInstance phase_constitutionInstance
use lattice, only: lattice_Sslip, & use lattice, only: lattice_sd, &
lattice_Sslip_v, &
lattice_maxNslipFamily, &
lattice_NslipSystem, &
lattice_maxNslip, &
lattice_sd, &
lattice_sn, & lattice_sn, &
lattice_st lattice_st
@ -850,197 +841,289 @@ integer(pInt), intent(in) :: g, & ! current grain ID
el ! current element el ! current element
real(pReal), intent(in) :: Temperature ! temperature real(pReal), intent(in) :: Temperature ! temperature
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
Fe, & ! elastic deformation gradient 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
!*** input/output variables
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(inout) :: & type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(inout) :: &
state ! microstructural state state ! microstructural state
!*** output variables !*** output variables
!*** local variables !*** local variables
integer(pInt) myInstance, & ! current instance of this constitution integer(pInt) neighboring_el, & ! element number of neighboring material point
myStructure, & ! current lattice structure neighboring_ip, & ! integration point of neighboring material point
myPhase, & instance, & ! my instance of this constitution
ns, & ! short notation for the total number of active slip systems neighboring_instance, & ! instance of this constitution of neighboring material point
neighboring_el, & ! element number of my neighbor latticeStruct, & ! my lattice structure
neighboring_ip, & ! integration point of my neighbor neighboring_latticeStruct, & ! lattice structure of neighboring material point
opposite_el, & ! element number of my opposite neighbor phase, &
opposite_ip, & ! integration point of my opposite neighbor neighboring_phase, &
ns, & ! total number of active slip systems at my material point
neighboring_ns, & ! total number of active slip systems at neighboring material point
c, & ! index of dilsocation character (edge, screw) c, & ! index of dilsocation character (edge, screw)
n, & ! index of my current neighbor s, & ! slip system index
opposite_n, & ! index of my opposite neighbor s2, & ! slip system index according to ordering in "lattice.f90"
s, & ! index of my current slip system
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-)
sLattice, & ! index of my current slip system according to lattice order dir, &
i, & deltaX, deltaY, deltaZ, &
j side
real(pReal) nu ! poisson's ratio integer(pInt), dimension(2,3) :: periodicImages
real(pReal), dimension(3) :: deltaCoG ! difference of two centers of gravities (distance vector to my neighbor in reference configuration) real(pReal) nu, & ! poisson's ratio
real(pReal), dimension(3,2) :: rhoExcessDifference, & ! finite differences of excess density (in 3 directions for edge and screw) x, y, z, & ! coordinates of connection vector in neighboring lattice frame
disloGradients ! spatial gradient in excess dislocation density (in 3 directions for edge and screw) xsquare, ysquare, zsquare, & ! squares of respective coordinates
real(pReal), dimension(3,3) :: sigma, & ! dislocation stress for one slip system in its slip system frame distance, & ! length of connection vector
lattice2slip, & ! orthogonal transformation matrix from lattice coordinate system to slip coordinate system with e1=bxn, e2=b, e3=n (passive rotation!!!) segmentLength, & ! segment length of dislocations
F, & ! total deformation gradient neighboring_Nexcess, & ! excess number of dislocation segments at neighboring material point for specific slip system and dislocation character
neighboring_F, & ! total deformation gradient of neighbor lambda, &
invFe, & ! inverse elastic deformation gradient R, Rsquare, Rcube, &
Q ! inverse transpose of 3x3 matrix with finite differences of opposing position vectors denominator, &
real(pReal), dimension(6) :: Tdislocation_v ! dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress in Mandel notation flipSign
real(pReal), dimension(2,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: & real(pReal), dimension(3) :: connection, & ! connection vector between me and my neighbor in the deformed configuration
rhoExcess ! central excess density connection_neighboringLattice, & ! connection vector between me and my neighbor in the lattice configuration of my neighbor
real(pReal), dimension(6,2,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: & connection_neighboringSlip, & ! connection vector between me and my neighbor in the slip system frame of my neighbor
neighboring_rhoExcess ! excess density for each neighbor, dislo character and slip system maxCoord, minCoord, &
real(pReal), dimension(3,6) :: neighboring_position ! position vector of each neighbor when seen from the centreal material point's lattice frame meshSize, &
ipCoords, &
neighboring_ipCoords
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 ! ! !)
Tdislo_neighboringLattice, & ! dislocation stress as 2nd Piola-Kirchhoff stress at neighboring material point
Tdislo, & ! dislocation stress as 2nd Piola-Kirchhoff stress at my material point
invFe, & ! inverse of my elastic deformation gradient
neighboringLattice2myLattice ! mapping from neighboring MPs lattice configuration to my lattice configuration
real(pReal), dimension(2,maxval(constitutive_nonlocal_totalNslip)) :: &
neighboring_rhoExcess ! excess density at neighboring material point
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) :: &
rhoDip ! dipole dislocation density (edge, screw) rhoDip ! dipole dislocation density (edge, screw)
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: & real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: &
transmissivity, & ! transmissivity
rhoForest, & ! forest dislocation density rhoForest, & ! forest dislocation density
tauThreshold, & ! threshold shear stress tauThreshold, & ! threshold shear stress
tau ! resolved shear stress tau ! resolved shear stress
myPhase = material_phase(g,ip,el) phase = material_phase(g,ip,el)
myInstance = phase_constitutionInstance(myPhase) instance = phase_constitutionInstance(phase)
myStructure = constitutive_nonlocal_structure(myInstance) latticeStruct = constitutive_nonlocal_structure(instance)
ns = constitutive_nonlocal_totalNslip(myInstance) ns = constitutive_nonlocal_totalNslip(instance)
!**********************************************************************
!*** get basic states !*** 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 = 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(myInstance))) & 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 rhoSgl(1:ns,1:4) = 0.0_pReal ! delete non-significant single density
!**********************************************************************
!*** calculate dependent states
!*** calculate the forest dislocation density !*** calculate the forest dislocation density
!*** (= projection of screw and edge dislocations)
forall (s = 1:ns) & forall (s = 1:ns) &
rhoForest(s) = dot_product( ( sum(abs(rhoSgl(1:ns,(/1,2,5,6/))),2) + rhoDip(1:ns,1) ), & 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) ) & constitutive_nonlocal_forestProjectionEdge(s,1:ns,instance)) &
+ dot_product( ( sum(abs(rhoSgl(1:ns,(/3,4,7,8/))),2) + rhoDip(1:ns,2) ), & + 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 constitutive_nonlocal_forestProjectionScrew(s,1:ns,instance))
!*** calculate the threshold shear stress for dislocation slip !*** calculate the threshold shear stress for dislocation slip
forall (s = 1:ns) & forall (s = 1:ns) &
tauThreshold(s) = constitutive_nonlocal_Gmod(myInstance) & tauThreshold(s) = constitutive_nonlocal_Gmod(instance) &
* constitutive_nonlocal_burgersPerSlipSystem(s, myInstance) & * constitutive_nonlocal_burgersPerSlipSystem(s,instance) &
* sqrt( dot_product( ( sum(abs(rhoSgl),2) + sum(abs(rhoDip),2) ), & * sqrt(dot_product((sum(abs(rhoSgl),2) + sum(abs(rhoDip),2)), &
constitutive_nonlocal_interactionMatrixSlipSlip(s, 1:ns, myInstance) ) ) constitutive_nonlocal_interactionMatrixSlipSlip(s,1:ns,instance)))
!*** calculate the dislocation stress of the neighboring excess dislocation densities !*** calculate the dislocation stress of the neighboring excess dislocation densities
!*** zero for material points of local constitution
Tdislocation_v = 0.0_pReal Tdislo = 0.0_pReal
if (.not. phase_localConstitution(myPhase)) then ! only calculate dislocation stress for nonlocal material if (.not. phase_localConstitution(phase)) then
invFe = math_inv3x3(Fe(1:3,1:3,1,ip,el))
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) & !* in case of periodic surfaces we have to find out how many periodic images in each direction we need
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 dir = 1,3
do n = 1,6 maxCoord(dir) = maxval(mesh_node(dir,:))
neighboring_el = mesh_ipNeighborhood(1,n,ip,el) minCoord(dir) = minval(mesh_node(dir,:))
neighboring_ip = mesh_ipNeighborhood(2,n,ip,el) enddo
if ( neighboring_ip == 0 .or. neighboring_el == 0 ) then ! at free surfaces ... meshSize = maxCoord - minCoord
neighboring_el = el ! ... use central values instead of neighboring values ipCoords = mesh_ipCenterOfGravity(1:3,ip,el)
neighboring_ip = ip periodicImages = 0_pInt
neighboring_position(1:3,n) = 0.0_pReal do dir = 1,3
neighboring_rhoExcess(n,1:2,1:ns) = rhoExcess if (mesh_periodicSurface(dir)) then
elseif (phase_localConstitution(material_phase(1,neighboring_ip,neighboring_el))) then ! for neighbors with local constitution periodicImages(1,dir) = floor((ipCoords(dir) - constitutive_nonlocal_R(instance) - minCoord(dir)) / meshSize(dir), pInt)
neighboring_el = el ! ... use central values instead of neighboring values periodicImages(2,dir) = ceiling((ipCoords(dir) + constitutive_nonlocal_R(instance) - maxCoord(dir)) / meshSize(dir), pInt)
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))
deltaCoG = mesh_ipCenterOfGravity(1:3,neighboring_ip,neighboring_el) - mesh_ipCenterOfGravity(1:3,ip,el)
if (math_mul3x3(deltaCoG, mesh_ipAreaNormal(1:3,n,neighboring_ip,neighboring_el)) < 0.0_pReal) then ! periodic surface, so this is a twin
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)
deltaCoG = mesh_ipCenterOfGravity(1:3,ip,el) - mesh_ipCenterOfGravity(1:3,opposite_ip,opposite_el) ! assume that the twin neighbor has the same size as the opposite neighbor
endif
neighboring_position(1:3,n) = 0.5_pReal * math_mul33x3(math_mul33x33(invFe,neighboring_F) + Fp(1:3,1:3,g,ip,el), deltaCoG)
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 endif
enddo enddo
Q = math_inv3x3(math_transpose3x3(neighboring_position(1:3,(/1,3,5/)) - neighboring_position(1:3,(/2,4,6/))))
do s = 1,ns !* loop through all material points (also through their periodic images if present),
!* but only consider nonlocal neighbors within a certain cutoff radius R
lattice2slip = math_transpose3x3(reshape((/ & do neighboring_el = 1,mesh_NcpElems
lattice_sd(1:3, constitutive_nonlocal_slipSystemLattice(s,myInstance), myStructure), & do neighboring_ip = 1,FE_Nips(mesh_element(2,neighboring_el))
-lattice_st(1:3, constitutive_nonlocal_slipSystemLattice(s,myInstance), myStructure), & neighboring_phase = material_phase(g,neighboring_ip,neighboring_el)
lattice_sn(1:3, constitutive_nonlocal_slipSystemLattice(s,myInstance), myStructure)/), (/3,3/))) neighboring_instance = phase_constitutionInstance(neighboring_phase)
neighboring_latticeStruct = constitutive_nonlocal_structure(neighboring_instance)
neighboring_ns = constitutive_nonlocal_totalNslip(neighboring_instance)
nu = constitutive_nonlocal_nu(neighboring_instance)
Tdislo_neighboringLattice = 0.0_pReal
do deltaX = periodicImages(1,1),periodicImages(2,1)
do deltaY = periodicImages(1,2),periodicImages(2,2)
do deltaZ = periodicImages(1,3),periodicImages(2,3)
if (neighboring_el == el .and. neighboring_ip == ip &
.and. deltaX == 0 .and. deltaY == 0 .and. deltaZ == 0) then
cycle ! this is myself
endif
neighboring_ipCoords = mesh_ipCenterOfGravity(1:3,neighboring_ip,neighboring_el) + dble(deltaX) * meshSize(1) &
+ dble(deltaY) * meshSize(2) &
+ dble(deltaZ) * meshSize(3)
connection = neighboring_ipCoords - ipCoords
distance = sqrt(sum(connection ** 2.0_pReal))
if (.not. phase_localConstitution(neighboring_phase) &
.and. distance <= constitutive_nonlocal_R(instance)) then
rhoExcessDifference = neighboring_rhoExcess((/1,3,5/),1:2,s) - neighboring_rhoExcess((/2,4,6/),1:2,s)
forall (c = 1:2) & !* determine the effective number of dislocations
disloGradients(1:3,c) = math_mul33x3(lattice2slip, math_mul33x3(Q, rhoExcessDifference(1:3,c))) !* the segment length is the minimum of the third root of the control volume and the ip distance
!* this ensures, that the central MP never sits on a neighboring dislocation segment
connection_neighboringLattice = math_mul33x3(math_inv3x3(Fe(1:3,1:3,1,neighboring_ip,neighboring_el)), connection)
forall (s = 1:neighboring_ns, c = 1:2) &
neighboring_rhoExcess(c,s) = state(g,neighboring_ip,neighboring_el)%p((2*c-2)*neighboring_ns+s) &
+ abs(state(g,neighboring_ip,neighboring_el)%p((2*c+2)*neighboring_ns+s)) &
- state(g,neighboring_ip,neighboring_el)%p((2*c-1)*neighboring_ns+s) &
- abs(state(g,neighboring_ip,neighboring_el)%p((2*c+3)*neighboring_ns+s))
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
!* and add up the stress contributions from egde and screw excess on these slip systems
do s = 1,neighboring_ns
if (all(abs(neighboring_rhoExcess(:,s)) < 1.0_pReal)) then
cycle ! not significant
endif
!* map the connection vector from the lattice into the slip system frame
s2 = constitutive_nonlocal_slipSystemLattice(s,neighboring_instance)
neighboringLattice2neighboringSlip = math_transpose3x3( &
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
sigma = 0.0_pReal sigma = 0.0_pReal
sigma(1,1) = + 0.375_pReal / (1.0_pReal-nu) * disloGradients(3,1) neighboring_Nexcess = neighboring_rhoExcess(1,s) * mesh_ipVolume(neighboring_ip,neighboring_el) / segmentLength
sigma(2,2) = + 0.5_pReal * nu / (1.0_pReal-nu) * disloGradients(3,1) flipSign = sign(1.0_pReal, -y)
sigma(3,3) = + 0.125_pReal / (1.0_pReal-nu) * disloGradients(3,1) do side = 1,-1,-2
sigma(1,2) = + 0.25_pReal * disloGradients(3,2) lambda = dble(side) * 0.5_pReal * segmentLength - y
sigma(1,3) = - 0.125_pReal / (1.0_pReal-nu) * disloGradients(1,1) - 0.25_pReal * disloGradients(2,2) 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) - dble(side) * flipSign * z / denominator &
* (1.0_pReal + xsquare / Rsquare + xsquare / denominator) &
* neighboring_Nexcess
sigma(2,2) = sigma(2,2) - dble(side) * (flipSign * 2.0_pReal * nu * z / denominator + z * lambda / Rcube) &
* neighboring_Nexcess
sigma(3,3) = sigma(3,3) + dble(side) * flipSign * z / denominator &
* (1.0_pReal - zsquare / Rsquare - zsquare / denominator) &
* neighboring_Nexcess
sigma(1,2) = sigma(1,2) + dble(side) * x * z / Rcube * neighboring_Nexcess
sigma(1,3) = sigma(1,3) + dble(side) * flipSign * x / denominator &
* (1.0_pReal - zsquare / Rsquare - zsquare / denominator) &
* neighboring_Nexcess
sigma(2,3) = sigma(2,3) - dble(side) * (nu / R - zsquare / Rcube) * neighboring_Nexcess
enddo
!* screw contribution to stress
neighboring_Nexcess = neighboring_rhoExcess(2,s) * mesh_ipVolume(neighboring_ip,neighboring_el) / segmentLength
flipSign = sign(1.0_pReal, x)
do side = 1,-1,-2
lambda = x + dble(side) * 0.5_pReal * segmentLength
R = sqrt(ysquare + 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,2) = sigma(1,2) - dble(side) * flipSign * z * (1.0_pReal - nu) / denominator * neighboring_Nexcess
sigma(1,3) = sigma(1,3) + dble(side) * flipSign * y * (1.0_pReal - nu) / denominator * neighboring_Nexcess
enddo
!* copy symmetric parts
sigma(2,1) = sigma(1,2) sigma(2,1) = sigma(1,2)
sigma(3,1) = sigma(1,3) sigma(3,1) = sigma(1,3)
sigma(3,2) = sigma(2,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), & !* scale stresses and map them into the neighboring material point's lattice configuration
math_mul33x33(sigma, lattice2slip)))
enddo sigma = sigma * constitutive_nonlocal_Gmod(neighboring_instance) &
* constitutive_nonlocal_burgersPerSlipSystem(s,neighboring_instance) &
/ (4.0_pReal * pi * (1.0_pReal - nu))
Tdislo_neighboringLattice = Tdislo_neighboringLattice &
+ math_mul33x33(math_transpose3x3(neighboringLattice2neighboringSlip), &
math_mul33x33(sigma, neighboringLattice2neighboringSlip))
enddo ! slip system loop
endif
enddo ! deltaZ loop
enddo ! deltaY loop
enddo ! deltaX loop
!* map the stress from the neighboring MP's lattice configuration into the deformed configuration
!* and back into my lattice configuration
neighboringLattice2myLattice = math_mul33x33(invFe, Fe(1:3,1:3,1,neighboring_ip,neighboring_el))
Tdislo = Tdislo + math_mul33x33(neighboringLattice2myLattice, &
math_mul33x33(Tdislo_neighboringLattice, math_transpose3x3(neighboringLattice2myLattice)))
enddo ! ip loop
enddo ! element loop
endif endif
!**********************************************************************
!*** set states !*** set states
state(g,ip,el)%p(1:8*ns) = reshape(rhoSgl,(/8*ns/)) ! ensure positive single mobile densities 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(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) = Tdislocation_v state(g,ip,el)%p(12*ns+1:12*ns+6) = math_Mandel33to6(Tdislo)
#ifndef _OPENMP #ifndef _OPENMP
@ -1050,7 +1133,7 @@ state(g,ip,el)%p(12*ns+1:12*ns+6) = Tdislocation_v
write(6,*) write(6,*)
write(6,'(a,/,12(x),12(e10.3,x))') '<< CONST >> rhoForest', rhoForest write(6,'(a,/,12(x),12(e10.3,x))') '<< CONST >> rhoForest', rhoForest
write(6,'(a,/,12(x),12(f10.5,x))') '<< CONST >> tauThreshold / MPa', tauThreshold/1e6 write(6,'(a,/,12(x),12(f10.5,x))') '<< CONST >> tauThreshold / MPa', tauThreshold/1e6
write(6,'(a,/,3(12(x),3(f10.5,x),/))') '<< CONST >> Tdislocation / MPa', math_Mandel6to33(Tdislocation_v)/1e6 write(6,'(a,/,3(12(x),3(f10.5,x),/))') '<< CONST >> Tdislocation / MPa', Tdislo/1e6
endif endif
#endif #endif

25
code/crystallite.f90
View File

@ -209,7 +209,6 @@ allocate(crystallite_localConstitution(gMax,iMax,eMax)); crystallite_loc
allocate(crystallite_requested(gMax,iMax,eMax)); crystallite_requested = .false. allocate(crystallite_requested(gMax,iMax,eMax)); crystallite_requested = .false.
allocate(crystallite_todo(gMax,iMax,eMax)); crystallite_todo = .false. allocate(crystallite_todo(gMax,iMax,eMax)); crystallite_todo = .false.
allocate(crystallite_converged(gMax,iMax,eMax)); crystallite_converged = .true. allocate(crystallite_converged(gMax,iMax,eMax)); crystallite_converged = .true.
allocate(crystallite_output(maxval(crystallite_Noutput), & allocate(crystallite_output(maxval(crystallite_Noutput), &
material_Ncrystallite)) ; crystallite_output = '' material_Ncrystallite)) ; crystallite_output = ''
allocate(crystallite_sizePostResults(material_Ncrystallite)) ; crystallite_sizePostResults = 0_pInt allocate(crystallite_sizePostResults(material_Ncrystallite)) ; crystallite_sizePostResults = 0_pInt
@ -349,8 +348,7 @@ crystallite_orientation0 = crystallite_orientation ! Store initial o
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do g = 1,myNgrains do g = 1,myNgrains
call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar0_v(1:6,g,i,e), & call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe, g, i, e) ! update dependent state variables to be consistent with basic states
crystallite_Fe, crystallite_Fp0, g, i, e) ! update dependent state variables to be consistent with basic states
enddo enddo
enddo enddo
enddo enddo
@ -986,8 +984,7 @@ do n = 1,4
!$OMP DO !$OMP DO
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(1:6,g,i,e), & call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe, g, i, e) ! update dependent state variables to be consistent with basic states
crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMP ENDDO !$OMP ENDDO
@ -1339,8 +1336,7 @@ do n = 1,5
!$OMP DO !$OMP DO
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(1:6,g,i,e), & call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe, g, i, e) ! update dependent state variables to be consistent with basic states
crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states
endif endif
constitutive_dotState(g,i,e)%p = 0.0_pReal ! reset dotState to zero constitutive_dotState(g,i,e)%p = 0.0_pReal ! reset dotState to zero
enddo; enddo; enddo enddo; enddo; enddo
@ -1518,8 +1514,7 @@ relTemperatureResiduum = 0.0_pReal
!$OMP DO !$OMP DO
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(1:6,g,i,e), & call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe, g, i, e) ! update dependent state variables to be consistent with basic states
crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMP ENDDO !$OMP ENDDO
@ -1717,8 +1712,7 @@ stateResiduum = 0.0_pReal
!$OMP DO !$OMP DO
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(1:6,g,i,e), & call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe, g, i, e) ! update dependent state variables to be consistent with basic states
crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states
endif endif
constitutive_dotState(g,i,e)%p = 0.0_pReal ! reset dotState to zero constitutive_dotState(g,i,e)%p = 0.0_pReal ! reset dotState to zero
enddo; enddo; enddo enddo; enddo; enddo
@ -1992,8 +1986,7 @@ endif
!$OMP DO !$OMP DO
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(1:6,g,i,e), & call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe, g, i, e) ! update dependent state variables to be consistent with basic states
crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMP ENDDO !$OMP ENDDO
@ -2165,8 +2158,7 @@ endif
!$OMP DO !$OMP DO
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(1:6,g,i,e), & call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe, g, i, e) ! update dependent state variables to be consistent with basic states
crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states
endif endif
constitutive_previousDotState2(g,i,e)%p = constitutive_previousDotState(g,i,e)%p ! age previous dotState constitutive_previousDotState2(g,i,e)%p = constitutive_previousDotState(g,i,e)%p ! age previous dotState
constitutive_previousDotState(g,i,e)%p = constitutive_dotState(g,i,e)%p ! age previous dotState constitutive_previousDotState(g,i,e)%p = constitutive_dotState(g,i,e)%p ! age previous dotState
@ -2273,8 +2265,7 @@ do while (any(crystallite_todo) .and. NiterationState < nState )
!$OMP DO !$OMP DO
do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains do e = eIter(1),eIter(2); do i = iIter(1,e),iIter(2,e); do g = gIter(1,e),gIter(2,e) ! iterate over elements, ips and grains
if (crystallite_todo(g,i,e)) then if (crystallite_todo(g,i,e)) then
call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Tstar_v(1:6,g,i,e), & call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe, g, i, e) ! update dependent state variables to be consistent with basic states
crystallite_Fe, crystallite_Fp, g, i, e) ! update dependent state variables to be consistent with basic states
endif endif
constitutive_previousDotState2(g,i,e)%p = constitutive_previousDotState(g,i,e)%p ! age previous dotState constitutive_previousDotState2(g,i,e)%p = constitutive_previousDotState(g,i,e)%p ! age previous dotState
constitutive_previousDotState(g,i,e)%p = constitutive_dotState(g,i,e)%p ! age previous dotState constitutive_previousDotState(g,i,e)%p = constitutive_dotState(g,i,e)%p ! age previous dotState