constitutive_nonlocal:

- read in activation energy for dislocation glide from material.config
- changed naming of dDipMin/Max to dLower/dUpper
- added new outputs: rho_dot, rho_dot_dip, rho_dot_gen, rho_dot_sgl2dip, rho_dot_dip2sgl, rho_dot_ann_ath, rho_dot_ann_the, rho_dot_flux, d_upper_edge, d_upper_screw, d_upper_dot_edge, d_upper_dot_screw
- poisson's ratio is now calculated from elastic constants
- microstrucutre has state as first argument, since this is our output variable
- periodic boundary conditions are taken into account for fluxes and internal stresses. for the moment, flag has to be set in constitutive_nonlocal. 
- corrected calculation for dipole formation by glide
- added terms for dipole formation/annihilation by stress decrease/increase 

constitutive:
- passing of arguments is adapted for constitutive_nonlocal model

crystallite:
- in stiffness calculation: call to collect_dotState used wrong arguments
- crystallite_postResults uses own Tstar_v and temperature, no need for passing them from materialpoint_postResults

homogenization:
- crystallite_postResults uses own Tstar_v and temperature, no need for passing them from materialpoint_postResults

IO:
- changed error message 229

material.config:
- changed example for nonlocal constitution according to constitutive_nonlocal

all:
- added some flush statements
This commit is contained in:
Christoph Kords 2009-10-20 14:36:03 +00:00
parent 48e642bb15
commit c1ee34d235
7 changed files with 490 additions and 197 deletions

View File

@ -994,7 +994,7 @@ endfunction
case (228)
msg = 'Non-positive minimum stable dipole distance'
case (229)
msg = 'Hardening interaction coefficients below one'
msg = 'Non-positive hardening interaction coefficients'
case (230)
msg = 'Non-positive atomic volume'
case (231)
@ -1129,6 +1129,7 @@ endfunction
endif
endif
write(6,'(a38)') '+------------------------------------+'
call flush(6)
endsubroutine

View File

@ -205,6 +205,7 @@ subroutine constitutive_init()
write(6,'(a32,x,7(i5,x))') 'constitutive_sizePostResults: ', shape(constitutive_sizePostResults)
write(6,*)
write(6,'(a32,x,7(i5,x))') 'maxSizeState: ', constitutive_maxSizeState
write(6,'(a32,x,7(i5,x))') 'maxSizeDotState: ', constitutive_maxSizeDotState
write(6,'(a32,x,7(i5,x))') 'maxSizePostResults: ', constitutive_maxSizePostResults
return
@ -293,7 +294,7 @@ real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)
call constitutive_dislotwin_microstructure(Temperature,constitutive_state,ipc,ip,el)
case (constitutive_nonlocal_label)
call constitutive_nonlocal_microstructure(Temperature, Fe, Fp, constitutive_state, ipc, ip, el)
call constitutive_nonlocal_microstructure(constitutive_state, Temperature, Fe, Fp, ipc, ip, el)
end select
@ -374,7 +375,7 @@ subroutine constitutive_collectDotState(Tstar_v, subTstar0_v, Fe, Fp, Temperatur
!* Definition of variables
integer(pInt), intent(in) :: ipc,ip,el
real(pReal), intent(in) :: Temperature, subdt
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: Fe, Fp
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: Fe, Fp
real(pReal), dimension(6), intent(in) :: Tstar_v, subTstar0_v
select case (phase_constitution(material_phase(ipc,ip,el)))
@ -443,7 +444,7 @@ function constitutive_dotTemperature(Tstar_v,Temperature,ipc,ip,el)
endfunction
pure function constitutive_postResults(Tstar_v,Temperature,dt,ipc,ip,el)
pure function constitutive_postResults(Tstar_v,subTstar0_v,Temperature,dt,subdt,ipc,ip,el)
!*********************************************************************
!* return array of constitutive results *
!* INPUT: *
@ -463,8 +464,8 @@ pure function constitutive_postResults(Tstar_v,Temperature,dt,ipc,ip,el)
!* Definition of variables
integer(pInt), intent(in) :: ipc,ip,el
real(pReal), intent(in) :: dt,Temperature
real(pReal), dimension(6), intent(in) :: Tstar_v
real(pReal), intent(in) :: dt, Temperature, subdt
real(pReal), dimension(6), intent(in) :: Tstar_v, subTstar0_v
real(pReal), dimension(constitutive_sizePostResults(ipc,ip,el)) :: constitutive_postResults
constitutive_postResults = 0.0_pReal
@ -480,10 +481,10 @@ pure function constitutive_postResults(Tstar_v,Temperature,dt,ipc,ip,el)
constitutive_postResults = constitutive_dislotwin_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el)
case (constitutive_nonlocal_label)
constitutive_postResults = constitutive_nonlocal_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el)
constitutive_postResults = constitutive_nonlocal_postResults(Tstar_v, subTstar0_v, Temperature, dt, subdt, constitutive_state,&
constitutive_subState0, constitutive_dotstate, ipc, ip, el)
end select
return
endfunction

View File

