From e6d88e44588a57f74760cccaf2a92062ffa0ccb4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 22 Jan 2014 15:47:49 +0000 Subject: [PATCH] prevent twin volume fraction from going above 1.0 --- code/constitutive_phenopowerlaw.f90 | 57 ++++++++++++++--------------- 1 file changed, 27 insertions(+), 30 deletions(-) diff --git a/code/constitutive_phenopowerlaw.f90 b/code/constitutive_phenopowerlaw.f90 index 96fdb1d3f..b156c6992 100644 --- a/code/constitutive_phenopowerlaw.f90 +++ b/code/constitutive_phenopowerlaw.f90 @@ -720,8 +720,7 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar p_vec use math, only: & math_Plain3333to99, & - math_Mandel6to33, & - math_mul33xx33 + math_Mandel6to33 use lattice, only: & lattice_Sslip, & lattice_Sslip_v, & @@ -790,30 +789,28 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar !-------------------------------------------------------------------------------------------------- ! Calculation of Lp + 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) = lattice_Sslip(1:3,1:3,1,index_myFamily+i,structID) nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,1) do k = 1,lattice_NnonSchmid(structID) - nonSchmid_tensor(1:3,1:3,1) = nonSchmid_tensor(1:3,1:3,1) & - + constitutive_phenopowerlaw_nonSchmidCoeff(k,matID) & - * lattice_Sslip(1:3,1:3,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) & - * lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+i,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)*& + lattice_Sslip(1:3,1:3,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)*& + lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+i,structID) enddo - tau_slip_pos(j) = math_mul33xx33(math_Mandel6to33(Tstar_v), nonSchmid_tensor(1:3,1:3,1)) - tau_slip_neg(j) = math_mul33xx33(math_Mandel6to33(Tstar_v), nonSchmid_tensor(1:3,1:3,2)) - gdot_slip_pos(j) = 0.5_pReal * (1.0_pReal - state(ipc,ip,el)%p(index_F)) & ! 1-F - * 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 * (1.0_pReal - state(ipc,ip,el)%p(index_F)) & ! 1-F - * 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 + (gdot_slip_pos(j) + gdot_slip_neg(j)) & - * lattice_Sslip(1:3,1:3,1,index_myFamily+i,structID) + 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,1,index_myFamily+i,structID) !-------------------------------------------------------------------------------------------------- ! Calculation of the tangent of Lp @@ -844,12 +841,11 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar !-------------------------------------------------------------------------------------------------- ! Calculation of Lp tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,structID)) - gdot_twin(j) = (1.0_pReal - state(ipc,ip,el)%p(index_F)) & ! 1-F - * constitutive_phenopowerlaw_gdot0_twin(matID) & - * (abs(tau_twin(j)) / state(ipc,ip,el)%p(nSlip+j)) & - **constitutive_phenopowerlaw_n_twin(matID) & - * max(0.0_pReal,sign(1.0_pReal,tau_twin(j))) - Lp = Lp + gdot_twin(j) * lattice_Stwin(1:3,1:3,index_myFamily+i,structID) + gdot_twin(j) = (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F + constitutive_phenopowerlaw_gdot0_twin(matID)*& + (abs(tau_twin(j))/state(ipc,ip,el)%p(nSlip+j))**& + constitutive_phenopowerlaw_n_twin(matID)*max(0.0_pReal,sign(1.0_pReal,tau_twin(j))) + Lp = Lp + gdot_twin(j)*lattice_Stwin(1:3,1:3,index_myFamily+i,structID) !-------------------------------------------------------------------------------------------------- ! Calculation of the tangent of Lp @@ -1029,8 +1025,9 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,state,ipc,ip,el) c_TwinTwin * left_TwinTwin(j) * & dot_product(constitutive_phenopowerlaw_hardeningMatrix_TwinTwin(j,1:nTwin,matID), & right_TwinTwin*gdot_twin) ! dot gamma_twin modulated by right-side twin factor - constitutive_phenopowerlaw_dotState(index_F) = constitutive_phenopowerlaw_dotState(index_F) + & - gdot_twin(j)/lattice_shearTwin(index_myFamily+i,structID) + if (state(ipc,ip,el)%p(index_F) < 0.98_pReal) & ! ensure twin volume fractions stays below 1.0 + constitutive_phenopowerlaw_dotState(index_F) = constitutive_phenopowerlaw_dotState(index_F) + & + gdot_twin(j)/lattice_shearTwin(index_myFamily+i,structID) constitutive_phenopowerlaw_dotState(offset_accshear_twin+j) = abs(gdot_twin(j)) enddo enddo twinFamiliesLoop2