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

use parameter structure and avoid conversion 33<->6

This commit is contained in:
Martin Diehl 2019-02-17 12:15:46 +01:00
parent 0f319e2cf6
commit cf32e7d1f5
3 changed files with 90 additions and 146 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
src/constitutive.f90
View File

@ -402,15 +402,15 @@ end subroutine constitutive_microstructure
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient
! ToDo: Discuss wheter it makes sense if crystallite handles the configuration conversion, i.e.
! Mp in, dLp_dMp out
!--------------------------------------------------------------------------------------------------
subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, el)
subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
S, Fi, ipc, ip, el)
use prec, only: &
pReal
use math, only: &
math_mul33x33, &
math_6toSym33, &
math_sym33to6, &
math_99to3333
math_mul33x33
use material, only: &
phasememberAt, &
phase_plasticity, &
@ -444,9 +444,8 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(6) :: &
S6 !< 2nd Piola-Kirchhoff stress (vector notation)
real(pReal), intent(in), dimension(3,3) :: &
S, & !< 2nd Piola-Kirchhoff stress
Fi !< intermediate deformation gradient
real(pReal), intent(out), dimension(3,3) :: &
Lp !< plastic velocity gradient
@ -455,11 +454,8 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e
dLp_dFi !< derivative of Lp with respect to Fi
real(pReal), dimension(3,3,3,3) :: &
dLp_dMp !< derivative of Lp with respect to Mandel stress
real(pReal), dimension(9,9) :: &
dLp_dMp99 !< derivative of Lp with respect to Mstar (matrix notation)
real(pReal), dimension(3,3) :: &
Mp, & !< Mandel stress work conjugate with Lp
S !< 2nd Piola-Kirchhoff stress
Mp !< Mandel stress work conjugate with Lp
integer(pInt) :: &
ho, & !< homogenization
tme !< thermal member position
@ -469,7 +465,6 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e
ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el)
S = math_6toSym33(S6)
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S)
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
@ -491,12 +486,11 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, e
case (PLASTICITY_KINEHARDENING_ID) plasticityType
of = phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMp, Mp,instance,of)
call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMp,Mp,instance,of)
case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMp99, math_sym33to6(Mp), &
call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMp,Mp, &
temperature(ho)%p(tme),ip,el)
dLp_dMp = math_99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
case (PLASTICITY_DISLOTWIN_ID) plasticityType
of = phasememberAt(ipc,ip,el)
@ -993,7 +987,7 @@ end subroutine constitutive_collectDeltaState
!--------------------------------------------------------------------------------------------------
!> @brief returns array of constitutive results
!--------------------------------------------------------------------------------------------------
function constitutive_postResults(S6, Fi, FeArray, ipc, ip, el)
function constitutive_postResults(S, Fi, FeArray, ipc, ip, el)
use prec, only: &
pReal
use math, only: &
@ -1058,8 +1052,8 @@ function constitutive_postResults(S6, Fi, FeArray, ipc, ip, el)
Fi !< intermediate deformation gradient
real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%Nelems) :: &
FeArray !< elastic deformation gradient
real(pReal), intent(in), dimension(6) :: &
S6 !< 2nd Piola Kirchhoff stress (vector notation)
real(pReal), intent(in), dimension(3,3) :: &
S !< 2nd Piola Kirchhoff stress
real(pReal), dimension(3,3) :: &
Mp !< Mandel stress
integer(pInt) :: &
@ -1067,11 +1061,11 @@ function constitutive_postResults(S6, Fi, FeArray, ipc, ip, el)
integer(pInt) :: &
ho, & !< homogenization
tme, & !< thermal member position
s, of, instance !< counter in source loop
i, of, instance !< counter in source loop
constitutive_postResults = 0.0_pReal
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_6toSym33(S6))
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S)
ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el)
@ -1112,13 +1106,13 @@ function constitutive_postResults(S6, Fi, FeArray, ipc, ip, el)
case (PLASTICITY_NONLOCAL_ID) plasticityType
constitutive_postResults(startPos:endPos) = &
plastic_nonlocal_postResults (S6,FeArray,ip,el)
plastic_nonlocal_postResults (Mp,FeArray,ip,el)
end select plasticityType
SourceLoop: do s = 1_pInt, phase_Nsources(material_phase(ipc,ip,el))
SourceLoop: do i = 1_pInt, phase_Nsources(material_phase(ipc,ip,el))
startPos = endPos + 1_pInt
endPos = endPos + sourceState(material_phase(ipc,ip,el))%p(s)%sizePostResults
sourceType: select case (phase_source(s,material_phase(ipc,ip,el)))
endPos = endPos + sourceState(material_phase(ipc,ip,el))%p(i)%sizePostResults
sourceType: select case (phase_source(i,material_phase(ipc,ip,el)))
case (SOURCE_damage_isoBrittle_ID) sourceType
constitutive_postResults(startPos:endPos) = source_damage_isoBrittle_postResults(ipc, ip, el)
case (SOURCE_damage_isoDuctile_ID) sourceType

