diff --git a/PRIVATE b/PRIVATE index 591964dcf..7846c7112 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 591964dcf8521d95f6cccbfe840d462c430e63d9 +Subproject commit 7846c71126705cc5d41dd79f2d595f4864434068 diff --git a/src/DAMASK_marc.f90 b/src/DAMASK_marc.f90 index 0ad68445c..0f9d37ddb 100644 --- a/src/DAMASK_marc.f90 +++ b/src/DAMASK_marc.f90 @@ -364,7 +364,8 @@ subroutine flux(f,ts,n,time) real(pReal), dimension(2), intent(out) :: & f - call thermal_conduction_getSourceAndItsTangent(f(1), f(2), ts(3), n(3),mesh_FEM2DAMASK_elem(n(1))) + f(2) = 0.0_pReal + call thermal_conduction_getSource(f(1), ts(3), n(3),mesh_FEM2DAMASK_elem(n(1))) end subroutine flux diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index 674ec4c43..e04936a3b 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -32,8 +32,8 @@ #include "constitutive_plastic_disloTungsten.f90" #include "constitutive_plastic_nonlocal.f90" #include "constitutive_thermal.f90" -#include "source_thermal_dissipation.f90" -#include "source_thermal_externalheat.f90" +#include "constitutive_thermal_dissipation.f90" +#include "constitutive_thermal_externalheat.f90" #include "kinematics_thermal_expansion.f90" #include "constitutive_damage.f90" #include "source_damage_isoBrittle.f90" @@ -51,4 +51,5 @@ #include "homogenization_mech_none.f90" #include "homogenization_mech_isostrain.f90" #include "homogenization_mech_RGC.f90" +#include "homogenization_thermal.f90" #include "CPFEM.f90" diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 7dac3e129..c6415a883 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -21,21 +21,6 @@ module constitutive private enum, bind(c); enumerator :: & - PLASTICITY_UNDEFINED_ID, & - PLASTICITY_NONE_ID, & - PLASTICITY_ISOTROPIC_ID, & - PLASTICITY_PHENOPOWERLAW_ID, & - PLASTICITY_KINEHARDENING_ID, & - PLASTICITY_DISLOTWIN_ID, & - PLASTICITY_DISLOTUNGSTEN_ID, & - PLASTICITY_NONLOCAL_ID, & - SOURCE_UNDEFINED_ID ,& - SOURCE_THERMAL_DISSIPATION_ID, & - SOURCE_THERMAL_EXTERNALHEAT_ID, & - SOURCE_DAMAGE_ISOBRITTLE_ID, & - SOURCE_DAMAGE_ISODUCTILE_ID, & - SOURCE_DAMAGE_ANISOBRITTLE_ID, & - SOURCE_DAMAGE_ANISODUCTILE_ID, & KINEMATICS_UNDEFINED_ID ,& KINEMATICS_CLEAVAGE_OPENING_ID, & KINEMATICS_SLIPPLANE_OPENING_ID, & @@ -81,15 +66,11 @@ module constitutive type(tDebugOptions) :: debugCrystallite - - integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable, public :: & - phase_plasticity !< plasticity of each phase - - integer(kind(SOURCE_undefined_ID)), dimension(:,:), allocatable :: & - phase_source, & !< active sources mechanisms of each phase + integer(kind(KINEMATICS_UNDEFINED_ID)), dimension(:,:), allocatable :: & phase_kinematics !< active kinematic mechanisms of each phase integer, dimension(:), allocatable, public :: & !< ToDo: should be protected (bug in Intel compiler) + thermal_Nsources, & phase_Nsources, & !< number of source mechanisms active in each phase phase_Nkinematics, & !< number of kinematic mechanisms active in each phase phase_NstiffnessDegradations, & !< number of stiffness degradation mechanisms active in each phase @@ -102,7 +83,7 @@ module constitutive type(tPlasticState), allocatable, dimension(:), public :: & plasticState type(tSourceState), allocatable, dimension(:), public :: & - sourceState + damageState, thermalState integer, public, protected :: & @@ -112,13 +93,15 @@ module constitutive interface ! == cleaned:begin ================================================================================= - module subroutine mech_init + module subroutine mech_init(phases) + class(tNode), pointer :: phases end subroutine mech_init module subroutine damage_init end subroutine damage_init - module subroutine thermal_init + module subroutine thermal_init(phases) + class(tNode), pointer :: phases end subroutine thermal_init @@ -137,21 +120,37 @@ module constitutive integer, intent(in) :: ph, me end subroutine mech_initializeRestorationPoints - module subroutine constitutive_mech_windForward(ph,me) + module subroutine thermal_initializeRestorationPoints(ph,me) integer, intent(in) :: ph, me - end subroutine constitutive_mech_windForward + end subroutine thermal_initializeRestorationPoints + + + module subroutine mech_windForward(ph,me) + integer, intent(in) :: ph, me + end subroutine mech_windForward + + module subroutine thermal_windForward(ph,me) + integer, intent(in) :: ph, me + end subroutine thermal_windForward + + + module subroutine mech_forward() + end subroutine mech_forward + + module subroutine thermal_forward() + end subroutine thermal_forward - module subroutine constitutive_mech_forward - end subroutine constitutive_mech_forward module subroutine mech_restore(ip,el,includeL) - integer, intent(in) :: & - ip, & - el - logical, intent(in) :: & - includeL + integer, intent(in) :: ip, el + logical, intent(in) :: includeL end subroutine mech_restore + module subroutine thermal_restore(ip,el) + integer, intent(in) :: ip, el + end subroutine thermal_restore + + module function constitutive_mech_dPdF(dt,co,ip,el) result(dPdF) real(pReal), intent(in) :: dt integer, intent(in) :: & @@ -172,43 +171,67 @@ module constitutive end subroutine mech_restartRead - module function constitutive_mech_getS(co,ip,el) result(S) - integer, intent(in) :: co, ip, el + module function mech_S(ph,me) result(S) + integer, intent(in) :: ph,me real(pReal), dimension(3,3) :: S - end function constitutive_mech_getS + end function mech_S - module function constitutive_mech_getLp(co,ip,el) result(Lp) - integer, intent(in) :: co, ip, el - real(pReal), dimension(3,3) :: Lp - end function constitutive_mech_getLp + module function mech_L_p(ph,me) result(L_p) + integer, intent(in) :: ph,me + real(pReal), dimension(3,3) :: L_p + end function mech_L_p module function constitutive_mech_getF(co,ip,el) result(F) integer, intent(in) :: co, ip, el real(pReal), dimension(3,3) :: F end function constitutive_mech_getF - module function constitutive_mech_getF_e(co,ip,el) result(F_e) - integer, intent(in) :: co, ip, el + module function mech_F_e(ph,me) result(F_e) + integer, intent(in) :: ph,me real(pReal), dimension(3,3) :: F_e - end function constitutive_mech_getF_e + end function mech_F_e module function constitutive_mech_getP(co,ip,el) result(P) integer, intent(in) :: co, ip, el real(pReal), dimension(3,3) :: P end function constitutive_mech_getP - module function constitutive_thermal_T(co,ip,el) result(T) - integer, intent(in) :: co, ip, el + module function thermal_T(ph,me) result(T) + integer, intent(in) :: ph,me real(pReal) :: T - end function constitutive_thermal_T + end function thermal_T + module subroutine constitutive_mech_setF(F,co,ip,el) real(pReal), dimension(3,3), intent(in) :: F integer, intent(in) :: co, ip, el end subroutine constitutive_mech_setF + module subroutine constitutive_thermal_setT(T,co,ip,el) + real(pReal), intent(in) :: T + integer, intent(in) :: co, ip, el + end subroutine constitutive_thermal_setT + ! == cleaned:end =================================================================================== + module function integrateThermalState(Delta_t,co,ip,el) result(broken) + real(pReal), intent(in) :: Delta_t + integer, intent(in) :: & + el, & !< element index in element loop + ip, & !< integration point index in ip loop + co !< grain index in grain loop + logical :: broken + end function integrateThermalState + + module function integrateDamageState(dt,co,ip,el) result(broken) + real(pReal), intent(in) :: dt + integer, intent(in) :: & + el, & !< element index in element loop + ip, & !< integration point index in ip loop + co !< grain index in grain loop + logical :: broken + end function integrateDamageState + module function crystallite_stress(dt,co,ip,el) result(converged_) real(pReal), intent(in) :: dt integer, intent(in) :: co, ip, el @@ -260,16 +283,15 @@ module constitutive dPhiDot_dPhi end subroutine constitutive_damage_getRateAndItsTangents - module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, ip, el) + module subroutine constitutive_thermal_getRate(TDot, T,ip,el) integer, intent(in) :: & ip, & !< integration point number el !< element number real(pReal), intent(in) :: & T - real(pReal), intent(inout) :: & - TDot, & - dTDot_dT - end subroutine constitutive_thermal_getRateAndItsTangents + real(pReal), intent(out) :: & + TDot + end subroutine constitutive_thermal_getRate @@ -282,6 +304,7 @@ module constitutive orientation !< crystal orientation end subroutine plastic_nonlocal_updateCompatibility + module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of) real(pReal), dimension(3,3), intent(out) :: & Li !< inleastic velocity gradient @@ -344,47 +367,13 @@ module constitutive end subroutine source_damage_isoBrittle_deltaState - module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & - S, Fi, co, ip, el) - integer, intent(in) :: & - co, & !< 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 constitutive_plastic_LpAndItsTangents - - module subroutine constitutive_plastic_dependentState(co,ip,el) integer, intent(in) :: & - co, & !< component-ID of integration point + co, & !< component-ID of integration point ip, & !< integration point el !< element end subroutine constitutive_plastic_dependentState - - - module subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, co, ip, el) - integer, intent(in) :: & - co, & !< 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_hooke_SandItsTangents - end interface @@ -394,15 +383,13 @@ module constitutive public :: & constitutive_init, & constitutive_homogenizedC, & - constitutive_LiAndItsTangents, & constitutive_damage_getRateAndItsTangents, & - constitutive_thermal_getRateAndItsTangents, & + constitutive_thermal_getRate, & constitutive_results, & constitutive_allocateState, & constitutive_forward, & constitutive_restore, & plastic_nonlocal_updateCompatibility, & - source_active, & kinematics_active, & converged, & crystallite_init, & @@ -412,29 +399,14 @@ module constitutive crystallite_push33ToRef, & constitutive_restartWrite, & constitutive_restartRead, & - integrateSourceState, & - constitutive_mech_setF, & + integrateThermalState, & + integrateDamageState, & + constitutive_thermal_setT, & constitutive_mech_getP, & - constitutive_mech_getLp, & + constitutive_mech_setF, & constitutive_mech_getF, & - constitutive_mech_getS, & constitutive_initializeRestorationPoints, & constitutive_windForward, & - PLASTICITY_UNDEFINED_ID, & - PLASTICITY_NONE_ID, & - PLASTICITY_ISOTROPIC_ID, & - PLASTICITY_PHENOPOWERLAW_ID, & - PLASTICITY_KINEHARDENING_ID, & - PLASTICITY_DISLOTWIN_ID, & - PLASTICITY_DISLOTUNGSTEN_ID, & - PLASTICITY_NONLOCAL_ID, & - SOURCE_UNDEFINED_ID ,& - SOURCE_THERMAL_DISSIPATION_ID, & - SOURCE_THERMAL_EXTERNALHEAT_ID, & - SOURCE_DAMAGE_ISOBRITTLE_ID, & - SOURCE_DAMAGE_ISODUCTILE_ID, & - SOURCE_DAMAGE_ANISOBRITTLE_ID, & - SOURCE_DAMAGE_ANISODUCTILE_ID, & KINEMATICS_UNDEFINED_ID ,& KINEMATICS_CLEAVAGE_OPENING_ID, & KINEMATICS_SLIPPLANE_OPENING_ID, & @@ -456,6 +428,8 @@ subroutine constitutive_init phases + print'(/,a)', ' <<<+- constitutive init -+>>>'; flush(IO_STDOUT) + debug_constitutive => config_debug%get('constitutive', defaultVal=emptyList) debugConstitutive%basic = debug_constitutive%contains('basic') debugConstitutive%extensive = debug_constitutive%contains('extensive') @@ -464,15 +438,14 @@ subroutine constitutive_init debugConstitutive%ip = config_debug%get_asInt('integrationpoint',defaultVal = 1) debugConstitutive%grain = config_debug%get_asInt('grain',defaultVal = 1) -!-------------------------------------------------------------------------------------------------- -! initialize constitutive laws - print'(/,a)', ' <<<+- constitutive init -+>>>'; flush(IO_STDOUT) - call mech_init - call damage_init - call thermal_init - phases => config_material%get('phase') + + call mech_init(phases) + call damage_init + call thermal_init(phases) + + constitutive_source_maxSizeDotState = 0 PhaseLoop2:do ph = 1,phases%length !-------------------------------------------------------------------------------------------------- @@ -480,48 +453,18 @@ subroutine constitutive_init plasticState(ph)%partitionedState0 = plasticState(ph)%state0 plasticState(ph)%state = plasticState(ph)%partitionedState0 forall(so = 1:phase_Nsources(ph)) - sourceState(ph)%p(so)%partitionedState0 = sourceState(ph)%p(so)%state0 - sourceState(ph)%p(so)%state = sourceState(ph)%p(so)%partitionedState0 + damageState(ph)%p(so)%partitionedState0 = damageState(ph)%p(so)%state0 + damageState(ph)%p(so)%state = damageState(ph)%p(so)%partitionedState0 end forall constitutive_source_maxSizeDotState = max(constitutive_source_maxSizeDotState, & - maxval(sourceState(ph)%p%sizeDotState)) + maxval(damageState(ph)%p%sizeDotState)) enddo PhaseLoop2 constitutive_plasticity_maxSizeDotState = maxval(plasticState%sizeDotState) end subroutine constitutive_init -!-------------------------------------------------------------------------------------------------- -!> @brief checks if a source mechanism is active or not -!-------------------------------------------------------------------------------------------------- -function source_active(source_label,src_length) result(active_source) - - character(len=*), intent(in) :: source_label !< name of source mechanism - integer, intent(in) :: src_length !< max. number of sources in system - logical, dimension(:,:), allocatable :: active_source - - class(tNode), pointer :: & - phases, & - phase, & - sources, & - src - integer :: p,s - - phases => config_material%get('phase') - allocate(active_source(src_length,phases%length), source = .false. ) - do p = 1, phases%length - phase => phases%get(p) - sources => phase%get('source',defaultVal=emptyList) - do s = 1, sources%length - src => sources%get(s) - if(src%get_asString('type') == source_label) active_source(s,p) = .true. - enddo - enddo - - -end function source_active - !-------------------------------------------------------------------------------------------------- !> @brief checks if a kinematic mechanism is active or not @@ -554,197 +497,6 @@ function kinematics_active(kinematics_label,kinematics_length) result(active_ki end function kinematics_active -!-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the velocity gradient -! ToDo: MD: S is Mi? -!-------------------------------------------------------------------------------------------------- -subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & - S, Fi, co, ip, el) - - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el !< element - real(pReal), intent(in), dimension(3,3) :: & - S !< 2nd Piola-Kirchhoff stress - real(pReal), intent(in), dimension(3,3) :: & - Fi !< intermediate deformation gradient - real(pReal), intent(out), dimension(3,3) :: & - Li !< intermediate velocity gradient - real(pReal), intent(out), dimension(3,3,3,3) :: & - dLi_dS, & !< derivative of Li with respect to S - dLi_dFi - - real(pReal), dimension(3,3) :: & - my_Li, & !< intermediate velocity gradient - FiInv, & - temp_33 - real(pReal), dimension(3,3,3,3) :: & - my_dLi_dS - real(pReal) :: & - detFi - integer :: & - k, i, j, & - instance, of - - Li = 0.0_pReal - dLi_dS = 0.0_pReal - dLi_dFi = 0.0_pReal - - plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) - case (PLASTICITY_isotropic_ID) plasticityType - of = material_phasememberAt(co,ip,el) - instance = phase_plasticityInstance(material_phaseAt(co,el)) - call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,instance,of) - case default plasticityType - my_Li = 0.0_pReal - my_dLi_dS = 0.0_pReal - end select plasticityType - - Li = Li + my_Li - dLi_dS = dLi_dS + my_dLi_dS - - KinematicsLoop: do k = 1, phase_Nkinematics(material_phaseAt(co,el)) - kinematicsType: select case (phase_kinematics(k,material_phaseAt(co,el))) - case (KINEMATICS_cleavage_opening_ID) kinematicsType - call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el) - case (KINEMATICS_slipplane_opening_ID) kinematicsType - call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el) - case (KINEMATICS_thermal_expansion_ID) kinematicsType - call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dS, co, ip, el) - case default kinematicsType - my_Li = 0.0_pReal - my_dLi_dS = 0.0_pReal - end select kinematicsType - Li = Li + my_Li - dLi_dS = dLi_dS + my_dLi_dS - enddo KinematicsLoop - - FiInv = math_inv33(Fi) - detFi = math_det33(Fi) - Li = matmul(matmul(Fi,Li),FiInv)*detFi !< push forward to intermediate configuration - temp_33 = matmul(FiInv,Li) - - do i = 1,3; do j = 1,3 - dLi_dS(1:3,1:3,i,j) = matmul(matmul(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi - dLi_dFi(1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i) - dLi_dFi(1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i) - enddo; enddo - -end subroutine constitutive_LiAndItsTangents - - -!-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the rate of change of microstructure -!-------------------------------------------------------------------------------------------------- -function constitutive_damage_collectDotState(co,ip,el,ph,of) result(broken) - - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el, & !< element - ph, & - of - integer :: & - so !< counter in source loop - logical :: broken - - - broken = .false. - - SourceLoop: do so = 1, phase_Nsources(ph) - - sourceType: select case (phase_source(so,ph)) - - case (SOURCE_damage_anisoBrittle_ID) sourceType - call source_damage_anisoBrittle_dotState(constitutive_mech_getS(co,ip,el), co, ip, el) ! correct stress? - - case (SOURCE_damage_isoDuctile_ID) sourceType - call source_damage_isoDuctile_dotState(co, ip, el) - - case (SOURCE_damage_anisoDuctile_ID) sourceType - call source_damage_anisoDuctile_dotState(co, ip, el) - - end select sourceType - - broken = broken .or. any(IEEE_is_NaN(sourceState(ph)%p(so)%dotState(:,of))) - - enddo SourceLoop - -end function constitutive_damage_collectDotState - - -!-------------------------------------------------------------------------------------------------- -!> @brief contains the constitutive equation for calculating the rate of change of microstructure -!-------------------------------------------------------------------------------------------------- -function constitutive_thermal_collectDotState(ph,me) result(broken) - - integer, intent(in) :: ph, me - logical :: broken - - integer :: i - - - broken = .false. - - SourceLoop: do i = 1, phase_Nsources(ph) - - if (phase_source(i,ph) == SOURCE_thermal_externalheat_ID) & - call source_thermal_externalheat_dotState(ph,me) - - broken = broken .or. any(IEEE_is_NaN(sourceState(ph)%p(i)%dotState(:,me))) - - enddo SourceLoop - -end function constitutive_thermal_collectDotState - - -!-------------------------------------------------------------------------------------------------- -!> @brief for constitutive models having an instantaneous change of state -!> will return false if delta state is not needed/supported by the constitutive model -!-------------------------------------------------------------------------------------------------- -function constitutive_damage_deltaState(Fe, co, ip, el, ph, of) result(broken) - - integer, intent(in) :: & - co, & !< component-ID of integration point - ip, & !< integration point - el, & !< element - ph, & - of - real(pReal), intent(in), dimension(3,3) :: & - Fe !< elastic deformation gradient - integer :: & - so, & - myOffset, & - mySize - logical :: & - broken - - - broken = .false. - - sourceLoop: do so = 1, phase_Nsources(ph) - - sourceType: select case (phase_source(so,ph)) - - case (SOURCE_damage_isoBrittle_ID) sourceType - call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(co,ip,el), Fe, & - co, ip, el) - broken = any(IEEE_is_NaN(sourceState(ph)%p(so)%deltaState(:,of))) - if(.not. broken) then - myOffset = sourceState(ph)%p(so)%offsetDeltaState - mySize = sourceState(ph)%p(so)%sizeDeltaState - sourceState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,of) = & - sourceState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,of) + sourceState(ph)%p(so)%deltaState(1:mySize,of) - endif - - end select sourceType - - enddo SourceLoop - -end function constitutive_damage_deltaState - - !-------------------------------------------------------------------------------------------------- !> @brief Allocate the components of the state structure for a given phase !-------------------------------------------------------------------------------------------------- @@ -795,12 +547,13 @@ subroutine constitutive_restore(ip,el,includeL) do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) do so = 1, phase_Nsources(material_phaseAt(co,el)) - sourceState(material_phaseAt(co,el))%p(so)%state( :,material_phasememberAt(co,ip,el)) = & - sourceState(material_phaseAt(co,el))%p(so)%partitionedState0(:,material_phasememberAt(co,ip,el)) + damageState(material_phaseAt(co,el))%p(so)%state( :,material_phasememberAt(co,ip,el)) = & + damageState(material_phaseAt(co,el))%p(so)%partitionedState0(:,material_phasememberAt(co,ip,el)) enddo enddo call mech_restore(ip,el,includeL) + call thermal_restore(ip,el) end subroutine constitutive_restore @@ -809,16 +562,17 @@ end subroutine constitutive_restore !> @brief Forward data after successful increment. ! ToDo: Any guessing for the current states possible? !-------------------------------------------------------------------------------------------------- -subroutine constitutive_forward +subroutine constitutive_forward() integer :: ph, so - call constitutive_mech_forward() + call mech_forward() + call thermal_forward() - do ph = 1, size(sourceState) + do ph = 1, size(damageState) do so = 1,phase_Nsources(ph) - sourceState(ph)%p(so)%state0 = sourceState(ph)%p(so)%state + damageState(ph)%p(so)%state0 = damageState(ph)%p(so)%state enddo; enddo end subroutine constitutive_forward @@ -827,7 +581,7 @@ end subroutine constitutive_forward !-------------------------------------------------------------------------------------------------- !> @brief writes constitutive results to HDF5 output file !-------------------------------------------------------------------------------------------------- -subroutine constitutive_results +subroutine constitutive_results() integer :: ph character(len=:), allocatable :: group @@ -851,7 +605,7 @@ end subroutine constitutive_results !-------------------------------------------------------------------------------------------------- !> @brief allocates and initialize per grain variables !-------------------------------------------------------------------------------------------------- -subroutine crystallite_init +subroutine crystallite_init() integer :: & ph, & @@ -863,7 +617,7 @@ subroutine crystallite_init cMax, & !< maximum number of integration point components iMax, & !< maximum number of integration points eMax !< maximum number of elements - + class(tNode), pointer :: & num_crystallite, & @@ -919,7 +673,10 @@ subroutine crystallite_init do ph = 1, phases%length do so = 1, phase_Nsources(ph) - allocate(sourceState(ph)%p(so)%subState0,source=sourceState(ph)%p(so)%state0) ! ToDo: hack + allocate(damageState(ph)%p(so)%subState0,source=damageState(ph)%p(so)%state0) ! ToDo: hack + enddo + do so = 1, thermal_Nsources(ph) + allocate(thermalState(ph)%p(so)%subState0,source=thermalState(ph)%p(so)%state0) ! ToDo: hack enddo enddo @@ -963,10 +720,12 @@ subroutine constitutive_initializeRestorationPoints(ip,el) me = material_phaseMemberAt(co,ip,el) call mech_initializeRestorationPoints(ph,me) + call thermal_initializeRestorationPoints(ph,me) - do so = 1, phase_Nsources(material_phaseAt(co,el)) - sourceState(ph)%p(so)%partitionedState0(:,me) = sourceState(ph)%p(so)%state0(:,me) + do so = 1, size(damageState(ph)%p) + damageState(ph)%p(so)%partitionedState0(:,me) = damageState(ph)%p(so)%state0(:,me) enddo + enddo end subroutine constitutive_initializeRestorationPoints @@ -990,10 +749,13 @@ subroutine constitutive_windForward(ip,el) ph = material_phaseAt(co,el) me = material_phaseMemberAt(co,ip,el) - call constitutive_mech_windForward(ph,me) + call mech_windForward(ph,me) + call thermal_windForward(ph,me) + do so = 1, phase_Nsources(material_phaseAt(co,el)) - sourceState(ph)%p(so)%partitionedState0(:,me) = sourceState(ph)%p(so)%state(:,me) + damageState(ph)%p(so)%partitionedState0(:,me) = damageState(ph)%p(so)%state(:,me) enddo + enddo end subroutine constitutive_windForward @@ -1011,7 +773,7 @@ subroutine crystallite_orientations(co,ip,el) call crystallite_orientation(co,ip,el)%fromMatrix(transpose(math_rotationalPart(& - constitutive_mech_getF_e(co,ip,el)))) + mech_F_e(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el))))) if (plasticState(material_phaseAt(1,el))%nonlocal) & call plastic_nonlocal_updateCompatibility(crystallite_orientation, & @@ -1043,109 +805,7 @@ function crystallite_push33ToRef(co,ip,el, tensor33) end function crystallite_push33ToRef -!-------------------------------------------------------------------------------------------------- -!> @brief integrate stress, state with adaptive 1st order explicit Euler method -!> using Fixed Point Iteration to adapt the stepsize -!-------------------------------------------------------------------------------------------------- -function integrateSourceState(dt,co,ip,el) result(broken) - real(pReal), intent(in) :: dt - integer, intent(in) :: & - el, & !< element index in element loop - ip, & !< integration point index in ip loop - co !< grain index in grain loop - - integer :: & - NiterationState, & !< number of iterations in state loop - ph, & - me, & - so - integer, dimension(maxval(phase_Nsources)) :: & - size_so - real(pReal) :: & - zeta - real(pReal), dimension(constitutive_source_maxSizeDotState) :: & - r ! state residuum - real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState - logical :: & - broken, converged_ - - - ph = material_phaseAt(co,el) - me = material_phaseMemberAt(co,ip,el) - - converged_ = .true. - broken = constitutive_thermal_collectDotState(ph,me) - broken = broken .or. constitutive_damage_collectDotState(co,ip,el,ph,me) - if(broken) return - - do so = 1, phase_Nsources(ph) - size_so(so) = sourceState(ph)%p(so)%sizeDotState - sourceState(ph)%p(so)%state(1:size_so(so),me) = sourceState(ph)%p(so)%subState0(1:size_so(so),me) & - + sourceState(ph)%p(so)%dotState (1:size_so(so),me) * dt - source_dotState(1:size_so(so),2,so) = 0.0_pReal - enddo - - iteration: do NiterationState = 1, num%nState - - do so = 1, phase_Nsources(ph) - if(nIterationState > 1) source_dotState(1:size_so(so),2,so) = source_dotState(1:size_so(so),1,so) - source_dotState(1:size_so(so),1,so) = sourceState(ph)%p(so)%dotState(:,me) - enddo - - broken = constitutive_thermal_collectDotState(ph,me) - broken = broken .or. constitutive_damage_collectDotState(co,ip,el,ph,me) - if(broken) exit iteration - - do so = 1, phase_Nsources(ph) - zeta = damper(sourceState(ph)%p(so)%dotState(:,me), & - source_dotState(1:size_so(so),1,so),& - source_dotState(1:size_so(so),2,so)) - sourceState(ph)%p(so)%dotState(:,me) = sourceState(ph)%p(so)%dotState(:,me) * zeta & - + source_dotState(1:size_so(so),1,so)* (1.0_pReal - zeta) - r(1:size_so(so)) = sourceState(ph)%p(so)%state (1:size_so(so),me) & - - sourceState(ph)%p(so)%subState0(1:size_so(so),me) & - - sourceState(ph)%p(so)%dotState (1:size_so(so),me) * dt - sourceState(ph)%p(so)%state(1:size_so(so),me) = sourceState(ph)%p(so)%state(1:size_so(so),me) & - - r(1:size_so(so)) - converged_ = converged_ .and. converged(r(1:size_so(so)), & - sourceState(ph)%p(so)%state(1:size_so(so),me), & - sourceState(ph)%p(so)%atol(1:size_so(so))) - enddo - - if(converged_) then - broken = constitutive_damage_deltaState(constitutive_mech_getF_e(co,ip,el),co,ip,el,ph,me) - exit iteration - endif - - enddo iteration - - broken = broken .or. .not. converged_ - - - contains - - !-------------------------------------------------------------------------------------------------- - !> @brief calculate the damping for correction of state and dot state - !-------------------------------------------------------------------------------------------------- - real(pReal) pure function damper(current,previous,previous2) - - real(pReal), dimension(:), intent(in) ::& - current, previous, previous2 - - real(pReal) :: dot_prod12, dot_prod22 - - dot_prod12 = dot_product(current - previous, previous - previous2) - dot_prod22 = dot_product(previous - previous2, previous - previous2) - if ((dot_product(current,previous) < 0.0_pReal .or. dot_prod12 < 0.0_pReal) .and. dot_prod22 > 0.0_pReal) then - damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) - else - damper = 1.0_pReal - endif - - end function damper - -end function integrateSourceState !-------------------------------------------------------------------------------------------------- diff --git a/src/constitutive_damage.f90 b/src/constitutive_damage.f90 index 3ce614666..8c9104946 100644 --- a/src/constitutive_damage.f90 +++ b/src/constitutive_damage.f90 @@ -2,6 +2,16 @@ !> @brief internal microstructure state for all damage sources and kinematics constitutive models !---------------------------------------------------------------------------------------------------- submodule(constitutive) constitutive_damage + enum, bind(c); enumerator :: & + DAMAGE_UNDEFINED_ID, & + DAMAGE_ISOBRITTLE_ID, & + DAMAGE_ISODUCTILE_ID, & + DAMAGE_ANISOBRITTLE_ID, & + DAMAGE_ANISODUCTILE_ID + end enum + + integer(kind(DAMAGE_UNDEFINED_ID)), dimension(:,:), allocatable :: & + phase_source !< active sources mechanisms of each phase interface @@ -119,24 +129,24 @@ module subroutine damage_init phases => config_material%get('phase') - allocate(sourceState (phases%length)) + allocate(damageState (phases%length)) allocate(phase_Nsources(phases%length),source = 0) ! same for kinematics do ph = 1,phases%length phase => phases%get(ph) sources => phase%get('source',defaultVal=emptyList) phase_Nsources(ph) = sources%length - allocate(sourceState(ph)%p(phase_Nsources(ph))) + allocate(damageState(ph)%p(phase_Nsources(ph))) enddo - allocate(phase_source(maxval(phase_Nsources),phases%length), source = SOURCE_undefined_ID) + allocate(phase_source(maxval(phase_Nsources),phases%length), source = DAMAGE_UNDEFINED_ID) ! initialize source mechanisms if(maxval(phase_Nsources) /= 0) then - where(source_damage_isoBrittle_init (maxval(phase_Nsources))) phase_source = SOURCE_damage_isoBrittle_ID - where(source_damage_isoDuctile_init (maxval(phase_Nsources))) phase_source = SOURCE_damage_isoDuctile_ID - where(source_damage_anisoBrittle_init (maxval(phase_Nsources))) phase_source = SOURCE_damage_anisoBrittle_ID - where(source_damage_anisoDuctile_init (maxval(phase_Nsources))) phase_source = SOURCE_damage_anisoDuctile_ID + where(source_damage_isoBrittle_init (maxval(phase_Nsources))) phase_source = DAMAGE_ISOBRITTLE_ID + where(source_damage_isoDuctile_init (maxval(phase_Nsources))) phase_source = DAMAGE_ISODUCTILE_ID + where(source_damage_anisoBrittle_init (maxval(phase_Nsources))) phase_source = DAMAGE_ANISOBRITTLE_ID + where(source_damage_anisoDuctile_init (maxval(phase_Nsources))) phase_source = DAMAGE_ANISODUCTILE_ID endif !-------------------------------------------------------------------------------------------------- @@ -189,16 +199,16 @@ module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi constituent = material_phasememberAt(grain,ip,el) do source = 1, phase_Nsources(phase) select case(phase_source(source,phase)) - case (SOURCE_damage_isoBrittle_ID) + case (DAMAGE_ISOBRITTLE_ID) call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - case (SOURCE_damage_isoDuctile_ID) + case (DAMAGE_ISODUCTILE_ID) call source_damage_isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - case (SOURCE_damage_anisoBrittle_ID) + case (DAMAGE_ANISOBRITTLE_ID) call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) - case (SOURCE_damage_anisoDuctile_ID) + case (DAMAGE_ANISODUCTILE_ID) call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) case default @@ -214,6 +224,111 @@ module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi end subroutine constitutive_damage_getRateAndItsTangents + +!-------------------------------------------------------------------------------------------------- +!> @brief integrate stress, state with adaptive 1st order explicit Euler method +!> using Fixed Point Iteration to adapt the stepsize +!-------------------------------------------------------------------------------------------------- +module function integrateDamageState(dt,co,ip,el) result(broken) + + real(pReal), intent(in) :: dt + integer, intent(in) :: & + el, & !< element index in element loop + ip, & !< integration point index in ip loop + co !< grain index in grain loop + logical :: broken + + integer :: & + NiterationState, & !< number of iterations in state loop + ph, & + me, & + so + integer, dimension(maxval(phase_Nsources)) :: & + size_so + real(pReal) :: & + zeta + real(pReal), dimension(constitutive_source_maxSizeDotState) :: & + r ! state residuum + real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState + logical :: & + converged_ + + + ph = material_phaseAt(co,el) + me = material_phaseMemberAt(co,ip,el) + + converged_ = .true. + broken = constitutive_damage_collectDotState(co,ip,el,ph,me) + if(broken) return + + do so = 1, phase_Nsources(ph) + size_so(so) = damageState(ph)%p(so)%sizeDotState + damageState(ph)%p(so)%state(1:size_so(so),me) = damageState(ph)%p(so)%subState0(1:size_so(so),me) & + + damageState(ph)%p(so)%dotState (1:size_so(so),me) * dt + source_dotState(1:size_so(so),2,so) = 0.0_pReal + enddo + + iteration: do NiterationState = 1, num%nState + + do so = 1, phase_Nsources(ph) + if(nIterationState > 1) source_dotState(1:size_so(so),2,so) = source_dotState(1:size_so(so),1,so) + source_dotState(1:size_so(so),1,so) = damageState(ph)%p(so)%dotState(:,me) + enddo + + broken = constitutive_damage_collectDotState(co,ip,el,ph,me) + if(broken) exit iteration + + do so = 1, phase_Nsources(ph) + zeta = damper(damageState(ph)%p(so)%dotState(:,me), & + source_dotState(1:size_so(so),1,so),& + source_dotState(1:size_so(so),2,so)) + damageState(ph)%p(so)%dotState(:,me) = damageState(ph)%p(so)%dotState(:,me) * zeta & + + source_dotState(1:size_so(so),1,so)* (1.0_pReal - zeta) + r(1:size_so(so)) = damageState(ph)%p(so)%state (1:size_so(so),me) & + - damageState(ph)%p(so)%subState0(1:size_so(so),me) & + - damageState(ph)%p(so)%dotState (1:size_so(so),me) * dt + damageState(ph)%p(so)%state(1:size_so(so),me) = damageState(ph)%p(so)%state(1:size_so(so),me) & + - r(1:size_so(so)) + converged_ = converged_ .and. converged(r(1:size_so(so)), & + damageState(ph)%p(so)%state(1:size_so(so),me), & + damageState(ph)%p(so)%atol(1:size_so(so))) + enddo + + if(converged_) then + broken = constitutive_damage_deltaState(mech_F_e(ph,me),co,ip,el,ph,me) + exit iteration + endif + + enddo iteration + + broken = broken .or. .not. converged_ + + + contains + + !-------------------------------------------------------------------------------------------------- + !> @brief calculate the damping for correction of state and dot state + !-------------------------------------------------------------------------------------------------- + real(pReal) pure function damper(current,previous,previous2) + + real(pReal), dimension(:), intent(in) ::& + current, previous, previous2 + + real(pReal) :: dot_prod12, dot_prod22 + + dot_prod12 = dot_product(current - previous, previous - previous2) + dot_prod22 = dot_product(previous - previous2, previous - previous2) + if ((dot_product(current,previous) < 0.0_pReal .or. dot_prod12 < 0.0_pReal) .and. dot_prod22 > 0.0_pReal) then + damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22) + else + damper = 1.0_pReal + endif + + end function damper + +end function integrateDamageState + + !---------------------------------------------------------------------------------------------- !< @brief writes damage sources results to HDF5 output file !---------------------------------------------------------------------------------------------- @@ -226,23 +341,23 @@ module subroutine damage_results(group,ph) sourceLoop: do so = 1, phase_Nsources(ph) - if (phase_source(so,ph) /= SOURCE_UNDEFINED_ID) & + if (phase_source(so,ph) /= DAMAGE_UNDEFINED_ID) & call results_closeGroup(results_addGroup(group//'sources/')) ! should be 'damage' sourceType: select case (phase_source(so,ph)) - case (SOURCE_damage_anisoBrittle_ID) sourceType - call source_damage_anisoBrittle_results(ph,group//'sources/') - - case (SOURCE_damage_anisoDuctile_ID) sourceType - call source_damage_anisoDuctile_results(ph,group//'sources/') - - case (SOURCE_damage_isoBrittle_ID) sourceType + case (DAMAGE_ISOBRITTLE_ID) sourceType call source_damage_isoBrittle_results(ph,group//'sources/') - case (SOURCE_damage_isoDuctile_ID) sourceType + case (DAMAGE_ISODUCTILE_ID) sourceType call source_damage_isoDuctile_results(ph,group//'sources/') + case (DAMAGE_ANISOBRITTLE_ID) sourceType + call source_damage_anisoBrittle_results(ph,group//'sources/') + + case (DAMAGE_ANISODUCTILE_ID) sourceType + call source_damage_anisoDuctile_results(ph,group//'sources/') + end select sourceType enddo SourceLoop @@ -250,4 +365,123 @@ module subroutine damage_results(group,ph) end subroutine damage_results +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the rate of change of microstructure +!-------------------------------------------------------------------------------------------------- +function constitutive_damage_collectDotState(co,ip,el,ph,of) result(broken) + + integer, intent(in) :: & + co, & !< component-ID of integration point + ip, & !< integration point + el, & !< element + ph, & + of + integer :: & + so !< counter in source loop + logical :: broken + + + broken = .false. + + SourceLoop: do so = 1, phase_Nsources(ph) + + sourceType: select case (phase_source(so,ph)) + + case (DAMAGE_ISODUCTILE_ID) sourceType + call source_damage_isoDuctile_dotState(co, ip, el) + + case (DAMAGE_ANISODUCTILE_ID) sourceType + call source_damage_anisoDuctile_dotState(co, ip, el) + + case (DAMAGE_ANISOBRITTLE_ID) sourceType + call source_damage_anisoBrittle_dotState(mech_S(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el)),& + co, ip, el) ! correct stress? + + end select sourceType + + broken = broken .or. any(IEEE_is_NaN(damageState(ph)%p(so)%dotState(:,of))) + + enddo SourceLoop + +end function constitutive_damage_collectDotState + + + +!-------------------------------------------------------------------------------------------------- +!> @brief for constitutive models having an instantaneous change of state +!> will return false if delta state is not needed/supported by the constitutive model +!-------------------------------------------------------------------------------------------------- +function constitutive_damage_deltaState(Fe, co, ip, el, ph, of) result(broken) + + integer, intent(in) :: & + co, & !< component-ID of integration point + ip, & !< integration point + el, & !< element + ph, & + of + real(pReal), intent(in), dimension(3,3) :: & + Fe !< elastic deformation gradient + integer :: & + so, & + myOffset, & + mySize + logical :: & + broken + + + broken = .false. + + sourceLoop: do so = 1, phase_Nsources(ph) + + sourceType: select case (phase_source(so,ph)) + + case (DAMAGE_ISOBRITTLE_ID) sourceType + call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(co,ip,el), Fe, & + co, ip, el) + broken = any(IEEE_is_NaN(damageState(ph)%p(so)%deltaState(:,of))) + if(.not. broken) then + myOffset = damageState(ph)%p(so)%offsetDeltaState + mySize = damageState(ph)%p(so)%sizeDeltaState + damageState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,of) = & + damageState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,of) + damageState(ph)%p(so)%deltaState(1:mySize,of) + endif + + end select sourceType + + enddo SourceLoop + +end function constitutive_damage_deltaState + + +!-------------------------------------------------------------------------------------------------- +!> @brief checks if a source mechanism is active or not +!-------------------------------------------------------------------------------------------------- +function source_active(source_label,src_length) result(active_source) + + character(len=*), intent(in) :: source_label !< name of source mechanism + integer, intent(in) :: src_length !< max. number of sources in system + logical, dimension(:,:), allocatable :: active_source + + class(tNode), pointer :: & + phases, & + phase, & + sources, & + src + integer :: p,s + + phases => config_material%get('phase') + allocate(active_source(src_length,phases%length), source = .false. ) + do p = 1, phases%length + phase => phases%get(p) + sources => phase%get('source',defaultVal=emptyList) + do s = 1, sources%length + src => sources%get(s) + if(src%get_asString('type') == source_label) active_source(s,p) = .true. + enddo + enddo + + +end function source_active + + end submodule constitutive_damage diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index 44e777c80..8bc85354f 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -7,12 +7,20 @@ submodule(constitutive) constitutive_mech ELASTICITY_UNDEFINED_ID, & ELASTICITY_HOOKE_ID, & STIFFNESS_DEGRADATION_UNDEFINED_ID, & - STIFFNESS_DEGRADATION_DAMAGE_ID + STIFFNESS_DEGRADATION_DAMAGE_ID, & + PLASTICITY_UNDEFINED_ID, & + PLASTICITY_NONE_ID, & + PLASTICITY_ISOTROPIC_ID, & + PLASTICITY_PHENOPOWERLAW_ID, & + PLASTICITY_KINEHARDENING_ID, & + PLASTICITY_DISLOTWIN_ID, & + PLASTICITY_DISLOTUNGSTEN_ID, & + PLASTICITY_NONLOCAL_ID end enum - integer(kind(ELASTICITY_undefined_ID)), dimension(:), allocatable :: & + integer(kind(ELASTICITY_UNDEFINED_ID)), dimension(:), allocatable :: & phase_elasticity !< elasticity of each phase - integer(kind(SOURCE_undefined_ID)), dimension(:,:), allocatable :: & + integer(kind(STIFFNESS_DEGRADATION_UNDEFINED_ID)), dimension(:,:), allocatable :: & phase_stiffnessDegradation !< active stiffness degradation mechanisms of each phase type(tTensorContainer), dimension(:), allocatable :: & @@ -41,6 +49,12 @@ submodule(constitutive) constitutive_mech constitutive_mech_partitionedS0 + + + integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable :: & + phase_plasticity !< plasticity of each phase + + interface module function plastic_none_init() result(myPlasticity) @@ -319,7 +333,10 @@ contains !> @brief Initialize mechanical field related constitutive models !> @details Initialize elasticity, plasticity and stiffness degradation models. !-------------------------------------------------------------------------------------------------- -module subroutine mech_init +module subroutine mech_init(phases) + + class(tNode), pointer :: & + phases integer :: & el, & @@ -331,7 +348,6 @@ module subroutine mech_init Nconstituents class(tNode), pointer :: & num_crystallite, & - phases, & phase, & mech, & elastic, & @@ -341,7 +357,6 @@ module subroutine mech_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) @@ -529,7 +544,7 @@ end function plastic_active !> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to !> the elastic and intermediate deformation gradients using Hooke's law !-------------------------------------------------------------------------------------------------- -module subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & +subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & Fe, Fi, co, ip, el) integer, intent(in) :: & @@ -615,7 +630,7 @@ end subroutine constitutive_plastic_dependentState ! ToDo: Discuss whether it makes sense if crystallite handles the configuration conversion, i.e. ! Mp in, dLp_dMp out !-------------------------------------------------------------------------------------------------- -module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & +subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, & S, Fi, co, ip, el) integer, intent(in) :: & co, & !< component-ID of integration point @@ -1485,7 +1500,7 @@ end subroutine mech_initializeRestorationPoints !-------------------------------------------------------------------------------------------------- !> @brief Wind homog inc forward. !-------------------------------------------------------------------------------------------------- -module subroutine constitutive_mech_windForward(ph,me) +module subroutine mech_windForward(ph,me) integer, intent(in) :: ph, me @@ -1499,14 +1514,14 @@ module subroutine constitutive_mech_windForward(ph,me) plasticState(ph)%partitionedState0(:,me) = plasticState(ph)%state(:,me) -end subroutine constitutive_mech_windForward +end subroutine mech_windForward !-------------------------------------------------------------------------------------------------- !> @brief Forward data after successful increment. ! ToDo: Any guessing for the current states possible? !-------------------------------------------------------------------------------------------------- -module subroutine constitutive_mech_forward() +module subroutine mech_forward() integer :: ph @@ -1521,7 +1536,7 @@ module subroutine constitutive_mech_forward() plasticState(ph)%state0 = plasticState(ph)%state enddo -end subroutine constitutive_mech_forward +end subroutine mech_forward @@ -1585,7 +1600,10 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) do so = 1, phase_Nsources(ph) - sourceState(ph)%p(so)%subState0(:,me) = sourceState(ph)%p(so)%partitionedState0(:,me) + damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%partitionedState0(:,me) + enddo + do so = 1, thermal_Nsources(ph) + thermalState(ph)%p(so)%subState0(:,me) = thermalState(ph)%p(so)%partitionedState0(:,me) enddo subFp0 = constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) subFi0 = constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) @@ -1613,7 +1631,10 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) subFi0 = constitutive_mech_Fi(ph)%data(1:3,1:3,me) subState0 = plasticState(ph)%state(:,me) do so = 1, phase_Nsources(ph) - sourceState(ph)%p(so)%subState0(:,me) = sourceState(ph)%p(so)%state(:,me) + damageState(ph)%p(so)%subState0(:,me) = damageState(ph)%p(so)%state(:,me) + enddo + do so = 1, thermal_Nsources(ph) + thermalState(ph)%p(so)%subState0(:,me) = thermalState(ph)%p(so)%state(:,me) enddo endif !-------------------------------------------------------------------------------------------------- @@ -1629,7 +1650,10 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) endif plasticState(ph)%state(:,me) = subState0 do so = 1, phase_Nsources(ph) - sourceState(ph)%p(so)%state(:,me) = sourceState(ph)%p(so)%subState0(:,me) + damageState(ph)%p(so)%state(:,me) = damageState(ph)%p(so)%subState0(:,me) + enddo + do so = 1, thermal_Nsources(ph) + thermalState(ph)%p(so)%state(:,me) = thermalState(ph)%p(so)%subState0(:,me) enddo todo = subStep > num%subStepMinCryst ! still on track or already done (beyond repair) @@ -1643,7 +1667,8 @@ module function crystallite_stress(dt,co,ip,el) result(converged_) constitutive_mech_Fe(ph)%data(1:3,1:3,me) = matmul(subF,math_inv33(matmul(constitutive_mech_Fi(ph)%data(1:3,1:3,me), & constitutive_mech_Fp(ph)%data(1:3,1:3,me)))) converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * dt,co,ip,el) - converged_ = converged_ .and. .not. integrateSourceState(subStep * dt,co,ip,el) + converged_ = converged_ .and. .not. integrateDamageState(subStep * dt,co,ip,el) + converged_ = converged_ .and. .not. integrateThermalState(subStep * dt,co,ip,el) endif enddo cutbackLooping @@ -1678,8 +1703,7 @@ module subroutine mech_restore(ip,el,includeL) constitutive_mech_Fi(ph)%data(1:3,1:3,me) = constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) constitutive_mech_S(ph)%data(1:3,1:3,me) = constitutive_mech_partitionedS0(ph)%data(1:3,1:3,me) - plasticState (material_phaseAt(co,el))%state( :,material_phasememberAt(co,ip,el)) = & - plasticState (material_phaseAt(co,el))%partitionedState0(:,material_phasememberAt(co,ip,el)) + plasticState(ph)%state(:,me) = plasticState(ph)%partitionedState0(:,me) enddo end subroutine mech_restore @@ -1846,31 +1870,37 @@ module subroutine mech_restartRead(groupHandle,ph) end subroutine mech_restartRead -! getter for non-mech (e.g. thermal) -module function constitutive_mech_getS(co,ip,el) result(S) +!---------------------------------------------------------------------------------------------- +!< @brief Get first Piola-Kichhoff stress (for use by non-mech physics) +!---------------------------------------------------------------------------------------------- +module function mech_S(ph,me) result(S) - integer, intent(in) :: co, ip, el + integer, intent(in) :: ph,me real(pReal), dimension(3,3) :: S - S = constitutive_mech_S(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el)) + S = constitutive_mech_S(ph)%data(1:3,1:3,me) -end function constitutive_mech_getS +end function mech_S -! getter for non-mech (e.g. thermal) -module function constitutive_mech_getLp(co,ip,el) result(Lp) +!---------------------------------------------------------------------------------------------- +!< @brief Get plastic velocity gradient (for use by non-mech physics) +!---------------------------------------------------------------------------------------------- +module function mech_L_p(ph,me) result(L_p) - integer, intent(in) :: co, ip, el - real(pReal), dimension(3,3) :: Lp + integer, intent(in) :: ph,me + real(pReal), dimension(3,3) :: L_p - Lp = constitutive_mech_Lp(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el)) + L_p = constitutive_mech_Lp(ph)%data(1:3,1:3,me) -end function constitutive_mech_getLp +end function mech_L_p -! getter for non-mech (e.g. thermal) +!---------------------------------------------------------------------------------------------- +!< @brief Get deformation gradient (for use by homogenization) +!---------------------------------------------------------------------------------------------- module function constitutive_mech_getF(co,ip,el) result(F) integer, intent(in) :: co, ip, el @@ -1882,20 +1912,24 @@ module function constitutive_mech_getF(co,ip,el) result(F) end function constitutive_mech_getF -! getter for non-mech (e.g. thermal) -module function constitutive_mech_getF_e(co,ip,el) result(F_e) +!---------------------------------------------------------------------------------------------- +!< @brief Get elastic deformation gradient (for use by non-mech physics) +!---------------------------------------------------------------------------------------------- +module function mech_F_e(ph,me) result(F_e) - integer, intent(in) :: co, ip, el + integer, intent(in) :: ph,me real(pReal), dimension(3,3) :: F_e - F_e = constitutive_mech_Fe(material_phaseAt(co,el))%data(1:3,1:3,material_phaseMemberAt(co,ip,el)) + F_e = constitutive_mech_Fe(ph)%data(1:3,1:3,me) -end function constitutive_mech_getF_e +end function mech_F_e -! getter for non-mech (e.g. thermal) +!---------------------------------------------------------------------------------------------- +!< @brief Get second Piola-Kichhoff stress (for use by homogenization) +!---------------------------------------------------------------------------------------------- module function constitutive_mech_getP(co,ip,el) result(P) integer, intent(in) :: co, ip, el @@ -1918,5 +1952,85 @@ module subroutine constitutive_mech_setF(F,co,ip,el) end subroutine constitutive_mech_setF + +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the velocity gradient +! ToDo: MD: S is Mi? +!-------------------------------------------------------------------------------------------------- +subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & + S, Fi, co, ip, el) + + integer, intent(in) :: & + co, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in), dimension(3,3) :: & + S !< 2nd Piola-Kirchhoff stress + real(pReal), intent(in), dimension(3,3) :: & + Fi !< intermediate deformation gradient + real(pReal), intent(out), dimension(3,3) :: & + Li !< intermediate velocity gradient + real(pReal), intent(out), dimension(3,3,3,3) :: & + dLi_dS, & !< derivative of Li with respect to S + dLi_dFi + + real(pReal), dimension(3,3) :: & + my_Li, & !< intermediate velocity gradient + FiInv, & + temp_33 + real(pReal), dimension(3,3,3,3) :: & + my_dLi_dS + real(pReal) :: & + detFi + integer :: & + k, i, j, & + instance, of + + Li = 0.0_pReal + dLi_dS = 0.0_pReal + dLi_dFi = 0.0_pReal + + plasticityType: select case (phase_plasticity(material_phaseAt(co,el))) + case (PLASTICITY_isotropic_ID) plasticityType + of = material_phasememberAt(co,ip,el) + instance = phase_plasticityInstance(material_phaseAt(co,el)) + call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S ,instance,of) + case default plasticityType + my_Li = 0.0_pReal + my_dLi_dS = 0.0_pReal + end select plasticityType + + Li = Li + my_Li + dLi_dS = dLi_dS + my_dLi_dS + + KinematicsLoop: do k = 1, phase_Nkinematics(material_phaseAt(co,el)) + kinematicsType: select case (phase_kinematics(k,material_phaseAt(co,el))) + case (KINEMATICS_cleavage_opening_ID) kinematicsType + call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el) + case (KINEMATICS_slipplane_opening_ID) kinematicsType + call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el) + case (KINEMATICS_thermal_expansion_ID) kinematicsType + call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dS, co, ip, el) + case default kinematicsType + my_Li = 0.0_pReal + my_dLi_dS = 0.0_pReal + end select kinematicsType + Li = Li + my_Li + dLi_dS = dLi_dS + my_dLi_dS + enddo KinematicsLoop + + FiInv = math_inv33(Fi) + detFi = math_det33(Fi) + Li = matmul(matmul(Fi,Li),FiInv)*detFi !< push forward to intermediate configuration + temp_33 = matmul(FiInv,Li) + + do i = 1,3; do j = 1,3 + dLi_dS(1:3,1:3,i,j) = matmul(matmul(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi + dLi_dFi(1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i) + dLi_dFi(1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i) + enddo; enddo + +end subroutine constitutive_LiAndItsTangents + end submodule constitutive_mech diff --git a/src/constitutive_thermal.f90 b/src/constitutive_thermal.f90 index 1a05b983f..9787cb0e4 100644 --- a/src/constitutive_thermal.f90 +++ b/src/constitutive_thermal.f90 @@ -3,6 +3,21 @@ !---------------------------------------------------------------------------------------------------- submodule(constitutive) constitutive_thermal + enum, bind(c); enumerator :: & + THERMAL_UNDEFINED_ID ,& + THERMAL_DISSIPATION_ID, & + THERMAL_EXTERNALHEAT_ID + end enum + + type :: tDataContainer + real(pReal), dimension(:), allocatable :: T + end type tDataContainer + integer(kind(THERMAL_UNDEFINED_ID)), dimension(:,:), allocatable :: & + thermal_source + + type(tDataContainer), dimension(:), allocatable :: current + + integer :: thermal_source_maxSizeDotState interface module function source_thermal_dissipation_init(source_length) result(mySources) @@ -21,7 +36,7 @@ submodule(constitutive) constitutive_thermal end function kinematics_thermal_expansion_init - module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase) + module subroutine thermal_dissipation_getRate(TDot, Tstar,Lp,phase) integer, intent(in) :: & phase !< phase ID of element real(pReal), intent(in), dimension(3,3) :: & @@ -29,18 +44,16 @@ submodule(constitutive) constitutive_thermal real(pReal), intent(in), dimension(3,3) :: & Lp !< plastic velocuty gradient for a given element real(pReal), intent(out) :: & - TDot, & - dTDot_dT - end subroutine source_thermal_dissipation_getRateAndItsTangent + TDot + end subroutine thermal_dissipation_getRate - module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phase, of) + module subroutine thermal_externalheat_getRate(TDot, phase,of) integer, intent(in) :: & phase, & of real(pReal), intent(out) :: & - TDot, & - dTDot_dT - end subroutine source_thermal_externalheat_getRateAndItsTangent + TDot + end subroutine thermal_externalheat_getRate end interface @@ -49,14 +62,60 @@ contains !---------------------------------------------------------------------------------------------- !< @brief initializes thermal sources and kinematics mechanism !---------------------------------------------------------------------------------------------- -module subroutine thermal_init +module subroutine thermal_init(phases) -! initialize source mechanisms - if(maxval(phase_Nsources) /= 0) then - where(source_thermal_dissipation_init (maxval(phase_Nsources))) phase_source = SOURCE_thermal_dissipation_ID - where(source_thermal_externalheat_init(maxval(phase_Nsources))) phase_source = SOURCE_thermal_externalheat_ID + class(tNode), pointer :: & + phases + + class(tNode), pointer :: & + phase, thermal, sources + + integer :: & + ph, so, & + Nconstituents + + + print'(/,a)', ' <<<+- constitutive_thermal init -+>>>' + + allocate(current(phases%length)) + + allocate(thermalState (phases%length)) + allocate(thermal_Nsources(phases%length),source = 0) + + do ph = 1, phases%length + + Nconstituents = count(material_phaseAt == ph) * discretization_nIPs + + allocate(current(ph)%T(Nconstituents)) + phase => phases%get(ph) + if(phase%contains('thermal')) then + thermal => phase%get('thermal') + sources => thermal%get('source',defaultVal=emptyList) + + thermal_Nsources(ph) = sources%length + endif + allocate(thermalstate(ph)%p(thermal_Nsources(ph))) + enddo + + allocate(thermal_source(maxval(thermal_Nsources),phases%length), source = THERMAL_UNDEFINED_ID) + + if(maxval(thermal_Nsources) /= 0) then + where(source_thermal_dissipation_init (maxval(thermal_Nsources))) thermal_source = THERMAL_DISSIPATION_ID + where(source_thermal_externalheat_init(maxval(thermal_Nsources))) thermal_source = THERMAL_EXTERNALHEAT_ID endif + thermal_source_maxSizeDotState = 0 + PhaseLoop2:do ph = 1,phases%length + + do so = 1,thermal_Nsources(ph) + thermalState(ph)%p(so)%partitionedState0 = thermalState(ph)%p(so)%state0 + thermalState(ph)%p(so)%state = thermalState(ph)%p(so)%partitionedState0 + enddo + + thermal_source_maxSizeDotState = max(thermal_source_maxSizeDotState, & + maxval(thermalState(ph)%p%sizeDotState)) + enddo PhaseLoop2 + !-------------------------------------------------------------------------------------------------- !initialize kinematic mechanisms if(maxval(phase_Nkinematics) /= 0) where(kinematics_thermal_expansion_init(maxval(phase_Nkinematics))) & @@ -68,75 +127,236 @@ end subroutine thermal_init !---------------------------------------------------------------------------------------------- !< @brief calculates thermal dissipation rate !---------------------------------------------------------------------------------------------- -module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, ip, el) +module subroutine constitutive_thermal_getRate(TDot, T, ip, el) integer, intent(in) :: & ip, & !< integration point number el !< element number real(pReal), intent(in) :: & T !< plastic velocity gradient - real(pReal), intent(inout) :: & - TDot, & - dTDot_dT + real(pReal), intent(out) :: & + TDot real(pReal) :: & - my_Tdot, & - my_dTdot_dT - real(pReal), dimension(3,3) :: Lp, S + my_Tdot integer :: & - phase, & + ph, & homog, & instance, & - grain, & - source, & - constituent + me, & + so, & + co homog = material_homogenizationAt(el) instance = thermal_typeInstance(homog) - do grain = 1, homogenization_Nconstituents(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) - Lp = constitutive_mech_getLp(grain,ip,el) - S = constitutive_mech_getS(grain,ip,el) - call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & - S, Lp, phase) + TDot = 0.0_pReal + do co = 1, homogenization_Nconstituents(homog) + ph = material_phaseAt(co,el) + me = material_phasememberAt(co,ip,el) + do so = 1, thermal_Nsources(ph) + select case(thermal_source(so,ph)) + case (THERMAL_DISSIPATION_ID) + call thermal_dissipation_getRate(my_Tdot, mech_S(ph,me),mech_L_p(ph,me),ph) - case (SOURCE_thermal_externalheat_ID) - call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & - phase, constituent) + case (THERMAL_EXTERNALHEAT_ID) + call thermal_externalheat_getRate(my_Tdot, ph,me) 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 constitutive_thermal_getRateAndItsTangents +end subroutine constitutive_thermal_getRate + + +!-------------------------------------------------------------------------------------------------- +!> @brief contains the constitutive equation for calculating the rate of change of microstructure +!-------------------------------------------------------------------------------------------------- +function constitutive_thermal_collectDotState(ph,me) result(broken) + + integer, intent(in) :: ph, me + logical :: broken + + integer :: i + + + broken = .false. + + SourceLoop: do i = 1, thermal_Nsources(ph) + + if (thermal_source(i,ph) == THERMAL_EXTERNALHEAT_ID) & + call source_thermal_externalheat_dotState(ph,me) + + broken = broken .or. any(IEEE_is_NaN(thermalState(ph)%p(i)%dotState(:,me))) + + enddo SourceLoop + +end function constitutive_thermal_collectDotState +!-------------------------------------------------------------------------------------------------- +!> @brief integrate state with 1st order explicit Euler method +!-------------------------------------------------------------------------------------------------- +module function integrateThermalState(Delta_t,co,ip,el) result(broken) -! getter for non-thermal (e.g. mech) -module function constitutive_thermal_T(co,ip,el) result(T) + real(pReal), intent(in) :: Delta_t + integer, intent(in) :: & + el, & !< element index in element loop + ip, & !< integration point index in ip loop + co !< grain index in grain loop + logical :: & + broken - integer, intent(in) :: co, ip, el + integer :: & + ph, & + me, & + so, & + sizeDotState + + + ph = material_phaseAt(co,el) + me = material_phaseMemberAt(co,ip,el) + + broken = constitutive_thermal_collectDotState(ph,me) + if(broken) return + + do so = 1, thermal_Nsources(ph) + sizeDotState = thermalState(ph)%p(so)%sizeDotState + thermalState(ph)%p(so)%state(1:sizeDotState,me) = thermalState(ph)%p(so)%subState0(1:sizeDotState,me) & + + thermalState(ph)%p(so)%dotState(1:sizeDotState,me) * Delta_t + enddo + +end function integrateThermalState + + +module subroutine thermal_initializeRestorationPoints(ph,me) + + integer, intent(in) :: ph, me + + integer :: so + + + do so = 1, size(thermalState(ph)%p) + thermalState(ph)%p(so)%partitionedState0(:,me) = thermalState(ph)%p(so)%state0(:,me) + enddo + +end subroutine thermal_initializeRestorationPoints + + + +module subroutine thermal_windForward(ph,me) + + integer, intent(in) :: ph, me + + integer :: so + + + do so = 1, size(thermalState(ph)%p) + thermalState(ph)%p(so)%partitionedState0(:,me) = thermalState(ph)%p(so)%state(:,me) + enddo + +end subroutine thermal_windForward + + +module subroutine thermal_forward() + + integer :: ph, so + + + do ph = 1, size(thermalState) + do so = 1, size(thermalState(ph)%p) + thermalState(ph)%p(so)%state0 = thermalState(ph)%p(so)%state + enddo + enddo + +end subroutine thermal_forward + + +module subroutine thermal_restore(ip,el) + + integer, intent(in) :: ip, el + + integer :: co, ph, me, so + + + do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) + ph = material_phaseAt(co,el) + me = material_phaseMemberAt(co,ip,el) + + do so = 1, size(thermalState(ph)%p) + thermalState(ph)%p(so)%state(:,me) = thermalState(ph)%p(so)%partitionedState0(:,me) + enddo + + enddo + +end subroutine thermal_restore + + +!---------------------------------------------------------------------------------------------- +!< @brief Get temperature (for use by non-thermal physics) +!---------------------------------------------------------------------------------------------- +module function thermal_T(ph,me) result(T) + + integer, intent(in) :: ph, me real(pReal) :: T - integer :: ho, tme - ho = material_homogenizationAt(el) - tme = material_homogenizationMemberAt(ip,el) + T = current(ph)%T(me) - T = temperature(ho)%p(tme) +end function thermal_T -end function constitutive_thermal_T + +!---------------------------------------------------------------------------------------------- +!< @brief Set temperature +!---------------------------------------------------------------------------------------------- +module subroutine constitutive_thermal_setT(T,co,ip,el) + + real(pReal), intent(in) :: T + integer, intent(in) :: co, ip, el + + + current(material_phaseAt(co,el))%T(material_phaseMemberAt(co,ip,el)) = T + +end subroutine constitutive_thermal_setT + + + +!-------------------------------------------------------------------------------------------------- +!> @brief checks if a source mechanism is active or not +!-------------------------------------------------------------------------------------------------- +function thermal_active(source_label,src_length) result(active_source) + + character(len=*), intent(in) :: source_label !< name of source mechanism + integer, intent(in) :: src_length !< max. number of sources in system + logical, dimension(:,:), allocatable :: active_source + + class(tNode), pointer :: & + phases, & + phase, & + sources, thermal, & + src + integer :: p,s + + phases => config_material%get('phase') + allocate(active_source(src_length,phases%length), source = .false. ) + do p = 1, phases%length + phase => phases%get(p) + if (phase%contains('thermal')) then + thermal => phase%get('thermal',defaultVal=emptyList) + sources => thermal%get('source',defaultVal=emptyList) + do s = 1, sources%length + src => sources%get(s) + if(src%get_asString('type') == source_label) active_source(s,p) = .true. + enddo + endif + enddo + + +end function thermal_active end submodule constitutive_thermal diff --git a/src/source_thermal_dissipation.f90 b/src/constitutive_thermal_dissipation.f90 similarity index 82% rename from src/source_thermal_dissipation.f90 rename to src/constitutive_thermal_dissipation.f90 index f28567aa7..ae2d5735e 100644 --- a/src/source_thermal_dissipation.f90 +++ b/src/constitutive_thermal_dissipation.f90 @@ -4,7 +4,7 @@ !> @brief material subroutine for thermal source due to plastic dissipation !> @details to be done !-------------------------------------------------------------------------------------------------- -submodule(constitutive:constitutive_thermal) source_thermal_dissipation +submodule(constitutive:constitutive_thermal) source_dissipation integer, dimension(:), allocatable :: & source_thermal_dissipation_offset, & !< which source is my current thermal dissipation mechanism? @@ -27,19 +27,20 @@ contains !-------------------------------------------------------------------------------------------------- module function source_thermal_dissipation_init(source_length) result(mySources) - integer, intent(in) :: source_length + integer, intent(in) :: source_length logical, dimension(:,:), allocatable :: mySources class(tNode), pointer :: & phases, & phase, & - sources, & - src + sources, thermal, & + src integer :: Ninstances,sourceOffset,Nconstituents,p - print'(/,a)', ' <<<+- source_thermal_dissipation init -+>>>' + print'(/,a)', ' <<<+- thermal_dissipation init -+>>>' + + mySources = thermal_active('dissipation',source_length) - mySources = source_active('thermal_dissipation',source_length) Ninstances = count(mySources) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) if(Ninstances == 0) return @@ -50,19 +51,20 @@ module function source_thermal_dissipation_init(source_length) result(mySources) allocate(source_thermal_dissipation_instance(phases%length), source=0) do p = 1, phases%length - phase => phases%get(p) - if(count(mySources(:,p)) == 0) cycle + phase => phases%get(p) if(any(mySources(:,p))) source_thermal_dissipation_instance(p) = count(mySources(:,1:p)) - sources => phase%get('source') + if(count(mySources(:,p)) == 0) cycle + thermal => phase%get('thermal') + sources => thermal%get('source') do sourceOffset = 1, sources%length if(mySources(sourceOffset,p)) then source_thermal_dissipation_offset(p) = sourceOffset associate(prm => param(source_thermal_dissipation_instance(p))) + src => sources%get(sourceOffset) - src => sources%get(sourceOffset) prm%kappa = src%get_asFloat('kappa') Nconstituents = count(material_phaseAt==p) * discretization_nIPs - call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,0,0,0) + call constitutive_allocateState(thermalState(p)%p(sourceOffset),Nconstituents,0,0,0) end associate endif @@ -76,7 +78,7 @@ end function source_thermal_dissipation_init !-------------------------------------------------------------------------------------------------- !> @brief Ninstancess dissipation rate !-------------------------------------------------------------------------------------------------- -module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase) +module subroutine thermal_dissipation_getRate(TDot, Tstar, Lp, phase) integer, intent(in) :: & phase @@ -86,14 +88,12 @@ module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT Lp real(pReal), intent(out) :: & - TDot, & - dTDot_dT + TDot associate(prm => param(source_thermal_dissipation_instance(phase))) TDot = prm%kappa*sum(abs(Tstar*Lp)) - dTDot_dT = 0.0_pReal end associate -end subroutine source_thermal_dissipation_getRateAndItsTangent +end subroutine thermal_dissipation_getRate -end submodule source_thermal_dissipation +end submodule source_dissipation diff --git a/src/source_thermal_externalheat.f90 b/src/constitutive_thermal_externalheat.f90 similarity index 83% rename from src/source_thermal_externalheat.f90 rename to src/constitutive_thermal_externalheat.f90 index 9ba4a051b..2e8c02f8c 100644 --- a/src/source_thermal_externalheat.f90 +++ b/src/constitutive_thermal_externalheat.f90 @@ -4,7 +4,7 @@ !> @author Philip Eisenlohr, Michigan State University !> @brief material subroutine for variable heat source !-------------------------------------------------------------------------------------------------- -submodule(constitutive:constitutive_thermal) source_thermal_externalheat +submodule(constitutive:constitutive_thermal) source_externalheat integer, dimension(:), allocatable :: & @@ -13,7 +13,7 @@ submodule(constitutive:constitutive_thermal) source_thermal_externalheat type :: tParameters !< container type for internal constitutive parameters real(pReal), dimension(:), allocatable :: & - t_n, & + t_n, & f_T integer :: & nIntervals @@ -31,19 +31,20 @@ contains !-------------------------------------------------------------------------------------------------- module function source_thermal_externalheat_init(source_length) result(mySources) - integer, intent(in) :: source_length + integer, intent(in) :: source_length logical, dimension(:,:), allocatable :: mySources class(tNode), pointer :: & phases, & phase, & - sources, & - src + sources, thermal, & + src integer :: Ninstances,sourceOffset,Nconstituents,p - print'(/,a)', ' <<<+- source_thermal_externalHeat init -+>>>' + print'(/,a)', ' <<<+- thermal_externalheat init -+>>>' + + mySources = thermal_active('externalheat',source_length) - mySources = source_active('thermal_externalheat',source_length) Ninstances = count(mySources) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT) if(Ninstances == 0) return @@ -54,15 +55,16 @@ module function source_thermal_externalheat_init(source_length) result(mySources allocate(source_thermal_externalheat_instance(phases%length), source=0) do p = 1, phases%length - phase => phases%get(p) + phase => phases%get(p) if(any(mySources(:,p))) source_thermal_externalheat_instance(p) = count(mySources(:,1:p)) if(count(mySources(:,p)) == 0) cycle - sources => phase%get('source') + thermal => phase%get('thermal') + sources => thermal%get('source') do sourceOffset = 1, sources%length if(mySources(sourceOffset,p)) then source_thermal_externalheat_offset(p) = sourceOffset associate(prm => param(source_thermal_externalheat_instance(p))) - src => sources%get(sourceOffset) + src => sources%get(sourceOffset) prm%t_n = src%get_asFloats('t_n') prm%nIntervals = size(prm%t_n) - 1 @@ -70,9 +72,8 @@ module function source_thermal_externalheat_init(source_length) result(mySources prm%f_T = src%get_asFloats('f_T',requiredSize = size(prm%t_n)) Nconstituents = count(material_phaseAt==p) * discretization_nIPs - call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,0) + call constitutive_allocateState(thermalState(p)%p(sourceOffset),Nconstituents,1,1,0) end associate - endif enddo enddo @@ -95,7 +96,7 @@ module subroutine source_thermal_externalheat_dotState(phase, of) sourceOffset = source_thermal_externalheat_offset(phase) - sourceState(phase)%p(sourceOffset)%dotState(1,of) = 1.0_pReal ! state is current time + thermalState(phase)%p(sourceOffset)%dotState(1,of) = 1.0_pReal ! state is current time end subroutine source_thermal_externalheat_dotState @@ -103,14 +104,13 @@ end subroutine source_thermal_externalheat_dotState !-------------------------------------------------------------------------------------------------- !> @brief returns local heat generation rate !-------------------------------------------------------------------------------------------------- -module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phase, of) +module subroutine thermal_externalheat_getRate(TDot, phase, of) integer, intent(in) :: & phase, & of real(pReal), intent(out) :: & - TDot, & - dTDot_dT + TDot integer :: & sourceOffset, interval @@ -121,7 +121,7 @@ module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_d associate(prm => param(source_thermal_externalheat_instance(phase))) do interval = 1, prm%nIntervals ! scan through all rate segments - frac_time = (sourceState(phase)%p(sourceOffset)%state(1,of) - prm%t_n(interval)) & + frac_time = (thermalState(phase)%p(sourceOffset)%state(1,of) - prm%t_n(interval)) & / (prm%t_n(interval+1) - prm%t_n(interval)) ! fractional time within segment if ( (frac_time < 0.0_pReal .and. interval == 1) & .or. (frac_time >= 1.0_pReal .and. interval == prm%nIntervals) & @@ -130,9 +130,8 @@ module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_d prm%f_T(interval+1) * frac_time ! interpolate heat rate between segment boundaries... ! ...or extrapolate if outside of bounds enddo - dTDot_dT = 0.0 end associate -end subroutine source_thermal_externalheat_getRateAndItsTangent +end subroutine thermal_externalheat_getRate -end submodule source_thermal_externalheat +end submodule source_externalheat diff --git a/src/damage_none.f90 b/src/damage_none.f90 index 3f1144833..078d42af7 100644 --- a/src/damage_none.f90 +++ b/src/damage_none.f90 @@ -25,10 +25,10 @@ subroutine damage_none_init if (damage_type(h) /= DAMAGE_NONE_ID) cycle Nmaterialpoints = count(material_homogenizationAt == h) - damageState(h)%sizeState = 0 - allocate(damageState(h)%state0 (0,Nmaterialpoints)) - allocate(damageState(h)%subState0(0,Nmaterialpoints)) - allocate(damageState(h)%state (0,Nmaterialpoints)) + damageState_h(h)%sizeState = 0 + allocate(damageState_h(h)%state0 (0,Nmaterialpoints)) + allocate(damageState_h(h)%subState0(0,Nmaterialpoints)) + allocate(damageState_h(h)%state (0,Nmaterialpoints)) allocate (damage(h)%p(Nmaterialpoints), source=1.0_pReal) diff --git a/src/damage_nonlocal.f90 b/src/damage_nonlocal.f90 index 3db63cab2..4423c1e3a 100644 --- a/src/damage_nonlocal.f90 +++ b/src/damage_nonlocal.f90 @@ -76,12 +76,12 @@ subroutine damage_nonlocal_init #endif Nmaterialpoints = count(material_homogenizationAt == h) - damageState(h)%sizeState = 1 - allocate(damageState(h)%state0 (1,Nmaterialpoints), source=1.0_pReal) - allocate(damageState(h)%subState0(1,Nmaterialpoints), source=1.0_pReal) - allocate(damageState(h)%state (1,Nmaterialpoints), source=1.0_pReal) + damageState_h(h)%sizeState = 1 + allocate(damageState_h(h)%state0 (1,Nmaterialpoints), source=1.0_pReal) + allocate(damageState_h(h)%subState0(1,Nmaterialpoints), source=1.0_pReal) + allocate(damageState_h(h)%state (1,Nmaterialpoints), source=1.0_pReal) - damage(h)%p => damageState(h)%state(1,:) + damage(h)%p => damageState_h(h)%state(1,:) end associate enddo diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index 259b45f33..9d804ec56 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -256,7 +256,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) PetscObject :: dummy PetscErrorCode :: ierr integer :: i, j, k, cell - real(pReal) :: Tdot, dTdot_dT + real(pReal) :: Tdot T_current = x_scal !-------------------------------------------------------------------------------------------------- @@ -278,7 +278,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr) cell = 0 do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) cell = cell + 1 - call thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T_current(i,j,k), 1, cell) + call thermal_conduction_getSource(Tdot, T_current(i,j,k), 1, cell) scalarField_real(i,j,k) = params%timeinc*(scalarField_real(i,j,k) + Tdot) & + thermal_conduction_getMassDensity (1,cell)* & thermal_conduction_getSpecificHeat(1,cell)*(T_lastInc(i,j,k) - & diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 686bb9885..9112562b9 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -27,6 +27,8 @@ module homogenization !-------------------------------------------------------------------------------------------------- ! General variables for the homogenization at a material point + real(pReal), dimension(:), allocatable, public :: & + homogenization_T real(pReal), dimension(:,:,:), allocatable, public :: & homogenization_F0, & !< def grad of IP at start of FE increment homogenization_F !< def grad of IP to be reached at end of FE increment @@ -56,6 +58,9 @@ module homogenization num_homog !< pointer to mechanical homogenization numerics data end subroutine mech_init + module subroutine thermal_init + end subroutine thermal_init + module subroutine mech_partition(subF,ip,el) real(pReal), intent(in), dimension(3,3) :: & subF @@ -64,6 +69,13 @@ module homogenization el !< element number end subroutine mech_partition + module subroutine thermal_partition(T,ip,el) + real(pReal), intent(in) :: T + integer, intent(in) :: & + ip, & !< integration point + el !< element number + end subroutine thermal_partition + module subroutine mech_homogenize(dt,ip,el) real(pReal), intent(in) :: dt integer, intent(in) :: & @@ -126,9 +138,10 @@ subroutine homogenization_init call mech_init(num_homog) + call thermal_init() - if (any(thermal_type == THERMAL_isothermal_ID)) call thermal_isothermal_init - if (any(thermal_type == THERMAL_conduction_ID)) call thermal_conduction_init + if (any(thermal_type == THERMAL_isothermal_ID)) call thermal_isothermal_init(homogenization_T) + if (any(thermal_type == THERMAL_conduction_ID)) call thermal_conduction_init(homogenization_T) if (any(damage_type == DAMAGE_none_ID)) call damage_none_init if (any(damage_type == DAMAGE_nonlocal_ID)) call damage_nonlocal_init @@ -173,7 +186,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE subStep = 1.0_pReal/num%subStepSizeHomog ! ... larger then the requested calculation if (homogState(ho)%sizeState > 0) homogState(ho)%subState0(:,me) = homogState(ho)%State0(:,me) - if (damageState(ho)%sizeState > 0) damageState(ho)%subState0(:,me) = damageState(ho)%State0(:,me) + if (damageState_h(ho)%sizeState > 0) damageState_h(ho)%subState0(:,me) = damageState_h(ho)%State0(:,me) cutBackLooping: do while (.not. terminallyIll .and. subStep > num%subStepMinHomog) @@ -187,7 +200,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE call constitutive_windForward(ip,el) if(homogState(ho)%sizeState > 0) homogState(ho)%subState0(:,me) = homogState(ho)%State(:,me) - if(damageState(ho)%sizeState > 0) damageState(ho)%subState0(:,me) = damageState(ho)%State(:,me) + if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%subState0(:,me) = damageState_h(ho)%State(:,me) endif steppingNeeded elseif ( (myNgrains == 1 .and. subStep <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite @@ -202,7 +215,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE call constitutive_restore(ip,el,subStep < 1.0_pReal) if(homogState(ho)%sizeState > 0) homogState(ho)%State(:,me) = homogState(ho)%subState0(:,me) - if(damageState(ho)%sizeState > 0) damageState(ho)%State(:,me) = damageState(ho)%subState0(:,me) + if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%State(:,me) = damageState_h(ho)%subState0(:,me) endif if (subStep > num%subStepMinHomog) doneAndHappy = [.false.,.true.] @@ -313,7 +326,7 @@ subroutine homogenization_forward do ho = 1, size(material_name_homogenization) homogState (ho)%state0 = homogState (ho)%state - damageState(ho)%state0 = damageState(ho)%state + damageState_h(ho)%state0 = damageState_h(ho)%state enddo end subroutine homogenization_forward diff --git a/src/homogenization_mech.f90 b/src/homogenization_mech.f90 index 783d08dd1..dd3e47c59 100644 --- a/src/homogenization_mech.f90 +++ b/src/homogenization_mech.f90 @@ -113,24 +113,24 @@ module subroutine mech_partition(subF,ip,el) el !< element number integer :: co - real(pReal) :: F(3,3,homogenization_Nconstituents(material_homogenizationAt(el))) + real(pReal), dimension (3,3,homogenization_Nconstituents(material_homogenizationAt(el))) :: Fs chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) case (HOMOGENIZATION_NONE_ID) chosenHomogenization - F(1:3,1:3,1) = subF + Fs(1:3,1:3,1) = subF case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization - call mech_isostrain_partitionDeformation(F,subF) + call mech_isostrain_partitionDeformation(Fs,subF) case (HOMOGENIZATION_RGC_ID) chosenHomogenization - call mech_RGC_partitionDeformation(F,subF,ip,el) + call mech_RGC_partitionDeformation(Fs,subF,ip,el) end select chosenHomogenization do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) - call constitutive_mech_setF(F(1:3,1:3,co),co,ip,el) + call constitutive_mech_setF(Fs(1:3,1:3,co),co,ip,el) enddo diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90 new file mode 100644 index 000000000..f181d97fa --- /dev/null +++ b/src/homogenization_thermal.f90 @@ -0,0 +1,39 @@ +!-------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, KU Leuven +!-------------------------------------------------------------------------------------------------- +submodule(homogenization) homogenization_thermal + + +contains + +!-------------------------------------------------------------------------------------------------- +!> @brief Allocate variables and set parameters. +!-------------------------------------------------------------------------------------------------- +module subroutine thermal_init() + + + print'(/,a)', ' <<<+- homogenization_thermal init -+>>>' + + allocate(homogenization_T(discretization_nIPs*discretization_Nelems)) + + +end subroutine thermal_init + + +!-------------------------------------------------------------------------------------------------- +!> @brief Partition T onto the individual constituents. +!-------------------------------------------------------------------------------------------------- +module subroutine thermal_partition(T,ip,el) + + real(pReal), intent(in) :: T + integer, intent(in) :: & + ip, & !< integration point + el !< element number + + + call constitutive_thermal_setT(T,1,ip,el) + +end subroutine thermal_partition + + +end submodule homogenization_thermal diff --git a/src/lattice.f90 b/src/lattice.f90 index 6af135e4e..c9b6d99ef 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -453,12 +453,13 @@ contains !-------------------------------------------------------------------------------------------------- subroutine lattice_init - integer :: Nphases, p,i + integer :: Nphases, ph,i class(tNode), pointer :: & phases, & phase, & mech, & - elasticity + elasticity, & + thermal print'(/,a)', ' <<<+- lattice init -+>>>'; flush(IO_STDOUT) @@ -476,67 +477,71 @@ subroutine lattice_init lattice_mu, lattice_nu,& source=[(0.0_pReal,i=1,Nphases)]) - do p = 1, phases%length - phase => phases%get(p) + do ph = 1, phases%length + phase => phases%get(ph) mech => phase%get('mechanics') elasticity => mech%get('elasticity') - lattice_C66(1,1,p) = elasticity%get_asFloat('C_11') - lattice_C66(1,2,p) = elasticity%get_asFloat('C_12') + lattice_C66(1,1,ph) = elasticity%get_asFloat('C_11') + lattice_C66(1,2,ph) = elasticity%get_asFloat('C_12') - lattice_C66(1,3,p) = elasticity%get_asFloat('C_13',defaultVal=0.0_pReal) - lattice_C66(2,2,p) = elasticity%get_asFloat('C_22',defaultVal=0.0_pReal) - lattice_C66(2,3,p) = elasticity%get_asFloat('C_23',defaultVal=0.0_pReal) - lattice_C66(3,3,p) = elasticity%get_asFloat('C_33',defaultVal=0.0_pReal) - lattice_C66(4,4,p) = elasticity%get_asFloat('C_44',defaultVal=0.0_pReal) - lattice_C66(5,5,p) = elasticity%get_asFloat('C_55',defaultVal=0.0_pReal) - lattice_C66(6,6,p) = elasticity%get_asFloat('C_66',defaultVal=0.0_pReal) + lattice_C66(1,3,ph) = elasticity%get_asFloat('C_13',defaultVal=0.0_pReal) + lattice_C66(2,2,ph) = elasticity%get_asFloat('C_22',defaultVal=0.0_pReal) + lattice_C66(2,3,ph) = elasticity%get_asFloat('C_23',defaultVal=0.0_pReal) + lattice_C66(3,3,ph) = elasticity%get_asFloat('C_33',defaultVal=0.0_pReal) + lattice_C66(4,4,ph) = elasticity%get_asFloat('C_44',defaultVal=0.0_pReal) + lattice_C66(5,5,ph) = elasticity%get_asFloat('C_55',defaultVal=0.0_pReal) + lattice_C66(6,6,ph) = elasticity%get_asFloat('C_66',defaultVal=0.0_pReal) select case(phase%get_asString('lattice')) case('cF') - lattice_structure(p) = lattice_FCC_ID + lattice_structure(ph) = lattice_FCC_ID case('cI') - lattice_structure(p) = lattice_BCC_ID + lattice_structure(ph) = lattice_BCC_ID case('hP') - lattice_structure(p) = lattice_HEX_ID + lattice_structure(ph) = lattice_HEX_ID case('tI') - lattice_structure(p) = lattice_BCT_ID + lattice_structure(ph) = lattice_BCT_ID case('oP') - lattice_structure(p) = lattice_ORT_ID + lattice_structure(ph) = lattice_ORT_ID case('aP') - lattice_structure(p) = lattice_ISO_ID + lattice_structure(ph) = lattice_ISO_ID case default call IO_error(130,ext_msg='lattice_init: '//phase%get_asString('lattice')) end select - lattice_C66(1:6,1:6,p) = applyLatticeSymmetryC66(lattice_C66(1:6,1:6,p),phase%get_asString('lattice')) + lattice_C66(1:6,1:6,ph) = applyLatticeSymmetryC66(lattice_C66(1:6,1:6,ph),phase%get_asString('lattice')) - lattice_nu(p) = lattice_equivalent_nu(lattice_C66(1:6,1:6,p),'voigt') - lattice_mu(p) = lattice_equivalent_mu(lattice_C66(1:6,1:6,p),'voigt') + lattice_nu(ph) = lattice_equivalent_nu(lattice_C66(1:6,1:6,ph),'voigt') + lattice_mu(ph) = lattice_equivalent_mu(lattice_C66(1:6,1:6,ph),'voigt') - lattice_C66(1:6,1:6,p) = math_sym3333to66(math_Voigt66to3333(lattice_C66(1:6,1:6,p))) ! Literature data is in Voigt notation + lattice_C66(1:6,1:6,ph) = math_sym3333to66(math_Voigt66to3333(lattice_C66(1:6,1:6,ph))) ! Literature data is in Voigt notation do i = 1, 6 - if (abs(lattice_C66(i,i,p)) phase%get('thermal') + lattice_K(1,1,ph) = thermal%get_asFloat('K_11',defaultVal=0.0_pReal) + lattice_K(2,2,ph) = thermal%get_asFloat('K_22',defaultVal=0.0_pReal) + lattice_K(3,3,ph) = thermal%get_asFloat('K_33',defaultVal=0.0_pReal) + lattice_K(1:3,1:3,ph) = lattice_applyLatticeSymmetry33(lattice_K(1:3,1:3,ph), & + phase%get_asString('lattice')) + lattice_c_p(ph) = thermal%get_asFloat('c_p', defaultVal=0.0_pReal) + endif + + + lattice_D(1,1,ph) = phase%get_asFloat('D_11',defaultVal=0.0_pReal) + lattice_D(2,2,ph) = phase%get_asFloat('D_22',defaultVal=0.0_pReal) + lattice_D(3,3,ph) = phase%get_asFloat('D_33',defaultVal=0.0_pReal) + lattice_D(1:3,1:3,ph) = lattice_applyLatticeSymmetry33(lattice_D(1:3,1:3,ph), & phase%get_asString('lattice')) - lattice_c_p(p) = phase%get_asFloat('c_p', defaultVal=0.0_pReal) - lattice_rho(p) = phase%get_asFloat('rho', defaultVal=0.0_pReal) - - lattice_D(1,1,p) = phase%get_asFloat('D_11',defaultVal=0.0_pReal) - lattice_D(2,2,p) = phase%get_asFloat('D_22',defaultVal=0.0_pReal) - lattice_D(3,3,p) = phase%get_asFloat('D_33',defaultVal=0.0_pReal) - lattice_D(1:3,1:3,p) = lattice_applyLatticeSymmetry33(lattice_D(1:3,1:3,p), & - phase%get_asString('lattice')) - - lattice_M(p) = phase%get_asFloat('M',defaultVal=0.0_pReal) + lattice_M(ph) = phase%get_asFloat('M',defaultVal=0.0_pReal) ! SHOULD NOT BE PART OF LATTICE END call selfTest diff --git a/src/material.f90 b/src/material.f90 index 581182d22..16116ca91 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -61,7 +61,7 @@ module material type(tState), allocatable, dimension(:), public :: & homogState, & - damageState + damageState_h type(Rotation), dimension(:,:,:), allocatable, public, protected :: & material_orientation0 !< initial orientation of each grain,IP,element @@ -101,7 +101,7 @@ subroutine material_init(restart) allocate(homogState (size(material_name_homogenization))) - allocate(damageState (size(material_name_homogenization))) + allocate(damageState_h (size(material_name_homogenization))) allocate(temperature (size(material_name_homogenization))) allocate(damage (size(material_name_homogenization))) diff --git a/src/source_damage_anisoBrittle.f90 b/src/source_damage_anisoBrittle.f90 index 0f923ceba..7c00c6580 100644 --- a/src/source_damage_anisoBrittle.f90 +++ b/src/source_damage_anisoBrittle.f90 @@ -101,9 +101,9 @@ module function source_damage_anisoBrittle_init(source_length) result(mySources) if (any(prm%s_crit < 0.0_pReal)) extmsg = trim(extmsg)//' s_crit' Nconstituents = count(material_phaseAt==p) * discretization_nIPs - call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,0) - sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal) - if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol' + call constitutive_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0) + damageState(p)%p(sourceOffset)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal) + if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol' end associate @@ -146,7 +146,7 @@ module subroutine source_damage_anisoBrittle_dotState(S, co, ip, el) damageOffset = material_homogenizationMemberAt(ip,el) associate(prm => param(source_damage_anisoBrittle_instance(phase))) - sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal + damageState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal do i = 1, prm%sum_N_cl traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i)) traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i)) @@ -154,8 +154,8 @@ module subroutine source_damage_anisoBrittle_dotState(S, co, ip, el) traction_crit = prm%g_crit(i)*damage(homog)%p(damageOffset)**2.0_pReal - sourceState(phase)%p(sourceOffset)%dotState(1,constituent) & - = sourceState(phase)%p(sourceOffset)%dotState(1,constituent) & + damageState(phase)%p(sourceOffset)%dotState(1,constituent) & + = damageState(phase)%p(sourceOffset)%dotState(1,constituent) & + prm%dot_o / prm%s_crit(i) & * ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%q + & (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%q + & @@ -185,7 +185,7 @@ module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, d sourceOffset = source_damage_anisoBrittle_offset(phase) - dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent) + dLocalphiDot_dPhi = -damageState(phase)%p(sourceOffset)%state(1,constituent) localphiDot = 1.0_pReal & + dLocalphiDot_dPhi*phi @@ -204,7 +204,7 @@ module subroutine source_damage_anisoBrittle_results(phase,group) integer :: o associate(prm => param(source_damage_anisoBrittle_instance(phase)), & - stt => sourceState(phase)%p(source_damage_anisoBrittle_offset(phase))%state) + stt => damageState(phase)%p(source_damage_anisoBrittle_offset(phase))%state) outputsLoop: do o = 1,size(prm%output) select case(trim(prm%output(o))) case ('f_phi') diff --git a/src/source_damage_anisoDuctile.f90 b/src/source_damage_anisoDuctile.f90 index 6f71fc145..7ec06cb62 100644 --- a/src/source_damage_anisoDuctile.f90 +++ b/src/source_damage_anisoDuctile.f90 @@ -87,9 +87,9 @@ module function source_damage_anisoDuctile_init(source_length) result(mySources) if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit' Nconstituents=count(material_phaseAt==p) * discretization_nIPs - call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,0) - sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal) - if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol' + call constitutive_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0) + damageState(p)%p(sourceOffset)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal) + if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol' end associate @@ -128,7 +128,7 @@ module subroutine source_damage_anisoDuctile_dotState(co, ip, el) damageOffset = material_homogenizationMemberAt(ip,el) associate(prm => param(source_damage_anisoDuctile_instance(phase))) - sourceState(phase)%p(sourceOffset)%dotState(1,constituent) & + damageState(phase)%p(sourceOffset)%dotState(1,constituent) & = sum(plasticState(phase)%slipRate(:,constituent)/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit) end associate @@ -154,7 +154,7 @@ module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, d sourceOffset = source_damage_anisoDuctile_offset(phase) - dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent) + dLocalphiDot_dPhi = -damageState(phase)%p(sourceOffset)%state(1,constituent) localphiDot = 1.0_pReal & + dLocalphiDot_dPhi*phi @@ -173,7 +173,7 @@ module subroutine source_damage_anisoDuctile_results(phase,group) integer :: o associate(prm => param(source_damage_anisoDuctile_instance(phase)), & - stt => sourceState(phase)%p(source_damage_anisoDuctile_offset(phase))%state) + stt => damageState(phase)%p(source_damage_anisoDuctile_offset(phase))%state) outputsLoop: do o = 1,size(prm%output) select case(trim(prm%output(o))) case ('f_phi') diff --git a/src/source_damage_isoBrittle.f90 b/src/source_damage_isoBrittle.f90 index 8c768b08d..1721b0201 100644 --- a/src/source_damage_isoBrittle.f90 +++ b/src/source_damage_isoBrittle.f90 @@ -74,9 +74,9 @@ module function source_damage_isoBrittle_init(source_length) result(mySources) if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit' Nconstituents = count(material_phaseAt==p) * discretization_nIPs - call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,1) - sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal) - if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' + call constitutive_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,1) + damageState(p)%p(sourceOffset)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal) + if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' end associate @@ -124,13 +124,13 @@ module subroutine source_damage_isoBrittle_deltaState(C, Fe, co, ip, el) strainenergy = 2.0_pReal*sum(strain*matmul(C,strain))/prm%W_crit ! ToDo: check strainenergy = 2.0_pReal*dot_product(strain,matmul(C,strain))/prm%W_crit - if (strainenergy > sourceState(phase)%p(sourceOffset)%subState0(1,constituent)) then - sourceState(phase)%p(sourceOffset)%deltaState(1,constituent) = & - strainenergy - sourceState(phase)%p(sourceOffset)%state(1,constituent) + if (strainenergy > damageState(phase)%p(sourceOffset)%subState0(1,constituent)) then + damageState(phase)%p(sourceOffset)%deltaState(1,constituent) = & + strainenergy - damageState(phase)%p(sourceOffset)%state(1,constituent) else - sourceState(phase)%p(sourceOffset)%deltaState(1,constituent) = & - sourceState(phase)%p(sourceOffset)%subState0(1,constituent) - & - sourceState(phase)%p(sourceOffset)%state(1,constituent) + damageState(phase)%p(sourceOffset)%deltaState(1,constituent) = & + damageState(phase)%p(sourceOffset)%subState0(1,constituent) - & + damageState(phase)%p(sourceOffset)%state(1,constituent) endif end associate @@ -158,8 +158,8 @@ module subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLo associate(prm => param(source_damage_isoBrittle_instance(phase))) localphiDot = 1.0_pReal & - - phi*sourceState(phase)%p(sourceOffset)%state(1,constituent) - dLocalphiDot_dPhi = - sourceState(phase)%p(sourceOffset)%state(1,constituent) + - phi*damageState(phase)%p(sourceOffset)%state(1,constituent) + dLocalphiDot_dPhi = - damageState(phase)%p(sourceOffset)%state(1,constituent) end associate end subroutine source_damage_isoBrittle_getRateAndItsTangent @@ -176,7 +176,7 @@ module subroutine source_damage_isoBrittle_results(phase,group) integer :: o associate(prm => param(source_damage_isoBrittle_instance(phase)), & - stt => sourceState(phase)%p(source_damage_isoBrittle_offset(phase))%state) + stt => damageState(phase)%p(source_damage_isoBrittle_offset(phase))%state) outputsLoop: do o = 1,size(prm%output) select case(trim(prm%output(o))) case ('f_phi') diff --git a/src/source_damage_isoDuctile.f90 b/src/source_damage_isoDuctile.f90 index 86222bbf9..dd2910182 100644 --- a/src/source_damage_isoDuctile.f90 +++ b/src/source_damage_isoDuctile.f90 @@ -78,9 +78,9 @@ module function source_damage_isoDuctile_init(source_length) result(mySources) if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit' Nconstituents=count(material_phaseAt==p) * discretization_nIPs - call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,0) - sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal) - if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol' + call constitutive_allocateState(damageState(p)%p(sourceOffset),Nconstituents,1,1,0) + damageState(p)%p(sourceOffset)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal) + if(any(damageState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol' end associate @@ -119,7 +119,7 @@ module subroutine source_damage_isoDuctile_dotState(co, ip, el) damageOffset = material_homogenizationMemberAt(ip,el) associate(prm => param(source_damage_isoDuctile_instance(phase))) - sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = & + damageState(phase)%p(sourceOffset)%dotState(1,constituent) = & sum(plasticState(phase)%slipRate(:,constituent))/(damage(homog)%p(damageOffset)**prm%q)/prm%gamma_crit end associate @@ -145,7 +145,7 @@ module subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLo sourceOffset = source_damage_isoDuctile_offset(phase) - dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent) + dLocalphiDot_dPhi = -damageState(phase)%p(sourceOffset)%state(1,constituent) localphiDot = 1.0_pReal & + dLocalphiDot_dPhi*phi @@ -164,7 +164,7 @@ module subroutine source_damage_isoDuctile_results(phase,group) integer :: o associate(prm => param(source_damage_isoDuctile_instance(phase)), & - stt => sourceState(phase)%p(source_damage_isoDuctile_offset(phase))%state) + stt => damageState(phase)%p(source_damage_isoDuctile_offset(phase))%state) outputsLoop: do o = 1,size(prm%output) select case(trim(prm%output(o))) case ('f_phi') diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index f98d36d3b..79fe0d6cd 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -10,6 +10,7 @@ module thermal_conduction use results use constitutive use YAML_types + use discretization implicit none private @@ -24,7 +25,7 @@ module thermal_conduction public :: & thermal_conduction_init, & - thermal_conduction_getSourceAndItsTangent, & + thermal_conduction_getSource, & thermal_conduction_getConductivity, & thermal_conduction_getSpecificHeat, & thermal_conduction_getMassDensity, & @@ -38,25 +39,28 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -subroutine thermal_conduction_init +subroutine thermal_conduction_init(T) - integer :: Ninstances,Nmaterialpoints,h + real(pReal), dimension(:), intent(inout) :: T + + integer :: Ninstances,Nmaterialpoints,ho,ip,el,ce class(tNode), pointer :: & material_homogenization, & homog, & homogThermal - + + print'(/,a)', ' <<<+- thermal_conduction init -+>>>'; flush(6) Ninstances = count(thermal_type == THERMAL_conduction_ID) allocate(param(Ninstances)) material_homogenization => config_material%get('homogenization') - do h = 1, size(material_name_homogenization) - if (thermal_type(h) /= THERMAL_conduction_ID) cycle - homog => material_homogenization%get(h) + do ho = 1, size(material_name_homogenization) + if (thermal_type(ho) /= THERMAL_conduction_ID) cycle + homog => material_homogenization%get(ho) homogThermal => homog%get('thermal') - associate(prm => param(thermal_typeInstance(h))) + associate(prm => param(thermal_typeInstance(ho))) #if defined (__GFORTRAN__) prm%output = output_asStrings(homogThermal) @@ -64,21 +68,30 @@ subroutine thermal_conduction_init prm%output = homogThermal%get_asStrings('output',defaultVal=emptyStringArray) #endif - Nmaterialpoints=count(material_homogenizationAt==h) + Nmaterialpoints=count(material_homogenizationAt==ho) - allocate (temperature (h)%p(Nmaterialpoints), source=thermal_initialT(h)) - allocate (temperatureRate(h)%p(Nmaterialpoints), source=0.0_pReal) + allocate (temperature (ho)%p(Nmaterialpoints), source=thermal_initialT(ho)) + allocate (temperatureRate(ho)%p(Nmaterialpoints), source=0.0_pReal) end associate enddo + ce = 0 + do el = 1, discretization_Nelems + do ip = 1, discretization_nIPs + ce = ce + 1 + ho = material_homogenizationAt(el) + if (thermal_type(ho) == THERMAL_conduction_ID) T(ce) = thermal_initialT(ho) + enddo + enddo + end subroutine thermal_conduction_init !-------------------------------------------------------------------------------------------------- !> @brief return heat generation rate !-------------------------------------------------------------------------------------------------- -subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) +subroutine thermal_conduction_getSource(Tdot, T,ip,el) integer, intent(in) :: & ip, & !< integration point number @@ -86,20 +99,17 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el) real(pReal), intent(in) :: & T real(pReal), intent(out) :: & - Tdot, dTdot_dT - integer :: & - homog - - Tdot = 0.0_pReal - dTdot_dT = 0.0_pReal + Tdot - homog = material_homogenizationAt(el) - call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, ip, el) + integer :: & + homog + + homog = material_homogenizationAt(el) + call constitutive_thermal_getRate(TDot, T,ip,el) Tdot = Tdot/real(homogenization_Nconstituents(homog),pReal) - dTdot_dT = dTdot_dT/real(homogenization_Nconstituents(homog),pReal) -end subroutine thermal_conduction_getSourceAndItsTangent +end subroutine thermal_conduction_getSource !-------------------------------------------------------------------------------------------------- @@ -112,13 +122,13 @@ function thermal_conduction_getConductivity(ip,el) el !< element number real(pReal), dimension(3,3) :: & thermal_conduction_getConductivity - + integer :: & co thermal_conduction_getConductivity = 0.0_pReal - + do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) thermal_conduction_getConductivity = thermal_conduction_getConductivity + & crystallite_push33ToRef(co,ip,el,lattice_K(:,:,material_phaseAt(co,el))) @@ -168,7 +178,7 @@ function thermal_conduction_getMassDensity(ip,el) el !< element number real(pReal) :: & thermal_conduction_getMassDensity - + integer :: & co diff --git a/src/thermal_isothermal.f90 b/src/thermal_isothermal.f90 index 2a41ada49..09e35931e 100644 --- a/src/thermal_isothermal.f90 +++ b/src/thermal_isothermal.f90 @@ -6,6 +6,7 @@ module thermal_isothermal use prec use config use material + use discretization implicit none public @@ -15,22 +16,33 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief allocates fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- -subroutine thermal_isothermal_init +subroutine thermal_isothermal_init(T) - integer :: h,Nmaterialpoints + real(pReal), dimension(:), intent(inout) :: T + + integer :: Ninstances,Nmaterialpoints,ho,ip,el,ce print'(/,a)', ' <<<+- thermal_isothermal init -+>>>'; flush(6) - do h = 1, size(material_name_homogenization) - if (thermal_type(h) /= THERMAL_isothermal_ID) cycle + do ho = 1, size(thermal_type) + if (thermal_type(ho) /= THERMAL_isothermal_ID) cycle - Nmaterialpoints = count(material_homogenizationAt == h) + Nmaterialpoints = count(material_homogenizationAt == ho) - allocate(temperature (h)%p(Nmaterialpoints),source=thermal_initialT(h)) - allocate(temperatureRate(h)%p(Nmaterialpoints),source = 0.0_pReal) + allocate(temperature (ho)%p(Nmaterialpoints),source=thermal_initialT(ho)) + allocate(temperatureRate(ho)%p(Nmaterialpoints),source = 0.0_pReal) enddo + ce = 0 + do el = 1, discretization_Nelems + do ip = 1, discretization_nIPs + ce = ce + 1 + ho = material_homogenizationAt(el) + if (thermal_type(ho) == THERMAL_isothermal_ID) T(ce) = thermal_initialT(ho) + enddo + enddo + end subroutine thermal_isothermal_init end module thermal_isothermal