From b5a10f23868f00a7fefc0914b4d97b9de49b5639 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Thu, 9 Jul 2020 01:01:08 +0200 Subject: [PATCH] sources and kinematics modules under submodules --- src/constitutive.f90 | 592 ++++++++++----------- src/constitutive_damage.f90 | 163 ++++++ src/constitutive_plastic.f90 | 274 ++++++++++ src/constitutive_plastic_disloUCLA.f90 | 2 +- src/constitutive_plastic_dislotwin.f90 | 2 +- src/constitutive_plastic_isotropic.f90 | 2 +- src/constitutive_plastic_kinehardening.f90 | 2 +- src/constitutive_plastic_none.f90 | 2 +- src/constitutive_plastic_nonlocal.f90 | 2 +- src/constitutive_plastic_phenopowerlaw.f90 | 2 +- src/constitutive_thermal.f90 | 112 ++++ src/damage_local.f90 | 40 +- src/damage_nonlocal.f90 | 41 +- src/kinematics_cleavage_opening.f90 | 20 +- src/kinematics_slipplane_opening.f90 | 20 +- src/kinematics_thermal_expansion.f90 | 57 +- src/source_damage_anisoBrittle.f90 | 29 +- src/source_damage_anisoDuctile.f90 | 33 +- src/source_damage_isoBrittle.f90 | 31 +- src/source_damage_isoDuctile.f90 | 31 +- src/source_thermal_dissipation.f90 | 19 +- src/source_thermal_externalheat.f90 | 21 +- src/thermal_adiabatic.f90 | 44 +- src/thermal_conduction.f90 | 43 +- 24 files changed, 937 insertions(+), 647 deletions(-) create mode 100644 src/constitutive_damage.f90 create mode 100644 src/constitutive_plastic.f90 create mode 100644 src/constitutive_thermal.f90 diff --git a/src/constitutive.f90 b/src/constitutive.f90 index d920c481f..5f509b0ba 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -16,15 +16,6 @@ module constitutive use lattice use discretization use geometry_plastic_nonlocal - use source_thermal_dissipation - use source_thermal_externalheat - use source_damage_isoBrittle - use source_damage_isoDuctile - use source_damage_anisoBrittle - use source_damage_anisoDuctile - use kinematics_cleavage_opening - use kinematics_slipplane_opening - use kinematics_thermal_expansion implicit none private @@ -35,114 +26,31 @@ module constitutive interface - module subroutine plastic_none_init - end subroutine plastic_none_init + module subroutine plastic_init + end subroutine plastic_init - module subroutine plastic_isotropic_init - end subroutine plastic_isotropic_init + module subroutine damage_init + end subroutine damage_init - module subroutine plastic_phenopowerlaw_init - end subroutine plastic_phenopowerlaw_init + module subroutine thermal_init + end subroutine thermal_init - module subroutine plastic_kinehardening_init - end subroutine plastic_kinehardening_init +! module pure function kinematics_thermal_expansion_initialStrain(homog,phase,offset) +! +! integer, intent(in) :: & +! phase, & +! homog, & +! offset +! end function kinematics_thermal_expansion_initialStrain - module subroutine plastic_dislotwin_init - end subroutine plastic_dislotwin_init - - module subroutine plastic_disloUCLA_init - end subroutine plastic_disloUCLA_init - - module subroutine plastic_nonlocal_init - end subroutine plastic_nonlocal_init - - - module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) - real(pReal), dimension(3,3), intent(out) :: & - Lp !< plastic velocity gradient - real(pReal), dimension(3,3,3,3), intent(out) :: & - dLp_dMp !< derivative of Lp with respect to the Mandel stress - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - end subroutine plastic_isotropic_LpAndItsTangent - - pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) - real(pReal), dimension(3,3), intent(out) :: & - Lp !< plastic velocity gradient - real(pReal), dimension(3,3,3,3), intent(out) :: & - dLp_dMp !< derivative of Lp with respect to the Mandel stress - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - end subroutine plastic_phenopowerlaw_LpAndItsTangent - - pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) - real(pReal), dimension(3,3), intent(out) :: & - Lp !< plastic velocity gradient - real(pReal), dimension(3,3,3,3), intent(out) :: & - dLp_dMp !< derivative of Lp with respect to the Mandel stress - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - integer, intent(in) :: & - instance, & - of - end subroutine plastic_kinehardening_LpAndItsTangent - - module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) - real(pReal), dimension(3,3), intent(out) :: & - Lp !< plastic velocity gradient - real(pReal), dimension(3,3,3,3), intent(out) :: & - dLp_dMp !< derivative of Lp with respect to the Mandel stress - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - T - integer, intent(in) :: & - instance, & - of - end subroutine plastic_dislotwin_LpAndItsTangent - - pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) - real(pReal), dimension(3,3), intent(out) :: & - Lp !< plastic velocity gradient - real(pReal), dimension(3,3,3,3), intent(out) :: & - dLp_dMp !< derivative of Lp with respect to the Mandel stress - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - T - integer, intent(in) :: & - instance, & - of - end subroutine plastic_disloUCLA_LpAndItsTangent - - module subroutine plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp, & - Mp,Temperature,instance,of,ip,el) - real(pReal), dimension(3,3), intent(out) :: & - Lp !< plastic velocity gradient - real(pReal), dimension(3,3,3,3), intent(out) :: & - dLp_dMp !< derivative of Lp with respect to the Mandel stress - - real(pReal), dimension(3,3), intent(in) :: & - Mp !< Mandel stress - real(pReal), intent(in) :: & - Temperature - integer, intent(in) :: & - instance, & - of, & - ip, & !< current integration point - el !< current element number - end subroutine plastic_nonlocal_LpAndItsTangent + module function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC) + real(pReal), dimension(6,6) :: & + homogenizedC + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + end function plastic_dislotwin_homogenizedC module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) @@ -158,6 +66,45 @@ module constitutive of end subroutine plastic_isotropic_LiAndItsTangent + module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) + + integer, intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal), intent(in), dimension(3,3) :: & + S + real(pReal), intent(out), dimension(3,3) :: & + Ld !< damage velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) + end subroutine kinematics_cleavage_opening_LiAndItsTangent + + module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) + + integer, intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal), intent(in), dimension(3,3) :: & + S + real(pReal), intent(out), dimension(3,3) :: & + Ld !< damage velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) + end subroutine kinematics_slipplane_opening_LiAndItsTangent + + module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, el) + + integer, intent(in) :: & + ipc, & !< grain number + ip, & !< integration point number + el !< element number + real(pReal), intent(out), dimension(3,3) :: & + Li !< thermal velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero) + end subroutine kinematics_thermal_expansion_LiAndItsTangent module subroutine plastic_isotropic_dotState(Mp,instance,of) real(pReal), dimension(3,3), intent(in) :: & @@ -220,31 +167,39 @@ module constitutive el !< current element number end subroutine plastic_nonlocal_dotState + + module subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el) - module subroutine plastic_dislotwin_dependentState(T,instance,of) - integer, intent(in) :: & - instance, & - of - real(pReal), intent(in) :: & - T - end subroutine plastic_dislotwin_dependentState - - module subroutine plastic_disloUCLA_dependentState(instance,of) - integer, intent(in) :: & - instance, & - of - end subroutine plastic_disloUCLA_dependentState - - module subroutine plastic_nonlocal_dependentState(F, Fp, instance, of, ip, el) - real(pReal), dimension(3,3), intent(in) :: & - F, & - Fp integer, intent(in) :: & - instance, & - of, & - ip, & - el - end subroutine plastic_nonlocal_dependentState + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + S + end subroutine source_damage_anisoBrittle_dotState + + module subroutine source_damage_anisoDuctile_dotState(ipc, ip, el) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + end subroutine source_damage_anisoDuctile_dotState + + module subroutine source_damage_isoDuctile_dotState(ipc, ip, el) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + end subroutine source_damage_isoDuctile_dotState + + module subroutine source_thermal_externalheat_dotState(phase, of) + + integer, intent(in) :: & + phase, & + of + end subroutine source_thermal_externalheat_dotState module subroutine plastic_kinehardening_deltaState(Mp,instance,of) @@ -265,24 +220,17 @@ module constitutive el end subroutine plastic_nonlocal_deltaState - - module function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC) - real(pReal), dimension(6,6) :: & - homogenizedC - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - end function plastic_dislotwin_homogenizedC - - module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e) + module subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el) integer, intent(in) :: & - instance, & - i, & - e - type(rotation), dimension(1,discretization_nIP,discretization_nElem), intent(in) :: & - orientation !< crystal orientation - end subroutine plastic_nonlocal_updateCompatibility + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + Fe + real(pReal), intent(in), dimension(6,6) :: & + C + end subroutine source_damage_isoBrittle_deltaState + module subroutine plastic_isotropic_results(instance,group) @@ -315,6 +263,88 @@ module constitutive character(len=*), intent(in) :: group end subroutine plastic_nonlocal_results + module subroutine source_damage_anisoBrittle_results(phase,group) + integer, intent(in) :: phase + character(len=*), intent(in) :: group + end subroutine source_damage_anisoBrittle_results + + module subroutine source_damage_anisoDuctile_results(phase,group) + integer, intent(in) :: phase + character(len=*), intent(in) :: group + end subroutine source_damage_anisoDuctile_results + + module subroutine source_damage_isoBrittle_results(phase,group) + integer, intent(in) :: phase + character(len=*), intent(in) :: group + end subroutine source_damage_isoBrittle_results + + module subroutine source_damage_isoDuctile_results(phase,group) + integer, intent(in) :: phase + character(len=*), intent(in) :: group + end subroutine source_damage_isoDuctile_results + + module subroutine plastic_dependentState(F, Fp, ipc, ip, el) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + F, & !< elastic deformation gradient + Fp !< plastic deformation gradient + end subroutine plastic_dependentState + + module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & + S, 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) :: & + S, & !< 2nd Piola-Kirchhoff stress + Fi !< intermediate deformation gradient + real(pReal), intent(out), dimension(3,3) :: & + Lp !< plastic velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLp_dS, & + dLp_dFi !< derivative of Lp with respect to Fi + + end subroutine plastic_LpAndItsTangents + + module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e) + integer, intent(in) :: & + instance, & + i, & + e + type(rotation), dimension(1,discretization_nIP,discretization_nElem), intent(in) :: & + orientation !< crystal orientation + end subroutine plastic_nonlocal_updateCompatibility + + module subroutine damage_source_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + phi + real(pReal), intent(inout) :: & + phiDot, & + dPhiDot_dPhi + end subroutine damage_source_getRateAndItsTangents + + module subroutine thermal_source_getRateAndItsTangents(TDot, dTDot_dT, T, Tstar, Lp, ip, el) + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + T + real(pReal), intent(in), dimension(:,:,:,:,:) :: & + Tstar, & + Lp + real(pReal), intent(inout) :: & + TDot, & + dTDot_dT + end subroutine thermal_source_getRateAndItsTangents + end interface @@ -332,16 +362,18 @@ module constitutive type(tDebugOptions) :: debugConstitutive public :: & - plastic_nonlocal_updateCompatibility, & constitutive_init, & constitutive_homogenizedC, & - constitutive_dependentState, & constitutive_LpAndItsTangents, & + constitutive_dependentState, & constitutive_LiAndItsTangents, & constitutive_initialFi, & constitutive_SandItsTangents, & constitutive_collectDotState, & constitutive_deltaState, & + plastic_nonlocal_updateCompatibility, & + constitutive_damage_getRateAndItsTangents, & + constitutive_thermal_getRateAndItsTangents, & constitutive_results contains @@ -368,31 +400,9 @@ subroutine constitutive_init !-------------------------------------------------------------------------------------------------- ! initialized plasticity - if (any(phase_plasticity == PLASTICITY_NONE_ID)) call plastic_none_init - if (any(phase_plasticity == PLASTICITY_ISOTROPIC_ID)) call plastic_isotropic_init - if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_init - if (any(phase_plasticity == PLASTICITY_KINEHARDENING_ID)) call plastic_kinehardening_init - if (any(phase_plasticity == PLASTICITY_DISLOTWIN_ID)) call plastic_dislotwin_init - if (any(phase_plasticity == PLASTICITY_DISLOUCLA_ID)) call plastic_disloucla_init - if (any(phase_plasticity == PLASTICITY_NONLOCAL_ID)) then - call plastic_nonlocal_init - else - call geometry_plastic_nonlocal_disable - endif -!-------------------------------------------------------------------------------------------------- -! initialize source mechanisms - if (any(phase_source == SOURCE_thermal_dissipation_ID)) call source_thermal_dissipation_init - if (any(phase_source == SOURCE_thermal_externalheat_ID)) call source_thermal_externalheat_init - if (any(phase_source == SOURCE_damage_isoBrittle_ID)) call source_damage_isoBrittle_init - if (any(phase_source == SOURCE_damage_isoDuctile_ID)) call source_damage_isoDuctile_init - if (any(phase_source == SOURCE_damage_anisoBrittle_ID)) call source_damage_anisoBrittle_init - if (any(phase_source == SOURCE_damage_anisoDuctile_ID)) call source_damage_anisoDuctile_init - -!-------------------------------------------------------------------------------------------------- -! initialize kinematic mechanisms - if (any(phase_kinematics == KINEMATICS_cleavage_opening_ID)) call kinematics_cleavage_opening_init - if (any(phase_kinematics == KINEMATICS_slipplane_opening_ID)) call kinematics_slipplane_opening_init - if (any(phase_kinematics == KINEMATICS_thermal_expansion_ID)) call kinematics_thermal_expansion_init + call plastic_init + call damage_init + call thermal_init write(6,'(/,a)') ' <<<+- constitutive init -+>>>'; flush(6) @@ -416,6 +426,43 @@ subroutine constitutive_init end subroutine constitutive_init +subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & + S, 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) :: & + S, & !< 2nd Piola-Kirchhoff stress + Fi !< intermediate deformation gradient + real(pReal), intent(out), dimension(3,3) :: & + Lp !< plastic velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLp_dS, & + dLp_dFi !< derivative of Lp with respect to Fi + + call plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S, Fi, ipc, ip, el) + +end subroutine constitutive_LpAndItsTangents + +!-------------------------------------------------------------------------------------------------- +!> @brief calls microstructure function of the different constitutive models +!-------------------------------------------------------------------------------------------------- +subroutine constitutive_dependentState(F, Fp, ipc, ip, el) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + F, & !< elastic deformation gradient + Fp !< plastic deformation gradient + + call plastic_dependentState(F,Fp,ipc,ip,el) + +end subroutine constitutive_dependentState + + !-------------------------------------------------------------------------------------------------- !> @brief returns the homogenize elasticity matrix !> ToDo: homogenizedC66 would be more consistent @@ -438,111 +485,6 @@ function constitutive_homogenizedC(ipc,ip,el) end function constitutive_homogenizedC -!-------------------------------------------------------------------------------------------------- -!> @brief calls microstructure function of the different constitutive models -!-------------------------------------------------------------------------------------------------- -subroutine constitutive_dependentState(F, Fp, ipc, ip, el) - - integer, intent(in) :: & - ipc, & !< component-ID of integration point - ip, & !< integration point - el !< element - real(pReal), intent(in), dimension(3,3) :: & - F, & !< elastic deformation gradient - Fp !< plastic deformation gradient - integer :: & - ho, & !< homogenization - tme, & !< thermal member position - instance, of - - ho = material_homogenizationAt(el) - tme = thermalMapping(ho)%p(ip,el) - of = material_phasememberAt(ipc,ip,el) - instance = phase_plasticityInstance(material_phaseAt(ipc,el)) - - plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) - case (PLASTICITY_DISLOTWIN_ID) plasticityType - call plastic_dislotwin_dependentState(temperature(ho)%p(tme),instance,of) - case (PLASTICITY_DISLOUCLA_ID) plasticityType - call plastic_disloUCLA_dependentState(instance,of) - case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_dependentState (F,Fp,instance,of,ip,el) - end select plasticityType - -end subroutine constitutive_dependentState - - -!-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the velocity gradient -! ToDo: Discuss whether it makes sense if crystallite handles the configuration conversion, i.e. -! Mp in, dLp_dMp out -!-------------------------------------------------------------------------------------------------- -subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & - S, 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) :: & - S, & !< 2nd Piola-Kirchhoff stress - Fi !< intermediate deformation gradient - real(pReal), intent(out), dimension(3,3) :: & - Lp !< plastic velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLp_dS, & - dLp_dFi !< derivative of Lp with respect to Fi - real(pReal), dimension(3,3,3,3) :: & - dLp_dMp !< derivative of Lp with respect to Mandel stress - real(pReal), dimension(3,3) :: & - Mp !< Mandel stress work conjugate with Lp - integer :: & - ho, & !< homogenization - tme !< thermal member position - integer :: & - i, j, instance, of - - ho = material_homogenizationAt(el) - tme = thermalMapping(ho)%p(ip,el) - - Mp = matmul(matmul(transpose(Fi),Fi),S) - of = material_phasememberAt(ipc,ip,el) - instance = phase_plasticityInstance(material_phaseAt(ipc,el)) - - plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) - - case (PLASTICITY_NONE_ID) plasticityType - Lp = 0.0_pReal - dLp_dMp = 0.0_pReal - - case (PLASTICITY_ISOTROPIC_ID) plasticityType - call plastic_isotropic_LpAndItsTangent (Lp,dLp_dMp,Mp,instance,of) - - case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType - call plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) - - case (PLASTICITY_KINEHARDENING_ID) plasticityType - call plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) - - case (PLASTICITY_NONLOCAL_ID) plasticityType - call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMp,Mp, temperature(ho)%p(tme),instance,of,ip,el) - - case (PLASTICITY_DISLOTWIN_ID) plasticityType - call plastic_dislotwin_LpAndItsTangent (Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,of) - - case (PLASTICITY_DISLOUCLA_ID) plasticityType - call plastic_disloucla_LpAndItsTangent (Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,of) - - end select plasticityType - - do i=1,3; do j=1,3 - dLp_dFi(i,j,1:3,1:3) = matmul(matmul(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + & - matmul(matmul(Fi,dLp_dMp(i,j,1:3,1:3)),S) - dLp_dS(i,j,1:3,1:3) = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi) - enddo; enddo - -end subroutine constitutive_LpAndItsTangents - - !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient ! ToDo: MD: S is Mi? @@ -649,7 +591,7 @@ pure function constitutive_initialFi(ipc, ip, el) homog = material_homogenizationAt(el) offset = thermalMapping(homog)%p(ip,el) constitutive_initialFi = & - constitutive_initialFi + kinematics_thermal_expansion_initialStrain(homog,phase,offset) + constitutive_initialFi !+ kinematics_thermal_expansion_initialStrain(homog,phase,offset) end select kinematicsType enddo KinematicsLoop @@ -890,43 +832,97 @@ function constitutive_deltaState(S, Fe, Fi, ipc, ip, el, phase, of) result(broke end function constitutive_deltaState +subroutine constitutive_thermal_getRateAndItsTangents(Tdot, dTDot_dT, T, Tstar, Lp, ip, el) + + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + T + real(pReal), intent(in), dimension(:,:,:,:,:) :: & + Tstar, & + Lp + real(pReal), intent(inout) :: & + Tdot, & + dTdot_dT + + call thermal_source_getRateAndItsTangents(Tdot, dTdot_dT, T, Tstar, Lp, ip, el) + +end subroutine constitutive_thermal_getRateAndItsTangents + +subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) + + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + phi + real(pReal), intent(inout) :: & + phiDot, & + dPhiDot_dPhi + + call damage_source_getRateAndItsTangents(phiDot,dPhiDot_dPhi,phi,ip,el) + +end subroutine constitutive_damage_getRateAndItsTangents + !-------------------------------------------------------------------------------------------------- !> @brief writes constitutive results to HDF5 output file !-------------------------------------------------------------------------------------------------- subroutine constitutive_results - integer :: p - character(len=pStringLen) :: group - do p=1,size(config_name_phase) + integer :: p,i + character(len=pStringLen) :: group,group_plastic,group_sources + + plasticityLoop: do p=1,size(config_name_phase) group = trim('current/constituent')//'/'//trim(config_name_phase(p)) call results_closeGroup(results_addGroup(group)) - group = trim(group)//'/plastic' + group_plastic = trim(group)//'/plastic' - call results_closeGroup(results_addGroup(group)) + call results_closeGroup(results_addGroup(group_plastic)) select case(phase_plasticity(p)) case(PLASTICITY_ISOTROPIC_ID) - call plastic_isotropic_results(phase_plasticityInstance(p),group) + call plastic_isotropic_results(phase_plasticityInstance(p),group_plastic) case(PLASTICITY_PHENOPOWERLAW_ID) - call plastic_phenopowerlaw_results(phase_plasticityInstance(p),group) + call plastic_phenopowerlaw_results(phase_plasticityInstance(p),group_plastic) case(PLASTICITY_KINEHARDENING_ID) - call plastic_kinehardening_results(phase_plasticityInstance(p),group) + call plastic_kinehardening_results(phase_plasticityInstance(p),group_plastic) case(PLASTICITY_DISLOTWIN_ID) - call plastic_dislotwin_results(phase_plasticityInstance(p),group) + call plastic_dislotwin_results(phase_plasticityInstance(p),group_plastic) case(PLASTICITY_DISLOUCLA_ID) - call plastic_disloUCLA_results(phase_plasticityInstance(p),group) + call plastic_disloUCLA_results(phase_plasticityInstance(p),group_plastic) case(PLASTICITY_NONLOCAL_ID) - call plastic_nonlocal_results(phase_plasticityInstance(p),group) + call plastic_nonlocal_results(phase_plasticityInstance(p),group_plastic) end select - enddo + sourceLoop: do i = 1, phase_Nsources(p) + group_sources = trim(group)//'/sources' + + call results_closeGroup(results_addGroup(group_sources)) + + sourceType: select case (phase_source(i,p)) + + case (SOURCE_damage_anisoBrittle_ID) sourceType + call source_damage_anisoBrittle_results(p,group_sources) + case (SOURCE_damage_anisoDuctile_ID) sourceType + call source_damage_anisoDuctile_results(p,group_sources) + case (SOURCE_damage_isoBrittle_ID) sourceType + call source_damage_isoBrittle_results(p,group_sources) + case (SOURCE_damage_isoDuctile_ID) sourceType + call source_damage_isoDuctile_results(p,group_sources) + end select sourceType + + enddo SourceLoop + + enddo plasticityLoop + end subroutine constitutive_results diff --git a/src/constitutive_damage.f90 b/src/constitutive_damage.f90 new file mode 100644 index 000000000..cb81c252d --- /dev/null +++ b/src/constitutive_damage.f90 @@ -0,0 +1,163 @@ +submodule(constitutive) constitutive_damage + + implicit none + + interface + + module subroutine source_damage_anisoBrittle_init + end subroutine source_damage_anisoBrittle_init + + module subroutine source_damage_anisoDuctile_init + end subroutine source_damage_anisoDuctile_init + + module subroutine source_damage_isoBrittle_init + end subroutine source_damage_isoBrittle_init + + module subroutine source_damage_isoDuctile_init + end subroutine source_damage_isoDuctile_init + + module subroutine kinematics_cleavage_opening_init + end subroutine kinematics_cleavage_opening_init + + module subroutine kinematics_slipplane_opening_init + end subroutine kinematics_slipplane_opening_init + + module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) + + integer, intent(in) :: & + phase, & + constituent + real(pReal), intent(in) :: & + phi + real(pReal), intent(out) :: & + localphiDot, & + dLocalphiDot_dPhi + + end subroutine source_damage_anisoBrittle_getRateAndItsTangent + + module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) + + integer, intent(in) :: & + phase, & + constituent + real(pReal), intent(in) :: & + phi + real(pReal), intent(out) :: & + localphiDot, & + dLocalphiDot_dPhi + + end subroutine source_damage_anisoDuctile_getRateAndItsTangent + + module subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) + + integer, intent(in) :: & + phase, & + constituent + real(pReal), intent(in) :: & + phi + real(pReal), intent(out) :: & + localphiDot, & + dLocalphiDot_dPhi + + end subroutine source_damage_isoBrittle_getRateAndItsTangent + + module subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) + + integer, intent(in) :: & + phase, & + constituent + real(pReal), intent(in) :: & + phi + real(pReal), intent(out) :: & + localphiDot, & + dLocalphiDot_dPhi + + end subroutine source_damage_isoDuctile_getRateAndItsTangent + + module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase) + + integer, intent(in) :: & + phase + real(pReal), intent(in), dimension(3,3) :: & + Tstar + real(pReal), intent(in), dimension(3,3) :: & + Lp + + real(pReal), intent(out) :: & + TDot, & + dTDot_dT + + end subroutine source_thermal_dissipation_getRateAndItsTangent + end interface + +contains + +module subroutine damage_init + +! initialize source mechanisms + if (any(phase_source == SOURCE_damage_isoBrittle_ID)) call source_damage_isoBrittle_init + if (any(phase_source == SOURCE_damage_isoDuctile_ID)) call source_damage_isoDuctile_init + if (any(phase_source == SOURCE_damage_anisoBrittle_ID)) call source_damage_anisoBrittle_init + if (any(phase_source == SOURCE_damage_anisoDuctile_ID)) call source_damage_anisoDuctile_init + +!-------------------------------------------------------------------------------------------------- +! initialize kinematic mechanisms + if (any(phase_kinematics == KINEMATICS_cleavage_opening_ID)) call kinematics_cleavage_opening_init + if (any(phase_kinematics == KINEMATICS_slipplane_opening_ID)) call kinematics_slipplane_opening_init + +end subroutine damage_init + + +module subroutine damage_source_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) + + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + phi + real(pReal), intent(inout) :: & + phiDot, & + dPhiDot_dPhi + + real(pReal) :: & + localphiDot, & + dLocalphiDot_dPhi + integer :: & + phase, & + grain, & + source, & + constituent + + phiDot = 0.0_pReal + dPhiDot_dPhi = 0.0_pReal + + do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) + phase = material_phaseAt(grain,el) + constituent = material_phasememberAt(grain,ip,el) + do source = 1, phase_Nsources(phase) + select case(phase_source(source,phase)) + case (SOURCE_damage_isoBrittle_ID) + call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) + + case (SOURCE_damage_isoDuctile_ID) + call source_damage_isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) + + case (SOURCE_damage_anisoBrittle_ID) + call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) + + case (SOURCE_damage_anisoDuctile_ID) + call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) + + case default + localphiDot = 0.0_pReal + dLocalphiDot_dPhi = 0.0_pReal + + end select + phiDot = phiDot + localphiDot + dPhiDot_dPhi = dPhiDot_dPhi + dLocalphiDot_dPhi + enddo + enddo + +end subroutine damage_source_getRateAndItsTangents + +end submodule diff --git a/src/constitutive_plastic.f90 b/src/constitutive_plastic.f90 new file mode 100644 index 000000000..d4d0f19a2 --- /dev/null +++ b/src/constitutive_plastic.f90 @@ -0,0 +1,274 @@ +submodule(constitutive) constitutive_plastic + + implicit none + + interface + + module subroutine plastic_none_init + end subroutine plastic_none_init + + module subroutine plastic_isotropic_init + end subroutine plastic_isotropic_init + + module subroutine plastic_phenopowerlaw_init + end subroutine plastic_phenopowerlaw_init + + module subroutine plastic_kinehardening_init + end subroutine plastic_kinehardening_init + + module subroutine plastic_dislotwin_init + end subroutine plastic_dislotwin_init + + module subroutine plastic_disloUCLA_init + end subroutine plastic_disloUCLA_init + + module subroutine plastic_nonlocal_init + end subroutine plastic_nonlocal_init + + + module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) + real(pReal), dimension(3,3), intent(out) :: & + Lp !< plastic velocity gradient + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp !< derivative of Lp with respect to the Mandel stress + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + end subroutine plastic_isotropic_LpAndItsTangent + + pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) + real(pReal), dimension(3,3), intent(out) :: & + Lp !< plastic velocity gradient + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp !< derivative of Lp with respect to the Mandel stress + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + end subroutine plastic_phenopowerlaw_LpAndItsTangent + + pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) + real(pReal), dimension(3,3), intent(out) :: & + Lp !< plastic velocity gradient + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp !< derivative of Lp with respect to the Mandel stress + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + integer, intent(in) :: & + instance, & + of + end subroutine plastic_kinehardening_LpAndItsTangent + + module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) + real(pReal), dimension(3,3), intent(out) :: & + Lp !< plastic velocity gradient + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp !< derivative of Lp with respect to the Mandel stress + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T + integer, intent(in) :: & + instance, & + of + end subroutine plastic_dislotwin_LpAndItsTangent + + pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of) + real(pReal), dimension(3,3), intent(out) :: & + Lp !< plastic velocity gradient + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp !< derivative of Lp with respect to the Mandel stress + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + T + integer, intent(in) :: & + instance, & + of + end subroutine plastic_disloUCLA_LpAndItsTangent + + module subroutine plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp, & + Mp,Temperature,instance,of,ip,el) + real(pReal), dimension(3,3), intent(out) :: & + Lp !< plastic velocity gradient + real(pReal), dimension(3,3,3,3), intent(out) :: & + dLp_dMp !< derivative of Lp with respect to the Mandel stress + + real(pReal), dimension(3,3), intent(in) :: & + Mp !< Mandel stress + real(pReal), intent(in) :: & + Temperature + integer, intent(in) :: & + instance, & + of, & + ip, & !< current integration point + el !< current element number + end subroutine plastic_nonlocal_LpAndItsTangent + + + module subroutine plastic_dislotwin_dependentState(T,instance,of) + integer, intent(in) :: & + instance, & + of + real(pReal), intent(in) :: & + T + end subroutine plastic_dislotwin_dependentState + + module subroutine plastic_disloUCLA_dependentState(instance,of) + integer, intent(in) :: & + instance, & + of + end subroutine plastic_disloUCLA_dependentState + + module subroutine plastic_nonlocal_dependentState(F, Fp, instance, of, ip, el) + real(pReal), dimension(3,3), intent(in) :: & + F, & + Fp + integer, intent(in) :: & + instance, & + of, & + ip, & + el + end subroutine plastic_nonlocal_dependentState + + end interface + + +contains + + +!-------------------------------------------------------------------------------------------------- +!> @brief allocates arrays pointing to array of the various constitutive modules +!-------------------------------------------------------------------------------------------------- +module subroutine plastic_init + +!-------------------------------------------------------------------------------------------------- +! initialized plasticity + if (any(phase_plasticity == PLASTICITY_NONE_ID)) call plastic_none_init + if (any(phase_plasticity == PLASTICITY_ISOTROPIC_ID)) call plastic_isotropic_init + if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_init + if (any(phase_plasticity == PLASTICITY_KINEHARDENING_ID)) call plastic_kinehardening_init + if (any(phase_plasticity == PLASTICITY_DISLOTWIN_ID)) call plastic_dislotwin_init + if (any(phase_plasticity == PLASTICITY_DISLOUCLA_ID)) call plastic_disloucla_init + if (any(phase_plasticity == PLASTICITY_NONLOCAL_ID)) then + call plastic_nonlocal_init + else + call geometry_plastic_nonlocal_disable + endif + +end subroutine plastic_init + + +!-------------------------------------------------------------------------------------------------- +!> @brief calls microstructure function of the different constitutive models +!-------------------------------------------------------------------------------------------------- +module subroutine plastic_dependentState(F, Fp, ipc, ip, el) + + integer, intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + F, & !< elastic deformation gradient + Fp !< plastic deformation gradient + integer :: & + ho, & !< homogenization + tme, & !< thermal member position + instance, of + + ho = material_homogenizationAt(el) + tme = thermalMapping(ho)%p(ip,el) + of = material_phasememberAt(ipc,ip,el) + instance = phase_plasticityInstance(material_phaseAt(ipc,el)) + + plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) + case (PLASTICITY_DISLOTWIN_ID) plasticityType + call plastic_dislotwin_dependentState(temperature(ho)%p(tme),instance,of) + case (PLASTICITY_DISLOUCLA_ID) plasticityType + call plastic_disloUCLA_dependentState(instance,of) + case (PLASTICITY_NONLOCAL_ID) plasticityType + call plastic_nonlocal_dependentState (F,Fp,instance,of,ip,el) + end select plasticityType + +end subroutine plastic_dependentState + +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the velocity gradient +! ToDo: Discuss whether it makes sense if crystallite handles the configuration conversion, i.e. +! Mp in, dLp_dMp out +!-------------------------------------------------------------------------------------------------- +module subroutine plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & + S, 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) :: & + S, & !< 2nd Piola-Kirchhoff stress + Fi !< intermediate deformation gradient + real(pReal), intent(out), dimension(3,3) :: & + Lp !< plastic velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLp_dS, & + dLp_dFi !< derivative of Lp with respect to Fi + real(pReal), dimension(3,3,3,3) :: & + dLp_dMp !< derivative of Lp with respect to Mandel stress + real(pReal), dimension(3,3) :: & + Mp !< Mandel stress work conjugate with Lp + integer :: & + ho, & !< homogenization + tme !< thermal member position + integer :: & + i, j, instance, of + + ho = material_homogenizationAt(el) + tme = thermalMapping(ho)%p(ip,el) + + Mp = matmul(matmul(transpose(Fi),Fi),S) + of = material_phasememberAt(ipc,ip,el) + instance = phase_plasticityInstance(material_phaseAt(ipc,el)) + + plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) + + case (PLASTICITY_NONE_ID) plasticityType + Lp = 0.0_pReal + dLp_dMp = 0.0_pReal + + case (PLASTICITY_ISOTROPIC_ID) plasticityType + call plastic_isotropic_LpAndItsTangent (Lp,dLp_dMp,Mp,instance,of) + + case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType + call plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) + + case (PLASTICITY_KINEHARDENING_ID) plasticityType + call plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of) + + case (PLASTICITY_NONLOCAL_ID) plasticityType + call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMp,Mp, temperature(ho)%p(tme),instance,of,ip,el) + + case (PLASTICITY_DISLOTWIN_ID) plasticityType + call plastic_dislotwin_LpAndItsTangent (Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,of) + + case (PLASTICITY_DISLOUCLA_ID) plasticityType + call plastic_disloucla_LpAndItsTangent (Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,of) + + end select plasticityType + + do i=1,3; do j=1,3 + dLp_dFi(i,j,1:3,1:3) = matmul(matmul(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + & + matmul(matmul(Fi,dLp_dMp(i,j,1:3,1:3)),S) + dLp_dS(i,j,1:3,1:3) = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi) + enddo; enddo + +end subroutine plastic_LpAndItsTangents + +end submodule constitutive_plastic + diff --git a/src/constitutive_plastic_disloUCLA.f90 b/src/constitutive_plastic_disloUCLA.f90 index 8c94817b9..71f8f4a27 100644 --- a/src/constitutive_plastic_disloUCLA.f90 +++ b/src/constitutive_plastic_disloUCLA.f90 @@ -5,7 +5,7 @@ !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief crystal plasticity model for bcc metals, especially Tungsten !-------------------------------------------------------------------------------------------------- -submodule(constitutive) plastic_disloUCLA +submodule(constitutive:constitutive_plastic) plastic_disloUCLA real(pReal), parameter :: & kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin diff --git a/src/constitutive_plastic_dislotwin.f90 b/src/constitutive_plastic_dislotwin.f90 index 09294ecd9..004238817 100644 --- a/src/constitutive_plastic_dislotwin.f90 +++ b/src/constitutive_plastic_dislotwin.f90 @@ -7,7 +7,7 @@ !> @brief material subroutine incoprorating dislocation and twinning physics !> @details to be done !-------------------------------------------------------------------------------------------------- -submodule(constitutive) plastic_dislotwin +submodule(constitutive:constitutive_plastic) plastic_dislotwin real(pReal), parameter :: & kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin diff --git a/src/constitutive_plastic_isotropic.f90 b/src/constitutive_plastic_isotropic.f90 index 60f40fe37..3484b6f64 100644 --- a/src/constitutive_plastic_isotropic.f90 +++ b/src/constitutive_plastic_isotropic.f90 @@ -7,7 +7,7 @@ !! resolving the stress on the slip systems. Will give the response of phenopowerlaw for an !! untextured polycrystal !-------------------------------------------------------------------------------------------------- -submodule(constitutive) plastic_isotropic +submodule(constitutive:constitutive_plastic) plastic_isotropic type :: tParameters real(pReal) :: & diff --git a/src/constitutive_plastic_kinehardening.f90 b/src/constitutive_plastic_kinehardening.f90 index d3e84cfe2..163b5794a 100644 --- a/src/constitutive_plastic_kinehardening.f90 +++ b/src/constitutive_plastic_kinehardening.f90 @@ -5,7 +5,7 @@ !> @brief Phenomenological crystal plasticity using a power law formulation for the shear rates !! and a Voce-type kinematic hardening rule !-------------------------------------------------------------------------------------------------- -submodule(constitutive) plastic_kinehardening +submodule(constitutive:constitutive_plastic) plastic_kinehardening type :: tParameters real(pReal) :: & diff --git a/src/constitutive_plastic_none.f90 b/src/constitutive_plastic_none.f90 index da3ee9796..4e6033499 100644 --- a/src/constitutive_plastic_none.f90 +++ b/src/constitutive_plastic_none.f90 @@ -4,7 +4,7 @@ !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief Dummy plasticity for purely elastic material !-------------------------------------------------------------------------------------------------- -submodule(constitutive) plastic_none +submodule(constitutive:constitutive_plastic) plastic_none contains diff --git a/src/constitutive_plastic_nonlocal.f90 b/src/constitutive_plastic_nonlocal.f90 index ea55024b9..53b173db3 100644 --- a/src/constitutive_plastic_nonlocal.f90 +++ b/src/constitutive_plastic_nonlocal.f90 @@ -4,7 +4,7 @@ !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @brief material subroutine for plasticity including dislocation flux !-------------------------------------------------------------------------------------------------- -submodule(constitutive) plastic_nonlocal +submodule(constitutive:constitutive_plastic) plastic_nonlocal use geometry_plastic_nonlocal, only: & nIPneighbors => geometry_plastic_nonlocal_nIPneighbors, & IPneighborhood => geometry_plastic_nonlocal_IPneighborhood, & diff --git a/src/constitutive_plastic_phenopowerlaw.f90 b/src/constitutive_plastic_phenopowerlaw.f90 index 53e55b319..3f6ade8d2 100644 --- a/src/constitutive_plastic_phenopowerlaw.f90 +++ b/src/constitutive_plastic_phenopowerlaw.f90 @@ -4,7 +4,7 @@ !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief phenomenological crystal plasticity formulation using a powerlaw fitting !-------------------------------------------------------------------------------------------------- -submodule(constitutive) plastic_phenopowerlaw +submodule(constitutive:constitutive_plastic) plastic_phenopowerlaw type :: tParameters real(pReal) :: & diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 new file mode 100644 index 000000000..f1cd12d7f --- /dev/null +++ b/src/constitutive_thermal.f90 @@ -0,0 +1,112 @@ +submodule(constitutive) constitutive_thermal + + implicit none + + interface + + module subroutine source_thermal_dissipation_init + end subroutine source_thermal_dissipation_init + + module subroutine source_thermal_externalheat_init + end subroutine source_thermal_externalheat_init + + module subroutine kinematics_thermal_expansion_init + end subroutine kinematics_thermal_expansion_init + + module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase) + + integer, intent(in) :: & + phase + real(pReal), intent(in), dimension(3,3) :: & + Tstar + real(pReal), intent(in), dimension(3,3) :: & + Lp + + real(pReal), intent(out) :: & + TDot, & + dTDot_dT + + end subroutine source_thermal_dissipation_getRateAndItsTangent + + module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phase, of) + + integer, intent(in) :: & + phase, & + of + real(pReal), intent(out) :: & + TDot, & + dTDot_dT + + end subroutine source_thermal_externalheat_getRateAndItsTangent + + end interface +contains + +module subroutine thermal_init + +! initialize source mechanisms + if (any(phase_source == SOURCE_thermal_dissipation_ID)) call source_thermal_dissipation_init + if (any(phase_source == SOURCE_thermal_externalheat_ID)) call source_thermal_externalheat_init + +!-------------------------------------------------------------------------------------------------- +!initialize kinematic mechanisms + if (any(phase_kinematics == KINEMATICS_thermal_expansion_ID)) call kinematics_thermal_expansion_init + +end subroutine thermal_init + + +module subroutine thermal_source_getRateAndItsTangents(Tdot, dTdot_dT, T, Tstar, Lp, ip, el) + + integer, intent(in) :: & + ip, & !< integration point number + el !< element number + real(pReal), intent(in) :: & + T + real(pReal), intent(in), dimension(:,:,:,:,:) :: & + Tstar, & + Lp + real(pReal), intent(inout) :: & + Tdot, & + dTdot_dT + + real(pReal) :: & + my_Tdot, & + my_dTdot_dT + integer :: & + phase, & + homog, & + instance, & + grain, & + source, & + constituent + + homog = material_homogenizationAt(el) + instance = thermal_typeInstance(homog) + + do grain = 1, homogenization_Ngrains(homog) + phase = material_phaseAt(grain,el) + constituent = material_phasememberAt(grain,ip,el) + do source = 1, phase_Nsources(phase) + select case(phase_source(source,phase)) + case (SOURCE_thermal_dissipation_ID) + call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & + Tstar(1:3,1:3,grain,ip,el), & + Lp(1:3,1:3,grain,ip,el), & + phase) + + case (SOURCE_thermal_externalheat_ID) + call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & + phase, constituent) + + case default + my_Tdot = 0.0_pReal + my_dTdot_dT = 0.0_pReal + end select + Tdot = Tdot + my_Tdot + dTdot_dT = dTdot_dT + my_dTdot_dT + enddo + enddo + +end subroutine thermal_source_getRateAndItsTangents + +end submodule diff --git a/src/damage_local.f90 b/src/damage_local.f90 index 85701acee..66b50064b 100644 --- a/src/damage_local.f90 +++ b/src/damage_local.f90 @@ -9,10 +9,7 @@ module damage_local use config use numerics use YAML_types - use source_damage_isoBrittle - use source_damage_isoDuctile - use source_damage_anisoBrittle - use source_damage_anisoDuctile + use constitutive use results implicit none @@ -128,42 +125,13 @@ subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el el !< element number real(pReal), intent(in) :: & phi - integer :: & - phase, & - grain, & - source, & - constituent real(pReal) :: & - phiDot, dPhiDot_dPhi, localphiDot, dLocalphiDot_dPhi + phiDot, dPhiDot_dPhi phiDot = 0.0_pReal dPhiDot_dPhi = 0.0_pReal - do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) - phase = material_phaseAt(grain,el) - constituent = material_phasememberAt(grain,ip,el) - do source = 1, phase_Nsources(phase) - select case(phase_source(source,phase)) - case (SOURCE_damage_isoBrittle_ID) - call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - - case (SOURCE_damage_isoDuctile_ID) - call source_damage_isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - - case (SOURCE_damage_anisoBrittle_ID) - call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - - case (SOURCE_damage_anisoDuctile_ID) - call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - - case default - localphiDot = 0.0_pReal - dLocalphiDot_dPhi = 0.0_pReal - - end select - phiDot = phiDot + localphiDot - dPhiDot_dPhi = dPhiDot_dPhi + dLocalphiDot_dPhi - enddo - enddo + + call constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) phiDot = phiDot/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) diff --git a/src/damage_nonlocal.f90 b/src/damage_nonlocal.f90 index d3c3c07a2..914133de7 100644 --- a/src/damage_nonlocal.f90 +++ b/src/damage_nonlocal.f90 @@ -10,10 +10,7 @@ module damage_nonlocal use YAML_types use crystallite use lattice - use source_damage_isoBrittle - use source_damage_isoDuctile - use source_damage_anisoBrittle - use source_damage_anisoDuctile + use constitutive use results implicit none @@ -97,43 +94,13 @@ subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el !< element number real(pReal), intent(in) :: & phi - integer :: & - phase, & - grain, & - source, & - constituent real(pReal) :: & - phiDot, dPhiDot_dPhi, localphiDot, dLocalphiDot_dPhi + phiDot, dPhiDot_dPhi phiDot = 0.0_pReal dPhiDot_dPhi = 0.0_pReal - do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) - phase = material_phaseAt(grain,el) - constituent = material_phasememberAt(grain,ip,el) - do source = 1, phase_Nsources(phase) - select case(phase_source(source,phase)) - case (SOURCE_damage_isoBrittle_ID) - call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - - case (SOURCE_damage_isoDuctile_ID) - call source_damage_isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - - case (SOURCE_damage_anisoBrittle_ID) - call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - - case (SOURCE_damage_anisoDuctile_ID) - call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - - case default - localphiDot = 0.0_pReal - dLocalphiDot_dPhi = 0.0_pReal - - end select - phiDot = phiDot + localphiDot - dPhiDot_dPhi = dPhiDot_dPhi + dLocalphiDot_dPhi - enddo - enddo - + + call constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) phiDot = phiDot/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) diff --git a/src/kinematics_cleavage_opening.f90 b/src/kinematics_cleavage_opening.f90 index e35f37e0e..2e314024c 100644 --- a/src/kinematics_cleavage_opening.f90 +++ b/src/kinematics_cleavage_opening.f90 @@ -4,16 +4,7 @@ !> @brief material subroutine incorporating kinematics resulting from opening of cleavage planes !> @details to be done !-------------------------------------------------------------------------------------------------- -module kinematics_cleavage_opening - use prec - use IO - use config - use math - use lattice - use material - - implicit none - private +submodule(constitutive:constitutive_damage) kinematics_cleavage_opening integer, dimension(:), allocatable :: kinematics_cleavage_opening_instance @@ -31,9 +22,6 @@ module kinematics_cleavage_opening type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) - public :: & - kinematics_cleavage_opening_init, & - kinematics_cleavage_opening_LiAndItsTangent contains @@ -42,7 +30,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine kinematics_cleavage_opening_init +module subroutine kinematics_cleavage_opening_init integer :: Ninstance,p integer, dimension(:), allocatable :: N_cl !< active number of cleavage systems per family @@ -95,7 +83,7 @@ end subroutine kinematics_cleavage_opening_init !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) +module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) integer, intent(in) :: & ipc, & !< grain number @@ -158,4 +146,4 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, i end subroutine kinematics_cleavage_opening_LiAndItsTangent -end module kinematics_cleavage_opening +end submodule kinematics_cleavage_opening diff --git a/src/kinematics_slipplane_opening.f90 b/src/kinematics_slipplane_opening.f90 index 847dc6c72..b15d206d3 100644 --- a/src/kinematics_slipplane_opening.f90 +++ b/src/kinematics_slipplane_opening.f90 @@ -4,16 +4,7 @@ !> @brief material subroutine incorporating kinematics resulting from opening of slip planes !> @details to be done !-------------------------------------------------------------------------------------------------- -module kinematics_slipplane_opening - use prec - use config - use IO - use math - use lattice - use material - - implicit none - private +submodule(constitutive:constitutive_damage) kinematics_slipplane_opening integer, dimension(:), allocatable :: kinematics_slipplane_opening_instance @@ -33,9 +24,6 @@ module kinematics_slipplane_opening type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) - public :: & - kinematics_slipplane_opening_init, & - kinematics_slipplane_opening_LiAndItsTangent contains @@ -44,7 +32,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine kinematics_slipplane_opening_init +module subroutine kinematics_slipplane_opening_init integer :: Ninstance,p,i character(len=pStringLen) :: extmsg = '' @@ -107,7 +95,7 @@ end subroutine kinematics_slipplane_opening_init !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) +module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) integer, intent(in) :: & ipc, & !< grain number @@ -183,4 +171,4 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, end subroutine kinematics_slipplane_opening_LiAndItsTangent -end module kinematics_slipplane_opening +end submodule kinematics_slipplane_opening diff --git a/src/kinematics_thermal_expansion.f90 b/src/kinematics_thermal_expansion.f90 index 63da6eb51..00871e94f 100644 --- a/src/kinematics_thermal_expansion.f90 +++ b/src/kinematics_thermal_expansion.f90 @@ -3,16 +3,7 @@ !> @brief material subroutine incorporating kinematics resulting from thermal expansion !> @details to be done !-------------------------------------------------------------------------------------------------- -module kinematics_thermal_expansion - use prec - use IO - use config - use math - use lattice - use material - - implicit none - private +submodule(constitutive:constitutive_thermal) kinematics_thermal_expansion integer, dimension(:), allocatable :: kinematics_thermal_expansion_instance @@ -25,10 +16,6 @@ module kinematics_thermal_expansion type(tParameters), dimension(:), allocatable :: param - public :: & - kinematics_thermal_expansion_init, & - kinematics_thermal_expansion_initialStrain, & - kinematics_thermal_expansion_LiAndItsTangent contains @@ -37,7 +24,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine kinematics_thermal_expansion_init +module subroutine kinematics_thermal_expansion_init integer :: Ninstance,p,i real(pReal), dimension(:), allocatable :: temp @@ -79,30 +66,30 @@ end subroutine kinematics_thermal_expansion_init !-------------------------------------------------------------------------------------------------- !> @brief report initial thermal strain based on current temperature deviation from reference !-------------------------------------------------------------------------------------------------- -pure function kinematics_thermal_expansion_initialStrain(homog,phase,offset) - - integer, intent(in) :: & - phase, & - homog, & - offset - - real(pReal), dimension(3,3) :: & - kinematics_thermal_expansion_initialStrain !< initial thermal strain (should be small strain, though) - - associate(prm => param(kinematics_thermal_expansion_instance(phase))) - kinematics_thermal_expansion_initialStrain = & - (temperature(homog)%p(offset) - prm%T_ref)**1 / 1. * prm%expansion(1:3,1:3,1) + & ! constant coefficient - (temperature(homog)%p(offset) - prm%T_ref)**2 / 2. * prm%expansion(1:3,1:3,2) + & ! linear coefficient - (temperature(homog)%p(offset) - prm%T_ref)**3 / 3. * prm%expansion(1:3,1:3,3) ! quadratic coefficient - end associate - -end function kinematics_thermal_expansion_initialStrain +!pure module function kinematics_thermal_expansion_initialStrain(homog,phase,offset) +! +! integer, intent(in) :: & +! phase, & +! homog, & +! offset +! +! real(pReal), dimension(3,3) :: & +! kinematics_thermal_expansion_initialStrain !< initial thermal strain (should be small strain, though) +! +! associate(prm => param(kinematics_thermal_expansion_instance(phase))) +! kinematics_thermal_expansion_initialStrain = & +! (temperature(homog)%p(offset) - prm%T_ref)**1 / 1. * prm%expansion(1:3,1:3,1) + & ! constant coefficient +! (temperature(homog)%p(offset) - prm%T_ref)**2 / 2. * prm%expansion(1:3,1:3,2) + & ! linear coefficient +! (temperature(homog)%p(offset) - prm%T_ref)**3 / 3. * prm%expansion(1:3,1:3,3) ! quadratic coefficient +! end associate +! +!end function kinematics_thermal_expansion_initialStrain !-------------------------------------------------------------------------------------------------- !> @brief constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, el) +module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, el) integer, intent(in) :: & ipc, & !< grain number @@ -140,4 +127,4 @@ subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, end subroutine kinematics_thermal_expansion_LiAndItsTangent -end module kinematics_thermal_expansion +end submodule kinematics_thermal_expansion diff --git a/src/source_damage_anisoBrittle.f90 b/src/source_damage_anisoBrittle.f90 index c5ff1f906..20cf3a914 100644 --- a/src/source_damage_anisoBrittle.f90 +++ b/src/source_damage_anisoBrittle.f90 @@ -4,18 +4,7 @@ !> @brief material subroutine incorporating anisotropic brittle damage source mechanism !> @details to be done !-------------------------------------------------------------------------------------------------- -module source_damage_anisoBrittle - use prec - use IO - use math - use discretization - use material - use config - use lattice - use results - - implicit none - private +submodule (constitutive:constitutive_damage) source_damage_anisoBrittle integer, dimension(:), allocatable :: & source_damage_anisoBrittle_offset, & !< which source is my current source mechanism? @@ -39,12 +28,6 @@ module source_damage_anisoBrittle type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) - public :: & - source_damage_anisoBrittle_init, & - source_damage_anisoBrittle_dotState, & - source_damage_anisobrittle_getRateAndItsTangent, & - source_damage_anisoBrittle_results - contains @@ -52,7 +35,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine source_damage_anisoBrittle_init +module subroutine source_damage_anisoBrittle_init integer :: Ninstance,sourceOffset,NipcMyPhase,p integer, dimension(:), allocatable :: N_cl @@ -123,7 +106,7 @@ end subroutine source_damage_anisoBrittle_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el) +module subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el) integer, intent(in) :: & ipc, & !< component-ID of integration point @@ -172,7 +155,7 @@ end subroutine source_damage_anisoBrittle_dotState !-------------------------------------------------------------------------------------------------- !> @brief returns local part of nonlocal damage driving force !-------------------------------------------------------------------------------------------------- -subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) +module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) integer, intent(in) :: & phase, & @@ -199,7 +182,7 @@ end subroutine source_damage_anisoBrittle_getRateAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- -subroutine source_damage_anisoBrittle_results(phase,group) +module subroutine source_damage_anisoBrittle_results(phase,group) integer, intent(in) :: phase character(len=*), intent(in) :: group @@ -218,4 +201,4 @@ subroutine source_damage_anisoBrittle_results(phase,group) end subroutine source_damage_anisoBrittle_results -end module source_damage_anisoBrittle +end submodule source_damage_anisoBrittle diff --git a/src/source_damage_anisoDuctile.f90 b/src/source_damage_anisoDuctile.f90 index e8a76dc3a..3722f8b1c 100644 --- a/src/source_damage_anisoDuctile.f90 +++ b/src/source_damage_anisoDuctile.f90 @@ -4,23 +4,13 @@ !> @brief material subroutine incorporating anisotropic ductile damage source mechanism !> @details to be done !-------------------------------------------------------------------------------------------------- -module source_damage_anisoDuctile - use prec - use IO - use math - use discretization - use material - use config - use results - - implicit none - private +submodule(constitutive:constitutive_damage) source_damage_anisoDuctile integer, dimension(:), allocatable :: & source_damage_anisoDuctile_offset, & !< which source is my current damage mechanism? source_damage_anisoDuctile_instance !< instance of damage source mechanism - type, private :: tParameters !< container type for internal constitutive parameters + type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & n real(pReal), dimension(:), allocatable :: & @@ -29,14 +19,7 @@ module source_damage_anisoDuctile output end type tParameters - type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) - - - public :: & - source_damage_anisoDuctile_init, & - source_damage_anisoDuctile_dotState, & - source_damage_anisoDuctile_getRateAndItsTangent, & - source_damage_anisoDuctile_results + type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) contains @@ -45,7 +28,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine source_damage_anisoDuctile_init +module subroutine source_damage_anisoDuctile_init integer :: Ninstance,sourceOffset,NipcMyPhase,p integer, dimension(:), allocatable :: N_sl @@ -105,7 +88,7 @@ end subroutine source_damage_anisoDuctile_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -subroutine source_damage_anisoDuctile_dotState(ipc, ip, el) +module subroutine source_damage_anisoDuctile_dotState(ipc, ip, el) integer, intent(in) :: & ipc, & !< component-ID of integration point @@ -136,7 +119,7 @@ end subroutine source_damage_anisoDuctile_dotState !-------------------------------------------------------------------------------------------------- !> @brief returns local part of nonlocal damage driving force !-------------------------------------------------------------------------------------------------- -subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) +module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) integer, intent(in) :: & phase, & @@ -163,7 +146,7 @@ end subroutine source_damage_anisoDuctile_getRateAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- -subroutine source_damage_anisoDuctile_results(phase,group) +module subroutine source_damage_anisoDuctile_results(phase,group) integer, intent(in) :: phase character(len=*), intent(in) :: group @@ -182,4 +165,4 @@ subroutine source_damage_anisoDuctile_results(phase,group) end subroutine source_damage_anisoDuctile_results -end module source_damage_anisoDuctile +end submodule source_damage_anisoDuctile diff --git a/src/source_damage_isoBrittle.f90 b/src/source_damage_isoBrittle.f90 index c4c4c72a4..00704fe26 100644 --- a/src/source_damage_isoBrittle.f90 +++ b/src/source_damage_isoBrittle.f90 @@ -4,23 +4,13 @@ !> @brief material subroutine incoprorating isotropic brittle damage source mechanism !> @details to be done !-------------------------------------------------------------------------------------------------- -module source_damage_isoBrittle - use prec - use IO - use math - use discretization - use material - use config - use results - - implicit none - private +submodule(constitutive:constitutive_damage) source_damage_isoBrittle integer, dimension(:), allocatable :: & source_damage_isoBrittle_offset, & source_damage_isoBrittle_instance - type, private :: tParameters !< container type for internal constitutive parameters + type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & critStrainEnergy, & N @@ -30,13 +20,6 @@ module source_damage_isoBrittle type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) - - public :: & - source_damage_isoBrittle_init, & - source_damage_isoBrittle_deltaState, & - source_damage_isoBrittle_getRateAndItsTangent, & - source_damage_isoBrittle_results - contains @@ -44,7 +27,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine source_damage_isoBrittle_init +module subroutine source_damage_isoBrittle_init integer :: Ninstance,sourceOffset,NipcMyPhase,p character(len=pStringLen) :: extmsg = '' @@ -99,7 +82,7 @@ end subroutine source_damage_isoBrittle_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el) +module subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el) integer, intent(in) :: & ipc, & !< component-ID of integration point @@ -145,7 +128,7 @@ end subroutine source_damage_isoBrittle_deltaState !-------------------------------------------------------------------------------------------------- !> @brief returns local part of nonlocal damage driving force !-------------------------------------------------------------------------------------------------- -subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) +module subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) integer, intent(in) :: & phase, & @@ -174,7 +157,7 @@ end subroutine source_damage_isoBrittle_getRateAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- -subroutine source_damage_isoBrittle_results(phase,group) +module subroutine source_damage_isoBrittle_results(phase,group) integer, intent(in) :: phase character(len=*), intent(in) :: group @@ -193,4 +176,4 @@ subroutine source_damage_isoBrittle_results(phase,group) end subroutine source_damage_isoBrittle_results -end module source_damage_isoBrittle +end submodule source_damage_isoBrittle diff --git a/src/source_damage_isoDuctile.f90 b/src/source_damage_isoDuctile.f90 index cf392fdc4..517332316 100644 --- a/src/source_damage_isoDuctile.f90 +++ b/src/source_damage_isoDuctile.f90 @@ -4,22 +4,13 @@ !> @brief material subroutine incorporating isotropic ductile damage source mechanism !> @details to be done !-------------------------------------------------------------------------------------------------- -module source_damage_isoDuctile - use prec - use IO - use discretization - use material - use config - use results - - implicit none - private +submodule (constitutive:constitutive_damage) source_damage_isoDuctile integer, dimension(:), allocatable :: & source_damage_isoDuctile_offset, & !< which source is my current damage mechanism? source_damage_isoDuctile_instance !< instance of damage source mechanism - type, private :: tParameters !< container type for internal constitutive parameters + type:: tParameters !< container type for internal constitutive parameters real(pReal) :: & critPlasticStrain, & N @@ -27,15 +18,9 @@ module source_damage_isoDuctile output end type tParameters - type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) + type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) - public :: & - source_damage_isoDuctile_init, & - source_damage_isoDuctile_dotState, & - source_damage_isoDuctile_getRateAndItsTangent, & - source_damage_isoDuctile_Results - contains @@ -43,7 +28,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine source_damage_isoDuctile_init +module subroutine source_damage_isoDuctile_init integer :: Ninstance,sourceOffset,NipcMyPhase,p character(len=pStringLen) :: extmsg = '' @@ -98,7 +83,7 @@ end subroutine source_damage_isoDuctile_init !-------------------------------------------------------------------------------------------------- !> @brief calculates derived quantities from state !-------------------------------------------------------------------------------------------------- -subroutine source_damage_isoDuctile_dotState(ipc, ip, el) +module subroutine source_damage_isoDuctile_dotState(ipc, ip, el) integer, intent(in) :: & ipc, & !< component-ID of integration point @@ -129,7 +114,7 @@ end subroutine source_damage_isoDuctile_dotState !-------------------------------------------------------------------------------------------------- !> @brief returns local part of nonlocal damage driving force !-------------------------------------------------------------------------------------------------- -subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) +module subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) integer, intent(in) :: & phase, & @@ -156,7 +141,7 @@ end subroutine source_damage_isoDuctile_getRateAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- -subroutine source_damage_isoDuctile_results(phase,group) +module subroutine source_damage_isoDuctile_results(phase,group) integer, intent(in) :: phase character(len=*), intent(in) :: group @@ -175,4 +160,4 @@ subroutine source_damage_isoDuctile_results(phase,group) end subroutine source_damage_isoDuctile_results -end module source_damage_isoDuctile +end submodule source_damage_isoDuctile diff --git a/src/source_thermal_dissipation.f90 b/src/source_thermal_dissipation.f90 index 8fa5ae0c7..58a0c6b3c 100644 --- a/src/source_thermal_dissipation.f90 +++ b/src/source_thermal_dissipation.f90 @@ -4,14 +4,7 @@ !> @brief material subroutine for thermal source due to plastic dissipation !> @details to be done !-------------------------------------------------------------------------------------------------- -module source_thermal_dissipation - use prec - use discretization - use material - use config - - implicit none - private +submodule(constitutive:constitutive_thermal) source_thermal_dissipation integer, dimension(:), allocatable :: & source_thermal_dissipation_offset, & !< which source is my current thermal dissipation mechanism? @@ -25,10 +18,6 @@ module source_thermal_dissipation type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) - public :: & - source_thermal_dissipation_init, & - source_thermal_dissipation_getRateAndItsTangent - contains @@ -36,7 +25,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine source_thermal_dissipation_init +module subroutine source_thermal_dissipation_init integer :: Ninstance,sourceOffset,NipcMyPhase,p @@ -76,7 +65,7 @@ end subroutine source_thermal_dissipation_init !-------------------------------------------------------------------------------------------------- !> @brief Ninstances dissipation rate !-------------------------------------------------------------------------------------------------- -subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase) +module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase) integer, intent(in) :: & phase @@ -96,4 +85,4 @@ subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar end subroutine source_thermal_dissipation_getRateAndItsTangent -end module source_thermal_dissipation +end submodule source_thermal_dissipation diff --git a/src/source_thermal_externalheat.f90 b/src/source_thermal_externalheat.f90 index b83b29596..8ffc7a4fb 100644 --- a/src/source_thermal_externalheat.f90 +++ b/src/source_thermal_externalheat.f90 @@ -4,14 +4,8 @@ !> @author Philip Eisenlohr, Michigan State University !> @brief material subroutine for variable heat source !-------------------------------------------------------------------------------------------------- -module source_thermal_externalheat - use prec - use discretization - use material - use config +submodule(constitutive:constitutive_thermal) source_thermal_externalheat - implicit none - private integer, dimension(:), allocatable :: & source_thermal_externalheat_offset, & !< which source is my current thermal dissipation mechanism? @@ -28,11 +22,6 @@ module source_thermal_externalheat type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) - public :: & - source_thermal_externalheat_init, & - source_thermal_externalheat_dotState, & - source_thermal_externalheat_getRateAndItsTangent - contains @@ -40,7 +29,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine source_thermal_externalheat_init +module subroutine source_thermal_externalheat_init integer :: Ninstance,sourceOffset,NipcMyPhase,p @@ -84,7 +73,7 @@ end subroutine source_thermal_externalheat_init !> @brief rate of change of state !> @details state only contains current time to linearly interpolate given heat powers !-------------------------------------------------------------------------------------------------- -subroutine source_thermal_externalheat_dotState(phase, of) +module subroutine source_thermal_externalheat_dotState(phase, of) integer, intent(in) :: & phase, & @@ -103,7 +92,7 @@ end subroutine source_thermal_externalheat_dotState !-------------------------------------------------------------------------------------------------- !> @brief returns local heat generation rate !-------------------------------------------------------------------------------------------------- -subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phase, of) +module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phase, of) integer, intent(in) :: & phase, & @@ -135,4 +124,4 @@ subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phas end subroutine source_thermal_externalheat_getRateAndItsTangent -end module source_thermal_externalheat +end submodule source_thermal_externalheat diff --git a/src/thermal_adiabatic.f90 b/src/thermal_adiabatic.f90 index a11e7b0a1..c52f0a3d0 100644 --- a/src/thermal_adiabatic.f90 +++ b/src/thermal_adiabatic.f90 @@ -8,8 +8,7 @@ module thermal_adiabatic use numerics use material use results - use source_thermal_dissipation - use source_thermal_externalheat + use constitutive use crystallite use lattice @@ -125,46 +124,15 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) T real(pReal), intent(out) :: & Tdot, dTdot_dT - - real(pReal) :: & - my_Tdot, my_dTdot_dT integer :: & - phase, & - homog, & - instance, & - grain, & - source, & - constituent - - homog = material_homogenizationAt(el) - instance = thermal_typeInstance(homog) - + homog + Tdot = 0.0_pReal dTdot_dT = 0.0_pReal - do grain = 1, homogenization_Ngrains(homog) - phase = material_phaseAt(grain,el) - constituent = material_phasememberAt(grain,ip,el) - do source = 1, phase_Nsources(phase) - select case(phase_source(source,phase)) - case (SOURCE_thermal_dissipation_ID) - call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & - crystallite_S(1:3,1:3,grain,ip,el), & - crystallite_Lp(1:3,1:3,grain,ip,el), & - phase) - case (SOURCE_thermal_externalheat_ID) - call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & - phase, constituent) - - case default - my_Tdot = 0.0_pReal - my_dTdot_dT = 0.0_pReal - end select - Tdot = Tdot + my_Tdot - dTdot_dT = dTdot_dT + my_dTdot_dT - enddo - enddo - + homog = material_homogenizationAt(el) + call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, crystallite_S, crystallite_Lp, ip, el) + Tdot = Tdot/real(homogenization_Ngrains(homog),pReal) dTdot_dT = dTdot_dT/real(homogenization_Ngrains(homog),pReal) diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index 82eed4cc4..60766710d 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -9,8 +9,7 @@ module thermal_conduction use lattice use results use crystallite - use source_thermal_dissipation - use source_thermal_externalheat + use constitutive implicit none private @@ -84,46 +83,14 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) T real(pReal), intent(out) :: & Tdot, dTdot_dT - real(pReal) :: & - my_Tdot, my_dTdot_dT integer :: & - phase, & - homog, & - offset, & - instance, & - grain, & - source, & - constituent - - homog = material_homogenizationAt(el) - offset = material_homogenizationMemberAt(ip,el) - instance = thermal_typeInstance(homog) - + homog + Tdot = 0.0_pReal dTdot_dT = 0.0_pReal - do grain = 1, homogenization_Ngrains(homog) - phase = material_phaseAt(grain,el) - constituent = material_phasememberAt(grain,ip,el) - do source = 1, phase_Nsources(phase) - select case(phase_source(source,phase)) - case (SOURCE_thermal_dissipation_ID) - call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & - crystallite_S(1:3,1:3,grain,ip,el), & - crystallite_Lp(1:3,1:3,grain,ip,el), & - phase) - case (SOURCE_thermal_externalheat_ID) - call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & - phase, constituent) - case default - my_Tdot = 0.0_pReal - my_dTdot_dT = 0.0_pReal - - end select - Tdot = Tdot + my_Tdot - dTdot_dT = dTdot_dT + my_dTdot_dT - enddo - enddo + homog = material_homogenizationAt(el) + call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, crystallite_S,crystallite_Lp ,ip, el) Tdot = Tdot/real(homogenization_Ngrains(homog),pReal) dTdot_dT = dTdot_dT/real(homogenization_Ngrains(homog),pReal)