@ -27,7 +27,6 @@ character(len=16), dimension(3), parameter :: constitutive_nonlocal_listDependen
'tauSlipThreshold', &
'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 :: constitutive_nonlocal_G = 2.5e-19_pReal
!* Definition of global variables
integer(pInt), dimension(:), allocatable :: constitutive_nonlocal_sizeDotState, & ! number of dotStates
@ -51,9 +50,10 @@ real(pReal), dimension(:), allocatable :: constitutive_nonlocal_
constitutive_nonlocal_C44, & ! C44 element in elasticity matrix
constitutive_nonlocal_Gmod, & ! shear modulus
constitutive_nonlocal_nu, & ! poisson's ratio
constitutive_nonlocal_Q0, & ! activation energy for dislocation glide
constitutive_nonlocal_atomicVolume, & ! atomic volume
constitutive_nonlocal_D0, & ! prefactor for self-diffusion coefficient
constitutive_nonlocal_Qsd, & ! activation energy for dislocation climb
constitutive_nonlocal_Qsd, & ! activation enthalpy for diffusion
constitutive_nonlocal_relevantRho ! dislocation density considered relevant
real(pReal), dimension(:,:,:), allocatable :: constitutive_nonlocal_Cslip_66 ! elasticity matrix in Mandel notation for each instance
real(pReal), dimension(:,:,:,:,:), allocatable :: constitutive_nonlocal_Cslip_3333 ! elasticity matrix for each instance
@ -69,16 +69,18 @@ real(pReal), dimension(:,:), allocatable :: constitutive_nonlocal_
constitutive_nonlocal_lambda0PerSlipSystem, & ! mean free path prefactor for each slip system and instance
constitutive_nonlocal_burgersPerSlipFamily, & ! absolute length of burgers vector [m] for each family and instance
constitutive_nonlocal_burgersPerSlipSystem, & ! absolute length of burgers vector [m] for each slip system and instance
constitutive_nonlocal_dDipMinEdgePerSlipFamily, & ! minimum stable edge dipole height for each family and instance
constitutive_nonlocal_dDipMinEdgePerSlipSystem, & ! minimum stable edge dipole height for each slip system and instance
constitutive_nonlocal_dDipMinScrewPerSlipFamily, & ! minimum stable screw dipole height for each family and instance
constitutive_nonlocal_dDipMinScrewPerSlipSystem, & ! minimum stable screw dipole height for each slip system and instance
constitutive_nonlocal_dLowerEdgePerSlipFamily, & ! minimum stable edge dipole height for each family and instance
constitutive_nonlocal_dLowerEdgePerSlipSystem, & ! minimum stable edge dipole height for each slip system and instance
constitutive_nonlocal_dLowerScrewPerSlipFamily, & ! minimum stable screw dipole height for each family and instance
constitutive_nonlocal_dLowerScrewPerSlipSystem, & ! minimum stable screw dipole height for each slip system and instance
constitutive_nonlocal_interactionSlipSlip ! coefficients for slip-slip interaction for each interaction type and instance
real(pReal), dimension(:,:,:), allocatable :: constitutive_nonlocal_forestProjectionEdge, & ! matrix of forest projections of edge dislocations for each instance
constitutive_nonlocal_forestProjectionScrew, & ! matrix of forest projections of screw dislocations for each instance
constitutive_nonlocal_interactionMatrixSlipSlip ! interaction matrix of the different slip systems for each instance
logical periodicBC = .false.
CONTAINS
!****************************************
@ -158,7 +160,7 @@ character(len=1024) line
write(6,*)
write(6,'(a20,a20,a12)') '<<<+- constitutive_',constitutive_nonlocal_label,' init -+>>>'
write(6,*) '$Id$'
write(6,*) '$Id$'
write(6,*)
maxNinstance = count(phase_constitution == constitutive_nonlocal_label)
@ -199,6 +201,7 @@ allocate(constitutive_nonlocal_C33(maxNinstance))
allocate(constitutive_nonlocal_C44(maxNinstance))
allocate(constitutive_nonlocal_Gmod(maxNinstance))
allocate(constitutive_nonlocal_nu(maxNinstance))
allocate(constitutive_nonlocal_Q0(maxNinstance))
allocate(constitutive_nonlocal_atomicVolume(maxNinstance))
allocate(constitutive_nonlocal_D0(maxNinstance))
allocate(constitutive_nonlocal_Qsd(maxNinstance))
@ -212,6 +215,7 @@ constitutive_nonlocal_C13 = 0.0_pReal
constitutive_nonlocal_C33 = 0.0_pReal
constitutive_nonlocal_C44 = 0.0_pReal
constitutive_nonlocal_Gmod = 0.0_pReal
constitutive_nonlocal_Q0 = 0.0_pReal
constitutive_nonlocal_atomicVolume = 0.0_pReal
constitutive_nonlocal_D0 = 0.0_pReal
constitutive_nonlocal_Qsd = 0.0_pReal
@ -230,8 +234,8 @@ allocate(constitutive_nonlocal_v0PerSlipFamily(lattice_maxNslipFamily, maxNinsta
allocate(constitutive_nonlocal_burgersPerSlipFamily(lattice_maxNslipFamily, maxNinstance))
allocate(constitutive_nonlocal_Lambda0PerSlipFamily(lattice_maxNslipFamily, maxNinstance))
allocate(constitutive_nonlocal_interactionSlipSlip(lattice_maxNinteraction, maxNinstance))
allocate(constitutive_nonlocal_dDipMinEdgePerSlipFamily(lattice_maxNslipFamily, maxNinstance))
allocate(constitutive_nonlocal_dDipMinScrewPerSlipFamily(lattice_maxNslipFamily, maxNinstance))
allocate(constitutive_nonlocal_dLowerEdgePerSlipFamily(lattice_maxNslipFamily, maxNinstance))
allocate(constitutive_nonlocal_dLowerScrewPerSlipFamily(lattice_maxNslipFamily, maxNinstance))
constitutive_nonlocal_rhoEdgePos0 = 0.0_pReal
constitutive_nonlocal_rhoEdgeNeg0 = 0.0_pReal
constitutive_nonlocal_rhoScrewPos0 = 0.0_pReal
@ -242,8 +246,8 @@ constitutive_nonlocal_v0PerSlipFamily = 0.0_pReal
constitutive_nonlocal_burgersPerSlipFamily = 0.0_pReal
constitutive_nonlocal_lambda0PerSlipFamily = 0.0_pReal
constitutive_nonlocal_interactionSlipSlip = 0.0_pReal
constitutive_nonlocal_dDipMinEdgePerSlipFamily = 0.0_pReal
constitutive_nonlocal_dDipMinScrewPerSlipFamily = 0.0_pReal
constitutive_nonlocal_dLowerEdgePerSlipFamily = 0.0_pReal
constitutive_nonlocal_dLowerScrewPerSlipFamily = 0.0_pReal
!*** readout data from material.config file
@ -308,10 +312,12 @@ do
forall (f = 1:lattice_maxNslipFamily) constitutive_nonlocal_burgersPerSlipFamily(f,i) = IO_floatValue(line,positions,1+f)
case('ddipminedge')
forall (f = 1:lattice_maxNslipFamily) &
constitutive_nonlocal_dDipMinEdgePerSlipFamily(f,i) = IO_floatValue(line,positions,1+f)
constitutive_nonlocal_dLowerEdgePerSlipFamily(f,i) = IO_floatValue(line,positions,1+f)
case('ddipminscrew')
forall (f = 1:lattice_maxNslipFamily) &
constitutive_nonlocal_dDipMinScrewPerSlipFamily(f,i) = IO_floatValue(line,positions,1+f)
constitutive_nonlocal_dLowerScrewPerSlipFamily(f,i) = IO_floatValue(line,positions,1+f)
case('q0')
constitutive_nonlocal_Q0(i) = IO_floatValue(line,positions,2)
case('atomicvolume')
constitutive_nonlocal_atomicVolume(i) = IO_floatValue(line,positions,2)
case('d0')
@ -337,6 +343,9 @@ enddo
if (myStructure < 1 .or. myStructure > 3) call IO_error(205)
if (sum(constitutive_nonlocal_Nslip(:,i)) <= 0_pInt) call IO_error(225)
do o = 1,maxval(phase_Noutput)
if(len(constitutive_nonlocal_output(o,i)) > 64) call IO_error(666)
enddo
do f = 1,lattice_maxNslipFamily
if (constitutive_nonlocal_Nslip(f,i) > 0_pInt) then
if (constitutive_nonlocal_rhoEdgePos0(f,i) < 0.0_pReal) call IO_error(220)
@ -348,12 +357,13 @@ enddo
if (constitutive_nonlocal_burgersPerSlipFamily(f,i) <= 0.0_pReal) call IO_error(221)
if (constitutive_nonlocal_v0PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(226)
if (constitutive_nonlocal_lambda0PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(227)
if (constitutive_nonlocal_dDipMinEdgePerSlipFamily(f,i) <= 0.0_pReal) call IO_error(228)
if (constitutive_nonlocal_dDipMinScrewPerSlipFamily(f,i) <= 0.0_pReal) call IO_error(228)
if (constitutive_nonlocal_dLowerEdgePerSlipFamily(f,i) <= 0.0_pReal) call IO_error(228)
if (constitutive_nonlocal_dLowerScrewPerSlipFamily(f,i) <= 0.0_pReal) call IO_error(228)
endif
enddo
if (any(constitutive_nonlocal_interactionSlipSlip(1:maxval(lattice_interactionSlipSlip(:,:,myStructure)),i) < 1.0_pReal)) &
if (any(constitutive_nonlocal_interactionSlipSlip(1:maxval(lattice_interactionSlipSlip(:,:,myStructure)),i) < 0.0_pReal)) &
call IO_error(229)
if (constitutive_nonlocal_Q0(i) <= 0.0_pReal) call IO_error(-1)
if (constitutive_nonlocal_atomicVolume(i) <= 0.0_pReal) call IO_error(230)
if (constitutive_nonlocal_D0(i) <= 0.0_pReal) call IO_error(231)
if (constitutive_nonlocal_Qsd(i) <= 0.0_pReal) call IO_error(232)
@ -362,7 +372,7 @@ enddo
!*** determine total number of active slip systems
constitutive_nonlocal_Nslip(:,i) = min( lattice_NslipSystem(:, myStructure), constitutive_nonlocal_Nslip(:,i) ) ! we can't use more slip systems per family than specified in lattice
constitutive_nonlocal_Nslip(:,i) = min( lattice_NslipSystem(:, myStructure), constitutive_nonlocal_Nslip(:,i) ) ! we can't use more slip systems per family than specified in lattice
constitutive_nonlocal_totalNslip(i) = sum(constitutive_nonlocal_Nslip(:,i))
enddo
@ -381,11 +391,11 @@ constitutive_nonlocal_v0PerSlipSystem = 0.0_pReal
allocate(constitutive_nonlocal_lambda0PerSlipSystem(maxTotalNslip, maxNinstance))
constitutive_nonlocal_lambda0PerSlipSystem = 0.0_pReal
allocate(constitutive_nonlocal_dDipMinEdgePerSlipSystem(maxTotalNslip, maxNinstance))
constitutive_nonlocal_dDipMinEdgePerSlipSystem = 0.0_pReal
allocate(constitutive_nonlocal_dLowerEdgePerSlipSystem(maxTotalNslip, maxNinstance))
constitutive_nonlocal_dLowerEdgePerSlipSystem = 0.0_pReal
allocate(constitutive_nonlocal_dDipMinScrewPerSlipSystem(maxTotalNslip, maxNinstance))
constitutive_nonlocal_dDipMinScrewPerSlipSystem = 0.0_pReal
allocate(constitutive_nonlocal_dLowerScrewPerSlipSystem(maxTotalNslip, maxNinstance))
constitutive_nonlocal_dLowerScrewPerSlipSystem = 0.0_pReal
allocate(constitutive_nonlocal_forestProjectionEdge(maxTotalNslip, maxTotalNslip, maxNinstance))
constitutive_nonlocal_forestProjectionEdge = 0.0_pReal
@ -422,9 +432,10 @@ do i = 1,maxNinstance
!*** determine size of postResults array
do o = 1,maxval(phase_Noutput)
do o = 1,maxval(phase_Noutput)
select case(constitutive_nonlocal_output(o,i))
case( 'rho', &
'delta', &
'rho_edge', &
'rho_screw', &
'excess_rho', &
@ -432,11 +443,24 @@ do i = 1,maxNinstance
'excess_rho_screw', &
'rho_forest', &
'rho_dip', &
'delta_dip', &
'rho_edge_dip', &
'rho_screw_dip', &
'shearrate', &
'resolvedstress', &
'resistance')
'resistance', &
'rho_dot', &
'rho_dot_dip', &
'rho_dot_gen', &
'rho_dot_sgl2dip', &
'rho_dot_dip2sgl', &
'rho_dot_ann_ath', &
'rho_dot_ann_the', &
'rho_dot_flux', &
'd_upper_edge', &
'd_upper_screw', &
'd_upper_dot_edge', &
'd_upper_dot_screw' )
mySize = constitutive_nonlocal_totalNslip(i)
case default
mySize = 0_pInt
@ -476,8 +500,8 @@ do i = 1,maxNinstance
constitutive_nonlocal_Cslip_3333(:,:,:,:,i) = math_Voigt66to3333(constitutive_nonlocal_Cslip_66(:,:,i))
constitutive_nonlocal_Gmod(i) = constitutive_nonlocal_C44(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)
constitutive_nonlocal_nu(i) = 0.5_pReal * constitutive_nonlocal_C12(i) &
/ (constitutive_nonlocal_C12(i) + constitutive_nonlocal_C44(i))
do s1 = 1,ns
@ -489,8 +513,8 @@ do i = 1,maxNinstance
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)
constitutive_nonlocal_dLowerEdgePerSlipSystem(s1,i) = constitutive_nonlocal_dLowerEdgePerSlipFamily(f,i)
constitutive_nonlocal_dLowerScrewPerSlipSystem(s1,i) = constitutive_nonlocal_dLowerScrewPerSlipFamily(f,i)
do s2 = 1,ns
@ -665,7 +689,7 @@ endfunction
!*********************************************************************
!* calculates quantities characterizing the microstructure *
!*********************************************************************
subroutine constitutive_nonlocal_microstructure(Temperature, Fe, Fp, state, g, ip, el)
subroutine constitutive_nonlocal_microstructure(state, Temperature, Fe, Fp, g, ip, el)
use prec, only: pReal, &
pInt, &
@ -724,7 +748,7 @@ integer(pInt) myInstance, & ! current instance
neighboring_ip, & ! integration point of my neighbor
n, & ! index of my current neighbor
s, & ! index of my current slip system
sLattice ! index of my current slip system as specified by lattice
sLattice ! index of my current slip system according to lattice order
real(pReal) gb, & ! short notation for G*b/2/pi
x, & ! coordinate in direction of lvec
y, & ! coordinate in direction of bvec
@ -758,6 +782,7 @@ real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstan
neighboring_rhoScrewExcess,&! screw excess dislocation density of my neighbor
neighboring_Nedge, & ! total number of edge excess dislocations in my neighbor
neighboring_Nscrew
logical flipConnectingVector
myInstance = phase_constitutionInstance(material_phase(g,ip,el))
@ -810,8 +835,35 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el))
neighboring_el = mesh_ipNeighborhood(1,n,ip,el)
neighboring_ip = mesh_ipNeighborhood(2,n,ip,el)
flipConnectingVector = .false.
if ( neighboring_el == 0 .or. neighboring_ip == 0 ) cycle
if ( neighboring_el == 0 .or. neighboring_ip == 0 ) then
if (.not. periodicBC) then
cycle
else
flipConnectingVector = .true.
select case (n)
case (1)
neighboring_el = mesh_ipNeighborhood(1,2,ip,el)
neighboring_ip = mesh_ipNeighborhood(2,2,ip,el)
case (2)
neighboring_el = mesh_ipNeighborhood(1,1,ip,el)
neighboring_ip = mesh_ipNeighborhood(2,1,ip,el)
case (3)
neighboring_el = mesh_ipNeighborhood(1,4,ip,el)
neighboring_ip = mesh_ipNeighborhood(2,4,ip,el)
case (4)
neighboring_el = mesh_ipNeighborhood(1,3,ip,el)
neighboring_ip = mesh_ipNeighborhood(2,3,ip,el)
case (5)
neighboring_el = mesh_ipNeighborhood(1,6,ip,el)
neighboring_ip = mesh_ipNeighborhood(2,6,ip,el)
case (6)
neighboring_el = mesh_ipNeighborhood(1,5,ip,el)
neighboring_ip = mesh_ipNeighborhood(2,5,ip,el)
endselect
endif
endif
! deformation gradients needed for mapping between configurations
neighboring_F = math_mul33x33(Fe(:,:,g,neighboring_ip,neighboring_el), Fp(:,:,g,neighboring_ip,neighboring_el))
@ -820,8 +872,13 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el))
neighboring_invFe = math_inv3x3(Fe(:,:,g,neighboring_ip,neighboring_el))
! calculate connection vector between me and my neighbor in its lattice configuration
connectingVector = math_mul33x3(neighboring_invFe, math_mul33x3(Favg, &
(mesh_ipCenterOfGravity(:,neighboring_ip,neighboring_el) - mesh_ipCenterOfGravity(:,ip,el)) ) )
if (flipConnectingVector) then
connectingVector = math_mul33x3(neighboring_invFe, math_mul33x3(Favg, &
-(mesh_ipCenterOfGravity(:,neighboring_ip,neighboring_el) - mesh_ipCenterOfGravity(:,ip,el)) ) )
else
connectingVector = math_mul33x3(neighboring_invFe, math_mul33x3(Favg, &
(mesh_ipCenterOfGravity(:,neighboring_ip,neighboring_el) - mesh_ipCenterOfGravity(:,ip,el)) ) )
endif
! neighboring dislocation densities
neighboring_rhoEdgePos = state(1, neighboring_ip, neighboring_el)%p( 1: ns)
@ -940,7 +997,7 @@ integer(pInt) myInstance, & ! curren
l, &
t, & ! dislocation type
s, & ! index of my current slip system
sLattice ! index of my current slip system as specified by lattice
sLattice ! index of my current slip system according to lattice order
real(pReal), dimension(6) :: Tdislocation_v ! dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress
real(pReal), dimension(3,3,3,3) :: dLp_dTstar3333 ! derivative of Lp with respect to Tstar (3x3x3x3 matrix)
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: &
@ -982,17 +1039,16 @@ forall (s =1:ns) &
!*** Calculation of gdot and its tangent
v = constitutive_nonlocal_v0PerSlipSystem(:,myInstance) &
* exp( - constitutive_nonlocal_G / ( kB * Temperature) * (1.0_pReal - (abs(tauSlip)/tauSlipThreshold) ) ) &
* exp( - constitutive_nonlocal_Q0(myInstance) / ( kB * Temperature) * (1.0_pReal - (abs(tauSlip)/tauSlipThreshold) ) ) &
* sign(1.0_pReal,tauSlip)
gdotSlip = sum(rho,2) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) * v
dgdot_dtauSlip = gdotSlip * constitutive_nonlocal_G / ( kB * Temperature * tauSlipThreshold )
dgdot_dtauSlip = gdotSlip * constitutive_nonlocal_Q0(myInstance) / ( kB * Temperature * tauSlipThreshold )
!*** Calculation of Lp and its tangent
do s = 1,ns
sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance)
Lp = Lp + gdotSlip(s) * lattice_Sslip(:,:,sLattice,myStructure)
@ -1004,17 +1060,20 @@ enddo
dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333)
! if (debugger) then
! !$OMP CRITICAL (write2out)
! write(6,*) '::: LpandItsTangent',g,ip,el
! write(6,*)
! 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),/)') '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
! !$OMPEND CRITICAL (write2out)
! endif
if (debugger) then
!$OMP CRITICAL (write2out)
write(6,*) '::: LpandItsTangent',g,ip,el
write(6,*)
write(6,'(a,/,3(3(f12.3,x)/))') 'Tdislocation / MPa', math_Mandel6to33(Tdislocation_v/1e6)
write(6,'(a,/,3(3(f12.3,x)/))') 'Tstar / MPa', math_Mandel6to33(Tstar_v/1e6)
write(6,'(a,/,12(f12.5,x),/)') 'tauSlip / MPa', tauSlip/1e6_pReal
! 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
call flush(6)
!$OMPEND CRITICAL (write2out)
endif
endsubroutine
@ -1087,7 +1146,8 @@ integer(pInt) myInstance, & ! current
c, & ! character of dislocation
n, & ! index of my current neighbor
t, & ! type of dislocation
s ! index of my current slip system
s, & ! index of my current slip system
sLattice ! index of my current slip system according to lattice order
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: &
rho, & ! dislocation densities (positive/negative screw and edge without dipoles)
totalRhoDot, & ! total rate of change of dislocation densities
@ -1106,10 +1166,10 @@ real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstan
rhoDip, & ! dipole dislocation densities (screw and edge dipoles)
totalRhoDipDot, & ! total rate of change of dipole dislocation densities
thisRhoDipDot, & ! rate of change of dipole dislocation densities for this mechanism
dDipMin, & ! minimum 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
dDipMaxDot ! rate of change of the maximum stable dipole distance for edges and screws
dLower, & ! minimum stable dipole distance for edges and screws
dUpper, & ! current maximum stable dipole distance for edges and screws
dUpper0, & ! maximum stable dipole distance for edges and screws at start of crystallite increment
dUpperDot ! rate of change of the maximum stable dipole distance for edges and screws
real(pReal), dimension(3,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: &
m ! direction of dislocation motion
real(pReal), dimension(3,3) :: F, & ! total deformation gradient
@ -1122,6 +1182,7 @@ real(pReal) norm_surfaceNormal, & ! euclidic
area, & ! area of the current interface
detFe, & ! determinant of elastic defornmation gradient
D ! self diffusion
logical flipAreaNormal
myInstance = phase_constitutionInstance(material_phase(g,ip,el))
myStructure = constitutive_nonlocal_structure(myInstance)
@ -1131,10 +1192,10 @@ tauSlip = 0.0_pReal
subTauSlip0 = 0.0_pReal
v = 0.0_pReal
gdot = 0.0_pReal
dDipMin = 0.0_pReal
dDipMax = 0.0_pReal
dDipMax0 = 0.0_pReal
dDipMaxDot = 0.0_pReal
dLower = 0.0_pReal
dUpper = 0.0_pReal
dUpper0 = 0.0_pReal
dUpperDot = 0.0_pReal
totalRhoDot = 0.0_pReal
thisRhoDot = 0.0_pReal
totalRhoDipDot = 0.0_pReal
@ -1154,16 +1215,15 @@ subTdislocation0_v = subState0(g,ip,el)%p(8*ns+1:8*ns+6)
!*** Calculate shear rate
do s = 1,ns ! loop over slip systems
tauSlip(s) = math_mul6x6( Tstar_v + Tdislocation_v, &
lattice_Sslip_v(:,constitutive_nonlocal_slipSystemLattice(s,myInstance),myStructure) )
sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance)
tauSlip(s) = math_mul6x6( Tstar_v + Tdislocation_v, lattice_Sslip_v(:,sLattice,myStructure) )
subTauSlip0(s) = math_mul6x6( subTstar0_v + subTdislocation0_v, lattice_Sslip_v(:,sLattice,myStructure) )
subTauSlip0(s) = math_mul6x6( subTstar0_v - subTdislocation0_v, &
lattice_Sslip_v(:,constitutive_nonlocal_slipSystemLattice(s,myInstance),myStructure) )
enddo
v = constitutive_nonlocal_v0PerSlipSystem(:,myInstance) &
* exp( - constitutive_nonlocal_G / ( kB * Temperature) * (1.0_pReal - (abs(tauSlip)/tauSlipThreshold) ) ) &
* exp( - constitutive_nonlocal_Q0(myInstance) / ( kB * Temperature) * (1.0_pReal - (abs(tauSlip)/tauSlipThreshold) ) ) &
* sign(1.0_pReal,tauSlip)
forall (t = 1:4) &
@ -1178,6 +1238,32 @@ if (debugger) then
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,/,(12(f12.5,x),/))') 'gdot total/ 1e-3', sum(gdot,2)*1e3_pReal
call flush(6)
!$OMP CRITICAL (write2out)
endif
!****************************************************************************
!*** calculate limits for stable dipole height and its rate of change
dLower(:,1) = constitutive_nonlocal_dLowerEdgePerSlipSystem(:,myInstance)
dLower(:,2) = constitutive_nonlocal_dLowerScrewPerSlipSystem(:,myInstance)
dUpper(:,2) = min( constitutive_nonlocal_Gmod(myInstance) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) &
/ ( 8.0_pReal * pi * abs(tauSlip) ), &
1.0_pReal / sqrt( sum(rho,2)+sum(rhoDip,2) ) )
dUpper(:,1) = dUpper(:,2) / ( 1.0_pReal - constitutive_nonlocal_nu(myInstance) )
dUpper0(:,2) = min( constitutive_nonlocal_Gmod(myInstance) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) &
/ ( 8.0_pReal * pi * abs(subTauSlip0) ), &
1.0_pReal / sqrt( sum(rho,2)+sum(rhoDip,2) ) )
dUpper0(:,1) = dUpper0(:,2) / ( 1.0_pReal - constitutive_nonlocal_nu(myInstance) )
if (subdt > 0) dUpperDot = (dUpper - dUpper0) / subdt
if (debugger) then
!$OMP CRITICAL (write2out)
write(6,'(a,/,2(12(e12.5,x),/))') 'dUpper:',dUpper
write(6,'(a,/,2(12(e12.5,x),/))') 'dUpperDot:',dUpperDot
call flush(6)
!$OMP CRITICAL (write2out)
endif
@ -1185,12 +1271,10 @@ endif
!****************************************************************************
!*** 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)
thisRhoDipDot = 0.0_pReal ! dipoles don't multiplicate
totalRhoDot = totalRhoDot + thisRhoDot
totalRhoDipDot = totalRhoDipDot + thisRhoDipDot
@ -1198,6 +1282,7 @@ totalRhoDipDot = totalRhoDipDot + thisRhoDipDot
if (debugger) then
!$OMP CRITICAL (write2out)
write(6,'(a,/,6(12(e12.5,x),/))') 'dislocation multiplication', thisRhoDot * subdt, thisRhoDipDot * subdt
call flush(6)
!$OMP CRITICAL (write2out)
endif
@ -1214,14 +1299,38 @@ m(:,:,3) = lattice_st(:, constitutive_nonlocal_slipSystemLattice(:,myInstance),
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))
detFe = math_det3x3(Fe(:,:,g,ip,el))
do n = 1,FE_NipNeighbors(mesh_element(2,el)) ! loop through my neighbors
neighboring_el = mesh_ipNeighborhood(1,n,ip,el)
neighboring_ip = mesh_ipNeighborhood(2,n,ip,el)
flipAreaNormal = .false.
! if neighbor exists, total deformation gradient is averaged over me and my neighbor
if (periodicBC .and. (neighboring_el == 0 .or. neighboring_ip == 0) ) then
flipAreaNormal = .true.
select case (n)
case (1)
neighboring_el = mesh_ipNeighborhood(1,2,ip,el)
neighboring_ip = mesh_ipNeighborhood(2,2,ip,el)
case (2)
neighboring_el = mesh_ipNeighborhood(1,1,ip,el)
neighboring_ip = mesh_ipNeighborhood(2,1,ip,el)
case (3)
neighboring_el = mesh_ipNeighborhood(1,4,ip,el)
neighboring_ip = mesh_ipNeighborhood(2,4,ip,el)
case (4)
neighboring_el = mesh_ipNeighborhood(1,3,ip,el)
neighboring_ip = mesh_ipNeighborhood(2,3,ip,el)
case (5)
neighboring_el = mesh_ipNeighborhood(1,6,ip,el)
neighboring_ip = mesh_ipNeighborhood(2,6,ip,el)
case (6)
neighboring_el = mesh_ipNeighborhood(1,5,ip,el)
neighboring_ip = mesh_ipNeighborhood(2,5,ip,el)
endselect
endif
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)
@ -1230,8 +1339,13 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el))
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)))
if (flipAreaNormal) then
surfaceNormal = math_det3x3(Favg) / detFe &
* math_mul33x3(transpose(Fe(:,:,g,ip,el)), math_mul33x3(Favg,-mesh_ipAreaNormal(:,n,ip,el)))
else
surfaceNormal = math_det3x3(Favg) / detFe &
* math_mul33x3(transpose(Fe(:,:,g,ip,el)), math_mul33x3(Favg,mesh_ipAreaNormal(:,n,ip,el)))
endif
norm_surfaceNormal = math_norm3(surfaceNormal)
surfaceNormal = surfaceNormal / norm_surfaceNormal
area = mesh_ipArea(n,ip,el) * norm_surfaceNormal
@ -1245,7 +1359,7 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el))
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
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 ...
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
!*****************************************************************************************************
@ -1265,6 +1379,7 @@ totalRhoDipDot = totalRhoDipDot + thisRhoDipDot
if (debugger) then
!$OMP CRITICAL (write2out)
write(6,'(a,/,6(12(e12.5,x),/))') 'dislocation flux', thisRhoDot * subdt, thisRhoDipDot * subdt
call flush(6)
!$OMP CRITICAL (write2out)
endif
@ -1272,36 +1387,14 @@ endif
!****************************************************************************
!*** calculate dipole formation and annihilation
!*** limits for stable dipole height and its tate of change
dDipMin(:,1) = constitutive_nonlocal_dDipMinEdgePerSlipSystem(:,myInstance)
dDipMin(:,2) = constitutive_nonlocal_dDipMinScrewPerSlipSystem(:,myInstance)
dDipMax(:,2) = min( constitutive_nonlocal_Gmod(myInstance) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) &
/ ( 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) )
dDipMax0(:,2) = min( constitutive_nonlocal_Gmod(myInstance) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) &
/ ( 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) )
dDipMaxDot(:,1) = (dDipMax(:,1) - dDipMax0(:,1)) / subdt
dDipMaxDot(:,2) = (dDipMax(:,2) - dDipMax0(:,2)) / subdt
! if (debugger) write(6,'(a,/,2(12(e12.5,x),/))') 'dDipMax:',dDipMax
! if (debugger) write(6,'(a,/,2(12(e12.5,x),/))') 'dDipMaxDot:',dDipMaxDot
!*** formation by glide
thisRhoDot = 0.0_pReal
thisRhoDipDot = 0.0_pReal
forall (c=1:2) &
thisRhoDipDot(:,c) = 4.0_pReal * dDipMax(:,c) / constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) &
thisRhoDipDot(:,c) = 4.0_pReal * dUpper(:,c) / constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) &
* ( 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)
thisRhoDot(:,t) = -0.5_pReal * thisRhoDipDot(:,(t-1)/2+1)
totalRhoDot = totalRhoDot + thisRhoDot
totalRhoDipDot = totalRhoDipDot + thisRhoDipDot
@ -1309,19 +1402,18 @@ 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
call flush(6)
!$OMP CRITICAL (write2out)
endif
!*** athermal annihilation
thisRhoDot = 0.0_pReal
thisRhoDipDot = 0.0_pReal
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
thisRhoDipDot(:,c) = - 4.0_pReal * dLower(:,c) / constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) &
* ( rho(:,2*c-1)*abs(gdot(:,2*c)) + rho(:,2*c)*abs(gdot(:,2*c-1)) & ! was single hitting single
+ 0.5_pReal * rhoDip(:,c) * (abs(gdot(:,2*c-1))+abs(gdot(:,2*c))) ) ! single knocks dipole constituent
thisRhoDot = 0.0_pReal ! singles themselves don't annihilate
totalRhoDot = totalRhoDot + thisRhoDot
totalRhoDipDot = totalRhoDipDot + thisRhoDipDot
@ -1329,22 +1421,22 @@ totalRhoDipDot = totalRhoDipDot + thisRhoDipDot
if (debugger) then
!$OMP CRITICAL (write2out)
write(6,'(a,/,6(12(e12.5,x),/))') 'athermal dipole annihilation', thisRhoDot * subdt, thisRhoDipDot * subdt
call flush(6)
!$OMP CRITICAL (write2out)
endif
!*** thermally activated annihilation
thisRhoDot = 0.0_pReal
thisRhoDipDot = 0.0_pReal
!*** thermally activated annihilation of dipoles
D = constitutive_nonlocal_D0(myInstance) * exp(-constitutive_nonlocal_Qsd(myInstance) / (kB * Temperature))
vClimb = constitutive_nonlocal_atomicVolume(myInstance) * D / ( kB * Temperature ) &
* constitutive_nonlocal_Gmod(myInstance) / ( 2.0_pReal * pi * (1.0_pReal-constitutive_nonlocal_nu(myInstance)) ) &
* 2.0_pReal / ( dDipMax(:,1) + dDipMin(:,1) )
* 2.0_pReal / ( dUpper(:,1) + dLower(:,1) )
thisRhoDipDot(:,1) = - 4.0_pReal * rho(:,1) * vClimb / ( dDipMax(:,1) - dDipMin(:,1) )
thisRhoDipDot(:,1) = - 4.0_pReal * rhoDip(:,1) * vClimb / ( dUpper(:,1) - dLower(:,1) ) ! edge climb
thisRhoDipDot(:,2) = 0.0_pReal !!! cross slipping still has to be implemented !!!
thisRhoDot = 0.0_pReal
totalRhoDot = totalRhoDot + thisRhoDot
totalRhoDipDot = totalRhoDipDot + thisRhoDipDot
@ -1352,30 +1444,46 @@ 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
call flush(6)
!$OMP CRITICAL (write2out)
endif
!*** formation by stress decrease = increase in dDipMax
! !!! MISSING !!!
!*** formation by stress change = alteration in dUpper
thisRhoDipDot = 0.0_pReal
!*** dipole dissociation by increased stress = decrease in dDipMax
forall (c=1:2, s=1:ns, dUpperDot(s,c) > 0) & ! stress decrease
thisRhoDipDot(s,c) = 8.0_pReal * rho(s,2*c-1) * rho(s,2*c) * dUpper(s,c) * dUpperDot(s,c)
forall (c=1:2, s=1:ns, dUpperDot(s,c) < 0) & ! increased stress
thisRhoDipDot(s,c) = rhoDip(s,c) * dUpperDot(s,c) / (dUpper(s,c) - dLower(s,c))
! !!! MISSING !!!
forall (t=1:4) &
thisRhoDot(:,t) = -0.5_pReal * thisRhoDipDot(:,(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 stability by stress change', thisRhoDot * subdt, thisRhoDipDot * subdt
call flush(6)
!$OMP CRITICAL (write2out)
endif
!****************************************************************************
!*** assign the rates of dislocation densities to my dotState
dotState(1,ip,el)%p(1:4*ns) = reshape(totalRhoDot,(/4*ns/))
dotState(1,ip,el)%p(1:4*ns) = reshape(totalRhoDot,(/4*ns/)) ! one-dimension only (linear list)
dotState(1,ip,el)%p(4*ns+1:6*ns) = reshape(totalRhoDipDot,(/2*ns/))
if (debugger) then
!$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
call flush(6)
!$OMP CRITICAL (write2out)
endif
@ -1419,12 +1527,13 @@ endfunction
!*********************************************************************
!* return array of constitutive results *
!*********************************************************************
pure function constitutive_nonlocal_postResults(Tstar_v, Temperature, dt, state, g, ip, el)
pure function constitutive_nonlocal_postResults(Tstar_v, subTstar0_v, Temperature, dt, subdt, state, subState0, dotState, g, ip, el)
use prec, only: pReal, &
pInt, &
p_vec
use math, only: math_mul6x6
use math, only: math_mul6x6, &
pi
use mesh, only: mesh_NcpElems, &
mesh_maxNips
use material, only: homogenization_maxNgrains, &
@ -1440,10 +1549,14 @@ integer(pInt), intent(in) :: g, & ! current grain nu
ip, & ! current integration point
el ! current element number
real(pReal), intent(in) :: dt, & ! time increment
subdt, & ! time increment of crystallite substep
Temperature ! temperature
real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola-Kirchhoff stress in Mandel notation
real(pReal), dimension(6), intent(in) :: Tstar_v, & ! 2nd Piola-Kirchhoff stress in Mandel notation
subTstar0_v ! 2nd Piola-Kirchhoff stress in Mandel notation at start of crystallite inc
type(p_vec), dimension(homogenization_maxNgrains, mesh_maxNips, mesh_NcpElems), intent(in) :: &
state ! microstructural state
state, & ! microstructural state
subState0, & ! microstructural state at start of crystallite inc
dotState ! evolution rate of microstructural state
!*** output variables
real(pReal), dimension(constitutive_nonlocal_sizePostResults(phase_constitutionInstance(material_phase(g,ip,el)))) :: &
@ -1455,96 +1568,228 @@ integer(pInt) myInstance, & ! current instance
ns, & ! short notation for the total number of active slip systems
o, & ! index of current output
s, & ! index of current slip system
sLattice, & ! index of current slip system as specified by lattice
c
real(pReal) tau, & ! resolved shear stress on current slip system
v ! dislocation velocity on current slip system
sLattice, & ! index of my current slip system according to lattice order
cs, & ! constitutive result index
c, & ! character of dislocation
t ! type of dislocation
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: &
rho, & ! dislocation densities (positive/negative screw and edge without dipoles)
rhoDot, & ! evolution rate of dislocation densities (positive/negative screw and edge without dipoles)
gdot ! shear rates
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: &
rhoForest, & ! forest dislocation density
tauSlipThreshold, & ! threshold shear stress
tauSlip, & ! current resolved shear stress
subTauSlip0, & ! resolved shear stress at start of crystallite increment
v, & ! dislocation velocity
invLambda, & ! inverse of mean free path for dislocations
vClimb ! climb velocity of edge dipoles
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),2) :: &
rhoDip, & ! dipole dislocation densities (screw and edge dipoles)
rhoDipDot, & ! evolution rate of dipole dislocation densities (screw and edge dipoles)
dLower, & ! minimum stable dipole distance for edges and screws
dUpper, & ! current maximum stable dipole distance for edges and screws
dUpper0, & ! maximum stable dipole distance for edges and screws at start of crystallite increment
dUpperDot ! rate of change of the maximum stable dipole distance for edges and screws
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
real(pReal) D ! self diffusion
myInstance = phase_constitutionInstance(material_phase(g,ip,el))
myStructure = constitutive_nonlocal_structure(myInstance)
ns = constitutive_nonlocal_totalNslip(myInstance)
c = 0_pInt
cs = 0_pInt
constitutive_nonlocal_postResults = 0.0_pReal
! short hand notations for state variables
forall (t = 1:4) rho(:,t) = state(g,ip,el)%p((t-1)*ns+1:t*ns)
forall (c = 1:2) rhoDip(:,c) = state(g,ip,el)%p((3+c)*ns+1:(4+c)*ns)
rhoForest = state(g,ip,el)%p(6*ns+1:7*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)
subTdislocation0_v = subState0(g,ip,el)%p(8*ns+1:8*ns+6)
forall (t = 1:4) rhoDot(:,t) = dotState(g,ip,el)%p((t-1)*ns+1:t*ns)
forall (c = 1:2) rhoDipDot(:,c) = dotState(g,ip,el)%p((3+c)*ns+1:(4+c)*ns)
! Calculate shear rate
do s = 1,ns
sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance)
tauSlip(s) = math_mul6x6( Tstar_v + Tdislocation_v, lattice_Sslip_v(:,sLattice,myStructure) )
subTauSlip0(s) = math_mul6x6( subTstar0_v + subTdislocation0_v, lattice_Sslip_v(:,sLattice,myStructure) )
enddo
v = constitutive_nonlocal_v0PerSlipSystem(:,myInstance) &
* exp( - constitutive_nonlocal_Q0(myInstance) / ( kB * Temperature) * (1.0_pReal - (abs(tauSlip)/tauSlipThreshold) ) ) &
* sign(1.0_pReal,tauSlip)
forall (t = 1:4) &
gdot(:,t) = rho(:,t) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) * v
! calculate limits for stable dipole height and its rate of change
dLower(:,1) = constitutive_nonlocal_dLowerEdgePerSlipSystem(:,myInstance)
dLower(:,2) = constitutive_nonlocal_dLowerScrewPerSlipSystem(:,myInstance)
dUpper(:,2) = min( constitutive_nonlocal_Gmod(myInstance) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) &
/ ( 8.0_pReal * pi * abs(tauSlip) ), &
1.0_pReal / sqrt( sum(rho,2)+sum(rhoDip,2) ) )
dUpper(:,1) = dUpper(:,2) / ( 1.0_pReal - constitutive_nonlocal_nu(myInstance) )
dUpper0(:,2) = min( constitutive_nonlocal_Gmod(myInstance) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) &
/ ( 8.0_pReal * pi * abs(subTauSlip0) ), &
1.0_pReal / sqrt( sum(rho,2)+sum(rhoDip,2) ) )
dUpper0(:,1) = dUpper0(:,2) / ( 1.0_pReal - constitutive_nonlocal_nu(myInstance) )
if (subdt > 0) then
dUpperDot = (dUpper - dUpper0) / subdt
else
dUpperDot = 0.0_pReal
endif
do o = 1,phase_Noutput(material_phase(g,ip,el))
select case(constitutive_nonlocal_output(o,myInstance))
case ('rho')
constitutive_nonlocal_postResults(c+1:c+ns) = state(g,ip,el)%p(1:ns) + state(g,ip,el)%p(ns+1:2*ns) &
+ state(g,ip,el)%p(2*ns+1:3*ns) + state(g,ip,el)%p(3*ns+1:4*ns)
c = c + ns
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(rho,2)
cs = cs + ns
case ('delta')
constitutive_nonlocal_postResults(cs+1:cs+ns) = 1.0_pReal / sqrt( sum(rho,2) )
cs = cs + ns
case ('rho_edge')
constitutive_nonlocal_postResults(c+1:c+ns) = state(g,ip,el)%p(1:ns) + state(g,ip,el)%p(ns+1:2*ns)
c = c + ns
constitutive_nonlocal_postResults(cs+1:cs+ns) = rho(:,1) + rho(:,2)
cs = cs + ns
case ('rho_screw')
constitutive_nonlocal_postResults(c+1:c+ns) = state(g,ip,el)%p(2*ns+1:3*ns) + state(g,ip,el)%p(3*ns+1:4*ns)
c = c + ns
constitutive_nonlocal_postResults(cs+1:cs+ns) = rho(:,3) + rho(:,4)
cs = cs + ns
case ('excess_rho')
constitutive_nonlocal_postResults(c+1:c+ns) = state(g,ip,el)%p(1:ns) - state(g,ip,el)%p(ns+1:2*ns) &
+ state(g,ip,el)%p(2*ns+1:3*ns) - state(g,ip,el)%p(3*ns+1:4*ns)
c = c + ns
constitutive_nonlocal_postResults(cs+1:cs+ns) = rho(:,1) - rho(:,2) + rho(:,3) - rho(:,4)
cs = cs + ns
case ('excess_rho_edge')
constitutive_nonlocal_postResults(c+1:c+ns) = state(g,ip,el)%p(1:ns) - state(g,ip,el)%p(ns+1:2*ns)
c = c + ns
constitutive_nonlocal_postResults(cs+1:cs+ns) = rho(:,1) - rho(:,2)
cs = cs + ns
case ('excess_rho_screw')
constitutive_nonlocal_postResults(c+1:c+ns) = state(g,ip,el)%p(2*ns+1:3*ns) - state(g,ip,el)%p(3*ns+1:4*ns)
c = c + ns
constitutive_nonlocal_postResults(cs+1:cs+ns) = rho(:,3) - rho(:,4)
cs = cs + ns
case ('rho_forest')
constitutive_nonlocal_postResults(c+1:c+ns) = state(g,ip,el)%p(6*ns+1:7*ns)
c = c + ns
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoForest
cs = cs + ns
case ('rho_dip')
constitutive_nonlocal_postResults(c+1:c+ns) = state(g,ip,el)%p(4*ns+1:5*ns) + state(g,ip,el)%p(5*ns+1:6*ns)
c = c + ns
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(rhoDip,2)
cs = cs + ns
case ('delta_dip')
constitutive_nonlocal_postResults(cs+1:cs+ns) = 1.0_pReal / sqrt( sum(rhoDip,2) )
cs = cs + ns
case ('rho_edge_dip')
constitutive_nonlocal_postResults(c+1:c+ns) = state(g,ip,el)%p(4*ns+1:5*ns)
c = c + ns
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoDip(:,1)
cs = cs + ns
case ('rho_screw_dip')
constitutive_nonlocal_postResults(c+1:c+ns) = state(g,ip,el)%p(5*ns+1:6*ns)
c = c + ns
constitutive_nonlocal_postResults(cs+1:cs+ns) = rhoDip(:,2)
cs = cs + ns
case ('shearrate')
do s = 1,ns
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) )
if (state(g,ip,el)%p(7*ns+s) > 0.0_pReal) then
v = constitutive_nonlocal_v0PerSlipSystem(s,myInstance) &
* exp( - constitutive_nonlocal_G / ( kB * Temperature) * (1.0_pReal - (abs(tau)/state(g,ip,el)%p(7*ns+s)) ) ) &
* sign(1.0_pReal,tau)
else
v = 0.0_pReal
endif
constitutive_nonlocal_postResults(c+s) = ( state(g,ip,el)%p(s) + state(g,ip,el)%p(ns+s) &
+ state(g,ip,el)%p(2*ns+s) + state(g,ip,el)%p(3*ns+s) ) &
* constitutive_nonlocal_burgersPerSlipSystem(s,myInstance) * v
enddo
c = c + ns
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(abs(gdot),2)
cs = cs + ns
case ('resolvedstress')
do s = 1,ns
sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance)
constitutive_nonlocal_postResults(c+s) = math_mul6x6( Tstar_v + state(g,ip,el)%p(8*ns+1:8*ns+6), &
lattice_Sslip_v(:,sLattice,myStructure) )
constitutive_nonlocal_postResults(cs+s) = math_mul6x6( Tstar_v + Tdislocation_v, lattice_Sslip_v(:,sLattice,myStructure) )
enddo
c = c + ns
cs = cs + ns
case ('resistance')
constitutive_nonlocal_postResults(c+1:c+ns) = state(g,ip,el)%p(7*ns+1:8*ns)
c = c + ns
constitutive_nonlocal_postResults(cs+1:cs+ns) = tauSlipThreshold
cs = cs + ns
case ('rho_dot')
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(rhoDot,2)
cs = cs + ns
case ('rho_dot_dip')
constitutive_nonlocal_postResults(cs+1:cs+ns) = sum(rhoDipDot,2)
cs = cs + ns
case ('rho_dot_gen')
invLambda = sqrt(rhoForest) / constitutive_nonlocal_lambda0PerSlipSystem(:,myInstance)
constitutive_nonlocal_postResults(cs+1:cs+ns) = &
sum(abs(gdot),2) * invLambda / constitutive_nonlocal_burgersPerSlipSystem(:,myInstance)
cs = cs + ns
case ('rho_dot_sgl2dip')
do c=1,2 ! dipole formation by glide
constitutive_nonlocal_postResults(cs+1:cs+ns) = constitutive_nonlocal_postResults(cs+1:cs+ns) + &
4.0_pReal * dUpper(:,c) / constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) &
* ( rho(:,2*c-1)*abs(gdot(:,2*c)) + rho(:,2*c)*abs(gdot(:,2*c-1)) )
enddo
do c=1,2
forall (s=1:ns, dUpperDot(s,c) > 0) & ! dipole formation by stress decrease
constitutive_nonlocal_postResults(cs+s) = constitutive_nonlocal_postResults(cs+s) + &
8.0_pReal * rho(s,2*c-1) * rho(s,2*c) * dUpper(s,c) * dUpperDot(s,c)
enddo
cs = cs + ns
case ('rho_dot_dip2sgl')
do c=1,2
forall (s=1:ns, dUpperDot(s,c) < 0) &
constitutive_nonlocal_postResults(cs+s) = constitutive_nonlocal_postResults(cs+s) - &
rhoDip(s,c) * dUpperDot(s,c) / (dUpper(s,c) - dLower(s,c))
enddo
cs = cs + ns
case ('rho_dot_ann_ath')
do c=1,2
constitutive_nonlocal_postResults(cs+1:cs+ns) = constitutive_nonlocal_postResults(cs+1:cs+ns) + &
4.0_pReal * dLower(:,c) / constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) &
* ( rho(:,2*c-1)*abs(gdot(:,2*c)) + rho(:,2*c)*abs(gdot(:,2*c-1)) &
+ 0.5_pReal * rhoDip(:,c) * (abs(gdot(:,2*c-1))+abs(gdot(:,2*c))) )
enddo
cs = cs + ns
case ('rho_dot_ann_the')
D = constitutive_nonlocal_D0(myInstance) * exp(-constitutive_nonlocal_Qsd(myInstance) / (kB * Temperature))
vClimb = constitutive_nonlocal_atomicVolume(myInstance) * D / ( kB * Temperature ) &
* constitutive_nonlocal_Gmod(myInstance) / ( 2.0_pReal * pi * (1.0_pReal-constitutive_nonlocal_nu(myInstance)) ) &
* 2.0_pReal / ( dUpper(:,1) + dLower(:,1) )
constitutive_nonlocal_postResults(cs+1:cs+ns) = 4.0_pReal * rhoDip(:,1) * vClimb / ( dUpper(:,1) - dLower(:,1) )
! !!! cross-slip of screws missing !!!
cs = cs + ns
case ('rho_dot_flux')
! !!! still has to be implemented !!!
constitutive_nonlocal_postResults(cs+1:cs+ns) = 0.0_pReal
cs = cs + ns
case ('d_upper_edge')
constitutive_nonlocal_postResults(cs+1:cs+ns) = dUpper(:,1)
cs = cs + ns
case ('d_upper_screw')
constitutive_nonlocal_postResults(cs+1:cs+ns) = dUpper(:,2)
cs = cs + ns
case ('d_upper_dot_edge')
constitutive_nonlocal_postResults(cs+1:cs+ns) = dUpperDot(:,1)
cs = cs + ns
case ('d_upper_dot_screw')
constitutive_nonlocal_postResults(cs+1:cs+ns) = dUpperDot(:,2)
cs = cs + ns
end select
enddo

