fixed bug in constitutive_nonlocal_kinetics that was introduced in rev2088 along with the non-Schmid behavior: nonSchmid matrix of last slip system was used for all slip systems; led to extremely bad convergence due to flawed dLp_dT

This commit is contained in:
Christoph Kords 2013-05-10 22:29:12 +00:00
parent b83b557af2
commit 6fac7c347c
1 changed files with 35 additions and 27 deletions

View File

@ -1706,7 +1706,8 @@ integer(pInt) myInstance, & ! curren
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(3,3,2) :: nonSchmid_tensor
real(pReal), dimension(3,3,2,constitutive_nonlocal_totalNslip(phase_plasticityInstance(material_phase(g,ip,el)))) :: &
nonSchmidTensor
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),8) :: &
rhoSgl ! single dislocation densities (including blocked)
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_plasticityInstance(material_phase(g,ip,el))),4) :: &
@ -1724,6 +1725,7 @@ real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_plasticityInstance
Lp = 0.0_pReal
dLp_dTstar3333 = 0.0_pReal
nonSchmidTensor = 0.0_pReal
myInstance = phase_plasticityInstance(material_phase(g,ip,el))
myStructure = constitutive_nonlocal_structure(myInstance)
@ -1743,40 +1745,43 @@ where (abs(rhoSgl) * mesh_ipVolume(ip,el) ** 0.667_pReal < constitutive_nonlocal
!*** get effective resolved shear stress
!*** add non schmid contributions to ONLY screw components if present (i.e. if NnonSchmid(myStructure) > 0)
do s = 1_pInt,ns
tau(s,1:4) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,constitutive_nonlocal_slipSystemLattice(s,myInstance),myStructure)) &
+ tauBack(s)
!*** adding non schmid contributions to ONLY screw components if present (i.e. if NnonSchmid(myStructure) > 0)
nonSchmid_tensor(1:3,1:3,1) = &
math_Mandel6to33(lattice_Sslip_v(1:6,1,constitutive_nonlocal_slipSystemLattice(s,myInstance),myStructure))
nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,1)
sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance)
tau(s,1:4) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,myStructure)) + tauBack(s)
nonSchmidTensor(1:3,1:3,1,s) = lattice_Sslip(1:3,1:3,sLattice,myStructure)
nonSchmidTensor(1:3,1:3,2,s) = nonSchmidTensor(1:3,1:3,1,s)
do k = 1_pInt, NnonSchmid(myStructure)
tau(s,3) = tau(s,3) + constitutive_nonlocal_nonSchmidCoeff(k,myInstance)* &
math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,2*k,constitutive_nonlocal_slipSystemLattice(s,myInstance),myStructure))
tau(s,4) = tau(s,4) + constitutive_nonlocal_nonSchmidCoeff(k,myInstance)* &
math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,2*k+1,constitutive_nonlocal_slipSystemLattice(s,myInstance),myStructure))
nonSchmid_tensor(1:3,1:3,1) = nonSchmid_tensor(1:3,1:3,1) + constitutive_nonlocal_nonSchmidCoeff(k,myInstance)*&
math_Mandel6to33(lattice_Sslip_v(1:6,2*k,constitutive_nonlocal_slipSystemLattice(s,myInstance),myStructure))
nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,2) + constitutive_nonlocal_nonSchmidCoeff(k,myInstance)*&
math_Mandel6to33(lattice_Sslip_v(1:6,2*k+1,constitutive_nonlocal_slipSystemLattice(s,myInstance),myStructure))
tau(s,3) = tau(s,3) + constitutive_nonlocal_nonSchmidCoeff(k,myInstance) &
* math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,2*k,sLattice,myStructure))
tau(s,4) = tau(s,4) + constitutive_nonlocal_nonSchmidCoeff(k,myInstance) &
* math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,2*k+1,sLattice,myStructure))
nonSchmidTensor(1:3,1:3,1,s) = nonSchmidTensor(1:3,1:3,1,s) &
+ constitutive_nonlocal_nonSchmidCoeff(k,myInstance) &
* math_Mandel6to33(lattice_Sslip_v(1:6,2*k,sLattice,myStructure))
nonSchmidTensor(1:3,1:3,2,s) = nonSchmidTensor(1:3,1:3,2,s) &
+ constitutive_nonlocal_nonSchmidCoeff(k,myInstance) &
* math_Mandel6to33(lattice_Sslip_v(1:6,2*k+1,sLattice,myStructure))
enddo
enddo
!*** get dislocation velocity and its tangent and store the velocity in the state array
if (myStructure == 1_pInt .and. NnonSchmid(myStructure) == 0_pInt) then ! for fcc all velcities are equal
call constitutive_nonlocal_kinetics(v(1:ns,1), tau(1:ns,1), 1_pInt, Temperature, state, g, ip, el, dv_dtau(1:ns,1))
if (myStructure == 1_pInt .and. NnonSchmid(myStructure) == 0_pInt) then ! for fcc all velcities are equal
call constitutive_nonlocal_kinetics(v(1:ns,1), tau(1:ns,1), 1_pInt, Temperature, state, &
g, ip, el, dv_dtau(1:ns,1))
do t = 1_pInt,4_pInt
v(1:ns,t) = v(1:ns,1)
dv_dtau(1:ns,t) = dv_dtau(1:ns,1)
state%p((13_pInt+t)*ns+1:(14_pInt+t)*ns) = v(1:ns,1)
enddo
else ! for all other lattice structures the velcities may vary with character and sign
else ! for all other lattice structures the velocities may vary with character and sign
do t = 1_pInt,4_pInt
c = (t-1_pInt)/2_pInt+1_pInt
call constitutive_nonlocal_kinetics(v(1:ns,t), tau(1:ns,t), c, Temperature, state, g, ip, el, dv_dtau(1:ns,t))
call constitutive_nonlocal_kinetics(v(1:ns,t), tau(1:ns,t), c, Temperature, state, &
g, ip, el, dv_dtau(1:ns,t))
state%p((13+t)*ns+1:(14+t)*ns) = v(1:ns,t)
enddo
endif
@ -1795,22 +1800,25 @@ if (constitutive_nonlocal_deadZoneScaling(myInstance)) then
forall(s = 1_pInt:ns, sum(abs(rhoSgl(s,1:8))) > 0.0_pReal) &
deadZoneSize(s) = maxval(abs(rhoSgl(s,5:8)) / (rhoSgl(s,1:4) + abs(rhoSgl(s,5:8))))
endif
gdotTotal = sum(rhoSgl(1:ns,1:4) * v, 2) * constitutive_nonlocal_burgers(1:ns,myInstance) * (1.0_pReal - deadZoneSize)
gdotTotal = sum(rhoSgl(1:ns,1:4) * v, 2) * constitutive_nonlocal_burgers(1:ns,myInstance) &
* (1.0_pReal - deadZoneSize)
do t = 1_pInt,4_pInt
dgdot_dtau(:,t) = rhoSgl(1:ns,t) * dv_dtau(1:ns,t) * constitutive_nonlocal_burgers(1:ns,myInstance) * (1.0_pReal - deadZoneSize)
dgdot_dtau(:,t) = rhoSgl(1:ns,t) * dv_dtau(1:ns,t) * constitutive_nonlocal_burgers(1:ns,myInstance) &
* (1.0_pReal - deadZoneSize)
enddo
!*** Calculation of Lp and its tangent
do s = 1_pInt,ns
sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance)
Lp = Lp + gdotTotal(s) * lattice_Sslip(1:3,1:3,sLattice,myStructure)
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) + &
dgdot_dtau(s,1)*lattice_Sslip(i,j,sLattice,myStructure)*lattice_Sslip(k,l,sLattice,myStructure) +&
dgdot_dtau(s,2)*lattice_Sslip(i,j,sLattice,myStructure)*lattice_Sslip(k,l,sLattice,myStructure) +&
dgdot_dtau(s,3)*lattice_Sslip(i,j,sLattice,myStructure)*nonSchmid_tensor(k,l,1) +&
dgdot_dtau(s,4)*lattice_Sslip(i,j,sLattice,myStructure)*nonSchmid_tensor(k,l,2)
dLp_dTstar3333(i,j,k,l) = dLp_dTstar3333(i,j,k,l) &
+ dgdot_dtau(s,1) * lattice_Sslip(i,j,sLattice,myStructure) * lattice_Sslip(k,l,sLattice,myStructure) &
+ dgdot_dtau(s,2) * lattice_Sslip(i,j,sLattice,myStructure) * lattice_Sslip(k,l,sLattice,myStructure) &
+ dgdot_dtau(s,3) * lattice_Sslip(i,j,sLattice,myStructure) * nonSchmidTensor(k,l,1,s) &
+ dgdot_dtau(s,4) * lattice_Sslip(i,j,sLattice,myStructure) * nonSchmidTensor(k,l,2,s)
enddo
dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333)
@ -1823,7 +1831,7 @@ dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333)
write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_LpandItsTangent at el ip g ',el,ip,g
write(6,*)
write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> gdot total / 1e-3',gdotTotal*1e3_pReal
write(6,'(a,/,3(12x,3(f12.7,1x),/))') '<< CONST >> Lp',Lp
write(6,'(a,/,3(12x,3(f12.7,1x),/))') '<< CONST >> Lp',transpose(Lp)
endif
#endif