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

constitutive_nonlocal 

- take orientation gradients into account when calculating dislocation stress and dislocation fluxes
- hard coded value for nu
- changed kinetics (parameter G0 is currently defined as a parameter, needs to be read from material.config)
- added some output statements

constitutive: 
- some functions and subroutines needed additional input variables for passing to constitutive_nonlocal 

crystallite:
- some functions and subroutines needed additional input variables for passing to constitutive
- call microstructure with current temperature, Fp, Fe, not "sub0" values
- show number of IPs, that are "onTrack" instead of those not "onTrack"
- calculate Fe at beginning of substep, since we need it for state preguess
This commit is contained in:
Christoph Kords 2009-10-07 15:31:52 +00:00
parent 7d6729f1c5
commit 9b3a59646a
3 changed files with 284 additions and 209 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
code/constitutive.f90
View File

@ -252,7 +252,7 @@ return
endfunction endfunction
subroutine constitutive_microstructure(Temperature,Fp,ipc,ip,el) subroutine constitutive_microstructure(Temperature,Fe,Fp,ipc,ip,el)
!********************************************************************* !*********************************************************************
!* This function calculates from state needed variables * !* This function calculates from state needed variables *
!* INPUT: * !* INPUT: *
@ -277,6 +277,7 @@ subroutine constitutive_microstructure(Temperature,Fp,ipc,ip,el)
!* Definition of variables !* Definition of variables
integer(pInt), intent(in) :: ipc,ip,el integer(pInt), intent(in) :: ipc,ip,el
real(pReal), intent(in) :: Temperature real(pReal), intent(in) :: Temperature
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: Fe
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: Fp real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: Fp
select case (phase_constitution(material_phase(ipc,ip,el))) select case (phase_constitution(material_phase(ipc,ip,el)))
@ -291,7 +292,7 @@ real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)
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(Temperature, Fp, constitutive_state,ipc,ip,el) call constitutive_nonlocal_microstructure(Temperature, Fe, Fp, constitutive_state, ipc, ip, el)
end select end select
@ -346,7 +347,7 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, Temperature, ip
endsubroutine endsubroutine
subroutine constitutive_collectDotState(Tstar_v, subTstar0_v, Fp, invFp, Temperature, subdt, ipc, ip, el) subroutine constitutive_collectDotState(Tstar_v, subTstar0_v, Fe, Fp, Temperature, subdt, ipc, ip, el)
!********************************************************************* !*********************************************************************
!* This subroutine contains the constitutive equation for * !* This subroutine contains the constitutive equation for *
!* calculating the rate of change of microstructure * !* calculating the rate of change of microstructure *
@ -359,9 +360,10 @@ subroutine constitutive_collectDotState(Tstar_v, subTstar0_v, Fp, invFp, Tempera
!* OUTPUT: * !* OUTPUT: *
!* - constitutive_dotState : evolution of state variable * !* - constitutive_dotState : evolution of state variable *
!********************************************************************* !*********************************************************************
use prec, only: pReal,pInt use prec, only: pReal, pInt
use debug use debug
use material, only: phase_constitution,material_phase use mesh, only: mesh_NcpElems, mesh_maxNips
use material, only: phase_constitution, material_phase, homogenization_maxNgrains
use constitutive_j2 use constitutive_j2
use constitutive_phenopowerlaw use constitutive_phenopowerlaw
use constitutive_dislotwin use constitutive_dislotwin
@ -369,10 +371,10 @@ subroutine constitutive_collectDotState(Tstar_v, subTstar0_v, Fp, invFp, Tempera
implicit none implicit none
!* Definition of variables !* Definition of variables
integer(pInt) ipc,ip,el integer(pInt), intent(in) :: ipc,ip,el
real(pReal) Temperature, subdt real(pReal), intent(in) :: Temperature, subdt
real(pReal), dimension(3,3) :: Fp, invFp real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: Fe, Fp
real(pReal), dimension(6) :: Tstar_v, subTstar0_v real(pReal), dimension(6), intent(in) :: Tstar_v, subTstar0_v
select case (phase_constitution(material_phase(ipc,ip,el))) select case (phase_constitution(material_phase(ipc,ip,el)))
@ -386,7 +388,7 @@ subroutine constitutive_collectDotState(Tstar_v, subTstar0_v, Fp, invFp, Tempera
constitutive_dotState(ipc,ip,el)%p = constitutive_dislotwin_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el) constitutive_dotState(ipc,ip,el)%p = constitutive_dislotwin_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_nonlocal_label) case (constitutive_nonlocal_label)
call constitutive_nonlocal_dotState(constitutive_dotState, Tstar_v, subTstar0_v, Fp, invFp, Temperature, subdt, & call constitutive_nonlocal_dotState(constitutive_dotState, Tstar_v, subTstar0_v, Fe, Fp, Temperature, subdt, &
constitutive_state, constitutive_subState0, ipc, ip, el) constitutive_state, constitutive_subState0, ipc, ip, el)
end select end select

380
code/constitutive_nonlocal.f90
View File

@ -27,7 +27,7 @@ character(len=16), dimension(3), parameter :: constitutive_nonlocal_listDependen
'tauSlipThreshold', & 'tauSlipThreshold', &
'Tdislocation_v ' /) ! list of microstructural state variables that depend on other state variables '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 real(pReal), parameter :: kB = 1.38e-23_pReal ! Physical parameter, Boltzmann constant in J/Kelvin
real(pReal), parameter :: constitutive_nonlocal_G = 2.5e-19_pReal
!* Definition of global variables !* Definition of global variables
integer(pInt), dimension(:), allocatable :: constitutive_nonlocal_sizeDotState, & ! number of dotStates integer(pInt), dimension(:), allocatable :: constitutive_nonlocal_sizeDotState, & ! number of dotStates
@ -476,45 +476,36 @@ do i = 1,maxNinstance
constitutive_nonlocal_Cslip_3333(:,:,:,:,i) = math_Voigt66to3333(constitutive_nonlocal_Cslip_66(:,:,i)) constitutive_nonlocal_Cslip_3333(:,:,:,:,i) = math_Voigt66to3333(constitutive_nonlocal_Cslip_66(:,:,i))
constitutive_nonlocal_Gmod(i) = constitutive_nonlocal_C44(i) constitutive_nonlocal_Gmod(i) = constitutive_nonlocal_C44(i)
constitutive_nonlocal_nu(i) = constitutive_nonlocal_C12(i) / constitutive_nonlocal_C11(i) constitutive_nonlocal_nu(i) = 0.3_pReal ! harcoded version, since not clear how to exactly calculate that value from C11,C12,C44 etc
! constitutive_nonlocal_nu(i) = constitutive_nonlocal_C12(i) / constitutive_nonlocal_C11(i)
do s1 = 1,ns
f = constitutive_nonlocal_slipFamily(s1,i)
!*** burgers vector, dislocation velocity prefactor, mean free path prefactor and minimum dipole distance for each slip system !*** burgers vector, dislocation velocity prefactor, mean free path prefactor and minimum dipole distance for each slip system
do s = 1,constitutive_nonlocal_totalNslip(i) constitutive_nonlocal_burgersPerSlipSystem(s1,i) = constitutive_nonlocal_burgersPerSlipFamily(f,i)
constitutive_nonlocal_v0PerSlipSystem(s1,i) = constitutive_nonlocal_v0PerSlipFamily(f,i)
constitutive_nonlocal_lambda0PerSlipSystem(s1,i) = constitutive_nonlocal_lambda0PerSlipFamily(f,i)
constitutive_nonlocal_dDipMinEdgePerSlipSystem(s1,i) = constitutive_nonlocal_dDipMinEdgePerSlipFamily(f,i)
constitutive_nonlocal_dDipMinScrewPerSlipSystem(s1,i) = constitutive_nonlocal_dDipMinScrewPerSlipFamily(f,i)
f = constitutive_nonlocal_slipFamily(s,i) do s2 = 1,ns
constitutive_nonlocal_burgersPerSlipSystem(s,i) = constitutive_nonlocal_burgersPerSlipFamily(f,i) !*** calculation of forest projections for edge and screw dislocations. s2 acts as forest to s1
constitutive_nonlocal_v0PerSlipSystem(s,i) = constitutive_nonlocal_v0PerSlipFamily(f,i)
constitutive_nonlocal_lambda0PerSlipSystem(s,i) = constitutive_nonlocal_lambda0PerSlipFamily(f,i)
constitutive_nonlocal_dDipMinEdgePerSlipSystem(s,i) = constitutive_nonlocal_dDipMinEdgePerSlipFamily(f,i)
constitutive_nonlocal_dDipMinScrewPerSlipSystem(s,i) = constitutive_nonlocal_dDipMinScrewPerSlipFamily(f,i)
enddo
!*** calculation of forest projections for edge and screw dislocations
do s1 = 1,constitutive_nonlocal_totalNslip(i)
do s2 = 1,constitutive_nonlocal_totalNslip(i)
constitutive_nonlocal_forestProjectionEdge(s1, s2, i) & constitutive_nonlocal_forestProjectionEdge(s1, s2, i) &
= abs(math_mul3x3(lattice_sn(:, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure), & = abs(math_mul3x3(lattice_sn(:, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure), &
lattice_st(:, constitutive_nonlocal_slipSystemLattice(s2,i), myStructure))) ! forest projection of edge dislocations is the projection of (b x n) onto the slip normal of the respective splip plane lattice_st(:, 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) & constitutive_nonlocal_forestProjectionScrew(s1, s2, i) &
= abs(math_mul3x3(lattice_sn(:, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure), & = abs(math_mul3x3(lattice_sn(:, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure), &
lattice_sd(:, 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 lattice_sd(:, 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
enddo; enddo
!*** calculation of interaction matrices !*** calculation of interaction matrices
do s1 = 1,constitutive_nonlocal_totalNslip(i)
do s2 = 1,constitutive_nonlocal_totalNslip(i)
constitutive_nonlocal_interactionMatrixSlipSlip(s1, s2, i) & constitutive_nonlocal_interactionMatrixSlipSlip(s1, s2, i) &
= constitutive_nonlocal_interactionSlipSlip( lattice_interactionSlipSlip(constitutive_nonlocal_slipSystemLattice(s1,i), & = constitutive_nonlocal_interactionSlipSlip( lattice_interactionSlipSlip(constitutive_nonlocal_slipSystemLattice(s1,i), &
constitutive_nonlocal_slipSystemLattice(s2,i), & constitutive_nonlocal_slipSystemLattice(s2,i), &
@ -558,31 +549,25 @@ real(pReal), dimension(constitutive_nonlocal_totalNslip(myInstance)) :: &
tauSlipThreshold ! threshold shear stress for slip tauSlipThreshold ! threshold shear stress for slip
integer(pInt) ns, & ! short notation for total number of active slip systems integer(pInt) ns, & ! short notation for total number of active slip systems
f, & ! index of lattice family f, & ! index of lattice family
s0, & from, &
s1, & upto, &
s ! index of slip system s ! index of slip system
constitutive_nonlocal_stateInit = 0.0_pReal constitutive_nonlocal_stateInit = 0.0_pReal
ns = constitutive_nonlocal_totalNslip(myInstance) ns = constitutive_nonlocal_totalNslip(myInstance)
!*** set the basic state variables !*** set the basic state variables
s1 = 0_pInt
do f = 1,lattice_maxNslipFamily do f = 1,lattice_maxNslipFamily
from = 1+sum(constitutive_nonlocal_Nslip(1:f-1,myInstance))
s0 = s1 + 1_pInt upto = sum(constitutive_nonlocal_Nslip(1:f,myInstance))
s1 = s0 + constitutive_nonlocal_Nslip(f,myInstance) - 1_pInt rhoEdgePos(from:upto) = constitutive_nonlocal_rhoEdgePos0(f, myInstance)
rhoEdgeNeg(from:upto) = constitutive_nonlocal_rhoEdgeNeg0(f, myInstance)
do s = s0,s1 rhoScrewPos(from:upto) = constitutive_nonlocal_rhoScrewPos0(f, myInstance)
rhoEdgePos(s) = constitutive_nonlocal_rhoEdgePos0(f, myInstance) rhoScrewNeg(from:upto) = constitutive_nonlocal_rhoScrewNeg0(f, myInstance)
rhoEdgeNeg(s) = constitutive_nonlocal_rhoEdgeNeg0(f, myInstance) rhoEdgeDip(from:upto) = constitutive_nonlocal_rhoEdgeDip0(f, myInstance)
rhoScrewPos(s) = constitutive_nonlocal_rhoScrewPos0(f, myInstance) rhoScrewDip(from:upto) = constitutive_nonlocal_rhoScrewDip0(f, myInstance)
rhoScrewNeg(s) = constitutive_nonlocal_rhoScrewNeg0(f, myInstance)
rhoEdgeDip(s) = constitutive_nonlocal_rhoEdgeDip0(f, myInstance)
rhoScrewDip(s) = constitutive_nonlocal_rhoScrewDip0(f, myInstance)
enddo
enddo enddo
@ -591,16 +576,16 @@ enddo
! forest dislocation density ! forest dislocation density
forall (s = 1:ns) & forall (s = 1:ns) &
rhoForest(s) & rhoForest(s) &
= dot_product( (rhoEdgePos + rhoEdgeNeg + rhoEdgeDip), constitutive_nonlocal_forestProjectionEdge(1:ns, s, myInstance) ) & = dot_product( (rhoEdgePos + rhoEdgeNeg + rhoEdgeDip), constitutive_nonlocal_forestProjectionEdge(s, 1:ns, myInstance) ) &
+ dot_product( (rhoScrewPos + rhoScrewNeg + rhoScrewDip), constitutive_nonlocal_forestProjectionScrew(1:ns, s, myInstance) ) ! calculation of forest dislocation density as projection of screw and edge dislocations + dot_product( (rhoScrewPos + rhoScrewNeg + rhoScrewDip), constitutive_nonlocal_forestProjectionScrew(s, 1:ns, myInstance) ) ! calculation of forest dislocation density as projection of screw and edge dislocations
! threshold shear stress for dislocation slip ! threshold shear stress for dislocation slip
forall (s = 1:ns) & forall (s = 1:ns) &
tauSlipThreshold(s) = constitutive_nonlocal_Gmod(myInstance) & tauSlipThreshold(s) = constitutive_nonlocal_Gmod(myInstance) &
* constitutive_nonlocal_burgersPerSlipSystem(s, myInstance) & * constitutive_nonlocal_burgersPerSlipSystem(s, myInstance) &
* sqrt( dot_product( (rhoEdgePos + rhoEdgeNeg + rhoScrewPos + rhoScrewNeg), & * sqrt( dot_product( (rhoEdgePos + rhoEdgeNeg + rhoScrewPos + rhoScrewNeg + rhoEdgeDip + rhoScrewDip), &
constitutive_nonlocal_interactionMatrixSlipSlip(1:ns, s, myInstance) ) ) constitutive_nonlocal_interactionMatrixSlipSlip(s, 1:ns, myInstance) ) )
!*** put everything together and in right order !*** put everything together and in right order
@ -680,7 +665,7 @@ endfunction
!********************************************************************* !*********************************************************************
!* calculates quantities characterizing the microstructure * !* calculates quantities characterizing the microstructure *
!********************************************************************* !*********************************************************************
subroutine constitutive_nonlocal_microstructure(Temperature, Fp, state, g, ip, el) subroutine constitutive_nonlocal_microstructure(Temperature, Fe, Fp, state, g, ip, el)
use prec, only: pReal, & use prec, only: pReal, &
pInt, & pInt, &
@ -691,6 +676,8 @@ use math, only: math_Plain3333to99, &
math_mul33x33, & math_mul33x33, &
math_mul3x3, & math_mul3x3, &
math_mul33x3, & math_mul33x3, &
math_inv3x3, &
math_det3x3, &
pi pi
use debug, only: debugger use debug, only: debugger
use mesh, only: mesh_NcpElems, & use mesh, only: mesh_NcpElems, &
@ -720,6 +707,7 @@ 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
Fp ! plastic deformation gradient Fp ! plastic deformation gradient
!*** input/output variables !*** input/output variables
@ -740,11 +728,19 @@ integer(pInt) myInstance, & ! current instance
real(pReal) gb, & ! short notation for G*b/2/pi real(pReal) gb, & ! short notation for G*b/2/pi
x, & ! coordinate in direction of lvec x, & ! coordinate in direction of lvec
y, & ! coordinate in direction of bvec y, & ! coordinate in direction of bvec
z ! coordinate in direction of nvec z, & ! coordinate in direction of nvec
detFe, & ! determinant of elastic deformation gradient
neighboring_detFe ! determinant of my neighboring elastic deformation gradient
real(pReal), dimension(3) :: connectingVector ! vector connecting the centers of gravity of me and my neigbor real(pReal), dimension(3) :: connectingVector ! vector connecting the centers of gravity of me and my neigbor
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(6) :: Tdislocation_v ! dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress in Mandel notation
real(pReal), dimension(3,3) :: transform, & ! orthogonal transformation matrix from slip coordinate system with e1=bxn, e2=b, e3=n to lattice coordinate system real(pReal), dimension(3,3) :: lattice2slip, & ! orthogonal transformation matrix from lattice coordinate system to slip coordinate system with e1=bxn, e2=b, e3=n
sigma ! Tdislocation resulting from the excess dislocation density of a single slip system and a single neighbor calculated in the coordinate system of the slip system neighboringSlip2myLattice, &! mapping from my neighbors slip coordinate system to my lattice coordinate system
sigma, & ! Tdislocation resulting from the excess dislocation density of a single slip system and a single neighbor calculated in the coordinate system of the slip system
F, & ! total deformation gradient
neighboring_F, & ! total deformation gradient of my neighbor
Favg, & ! average total deformation gradient of me and my neighbor
invFe, & ! inverse of elastic deformation gradient
neighboring_invFe ! inverse of my neighboring elastic deformation gradient
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)))) :: &
rhoEdgePos, & ! positive edge dislocation density rhoEdgePos, & ! positive edge dislocation density
rhoEdgeNeg, & ! negative edge dislocation density rhoEdgeNeg, & ! negative edge dislocation density
@ -791,12 +787,13 @@ forall (s = 1:ns) &
+ dot_product( (rhoScrewPos + rhoScrewNeg + rhoScrewDip), constitutive_nonlocal_forestProjectionScrew(1:ns, s, myInstance) ) ! calculation of forest dislocation density as projection of screw and edge dislocations + dot_product( (rhoScrewPos + rhoScrewNeg + rhoScrewDip), constitutive_nonlocal_forestProjectionScrew(1:ns, s, myInstance) ) ! calculation of forest dislocation density as projection of screw and edge dislocations
! if (debugger) write(6,'(a23,3(i3,x),/,12(e10.3,x),/)') 'forest dislocation density at ',g,ip,el, rhoForest ! if (debugger) write(6,'(a23,3(i3,x),/,12(e10.3,x),/)') 'forest dislocation density at ',g,ip,el, rhoForest
!*** calculate the threshold shear stress for dislocation slip !*** calculate the threshold shear stress for dislocation slip
forall (s = 1:ns) & forall (s = 1:ns) &
tauSlipThreshold(s) = constitutive_nonlocal_Gmod(myInstance) & tauSlipThreshold(s) = constitutive_nonlocal_Gmod(myInstance) &
* constitutive_nonlocal_burgersPerSlipSystem(s, myInstance) & * constitutive_nonlocal_burgersPerSlipSystem(s, myInstance) &
* sqrt( dot_product( (rhoEdgePos + rhoEdgeNeg + rhoScrewPos + rhoScrewNeg), & * sqrt( dot_product( (rhoEdgePos + rhoEdgeNeg + rhoEdgeDip + rhoScrewPos + rhoScrewNeg + rhoScrewDip), &
constitutive_nonlocal_interactionMatrixSlipSlip(1:ns, s, myInstance) ) ) constitutive_nonlocal_interactionMatrixSlipSlip(1:ns, s, myInstance) ) )
! if (debugger) write(6,'(a26,3(i3,x),/,12(f10.5,x),/)') 'tauSlipThreshold / MPa at ',g,ip,el, tauSlipThreshold/1e6 ! if (debugger) write(6,'(a26,3(i3,x),/,12(f10.5,x),/)') 'tauSlipThreshold / MPa at ',g,ip,el, tauSlipThreshold/1e6
@ -804,10 +801,11 @@ forall (s = 1:ns) &
!*** calculate the dislocation stress of the neighboring excess dislocation densities !*** calculate the dislocation stress of the neighboring excess dislocation densities
Tdislocation_v = 0.0_pReal Tdislocation_v = 0.0_pReal
F = math_mul33x33(Fe(:,:,g,ip,el), Fp(:,:,g,ip,el))
detFe = math_det3x3(Fe)
invFe = math_inv3x3(Fe)
! loop through my neighbors (if existent!) ! loop through my neighbors (if existent!)
do n = 1,FE_NipNeighbors(mesh_element(2,el)) do n = 1,FE_NipNeighbors(mesh_element(2,el))
neighboring_el = mesh_ipNeighborhood(1,n,ip,el) neighboring_el = mesh_ipNeighborhood(1,n,ip,el)
@ -815,17 +813,21 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el))
if ( neighboring_el == 0 .or. neighboring_ip == 0 ) cycle if ( neighboring_el == 0 .or. neighboring_ip == 0 ) cycle
! deformation gradients needed for mapping between configurations
neighboring_F = math_mul33x33(Fe(:,:,g,neighboring_ip,neighboring_el), Fp(:,:,g,neighboring_ip,neighboring_el))
Favg = 0.5_pReal * (F + neighboring_F)
neighboring_detFe = math_det3x3(Fe(:,:,g,neighboring_ip,neighboring_el))
neighboring_invFe = math_inv3x3(Fe(:,:,g,neighboring_ip,neighboring_el))
! calculate the connecting vector between me and my neighbor and his excess dislocation density ! calculate connection vector between me and my neighbor in its lattice configuration
connectingVector = math_mul33x3(neighboring_invFe, math_mul33x3(Favg, &
connectingVector = math_mul33x3( Fp(:,:,g,neighboring_ip,neighboring_el), & (mesh_ipCenterOfGravity(:,neighboring_ip,neighboring_el) - mesh_ipCenterOfGravity(:,ip,el)) ) )
(mesh_ipCenterOfGravity(:,ip,el) - mesh_ipCenterOfGravity(:,neighboring_ip,neighboring_el)) )
! neighboring dislocation densities
neighboring_rhoEdgePos = state(1, neighboring_ip, neighboring_el)%p( 1: ns) neighboring_rhoEdgePos = state(1, neighboring_ip, neighboring_el)%p( 1: ns)
neighboring_rhoEdgeNeg = state(1, neighboring_ip, neighboring_el)%p( ns+1:2*ns) neighboring_rhoEdgeNeg = state(1, neighboring_ip, neighboring_el)%p( ns+1:2*ns)
neighboring_rhoScrewPos = state(1, neighboring_ip, neighboring_el)%p(2*ns+1:3*ns) neighboring_rhoScrewPos = state(1, neighboring_ip, neighboring_el)%p(2*ns+1:3*ns)
neighboring_rhoScrewNeg = state(1, neighboring_ip, neighboring_el)%p(3*ns+1:4*ns) neighboring_rhoScrewNeg = state(1, neighboring_ip, neighboring_el)%p(3*ns+1:4*ns)
neighboring_rhoEdgeExcess = neighboring_rhoEdgePos - neighboring_rhoEdgeNeg neighboring_rhoEdgeExcess = neighboring_rhoEdgePos - neighboring_rhoEdgeNeg
neighboring_rhoScrewExcess = neighboring_rhoScrewPos - neighboring_rhoScrewNeg neighboring_rhoScrewExcess = neighboring_rhoScrewPos - neighboring_rhoScrewNeg
@ -833,19 +835,16 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el))
neighboring_Nscrew = neighboring_rhoScrewExcess * mesh_ipVolume(neighboring_ip, neighboring_el) ** (2.0_pReal/3.0_pReal) neighboring_Nscrew = neighboring_rhoScrewExcess * mesh_ipVolume(neighboring_ip, neighboring_el) ** (2.0_pReal/3.0_pReal)
! loop over slip systems and get their slip directions, slip normals, and sd x sn ! loop over slip systems and get their slip directions, slip normals, and sd x sn
do s = 1,ns do s = 1,ns
transform = reshape( (/lattice_st(:, constitutive_nonlocal_slipSystemLattice(s,myInstance), myStructure), & lattice2slip = reshape( (/lattice_st(:, constitutive_nonlocal_slipSystemLattice(s,myInstance), myStructure), &
lattice_sd(:, constitutive_nonlocal_slipSystemLattice(s,myInstance), myStructure), & lattice_sd(:, constitutive_nonlocal_slipSystemLattice(s,myInstance), myStructure), &
lattice_sn(:, constitutive_nonlocal_slipSystemLattice(s,myInstance), myStructure)/), (/3,3/) ) lattice_sn(:, constitutive_nonlocal_slipSystemLattice(s,myInstance), myStructure)/), (/3,3/) )
! coordinate transformation of connecting vector from the lattice coordinate system to the slip coordinate system
! coordinate transformation of p from the lattice coordinate system to the slip coordinate system x = math_mul3x3(lattice2slip(1,:), -connectingVector)
x = math_mul3x3(connectingVector, transform(:,1)) y = math_mul3x3(lattice2slip(2,:), -connectingVector)
y = math_mul3x3(connectingVector, transform(:,2)) z = math_mul3x3(lattice2slip(3,:), -connectingVector)
z = math_mul3x3(connectingVector, transform(:,3))
! calculate the back stress in the slip coordinate system for this slip system ! calculate the back stress in the slip coordinate system for this slip system
gb = constitutive_nonlocal_Gmod(myInstance) * constitutive_nonlocal_burgersPerSlipSystem(s,myInstance) / (2.0_pReal*pi) gb = constitutive_nonlocal_Gmod(myInstance) * constitutive_nonlocal_burgersPerSlipSystem(s,myInstance) / (2.0_pReal*pi)
@ -861,7 +860,7 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el))
sigma(1,2) = gb * neighboring_Nscrew(s) * z / (x**2.0_pReal + z**2.0_pReal) sigma(1,2) = gb * neighboring_Nscrew(s) * z / (x**2.0_pReal + z**2.0_pReal)
sigma(2,3) = gb * ( neighboring_Nedge(s) / (1.0_pReal-constitutive_nonlocal_nu(myInstance)) & sigma(2,3) = gb * ( neighboring_Nedge(s) / (1.0_pReal-constitutive_nonlocal_nu(myInstance)) &
* (y**2.0_pReal - z**2.0_pReal) / (y**2.0_pReal + z**2.0_pReal) & * y * (y**2.0_pReal - z**2.0_pReal) / (y**2.0_pReal + z**2.0_pReal)**2.0_pReal &
- neighboring_Nscrew(s) * x / (x**2.0_pReal + z**2.0_pReal) ) - neighboring_Nscrew(s) * x / (x**2.0_pReal + z**2.0_pReal) )
sigma(2,1) = sigma(1,2) sigma(2,1) = sigma(1,2)
@ -870,11 +869,19 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el))
sigma(3,1) = 0.0_pReal sigma(3,1) = 0.0_pReal
! coordinate transformation from the slip coordinate system to the lattice coordinate system ! coordinate transformation from the slip coordinate system to the lattice coordinate system
Tdislocation_v = Tdislocation_v + math_Mandel33to6( math_mul33x33(transform, math_mul33x33(sigma, transpose(transform)) ) ) neighboringSlip2myLattice = math_mul33x33(invFe,math_mul33x33(Fe(:,:,g,neighboring_ip,neighboring_el),transpose(lattice2slip)))
! if (debugger) write(6,'(a15,3(i3,x),/,3(3(f12.3,x)/))') 'sigma / MPa at ',g,ip,el, sigma/1e6 Tdislocation_v = Tdislocation_v + math_Mandel33to6( detFe / neighboring_detFe &
! if (debugger) write(6,'(a15,3(i3,x),/,3(3(f12.3,x)/))') 'Tdislocation / MPa at ',g,ip,el, math_Mandel6to33(Tdislocation_v/1e6) * math_mul33x33(neighboringSlip2myLattice, math_mul33x33(sigma, transpose(neighboringSlip2myLattice))) )
enddo enddo
enddo enddo
! if (debugger) then
! !$OMP CRITICAL (write2out)
! write(6,*) '::: constitutive_nonlocal_microstructure at ',g,ip,el
! write(6,*)
! write(6,'(a,/,6(f12.5,x),/)') 'Tdislocation_v / MPa', Tdislocation_v/1e6_pReal
! !$OMP CRITICAL (write2out)
! endif
!********************************************************************** !**********************************************************************
@ -965,9 +972,6 @@ forall (t = 1:4) rho(:,t) = state(g,ip,el)%p((t-1)*ns+1:t*ns)
rhoForest = state(g,ip,el)%p(6*ns+1:7*ns) rhoForest = state(g,ip,el)%p(6*ns+1:7*ns)
tauSlipThreshold = state(g,ip,el)%p(7*ns+1:8*ns) tauSlipThreshold = state(g,ip,el)%p(7*ns+1:8*ns)
Tdislocation_v = state(g,ip,el)%p(8*ns+1:8*ns+6) Tdislocation_v = state(g,ip,el)%p(8*ns+1:8*ns+6)
! if (debugger) write(6,'(a20,3(i3,x),/,3(3(f12.3,x)/))') 'Tdislocation / MPa at ', g,ip,el, math_Mandel6to33(Tdislocation_v/1e6)
! if (debugger) write(6,'(a15,3(i3,x),/,3(3(f12.3,x)/))') 'Tstar / MPa at ',g,ip,el, math_Mandel6to33(Tstar_v/1e6)
!*** calculation of resolved stress !*** calculation of resolved stress
@ -975,19 +979,15 @@ forall (s =1:ns) &
tauSlip(s) = math_mul6x6( Tstar_v + Tdislocation_v, & tauSlip(s) = math_mul6x6( Tstar_v + Tdislocation_v, &
lattice_Sslip_v(:,constitutive_nonlocal_slipSystemLattice(s,myInstance),myStructure) ) lattice_Sslip_v(:,constitutive_nonlocal_slipSystemLattice(s,myInstance),myStructure) )
!*** Calculation of gdot and its tangent !*** Calculation of gdot and its tangent
v = constitutive_nonlocal_v0PerSlipSystem(:,myInstance) & v = constitutive_nonlocal_v0PerSlipSystem(:,myInstance) &
* exp( - ( tauSlipThreshold - abs(tauSlip) ) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance)**2.0_pReal & * exp( - constitutive_nonlocal_G / ( kB * Temperature) * (1.0_pReal - (abs(tauSlip)/tauSlipThreshold) ) ) &
/ ( kB * Temperature * sqrt(rhoForest) ) ) &
* sign(1.0_pReal,tauSlip) * sign(1.0_pReal,tauSlip)
gdotSlip = sum(rho,2) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) * v gdotSlip = sum(rho,2) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) * v
dgdot_dtauSlip = gdotSlip * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance)**2.0_pReal & dgdot_dtauSlip = gdotSlip * constitutive_nonlocal_G / ( kB * Temperature * tauSlipThreshold )
/ ( kB * Temperature * sqrt(rhoForest) )
!*** Calculation of Lp and its tangent !*** Calculation of Lp and its tangent
@ -996,7 +996,6 @@ do s = 1,ns
sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance) sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance)
Lp = Lp + gdotSlip(s) * lattice_Sslip(:,:,sLattice,myStructure) Lp = Lp + gdotSlip(s) * lattice_Sslip(:,:,sLattice,myStructure)
! if (debugger) write(6,'(a4,i2,a3,/,3(3(f15.7)/))') 'dLp(',s,'): ',gdotSlip(s) * lattice_Sslip(:,:,sLattice,myStructure)
forall (i=1:3,j=1:3,k=1:3,l=1:3) & 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) + dgdot_dtauSlip(s) * lattice_Sslip(i,j, sLattice,myStructure) & dLp_dTstar3333(i,j,k,l) = dLp_dTstar3333(i,j,k,l) + dgdot_dtauSlip(s) * lattice_Sslip(i,j, sLattice,myStructure) &
@ -1009,10 +1008,11 @@ dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333)
! !$OMP CRITICAL (write2out) ! !$OMP CRITICAL (write2out)
! write(6,*) '::: LpandItsTangent',g,ip,el ! write(6,*) '::: LpandItsTangent',g,ip,el
! write(6,*) ! write(6,*)
! write(6,'(a,/,12(f12.5,x))') 'gdot/1e-3',gdotSlip*1e3_pReal ! write(6,'(a,/,12(f12.5,x),/)') 'tauSlip / MPa', tauSlip/1e6_pReal
! write(6,*) ! write(6,'(a,/,12(f12.5,x),/)') 'tauSlipThreshold / MPa', tauSlipThreshold/1e6_pReal
! write(6,'(a,/,12(f12.5,x),/)') 'v', v
! write(6,'(a,/,12(f12.5,x),/)') 'gdot total /1e-3',gdotSlip*1e3_pReal
! write(6,'(a,/,3(3(f12.7,x)/))') 'Lp',Lp ! write(6,'(a,/,3(3(f12.7,x)/))') 'Lp',Lp
! write(6,*)
! !$OMPEND CRITICAL (write2out) ! !$OMPEND CRITICAL (write2out)
! endif ! endif
@ -1023,7 +1023,7 @@ endsubroutine
!********************************************************************* !*********************************************************************
!* rate of change of microstructure * !* rate of change of microstructure *
!********************************************************************* !*********************************************************************
subroutine constitutive_nonlocal_dotState(dotState, Tstar_v, subTstar0_v, Fp, invFp, Temperature, subdt, state, subState0, g,ip,el) subroutine constitutive_nonlocal_dotState(dotState, Tstar_v, subTstar0_v, Fe, Fp, Temperature, subdt, state, subState0, g,ip,el)
use prec, only: pReal, & use prec, only: pReal, &
pInt, & pInt, &
@ -1033,7 +1033,9 @@ use math, only: math_norm3, &
math_mul6x6, & math_mul6x6, &
math_mul3x3, & math_mul3x3, &
math_mul33x3, & math_mul33x3, &
math_transpose3x3, & math_mul33x33, &
math_inv3x3, &
math_det3x3, &
pi pi
use mesh, only: mesh_NcpElems, & use mesh, only: mesh_NcpElems, &
mesh_maxNips, & mesh_maxNips, &
@ -1063,8 +1065,9 @@ real(pReal), intent(in) :: Temperature, & ! temperat
subdt ! substepped crystallite time increment subdt ! substepped crystallite time increment
real(pReal), dimension(6), intent(in) :: Tstar_v, & ! current 2nd Piola-Kirchhoff stress in Mandel notation real(pReal), dimension(6), intent(in) :: Tstar_v, & ! current 2nd Piola-Kirchhoff stress in Mandel notation
subTstar0_v ! 2nd Piola-Kirchhoff stress in Mandel notation at start of crystallite increment subTstar0_v ! 2nd Piola-Kirchhoff stress in Mandel notation at start of crystallite increment
real(pReal), dimension(3,3), intent(in) :: Fp, & ! plastic deformation gradient real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
invFp ! inverse of plastic deformation gradient Fe, & ! elastic deformation gradient
Fp ! plastic deformation gradient
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
state, & ! current microstructural state state, & ! current microstructural state
subState0 ! microstructural state at start of crystallite increment subState0 ! microstructural state at start of crystallite increment
@ -1087,7 +1090,8 @@ integer(pInt) myInstance, & ! current
s ! index of my current slip system s ! index of my current slip system
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: & real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: &
rho, & ! dislocation densities (positive/negative screw and edge without dipoles) rho, & ! dislocation densities (positive/negative screw and edge without dipoles)
rhoDot, & ! rate of change of dislocation densities totalRhoDot, & ! total rate of change of dislocation densities
thisRhoDot, & ! rate of change of dislocation densities for this mechanism
gdot, & ! shear rates gdot, & ! shear rates
lineLength ! dislocation line length leaving the current interface lineLength ! dislocation line length leaving the current interface
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)))) :: &
@ -1100,22 +1104,25 @@ real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstan
vClimb ! climb velocity of edge dipoles vClimb ! climb velocity of edge dipoles
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 densities (screw and edge dipoles) rhoDip, & ! dipole dislocation densities (screw and edge dipoles)
rhoDipDot, & ! rate of change of dipole dislocation densities totalRhoDipDot, & ! total rate of change of dipole dislocation densities
rhoDotTransfer, & ! dislocation density rate that is transferred from single dislocation to dipole dislocation thisRhoDipDot, & ! rate of change of dipole dislocation densities for this mechanism
dDipMin, & ! minimum stable dipole distance for edges and screws dDipMin, & ! minimum stable dipole distance for edges and screws
dDipMax, & ! current maximum stable dipole distance for edges and screws dDipMax, & ! current maximum stable dipole distance for edges and screws
dDipMax0, & ! maximum stable dipole distance for edges and screws at start of crystallite increment dDipMax0, & ! maximum stable dipole distance for edges and screws at start of crystallite increment
dDipMaxDot ! rate of change of the maximum stable dipole distance for edges and screws dDipMaxDot ! 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) :: & real(pReal), dimension(3,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: &
m ! direction of dislocation motion m ! direction of dislocation motion
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 real(pReal), dimension(6) :: Tdislocation_v, & ! current dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress
subTdislocation0_v ! dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress at start of crystallite increment subTdislocation0_v ! dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress at start of crystallite increment
real(pReal), dimension(3) :: surfaceNormal ! surface normal of the current interface real(pReal), dimension(3) :: surfaceNormal ! surface normal of the current interface
real(pReal) norm_surfaceNormal, & ! euclidic norm of the surface normal real(pReal) norm_surfaceNormal, & ! euclidic norm of the surface normal
area, & ! area of the current interface area, & ! area of the current interface
detFe, & ! determinant of elastic defornmation gradient
D ! self diffusion D ! self diffusion
myInstance = phase_constitutionInstance(material_phase(g,ip,el)) myInstance = phase_constitutionInstance(material_phase(g,ip,el))
myStructure = constitutive_nonlocal_structure(myInstance) myStructure = constitutive_nonlocal_structure(myInstance)
ns = constitutive_nonlocal_totalNslip(myInstance) ns = constitutive_nonlocal_totalNslip(myInstance)
@ -1128,9 +1135,10 @@ dDipMin = 0.0_pReal
dDipMax = 0.0_pReal dDipMax = 0.0_pReal
dDipMax0 = 0.0_pReal dDipMax0 = 0.0_pReal
dDipMaxDot = 0.0_pReal dDipMaxDot = 0.0_pReal
rhoDot = 0.0_pReal totalRhoDot = 0.0_pReal
rhoDipDot = 0.0_pReal thisRhoDot = 0.0_pReal
rhoDotTransfer = 0.0_pReal totalRhoDipDot = 0.0_pReal
thisRhoDipDot = 0.0_pReal
!*** shortcut to state variables !*** shortcut to state variables
@ -1150,70 +1158,94 @@ do s = 1,ns ! loop over slip systems
tauSlip(s) = math_mul6x6( Tstar_v + Tdislocation_v, & tauSlip(s) = math_mul6x6( Tstar_v + Tdislocation_v, &
lattice_Sslip_v(:,constitutive_nonlocal_slipSystemLattice(s,myInstance),myStructure) ) lattice_Sslip_v(:,constitutive_nonlocal_slipSystemLattice(s,myInstance),myStructure) )
subTauSlip0(s) = math_mul6x6( subTstar0_v + subTdislocation0_v, & subTauSlip0(s) = math_mul6x6( subTstar0_v - subTdislocation0_v, &
lattice_Sslip_v(:,constitutive_nonlocal_slipSystemLattice(s,myInstance),myStructure) ) lattice_Sslip_v(:,constitutive_nonlocal_slipSystemLattice(s,myInstance),myStructure) )
enddo enddo
v = constitutive_nonlocal_v0PerSlipSystem(:,myInstance) & v = constitutive_nonlocal_v0PerSlipSystem(:,myInstance) &
* exp( - ( tauSlipThreshold - abs(tauSlip) ) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance)**2.0_pReal & * exp( - constitutive_nonlocal_G / ( kB * Temperature) * (1.0_pReal - (abs(tauSlip)/tauSlipThreshold) ) ) &
/ ( kB * Temperature * sqrt(rhoForest) ) ) &
* sign(1.0_pReal,tauSlip) * sign(1.0_pReal,tauSlip)
forall (t = 1:4) & forall (t = 1:4) &
gdot(:,t) = rho(:,t) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) * v gdot(:,t) = rho(:,t) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) * v
!****************************************************************************
!*** calculate dislocation multiplication
invLambda = sqrt(rhoForest) / constitutive_nonlocal_lambda0PerSlipSystem(:,myInstance)
rhoDot = rhoDot + spread(0.25_pReal * sum(abs(gdot),2) * invLambda / constitutive_nonlocal_burgersPerSlipSystem(:,myInstance), 2, 4)
if (debugger) then if (debugger) then
!$OMP CRITICAL (write2out)
write(6,*) '::: constitutive_nonlocal_dotState at ',g,ip,el write(6,*) '::: constitutive_nonlocal_dotState at ',g,ip,el
write(6,*) write(6,*)
write(6,'(a,/,12(f12.5,x),/)') 'tauSlip / MPa', tauSlip/1e6_pReal write(6,'(a,/,12(f12.5,x),/)') 'tauSlip / MPa', tauSlip/1e6_pReal
write(6,'(a,/,12(f12.5,x),/)') 'tauSlipThreshold / MPa', tauSlipThreshold/1e6_pReal write(6,'(a,/,12(f12.5,x),/)') 'tauSlipThreshold / MPa', tauSlipThreshold/1e6_pReal
write(6,'(a,/,12(e10.3,x),/)') 'v', v write(6,'(a,/,12(e12.3,x),/)') 'v', v
write(6,'(a,/,4(12(f12.5,x),/))') 'gdot / 1e-3', gdot*1e3_pReal write(6,'(a,/,4(12(f12.5,x),/))') 'gdot / 1e-3', gdot*1e3_pReal
write(6,'(a,/,(12(f12.5,x),/))') 'gdot total/ 1e-3', sum(gdot,2)*1e3_pReal write(6,'(a,/,(12(f12.5,x),/))') 'gdot total/ 1e-3', sum(gdot,2)*1e3_pReal
write(6,'(a,/,6(12(e12.5,x),/))') 'dislocation multiplication', & !$OMP CRITICAL (write2out)
spread(0.25_pReal * sum(abs(gdot),2) * invLambda / constitutive_nonlocal_burgersPerSlipSystem(:,myInstance), 2, 4)*subdt, & endif
0.0_pReal*rhoDotTransfer
!****************************************************************************
!*** calculate dislocation multiplication
thisRhoDot = 0.0_pReal
thisRhoDipDot = 0.0_pReal
invLambda = sqrt(rhoForest) / constitutive_nonlocal_lambda0PerSlipSystem(:,myInstance)
thisRhoDot = spread(0.25_pReal * sum(abs(gdot),2) * invLambda / constitutive_nonlocal_burgersPerSlipSystem(:,myInstance), 2, 4)
totalRhoDot = totalRhoDot + thisRhoDot
totalRhoDipDot = totalRhoDipDot + thisRhoDipDot
if (debugger) then
!$OMP CRITICAL (write2out)
write(6,'(a,/,6(12(e12.5,x),/))') 'dislocation multiplication', thisRhoDot * subdt, thisRhoDipDot * subdt
!$OMP CRITICAL (write2out)
endif endif
!**************************************************************************** !****************************************************************************
!*** calculate dislocation fluxes !*** calculate dislocation fluxes
thisRhoDot = 0.0_pReal
thisRhoDipDot = 0.0_pReal
m(:,:,1) = lattice_sd(:, constitutive_nonlocal_slipSystemLattice(:,myInstance), myStructure) m(:,:,1) = lattice_sd(:, constitutive_nonlocal_slipSystemLattice(:,myInstance), myStructure)
m(:,:,2) = -lattice_sd(:, constitutive_nonlocal_slipSystemLattice(:,myInstance), myStructure) m(:,:,2) = -lattice_sd(:, constitutive_nonlocal_slipSystemLattice(:,myInstance), myStructure)
m(:,:,3) = lattice_st(:, constitutive_nonlocal_slipSystemLattice(:,myInstance), myStructure) m(:,:,3) = lattice_st(:, constitutive_nonlocal_slipSystemLattice(:,myInstance), myStructure)
m(:,:,4) = -lattice_st(:, constitutive_nonlocal_slipSystemLattice(:,myInstance), myStructure) m(:,:,4) = -lattice_st(:, constitutive_nonlocal_slipSystemLattice(:,myInstance), myStructure)
F = math_mul33x33(Fe(:,:,g,ip,el), Fp(:,:,g,ip,el))
detFe = math_det3x3(Fe(:,:,g,ip,el))
do n = 1,FE_NipNeighbors(mesh_element(2,el)) ! loop through my neighbors do n = 1,FE_NipNeighbors(mesh_element(2,el)) ! loop through my neighbors
neighboring_el = mesh_ipNeighborhood(1,n,ip,el) neighboring_el = mesh_ipNeighborhood(1,n,ip,el)
neighboring_ip = mesh_ipNeighborhood(2,n,ip,el) neighboring_ip = mesh_ipNeighborhood(2,n,ip,el)
! calculate the area and the surface normal of the interface ! if neighbor exists, total deformation gradient is averaged over me and my neighbor
surfaceNormal = math_mul33x3(math_transpose3x3(invFp), mesh_ipAreaNormal(:,n,ip,el)) if ( neighboring_el > 0 .and. neighboring_ip > 0 ) then
neighboring_F = math_mul33x33(Fe(:,:,g,neighboring_ip,neighboring_el), Fp(:,:,g,neighboring_ip,neighboring_el))
Favg = 0.5_pReal * (F + neighboring_F)
else
Favg = F
endif
! calculate the area and the surface normal (of unit length) of the interface in lattice configuration
surfaceNormal = math_det3x3(Favg) / detFe &
* math_mul33x3(transpose(Fe(:,:,g,ip,el)), math_mul33x3(Favg,mesh_ipAreaNormal(:,n,ip,el)))
norm_surfaceNormal = math_norm3(surfaceNormal) norm_surfaceNormal = math_norm3(surfaceNormal)
surfaceNormal = surfaceNormal / norm_surfaceNormal surfaceNormal = surfaceNormal / norm_surfaceNormal
area = mesh_ipArea(n,ip,el) / norm_surfaceNormal area = mesh_ipArea(n,ip,el) * norm_surfaceNormal
lineLength = 0.0_pReal lineLength = 0.0_pReal
do s = 1,ns ! loop over slip systems do s = 1,ns ! loop over slip systems
do t = 1,4 ! loop over dislocation types do t = 1,4 ! loop over dislocation types
if ( sign(1.0_pReal,math_mul3x3(m(:,s,t),surfaceNormal)) == sign(1.0_pReal,gdot(s,t)) ) then if ( sign(1.0_pReal,math_mul3x3(m(:,s,t),surfaceNormal)) == sign(1.0_pReal,gdot(s,t)) ) then
lineLength(s,t) = gdot(s,t) / constitutive_nonlocal_burgersPerSlipSystem(s,myInstance) & lineLength(s,t) = gdot(s,t) / constitutive_nonlocal_burgersPerSlipSystem(s,myInstance) &
* math_mul3x3(m(:,s,t),surfaceNormal) * area ! dislocation line length that leaves this interface per second * math_mul3x3(m(:,s,t),surfaceNormal) * area ! dislocation line length that leaves this interface per second
rhoDot(s,t) = rhoDot(s,t) - lineLength(s,t) / mesh_ipVolume(ip,el) ! subtract dislocation density rate (= line length over volume) that leaves through an interface from my dotState ... thisRhoDot(s,t) = thisRhoDot(s,t) - lineLength(s,t) / mesh_ipVolume(ip,el) ! subtract dislocation density rate (= line length over volume) that leaves through an interface from my dotState ...
if ( neighboring_el > 0 .and. neighboring_ip > 0 ) then if ( neighboring_el > 0 .and. neighboring_ip > 0 ) then
!***************************************************************************************************** !*****************************************************************************************************
@ -1226,7 +1258,15 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el))
enddo enddo
enddo enddo
enddo enddo
if (debugger) write(6,'(a,/,6(12(e12.5,x),/))') 'dislocation flux', lineLength/mesh_ipVolume(ip,el)*subdt, 0.0_pReal*rhoDotTransfer
totalRhoDot = totalRhoDot + thisRhoDot
totalRhoDipDot = totalRhoDipDot + thisRhoDipDot
if (debugger) then
!$OMP CRITICAL (write2out)
write(6,'(a,/,6(12(e12.5,x),/))') 'dislocation flux', thisRhoDot * subdt, thisRhoDipDot * subdt
!$OMP CRITICAL (write2out)
endif
!**************************************************************************** !****************************************************************************
@ -1236,11 +1276,13 @@ if (debugger) write(6,'(a,/,6(12(e12.5,x),/))') 'dislocation flux', lineLength/m
dDipMin(:,1) = constitutive_nonlocal_dDipMinEdgePerSlipSystem(:,myInstance) dDipMin(:,1) = constitutive_nonlocal_dDipMinEdgePerSlipSystem(:,myInstance)
dDipMin(:,2) = constitutive_nonlocal_dDipMinScrewPerSlipSystem(:,myInstance) dDipMin(:,2) = constitutive_nonlocal_dDipMinScrewPerSlipSystem(:,myInstance)
dDipMax(:,2) = constitutive_nonlocal_Gmod(myInstance) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) & dDipMax(:,2) = min( constitutive_nonlocal_Gmod(myInstance) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) &
/ ( 8.0_pReal * pi * abs(tauSlip) ) / ( 8.0_pReal * pi * abs(tauSlip) ), &
1.0_pReal / sqrt( 0.5_pReal * (rhoDip(:,1)+rhoDip(:,2)) ) )
dDipMax(:,1) = dDipMax(:,2) / ( 1.0_pReal - constitutive_nonlocal_nu(myInstance) ) dDipMax(:,1) = dDipMax(:,2) / ( 1.0_pReal - constitutive_nonlocal_nu(myInstance) )
dDipMax0(:,2) = constitutive_nonlocal_Gmod(myInstance) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) & dDipMax0(:,2) = min( constitutive_nonlocal_Gmod(myInstance) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) &
/ ( 8.0_pReal * pi * abs(subTauSlip0) ) / ( 8.0_pReal * pi * abs(subTauSlip0) ), &
1.0_pReal / sqrt( 0.5_pReal * (rhoDip(:,1)+rhoDip(:,2)) ) )
dDipMax0(:,1) = dDipMax0(:,2) / ( 1.0_pReal - constitutive_nonlocal_nu(myInstance) ) dDipMax0(:,1) = dDipMax0(:,2) / ( 1.0_pReal - constitutive_nonlocal_nu(myInstance) )
dDipMaxDot(:,1) = (dDipMax(:,1) - dDipMax0(:,1)) / subdt dDipMaxDot(:,1) = (dDipMax(:,1) - dDipMax0(:,1)) / subdt
@ -1251,76 +1293,91 @@ dDipMaxDot(:,2) = (dDipMax(:,2) - dDipMax0(:,2)) / subdt
!*** formation by glide !*** formation by glide
forall (c=1:2) & thisRhoDot = 0.0_pReal
rhoDotTransfer(:,c) = 2.0_pReal * dDipMax(:,c) / constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) & thisRhoDipDot = 0.0_pReal
* ( rho(:,2*c-1)*abs(gdot(:,2*c)) + rho(:,2*c)*abs(gdot(:,2*c-1)) )
if (debugger) write(6,'(a,/,6(12(e12.5,x),/))') 'dipole formation by glide', &
-rhoDotTransfer*subdt,-rhoDotTransfer*subdt,2.0_pReal*rhoDotTransfer*subdt
rhoDot(:,(/1,3/)) = rhoDot(:,(/1,3/)) - rhoDotTransfer ! subtract from positive single dislocation density of this character forall (c=1:2) &
rhoDot(:,(/2,4/)) = rhoDot(:,(/2,4/)) - rhoDotTransfer ! subtract from negative single dislocation density of this character thisRhoDipDot(:,c) = 4.0_pReal * dDipMax(:,c) / constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) &
rhoDipDot = rhoDipDot + 2.0_pReal * rhoDotTransfer ! add twice to dipole dislocation density of this character * ( rho(:,2*c-1)*abs(gdot(:,2*c)) + rho(:,2*c)*abs(gdot(:,2*c-1)) )
forall (t=1:4) &
thisRhoDot(:,t) = 0.5_pReal * thisRhoDipDot(:,mod(t+1,2)+1)
totalRhoDot = totalRhoDot + thisRhoDot
totalRhoDipDot = totalRhoDipDot + thisRhoDipDot
if (debugger) then
!$OMP CRITICAL (write2out)
write(6,'(a,/,6(12(e12.5,x),/))') 'dipole formation by glide', thisRhoDot * subdt, thisRhoDipDot * subdt
!$OMP CRITICAL (write2out)
endif
!*** athermal annihilation !*** athermal annihilation
forall (c=1:2) & thisRhoDot = 0.0_pReal
rhoDotTransfer(:,c) = - 2.0_pReal * dDipMin(:,c) / constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) & thisRhoDipDot = 0.0_pReal
* ( rho(:,2*c-1)*abs(gdot(:,2*c)) + rho(:,2*c)*abs(gdot(:,2*c-1)) )
if (debugger) write(6,'(a,/,6(12(e12.5,x),/))') 'athermal dipole annihilation', &
0.0_pReal*rhoDotTransfer,0.0_pReal*rhoDotTransfer,2.0_pReal*rhoDotTransfer*subdt
rhoDipDot = rhoDipDot + 2.0_pReal * rhoDotTransfer ! add twice to dipole dislocation density of this character forall (c=1:2) &
thisRhoDipDot(:,c) = - 4.0_pReal * dDipMin(:,c) / constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) &
* ( rho(:,2*c-1)*abs(gdot(:,2*c)) + rho(:,2*c)*abs(gdot(:,2*c-1)) & ! single hits single
+ 0.5_pReal * rhoDip(:,c) * (abs(gdot(:,2*c-1))+abs(gdot(:,2*c))) ) ! single knocks dipole constituent
totalRhoDot = totalRhoDot + thisRhoDot
totalRhoDipDot = totalRhoDipDot + thisRhoDipDot
if (debugger) then
!$OMP CRITICAL (write2out)
write(6,'(a,/,6(12(e12.5,x),/))') 'athermal dipole annihilation', thisRhoDot * subdt, thisRhoDipDot * subdt
!$OMP CRITICAL (write2out)
endif
!*** thermally activated annihilation !*** thermally activated annihilation
thisRhoDot = 0.0_pReal
thisRhoDipDot = 0.0_pReal
D = constitutive_nonlocal_D0(myInstance) * exp(-constitutive_nonlocal_Qsd(myInstance) / (kB * Temperature)) D = constitutive_nonlocal_D0(myInstance) * exp(-constitutive_nonlocal_Qsd(myInstance) / (kB * Temperature))
vClimb = constitutive_nonlocal_atomicVolume(myInstance) * D / ( 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)) ) & * constitutive_nonlocal_Gmod(myInstance) / ( 2.0_pReal * pi * (1.0_pReal-constitutive_nonlocal_nu(myInstance)) ) &
* 2.0_pReal / ( dDipMax(:,1) + dDipMin(:,1) ) * 2.0_pReal / ( dDipMax(:,1) + dDipMin(:,1) )
rhoDipDot(:,1) = rhoDipDot(:,1) - 4.0_pReal * rho(:,1) * vClimb / ( dDipMax(:,1) - dDipMin(:,1) ) thisRhoDipDot(:,1) = - 4.0_pReal * rho(:,1) * vClimb / ( dDipMax(:,1) - dDipMin(:,1) )
if (debugger) write(6,'(a,/,6(12(e12.5,x),/))') 'thermally activated dipole annihilation', &
0.0_pReal*rhoDotTransfer, 0.0_pReal*rhoDotTransfer, - 4.0_pReal * rho(:,1) * vClimb / ( dDipMax(:,1) - dDipMin(:,1) )*subdt, & totalRhoDot = totalRhoDot + thisRhoDot
0.0_pReal*vClimb totalRhoDipDot = totalRhoDipDot + thisRhoDipDot
if (debugger) then
!$OMP CRITICAL (write2out)
write(6,'(a,/,6(12(e12.5,x),/))') 'thermally activated dipole annihilation', thisRhoDot * subdt, thisRhoDipDot * subdt
!$OMP CRITICAL (write2out)
endif
! !*** formation by stress decrease = increase in dDipMax !*** formation by stress decrease = increase in dDipMax
! forall (s=1:ns, dDipMaxDot(s,1) > 0.0_pReal) & ! !!! MISSING !!!
! rhoDotTransfer(s,:) = 4.0_pReal * rho(s,(/1,3/)) * rho(s,(/2,4/)) * dDipMax0(s,:) * dDipMaxDot(s,:)
! if (debugger) write(6,'(a,/,6(12(e12.5,x),/))') 'dipole formation by stress decrease',&
! -rhoDotTransfer*subdt,-rhoDotTransfer*subdt,2.0_pReal*rhoDotTransfer*subdt
! rhoDot(:,(/1,3/)) = rhoDot(:,(/1,3/)) - rhoDotTransfer ! subtract from positive single dislocation density of this character
! rhoDot(:,(/2,4/)) = rhoDot(:,(/2,4/)) - rhoDotTransfer ! subtract from negative single dislocation density of this character
! rhoDipDot = rhoDipDot + 2.0_pReal * rhoDotTransfer ! add twice to dipole dislocation density of this character
! !*** dipole dissociation by increased stress = decrease in dDipMax !*** dipole dissociation by increased stress = decrease in dDipMax
! forall (s=1:ns, dDipMaxDot(s,1) < 0.0_pReal) & ! !!! MISSING !!!
! rhoDotTransfer(s,:) = 0.5_pReal * rhoDip(s,:) * dDipMaxDot(s,:) / (dDipMax0(s,:) - dDipMin(s,:))
! if (debugger) write(6,'(a,/,6(12(e12.5,x),/))') 'dipole formation by stress decrease',&
! -rhoDotTransfer*subdt,-rhoDotTransfer*subdt,2.0_pReal*rhoDotTransfer*subdt
! rhoDot(:,(/1,3/)) = rhoDot(:,(/1,3/)) - rhoDotTransfer ! subtract from positive single dislocation density of this character
! rhoDot(:,(/2,4/)) = rhoDot(:,(/2,4/)) - rhoDotTransfer ! subtract from negative single dislocation density of this character
! rhoDipDot = rhoDipDot + 2.0_pReal * rhoDotTransfer ! add twice to dipole dislocation density of this character
!**************************************************************************** !****************************************************************************
!*** assign the rates of dislocation densities to my dotState !*** assign the rates of dislocation densities to my dotState
dotState(1,ip,el)%p(1:4*ns) = reshape(rhoDot,(/4*ns/)) dotState(1,ip,el)%p(1:4*ns) = reshape(totalRhoDot,(/4*ns/))
dotState(1,ip,el)%p(4*ns+1:6*ns) = reshape(rhoDipDot,(/2*ns/)) dotState(1,ip,el)%p(4*ns+1:6*ns) = reshape(totalRhoDipDot,(/2*ns/))
if (debugger) write(6,'(a,/,4(12(e12.5,x),/))') 'deltaRho:',rhoDot*subdt if (debugger) then
if (debugger) write(6,'(a,/,2(12(e12.5,x),/))') 'deltaRhoDip:',rhoDipDot*subdt !$OMP CRITICAL (write2out)
write(6,'(a,/,4(12(e12.5,x),/))') 'deltaRho:', totalRhoDot * subdt
write(6,'(a,/,2(12(e12.5,x),/))') 'deltaRhoDip:', totalRhoDipDot * subdt
!$OMP CRITICAL (write2out)
endif
endsubroutine endsubroutine
@ -1462,10 +1519,9 @@ do o = 1,phase_Noutput(material_phase(g,ip,el))
sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance) sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance)
tau = math_mul6x6( Tstar_v + state(g,ip,el)%p(8*ns+1:8*ns+6), lattice_Sslip_v(:,sLattice,myStructure) ) tau = math_mul6x6( Tstar_v + state(g,ip,el)%p(8*ns+1:8*ns+6), lattice_Sslip_v(:,sLattice,myStructure) )
if (state(g,ip,el)%p(4*ns+s) > 0.0_pReal) then if (state(g,ip,el)%p(7*ns+s) > 0.0_pReal) then
v = constitutive_nonlocal_v0PerSlipSystem(s,myInstance) & v = constitutive_nonlocal_v0PerSlipSystem(s,myInstance) &
* exp( - ( state(g,ip,el)%p(7*ns+s) - abs(tau) ) * constitutive_nonlocal_burgersPerSlipSystem(s,myInstance)**2.0_pReal & * exp( - constitutive_nonlocal_G / ( kB * Temperature) * (1.0_pReal - (abs(tau)/state(g,ip,el)%p(7*ns+s)) ) ) &
/ ( kB * Temperature * sqrt(state(g,ip,el)%p(4*ns+s)) ) ) &
* sign(1.0_pReal,tau) * sign(1.0_pReal,tau)
else else
v = 0.0_pReal v = 0.0_pReal