View File

@ -350,7 +350,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
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 g = 1,myNgrains
debugger = (e == 1 .and. i == 1 .and. g == 1)
! debugger = (e == 1 .and. i == 1 .and. g == 1)
if (crystallite_converged(g,i,e)) then
if (debugger) then
!$OMP CRITICAL (write2out)
@ -358,6 +358,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
crystallite_subFrac(g,i,e),' to current crystallite_subfrac ', &
crystallite_subFrac(g,i,e)+crystallite_subStep(g,i,e),' in crystallite_stressAndItsTangent'
write(6,*)
call flush(6)
!$OMPEND CRITICAL (write2out)
endif
crystallite_subFrac(g,i,e) = crystallite_subFrac(g,i,e) + crystallite_subStep(g,i,e)
@ -388,6 +389,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
write(6,'(a78,f10.8)') 'cutback step in crystallite_stressAndItsTangent with new crystallite_subStep: ',&
crystallite_subStep(g,i,e)
write(6,*)
call flush(6)
!$OMPEND CRITICAL (write2out)
endif
endif
@ -452,7 +454,12 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
endif
enddo; enddo; enddo
!$OMPEND PARALLEL DO
write(6,*) count(crystallite_onTrack(1,:,:)),'IPs onTrack after preguess for state'
if (debugger) then
!$OMP CRITICAL (write2out)
write(6,*) count(crystallite_onTrack(1,:,:)),'IPs onTrack after preguess for state'
call flush(6)
!$OMPEND CRITICAL (write2out)
endif
! --+>> state loop <<+--
@ -482,9 +489,13 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
enddo
enddo
!$OMPEND PARALLEL DO
write(6,*) count(crystallite_onTrack(1,:,:)),'IPs onTrack after stress integration'
if (debugger) then
!$OMP CRITICAL (write2out)
write(6,*) count(crystallite_onTrack(1,:,:)),'IPs onTrack after stress integration'
call flush(6)
!$OMPEND CRITICAL (write2out)
endif
crystallite_todo = crystallite_todo .and. crystallite_onTrack
if (any(.not. crystallite_onTrack .and. .not. crystallite_localConstitution)) &
crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! all nonlocal crystallites can be skipped
@ -544,10 +555,15 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
if (any(.not. crystallite_onTrack .and. .not. crystallite_localConstitution)) &
crystallite_todo = crystallite_todo .and. crystallite_localConstitution ! all nonlocal crystallites can be skipped
write(6,*) count(crystallite_onTrack(1,:,:)),'IPs onTrack after state update'
write(6,*) count(crystallite_converged(1,:,:)),'IPs converged'
write(6,*) count(crystallite_todo(1,:,:)),'IPs todo'
write(6,*)
if (debugger) then
!$OMP CRITICAL (write2out)
write(6,*) count(crystallite_onTrack(1,:,:)),'IPs onTrack after state update'
write(6,*) count(crystallite_converged(1,:,:)),'IPs converged'
write(6,*) count(crystallite_todo(1,:,:)),'IPs todo'
write(6,*)
call flush(6)
!$OMPEND CRITICAL (write2out)
endif
enddo ! crystallite convergence loop
@ -612,6 +628,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
write (6,'(a,/,3(3(f12.8,x)/))') ' Fp of 1 1 1',myFp(1:3,:)
write (6,'(a,/,3(3(f12.8,x)/))') ' Lp of 1 1 1',myLp(1:3,:)
write (6,'(a,/,16(6(e12.4,x)/),2(f12.4,x))') 'state of 1 1 1',myState/1e6
call flush(6)
!$OMPEND CRITICAL (write2out)
endif
do k = 1,3 ! perturbation...
@ -624,6 +641,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
write (6,'(i1,x,i1)') k,l
write (6,*) '============='
write (6,'(a,/,3(3(f12.6,x)/))') 'pertF of 1 1 1',crystallite_subF(1:3,:,g,i,e)
call flush(6)
!$OMPEND CRITICAL (write2out)
endif
onTrack = .true.
@ -634,8 +652,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
onTrack = crystallite_integrateStress(g,i,e) ! stress of perturbed situation (overwrites _P,_Tstar_v,_Fp,_Lp,_Fe)
if (onTrack) then
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_Temperature(g,i,e), crystallite_subdt(g,i,e), g, i, e)
crystallite_Fe, crystallite_Fp, crystallite_Temperature(g,i,e), &
crystallite_subdt(g,i,e), g, i, e)
stateConverged = crystallite_updateState(g,i,e) ! update state
temperatureConverged = crystallite_updateTemperature(g,i,e) ! update temperature
@ -649,6 +667,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
write (6,'(a,/,3(3(f12.4,x)/))') 'DP of 1 1 1',(crystallite_P(1:3,:,g,i,e)-myP(1:3,:))/1e6
write (6,'(a,/,16(6(e12.4,x)/),/,2(f12.4,x))') 'state of 1 1 1',constitutive_state(g,i,e)%p/1e6
write (6,'(a,/,16(6(e12.4,x)/),/,2(f12.4,x))') 'Dstate of 1 1 1',(constitutive_state(g,i,e)%p-myState)/1e6
call flush(6)
!$OMPEND CRITICAL (write2out)
endif
enddo
@ -743,6 +762,7 @@ endsubroutine
if (debugger) then
!$OMP CRITICAL (write2out)
write(6,*) '::: updateState encountered NaN',g,i,e
call flush(6)
!$OMPEND CRITICAL (write2out)
endif
return
@ -768,6 +788,7 @@ endsubroutine
write(6,*)
write(6,'(a,/,12(f12.5,x))') 'resid tolerance',abs(residuum/rTol_crystalliteState/constitutive_state(g,i,e)%p(1:mySize))
write(6,*)
call flush(6)
!$OMPEND CRITICAL (write2out)
endif
return
@ -825,6 +846,7 @@ endsubroutine
crystallite_updateTemperature = .false. ! indicate update failed
!$OMP CRITICAL (write2out)
write(6,*) '::: updateTemperature encountered NaN',g,i,e
call flush(6)
!$OMPEND CRITICAL (write2out)
return
endif
@ -952,6 +974,7 @@ endsubroutine
write(6,*) '::: integrateStress failed on invFp_current inversion',g,i,e
write(6,*)
write(6,'(a11,3(i3,x),/,3(3(f12.7,x)/))') 'invFp_new at ',g,i,e,invFp_new
call flush(6)
!$OMPEND CRITICAL (write2out)
endif
return
@ -981,6 +1004,7 @@ LpLoop: do
!$OMP CRITICAL (write2out)
write(6,*) '::: integrateStress reached loop limit',g,i,e
write(6,*)
call flush(6)
!$OMPEND CRITICAL (write2out)
endif
return
@ -1009,6 +1033,7 @@ LpLoop: do
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
call flush(6)
!$OMPEND CRITICAL (write2out)
endif
@ -1030,6 +1055,7 @@ LpLoop: do
if (debugger) then
!$OMP CRITICAL (write2out)
write(6,*) '::: integrateStress encountered NaN at iteration', NiterationStress,'at',g,i,e
call flush(6)
!$OMPEND CRITICAL (write2out)
endif
return
@ -1069,6 +1095,7 @@ LpLoop: do
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(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
call flush(6)
!$OMPEND CRITICAL (write2out)
endif
return
@ -1099,6 +1126,7 @@ LpLoop: do
write(6,*) '::: integrateStress failed on invFp_new inversion at iteration', NiterationStress
write(6,*)
write(6,'(a11,3(i3,x),/,3(3(f12.7,x)/))') 'invFp_new at ',g,i,e,invFp_new
call flush(6)
!$OMPEND CRITICAL (write2out)
endif
return
@ -1127,6 +1155,7 @@ LpLoop: do
write(6,'(a,/,3(3(f12.7,x)/))') 'P / MPa',crystallite_P(:,:,g,i,e)/1e6
write(6,'(a,/,3(3(f12.7,x)/))') 'Lp',crystallite_Lp(:,:,g,i,e)
write(6,'(a,/,3(3(f12.7,x)/))') 'Fp',crystallite_Fp(:,:,g,i,e)
call flush(6)
!$OMP CRITICAL (write2out)
endif
@ -1145,8 +1174,6 @@ LpLoop: do
! return results of particular grain
!********************************************************************
function crystallite_postResults(&
Tstar_v,& ! stress
Temperature, & ! temperature
dt,& ! time increment
g,& ! grain number
i,& ! integration point number
@ -1171,16 +1198,13 @@ function crystallite_postResults(&
integer(pInt), intent(in):: e, & ! element index
i, & ! integration point index
g ! grain index
real(pReal), intent(in):: Temperature, & ! temperature
dt ! time increment
real(pReal), dimension(6), intent(in):: Tstar_v ! 2nd Piola-Kirchhoff stress in Mandel notation
real(pReal), intent(in):: dt ! time increment
!*** output variables ***!
real(pReal), dimension(crystallite_Nresults + constitutive_sizePostResults(g,i,e)) :: crystallite_postResults
!*** local variables ***!
real(pReal), dimension(3,3) :: U, &
R
real(pReal), dimension(3,3) :: U, R
logical error
if (crystallite_Nresults >= 2) then
@ -1198,7 +1222,9 @@ function crystallite_postResults(&
endif
crystallite_postResults(crystallite_Nresults+1:crystallite_Nresults+constitutive_sizePostResults(g,i,e)) = &
constitutive_postResults(Tstar_v,Temperature,dt,g,i,e)
constitutive_postResults(crystallite_Tstar_v(:,g,i,e), crystallite_subTstar0_v(:,g,i,e), crystallite_Temperature(g,i,e), &
dt, crystallite_subdt(g,i,e), g, i, e)
return
endfunction

View File

@ -188,6 +188,7 @@ endsubroutine
enddo
write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_MaterialpointLoopDistribution)
write(6,*)
call flush(6)
endsubroutine

View File

@ -236,6 +236,7 @@ subroutine materialpoint_stressAndItsTangent(&
write (6,'(a,/,3(3(f12.7,x)/))') 'F of 1 1',materialpoint_F(1:3,:,1,1)
write (6,'(a,/,3(3(f12.7,x)/))') 'Fp0 of 1 1 1',crystallite_Fp0(1:3,:,1,1,1)
write (6,'(a,/,3(3(f12.7,x)/))') 'Lp0 of 1 1 1',crystallite_Lp0(1:3,:,1,1,1)
call flush(6)
!$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
@ -284,6 +285,7 @@ subroutine materialpoint_stressAndItsTangent(&
write(6,'(a21,f10.8,a34,f10.8,a37,/)') 'winding forward from ', &
materialpoint_subFrac(i,e), ' to current materialpoint_subFrac ', &
materialpoint_subFrac(i,e)+materialpoint_subStep(i,e),' in materialpoint_stressAndItsTangent'
call flush(6)
!$OMPEND CRITICAL (write2out)
endif
@ -320,6 +322,7 @@ subroutine materialpoint_stressAndItsTangent(&
!$OMP CRITICAL (write2out)
write(6,'(a82,f10.8,/)') 'cutback step in materialpoint_stressAndItsTangent with new materialpoint_subStep: ',&
materialpoint_subStep(i,e)
call flush(6)
!$OMPEND CRITICAL (write2out)
endif
@ -460,7 +463,7 @@ subroutine materialpoint_postResults(dt)
use mesh, only: mesh_element
use material, only: homogenization_Ngrains
use constitutive, only: constitutive_sizePostResults, constitutive_postResults
use crystallite
use crystallite, only: crystallite_Nresults, crystallite_postResults
implicit none
real(pReal), intent(in) :: dt
@ -476,13 +479,13 @@ subroutine materialpoint_postResults(dt)
materialpoint_results(c+1,i,e) = d; c = c+1_pInt ! tell size of homogenization results
if (d > 0_pInt) then ! any homogenization results to mention?
materialpoint_results(c+1:c+d,i,e) = & ! tell homogenization results
homogenization_postResults(i,e); c = c+d
homogenization_postResults(i,e); c = c+d
endif
do g = 1,myNgrains !
d = crystallite_Nresults+constitutive_sizePostResults(g,i,e)
d = crystallite_Nresults + constitutive_sizePostResults(g,i,e)
materialpoint_results(c+1,i,e) = d; c = c+1_pInt ! tell size of crystallite results
materialpoint_results(c+1:c+d,i,e) = & ! tell crystallite results
crystallite_postResults(crystallite_Tstar_v(:,g,i,e),crystallite_Temperature(g,i,e),dt,g,i,e); c = c+d
crystallite_postResults(dt,g,i,e); c = c+d
enddo
enddo
enddo

View File

@ -111,17 +111,32 @@ constitution nonlocal
/nonlocal/
(output) rho
(output) delta
(output) rho_edge
(output) rho_screw
(output) excess_rho_edge
(output) excess_rho_screw
(output) rho_forest
(output) rho_dip
(output) delta_dip
(output) rho_edge_dip
(output) rho_screw_dip
(output) shearrate
(output) resolvedstress
(output) resistance
(output) rho_dot
(output) rho_dot_dip
(output) rho_dot_gen
(output) rho_dot_sgl2dip
(output) rho_dot_dip2sgl
(output) rho_dot_ann_ath
(output) rho_dot_ann_the
(output) rho_dot_flux
(output) d_upper_edge
(output) d_upper_screw
(output) d_upper_dot_screw
(output) d_upper_dot_edge
lattice_structure fcc
Nslip 12 0 0 0 # number of slip systems per family
@ -138,12 +153,13 @@ rhoScrewNeg0 1e10 0 0 0 # Initial negative screw disloc
rhoEdgeDip0 1 0 0 0 # Initial edge dipole dislocation density in m/m**3
rhoScrewDip0 1 0 0 0 # Initial screw dipole dislocation density in m/m**3
v0 1e-4 0 0 0 # prefactor for dislocation velocity
Q0 3e-19 # activation energy for dislocation glide
dDipMinEdge 1e-9 0 0 0 # minimum distance for stable edge dipoles in m
dDipMinScrew 1e-9 0 0 0 # minimum distance for stable screw dipoles in m
lambda0 100 0 0 0 # prefactor for mean free path
atomicVolume 1.7e-29
D0 1e-4 # prefactor for self-diffusion coefficient
Qsd 2.3e-19 # activation energy for dislocation climb
Qsd 2.3e-19 # activation enthalpy for seld-diffusion
relevantRho 1e3 # dislocation density considered relevant
interaction_SlipSlip 1.0 2.2 3.0 1.6 3.8 4.5 # Dislocation interaction coefficient