From 60fec0e8ec861a3fa49f09920736e4eabbfc1b23 Mon Sep 17 00:00:00 2001
From: Pratheek Shanthraj
Date: Mon, 21 Jan 2013 23:11:16 +0000
Subject: [PATCH] =?UTF-8?q?added=20code=20structure=20for=20non-schmid=20m?=
=?UTF-8?q?echanics.=20work=20in=20progress=E2=80=A6?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
---
code/constitutive_dislotwin.f90 | 12 ++--
code/constitutive_nonlocal.f90 | 10 +--
code/constitutive_phenopowerlaw.f90 | 102 +++++++++++++++++++++-------
code/constitutive_titanmod.f90 | 2 +-
code/lattice.f90 | 56 +++++++++++++--
5 files changed, 139 insertions(+), 43 deletions(-)
diff --git a/code/constitutive_dislotwin.f90 b/code/constitutive_dislotwin.f90
index 598a77625..f87e14c34 100644
--- a/code/constitutive_dislotwin.f90
+++ b/code/constitutive_dislotwin.f90
@@ -1000,7 +1000,7 @@ do f = 1_pInt,lattice_maxNslipFamily ! loop over
!* Calculation of Lp
!* Resolved shear stress on slip system
- tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure))
+ tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,myStructure))
!* Stress ratios
StressRatio_p = (abs(tau_slip(j))/state(g,ip,el)%p(5*ns+3*nt+j))**constitutive_dislotwin_p(myInstance)
@@ -1197,7 +1197,7 @@ do f = 1_pInt,lattice_maxNslipFamily ! loop over
!* Resolved shear stress on slip system
- tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure))
+ tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,myStructure))
!* Stress ratios
StressRatio_p = (abs(tau_slip(j))/state(g,ip,el)%p(5_pInt*ns+3_pInt*nt+j))**&
constitutive_dislotwin_p(myInstance)
@@ -1436,7 +1436,7 @@ do o = 1_pInt,phase_Noutput(material_phase(g,ip,el))
j = j + 1_pInt
!* Resolved shear stress on slip system
- tau = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure))
+ tau = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,myStructure))
!* Stress ratios
StressRatio_p = (abs(tau)/state(g,ip,el)%p(5_pInt*ns+3_pInt*nt+j))**&
constitutive_dislotwin_p(myInstance)
@@ -1466,7 +1466,7 @@ do o = 1_pInt,phase_Noutput(material_phase(g,ip,el))
do i = 1_pInt,constitutive_dislotwin_Nslip(f,myInstance) ! process each (active) slip system in family
j = j + 1_pInt
constitutive_dislotwin_postResults(c+j) =&
- dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure))
+ dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,myStructure))
enddo; enddo
c = c + ns
case ('threshold_stress_slip')
@@ -1481,7 +1481,7 @@ do o = 1_pInt,phase_Noutput(material_phase(g,ip,el))
j = j + 1_pInt
constitutive_dislotwin_postResults(c+j) = &
(3.0_pReal*constitutive_dislotwin_Gmod(myInstance)*constitutive_dislotwin_burgersPerSlipSystem(f,myInstance))/&
- (16.0_pReal*pi*abs(dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure))))
+ (16.0_pReal*pi*abs(dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,myStructure))))
constitutive_dislotwin_postResults(c+j) = min(constitutive_dislotwin_postResults(c+j),state(g,ip,el)%p(4*ns+2*nt+j))
! constitutive_dislotwin_postResults(c+j) = max(constitutive_dislotwin_postResults(c+j),state(g,ip,el)%p(4*ns+2*nt+j))
enddo; enddo
@@ -1563,7 +1563,7 @@ do o = 1_pInt,phase_Noutput(material_phase(g,ip,el))
j = j + 1_pInt
!* Resolved shear stress on slip system
- tau = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure))
+ tau = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,myStructure))
!* Stress ratios
StressRatio_p = (abs(tau)/state(g,ip,el)%p(5_pInt*ns+3_pInt*nt+j))**&
constitutive_dislotwin_p(myInstance)
diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90
index 06fc31c3e..d185c3bdd 100644
--- a/code/constitutive_nonlocal.f90
+++ b/code/constitutive_nonlocal.f90
@@ -1714,7 +1714,7 @@ where (abs(rhoSgl) * mesh_ipVolume(ip,el) ** 0.667_pReal < constitutive_nonlocal
!*** get effective resolved shear stress
do s = 1_pInt,ns
- tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(:,constitutive_nonlocal_slipSystemLattice(s,myInstance),myStructure)) &
+ tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(:,1,constitutive_nonlocal_slipSystemLattice(s,myInstance),myStructure)) &
+ tauBack(s)
enddo
@@ -1912,7 +1912,7 @@ enddo
do s = 1_pInt,ns
sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance)
- tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,sLattice,myStructure)) + tauBack(s)
+ tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,myStructure)) + tauBack(s)
if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal
enddo
dLower = constitutive_nonlocal_minimumDipoleHeight(1:ns,1:2,myInstance)
@@ -2185,7 +2185,7 @@ forall (t = 1_pInt:4_pInt) &
do s = 1_pInt,ns ! loop over slip systems
sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance)
- tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,sLattice,myStructure)) + tauBack(s)
+ tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,myStructure)) + tauBack(s)
if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal
enddo
@@ -3246,7 +3246,7 @@ forall (t = 1_pInt:4_pInt) &
do s = 1_pInt,ns
sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance)
- tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,sLattice,myStructure)) + tauBack(s)
+ tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,myStructure)) + tauBack(s)
if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal
enddo
@@ -3429,7 +3429,7 @@ do o = 1_pInt,phase_Noutput(material_phase(g,ip,el))
case ('resolvedstress_external')
do s = 1_pInt,ns
sLattice = constitutive_nonlocal_slipSystemLattice(s,myInstance)
- constitutive_nonlocal_postResults(cs+s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,sLattice,myStructure))
+ constitutive_nonlocal_postResults(cs+s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,myStructure))
enddo
cs = cs + ns
diff --git a/code/constitutive_phenopowerlaw.f90 b/code/constitutive_phenopowerlaw.f90
index 6b3cadc07..717196968 100644
--- a/code/constitutive_phenopowerlaw.f90
+++ b/code/constitutive_phenopowerlaw.f90
@@ -96,6 +96,9 @@ module constitutive_phenopowerlaw
constitutive_phenopowerlaw_hardeningMatrix_TwinTwin, &
constitutive_phenopowerlaw_Cslip_66
+ real(pReal), dimension(:,:), allocatable, private :: &
+ constitutive_phenopowerlaw_nonSchmidCoeff
+
public :: &
constitutive_phenopowerlaw_init, &
constitutive_phenopowerlaw_homogenizedC, &
@@ -125,6 +128,7 @@ subroutine constitutive_phenopowerlaw_init(myFile)
use lattice, only: lattice_initializeStructure, lattice_symmetrizeC66, &
lattice_maxNslipFamily, lattice_maxNtwinFamily, &
lattice_maxNinteraction, lattice_NslipSystem, lattice_NtwinSystem, &
+ lattice_maxNonSchmid, &
lattice_interactionSlipSlip, &
lattice_interactionSlipTwin, &
lattice_interactionTwinSlip, &
@@ -227,6 +231,8 @@ subroutine constitutive_phenopowerlaw_init(myFile)
constitutive_phenopowerlaw_aTolShear = 0.0_pReal
allocate(constitutive_phenopowerlaw_aTolTwinfrac(maxNinstance))
constitutive_phenopowerlaw_aTolTwinfrac = 0.0_pReal
+ allocate(constitutive_phenopowerlaw_nonSchmidCoeff(lattice_maxNonSchmid,maxNinstance))
+ constitutive_phenopowerlaw_nonSchmidCoeff = 0.0_pReal
rewind(myFile)
section = 0_pInt
@@ -338,6 +344,9 @@ subroutine constitutive_phenopowerlaw_init(myFile)
case ('interaction_twintwin')
forall (j = 1_pInt:lattice_maxNinteraction) &
constitutive_phenopowerlaw_interaction_TwinTwin(j,i) = IO_floatValue(line,positions,1_pInt+j)
+ case ('nonschmid_coefficients')
+ forall (j = 1_pInt:lattice_maxNonSchmid) &
+ constitutive_phenopowerlaw_nonSchmidCoeff(j,i) = IO_floatValue(line,positions,1_pInt+j)
case default
call IO_error(210_pInt,ext_msg=tag//' ('//constitutive_phenopowerlaw_label//')')
end select
@@ -611,9 +620,9 @@ end subroutine constitutive_phenopowerlaw_microstructure
!--------------------------------------------------------------------------------------------------
subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,state,ipc,ip,el)
use prec, only: p_vec
- use math, only: math_Plain3333to99
+ use math, only: math_Plain3333to99,math_Mandel6to33
use lattice, only: lattice_Sslip,lattice_Sslip_v,lattice_Stwin,lattice_Stwin_v, lattice_maxNslipFamily, lattice_maxNtwinFamily, &
- lattice_NslipSystem,lattice_NtwinSystem
+ lattice_NslipSystem,lattice_NtwinSystem,NnonSchmid
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase, phase_plasticityInstance
@@ -629,8 +638,9 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temp
real(pReal), dimension(3,3), intent(out) :: Lp ! plastic velocity gradient
real(pReal), dimension(3,3,3,3) :: dLp_dTstar3333 ! derivative of Lp (4th-rank tensor)
real(pReal), dimension(9,9), intent(out) :: dLp_dTstar
+ real(pReal), dimension(3,3,2) :: nonSchmid_tensor
real(pReal), dimension(constitutive_phenopowerlaw_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
- gdot_slip,dgdot_dtauslip,tau_slip
+ gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg,tau_slip_pos,tau_slip_neg
real(pReal), dimension(constitutive_phenopowerlaw_totalNtwin(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
gdot_twin,dgdot_dtautwin,tau_twin
@@ -655,20 +665,43 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temp
!--------------------------------------------------------------------------------------------------
! Calculation of Lp
- tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(1:6,index_myFamily+i,structID))
- gdot_slip(j) = constitutive_phenopowerlaw_gdot0_slip(matID)*(abs(tau_slip(j))/state(ipc,ip,el)%p(j))**&
- constitutive_phenopowerlaw_n_slip(matID)*sign(1.0_pReal,tau_slip(j))
- Lp = Lp + (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F
- gdot_slip(j)*lattice_Sslip(1:3,1:3,index_myFamily+i,structID)
+ tau_slip_pos(j) = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,structID))
+ tau_slip_neg(j) = tau_slip_pos(j)
+ nonSchmid_tensor(1:3,1:3,1) = math_Mandel6to33(lattice_Sslip_v(1:6,1,index_myFamily+i,structID))
+ nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,1)
+ do k = 1, NnonSchmid(structID)
+ tau_slip_pos(j) = tau_slip_pos(j) + constitutive_phenopowerlaw_nonSchmidCoeff(k,matID)* &
+ dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k,index_myFamily+i,structID))
+ tau_slip_neg(j) = tau_slip_neg(j) + constitutive_phenopowerlaw_nonSchmidCoeff(k,matID)* &
+ dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,structID))
+ nonSchmid_tensor(1:3,1:3,1) = nonSchmid_tensor(1:3,1:3,1) + constitutive_phenopowerlaw_nonSchmidCoeff(k,matID)*&
+ math_Mandel6to33(lattice_Sslip_v(1:6,2*k,index_myFamily+i,structID))
+ nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,2) + constitutive_phenopowerlaw_nonSchmidCoeff(k,matID)*&
+ math_Mandel6to33(lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,structID))
+ enddo
+ gdot_slip_pos(j) = 0.5_pReal*constitutive_phenopowerlaw_gdot0_slip(matID)* &
+ ((abs(tau_slip_pos(j))/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(matID))*sign(1.0_pReal,tau_slip_pos(j))
+ gdot_slip_neg(j) = 0.5_pReal*constitutive_phenopowerlaw_gdot0_slip(matID)* &
+ ((abs(tau_slip_neg(j))/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(matID))*sign(1.0_pReal,tau_slip_neg(j))
+ Lp = Lp + (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F
+ (gdot_slip_pos(j)+gdot_slip_neg(j))*lattice_Sslip(1:3,1:3,index_myFamily+i,structID)
!--------------------------------------------------------------------------------------------------
! Calculation of the tangent of Lp
- if (gdot_slip(j) /= 0.0_pReal) then
- dgdot_dtauslip(j) = gdot_slip(j)*constitutive_phenopowerlaw_n_slip(matID)/tau_slip(j)
+ if (gdot_slip_pos(j) /= 0.0_pReal) then
+ dgdot_dtauslip_pos(j) = gdot_slip_pos(j)*constitutive_phenopowerlaw_n_slip(matID)/tau_slip_pos(j)
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
- dgdot_dtauslip(j)*lattice_Sslip(k,l,index_myFamily+i,structID)* &
- lattice_Sslip(m,n,index_myFamily+i,structID)
+ dgdot_dtauslip_pos(j)*lattice_Sslip(k,l,index_myFamily+i,structID)* &
+ nonSchmid_tensor(m,n,1)
+ endif
+
+ if (gdot_slip_neg(j) /= 0.0_pReal) then
+ dgdot_dtauslip_neg(j) = gdot_slip_neg(j)*constitutive_phenopowerlaw_n_slip(matID)/tau_slip_neg(j)
+ forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
+ dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
+ dgdot_dtauslip_neg(j)*lattice_Sslip(k,l,index_myFamily+i,structID)* &
+ nonSchmid_tensor(m,n,2)
endif
enddo
enddo
@@ -711,7 +744,7 @@ end subroutine constitutive_phenopowerlaw_LpAndItsTangent
function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el)
use prec, only: p_vec
use lattice, only: lattice_Sslip_v, lattice_Stwin_v, lattice_maxNslipFamily, lattice_maxNtwinFamily, &
- lattice_NslipSystem,lattice_NtwinSystem,lattice_shearTwin
+ lattice_NslipSystem,lattice_NtwinSystem,lattice_shearTwin,NnonSchmid
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase, phase_plasticityInstance
@@ -720,12 +753,12 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el
ipc, & !< component-ID at current integration point
ip, & !< current integration point
el !< current element
- integer(pInt) matID,nSlip,nTwin,f,i,j, structID,index_Gamma,index_F,index_myFamily
+ integer(pInt) matID,nSlip,nTwin,f,i,j,k,structID,index_Gamma,index_F,index_myFamily
real(pReal) Temperature,c_SlipSlip,c_SlipTwin,c_TwinSlip,c_TwinTwin, ssat_offset
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state
real(pReal), dimension(6), intent(in) :: Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel)
real(pReal), dimension(constitutive_phenopowerlaw_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
- gdot_slip,tau_slip,left_SlipSlip,left_SlipTwin,right_SlipSlip,right_TwinSlip
+ gdot_slip,tau_slip_pos,tau_slip_neg,left_SlipSlip,left_SlipTwin,right_SlipSlip,right_TwinSlip
real(pReal), dimension(constitutive_phenopowerlaw_totalNtwin(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
gdot_twin,tau_twin,left_TwinSlip,left_TwinTwin,right_SlipTwin,right_TwinTwin
real(pReal), dimension(constitutive_phenopowerlaw_sizeDotState(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
@@ -770,9 +803,18 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el
!--------------------------------------------------------------------------------------------------
! Calculation of dot gamma
- tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(1:6,index_myFamily+i,structID))
- gdot_slip(j) = constitutive_phenopowerlaw_gdot0_slip(matID)*(abs(tau_slip(j))/state(ipc,ip,el)%p(j))**&
- constitutive_phenopowerlaw_n_slip(matID)*sign(1.0_pReal,tau_slip(j))
+ tau_slip_pos(j) = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,structID))
+ tau_slip_neg(j) = tau_slip_pos(j)
+ do k = 1, NnonSchmid(structID)
+ tau_slip_pos(j) = tau_slip_pos(j) + constitutive_phenopowerlaw_nonSchmidCoeff(k,matID)* &
+ dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k,index_myFamily+i,structID))
+ tau_slip_neg(j) = tau_slip_neg(j) + constitutive_phenopowerlaw_nonSchmidCoeff(k,matID)* &
+ dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,structID))
+ enddo
+ gdot_slip(j) = constitutive_phenopowerlaw_gdot0_slip(matID)*0.5_pReal* &
+ ((abs(tau_slip_pos(j))/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(matID) &
+ +(abs(tau_slip_neg(j))/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(matID))&
+ *sign(1.0_pReal,tau_slip_pos(j))
enddo
enddo
@@ -896,7 +938,7 @@ end function constitutive_phenopowerlaw_dotTemperature
pure function constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,state,ipc,ip,el)
use prec, only: pReal,pInt,p_vec
use lattice, only: lattice_Sslip_v,lattice_Stwin_v, lattice_maxNslipFamily, lattice_maxNtwinFamily, &
- lattice_NslipSystem,lattice_NtwinSystem
+ lattice_NslipSystem,lattice_NtwinSystem,NnonSchmid
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase,phase_plasticityInstance,phase_Noutput
@@ -910,8 +952,8 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,stat
Temperature
real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola Kirchhoff stress tensor (Mandel)
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state
- integer(pInt) matID,o,f,i,c,nSlip,nTwin,j, structID,index_Gamma,index_F,index_myFamily
- real(pReal) tau
+ integer(pInt) matID,o,f,i,c,nSlip,nTwin,j,k,structID,index_Gamma,index_F,index_myFamily
+ real(pReal) tau_slip_pos,tau_slip_neg,tau
real(pReal), dimension(constitutive_phenopowerlaw_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
constitutive_phenopowerlaw_postResults
@@ -939,10 +981,18 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,stat
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,structID)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family
j = j + 1_pInt
- tau = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,structID))
- constitutive_phenopowerlaw_postResults(c+j) = constitutive_phenopowerlaw_gdot0_slip(matID)*&
- (abs(tau)/state(ipc,ip,el)%p(j))**&
- constitutive_phenopowerlaw_n_slip(matID)*sign(1.0_pReal,tau)
+ tau_slip_pos = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,structID))
+ tau_slip_neg = tau_slip_pos
+ do k = 1, NnonSchmid(structID)
+ tau_slip_pos = tau_slip_pos + constitutive_phenopowerlaw_nonSchmidCoeff(k,matID)* &
+ dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k,index_myFamily+i,structID))
+ tau_slip_neg = tau_slip_neg + constitutive_phenopowerlaw_nonSchmidCoeff(k,matID)* &
+ dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,structID))
+ enddo
+ constitutive_phenopowerlaw_postResults(c+j) = constitutive_phenopowerlaw_gdot0_slip(matID)*0.5_pReal* &
+ ((abs(tau_slip_pos)/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(matID) &
+ +(abs(tau_slip_neg)/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(matID))&
+ *sign(1.0_pReal,tau_slip_pos)
enddo; enddo
c = c + nSlip
@@ -952,7 +1002,7 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,stat
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,structID)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family
j = j + 1_pInt
- constitutive_phenopowerlaw_postResults(c+j) = dot_product(Tstar_v,lattice_Sslip_v(1:6,index_myFamily+i,structID))
+ constitutive_phenopowerlaw_postResults(c+j) = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,structID))
enddo; enddo
c = c + nSlip
diff --git a/code/constitutive_titanmod.f90 b/code/constitutive_titanmod.f90
index 83906f4fc..71791b6fb 100644
--- a/code/constitutive_titanmod.f90
+++ b/code/constitutive_titanmod.f90
@@ -1343,7 +1343,7 @@ do f = 1_pInt,lattice_maxNslipFamily ! loop over
!* Calculation of Lp
!* Resolved shear stress on slip system
- tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure))
+ tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,1,index_myFamily+i,myStructure))
!*************************************************
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! if(myStructure>=3.and.j>3) then ! for all non-basal slip systems
diff --git a/code/lattice.f90 b/code/lattice.f90
index 6d9bb1bf9..5d32d9e28 100644
--- a/code/lattice.f90
+++ b/code/lattice.f90
@@ -43,7 +43,8 @@ module lattice
lattice_maxNtwinFamily = 4_pInt, & !< max # of twin system families over lattice structures
lattice_maxNslip = 54_pInt, & !< max # of slip systems over lattice structures
lattice_maxNtwin = 24_pInt, & !< max # of twin systems over lattice structures
- lattice_maxNinteraction = 30_pInt !< max # of interaction types (in hardening matrix part)
+ lattice_maxNinteraction = 30_pInt, & !< max # of interaction types (in hardening matrix part)
+ lattice_maxNonSchmid = 6_pInt !< max # of non schmid contributions over lattice structures
integer(pInt), allocatable, dimension(:,:), protected, public :: &
lattice_NslipSystem, & !< # of slip systems in each family
@@ -57,10 +58,10 @@ module lattice
real(pReal), allocatable, dimension(:,:,:,:), protected, public :: &
+ lattice_Sslip_v, &
lattice_Sslip !< Schmid matrices, normal, shear direction and d x n of slip systems
real(pReal), allocatable, dimension(:,:,:), protected, public :: &
- lattice_Sslip_v, &
lattice_sn, &
lattice_sd, &
lattice_st
@@ -88,6 +89,9 @@ module lattice
interactionSlipTwin, &
interactionTwinSlip, &
interactionTwinTwin
+
+ integer(pInt), allocatable, dimension(:), protected, public :: &
+ NnonSchmid !< # of Non Schmid contributions for each structure
!--------------------------------------------------------------------------------------------------
! fcc (1)
@@ -214,6 +218,13 @@ module lattice
2,2,2,2,2,2,2,2,2,1,1,1, &
2,2,2,2,2,2,2,2,2,1,1,1 &
],pInt),[lattice_fcc_Ntwin,lattice_fcc_Ntwin],order=[2,1])
+
+! Number of Non Schmid contributions for FCC
+ integer(pInt), parameter, private :: NnonSchmid_fcc = 0_pInt
+
+! Tensor for each non schmid contribution
+ real(pReal), dimension(3,3,2,NnonSchmid_fcc,lattice_fcc_Nslip), parameter, private :: &
+ lattice_nonSchmid_fcc = 0.0_pReal ! reshape([],[3,3,2,NnonSchmid_fcc,lattice_fcc_Nslip])
!--------------------------------------------------------------------------------------------------
@@ -416,6 +427,13 @@ module lattice
!< 1 --- self interaction
!< 2 --- collinear interaction
!< 3 --- other interaction
+
+! Number of Non Schmid contributions for BCC
+ integer(pInt), parameter, private :: NnonSchmid_bcc = 0_pInt
+
+! Tensor for each non schmid contribution
+ real(pReal), dimension(3,3,2,NnonSchmid_bcc,lattice_bcc_Nslip), parameter, private :: &
+ lattice_nonSchmid_bcc = 0.0_pReal ! reshape([],[3,3,2,NnonSchmid_bcc,lattice_bcc_Nslip])
!--------------------------------------------------------------------------------------------------
! hex (3+)
@@ -677,6 +695,13 @@ module lattice
20,20,20,20,20,20, 19,19,19,19,19,19, 17,17,17,17,17,17, 8, 8, 8, 8, 4, 8, &
20,20,20,20,20,20, 19,19,19,19,19,19, 17,17,17,17,17,17, 8, 8, 8, 8, 8, 4 &
],pInt),[lattice_hex_Ntwin,lattice_hex_Ntwin],order=[2,1])
+
+! Number of Non Schmid contributions for hex
+ integer(pInt), parameter, private :: NnonSchmid_hex = 0_pInt
+
+! Tensor for each non schmid contribution
+ real(pReal), dimension(3,3,2,NnonSchmid_hex,lattice_hex_Nslip), parameter, private :: &
+ lattice_nonSchmid_hex = 0.0_pReal ! reshape([],[3,3,2,NnonSchmid_hex,lattice_hex_Nslip])
public :: &
lattice_init, &
@@ -819,8 +844,9 @@ subroutine lattice_init
!$OMP END CRITICAL (write2out)
endif
+ allocate(NnonSchmid(lattice_Nstructure)); NnonSchmid = 0_pInt
allocate(lattice_Sslip(3,3,lattice_maxNslip,lattice_Nstructure)); lattice_Sslip = 0.0_pReal
- allocate(lattice_Sslip_v(6,lattice_maxNslip,lattice_Nstructure)); lattice_Sslip_v = 0.0_pReal
+ allocate(lattice_Sslip_v(6,1+2*lattice_maxNonSchmid,lattice_maxNslip,lattice_Nstructure)); lattice_Sslip_v = 0.0_pReal
allocate(lattice_sd(3,lattice_maxNslip,lattice_Nstructure)); lattice_sd = 0.0_pReal
allocate(lattice_st(3,lattice_maxNslip,lattice_Nstructure)); lattice_st = 0.0_pReal
allocate(lattice_sn(3,lattice_maxNslip,lattice_Nstructure)); lattice_sn = 0.0_pReal
@@ -869,12 +895,13 @@ integer(pInt) function lattice_initializeStructure(struct,CoverA)
real(pReal) CoverA
real(pReal), dimension(3,lattice_maxNslip) :: sd = 0.0_pReal, &
sn = 0.0_pReal
+ real(pReal), dimension(12,lattice_maxNonSchmid,lattice_maxNslip) :: sns = 0.0_pReal
real(pReal), dimension(3,lattice_maxNtwin) :: td = 0.0_pReal, &
tn = 0.0_pReal
real(pReal), dimension(lattice_maxNtwin) :: ts = 0.0_pReal
integer(pInt), dimension(lattice_maxNslipFamily) :: myNslipSystem = 0_pInt
integer(pInt), dimension(lattice_maxNtwinFamily) :: myNtwinSystem = 0_pInt
- integer(pInt) :: i,myNslip,myNtwin,myStructure = 0_pInt
+ integer(pInt) :: i,j,myNslip,myNtwin,myStructure = 0_pInt
logical :: processMe
processMe = .false.
@@ -889,9 +916,14 @@ integer(pInt) function lattice_initializeStructure(struct,CoverA)
lattice_fcc_Nstructure = lattice_fcc_Nstructure + 1_pInt ! count fcc instances
if (lattice_fcc_Nstructure == 1_pInt) then ! me is first fcc structure
processMe = .true.
+ NnonSchmid(myStructure) = NnonSchmid_fcc ! Currently no known non schmid contributions for FCC (to be changed later)
do i = 1_pInt,myNslip ! assign slip system vectors
sd(1:3,i) = lattice_fcc_systemSlip(1:3,i)
sn(1:3,i) = lattice_fcc_systemSlip(4:6,i)
+ do j = 1_pInt, NnonSchmid_fcc
+ sns(1:6,j,i) = math_Mandel33to6(lattice_nonSchmid_fcc(1:3,1:3,1,j,i))
+ sns(7:12,j,i) = math_Mandel33to6(lattice_nonSchmid_fcc(1:3,1:3,2,j,i))
+ enddo
enddo
do i = 1_pInt,myNtwin ! assign twin system vectors and shears
td(1:3,i) = lattice_fcc_systemTwin(1:3,i)
@@ -913,9 +945,14 @@ integer(pInt) function lattice_initializeStructure(struct,CoverA)
lattice_bcc_Nstructure = lattice_bcc_Nstructure + 1_pInt ! count bcc instances
if (lattice_bcc_Nstructure == 1_pInt) then ! me is first bcc structure
processMe = .true.
+ NnonSchmid(myStructure) = NnonSchmid_BCC ! 5 known non schmid contributions for BCC (A. Koester, A. Ma, A. Hartmaier 2012)
do i = 1_pInt,myNslip ! assign slip system vectors
sd(1:3,i) = lattice_bcc_systemSlip(1:3,i)
sn(1:3,i) = lattice_bcc_systemSlip(4:6,i)
+ do j = 1_pInt, NnonSchmid_bcc
+ sns(1:6,j,i) = math_Mandel33to6(lattice_nonSchmid_bcc(1:3,1:3,1,j,i))
+ sns(7:12,j,i) = math_Mandel33to6(lattice_nonSchmid_bcc(1:3,1:3,2,j,i))
+ enddo
enddo
do i = 1_pInt,myNtwin ! assign twin system vectors and shears
td(1:3,i) = lattice_bcc_systemTwin(1:3,i)
@@ -937,6 +974,7 @@ integer(pInt) function lattice_initializeStructure(struct,CoverA)
myNslip = lattice_hex_Nslip ! overall number of slip systems
myNtwin = lattice_hex_Ntwin ! overall number of twin systems
processMe = .true.
+ NnonSchmid(myStructure) = NnonSchmid_hex ! Currently no known non schmid contributions for hex (to be changed later)
! converting from 4 axes coordinate system (a1=a2=a3=c) to ortho-hexgonal system (a, b, c)
do i = 1_pInt,myNslip
sd(1,i) = lattice_hex_systemSlip(1,i)*1.5_pReal ! direction [uvtw]->[3u/2 (u+2v)*sqrt(3)/2 w*(c/a)]
@@ -945,6 +983,10 @@ integer(pInt) function lattice_initializeStructure(struct,CoverA)
sn(1,i) = lattice_hex_systemSlip(5,i) ! plane (hkil)->(h (h+2k)/sqrt(3) l/(c/a))
sn(2,i) = (lattice_hex_systemSlip(5,i)+2.0_pReal*lattice_hex_systemSlip(6,i))/sqrt(3.0_pReal)
sn(3,i) = lattice_hex_systemSlip(8,i)/CoverA
+ do j = 1_pInt, NnonSchmid_hex
+ sns(1:6,j,i) = math_Mandel33to6(lattice_nonSchmid_hex(1:3,1:3,1,j,i))
+ sns(7:12,j,i) = math_Mandel33to6(lattice_nonSchmid_hex(1:3,1:3,2,j,i))
+ enddo
enddo
do i = 1_pInt,myNtwin
td(1,i) = lattice_hex_systemTwin(1,i)*1.5_pReal
@@ -983,7 +1025,11 @@ integer(pInt) function lattice_initializeStructure(struct,CoverA)
lattice_sn(1:3,i,myStructure))
lattice_Sslip(1:3,1:3,i,myStructure) = math_tensorproduct(lattice_sd(1:3,i,myStructure), &
lattice_sn(1:3,i,myStructure))
- lattice_Sslip_v(1:6,i,myStructure) = math_Mandel33to6(math_symmetric33(lattice_Sslip(1:3,1:3,i,myStructure)))
+ lattice_Sslip_v(1:6,1,i,myStructure) = math_Mandel33to6(math_symmetric33(lattice_Sslip(1:3,1:3,i,myStructure)))
+ do j = 1_pInt, NnonSchmid(myStructure)
+ lattice_Sslip_v(1:6,2*j,i,myStructure) = sns(1:6,j,i)
+ lattice_Sslip_v(1:6,2*j+1,i,myStructure) = sns(7:12,j,i)
+ enddo
if (abs(math_trace33(lattice_Sslip(1:3,1:3,i,myStructure))) > 1.0e-8_pReal) &
call IO_error(0_pInt,myStructure,i,0_pInt,ext_msg = 'dilatational slip Schmid matrix')
enddo