From f9cec3079e2a83a038f5e83b1d21223718d94359 Mon Sep 17 00:00:00 2001 From: Pratheek Shanthraj Date: Fri, 31 Oct 2014 19:03:08 +0000 Subject: [PATCH] overhaul of kinematics from FeFp to FeFiFp decomposition --- code/constitutive.f90 | 263 ++++++++++++++++++++++-------- code/crystallite.f90 | 307 +++++++++++++++++++++++++++-------- code/damage_anisoBrittle.f90 | 303 +++++++++++++++++++++++++++------- code/thermal_adiabatic.f90 | 158 ++++++++++++++++++ 4 files changed, 838 insertions(+), 193 deletions(-) diff --git a/code/constitutive.f90 b/code/constitutive.f90 index 1112ce77f..5a9ab361a 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -27,6 +27,10 @@ module constitutive constitutive_damagedC, & constitutive_microstructure, & constitutive_LpAndItsTangent, & + constitutive_LiAndItsTangent, & + constitutive_getFi, & + constitutive_getFi0, & + constitutive_getPartionedFi0, & constitutive_TandItsTangent, & constitutive_collectDotState, & constitutive_collectDeltaState, & @@ -34,12 +38,10 @@ module constitutive constitutive_putLocalDamage, & constitutive_getDamage, & constitutive_getSlipDamage, & - constitutive_getDamageStrain, & constitutive_getDamageDiffusion33, & constitutive_getAdiabaticTemperature, & constitutive_putAdiabaticTemperature, & constitutive_getTemperature, & - constitutive_getThermalStrain, & constitutive_getLocalVacancyConcentration, & constitutive_putLocalVacancyConcentration, & constitutive_getVacancyConcentration, & @@ -701,6 +703,199 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, ipc, ip, el) end subroutine constitutive_LpAndItsTangent +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the velocity gradient +!-------------------------------------------------------------------------------------------------- +subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar, Tstar_v, Lp, ipc, ip, el) + use prec, only: & + pReal + use material, only: & + phase_damage, & + phase_thermal, & + material_phase, & + LOCAL_DAMAGE_anisoBrittle_ID, & + LOCAL_THERMAL_adiabatic_ID + use damage_anisoBrittle, only: & + damage_anisoBrittle_LdAndItsTangent + use thermal_adiabatic, only: & + thermal_adiabatic_LTAndItsTangent + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal), intent(in), dimension(6) :: & + Tstar_v !< 2nd Piola-Kirchhoff stress + real(pReal), intent(in), dimension(3,3) :: & + Lp !< plastic velocity gradient + real(pReal), intent(out), dimension(3,3) :: & + Li !< intermediate velocity gradient + real(pReal), intent(out), dimension(9,9) :: & + dLi_dTstar !< derivative of Li with respect to Tstar (2nd-order tensor) + real(pReal), dimension(3,3) :: & + Li_temp !< intermediate velocity gradient + real(pReal), dimension(9,9) :: & + dLi_dTstar_temp !< derivative of Li with respect to Tstar (4th-order tensor) + + Li = 0.0_pReal + dLi_dTstar = 0.0_pReal + + select case (phase_damage(material_phase(ipc,ip,el))) + case (LOCAL_DAMAGE_anisoBrittle_ID) + call damage_anisoBrittle_LdAndItsTangent(Li_temp, dLi_dTstar_temp, Tstar_v, ipc, ip, el) + Li = Li + Li_temp + dLi_dTstar = dLi_dTstar + dLi_dTstar_temp + + end select + + select case (phase_thermal(material_phase(ipc,ip,el))) + case (LOCAL_THERMAL_adiabatic_ID) + call thermal_adiabatic_LTAndItsTangent(Li_temp, dLi_dTstar_temp, Tstar_v, Lp, ipc, ip, el) + Li = Li + Li_temp + dLi_dTstar = dLi_dTstar + dLi_dTstar_temp + + end select + +end subroutine constitutive_LiAndItsTangent + + +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the intermediate deformation gradient +!-------------------------------------------------------------------------------------------------- +pure function constitutive_getFi(ipc, ip, el) + use prec, only: & + pReal + use math, only: & + math_I3, & + math_mul33x33 + use material, only: & + phase_damage, & + phase_thermal, & + material_phase, & + LOCAL_DAMAGE_anisoBrittle_ID, & + LOCAL_THERMAL_adiabatic_ID + use damage_anisoBrittle, only: & + damage_anisoBrittle_getFd + use thermal_adiabatic, only: & + thermal_adiabatic_getFT + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal), dimension(3,3) :: & + constitutive_getFi !< intermediate deformation gradient + + constitutive_getFi = math_I3 + + select case (phase_damage(material_phase(ipc,ip,el))) + case (LOCAL_DAMAGE_anisoBrittle_ID) + constitutive_getFi = math_mul33x33(constitutive_getFi,damage_anisoBrittle_getFd (ipc, ip, el)) + + end select + + select case (phase_thermal(material_phase(ipc,ip,el))) + case (LOCAL_THERMAL_adiabatic_ID) + constitutive_getFi = math_mul33x33(constitutive_getFi,thermal_adiabatic_getFT (ipc, ip, el)) + + end select + +end function constitutive_getFi + + +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the intermediate deformation gradient +!-------------------------------------------------------------------------------------------------- +pure function constitutive_getFi0(ipc, ip, el) + use prec, only: & + pReal + use math, only: & + math_I3, & + math_mul33x33 + use material, only: & + phase_damage, & + phase_thermal, & + material_phase, & + LOCAL_DAMAGE_anisoBrittle_ID, & + LOCAL_THERMAL_adiabatic_ID + use damage_anisoBrittle, only: & + damage_anisoBrittle_getFd0 + use thermal_adiabatic, only: & + thermal_adiabatic_getFT0 + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal), dimension(3,3) :: & + constitutive_getFi0 !< intermediate deformation gradient + + constitutive_getFi0 = math_I3 + + select case (phase_damage(material_phase(ipc,ip,el))) + case (LOCAL_DAMAGE_anisoBrittle_ID) + constitutive_getFi0 = math_mul33x33(constitutive_getFi0,damage_anisoBrittle_getFd0 (ipc, ip, el)) + + end select + + select case (phase_thermal(material_phase(ipc,ip,el))) + case (LOCAL_THERMAL_adiabatic_ID) + constitutive_getFi0 = math_mul33x33(constitutive_getFi0,thermal_adiabatic_getFT0 (ipc, ip, el)) + + end select + +end function constitutive_getFi0 + + +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the intermediate deformation gradient +!-------------------------------------------------------------------------------------------------- +pure function constitutive_getPartionedFi0(ipc, ip, el) + use prec, only: & + pReal + use math, only: & + math_I3, & + math_mul33x33 + use material, only: & + phase_damage, & + phase_thermal, & + material_phase, & + LOCAL_DAMAGE_anisoBrittle_ID, & + LOCAL_THERMAL_adiabatic_ID + use damage_anisoBrittle, only: & + damage_anisoBrittle_getPartionedFd0 + use thermal_adiabatic, only: & + thermal_adiabatic_getPartionedFT0 + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal), dimension(3,3) :: & + constitutive_getPartionedFi0 !< intermediate deformation gradient + + constitutive_getPartionedFi0 = math_I3 + + select case (phase_damage(material_phase(ipc,ip,el))) + case (LOCAL_DAMAGE_anisoBrittle_ID) + constitutive_getPartionedFi0 = math_mul33x33(constitutive_getPartionedFi0, & + damage_anisoBrittle_getPartionedFd0(ipc, ip, el)) + + end select + + select case (phase_thermal(material_phase(ipc,ip,el))) + case (LOCAL_THERMAL_adiabatic_ID) + constitutive_getPartionedFi0 = math_mul33x33(constitutive_getPartionedFi0, & + thermal_adiabatic_getPartionedFT0(ipc, ip, el)) + + end select + +end function constitutive_getPartionedFi0 + !-------------------------------------------------------------------------------------------------- !> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to @@ -1160,41 +1355,6 @@ function constitutive_getSlipDamage(nSlip, Tstar_v, ipc, ip, el) end function constitutive_getSlipDamage -!-------------------------------------------------------------------------------------------------- -!> @brief returns damage deformation gradient -!-------------------------------------------------------------------------------------------------- -function constitutive_getDamageStrain(ipc, ip, el) - use prec, only: & - pReal - use math, only: & - math_I3 - use material, only: & - material_phase, & - LOCAL_DAMAGE_anisoBrittle_ID, & - phase_damage - use damage_anisoBrittle, only: & - damage_anisoBrittle_getDamageStrain - - - implicit none - integer(pInt), intent(in) :: & - ipc, & !< grain number - ip, & !< integration point number - el !< element number - real(pReal), dimension(3,3) :: & - constitutive_getDamageStrain - - - select case (phase_damage(material_phase(ipc,ip,el))) - case (LOCAL_DAMAGE_anisoBrittle_ID) - constitutive_getDamageStrain = damage_anisoBrittle_getDamageStrain(ipc, ip, el) - case default - constitutive_getDamageStrain = math_I3 - - end select - -end function constitutive_getDamageStrain - !-------------------------------------------------------------------------------------------------- !> @brief returns damage diffusion tensor !-------------------------------------------------------------------------------------------------- @@ -1328,35 +1488,6 @@ function constitutive_getTemperature(ipc, ip, el) end function constitutive_getTemperature -!-------------------------------------------------------------------------------------------------- -!> @brief returns thermal deformation gradient -!-------------------------------------------------------------------------------------------------- -function constitutive_getThermalStrain(ipc, ip, el) - use prec, only: & - pReal - use math, only: & - math_I3 - use lattice, only: & - lattice_referenceTemperature, & - lattice_thermalExpansion33 - use material, only: & - mappingConstitutive - - implicit none - integer(pInt), intent(in) :: & - ipc, & !< grain number - ip, & !< integration point number - el !< element number - real(pReal), dimension(3,3) :: & - constitutive_getThermalStrain - - constitutive_getThermalStrain = math_I3 + & - (constitutive_getTemperature(ipc, ip, el) - & - lattice_referenceTemperature(mappingConstitutive(2,ipc,ip,el)))* & - lattice_thermalExpansion33(1:3,1:3,mappingConstitutive(2,ipc,ip,el)) - -end function constitutive_getThermalStrain - !-------------------------------------------------------------------------------------------------- !> @brief returns local vacancy concentration !-------------------------------------------------------------------------------------------------- diff --git a/code/crystallite.f90 b/code/crystallite.f90 index edb946936..155b1fb3e 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -516,7 +516,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco) math_I3, & math_mul3333xx3333, & math_mul33xx33, & - math_invert + math_invert, & + math_det33 use FEsolving, only: & FEsolving_execElem, & FEsolving_execIP @@ -537,7 +538,11 @@ subroutine crystallite_stressAndItsTangent(updateJaco) homogenization_maxNgrains use constitutive, only: & constitutive_TandItsTangent, & - constitutive_LpAndItsTangent + constitutive_LpAndItsTangent, & + constitutive_LiAndItsTangent, & + constitutive_getFi, & + constitutive_getFi0, & + constitutive_getPartionedFi0 implicit none logical, intent(in) :: & @@ -579,12 +584,21 @@ subroutine crystallite_stressAndItsTangent(updateJaco) logical, dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & convergenceFlag_backup ! local variables used for calculating analytic Jacobian - real(pReal), dimension(3,3):: junk + real(pReal) :: detInvFi + real(pReal), dimension(3,3) :: junk, & + Fi, & + invFi, & + Fi0, & + invFi0 real(pReal), dimension(3,3,3,3) :: dSdFe, & - dFedF, & dSdF, & junk2, & - dLpdS,dFpinvdF,rhs_3333,lhs_3333,temp_3333 + dLidS, & + dLpdS, & + dFpinvdF, & + rhs_3333, & + lhs_3333, & + temp_3333 real(pReal), dimension(9,9):: temp_99 logical :: error @@ -624,8 +638,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco) crystallite_dPdF0(1:3,1:3,1:3,1:3,g,i,e) = crystallite_partioneddPdF0(1:3,1:3,1:3,1:3,g,i,e) ! ...stiffness crystallite_subF0(1:3,1:3,g,i,e) = crystallite_partionedF0(1:3,1:3,g,i,e) ! ...def grad crystallite_subTstar0_v(1:6,g,i,e) = crystallite_partionedTstar0_v(1:6,g,i,e) !...2nd PK stress - crystallite_subFe0(1:3,1:3,g,i,e) = math_mul33x33(crystallite_subF0(1:3,1:3,g,i,e), & - math_inv33(crystallite_subFp0(1:3,1:3,g,i,e))) ! only needed later on for stiffness calculation + crystallite_subFe0(1:3,1:3,g,i,e) = math_mul33x33(math_mul33x33(crystallite_subF0(1:3,1:3,g,i,e), & + math_inv33(crystallite_subFp0(1:3,1:3,g,i,e))), & + math_inv33(constitutive_getFi0(g,i,e))) ! only needed later on for stiffness calculation crystallite_subFrac(g,i,e) = 0.0_pReal crystallite_subStep(g,i,e) = 1.0_pReal/subStepSizeCryst crystallite_todo(g,i,e) = .true. @@ -893,8 +908,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco) !$OMP FLUSH(crystallite_subF0) crystallite_subLp0(1:3,1:3,g,i,e) = crystallite_Lp(1:3,1:3,g,i,e) ! ...plastic velocity gradient crystallite_subFp0(1:3,1:3,g,i,e) = crystallite_Fp(1:3,1:3,g,i,e) ! ...plastic def grad - crystallite_subFe0(1:3,1:3,g,i,e) = math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), & - crystallite_invFp(1:3,1:3,g,i,e)) ! only needed later on for stiffness calculation + crystallite_subFe0(1:3,1:3,g,i,e) = math_mul33x33(math_mul33x33(crystallite_subF (1:3,1:3,g,i,e), & + crystallite_invFp(1:3,1:3,g,i,e)), & + math_inv33(constitutive_getFi(g,i,e))) ! only needed later on for stiffness calculation !if abbrevation, make c and p private in omp plasticState(mappingConstitutive(2,g,i,e))%subState0(:,mappingConstitutive(1,g,i,e)) = & plasticState(mappingConstitutive(2,g,i,e))%state( :,mappingConstitutive(1,g,i,e)) @@ -979,6 +995,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco) - crystallite_partionedF0(1:3,1:3,g,i,e)) !$OMP FLUSH(crystallite_subF) crystallite_Fe(1:3,1:3,g,i,e) = math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), crystallite_invFp(1:3,1:3,g,i,e)) + crystallite_Fe(1:3,1:3,g,i,e) = math_mul33x33(math_mul33x33(crystallite_subF (1:3,1:3,g,i,e), & + crystallite_invFp(1:3,1:3,g,i,e)), & + math_inv33(constitutive_getFi(g,i,e))) crystallite_subdt(g,i,e) = crystallite_subStep(g,i,e) * crystallite_dt(g,i,e) crystallite_converged(g,i,e) = .false. ! start out non-converged endif @@ -1054,9 +1073,12 @@ subroutine crystallite_stressAndItsTangent(updateJaco) write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,/)') '<< CRYST >> no convergence: respond fully elastic at el (elFE) ip g ', & e,'(',mesh_element(1,e),')',i,g invFp = math_inv33(crystallite_partionedFp0(1:3,1:3,g,i,e)) - Fe_guess = math_mul33x33(crystallite_partionedF(1:3,1:3,g,i,e), invFp) + Fe_guess = math_mul33x33(math_mul33x33(crystallite_partionedF(1:3,1:3,g,i,e), invFp), & + math_inv33(constitutive_getPartionedFi0(g,i,e))) call constitutive_TandItsTangent(Tstar, junk2, Fe_guess,g,i,e) - crystallite_P(1:3,1:3,g,i,e) = math_mul33x33(Fe_guess,math_mul33x33(Tstar,transpose(invFp))) + crystallite_P(1:3,1:3,g,i,e) = math_mul33x33(math_mul33x33(Fe_guess,constitutive_getPartionedFi0(g,i,e)), & + math_mul33x33(Tstar,transpose(invFp)))* & + math_det33(constitutive_getPartionedFi0(g,i,e)) endif if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt & .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) & @@ -1082,7 +1104,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco) ! --- ANALYTIC JACOBIAN --- - !$OMP PARALLEL DO PRIVATE(dFedF,dSdF,dSdFe,dLpdS,dFpinvdF,rhs_3333,lhs_3333,temp_3333,temp_99,junk,myNgrains) + !$OMP PARALLEL DO PRIVATE(dFedF,dSdF,dSdFe,dLpdS,dFpinvdF,rhs_3333,lhs_3333,temp_99,junk,dLidS,Fi,invFi,Fi0,invFi0,detInvFi,temp_3333,myNgrains) elementLooping6: do e = FEsolving_execElem(1),FEsolving_execElem(2) myNgrains = homogenization_Ngrains(mesh_element(3,e)) do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed @@ -1091,46 +1113,71 @@ subroutine crystallite_stressAndItsTangent(updateJaco) call constitutive_LpAndItsTangent(junk,temp_99,crystallite_Tstar_v(1:6,g,i,e), & g,i,e) dLpdS = reshape(temp_99,shape=[3,3,3,3]) + call constitutive_LiAndItsTangent(junk,temp_99,crystallite_Tstar_v(1:6,g,i,e), & + crystallite_Lp(1:3,1:3,g,i,e), g,i,e) + dLidS = reshape(temp_99,shape=[3,3,3,3]) + Fi = constitutive_getFi(g,i,e) + invFi = math_inv33(Fi) + Fi0 = constitutive_getFi0(g,i,e) + invFi0 = math_inv33(Fi0) + detInvFi = math_det33(invFi) + + junk = math_transpose33(math_mul33x33(crystallite_invFp(1:3,1:3,g,i,e),invFi)) rhs_3333 = 0.0_pReal forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & - rhs_3333(p,o,1:3,1:3) = math_mul33x33(dSdFe(p,o,1:3,1:3), & - math_transpose33(crystallite_invFp(1:3,1:3,g,i,e))) - junk = math_mul33x33(math_mul33x33(crystallite_Fe(1:3,1:3,g,i,e), & - crystallite_Fp(1:3,1:3,g,i,e)), & - math_inv33(crystallite_Fp0(1:3,1:3,g,i,e))) + rhs_3333(p,o,1:3,1:3) = math_mul33x33(dSdFe(p,o,1:3,1:3),junk) + temp_3333 = 0.0_pReal + junk = math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), & + math_inv33(crystallite_subFp0(1:3,1:3,g,i,e))) forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & - temp_3333(1:3,1:3,p,o) = math_mul33x33(junk,dLpdS(1:3,1:3,p,o)) - lhs_3333 = crystallite_dt(g,i,e)*math_mul3333xx3333(dSdFe,temp_3333) + temp_3333(1:3,1:3,p,o) = math_mul33x33(math_mul33x33(junk,dLpdS(1:3,1:3,p,o)),invFi) + + junk = math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), & + crystallite_invFp(1:3,1:3,g,i,e)), invFi0) + forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & + temp_3333(1:3,1:3,p,o) = temp_3333(1:3,1:3,p,o) + math_mul33x33(junk,dLidS(1:3,1:3,p,o)) + + lhs_3333 = crystallite_subdt(g,i,e)*math_mul3333xx3333(dSdFe,temp_3333) + call math_invert(9_pInt,math_identity2nd(9_pInt)+reshape(lhs_3333,shape=[9,9]),temp_99,error) if (error) call IO_error(error_ID=400_pInt,ext_msg='analytic tangent inversion') dSdF = math_mul3333xx3333(reshape(temp_99,shape=[3,3,3,3]),rhs_3333) - temp_3333 = math_mul3333xx3333(dLpdS,dSdF) + dFpinvdF = 0.0_pReal + temp_3333 = math_mul3333xx3333(dLpdS,dSdF) forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & - dFpinvdF(1:3,1:3,p,o) = -crystallite_dt(g,i,e)* & - math_mul33x33(math_inv33(crystallite_Fp0(1:3,1:3,g,i,e)), & - temp_3333(1:3,1:3,p,o)) - temp_3333 = 0.0_pReal - forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & - temp_3333(o,p,o,1:3) = crystallite_invFp(1:3,p,g,i,e) ! dFe^T_ij/dF_kl = delta_jk * (Fp current^-1)_li - dFedF = 0.0_pReal - forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & - dFedF(1:3,1:3,p,o) = temp_3333(1:3,1:3,p,o) + & - math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), & - dFpinvdF(1:3,1:3,p,o)) + dFpinvdF(1:3,1:3,p,o) = -crystallite_subdt(g,i,e)* & + math_mul33x33(math_inv33(crystallite_subFp0(1:3,1:3,g,i,e)), & + math_mul33x33(temp_3333(1:3,1:3,p,o),invFi)) + + crystallite_dPdF(1:3,1:3,1:3,1:3,g,i,e) = 0.0_pReal + junk = math_mul33x33(crystallite_invFp(1:3,1:3,g,i,e), & + math_mul33x33(math_Mandel6to33(crystallite_Tstar_v(1:6,g,i,e)), & + math_transpose33(crystallite_invFp(1:3,1:3,g,i,e))))/detInvFi + forall(p=1_pInt:3_pInt) & + crystallite_dPdF(p,1:3,p,1:3,g,i,e) = math_transpose33(junk) + junk = math_mul33x33(math_Mandel6to33(crystallite_Tstar_v(1:6,g,i,e)), & + math_transpose33(crystallite_invFp(1:3,1:3,g,i,e)))/detInvFi forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & - crystallite_dPdF(1:3,1:3,o,p,g,i,e) = & - math_mul33x33(math_mul33x33(dFedF(1:3,1:3,o,p),& - math_Mandel6to33(crystallite_Tstar_v(1:6,g,i,e))), & - math_transpose33(crystallite_invFp(1:3,1:3,g,i,e))) + & ! dP/dF = dFe/dF * S * Fp^-T... - math_mul33x33(crystallite_Fe(1:3,1:3,g,i,e), & - math_mul33x33(dSdF(1:3,1:3,o,p), & - math_transpose33(crystallite_invFp(1:3,1:3,g,i,e)))) + & ! + Fe * dS/dF * Fp^-T - math_mul33x33(crystallite_Fe(1:3,1:3,g,i,e), & - math_mul33x33(math_Mandel6to33(crystallite_Tstar_v(1:6,g,i,e)), & - math_transpose33(dFpinvdF(1:3,1:3,p,o)))) ! + Fe * S * dFp^-T/dF + crystallite_dPdF(1:3,1:3,p,o,g,i,e) = crystallite_dPdF(1:3,1:3,p,o,g,i,e) + & + math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,g,i,e),dFpinvdF(1:3,1:3,p,o)),junk) + + junk = math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), & + crystallite_invFp(1:3,1:3,g,i,e))/detInvFi + forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & + crystallite_dPdF(1:3,1:3,p,o,g,i,e) = crystallite_dPdF(1:3,1:3,p,o,g,i,e) + & + math_mul33x33(math_mul33x33(junk,dSdF(1:3,1:3,p,o)), & + math_transpose33(crystallite_invFp(1:3,1:3,g,i,e))) + + junk = math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,g,i,e), & + crystallite_invFp(1:3,1:3,g,i,e)), & + math_Mandel6to33(crystallite_Tstar_v(1:6,g,i,e)))/detInvFi + forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) & + crystallite_dPdF(1:3,1:3,p,o,g,i,e) = crystallite_dPdF(1:3,1:3,p,o,g,i,e) + & + math_mul33x33(junk,math_transpose33(dFpinvdF(1:3,1:3,p,o))) + enddo; enddo enddo elementLooping6 !$OMP END PARALLEL DO @@ -3296,11 +3343,12 @@ logical function crystallite_integrateStress(& debug_cumLpTicks, & debug_StressLoopDistribution use constitutive, only: constitutive_LpAndItsTangent, & - constitutive_TandItsTangent, & - constitutive_getThermalStrain, & - constitutive_getDamageStrain + constitutive_LiAndItsTangent, & + constitutive_getFi0, & + constitutive_TandItsTangent use math, only: math_mul33x33, & math_mul33xx33, & + math_mul3333xx3333, & math_mul66x6, & math_mul99x99, & math_transpose33, & @@ -3328,45 +3376,65 @@ logical function crystallite_integrateStress(& !*** local variables ***! real(pReal), dimension(3,3):: Fg_new, & ! deformation gradient at end of timestep Fp_current, & ! plastic deformation gradient at start of timestep + Fi_current, & ! intermediate deformation gradient at start of timestep Fp_new, & ! plastic deformation gradient at end of timestep Fe_new, & ! elastic deformation gradient at end of timestep invFp_new, & ! inverse of Fp_new invFp_current, & ! inverse of Fp_current + invFi_current, & ! inverse of Fp_current Lpguess, & ! current guess for plastic velocity gradient Lpguess_old, & ! known last good guess for plastic velocity gradient Lp_constitutive, & ! plastic velocity gradient resulting from constitutive law residuum, & ! current residuum of plastic velocity gradient residuum_old, & ! last residuum of plastic velocity gradient deltaLp, & ! direction of next guess - Tstar,& ! 2nd Piola-Kirchhoff Stress - A,& + Liguess, & ! current guess for plastic velocity gradient + Liguess_old, & ! known last good guess for plastic velocity gradient + Li_constitutive, & ! plastic velocity gradient resulting from constitutive law + residuumI, & ! current residuum of plastic velocity gradient + residuumI_old, & ! last residuum of plastic velocity gradient + deltaLi, & ! direction of next guess + Tstar, & ! 2nd Piola-Kirchhoff Stress in plastic (lattice) configuration + Tstar_unloaded, & ! 2nd Piola-Kirchhoff Stress in unloaded configuration + A, & B, & Fe, & ! elastic deformation gradient Fi, & ! gradient of intermediate deformation stages - Ci, & ! stretch of intermediate deformation stages - invFi + invFi, & + temp_33 real(pReal), dimension(6):: Tstar_v ! 2nd Piola-Kirchhoff Stress in Mandel-Notation real(pReal), dimension(9):: work ! needed for matrix inversion by LAPACK integer(pInt), dimension(9) :: ipiv ! needed for matrix inversion by LAPACK real(pReal), dimension(9,9) :: dLp_dT_constitutive, & ! partial derivative of plastic velocity gradient calculated by constitutive law + dLi_dT_constitutive, & ! partial derivative of intermediate velocity gradient calculated by constitutive law dT_dFe_constitutive, & ! partial derivative of 2nd Piola-Kirchhoff stress calculated by constitutive law dFe_dLp, & ! partial derivative of elastic deformation gradient dR_dLp, & ! partial derivative of residuum (Jacobian for NEwton-Raphson scheme) - dR_dLp2 ! working copy of dRdLp + dR_dLp2, & ! working copy of dRdLp + dRI_dLp, & ! partial derivative of residuumI (Jacobian for NEwton-Raphson scheme) + dRI_dLp2 ! working copy of dRIdLp real(pReal), dimension(3,3,3,3):: dT_dFe3333, & ! partial derivative of 2nd Piola-Kirchhoff stress - dFe_dLp3333 ! partial derivative of elastic deformation gradient + dFe_dLp3333, & ! partial derivative of elastic deformation gradient + dInvFi_dLp3333, & + dFe_dInvFi3333, & + dT_dInvFi3333, & + temp_3333 real(pReal) det, & ! determinant + detInvFi, & steplength0, & steplength, & + steplengthI0, & + steplengthI, & dt, & ! time increment - aTol, & - detFi + aTol logical error ! flag indicating an error integer(pInt) NiterationStress, & ! number of stress integrations + NiterationStressI, & ! number of inner stress integrations ierr, & ! error indicator for LAPACK o, & p, & - jacoCounter ! counter to check for Jacobian update + jacoCounter, & + jacoICounter ! counters to check for Jacobian update integer(pLongInt) tick, & tock, & tickrate, & @@ -3424,20 +3492,28 @@ logical function crystallite_integrateStress(& return endif A = math_mul33x33(Fg_new,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp - Fi = constitutive_getDamageStrain(g,i,e) ! intermediate configuration, assume decomposition as F = Fe Fi Fp - Fi = math_mul33x33(Fi,constitutive_getThermalStrain(g,i,e)) ! Fi = Ft Fd - Ci = math_mul33x33(math_transpose33(Fi),Fi) ! non-plastic stretch tensor (neglecting elastic contributions) - invFi = math_inv33(Fi) - detFi = math_det33(Fi) + Fi_current = constitutive_getFi0(g,i,e) ! intermediate configuration, assume decomposition as F = Fe Fi Fp + invFi_current = math_inv33(Fi_current) + if (all(invFi_current == 0.0_pReal)) then ! ... failed? +#ifndef _OPENMP + if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then + write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on inversion of Fi_current at el (elFE) ip g ',& + e,'(',mesh_element(1,e),')',i,g + if (iand(debug_level(debug_crystallite), debug_levelExtensive) > 0_pInt) & + write(6,'(/,a,/,3(12x,3(f12.7,1x)/))') '<< CRYST >> Fp_current',math_transpose33(Fi_current(1:3,1:3)) + endif +#endif + return + endif !* start LpLoop with normal step length NiterationStress = 0_pInt - jacoCounter = 0_pInt - steplength0 = 1.0_pReal - steplength = steplength0 - residuum_old = 0.0_pReal + jacoCounter = 0_pInt + steplength0 = 1.0_pReal + steplength = steplength0 + residuum_old = 0.0_pReal LpLoop: do NiterationStress = NiterationStress + 1_pInt @@ -3449,17 +3525,106 @@ logical function crystallite_integrateStress(& #endif return endif loopsExeced + B = math_I3 - dt*Lpguess + NiterationStressI = 0_pInt + Liguess = 0.0_pReal + Liguess_old = Liguess + jacoICounter = 0_pInt + steplengthI0 = 1.0_pReal + steplengthI = steplengthI0 + residuumI_old = 0.0_pReal + LiLoop: do ! inner stress integration loop for consistency with Fi + NiterationStressI = NiterationStressI + 1_pInt + IloopsExeced: if (NiterationStressI > nStress) then +#ifndef _OPENMP + if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) & + write(6,'(a,i3,a,i8,1x,a,i8,a,1x,i2,1x,i3,/)') '<< CRYST >> integrateStress reached inelastic loop limit',nStress, & + ' at el (elFE) ip g ', e,mesh_element(1,e),i,g +#endif + return + endif IloopsExeced + + !* calculate (elastic) 2nd Piola--Kirchhoff stress tensor and its tangent from constitutive law + invFi = math_mul33x33(invFi_current,math_I3 - dt*Liguess) + detInvFi = math_det33(invFi) + Fi = math_inv33(invFi) + Fe = math_mul33x33(math_mul33x33(A,B),invFi) ! current elastic deformation tensor + call constitutive_TandItsTangent(Tstar_unloaded, dT_dFe3333, Fe, g,i,e) ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative in unloaded configuration + Tstar = math_mul33x33(invFi, & + math_mul33x33(Tstar_unloaded,math_transpose33(invFi)))/detInvFi ! push Tstar forward from unloaded to plastic (lattice) configuration + Tstar_v = math_Mandel33to6(Tstar) - !* calculate (elastic) 2nd Piola--Kirchhoff stress tensor and its tangent from constitutive law + !* calculate intermediate velocity gradient and its tangent from constitutive law - B = math_I3 - dt*Lpguess - Fe = math_mul33x33(math_mul33x33(A,B),invFi) ! current elastic deformation tensor - call constitutive_TandItsTangent(Tstar, dT_dFe3333, Fe, g,i,e) ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative in unloaded configuration - Tstar = math_mul33x33(Ci, math_mul33x33(invFi, & - math_mul33x33(Tstar,math_transpose33(invFi))))/detFi ! push Tstar forward to plastic (lattice) configuration - Tstar_v = math_Mandel33to6(Tstar) + call constitutive_LiAndItsTangent(Li_constitutive, dLi_dT_constitutive, Tstar_v, Lpguess, & + g, i, e) + !* update current residuum and check for convergence of loop + + aTol = max(rTol_crystalliteStress * max(math_norm33(Liguess),math_norm33(Li_constitutive)), & ! absolute tolerance from largest acceptable relative error + aTol_crystalliteStress) ! minimum lower cutoff + residuumI = Liguess - Li_constitutive + if (any(residuumI /= residuumI)) then ! NaN in residuum... + return ! ...me = .false. to inform integrator about problem + elseif (math_norm33(residuumI) < aTol) then ! converged if below absolute tolerance + exit LiLoop ! ...leave iteration loop + elseif (math_norm33(residuumI) < math_norm33(residuumI_old) .or. NiterationStress == 1_pInt ) then ! not converged, but improved norm of residuum (always proceed in first iteration)... + residuumI_old = residuumI ! ...remember old values and... + Liguess_old = Liguess + steplengthI = steplengthI0 ! ...proceed with normal step length (calculate new search direction) + else ! not converged and residuum not improved... + steplengthI = 0.5_pReal * steplengthI ! ...try with smaller step length in same direction + Liguess = Liguess_old + steplengthI * deltaLi + cycle LiLoop + endif + + !* calculate Jacobian for correction term + + if (mod(jacoICounter, iJacoLpresiduum) == 0_pInt) then + dInvFi_dLp3333 = 0.0_pReal + do o=1_pInt,3_pInt + dInvFi_dLp3333(1:3,o,1:3,o) = invFi_current + enddo + dInvFi_dLp3333 = -dt * dInvFi_dLp3333 + dFe_dInvFi3333 = 0.0_pReal + temp_33 = math_mul33x33(A,B) + do o=1_pInt,3_pInt + dFe_dInvFi3333(1:3,o,1:3,o) = temp_33 + enddo + dT_dInvFi3333 = 0.0_pReal + do o=1_pInt,3_pInt; do p=1_pInt,3_pInt + dT_dInvFi3333(1:3,1:3,p,o) = -Tstar * Fi(o,p) + enddo; enddo + temp_33 = math_mul33x33(invFi,Tstar_unloaded)/detInvFi + do o=1_pInt,3_pInt + dT_dInvFi3333(1:3,o,o,1:3) = dT_dInvFi3333(1:3,o,o,1:3) + temp_33 + dT_dInvFi3333(o,1:3,o,1:3) = dT_dInvFi3333(o,1:3,o,1:3) + temp_33 + enddo + temp_3333 = math_mul3333xx3333(dT_dFe3333,dFe_dInvFi3333) + do o=1_pInt,3_pInt; do p=1_pInt,3_pInt + dT_dInvFi3333(1:3,1:3,o,p) = dT_dInvFi3333(1:3,1:3,o,p) + & + math_mul33x33(math_mul33x33(invFi,temp_3333(1:3,1:3,o,p)), & + math_transpose33(invFi))/detInvFi + enddo; enddo + dRI_dLp = math_identity2nd(9_pInt) - & + math_mul99x99(dLi_dT_constitutive,math_Plain3333to99(math_mul3333xx3333(dT_dInvFi3333,dInvFi_dLp3333))) + dRI_dLp2 = dRI_dLp ! will be overwritten in first call to LAPACK routine + work = math_plain33to9(residuumI) +#if(FLOAT==8) + call dgesv(9,1,dRI_dLp2,9,ipiv,work,9,ierr) ! solve dR/dLp * delta Lp = -res for dR/dLp +#elif(FLOAT==4) + call sgesv(9,1,dRI_dLp2,9,ipiv,work,9,ierr) ! solve dR/dLp * delta Lp = -res for dR/dLp +#endif + if (ierr /= 0_pInt) then + return + endif + deltaLi = - math_plain9to33(work) + endif + jacoICounter = jacoICounter + 1_pInt ! increase counter for jaco update + + Liguess = Liguess + steplengthI * deltaLi + enddo LiLoop !* calculate plastic velocity gradient and its tangent from constitutive law @@ -3524,7 +3689,7 @@ logical function crystallite_integrateStress(& if (mod(jacoCounter, iJacoLpresiduum) == 0_pInt) then dFe_dLp3333 = 0.0_pReal do o=1_pInt,3_pInt; do p=1_pInt,3_pInt - dFe_dLp3333(o,p,1:3,p) = A(o,1:3) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) delta(l,j) + dFe_dLp3333(o,1:3,p,1:3) = A(o,p)*math_transpose33(invFi) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j) enddo; enddo dFe_dLp3333 = -dt * dFe_dLp3333 dFe_dLp = math_Plain3333to99(dFe_dLp3333) @@ -3594,7 +3759,7 @@ logical function crystallite_integrateStress(& crystallite_P(1:3,1:3,g,i,e) = math_mul33x33(math_mul33x33(Fe_new,Fi), & math_mul33x33(math_Mandel6to33(Tstar_v), & - math_transpose33(invFp_new)))*detFi + math_transpose33(invFp_new)))/detInvFi !* store local values in global variables diff --git a/code/damage_anisoBrittle.f90 b/code/damage_anisoBrittle.f90 index 370275d7a..2b5ab3b8d 100644 --- a/code/damage_anisoBrittle.f90 +++ b/code/damage_anisoBrittle.f90 @@ -32,7 +32,8 @@ module damage_anisoBrittle damage_anisoBrittle_Ncleavage !< number of cleavage systems per family real(pReal), dimension(:), allocatable, private :: & - damage_anisoBrittle_aTol, & + damage_anisoBrittle_aTol_damage, & + damage_anisoBrittle_aTol_disp, & damage_anisoBrittle_sdot_0, & damage_anisoBrittle_N @@ -54,10 +55,13 @@ module damage_anisoBrittle damage_anisoBrittle_stateInit, & damage_anisoBrittle_aTolState, & damage_anisoBrittle_dotState, & + damage_anisoBrittle_LdAndItsTangent, & + damage_anisoBrittle_getFd, & + damage_anisoBrittle_getFd0, & + damage_anisoBrittle_getPartionedFd0, & damage_anisoBrittle_getDamage, & damage_anisoBrittle_putLocalDamage, & damage_anisoBrittle_getLocalDamage, & - damage_anisoBrittle_getDamageStrain, & damage_anisoBrittle_postResults contains @@ -139,9 +143,10 @@ subroutine damage_anisoBrittle_init(fileUnit) allocate(damage_anisoBrittle_Noutput(maxNinstance), source=0_pInt) allocate(damage_anisoBrittle_critDisp(lattice_maxNcleavageFamily,maxNinstance), source=0.0_pReal) allocate(damage_anisoBrittle_critLoad(lattice_maxNcleavageFamily,maxNinstance), source=0.0_pReal) - allocate(damage_anisoBrittle_Ncleavage(lattice_maxNcleavageFamily,maxNinstance), source=0_pInt) - allocate(damage_anisoBrittle_totalNcleavage(maxNinstance), source=0_pInt) - allocate(damage_anisoBrittle_aTol(maxNinstance), source=0.0_pReal) + allocate(damage_anisoBrittle_Ncleavage(lattice_maxNcleavageFamily,maxNinstance),source=0_pInt) + allocate(damage_anisoBrittle_totalNcleavage(maxNinstance), source=0_pInt) + allocate(damage_anisoBrittle_aTol_damage(maxNinstance), source=0.0_pReal) + allocate(damage_anisoBrittle_aTol_disp(maxNinstance), source=0.0_pReal) allocate(damage_anisoBrittle_sdot_0(maxNinstance), source=0.0_pReal) allocate(damage_anisoBrittle_N(maxNinstance), source=0.0_pReal) @@ -177,7 +182,10 @@ subroutine damage_anisoBrittle_init(fileUnit) end select case ('atol_damage') - damage_anisoBrittle_aTol(instance) = IO_floatValue(line,positions,2_pInt) + damage_anisoBrittle_aTol_damage(instance) = IO_floatValue(line,positions,2_pInt) + + case ('atol_disp') + damage_anisoBrittle_aTol_disp(instance) = IO_floatValue(line,positions,2_pInt) case ('sdot_0') damage_anisoBrittle_sdot_0(instance) = IO_floatValue(line,positions,2_pInt) @@ -226,8 +234,8 @@ subroutine damage_anisoBrittle_init(fileUnit) enddo outputsLoop ! Determine size of state array sizeDotState = 2_pInt + & ! viscous and non-viscous damage values - 9_pInt + & ! damage deformation gradient - 3_pInt*damage_anisoBrittle_totalNcleavage(instance) ! opening on each damage system + 3_pInt*damage_anisoBrittle_totalNcleavage(instance) + & ! opening on each damage system + 9_pInt ! Fd sizeState = sizeDotState damageState(phase)%sizeState = sizeState @@ -252,7 +260,7 @@ subroutine damage_anisoBrittle_init(fileUnit) if (any(numerics_integrator == 5_pInt)) & allocate(damageState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal) - call damage_anisoBrittle_stateInit(phase) + call damage_anisoBrittle_stateInit(phase,instance) call damage_anisoBrittle_aTolState(phase,instance) endif @@ -262,27 +270,22 @@ end subroutine damage_anisoBrittle_init !-------------------------------------------------------------------------------------------------- !> @brief sets the relevant state values for a given instance of this damage !-------------------------------------------------------------------------------------------------- -subroutine damage_anisoBrittle_stateInit(phase) +subroutine damage_anisoBrittle_stateInit(phase,instance) use material, only: & damageState + use math, only: & + math_I3 implicit none - integer(pInt), intent(in) :: phase !< number specifying the phase of the damage + integer(pInt), intent(in) :: phase, instance !< number specifying the phase of the damage real(pReal), dimension(damageState(phase)%sizeState) :: tempState - tempState = 0.0_pReal tempState(1) = 1.0_pReal tempState(2) = 1.0_pReal - tempState(3) = 1.0_pReal - tempState(4) = 0.0_pReal - tempState(5) = 0.0_pReal - tempState(6) = 0.0_pReal - tempState(7) = 1.0_pReal - tempState(8) = 0.0_pReal - tempState(9) = 0.0_pReal - tempState(10) = 0.0_pReal - tempState(11) = 1.0_pReal + tempState(3:2+3*damage_anisoBrittle_totalNcleavage(instance)) = 0.0_pReal + tempState(3*damage_anisoBrittle_totalNcleavage(instance)+3: & + 3*damage_anisoBrittle_totalNcleavage(instance)+11) = reshape(math_I3, shape=[9]) damageState(phase)%state = spread(tempState,2,size(damageState(phase)%state(1,:))) damageState(phase)%state0 = damageState(phase)%state damageState(phase)%partionedState0 = damageState(phase)%state @@ -301,7 +304,11 @@ subroutine damage_anisoBrittle_aTolState(phase,instance) instance ! number specifying the current instance of the damage real(pReal), dimension(damageState(phase)%sizeState) :: tempTol - tempTol = damage_anisoBrittle_aTol(instance) + tempTol(1) = damage_anisoBrittle_aTol_damage(instance) + tempTol(2) = damage_anisoBrittle_aTol_damage(instance) + tempTol(3:2+3*damage_anisoBrittle_totalNcleavage(instance)) = damage_anisoBrittle_aTol_disp(instance) + tempTol(3*damage_anisoBrittle_totalNcleavage(instance)+3: & + 3*damage_anisoBrittle_totalNcleavage(instance)+11) = damage_anisoBrittle_aTol_disp(instance) damageState(phase)%aTolState = tempTol end subroutine damage_anisoBrittle_aTolState @@ -330,13 +337,13 @@ subroutine damage_anisoBrittle_dotState(Tstar_v,ipc, ip, el) el !< element real(pReal), intent(in), dimension(6) :: & Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel) + real(pReal), dimension(3,3) :: & + Ld, Fd integer(pInt) :: & phase, & constituent, & instance, & f, i, index, index_myFamily - real(pReal), dimension(3,3) :: & - Ld, Fd real(pReal) :: & traction_d, traction_t, traction_n, traction_crit, & opening, & @@ -350,20 +357,20 @@ subroutine damage_anisoBrittle_dotState(Tstar_v,ipc, ip, el) (1.0_pReal/lattice_DamageMobility(phase))* & (damageState(phase)%state(2,constituent) - & damageState(phase)%state(1,constituent)) - + damageState(phase)%dotState(2,constituent) = 0.0_pReal nonLocalFactor = 1.0_pReal + & - (damage_anisoBrittle_getDamage(ipc, ip, el) - & - damageState(phase)%state(1,constituent)) - Ld = 0.0_pReal - index = 11_pInt + (damage_anisoBrittle_getDamage (ipc, ip, el) - & + damage_anisoBrittle_getLocalDamage(ipc, ip, el)) + Ld = 0.0_pReal + index = 2_pInt do f = 1_pInt,lattice_maxNcleavageFamily index_myFamily = sum(lattice_NcleavageSystem(1:f-1_pInt,phase)) ! at which index starts my family do i = 1_pInt,damage_anisoBrittle_Ncleavage(f,instance) ! process each (active) cleavage system in family - opening = math_norm3(damageState(phase)%state(index+1:index+3,constituent)) + opening = math_norm3(damageState(phase)%state0(index+1:index+3,constituent)) traction_d = dot_product(Tstar_v,lattice_Scleavage_v(1:6,1,index_myFamily+i,phase)) traction_t = dot_product(Tstar_v,lattice_Scleavage_v(1:6,2,index_myFamily+i,phase)) traction_n = dot_product(Tstar_v,lattice_Scleavage_v(1:6,3,index_myFamily+i,phase)) - traction_crit = (1.0_pReal - opening/damage_anisoBrittle_critDisp(f,instance))* & + traction_crit = max(0.001_pReal,(1.0_pReal - opening/damage_anisoBrittle_critDisp(f,instance)))* & damage_anisoBrittle_critLoad(f,instance)*nonLocalFactor damageState(phase)%dotState(index+1,constituent) = & sign(1.0_pReal,traction_d)* & @@ -393,11 +400,218 @@ subroutine damage_anisoBrittle_dotState(Tstar_v,ipc, ip, el) index = index + 3_pInt enddo enddo - Fd = reshape(damageState(phase)%state(3:11,constituent),shape=[3,3]) - damageState(phase)%dotState(3:11,constituent) = reshape(math_mul33x33(Ld,Fd),shape=[9]) + Fd = reshape(damageState(phase)%state(3*damage_anisoBrittle_totalNcleavage(instance)+3: & + 3*damage_anisoBrittle_totalNcleavage(instance)+11,constituent), & + shape=[3,3]) + damageState(phase)%dotState(3*damage_anisoBrittle_totalNcleavage(instance)+3: & + 3*damage_anisoBrittle_totalNcleavage(instance)+11,constituent) = & + reshape(math_mul33x33(Ld,Fd), shape=[9]) end subroutine damage_anisoBrittle_dotState + +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the velocity gradient +!-------------------------------------------------------------------------------------------------- +subroutine damage_anisoBrittle_LdAndItsTangent(Ld, dLd_dTstar, Tstar_v, ipc, ip, el) + use material, only: & + mappingConstitutive, & + phase_damageInstance, & + damageState + use math, only: & + math_norm3, & + math_Plain3333to99 + use lattice, only: & + lattice_Scleavage, & + lattice_Scleavage_v, & + lattice_maxNcleavageFamily, & + lattice_NcleavageSystem + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal), intent(in), dimension(6) :: & + Tstar_v !< 2nd Piola-Kirchhoff stress + real(pReal), intent(out), dimension(3,3) :: & + Ld !< damage velocity gradient + real(pReal), intent(out), dimension(9,9) :: & + dLd_dTstar !< derivative of Ld with respect to Tstar (2nd-order tensor) + real(pReal), dimension(3,3,3,3) :: & + dLd_dTstar3333 !< derivative of Ld with respect to Tstar (4th-order tensor) + integer(pInt) :: & + phase, & + constituent, & + instance, & + f, i, index, index_myFamily, k, l, m, n + real(pReal) :: & + traction_d, traction_t, traction_n, traction_crit, & + udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt, & + opening, & + nonLocalFactor + + phase = mappingConstitutive(2,ipc,ip,el) + constituent = mappingConstitutive(1,ipc,ip,el) + instance = phase_damageInstance(phase) + + Ld = 0.0_pReal + dLd_dTstar3333 = 0.0_pReal + nonLocalFactor = 1.0_pReal + & + (damage_anisoBrittle_getDamage (ipc, ip, el) - & + damage_anisoBrittle_getLocalDamage(ipc, ip, el)) + index = 2_pInt + do f = 1_pInt,lattice_maxNcleavageFamily + index_myFamily = sum(lattice_NcleavageSystem(1:f-1_pInt,phase)) ! at which index starts my family + do i = 1_pInt,damage_anisoBrittle_Ncleavage(f,instance) ! process each (active) cleavage system in family + opening = math_norm3(damageState(phase)%state0(index+1:index+3,constituent)) + traction_d = dot_product(Tstar_v,lattice_Scleavage_v(1:6,1,index_myFamily+i,phase)) + traction_t = dot_product(Tstar_v,lattice_Scleavage_v(1:6,2,index_myFamily+i,phase)) + traction_n = dot_product(Tstar_v,lattice_Scleavage_v(1:6,3,index_myFamily+i,phase)) + traction_crit = max(0.001_pReal,(1.0_pReal - opening/damage_anisoBrittle_critDisp(f,instance)))* & + damage_anisoBrittle_critLoad(f,instance)*nonLocalFactor + udotd = & + sign(1.0_pReal,traction_d)* & + damage_anisoBrittle_sdot_0(instance)* & + (abs(traction_d)/traction_crit)**damage_anisoBrittle_N(instance) + if (udotd /= 0.0_pReal) then + Ld = Ld + udotd*lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase) + dudotd_dt = udotd*damage_anisoBrittle_N(instance)/traction_d + forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & + dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + & + dudotd_dt*lattice_Scleavage(k,l,1,index_myFamily+i,phase)* & + lattice_Scleavage(m,n,1,index_myFamily+i,phase) + endif + + udott = & + sign(1.0_pReal,traction_t)* & + damage_anisoBrittle_sdot_0(instance)* & + (abs(traction_t)/traction_crit)**damage_anisoBrittle_N(instance) + if (udott /= 0.0_pReal) then + Ld = Ld + udott*lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase) + dudott_dt = udott*damage_anisoBrittle_N(instance)/traction_t + forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & + dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + & + dudott_dt*lattice_Scleavage(k,l,2,index_myFamily+i,phase)* & + lattice_Scleavage(m,n,2,index_myFamily+i,phase) + endif + + udotn = & + sign(1.0_pReal,traction_n)* & + damage_anisoBrittle_sdot_0(instance)* & + (abs(traction_n)/traction_crit)**damage_anisoBrittle_N(instance) + if (udotn /= 0.0_pReal) then + Ld = Ld + udotn*lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase) + dudotn_dt = udotn*damage_anisoBrittle_N(instance)/traction_n + forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & + dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + & + dudotn_dt*lattice_Scleavage(k,l,3,index_myFamily+i,phase)* & + lattice_Scleavage(m,n,3,index_myFamily+i,phase) + endif + + index = index + 3_pInt + enddo + enddo + dLd_dTstar = math_Plain3333to99(dLd_dTstar3333) + +end subroutine damage_anisoBrittle_LdAndItsTangent + +!-------------------------------------------------------------------------------------------------- +!> @brief returns local damage deformation gradient +!-------------------------------------------------------------------------------------------------- +pure function damage_anisoBrittle_getFd(ipc, ip, el) + use material, only: & + mappingConstitutive, & + phase_damageInstance, & + damageState + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + integer(pInt) :: & + phase, & + constituent, & + instance + real(pReal), dimension(3,3) :: & + damage_anisoBrittle_getFd + + phase = mappingConstitutive(2,ipc,ip,el) + constituent = mappingConstitutive(1,ipc,ip,el) + instance = phase_damageInstance(phase) + + damage_anisoBrittle_getFd = & + reshape(damageState(phase)%state(3*damage_anisoBrittle_totalNcleavage(instance)+3: & + 3*damage_anisoBrittle_totalNcleavage(instance)+11,constituent), & + shape=[3,3]) + +end function damage_anisoBrittle_getFd + +!-------------------------------------------------------------------------------------------------- +!> @brief returns local damage deformation gradient +!-------------------------------------------------------------------------------------------------- +pure function damage_anisoBrittle_getFd0(ipc, ip, el) + use material, only: & + mappingConstitutive, & + phase_damageInstance, & + damageState + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + integer(pInt) :: & + phase, & + constituent, & + instance + real(pReal), dimension(3,3) :: & + damage_anisoBrittle_getFd0 + + phase = mappingConstitutive(2,ipc,ip,el) + constituent = mappingConstitutive(1,ipc,ip,el) + instance = phase_damageInstance(phase) + + damage_anisoBrittle_getFd0 = & + reshape(damageState(phase)%subState0(3*damage_anisoBrittle_totalNcleavage(instance)+3: & + 3*damage_anisoBrittle_totalNcleavage(instance)+11,constituent), & + shape=[3,3]) + +end function damage_anisoBrittle_getFd0 + +!-------------------------------------------------------------------------------------------------- +!> @brief returns local damage deformation gradient +!-------------------------------------------------------------------------------------------------- +pure function damage_anisoBrittle_getPartionedFd0(ipc, ip, el) + use material, only: & + mappingConstitutive, & + phase_damageInstance, & + damageState + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + integer(pInt) :: & + phase, & + constituent, & + instance + real(pReal), dimension(3,3) :: & + damage_anisoBrittle_getPartionedFd0 + + phase = mappingConstitutive(2,ipc,ip,el) + constituent = mappingConstitutive(1,ipc,ip,el) + instance = phase_damageInstance(phase) + + damage_anisoBrittle_getPartionedFd0 = & + reshape(damageState(phase)%partionedState0(3*damage_anisoBrittle_totalNcleavage(instance)+3: & + 3*damage_anisoBrittle_totalNcleavage(instance)+11,constituent), & + shape=[3,3]) + +end function damage_anisoBrittle_getPartionedFd0 + !-------------------------------------------------------------------------------------------------- !> @brief returns damage !-------------------------------------------------------------------------------------------------- @@ -473,29 +687,6 @@ function damage_anisoBrittle_getLocalDamage(ipc, ip, el) end function damage_anisoBrittle_getLocalDamage -!-------------------------------------------------------------------------------------------------- -!> @brief returns local damage deformation gradient -!-------------------------------------------------------------------------------------------------- -function damage_anisoBrittle_getDamageStrain(ipc, ip, el) - use material, only: & - mappingConstitutive, & - phase_damageInstance, & - damageState - - implicit none - integer(pInt), intent(in) :: & - ipc, & !< grain number - ip, & !< integration point number - el !< element number - real(pReal), dimension(3,3) :: & - damage_anisoBrittle_getDamageStrain - - damage_anisoBrittle_getDamageStrain = & - reshape(damageState(mappingConstitutive(2,ipc,ip,el))%state(3:11,mappingConstitutive(1,ipc,ip,el)), & - shape=[3,3]) - -end function damage_anisoBrittle_getDamageStrain - !-------------------------------------------------------------------------------------------------- !> @brief return array of constitutive results !-------------------------------------------------------------------------------------------------- diff --git a/code/thermal_adiabatic.f90 b/code/thermal_adiabatic.f90 index 703981a99..f3697a6b3 100644 --- a/code/thermal_adiabatic.f90 +++ b/code/thermal_adiabatic.f90 @@ -40,6 +40,10 @@ module thermal_adiabatic thermal_adiabatic_stateInit, & thermal_adiabatic_aTolState, & thermal_adiabatic_dotState, & + thermal_adiabatic_LTAndItsTangent, & + thermal_adiabatic_getFT, & + thermal_adiabatic_getFT0, & + thermal_adiabatic_getPartionedFT0, & thermal_adiabatic_getTemperature, & thermal_adiabatic_putTemperature, & thermal_adiabatic_postResults @@ -283,6 +287,160 @@ subroutine thermal_adiabatic_dotState(Tstar_v, Lp, ipc, ip, el) end subroutine thermal_adiabatic_dotState +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the velocity gradient +!-------------------------------------------------------------------------------------------------- +subroutine thermal_adiabatic_LTAndItsTangent(LT, dLT_dTstar, Tstar_v, Lp, ipc, ip, el) + use lattice, only: & + lattice_massDensity, & + lattice_specificHeat + use material, only: & + mappingConstitutive, & + phase_thermalInstance, & + thermalState + use math, only: & + math_Plain3333to99, & + math_Mandel6to33, & + math_I3 + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal), intent(in), dimension(6) :: & + Tstar_v !< 2nd Piola-Kirchhoff stress + real(pReal), intent(in), dimension(3,3) :: & + Lp !< plastic velocity gradient + real(pReal), intent(out), dimension(3,3) :: & + LT !< thermal velocity gradient + real(pReal), intent(out), dimension(9,9) :: & + dLT_dTstar !< derivative of LT with respect to Tstar (2nd-order tensor) + real(pReal), dimension(3,3,3,3) :: & + dLT_dTstar3333 !< derivative of LT with respect to Tstar (4th-order tensor) + integer(pInt) :: & + phase, & + constituent, & + i, j, k, l + real(pReal) :: & + Tdot + + phase = mappingConstitutive(2,ipc,ip,el) + constituent = mappingConstitutive(1,ipc,ip,el) + + Tdot = 0.95_pReal & + * sum(abs(math_Mandel6to33(Tstar_v))*Lp) & + / (lattice_massDensity(phase)*lattice_specificHeat(phase)) + + LT = Tdot*math_I3 + dLT_dTstar3333 = 0.0_pReal + forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt,k=1_pInt:3_pInt,l=1_pInt:3_pInt) & + dLT_dTstar3333(i,j,k,l) = dLT_dTstar3333(i,j,k,l) + Lp(k,l)*math_I3(i,j) + + dLT_dTstar3333 = 0.95_pReal*dLT_dTstar3333/(lattice_massDensity(phase)*lattice_specificHeat(phase)) + dLT_dTstar = math_Plain3333to99(dLT_dTstar3333) + +end subroutine thermal_adiabatic_LTAndItsTangent + +!-------------------------------------------------------------------------------------------------- +!> @brief returns local thermal deformation gradient +!-------------------------------------------------------------------------------------------------- +pure function thermal_adiabatic_getFT(ipc, ip, el) + use material, only: & + mappingConstitutive, & + thermalState + use math, only: & + math_I3 + use lattice, only: & + lattice_referenceTemperature, & + lattice_thermalExpansion33 + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + integer(pInt) :: & + phase, & + constituent + real(pReal), dimension(3,3) :: & + thermal_adiabatic_getFT + + phase = mappingConstitutive(2,ipc,ip,el) + constituent = mappingConstitutive(1,ipc,ip,el) + + thermal_adiabatic_getFT = math_I3 + & + lattice_thermalExpansion33(1:3,1:3,phase)* & + (thermalState(phase)%state(1,constituent) - lattice_referenceTemperature(phase)) + +end function thermal_adiabatic_getFT + +!-------------------------------------------------------------------------------------------------- +!> @brief returns local thermal deformation gradient +!-------------------------------------------------------------------------------------------------- +pure function thermal_adiabatic_getFT0(ipc, ip, el) + use material, only: & + mappingConstitutive, & + thermalState + use math, only: & + math_I3 + use lattice, only: & + lattice_referenceTemperature, & + lattice_thermalExpansion33 + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + integer(pInt) :: & + phase, & + constituent + real(pReal), dimension(3,3) :: & + thermal_adiabatic_getFT0 + + phase = mappingConstitutive(2,ipc,ip,el) + constituent = mappingConstitutive(1,ipc,ip,el) + + thermal_adiabatic_getFT0 = math_I3 + & + lattice_thermalExpansion33(1:3,1:3,phase)* & + (thermalState(phase)%subState0(1,constituent) - lattice_referenceTemperature(phase)) + +end function thermal_adiabatic_getFT0 + +!-------------------------------------------------------------------------------------------------- +!> @brief returns local thermal deformation gradient +!-------------------------------------------------------------------------------------------------- +pure function thermal_adiabatic_getPartionedFT0(ipc, ip, el) + use material, only: & + mappingConstitutive, & + thermalState + use math, only: & + math_I3 + use lattice, only: & + lattice_referenceTemperature, & + lattice_thermalExpansion33 + + implicit none + integer(pInt), intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + integer(pInt) :: & + phase, & + constituent + real(pReal), dimension(3,3) :: & + thermal_adiabatic_getPartionedFT0 + + phase = mappingConstitutive(2,ipc,ip,el) + constituent = mappingConstitutive(1,ipc,ip,el) + + thermal_adiabatic_getPartionedFT0 = math_I3 + & + lattice_thermalExpansion33(1:3,1:3,phase)* & + (thermalState(phase)%partionedState0(1,constituent) - lattice_referenceTemperature(phase)) + +end function thermal_adiabatic_getPartionedFT0 + !-------------------------------------------------------------------------------------------------- !> @brief returns temperature based on local damage model state layout !--------------------------------------------------------------------------------------------------