63
code/crystallite.f90
View File

@ -397,6 +397,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
crystallite_subF(:,:,g,i,e) = crystallite_subF0(:,:,g,i,e) + & crystallite_subF(:,:,g,i,e) = crystallite_subF0(:,:,g,i,e) + &
crystallite_subStep(g,i,e) * & crystallite_subStep(g,i,e) * &
(crystallite_partionedF(:,:,g,i,e) - crystallite_partionedF0(:,:,g,i,e)) (crystallite_partionedF(:,:,g,i,e) - crystallite_partionedF0(:,:,g,i,e))
crystallite_Fe(:,:,g,i,e) = math_mul33x33(crystallite_subF(:,:,g,i,e),crystallite_invFp(:,:,g,i,e))
crystallite_subdt(g,i,e) = crystallite_subStep(g,i,e) * crystallite_dt(g,i,e) crystallite_subdt(g,i,e) = crystallite_subStep(g,i,e) * crystallite_dt(g,i,e)
crystallite_converged(g,i,e) = .false. ! start out non-converged crystallite_converged(g,i,e) = .false. ! start out non-converged
endif endif
@ -431,11 +432,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains do g = 1,myNgrains
debugger = (e == 1 .and. i == 1 .and. g == 1) ! debugger = (e == 1 .and. i == 1 .and. g == 1)
if (crystallite_todo(g,i,e)) & ! all undone crystallites if (crystallite_todo(g,i,e)) & ! all undone crystallites
call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), &
crystallite_Fp(:,:,g,i,e), crystallite_invFp(:,:,g,i,e), & crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), &
crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), g, i, e) crystallite_subdt(g,i,e), g, i, e)
enddo; enddo; enddo enddo; enddo; enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
!$OMP PARALLEL DO !$OMP PARALLEL DO
@ -443,7 +444,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains do g = 1,myNgrains
debugger = (e == 1 .and. i == 1 .and. g == 1) ! debugger = (e == 1 .and. i == 1 .and. g == 1)
if (crystallite_todo(g,i,e)) then ! all undone crystallites if (crystallite_todo(g,i,e)) then ! all undone crystallites
crystallite_stateConverged(g,i,e) = crystallite_updateState(g,i,e) ! update state crystallite_stateConverged(g,i,e) = crystallite_updateState(g,i,e) ! update state
crystallite_temperatureConverged(g,i,e) = crystallite_updateTemperature(g,i,e) ! update temperature crystallite_temperatureConverged(g,i,e) = crystallite_updateTemperature(g,i,e) ! update temperature
@ -451,7 +452,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
endif endif
enddo; enddo; enddo enddo; enddo; enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
write(6,*) count(.not. crystallite_onTrack(1,:,:)),'IPs not onTrack after preguess for state' write(6,*) count(crystallite_onTrack(1,:,:)),'IPs onTrack after preguess for state'
! --+>> state loop <<+-- ! --+>> state loop <<+--
@ -474,7 +475,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains do g = 1,myNgrains
debugger = (e == 1 .and. i == 1 .and. g == 1) ! debugger = (e == 1 .and. i == 1 .and. g == 1)
if (crystallite_todo(g,i,e)) & ! all undone crystallites if (crystallite_todo(g,i,e)) & ! all undone crystallites
crystallite_onTrack(g,i,e) = crystallite_integrateStress(g,i,e) crystallite_onTrack(g,i,e) = crystallite_integrateStress(g,i,e)
enddo enddo
@ -482,7 +483,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
enddo enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
write(6,*) count(.not. crystallite_onTrack(1,:,:)),'IPs not onTrack after stress integration' write(6,*) count(crystallite_onTrack(1,:,:)),'IPs onTrack after stress integration'
crystallite_todo = crystallite_todo .and. crystallite_onTrack crystallite_todo = crystallite_todo .and. crystallite_onTrack
if (any(.not. crystallite_onTrack .and. .not. crystallite_localConstitution)) & if (any(.not. crystallite_onTrack .and. .not. crystallite_localConstitution)) &
@ -510,11 +511,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains do g = 1,myNgrains
debugger = (e == 1 .and. i == 1 .and. g == 1) ! debugger = (e == 1 .and. i == 1 .and. g == 1)
if (crystallite_todo(g,i,e)) & ! all undone crystallites if (crystallite_todo(g,i,e)) & ! all undone crystallites
call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), & call constitutive_collectDotState(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), &
crystallite_Fp(:,:,g,i,e), crystallite_invFp(:,:,g,i,e), & crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), &
crystallite_Temperature(g,i,e), crystallite_subdt(g,i,e), g, i, e) crystallite_subdt(g,i,e), g, i, e)
enddo; enddo; enddo enddo; enddo; enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
!$OMP PARALLEL DO !$OMP PARALLEL DO
@ -522,7 +523,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains do g = 1,myNgrains
debugger = (e == 1 .and. i == 1 .and. g == 1) ! debugger = (e == 1 .and. i == 1 .and. g == 1)
if (crystallite_todo(g,i,e)) then ! all undone crystallites if (crystallite_todo(g,i,e)) then ! all undone crystallites
crystallite_stateConverged(g,i,e) = crystallite_updateState(g,i,e) ! update state crystallite_stateConverged(g,i,e) = crystallite_updateState(g,i,e) ! update state
crystallite_temperatureConverged(g,i,e) = crystallite_updateTemperature(g,i,e) ! update temperature crystallite_temperatureConverged(g,i,e) = crystallite_updateTemperature(g,i,e) ! update temperature
@ -543,9 +544,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
if (any(.not. crystallite_onTrack .and. .not. crystallite_localConstitution)) & if (any(.not. crystallite_onTrack .and. .not. crystallite_localConstitution)) &
crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! all nonlocal crystallites can be skipped crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! all nonlocal crystallites can be skipped
write(6,*) count(.not. crystallite_onTrack(1,:,:)),'IPs not onTrack after state update' write(6,*) count(crystallite_onTrack(1,:,:)),'IPs onTrack after state update'
write(6,*) count(crystallite_converged(1,:,:)),'IPs converged' write(6,*) count(crystallite_converged(1,:,:)),'IPs converged'
write(6,*) count(crystallite_todo(1,:,:)),'IPs todo' write(6,*) count(crystallite_todo(1,:,:)),'IPs todo'
write(6,*)
enddo ! crystallite convergence loop enddo ! crystallite convergence loop
@ -748,7 +750,7 @@ endsubroutine
! update the microstructure ! update the microstructure
constitutive_state(g,i,e)%p(1:mySize) = constitutive_state(g,i,e)%p(1:mySize) - residuum constitutive_state(g,i,e)%p(1:mySize) = constitutive_state(g,i,e)%p(1:mySize) - residuum
call constitutive_microstructure(crystallite_subTemperature0(g,i,e), crystallite_subFp0, g, i, e) call constitutive_microstructure(crystallite_Temperature(g,i,e), crystallite_Fe, crystallite_Fp, g, i, e)
! setting flag to true if state is below relative tolerance, otherwise set it to false ! setting flag to true if state is below relative tolerance, otherwise set it to false
@ -906,7 +908,7 @@ endsubroutine
AB, & AB, &
BTA BTA
real(pReal), dimension(6):: Tstar_v ! 2nd Piola-Kirchhoff Stress in Mandel-Notation real(pReal), dimension(6):: Tstar_v ! 2nd Piola-Kirchhoff Stress in Mandel-Notation
real(pReal), dimension(9,9):: dLp_constitutive, & ! partial derivative of plastic velocity gradient calculated by constitutive law real(pReal), dimension(9,9):: dLpdT_constitutive, & ! partial derivative of plastic velocity gradient calculated by constitutive law
dTdLp, & ! partial derivative of 2nd Piola-Kirchhoff stress dTdLp, & ! partial derivative of 2nd Piola-Kirchhoff stress
dRdLp, & ! partial derivative of residuum (Jacobian for NEwton-Raphson scheme) dRdLp, & ! partial derivative of residuum (Jacobian for NEwton-Raphson scheme)
invdRdLp ! inverse of dRdLp invdRdLp ! inverse of dRdLp
@ -974,7 +976,15 @@ LpLoop: do
NiterationStress = NiterationStress + 1 NiterationStress = NiterationStress + 1
! too many loops required ? ! too many loops required ?
if (NiterationStress > nStress) return if (NiterationStress > nStress) then
if (debugger) then
!$OMP CRITICAL (write2out)
write(6,*) '::: integrateStress reached loop limit',g,i,e
write(6,*)
!$OMPEND CRITICAL (write2out)
endif
return
endif
B = math_I3 - crystallite_subdt(g,i,e)*Lpguess B = math_I3 - crystallite_subdt(g,i,e)*Lpguess
BT = transpose(B) BT = transpose(B)
@ -988,11 +998,19 @@ LpLoop: do
! calculate plastic velocity gradient and its tangent according to constitutive law ! calculate plastic velocity gradient and its tangent according to constitutive law
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks) call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
call constitutive_LpAndItsTangent(Lp_constitutive, dLp_constitutive, Tstar_v, crystallite_Temperature(g,i,e), g, i, e) call constitutive_LpAndItsTangent(Lp_constitutive, dLpdT_constitutive, Tstar_v, crystallite_Temperature(g,i,e), g, i, e)
call system_clock(count=tock,count_rate=tickrate,count_max=maxticks) call system_clock(count=tock,count_rate=tickrate,count_max=maxticks)
debug_cumLpCalls = debug_cumLpCalls + 1_pInt debug_cumLpCalls = debug_cumLpCalls + 1_pInt
debug_cumLpTicks = debug_cumLpTicks + tock-tick debug_cumLpTicks = debug_cumLpTicks + tock-tick
if (tock < tick) debug_cumLpTicks = debug_cumLpTicks + maxticks if (tock < tick) debug_cumLpTicks = debug_cumLpTicks + maxticks
if (debugger) then
!$OMP CRITICAL (write2out)
write(6,*) '::: integrateStress at iteration', NiterationStress
write(6,*)
write(6,'(a19,3(i3,x),/,3(3(f20.7,x)/))') 'Lp_constitutive at ',g,i,e,Lp_constitutive
write(6,'(a11,3(i3,x),/,3(3(f20.7,x)/))') 'Lpguess at ',g,i,e,Lpguess
!$OMPEND CRITICAL (write2out)
endif
! update current residuum ! update current residuum
residuum = Lpguess - Lp_constitutive residuum = Lpguess - Lp_constitutive
@ -1037,10 +1055,9 @@ LpLoop: do
if (mod(jacoCounter, iJacoLpresiduum) == 0_pInt) then if (mod(jacoCounter, iJacoLpresiduum) == 0_pInt) then
dTdLp = 0.0_pReal dTdLp = 0.0_pReal
forall (h=1:3,j=1:3,k=1:3,l=1:3,m=1:3) & forall (h=1:3,j=1:3,k=1:3,l=1:3,m=1:3) &
dTdLp(3*(h-1)+j,3*(k-1)+l) = dTdLp(3*(h-1)+j,3*(k-1)+l) + & dTdLp(3*(h-1)+j,3*(k-1)+l) = dTdLp(3*(h-1)+j,3*(k-1)+l) + C(h,j,l,m)*AB(k,m)+C(h,j,m,l)*BTA(m,k)
C(h,j,l,m)*AB(k,m)+C(h,j,m,l)*BTA(m,k)
dTdLp = -0.5_pReal*crystallite_subdt(g,i,e)*dTdLp dTdLp = -0.5_pReal*crystallite_subdt(g,i,e)*dTdLp
dRdLp = math_identity2nd(9) - math_mul99x99(dLp_constitutive,dTdLp) dRdLp = math_identity2nd(9) - math_mul99x99(dLpdT_constitutive,dTdLp)
invdRdLp = 0.0_pReal invdRdLp = 0.0_pReal
call math_invert(9,dRdLp,invdRdLp,dummy,error) ! invert dR/dLp --> dLp/dR call math_invert(9,dRdLp,invdRdLp,dummy,error) ! invert dR/dLp --> dLp/dR
if (error) then if (error) then
@ -1048,10 +1065,10 @@ LpLoop: do
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) '::: integrateStress failed on dR/dLp inversion at iteration', NiterationStress write(6,*) '::: integrateStress failed on dR/dLp inversion at iteration', NiterationStress
write(6,*) write(6,*)
write(6,'(a9,3(i3,x),/,9(9(f12.7,x)/))') 'dRdLp at ',g,i,e,dRdLp write(6,'(a9,3(i3,x),/,9(9(f15.3,x)/))') 'dRdLp at ',g,i,e,dRdLp
write(6,'(a20,3(i3,x),/,9(9(e12.2,x)/))') 'dLp_constitutive at ',g,i,e,dLp_constitutive write(6,'(a20,3(i3,x),/,9(9(f15.3,x)/))') 'dLpdT_constitutive at ',g,i,e,dLpdT_constitutive
write(6,'(a19,3(i3,x),/,3(3(f12.7,x)/))') 'Lp_constitutive at ',g,i,e,Lp_constitutive write(6,'(a19,3(i3,x),/,3(3(f20.7,x)/))') 'Lp_constitutive at ',g,i,e,Lp_constitutive
write(6,'(a11,3(i3,x),/,3(3(f12.7,x)/))') 'Lpguess at ',g,i,e,Lpguess write(6,'(a11,3(i3,x),/,3(3(f20.7,x)/))') 'Lpguess at ',g,i,e,Lpguess
!$OMPEND CRITICAL (write2out) !$OMPEND CRITICAL (write2out)
endif endif
return return