diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 358937e4b..71a386903 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -19,7 +19,7 @@ module constitutive implicit none private - integer(kind(ELASTICITY_undefined_ID)), dimension(:), allocatable, protected :: & + integer(kind(ELASTICITY_undefined_ID)), dimension(:), allocatable :: & phase_elasticity !< elasticity of each phase integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable :: & !ToDo: old intel compiler complains about protected @@ -52,8 +52,8 @@ module constitutive interface - module subroutine plastic_init - end subroutine plastic_init + module subroutine mech_init + end subroutine mech_init module subroutine damage_init end subroutine damage_init @@ -301,6 +301,21 @@ module constitutive C end subroutine source_damage_isoBrittle_deltaState + module subroutine constitutive_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el) + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + Fe, & !< elastic deformation gradient + Fi !< intermediate deformation gradient + real(pReal), intent(out), dimension(3,3) :: & + S !< 2nd Piola-Kirchhoff stress tensor + real(pReal), intent(out), dimension(3,3,3,3) :: & + dS_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient + dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient + end subroutine constitutive_SandItsTangents + module subroutine plastic_results end subroutine plastic_results @@ -385,14 +400,10 @@ subroutine constitutive_init integer :: & p, & !< counter in phase loop - s, & !< counter in source loop - stiffDegradationCtr + s !< counter in source loop class (tNode), pointer :: & debug_constitutive, & - phases, & - phase, & - elastic, & - stiffDegradation + phases debug_constitutive => config_debug%get('constitutive', defaultVal=emptyList) debugConstitutive%basic = debug_constitutive%contains('basic') @@ -402,52 +413,15 @@ subroutine constitutive_init debugConstitutive%ip = config_debug%get_asInt('integrationpoint',defaultVal = 1) debugConstitutive%grain = config_debug%get_asInt('grain',defaultVal = 1) -!------------------------------------------------------------------------------------------------- -! initialize elasticity (hooke) !ToDO: Maybe move to elastic submodule along with function homogenizedC? - phases => config_material%get('phase') - allocate(phase_elasticity(phases%length), source = ELASTICITY_undefined_ID) - allocate(phase_elasticityInstance(phases%length), source = 0) - allocate(phase_NstiffnessDegradations(phases%length),source=0) - - do p = 1, phases%length - phase => phases%get(p) - elastic => phase%get('elasticity') - if(elastic%get_asString('type') == 'hooke') then - phase_elasticity(p) = ELASTICITY_HOOKE_ID - else - call IO_error(200,ext_msg=elastic%get_asString('type')) - endif - stiffDegradation => phase%get('stiffness_degradation',defaultVal=emptyList) ! check for stiffness degradation mechanisms - phase_NstiffnessDegradations(p) = stiffDegradation%length - enddo - - allocate(phase_stiffnessDegradation(maxval(phase_NstiffnessDegradations),phases%length), & - source=STIFFNESS_DEGRADATION_undefined_ID) - - if(maxVal(phase_NstiffnessDegradations)/=0) then - do p = 1, phases%length - phase => phases%get(p) - stiffDegradation => phase%get('stiffness_degradation',defaultVal=emptyList) - do stiffDegradationCtr = 1, stiffDegradation%length - if(stiffDegradation%get_asString(stiffDegradationCtr) == 'damage') & - phase_stiffnessDegradation(stiffDegradationCtr,p) = STIFFNESS_DEGRADATION_damage_ID - enddo - enddo - endif - - do p = 1, phases%length - phase_elasticityInstance(p) = count(phase_elasticity(1:p) == phase_elasticity(p)) - enddo - - !-------------------------------------------------------------------------------------------------- ! initialize constitutive laws - call plastic_init + call mech_init call damage_init call thermal_init print'(/,a)', ' <<<+- constitutive init -+>>>'; flush(IO_STDOUT) + phases => config_material%get('phase') constitutive_source_maxSizeDotState = 0 PhaseLoop2:do p = 1,phases%length !-------------------------------------------------------------------------------------------------- @@ -666,80 +640,6 @@ pure function constitutive_initialFi(ipc, ip, el) end function constitutive_initialFi -!-------------------------------------------------------------------------------------------------- -!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to -!> the elastic/intermediate deformation gradients depending on the selected elastic law -!! (so far no case switch because only Hooke is implemented) -!-------------------------------------------------------------------------------------------------- -subroutine constitutive_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el) - - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - real(pReal), intent(in), dimension(3,3) :: & - Fe, & !< elastic deformation gradient - Fi !< intermediate deformation gradient - real(pReal), intent(out), dimension(3,3) :: & - S !< 2nd Piola-Kirchhoff stress tensor - real(pReal), intent(out), dimension(3,3,3,3) :: & - dS_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient - dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient - - call constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el) - - -end subroutine constitutive_SandItsTangents - - -!-------------------------------------------------------------------------------------------------- -!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to -!> the elastic and intermediate deformation gradients using Hooke's law -!-------------------------------------------------------------------------------------------------- -subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & - Fe, Fi, ipc, ip, el) - - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - real(pReal), intent(in), dimension(3,3) :: & - Fe, & !< elastic deformation gradient - Fi !< intermediate deformation gradient - real(pReal), intent(out), dimension(3,3) :: & - S !< 2nd Piola-Kirchhoff stress tensor in lattice configuration - real(pReal), intent(out), dimension(3,3,3,3) :: & - dS_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient - dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient - real(pReal), dimension(3,3) :: E - real(pReal), dimension(3,3,3,3) :: C - integer :: & - ho, & !< homogenization - d !< counter in degradation loop - integer :: & - i, j - - ho = material_homogenizationAt(el) - C = math_66toSym3333(constitutive_homogenizedC(ipc,ip,el)) - - DegradationLoop: do d = 1, phase_NstiffnessDegradations(material_phaseAt(ipc,el)) - degradationType: select case(phase_stiffnessDegradation(d,material_phaseAt(ipc,el))) - case (STIFFNESS_DEGRADATION_damage_ID) degradationType - C = C * damage(ho)%p(damageMapping(ho)%p(ip,el))**2 - end select degradationType - enddo DegradationLoop - - E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration - S = math_mul3333xx33(C,matmul(matmul(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration - - do i =1, 3;do j=1,3 - dS_dFe(i,j,1:3,1:3) = matmul(Fe,matmul(matmul(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko - dS_dFi(i,j,1:3,1:3) = 2.0_pReal*matmul(matmul(E,Fi),C(i,j,1:3,1:3)) !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn - enddo; enddo - -end subroutine constitutive_hooke_SandItsTangents - - !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- diff --git a/src/constitutive_plastic.f90 b/src/constitutive_plastic.f90 index bf6bc079e..6809443ca 100644 --- a/src/constitutive_plastic.f90 +++ b/src/constitutive_plastic.f90 @@ -193,14 +193,56 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief Initialize constitutive models for plasticity !-------------------------------------------------------------------------------------------------- -module subroutine plastic_init +module subroutine mech_init - integer :: p - class(tNode), pointer :: phases + integer :: & + p, & + stiffDegradationCtr + class(tNode), pointer :: & + phases, & + phase, & + elastic, & + stiffDegradation print'(/,a)', ' <<<+- constitutive_plastic init -+>>>' +!------------------------------------------------------------------------------------------------- +! initialize elasticity (hooke) !ToDO: Maybe move to elastic submodule along with function homogenizedC? phases => config_material%get('phase') + allocate(phase_elasticity(phases%length), source = ELASTICITY_undefined_ID) + allocate(phase_elasticityInstance(phases%length), source = 0) + allocate(phase_NstiffnessDegradations(phases%length),source=0) + + do p = 1, phases%length + phase => phases%get(p) + elastic => phase%get('elasticity') + if(elastic%get_asString('type') == 'hooke') then + phase_elasticity(p) = ELASTICITY_HOOKE_ID + else + call IO_error(200,ext_msg=elastic%get_asString('type')) + endif + stiffDegradation => phase%get('stiffness_degradation',defaultVal=emptyList) ! check for stiffness degradation mechanisms + phase_NstiffnessDegradations(p) = stiffDegradation%length + enddo + + allocate(phase_stiffnessDegradation(maxval(phase_NstiffnessDegradations),phases%length), & + source=STIFFNESS_DEGRADATION_undefined_ID) + + if(maxVal(phase_NstiffnessDegradations)/=0) then + do p = 1, phases%length + phase => phases%get(p) + stiffDegradation => phase%get('stiffness_degradation',defaultVal=emptyList) + do stiffDegradationCtr = 1, stiffDegradation%length + if(stiffDegradation%get_asString(stiffDegradationCtr) == 'damage') & + phase_stiffnessDegradation(stiffDegradationCtr,p) = STIFFNESS_DEGRADATION_damage_ID + enddo + enddo + endif + + do p = 1, phases%length + phase_elasticityInstance(p) = count(phase_elasticity(1:p) == phase_elasticity(p)) + enddo + allocate(plasticState(phases%length)) allocate(phase_plasticity(phases%length),source = PLASTICITY_undefined_ID) @@ -220,7 +262,7 @@ module subroutine plastic_init enddo -end subroutine plastic_init +end subroutine mech_init !-------------------------------------------------------------------------------------------------- @@ -248,6 +290,80 @@ module function plastic_active(plastic_label) result(active_plastic) end function plastic_active +!-------------------------------------------------------------------------------------------------- +!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to +!> the elastic/intermediate deformation gradients depending on the selected elastic law +!! (so far no case switch because only Hooke is implemented) +!-------------------------------------------------------------------------------------------------- +module subroutine constitutive_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + Fe, & !< elastic deformation gradient + Fi !< intermediate deformation gradient + real(pReal), intent(out), dimension(3,3) :: & + S !< 2nd Piola-Kirchhoff stress tensor + real(pReal), intent(out), dimension(3,3,3,3) :: & + dS_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient + dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient + + call constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el) + + +end subroutine constitutive_SandItsTangents + + +!-------------------------------------------------------------------------------------------------- +!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to +!> the elastic and intermediate deformation gradients using Hooke's law +!-------------------------------------------------------------------------------------------------- +subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & + Fe, Fi, ipc, ip, el) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + Fe, & !< elastic deformation gradient + Fi !< intermediate deformation gradient + real(pReal), intent(out), dimension(3,3) :: & + S !< 2nd Piola-Kirchhoff stress tensor in lattice configuration + real(pReal), intent(out), dimension(3,3,3,3) :: & + dS_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient + dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient + real(pReal), dimension(3,3) :: E + real(pReal), dimension(3,3,3,3) :: C + integer :: & + ho, & !< homogenization + d !< counter in degradation loop + integer :: & + i, j + + ho = material_homogenizationAt(el) + C = math_66toSym3333(constitutive_homogenizedC(ipc,ip,el)) + + DegradationLoop: do d = 1, phase_NstiffnessDegradations(material_phaseAt(ipc,el)) + degradationType: select case(phase_stiffnessDegradation(d,material_phaseAt(ipc,el))) + case (STIFFNESS_DEGRADATION_damage_ID) degradationType + C = C * damage(ho)%p(damageMapping(ho)%p(ip,el))**2 + end select degradationType + enddo DegradationLoop + + E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration + S = math_mul3333xx33(C,matmul(matmul(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration + + do i =1, 3;do j=1,3 + dS_dFe(i,j,1:3,1:3) = matmul(Fe,matmul(matmul(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko + dS_dFi(i,j,1:3,1:3) = 2.0_pReal*matmul(matmul(E,Fi),C(i,j,1:3,1:3)) !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn + enddo; enddo + +end subroutine constitutive_hooke_SandItsTangents + + !-------------------------------------------------------------------------------------------------- !> @brief calls microstructure function of the different plasticity constitutive models !--------------------------------------------------------------------------------------------------