From 6fac7c347c64550c03f70b2abb3173999d6a6552 Mon Sep 17 00:00:00 2001 From: Christoph Kords Date: Fri, 10 May 2013 22:29:12 +0000 Subject: [PATCH] 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 --- code/constitutive_nonlocal.f90 | 62 +++++++++++++++++++--------------- 1 file changed, 35 insertions(+), 27 deletions(-) diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index 504db8d6b..89fa645b3 100644 --- a/code/constitutive_nonlocal.f90 +++ b/code/constitutive_nonlocal.f90 @@ -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