From ef792a578ba8f9e5f1f26a03f6633a1fe87e3966 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Tue, 16 Mar 2021 17:20:24 +0100 Subject: [PATCH] separate elastic submodule --- src/commercialFEM_fileList.f90 | 1 + src/phase.f90 | 4 - src/phase_mechanical.f90 | 112 ++++++---------------------- src/phase_mechanical_elastic.f90 | 123 +++++++++++++++++++++++++++++++ 4 files changed, 146 insertions(+), 94 deletions(-) create mode 100644 src/phase_mechanical_elastic.f90 diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index 97e7520f9..ff7c67f15 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -21,6 +21,7 @@ #include "lattice.f90" #include "phase.f90" #include "phase_mechanical.f90" +#include "phase_mechanical_elastic.f90" #include "phase_mechanical_plastic.f90" #include "phase_mechanical_plastic_none.f90" #include "phase_mechanical_plastic_isotropic.f90" diff --git a/src/phase.f90 b/src/phase.f90 index add26ae0b..64245e1d1 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -60,10 +60,6 @@ module phase type(tDebugOptions) :: debugCrystallite - integer, dimension(:), allocatable, public :: & !< ToDo: should be protected (bug in Intel compiler) - phase_elasticityInstance, & - phase_NstiffnessDegradations - logical, dimension(:), allocatable, public :: & ! ToDo: should be protected (bug in Intel Compiler) phase_localPlasticity !< flags phases with local constitutive law diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index a7e0c3254..7d8f6e463 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -5,10 +5,6 @@ submodule(phase) mechanical enum, bind(c); enumerator :: & - ELASTICITY_UNDEFINED_ID, & - ELASTICITY_HOOKE_ID, & - STIFFNESS_DEGRADATION_UNDEFINED_ID, & - STIFFNESS_DEGRADATION_DAMAGE_ID, & PLASTICITY_UNDEFINED_ID, & PLASTICITY_NONE_ID, & PLASTICITY_ISOTROPIC_ID, & @@ -23,11 +19,6 @@ submodule(phase) mechanical KINEMATICS_THERMAL_EXPANSION_ID end enum - integer(kind(ELASTICITY_UNDEFINED_ID)), dimension(:), allocatable :: & - phase_elasticity !< elasticity of each phase - integer(kind(STIFFNESS_DEGRADATION_UNDEFINED_ID)), dimension(:,:), allocatable :: & - phase_stiffnessDegradation !< active stiffness degradation mechanisms of each phase - type(tTensorContainer), dimension(:), allocatable :: & ! current value phase_mechanical_Fe, & @@ -57,9 +48,27 @@ submodule(phase) mechanical class(tNode), pointer :: phases end subroutine eigendeformation_init + module subroutine elastic_init(phases) + class(tNode), pointer :: phases + end subroutine elastic_init + module subroutine plastic_init end subroutine plastic_init + module subroutine phase_hooke_SandItsTangents(S,dS_dFe,dS_dFi,Fe,Fi,ph,me) + integer, intent(in) :: & + ph, & + me + real(pReal), intent(in), dimension(3,3) :: & + Fe, & !< elastic deformation gradient + Fi !< intermediate deformation gradient + real(pReal), intent(out), dimension(3,3) :: & + S !< 2nd Piola-Kirchhoff stress tensor in lattice configuration + real(pReal), intent(out), dimension(3,3,3,3) :: & + dS_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient + dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient + end subroutine phase_hooke_SandItsTangents + module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,me) real(pReal), dimension(3,3), intent(out) :: & Li !< inleastic velocity gradient @@ -73,7 +82,6 @@ submodule(phase) mechanical end subroutine plastic_isotropic_LiAndItsTangent module function plastic_dotState(subdt,co,ip,el,ph,me) result(broken) - integer, intent(in) :: & co, & !< component-ID of integration point ip, & !< integration point @@ -190,7 +198,6 @@ module subroutine mechanical_init(materials,phases) ce, & ph, & me, & - stiffDegradationCtr, & Nmembers class(tNode), pointer :: & num_crystallite, & @@ -198,17 +205,11 @@ module subroutine mechanical_init(materials,phases) constituents, & constituent, & phase, & - mech, & - elastic, & - stiffDegradation + mech print'(/,a)', ' <<<+- phase:mechanical init -+>>>' !------------------------------------------------------------------------------------------------- -! initialize elasticity (hooke) !ToDO: Maybe move to elastic submodule along with function homogenizedC? - allocate(phase_elasticity(phases%length), source = ELASTICITY_undefined_ID) - allocate(phase_elasticityInstance(phases%length), source = 0) - allocate(phase_NstiffnessDegradations(phases%length),source=0) allocate(output_constituent(phases%length)) allocate(phase_mechanical_Fe(phases%length)) @@ -253,32 +254,8 @@ module subroutine mechanical_init(materials,phases) #else output_constituent(ph)%label = mech%get_asStrings('output',defaultVal=emptyStringArray) #endif - elastic => mech%get('elastic') - if(elastic%get_asString('type') == 'hooke') then - phase_elasticity(ph) = ELASTICITY_HOOKE_ID - else - call IO_error(200,ext_msg=elastic%get_asString('type')) - endif - stiffDegradation => mech%get('stiffness_degradation',defaultVal=emptyList) ! check for stiffness degradation mechanisms - phase_NstiffnessDegradations(ph) = stiffDegradation%length enddo - allocate(phase_stiffnessDegradation(maxval(phase_NstiffnessDegradations),phases%length), & - source=STIFFNESS_DEGRADATION_undefined_ID) - - if(maxVal(phase_NstiffnessDegradations)/=0) then - do ph = 1, phases%length - phase => phases%get(ph) - mech => phase%get('mechanical') - stiffDegradation => mech%get('stiffness_degradation',defaultVal=emptyList) - do stiffDegradationCtr = 1, stiffDegradation%length - if(stiffDegradation%get_asString(stiffDegradationCtr) == 'damage') & - phase_stiffnessDegradation(stiffDegradationCtr,ph) = STIFFNESS_DEGRADATION_damage_ID - enddo - enddo - endif - - do el = 1, size(material_phaseMemberAt,3); do ip = 1, size(material_phaseMemberAt,2) do co = 1, homogenization_Nconstituents(material_homogenizationAt(el)) material => materials%get(discretization_materialAt(el)) @@ -306,6 +283,9 @@ module subroutine mechanical_init(materials,phases) enddo; enddo +! initialize elasticity + call elastic_init(phases) + ! initialize plasticity allocate(plasticState(phases%length)) allocate(phase_plasticity(phases%length),source = PLASTICITY_undefined_ID) @@ -313,9 +293,6 @@ module subroutine mechanical_init(materials,phases) call plastic_init() - do ph = 1, phases%length - phase_elasticityInstance(ph) = count(phase_elasticity(1:ph) == phase_elasticity(ph)) - enddo num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict) @@ -348,51 +325,6 @@ module subroutine mechanical_init(materials,phases) end subroutine mechanical_init -!-------------------------------------------------------------------------------------------------- -!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to -!> the elastic and intermediate deformation gradients using Hooke's law -!-------------------------------------------------------------------------------------------------- -subroutine phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & - Fe, Fi, ph, me) - - integer, intent(in) :: & - ph, & - me - real(pReal), intent(in), dimension(3,3) :: & - Fe, & !< elastic deformation gradient - Fi !< intermediate deformation gradient - real(pReal), intent(out), dimension(3,3) :: & - S !< 2nd Piola-Kirchhoff stress tensor in lattice configuration - real(pReal), intent(out), dimension(3,3,3,3) :: & - dS_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient - dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient - - real(pReal), dimension(3,3) :: E - real(pReal), dimension(3,3,3,3) :: C - integer :: & - d, & !< counter in degradation loop - i, j - - C = math_66toSym3333(phase_homogenizedC(ph,me)) - - DegradationLoop: do d = 1, phase_NstiffnessDegradations(ph) - degradationType: select case(phase_stiffnessDegradation(d,ph)) - case (STIFFNESS_DEGRADATION_damage_ID) degradationType - C = C * damage_phi(ph,me)**2 - end select degradationType - enddo DegradationLoop - - E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration - S = math_mul3333xx33(C,matmul(matmul(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration - - do i =1, 3;do j=1,3 - dS_dFe(i,j,1:3,1:3) = matmul(Fe,matmul(matmul(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko - dS_dFi(i,j,1:3,1:3) = 2.0_pReal*matmul(matmul(E,Fi),C(i,j,1:3,1:3)) !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn - enddo; enddo - -end subroutine phase_hooke_SandItsTangents - - module subroutine mechanical_results(group,ph) character(len=*), intent(in) :: group diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 new file mode 100644 index 000000000..2bd4ed37b --- /dev/null +++ b/src/phase_mechanical_elastic.f90 @@ -0,0 +1,123 @@ +submodule(phase:mechanical) elastic + + enum, bind(c); enumerator :: & + ELASTICITY_UNDEFINED_ID, & + ELASTICITY_HOOKE_ID, & + STIFFNESS_DEGRADATION_UNDEFINED_ID, & + STIFFNESS_DEGRADATION_DAMAGE_ID + end enum + + integer, dimension(:), allocatable :: & !< ToDo: should be protected (bug in Intel compiler) + phase_NstiffnessDegradations + + integer(kind(ELASTICITY_UNDEFINED_ID)), dimension(:), allocatable :: & + phase_elasticity !< elasticity of each phase + integer(kind(STIFFNESS_DEGRADATION_UNDEFINED_ID)), dimension(:,:), allocatable :: & + phase_stiffnessDegradation !< active stiffness degradation mechanisms of each phase + +contains + +module subroutine elastic_init(phases) + + class(tNode), pointer :: & + phases + + integer :: & + ph, & + stiffDegradationCtr + class(tNode), pointer :: & + phase, & + mech, & + elastic, & + stiffDegradation + + + print'(/,a)', ' <<<+- phase:mechanical:elastic init -+>>>' + +! initialize elasticity (hooke) + allocate(phase_elasticity(phases%length), source = ELASTICITY_undefined_ID) + allocate(phase_NstiffnessDegradations(phases%length),source=0) + + do ph = 1, phases%length + phase => phases%get(ph) + mech => phase%get('mechanical') + elastic => mech%get('elastic') + if(elastic%get_asString('type') == 'hooke') then + phase_elasticity(ph) = ELASTICITY_HOOKE_ID + else + call IO_error(200,ext_msg=elastic%get_asString('type')) + endif + stiffDegradation => mech%get('stiffness_degradation',defaultVal=emptyList) ! check for stiffness degradation mechanisms + phase_NstiffnessDegradations(ph) = stiffDegradation%length + enddo + + allocate(phase_stiffnessDegradation(maxval(phase_NstiffnessDegradations),phases%length), & + source=STIFFNESS_DEGRADATION_undefined_ID) + + if(maxVal(phase_NstiffnessDegradations)/=0) then + do ph = 1, phases%length + phase => phases%get(ph) + mech => phase%get('mechanical') + stiffDegradation => mech%get('stiffness_degradation',defaultVal=emptyList) + do stiffDegradationCtr = 1, stiffDegradation%length + if(stiffDegradation%get_asString(stiffDegradationCtr) == 'damage') & + phase_stiffnessDegradation(stiffDegradationCtr,ph) = STIFFNESS_DEGRADATION_damage_ID + enddo + enddo + endif + +end subroutine elastic_init + + +!-------------------------------------------------------------------------------------------------- +!> @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 phase_hooke_SandItsTangents(S, dS_dFe, dS_dFi, & + Fe, Fi, ph, me) + + integer, intent(in) :: & + ph, & + me + real(pReal), intent(in), dimension(3,3) :: & + Fe, & !< elastic deformation gradient + Fi !< intermediate deformation gradient + real(pReal), intent(out), dimension(3,3) :: & + S !< 2nd Piola-Kirchhoff stress tensor in lattice configuration + real(pReal), intent(out), dimension(3,3,3,3) :: & + dS_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient + dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient + + real(pReal), dimension(3,3) :: E + real(pReal), dimension(3,3,3,3) :: C + integer :: & + d, & !< counter in degradation loop + i, j + + C = math_66toSym3333(phase_homogenizedC(ph,me)) + + DegradationLoop: do d = 1, phase_NstiffnessDegradations(ph) + degradationType: select case(phase_stiffnessDegradation(d,ph)) + case (STIFFNESS_DEGRADATION_damage_ID) degradationType + C = C * damage_phi(ph,me)**2 + end select degradationType + enddo DegradationLoop + + E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration + S = math_mul3333xx33(C,matmul(matmul(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration + + do i =1, 3;do j=1,3 + dS_dFe(i,j,1:3,1:3) = matmul(Fe,matmul(matmul(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko + dS_dFi(i,j,1:3,1:3) = 2.0_pReal*matmul(matmul(E,Fi),C(i,j,1:3,1:3)) !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn + enddo; enddo + +end subroutine phase_hooke_SandItsTangents + + +end submodule elastic + + + + + +