25
src/crystallite.f90
View File

@ -286,7 +286,7 @@ subroutine crystallite_init
crystallite_outputID(o,c) = grainrotation_ID
case ('eulerangles') outputName
crystallite_outputID(o,c) = eulerangles_ID
case ('defgrad','f') outputName
case ('defgrad','f') outputName ! ToDo: no alias (f only)
crystallite_outputID(o,c) = defgrad_ID
case ('fe') outputName
crystallite_outputID(o,c) = fe_ID
@ -298,15 +298,15 @@ subroutine crystallite_init
crystallite_outputID(o,c) = lp_ID
case ('li') outputName
crystallite_outputID(o,c) = li_ID
case ('p','firstpiola','1stpiola') outputName
case ('p','firstpiola','1stpiola') outputName ! ToDo: no alias (p only)
crystallite_outputID(o,c) = p_ID
case ('s','tstar','secondpiola','2ndpiola') outputName
case ('s','tstar','secondpiola','2ndpiola') outputName ! ToDo: no alias (s only)
crystallite_outputID(o,c) = s_ID
case ('elasmatrix') outputName
crystallite_outputID(o,c) = elasmatrix_ID
case ('neighboringip') outputName
case ('neighboringip') outputName ! ToDo: this is not a result, it is static. Should be written out by mesh
crystallite_outputID(o,c) = neighboringip_ID
case ('neighboringelement') outputName
case ('neighboringelement') outputName ! ToDo: this is not a result, it is static. Should be written out by mesh
crystallite_outputID(o,c) = neighboringelement_ID
case default outputName
call IO_error(105_pInt,ext_msg=trim(str(o))//' (Crystallite)')
@ -426,7 +426,7 @@ end subroutine crystallite_init
!--------------------------------------------------------------------------------------------------
!> @brief calculate stress (P)
!--------------------------------------------------------------------------------------------------
function crystallite_stress(a)
function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
use prec, only: &
tol_math_check, &
dNeq0
@ -462,14 +462,11 @@ function crystallite_stress(a)
sourceState, &
phase_Nsources, &
phaseAt, phasememberAt
use constitutive, only: &
constitutive_SandItsTangents, &
constitutive_LpAndItsTangents, &
constitutive_LiAndItsTangents
implicit none
logical, dimension(theMesh%elem%nIPs,theMesh%Nelems) :: crystallite_stress
real(pReal), intent(in), optional :: a !ToDo: for some reason this prevents an internal compiler error in GNU. Very strange
real(pReal), intent(in), optional :: &
dummyArgumentToPreventInternalCompilerErrorWithGCC
real(pReal) :: &
formerSubStep
integer(pInt) :: &
@ -793,7 +790,7 @@ subroutine crystallite_stressTangent()
endif
call constitutive_LpAndItsTangents(devNull,dLpdS,dLpdFi, &
crystallite_Tstar_v(1:6,c,i,e), &
math_6toSym33(crystallite_Tstar_v(1:6,c,i,e)), &
crystallite_Fi(1:3,1:3,c,i,e),c,i,e) ! call constitutive law to calculate Lp tangent in lattice configuration
dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS
@ -1078,7 +1075,7 @@ function crystallite_postResults(ipc, ip, el)
c = c + 1_pInt
if (size(crystallite_postResults)-c > 0_pInt) &
crystallite_postResults(c+1:size(crystallite_postResults)) = &
constitutive_postResults(crystallite_Tstar_v(1:6,ipc,ip,el), crystallite_Fi(1:3,1:3,ipc,ip,el), &
constitutive_postResults(math_6toSym33(crystallite_Tstar_v(1:6,ipc,ip,el)), crystallite_Fi(1:3,1:3,ipc,ip,el), &
crystallite_Fe, ipc, ip, el)
end function crystallite_postResults
@ -1289,7 +1286,7 @@ logical function integrateStress(&
!* calculate plastic velocity gradient and its tangent from constitutive law
call constitutive_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, &
math_sym33to6(S), Fi_new, ipc, ip, el)
S, Fi_new, ipc, ip, el)
#ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &

169
src/plastic_nonlocal.f90
View File

@ -107,14 +107,11 @@ module plastic_nonlocal
rhoDotMultiplicationOutput, &
rhoDotSingle2DipoleGlideOutput, &
rhoDotAthermalAnnihilationOutput, &
rhoDotThermalAnnihilationOutput, &
nonSchmidProjection !< combined projection of Schmid and non-Schmid contributions to the resolved shear stress (only for screws)
rhoDotThermalAnnihilationOutput !< combined projection of Schmid and non-Schmid contributions to the resolved shear stress (only for screws)
real(pReal), dimension(:,:,:,:,:,:), allocatable, private :: &
compatibility !< slip system compatibility between me and my neighbors
real(pReal), dimension(:,:), allocatable, private :: &
nonSchmidCoeff
logical, dimension(:), allocatable, private :: &
shortRangeStressCorrection, & !< flag indicating the use of the short range stress correction by a excess density gradient term
@ -275,7 +272,6 @@ use IO, only: IO_read, &
IO_intValue, &
IO_warning, &
IO_error, &
IO_timeStamp, &
IO_EOF
use debug, only: debug_level, &
debug_constitutive, &
@ -319,7 +315,6 @@ integer(pInt) :: phase, &
c, & ! index of dislocation character
Nchunks_SlipSlip = 0_pInt, &
Nchunks_SlipFamilies = 0_pInt, &
Nchunks_nonSchmid = 0_pInt, &
mySize = 0_pInt ! to suppress warnings, safe as init is called only once
character(len=65536) :: &
tag = '', &
@ -396,7 +391,6 @@ allocate(lambda0PerSlipFamily(lattice_maxNslipFamily,maxNinstances), s
allocate(interactionSlipSlip(lattice_maxNinteraction,maxNinstances), source=0.0_pReal)
allocate(minDipoleHeightPerSlipFamily(lattice_maxNslipFamily,2,maxNinstances), source=-1.0_pReal)
allocate(peierlsStressPerSlipFamily(lattice_maxNslipFamily,2,maxNinstances), source=0.0_pReal)
allocate(nonSchmidCoeff(lattice_maxNnonSchmid,maxNinstances), source=0.0_pReal)
rewind(fileUnit)
@ -417,7 +411,6 @@ allocate(nonSchmidCoeff(lattice_maxNnonSchmid,maxNinstances), s
if (phase_plasticity(phase) == PLASTICITY_NONLOCAL_ID) then
Nchunks_SlipFamilies = count(lattice_NslipSystem(:,phase) > 0_pInt)
Nchunks_SlipSlip = maxval(lattice_InteractionSlipSlip(:,:,phase))
Nchunks_nonSchmid = lattice_NnonSchmid(phase)
endif
cycle
endif
@ -539,12 +532,6 @@ allocate(nonSchmidCoeff(lattice_maxNnonSchmid,maxNinstances), s
fEdgeMultiplication(instance) = IO_floatValue(line,chunkPos,2_pInt)
case('shortrangestresscorrection')
shortRangeStressCorrection(instance) = IO_floatValue(line,chunkPos,2_pInt) > 0.0_pReal
case ('nonschmid_coefficients')
if (chunkPos(1) < 1_pInt + Nchunks_nonSchmid) &
call IO_warning(52_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_NONLOCAL_label//')')
do f = 1_pInt,Nchunks_nonSchmid
nonSchmidCoeff(f,instance) = IO_floatValue(line,chunkPos,1_pInt+f)
enddo
case('probabilisticmultiplication','randomsources','randommultiplication','discretesources')
probabilisticMultiplication(instance) = IO_floatValue(line,chunkPos,2_pInt) > 0.0_pReal
end select
@ -686,7 +673,6 @@ allocate(compatibility(2,maxTotalNslip,maxTotalNslip,theMesh%elem%nIPneighbors,t
source=0.0_pReal)
allocate(peierlsStress(maxTotalNslip,2,maxNinstances), source=0.0_pReal)
allocate(colinearSystem(maxTotalNslip,maxNinstances), source=0_pInt)
allocate(nonSchmidProjection(3,3,4,maxTotalNslip,maxNinstances), source=0.0_pReal)
initializeInstances: do phase = 1_pInt, size(phase_plasticity)
NofMyPhase=count(material_phase==phase)
@ -856,19 +842,6 @@ allocate(nonSchmidProjection(3,3,4,maxTotalNslip,maxNinstances),
!* 3) negative screw at positive resolved stress
!* 4) negative screw at negative resolved stress
do s = 1_pInt,ns
do l = 1_pInt,lattice_NnonSchmid(phase)
nonSchmidProjection(1:3,1:3,1,s,instance) = nonSchmidProjection(1:3,1:3,1,s,instance) &
+ nonSchmidCoeff(l,instance) * lattice_Sslip(1:3,1:3,2*l,slipSystemLattice(s,instance),phase)
nonSchmidProjection(1:3,1:3,2,s,instance) = nonSchmidProjection(1:3,1:3,2,s,instance) &
+ nonSchmidCoeff(l,instance) * lattice_Sslip(1:3,1:3,2*l+1,slipSystemLattice(s,instance),phase)
enddo
nonSchmidProjection(1:3,1:3,3,s,instance) = -nonSchmidProjection(1:3,1:3,2,s,instance)
nonSchmidProjection(1:3,1:3,4,s,instance) = -nonSchmidProjection(1:3,1:3,1,s,instance)
forall (t = 1:4) &
nonSchmidProjection(1:3,1:3,t,s,instance) = nonSchmidProjection(1:3,1:3,t,s,instance) &
+ lattice_Sslip(1:3,1:3,1,slipSystemLattice(s,instance),phase)
enddo
call plastic_nonlocal_aTolState(phase,instance)
endif myPhase2
@ -1322,7 +1295,9 @@ real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(1_pI
ph = phaseAt(1,ip,el)
of = phasememberAt(1,ip,el)
instance = phase_plasticityInstance(ph)
ns = totalNslip(instance)
associate(prm => param(instance))
ns = prm%totalNslip
!*** get basic states
@ -1334,11 +1309,11 @@ endforall
forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) &
rhoDip(s,c) = max(plasticState(ph)%state(iRhoD(s,c,instance),of), 0.0_pReal) ! ensure positive dipole densities
where (abs(rhoSgl) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance) &
.or. abs(rhoSgl) < significantRho(instance)) &
where (abs(rhoSgl) * mesh_ipVolume(ip,el) ** 0.667_pReal < prm%significantN &
.or. abs(rhoSgl) < prm%significantRho) &
rhoSgl = 0.0_pReal
where (abs(rhoDip) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance) &
.or. abs(rhoDip) < significantRho(instance)) &
where (abs(rhoDip) * mesh_ipVolume(ip,el) ** 0.667_pReal < prm%significantN &
.or. abs(rhoDip) < prm%significantRho) &
rhoDip = 0.0_pReal
!*** calculate the forest dislocation density
@ -1377,7 +1352,7 @@ forall (s = 1_pInt:ns) &
tauBack = 0.0_pReal
if (.not. phase_localPlasticity(ph) .and. shortRangeStressCorrection(instance)) then
if (.not. phase_localPlasticity(ph) .and. prm%shortRangeStressCorrection) then
invFe = math_inv33(Fe)
invFp = math_inv33(Fp)
rhoExcess(1,1:ns) = rhoSgl(1:ns,1) - rhoSgl(1:ns,2)
@ -1418,10 +1393,9 @@ if (.not. phase_localPlasticity(ph) .and. shortRangeStressCorrection(instance))
math_mul33x3(invFe, mesh_ipCoordinates(1:3,neighbor_ip,neighbor_el) &
- mesh_ipCoordinates(1:3,ip,el))
normal_latticeConf = math_mul33x3(transpose(invFp), mesh_ipAreaNormal(1:3,n,ip,el))
if (math_mul3x3(normal_latticeConf,connection_latticeConf(1:3,n)) < 0.0_pReal) then ! neighboring connection points in opposite direction to face normal: must be periodic image
if (math_mul3x3(normal_latticeConf,connection_latticeConf(1:3,n)) < 0.0_pReal) & ! neighboring connection points in opposite direction to face normal: must be periodic image
connection_latticeConf(1:3,n) = normal_latticeConf * mesh_ipVolume(ip,el) &
/ mesh_ipArea(n,ip,el) ! instead take the surface normal scaled with the diameter of the cell
endif
else
! different number of active slip systems
call IO_error(-1_pInt,ext_msg='different number of active slip systems in neighboring IPs of same crystal structure')
@ -1507,7 +1481,7 @@ plasticState(ph)%state(iTauB(1:ns,instance),of) = tauBack
write(6,'(a,/,12x,12(f10.5,1x),/)') '<< CONST >> tauBack / MPa', tauBack*1e-6
endif
#endif
end associate
end subroutine plastic_nonlocal_microstructure
@ -1671,7 +1645,7 @@ end subroutine plastic_nonlocal_kinetics
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dTstar99, Tstar_v, Temperature, ip, el)
subroutine plastic_nonlocal_LpAndItsTangent(Lp, dLp_dMp, Mp, Temperature, ip, el)
use math, only: math_3333to99, &
math_mul6x6, &
@ -1687,9 +1661,6 @@ use material, only: material_phase, &
plasticState, &
phaseAt, phasememberAt,&
phase_plasticityInstance
use lattice, only: lattice_Sslip, &
lattice_Sslip_v, &
lattice_NnonSchmid
use mesh, only: mesh_ipVolume
implicit none
@ -1698,12 +1669,12 @@ implicit none
integer(pInt), intent(in) :: ip, & !< current integration point
el !< current element number
real(pReal), intent(in) :: Temperature !< temperature
real(pReal), dimension(6), intent(in) :: Tstar_v !< 2nd Piola-Kirchhoff stress in Mandel notation
real(pReal), dimension(3,3), intent(in) :: Mp
!*** output variables
real(pReal), dimension(3,3), intent(out) :: Lp !< plastic velocity gradient
real(pReal), dimension(9,9), intent(out) :: dLp_dTstar99 !< derivative of Lp with respect to Tstar (9x9 matrix)
real(pReal), dimension(3,3,3,3), intent(out) :: dLp_dMp !< derivative of Lp with respect to Tstar (9x9 matrix)
!*** local variables
integer(pInt) instance, & !< current instance of this plasticity
@ -1717,7 +1688,6 @@ integer(pInt) instance, &
t, & !< dislocation type
s, & !< index of my current slip system
sLattice !< index of my current slip system according to lattice order
real(pReal), dimension(3,3,3,3) :: dLp_dTstar3333 !< derivative of Lp with respect to Tstar (3x3x3x3 matrix)
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el))),8) :: &
rhoSgl !< single dislocation densities (including blocked)
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el))),4) :: &
@ -1737,7 +1707,7 @@ of = phasememberAt(1_pInt,ip,el)
!*** initialize local variables
Lp = 0.0_pReal
dLp_dTstar3333 = 0.0_pReal
dLp_dMp = 0.0_pReal
instance = phase_plasticityInstance(ph)
associate(prm => param(instance))
@ -1762,16 +1732,15 @@ tauThreshold = plasticState(ph)%state(iTauF(1:ns,instance),of)
!*** for screws possible non-schmid contributions are also taken into account
do s = 1_pInt,ns
sLattice = slipSystemLattice(s,instance)
tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,ph))
tau(s) = math_mul33xx33(Mp, prm%Schmid(1:3,1:3,s))
tauNS(s,1) = tau(s)
tauNS(s,2) = tau(s)
if (tau(s) > 0.0_pReal) then
tauNS(s,3) = math_mul33xx33(math_6toSym33(Tstar_v), nonSchmidProjection(1:3,1:3,1,s,instance))
tauNS(s,4) = math_mul33xx33(math_6toSym33(Tstar_v), nonSchmidProjection(1:3,1:3,3,s,instance))
tauNS(s,3) = math_mul33xx33(Mp, +prm%nonSchmid_pos(1:3,1:3,s))
tauNS(s,4) = math_mul33xx33(Mp, -prm%nonSchmid_neg(1:3,1:3,s))
else
tauNS(s,3) = math_mul33xx33(math_6toSym33(Tstar_v), nonSchmidProjection(1:3,1:3,2,s,instance))
tauNS(s,4) = math_mul33xx33(math_6toSym33(Tstar_v), nonSchmidProjection(1:3,1:3,4,s,instance))
tauNS(s,3) = math_mul33xx33(Mp, +prm%nonSchmid_neg(1:3,1:3,s))
tauNS(s,4) = math_mul33xx33(Mp, -prm%nonSchmid_pos(1:3,1:3,s))
endif
enddo
forall (t = 1_pInt:4_pInt) &
@ -1790,7 +1759,7 @@ dv_dtau(1:ns,2) = dv_dtau(1:ns,1)
dv_dtauNS(1:ns,2) = dv_dtauNS(1:ns,1)
!screws
if (lattice_NnonSchmid(ph) == 0_pInt) then ! no non-Schmid contributions
if (size(prm%nonSchmidCoeff) == 0_pInt) then ! no non-Schmid contributions
forall(t = 3_pInt:4_pInt)
v(1:ns,t) = v(1:ns,1)
dv_dtau(1:ns,t) = dv_dtau(1:ns,1)
@ -1817,47 +1786,37 @@ forall (s = 1_pInt:ns, t = 5_pInt:8_pInt, rhoSgl(s,t) * v(s,t-4_pInt) < 0.0_pRea
!*** Calculation of Lp and its tangent
gdotTotal = sum(rhoSgl(1:ns,1:4) * v, 2) * burgers(1:ns,instance)
gdotTotal = sum(rhoSgl(1:ns,1:4) * v, 2) * prm%burgers(1:ns)
do s = 1_pInt,ns
sLattice = slipSystemLattice(s,instance)
Lp = Lp + gdotTotal(s) * lattice_Sslip(1:3,1:3,1,sLattice,ph)
Lp = Lp + gdotTotal(s) * prm%Schmid(1:3,1:3,s)
! Schmid contributions to tangent
forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt,k=1_pInt:3_pInt,l=1_pInt:3_pInt) &
dLp_dTstar3333(i,j,k,l) = dLp_dTstar3333(i,j,k,l) &
+ lattice_Sslip(i,j,1,sLattice,ph) * lattice_Sslip(k,l,1,sLattice,ph) &
* sum(rhoSgl(s,1:4) * dv_dtau(s,1:4)) * burgers(s,instance)
dLp_dMp(i,j,k,l) = dLp_dMp(i,j,k,l) &
+ prm%Schmid(i,j,s) * prm%Schmid(k,l,s) &
* sum(rhoSgl(s,1:4) * dv_dtau(s,1:4)) * prm%burgers(s)
! non Schmid contributions to tangent
if (tau(s) > 0.0_pReal) then
forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt,k=1_pInt:3_pInt,l=1_pInt:3_pInt) &
dLp_dTstar3333(i,j,k,l) = dLp_dTstar3333(i,j,k,l) &
+ lattice_Sslip(i,j,1,sLattice,ph) &
* ( nonSchmidProjection(k,l,1,s,instance) * rhoSgl(s,3) * dv_dtauNS(s,3) &
+ nonSchmidProjection(k,l,3,s,instance) * rhoSgl(s,4) * dv_dtauNS(s,4) ) &
* burgers(s,instance)
dLp_dMp(i,j,k,l) = dLp_dMp(i,j,k,l) &
+ prm%Schmid(i,j,s) &
* ( prm%nonSchmid_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) &
- prm%nonSchmid_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4) ) &
* prm%burgers(s)
else
forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt,k=1_pInt:3_pInt,l=1_pInt:3_pInt) &
dLp_dTstar3333(i,j,k,l) = dLp_dTstar3333(i,j,k,l) &
+ lattice_Sslip(i,j,1,sLattice,ph) &
* ( nonSchmidProjection(k,l,2,s,instance) * rhoSgl(s,3) * dv_dtauNS(s,3) &
+ nonSchmidProjection(k,l,4,s,instance) * rhoSgl(s,4) * dv_dtauNS(s,4) ) &
* burgers(s,instance)
dLp_dMp(i,j,k,l) = dLp_dMp(i,j,k,l) &
+ prm%Schmid(i,j,s) &
* ( prm%nonSchmid_neg(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) &
- prm%nonSchmid_pos(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4) ) &
* prm%burgers(s)
endif
enddo
dLp_dTstar99 = math_3333to99(dLp_dTstar3333)
#ifdef DEBUG
if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt &
.and. ((debug_e == el .and. debug_i == ip)&
.or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then
write(6,'(/,a,i8,1x,i2,1x,i1,/)') '<< CONST >> nonlocal_LpandItsTangent at el ip',el,ip
write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> gdot total',gdotTotal
write(6,'(a,/,3(12x,3(f12.7,1x),/))') '<< CONST >> Lp',transpose(Lp)
endif
#endif
end associate
end subroutine plastic_nonlocal_LpAndItsTangent
@ -1866,7 +1825,7 @@ end subroutine plastic_nonlocal_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief (instantaneous) incremental change of microstructure
!--------------------------------------------------------------------------------------------------
subroutine plastic_nonlocal_deltaState(Tstar_v,ip,el)
subroutine plastic_nonlocal_deltaState(Mp,ip,el)
use prec, only: &
dNeq0
use debug, only: debug_level, &
@ -1877,9 +1836,8 @@ use debug, only: debug_level, &
debug_i, &
debug_e
use math, only: pi, &
math_mul6x6
use lattice, only: lattice_Sslip_v ,&
lattice_mu, &
math_mul33xx33
use lattice, only: lattice_mu, &
lattice_nu
use mesh, only: mesh_ipVolume
use material, only: material_phase, &
@ -1890,7 +1848,7 @@ use material, only: material_phase, &
implicit none
integer(pInt), intent(in) :: ip, & ! current grain number
el ! current element number
real(pReal), dimension(6), intent(in) :: Tstar_v ! current 2nd Piola-Kirchhoff stress in Mandel notation
real(pReal), dimension(3,3), intent(in) :: Mp !< MandelStress
integer(pInt) :: &
@ -1931,6 +1889,7 @@ real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,e
ph = phaseAt(1,ip,el)
of = phasememberAt(1,ip,el)
instance = phase_plasticityInstance(ph)
associate(prm => param(instance))
ns = totalNslip(instance)
@ -1980,8 +1939,7 @@ enddo
!*** calculate limits for stable dipole height
do s = 1_pInt,ns
sLattice = slipSystemLattice(s,instance)
tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,ph)) + tauBack(s)
tau(s) = math_mul33xx33(Mp, prm%Schmid(1:3,1:3,s)) + tauBack(s)
if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal
enddo
dLower = minDipoleHeight(1:ns,1:2,instance)
@ -2042,13 +2000,14 @@ forall (s = 1:ns, c = 1_pInt:2_pInt) &
write(6,'(a,/,10(12x,12(e12.5,1x),/),/)') '<< CONST >> dipole dissociation by stress increase', deltaRhoDipole2SingleStress
endif
#endif
end associate
end subroutine plastic_nonlocal_deltaState
!---------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!---------------------------------------------------------------------------------------------------
subroutine plastic_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature, &
subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
timestep,subfrac, ip,el)
use, intrinsic :: &
IEEE_arithmetic
@ -2065,9 +2024,9 @@ use debug, only: debug_level, &
debug_g, &
debug_i, &
debug_e
use math, only: math_mul6x6, &
math_mul3x3, &
use math, only: math_mul3x3, &
math_mul33x3, &
math_mul33xx33, &
math_mul33x33, &
math_inv33, &
math_det33, &
@ -2086,8 +2045,7 @@ use material, only: homogenization_maxNgrains, &
phaseAt, phasememberAt, &
phase_plasticity ,&
PLASTICITY_NONLOCAL_ID
use lattice, only: lattice_Sslip_v, &
lattice_sd, &
use lattice, only: lattice_sd, &
lattice_st ,&
lattice_mu, &
lattice_nu, &
@ -2102,7 +2060,7 @@ integer(pInt), intent(in) :: ip, &
el !< current element number
real(pReal), intent(in) :: Temperature, & !< temperature
timestep !< substepped crystallite time increment
real(pReal), dimension(6), intent(in) :: Tstar_v !< current 2nd Piola-Kirchhoff stress in Mandel notation
real(pReal), dimension(3,3), intent(in) :: Mp !< MandelStress
real(pReal), dimension(homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%nElems), intent(in) :: &
subfrac !< fraction of timestep at the beginning of the substepped crystallite time increment
real(pReal), dimension(3,3,homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%nElems), intent(in) :: &
@ -2198,6 +2156,7 @@ logical considerEnteringFlux, &
ph = material_phase(1_pInt,ip,el)
instance = phase_plasticityInstance(ph)
associate(prm => param(instance))
ns = totalNslip(instance)
tau = 0.0_pReal
@ -2271,8 +2230,7 @@ forall (t = 1_pInt:4_pInt) &
!*** calculate limits for stable dipole height
do s = 1_pInt,ns ! loop over slip systems
sLattice = slipSystemLattice(s,instance)
tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,ph)) + tauBack(s)
tau(s) = math_mul33xx33(Mp, prm%Schmid(1:3,1:3,s)) + tauBack(s)
if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal
enddo
@ -2661,7 +2619,7 @@ else
forall (s = 1:ns) &
plasticState(p)%dotState(iGamma(s,instance),o) = sum(gdot(s,1:4))
endif
end associate
end subroutine plastic_nonlocal_dotState
@ -2831,13 +2789,12 @@ end subroutine plastic_nonlocal_updateCompatibility
!--------------------------------------------------------------------------------------------------
!> @brief return array of constitutive results
!--------------------------------------------------------------------------------------------------
function plastic_nonlocal_postResults(Tstar_v,Fe,ip,el)
function plastic_nonlocal_postResults(Mp,Fe,ip,el)
use prec, only: &
dNeq0
use math, only: &
math_mul6x6, &
math_mul33x3, &
math_mul33x33, &
math_mul33xx33, &
pi
use mesh, only: &
theMesh
@ -2848,7 +2805,6 @@ function plastic_nonlocal_postResults(Tstar_v,Fe,ip,el)
plasticState, &
phase_plasticityInstance
use lattice, only: &
lattice_Sslip_v, &
lattice_sd, &
lattice_st, &
lattice_sn, &
@ -2856,8 +2812,7 @@ function plastic_nonlocal_postResults(Tstar_v,Fe,ip,el)
lattice_nu
implicit none
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), dimension(3,3), intent(in) :: Mp !< MandelStress
real(pReal), dimension(3,3,homogenization_maxNgrains,theMesh%elem%nIPs,theMesh%nElems), intent(in) :: &
Fe !< elastic deformation gradient
integer(pInt), intent(in) :: &
@ -2910,7 +2865,7 @@ ns = totalNslip(instance)
cs = 0_pInt
plastic_nonlocal_postResults = 0.0_pReal
associate(prm => param(instance))
!* short hand notations for state variables
forall (s = 1_pInt:ns, t = 1_pInt:4_pInt)
@ -2937,8 +2892,7 @@ forall (t = 1_pInt:4_pInt) &
!* calculate limits for stable dipole height
do s = 1_pInt,ns
sLattice = slipSystemLattice(s,instance)
tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,ph)) + tauBack(s)
tau(s) = math_mul33xx33(Mp, prm%Schmid(1:3,1:3,s)) + tauBack(s)
if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal
enddo
@ -3029,8 +2983,7 @@ outputsLoop: do o = 1_pInt,size(param(instance)%outputID)
case (resolvedstress_external_ID)
do s = 1_pInt,ns
sLattice = slipSystemLattice(s,instance)
plastic_nonlocal_postResults(cs+s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,ph))
plastic_nonlocal_postResults(cs+s) = math_mul33xx33(Mp, prm%Schmid(1:3,1:3,s))
enddo
cs = cs + ns
@ -3053,7 +3006,7 @@ outputsLoop: do o = 1_pInt,size(param(instance)%outputID)
case (rho_dot_gen_ID) ! Obsolete
plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotMultiplicationOutput(1:ns,1,1_pInt,ip,el) &
+ rhoDotMultiplicationOutput(1:ns,2,1_pInt,ip,el)
+ rhoDotMultiplicationOutput(1:ns,2,1_pInt,ip,el)
cs = cs + ns
case (rho_dot_gen_edge_ID)
@ -3074,7 +3027,7 @@ outputsLoop: do o = 1_pInt,size(param(instance)%outputID)
case (rho_dot_ann_ath_ID)
plastic_nonlocal_postResults(cs+1_pInt:cs+ns) = rhoDotAthermalAnnihilationOutput(1:ns,1,1_pInt,ip,el) &
+ rhoDotAthermalAnnihilationOutput(1:ns,2,1_pInt,ip,el)
+ rhoDotAthermalAnnihilationOutput(1:ns,2,1_pInt,ip,el)
cs = cs + ns
case (rho_dot_ann_the_edge_ID)
@ -3133,7 +3086,7 @@ outputsLoop: do o = 1_pInt,size(param(instance)%outputID)
end select
enddo outputsLoop
end associate
end function plastic_nonlocal_postResults
end module plastic_nonlocal