223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

Merge branch 'development' into marc-python-modifications

This commit is contained in:
Martin Diehl 2020-07-22 14:10:56 +02:00
parent 888e7028f9 83504b417d
commit ed16cb087e
29 changed files with 1049 additions and 747 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
VERSION
View File

@ -1 +1 @@
v2.0.3-2881-gc07efe84 v2.0.3-2929-gfb02c69a

21
src/commercialFEM_fileList.f90
View File

@ -26,16 +26,8 @@
#endif #endif
#include "material.f90" #include "material.f90"
#include "lattice.f90" #include "lattice.f90"
#include "source_thermal_dissipation.f90"
#include "source_thermal_externalheat.f90"
#include "source_damage_isoBrittle.f90"
#include "source_damage_isoDuctile.f90"
#include "source_damage_anisoBrittle.f90"
#include "source_damage_anisoDuctile.f90"
#include "kinematics_cleavage_opening.f90"
#include "kinematics_slipplane_opening.f90"
#include "kinematics_thermal_expansion.f90"
#include "constitutive.f90" #include "constitutive.f90"
#include "constitutive_plastic.f90"
#include "constitutive_plastic_none.f90" #include "constitutive_plastic_none.f90"
#include "constitutive_plastic_isotropic.f90" #include "constitutive_plastic_isotropic.f90"
#include "constitutive_plastic_phenopowerlaw.f90" #include "constitutive_plastic_phenopowerlaw.f90"
@ -43,6 +35,17 @@
#include "constitutive_plastic_dislotwin.f90" #include "constitutive_plastic_dislotwin.f90"
#include "constitutive_plastic_disloUCLA.f90" #include "constitutive_plastic_disloUCLA.f90"
#include "constitutive_plastic_nonlocal.f90" #include "constitutive_plastic_nonlocal.f90"
#include "constitutive_thermal.f90"
#include "source_thermal_dissipation.f90"
#include "source_thermal_externalheat.f90"
#include "kinematics_thermal_expansion.f90"
#include "constitutive_damage.f90"
#include "source_damage_isoBrittle.f90"
#include "source_damage_isoDuctile.f90"
#include "source_damage_anisoBrittle.f90"
#include "source_damage_anisoDuctile.f90"
#include "kinematics_cleavage_opening.f90"
#include "kinematics_slipplane_opening.f90"
#include "crystallite.f90" #include "crystallite.f90"
#include "thermal_isothermal.f90" #include "thermal_isothermal.f90"
#include "thermal_adiabatic.f90" #include "thermal_adiabatic.f90"

564
src/constitutive.f90
View File

@ -1,7 +1,7 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief elasticity, plasticity, internal microstructure state !> @brief elasticity, plasticity, damage & thermal internal microstructure state
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module constitutive module constitutive
use prec use prec
@ -17,15 +17,6 @@ module constitutive
use discretization use discretization
use geometry_plastic_nonlocal, only: & use geometry_plastic_nonlocal, only: &
geometry_plastic_nonlocal_disable geometry_plastic_nonlocal_disable
use source_thermal_dissipation
use source_thermal_externalheat
use source_damage_isoBrittle
use source_damage_isoDuctile
use source_damage_anisoBrittle
use source_damage_anisoDuctile
use kinematics_cleavage_opening
use kinematics_slipplane_opening
use kinematics_thermal_expansion
implicit none implicit none
private private
@ -36,128 +27,14 @@ module constitutive
interface interface
module subroutine plastic_none_init module subroutine plastic_init
end subroutine plastic_none_init end subroutine plastic_init
module subroutine plastic_isotropic_init module subroutine damage_init
end subroutine plastic_isotropic_init end subroutine damage_init
module subroutine plastic_phenopowerlaw_init module subroutine thermal_init
end subroutine plastic_phenopowerlaw_init end subroutine thermal_init
module subroutine plastic_kinehardening_init
end subroutine plastic_kinehardening_init
module subroutine plastic_dislotwin_init
end subroutine plastic_dislotwin_init
module subroutine plastic_disloUCLA_init
end subroutine plastic_disloUCLA_init
module subroutine plastic_nonlocal_init
end subroutine plastic_nonlocal_init
module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
end subroutine plastic_isotropic_LpAndItsTangent
pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
end subroutine plastic_phenopowerlaw_LpAndItsTangent
pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
end subroutine plastic_kinehardening_LpAndItsTangent
module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of)
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
real(pReal), intent(in) :: &
T
integer, intent(in) :: &
instance, &
of
end subroutine plastic_dislotwin_LpAndItsTangent
pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of)
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
real(pReal), intent(in) :: &
T
integer, intent(in) :: &
instance, &
of
end subroutine plastic_disloUCLA_LpAndItsTangent
module subroutine plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
Mp,Temperature,instance,of,ip,el)
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
real(pReal), intent(in) :: &
Temperature
integer, intent(in) :: &
instance, &
of, &
ip, & !< current integration point
el !< current element number
end subroutine plastic_nonlocal_LpAndItsTangent
module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of)
real(pReal), dimension(3,3), intent(out) :: &
Li !< inleastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLi_dMi !< derivative of Li with respect to Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mi !< Mandel stress
integer, intent(in) :: &
instance, &
of
end subroutine plastic_isotropic_LiAndItsTangent
module subroutine plastic_isotropic_dotState(Mp,instance,of) module subroutine plastic_isotropic_dotState(Mp,instance,of)
@ -205,8 +82,8 @@ module constitutive
end subroutine plastic_disloUCLA_dotState end subroutine plastic_disloUCLA_dotState
module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, & module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, &
instance,of,ip,el) instance,of,ip,el)
real(pReal), dimension(3,3), intent(in) ::& real(pReal), dimension(3,3), intent(in) :: &
Mp !< MandelStress Mp !< MandelStress
real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: & real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: &
F, & !< deformation gradient F, & !< deformation gradient
@ -222,30 +99,135 @@ module constitutive
end subroutine plastic_nonlocal_dotState end subroutine plastic_nonlocal_dotState
module subroutine plastic_dislotwin_dependentState(T,instance,of) module subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
instance, & ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(3,3) :: &
S
end subroutine source_damage_anisoBrittle_dotState
module subroutine source_damage_anisoDuctile_dotState(ipc, ip, el)
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
end subroutine source_damage_anisoDuctile_dotState
module subroutine source_damage_isoDuctile_dotState(ipc, ip, el)
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
end subroutine source_damage_isoDuctile_dotState
module subroutine source_thermal_externalheat_dotState(phase, of)
integer, intent(in) :: &
phase, &
of of
real(pReal), intent(in) :: & end subroutine source_thermal_externalheat_dotState
module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: &
phi !< damage parameter
real(pReal), intent(inout) :: &
phiDot, &
dPhiDot_dPhi
end subroutine constitutive_damage_getRateAndItsTangents
module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, S, Lp, ip, el)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: &
T T
end subroutine plastic_dislotwin_dependentState real(pReal), intent(in), dimension(:,:,:,:,:) :: &
S, & !< current 2nd Piola Kitchoff stress vector
Lp !< plastic velocity gradient
real(pReal), intent(inout) :: &
TDot, &
dTDot_dT
end subroutine constitutive_thermal_getRateAndItsTangents
module subroutine plastic_disloUCLA_dependentState(instance,of) module function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC)
integer, intent(in) :: & real(pReal), dimension(6,6) :: &
instance, & homogenizedC
of integer, intent(in) :: &
end subroutine plastic_disloUCLA_dependentState ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
end function plastic_dislotwin_homogenizedC
module subroutine plastic_nonlocal_dependentState(F, Fp, instance, of, ip, el) pure module function kinematics_thermal_expansion_initialStrain(homog,phase,offset) result(initialStrain)
real(pReal), dimension(3,3), intent(in) :: & integer, intent(in) :: &
F, & phase, &
Fp homog, &
offset
real(pReal), dimension(3,3) :: &
initialStrain
end function kinematics_thermal_expansion_initialStrain
module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e)
integer, intent(in) :: & integer, intent(in) :: &
instance, & instance, &
of, & i, &
ip, & e
el type(rotation), dimension(1,discretization_nIP,discretization_nElem), intent(in) :: &
end subroutine plastic_nonlocal_dependentState 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
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLi_dMi !< derivative of Li with respect to Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mi !< Mandel stress
integer, intent(in) :: &
instance, &
of
end subroutine plastic_isotropic_LiAndItsTangent
module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el)
integer, intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), intent(in), dimension(3,3) :: &
S
real(pReal), intent(out), dimension(3,3) :: &
Ld !< damage velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor)
end subroutine kinematics_cleavage_opening_LiAndItsTangent
module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el)
integer, intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), intent(in), dimension(3,3) :: &
S
real(pReal), intent(out), dimension(3,3) :: &
Ld !< damage velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor)
end subroutine kinematics_slipplane_opening_LiAndItsTangent
module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, el)
integer, intent(in) :: &
ipc, & !< grain number
ip, & !< integration point number
el !< element number
real(pReal), intent(out), dimension(3,3) :: &
Li !< thermal velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLi_dTstar !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero)
end subroutine kinematics_thermal_expansion_LiAndItsTangent
module subroutine plastic_kinehardening_deltaState(Mp,instance,of) module subroutine plastic_kinehardening_deltaState(Mp,instance,of)
@ -266,58 +248,60 @@ module constitutive
el el
end subroutine plastic_nonlocal_deltaState end subroutine plastic_nonlocal_deltaState
module subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el)
module function plastic_dislotwin_homogenizedC(ipc,ip,el) result(homogenizedC) integer, intent(in) :: &
real(pReal), dimension(6,6) :: &
homogenizedC
integer, intent(in) :: &
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el !< element
end function plastic_dislotwin_homogenizedC real(pReal), intent(in), dimension(3,3) :: &
Fe
real(pReal), intent(in), dimension(6,6) :: &
C
end subroutine source_damage_isoBrittle_deltaState
module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e) module subroutine plastic_results
integer, intent(in) :: & end subroutine plastic_results
instance, &
i, &
e
type(rotation), dimension(1,discretization_nIP,discretization_nElem), intent(in) :: &
orientation !< crystal orientation
end subroutine plastic_nonlocal_updateCompatibility
module subroutine damage_results
module subroutine plastic_isotropic_results(instance,group) end subroutine damage_results
integer, intent(in) :: instance
character(len=*), intent(in) :: group
end subroutine plastic_isotropic_results
module subroutine plastic_phenopowerlaw_results(instance,group)
integer, intent(in) :: instance
character(len=*), intent(in) :: group
end subroutine plastic_phenopowerlaw_results
module subroutine plastic_kinehardening_results(instance,group)
integer, intent(in) :: instance
character(len=*), intent(in) :: group
end subroutine plastic_kinehardening_results
module subroutine plastic_dislotwin_results(instance,group)
integer, intent(in) :: instance
character(len=*), intent(in) :: group
end subroutine plastic_dislotwin_results
module subroutine plastic_disloUCLA_results(instance,group)
integer, intent(in) :: instance
character(len=*), intent(in) :: group
end subroutine plastic_disloUCLA_results
module subroutine plastic_nonlocal_results(instance,group)
integer, intent(in) :: instance
character(len=*), intent(in) :: group
end subroutine plastic_nonlocal_results
end interface end interface
interface constitutive_LpAndItsTangents
module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
S, Fi, ipc, ip, el)
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(3,3) :: &
S, & !< 2nd Piola-Kirchhoff stress
Fi !< intermediate deformation gradient
real(pReal), intent(out), dimension(3,3) :: &
Lp !< plastic velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLp_dS, &
dLp_dFi !< derivative of Lp with respect to Fi
end subroutine constitutive_plastic_LpAndItsTangents
end interface constitutive_LpAndItsTangents
interface constitutive_dependentState
module subroutine constitutive_plastic_dependentState(F, Fp, ipc, ip, el)
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(3,3) :: &
F, & !< elastic deformation gradient
Fp !< plastic deformation gradient
end subroutine constitutive_plastic_dependentState
end interface constitutive_dependentState
type :: tDebugOptions type :: tDebugOptions
logical :: & logical :: &
@ -333,23 +317,25 @@ module constitutive
type(tDebugOptions) :: debugConstitutive type(tDebugOptions) :: debugConstitutive
public :: & public :: &
plastic_nonlocal_updateCompatibility, &
constitutive_init, & constitutive_init, &
constitutive_homogenizedC, & constitutive_homogenizedC, &
constitutive_dependentState, &
constitutive_LpAndItsTangents, & constitutive_LpAndItsTangents, &
constitutive_dependentState, &
constitutive_LiAndItsTangents, & constitutive_LiAndItsTangents, &
constitutive_initialFi, & constitutive_initialFi, &
constitutive_SandItsTangents, & constitutive_SandItsTangents, &
constitutive_collectDotState, & constitutive_collectDotState, &
constitutive_deltaState, & constitutive_deltaState, &
plastic_nonlocal_updateCompatibility, &
constitutive_damage_getRateAndItsTangents, &
constitutive_thermal_getRateAndItsTangents, &
constitutive_results constitutive_results
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief allocates arrays pointing to array of the various constitutive modules !> @brief Initialze constitutive models for individual physics
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_init subroutine constitutive_init
@ -367,33 +353,10 @@ subroutine constitutive_init
debugConstitutive%ip = debug_root%get_asInt('integrationpoint',defaultVal = 1) debugConstitutive%ip = debug_root%get_asInt('integrationpoint',defaultVal = 1)
debugConstitutive%grain = debug_root%get_asInt('grain',defaultVal = 1) debugConstitutive%grain = debug_root%get_asInt('grain',defaultVal = 1)
!--------------------------------------------------------------------------------------------------
! initialized plasticity
if (any(phase_plasticity == PLASTICITY_NONE_ID)) call plastic_none_init
if (any(phase_plasticity == PLASTICITY_ISOTROPIC_ID)) call plastic_isotropic_init
if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_init
if (any(phase_plasticity == PLASTICITY_KINEHARDENING_ID)) call plastic_kinehardening_init
if (any(phase_plasticity == PLASTICITY_DISLOTWIN_ID)) call plastic_dislotwin_init
if (any(phase_plasticity == PLASTICITY_DISLOUCLA_ID)) call plastic_disloucla_init
if (any(phase_plasticity == PLASTICITY_NONLOCAL_ID)) then
call plastic_nonlocal_init
else
call geometry_plastic_nonlocal_disable
endif
!--------------------------------------------------------------------------------------------------
! initialize source mechanisms
if (any(phase_source == SOURCE_thermal_dissipation_ID)) call source_thermal_dissipation_init
if (any(phase_source == SOURCE_thermal_externalheat_ID)) call source_thermal_externalheat_init
if (any(phase_source == SOURCE_damage_isoBrittle_ID)) call source_damage_isoBrittle_init
if (any(phase_source == SOURCE_damage_isoDuctile_ID)) call source_damage_isoDuctile_init
if (any(phase_source == SOURCE_damage_anisoBrittle_ID)) call source_damage_anisoBrittle_init
if (any(phase_source == SOURCE_damage_anisoDuctile_ID)) call source_damage_anisoDuctile_init
!-------------------------------------------------------------------------------------------------- call plastic_init
! initialize kinematic mechanisms call damage_init
if (any(phase_kinematics == KINEMATICS_cleavage_opening_ID)) call kinematics_cleavage_opening_init call thermal_init
if (any(phase_kinematics == KINEMATICS_slipplane_opening_ID)) call kinematics_slipplane_opening_init
if (any(phase_kinematics == KINEMATICS_thermal_expansion_ID)) call kinematics_thermal_expansion_init
write(6,'(/,a)') ' <<<+- constitutive init -+>>>'; flush(6) write(6,'(/,a)') ' <<<+- constitutive init -+>>>'; flush(6)
@ -407,8 +370,7 @@ subroutine constitutive_init
sourceState(ph)%p(s)%partionedState0 = sourceState(ph)%p(s)%state0 sourceState(ph)%p(s)%partionedState0 = sourceState(ph)%p(s)%state0
sourceState(ph)%p(s)%state = sourceState(ph)%p(s)%partionedState0 sourceState(ph)%p(s)%state = sourceState(ph)%p(s)%partionedState0
end forall end forall
!--------------------------------------------------------------------------------------------------
! determine max size of source state
constitutive_source_maxSizeDotState = max(constitutive_source_maxSizeDotState, & constitutive_source_maxSizeDotState = max(constitutive_source_maxSizeDotState, &
maxval(sourceState(ph)%p%sizeDotState)) maxval(sourceState(ph)%p%sizeDotState))
enddo PhaseLoop2 enddo PhaseLoop2
@ -416,134 +378,29 @@ subroutine constitutive_init
end subroutine constitutive_init end subroutine constitutive_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns the homogenize elasticity matrix !> @brief returns the homogenize elasticity matrix
!> ToDo: homogenizedC66 would be more consistent !> ToDo: homogenizedC66 would be more consistent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function constitutive_homogenizedC(ipc,ip,el) function constitutive_homogenizedC(ipc,ip,el)
real(pReal), dimension(6,6) :: constitutive_homogenizedC real(pReal), dimension(6,6) :: &
integer, intent(in) :: & constitutive_homogenizedC
integer, intent(in) :: &
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el !< element el !< element
plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el))) plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el)))
case (PLASTICITY_DISLOTWIN_ID) plasticityType case (PLASTICITY_DISLOTWIN_ID) plasticityType
constitutive_homogenizedC = plastic_dislotwin_homogenizedC(ipc,ip,el) constitutive_homogenizedC = plastic_dislotwin_homogenizedC(ipc,ip,el)
case default plasticityType case default plasticityType
constitutive_homogenizedC = lattice_C66(1:6,1:6,material_phaseAt(ipc,el)) constitutive_homogenizedC = lattice_C66(1:6,1:6,material_phaseAt(ipc,el))
end select plasticityType end select plasticityType
end function constitutive_homogenizedC end function constitutive_homogenizedC
!--------------------------------------------------------------------------------------------------
!> @brief calls microstructure function of the different constitutive models
!--------------------------------------------------------------------------------------------------
subroutine constitutive_dependentState(F, Fp, ipc, ip, el)
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(3,3) :: &
F, & !< elastic deformation gradient
Fp !< plastic deformation gradient
integer :: &
ho, & !< homogenization
tme, & !< thermal member position
instance, of
ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el)
of = material_phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phaseAt(ipc,el))
plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el)))
case (PLASTICITY_DISLOTWIN_ID) plasticityType
call plastic_dislotwin_dependentState(temperature(ho)%p(tme),instance,of)
case (PLASTICITY_DISLOUCLA_ID) plasticityType
call plastic_disloUCLA_dependentState(instance,of)
case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_dependentState (F,Fp,instance,of,ip,el)
end select plasticityType
end subroutine constitutive_dependentState
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient
! ToDo: Discuss whether it makes sense if crystallite handles the configuration conversion, i.e.
! Mp in, dLp_dMp out
!--------------------------------------------------------------------------------------------------
subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
S, Fi, ipc, ip, el)
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(3,3) :: &
S, & !< 2nd Piola-Kirchhoff stress
Fi !< intermediate deformation gradient
real(pReal), intent(out), dimension(3,3) :: &
Lp !< plastic velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLp_dS, &
dLp_dFi !< derivative of Lp with respect to Fi
real(pReal), dimension(3,3,3,3) :: &
dLp_dMp !< derivative of Lp with respect to Mandel stress
real(pReal), dimension(3,3) :: &
Mp !< Mandel stress work conjugate with Lp
integer :: &
ho, & !< homogenization
tme !< thermal member position
integer :: &
i, j, instance, of
ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el)
Mp = matmul(matmul(transpose(Fi),Fi),S)
of = material_phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phaseAt(ipc,el))
plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el)))
case (PLASTICITY_NONE_ID) plasticityType
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
case (PLASTICITY_ISOTROPIC_ID) plasticityType
call plastic_isotropic_LpAndItsTangent (Lp,dLp_dMp,Mp,instance,of)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
call plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
case (PLASTICITY_KINEHARDENING_ID) plasticityType
call plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMp,Mp, temperature(ho)%p(tme),instance,of,ip,el)
case (PLASTICITY_DISLOTWIN_ID) plasticityType
call plastic_dislotwin_LpAndItsTangent (Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,of)
case (PLASTICITY_DISLOUCLA_ID) plasticityType
call plastic_disloucla_LpAndItsTangent (Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,of)
end select plasticityType
do i=1,3; do j=1,3
dLp_dFi(i,j,1:3,1:3) = matmul(matmul(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + &
matmul(matmul(Fi,dLp_dMp(i,j,1:3,1:3)),S)
dLp_dS(i,j,1:3,1:3) = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi)
enddo; enddo
end subroutine constitutive_LpAndItsTangents
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient !> @brief contains the constitutive equation for calculating the velocity gradient
! ToDo: MD: S is Mi? ! ToDo: MD: S is Mi?
@ -739,7 +596,7 @@ function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el
integer, intent(in) :: & integer, intent(in) :: &
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
el, & !< element el, & !< element
phase, & phase, &
of of
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
@ -794,7 +651,7 @@ function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el
sourceType: select case (phase_source(i,phase)) sourceType: select case (phase_source(i,phase))
case (SOURCE_damage_anisoBrittle_ID) sourceType case (SOURCE_damage_anisoBrittle_ID) sourceType
call source_damage_anisoBrittle_dotState (S, ipc, ip, el) !< correct stress? call source_damage_anisoBrittle_dotState (S, ipc, ip, el) ! correct stress?
case (SOURCE_damage_isoDuctile_ID) sourceType case (SOURCE_damage_isoDuctile_ID) sourceType
call source_damage_isoDuctile_dotState ( ipc, ip, el) call source_damage_isoDuctile_dotState ( ipc, ip, el)
@ -897,37 +754,8 @@ end function constitutive_deltaState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_results subroutine constitutive_results
integer :: p call plastic_results
character(len=pStringLen) :: group call damage_results
do p=1,size(config_name_phase)
group = trim('current/constituent')//'/'//trim(config_name_phase(p))
call results_closeGroup(results_addGroup(group))
group = trim(group)//'/plastic'
call results_closeGroup(results_addGroup(group))
select case(phase_plasticity(p))
case(PLASTICITY_ISOTROPIC_ID)
call plastic_isotropic_results(phase_plasticityInstance(p),group)
case(PLASTICITY_PHENOPOWERLAW_ID)
call plastic_phenopowerlaw_results(phase_plasticityInstance(p),group)
case(PLASTICITY_KINEHARDENING_ID)
call plastic_kinehardening_results(phase_plasticityInstance(p),group)
case(PLASTICITY_DISLOTWIN_ID)
call plastic_dislotwin_results(phase_plasticityInstance(p),group)
case(PLASTICITY_DISLOUCLA_ID)
call plastic_disloUCLA_results(phase_plasticityInstance(p),group)
case(PLASTICITY_NONLOCAL_ID)
call plastic_nonlocal_results(phase_plasticityInstance(p),group)
end select
enddo
end subroutine constitutive_results end subroutine constitutive_results

202
src/constitutive_damage.f90 Normal file
View File

@ -0,0 +1,202 @@
!----------------------------------------------------------------------------------------------------
!> @brief internal microstructure state for all damage sources and kinematics constitutive models
!----------------------------------------------------------------------------------------------------
submodule(constitutive) constitutive_damage
interface
module subroutine source_damage_anisoBrittle_init
end subroutine source_damage_anisoBrittle_init
module subroutine source_damage_anisoDuctile_init
end subroutine source_damage_anisoDuctile_init
module subroutine source_damage_isoBrittle_init
end subroutine source_damage_isoBrittle_init
module subroutine source_damage_isoDuctile_init
end subroutine source_damage_isoDuctile_init
module subroutine kinematics_cleavage_opening_init
end subroutine kinematics_cleavage_opening_init
module subroutine kinematics_slipplane_opening_init
end subroutine kinematics_slipplane_opening_init
module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
integer, intent(in) :: &
phase, & !< phase ID of element
constituent !< position of element within its phase instance
real(pReal), intent(in) :: &
phi !< damage parameter
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
end subroutine source_damage_anisoBrittle_getRateAndItsTangent
module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
integer, intent(in) :: &
phase, & !< phase ID of element
constituent !< position of element within its phase instance
real(pReal), intent(in) :: &
phi !< damage parameter
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
end subroutine source_damage_anisoDuctile_getRateAndItsTangent
module subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
integer, intent(in) :: &
phase, & !< phase ID of element
constituent !< position of element within its phase instance
real(pReal), intent(in) :: &
phi !< damage parameter
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
end subroutine source_damage_isoBrittle_getRateAndItsTangent
module subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
integer, intent(in) :: &
phase, & !< phase ID of element
constituent !< position of element within its phase instance
real(pReal), intent(in) :: &
phi !< damage parameter
real(pReal), intent(out) :: &
localphiDot, &
dLocalphiDot_dPhi
end subroutine source_damage_isoDuctile_getRateAndItsTangent
module subroutine source_damage_anisoBrittle_results(phase,group)
integer, intent(in) :: phase
character(len=*), intent(in) :: group
end subroutine source_damage_anisoBrittle_results
module subroutine source_damage_anisoDuctile_results(phase,group)
integer, intent(in) :: phase
character(len=*), intent(in) :: group
end subroutine source_damage_anisoDuctile_results
module subroutine source_damage_isoBrittle_results(phase,group)
integer, intent(in) :: phase
character(len=*), intent(in) :: group
end subroutine source_damage_isoBrittle_results
module subroutine source_damage_isoDuctile_results(phase,group)
integer, intent(in) :: phase
character(len=*), intent(in) :: group
end subroutine source_damage_isoDuctile_results
end interface
contains
!----------------------------------------------------------------------------------------------
!< @brief initialize damage sources and kinematics mechanism
!----------------------------------------------------------------------------------------------
module subroutine damage_init
! initialize source mechanisms
if (any(phase_source == SOURCE_damage_isoBrittle_ID)) call source_damage_isoBrittle_init
if (any(phase_source == SOURCE_damage_isoDuctile_ID)) call source_damage_isoDuctile_init
if (any(phase_source == SOURCE_damage_anisoBrittle_ID)) call source_damage_anisoBrittle_init
if (any(phase_source == SOURCE_damage_anisoDuctile_ID)) call source_damage_anisoDuctile_init
!--------------------------------------------------------------------------------------------------
! initialize kinematic mechanisms
if (any(phase_kinematics == KINEMATICS_cleavage_opening_ID)) call kinematics_cleavage_opening_init
if (any(phase_kinematics == KINEMATICS_slipplane_opening_ID)) call kinematics_slipplane_opening_init
end subroutine damage_init
!----------------------------------------------------------------------------------------------
!< @brief returns local part of nonlocal damage driving force
!----------------------------------------------------------------------------------------------
module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: &
phi !< damage parameter
real(pReal), intent(inout) :: &
phiDot, &
dPhiDot_dPhi
real(pReal) :: &
localphiDot, &
dLocalphiDot_dPhi
integer :: &
phase, &
grain, &
source, &
constituent
phiDot = 0.0_pReal
dPhiDot_dPhi = 0.0_pReal
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el))
phase = material_phaseAt(grain,el)
constituent = material_phasememberAt(grain,ip,el)
do source = 1, phase_Nsources(phase)
select case(phase_source(source,phase))
case (SOURCE_damage_isoBrittle_ID)
call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case (SOURCE_damage_isoDuctile_ID)
call source_damage_isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case (SOURCE_damage_anisoBrittle_ID)
call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case (SOURCE_damage_anisoDuctile_ID)
call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case default
localphiDot = 0.0_pReal
dLocalphiDot_dPhi = 0.0_pReal
end select
phiDot = phiDot + localphiDot
dPhiDot_dPhi = dPhiDot_dPhi + dLocalphiDot_dPhi
enddo
enddo
end subroutine constitutive_damage_getRateAndItsTangents
!----------------------------------------------------------------------------------------------
!< @brief writes damage sources resultsvto HDF5 output file
!----------------------------------------------------------------------------------------------
module subroutine damage_results
integer :: p,i
character(len=pStringLen) :: group
do p = 1, size(config_name_phase)
sourceLoop: do i = 1, phase_Nsources(p)
group = trim('current/constituent')//'/'//trim(config_name_phase(p))
group = trim(group)//'/sources'
call results_closeGroup(results_addGroup(group))
sourceType: select case (phase_source(i,p))
case (SOURCE_damage_anisoBrittle_ID) sourceType
call source_damage_anisoBrittle_results(p,group)
case (SOURCE_damage_anisoDuctile_ID) sourceType
call source_damage_anisoDuctile_results(p,group)
case (SOURCE_damage_isoBrittle_ID) sourceType
call source_damage_isoBrittle_results(p,group)
case (SOURCE_damage_isoDuctile_ID) sourceType
call source_damage_isoDuctile_results(p,group)
end select sourceType
enddo SourceLoop
enddo
end subroutine damage_results
end submodule constitutive_damage

348
src/constitutive_plastic.f90 Normal file
View File

@ -0,0 +1,348 @@
!----------------------------------------------------------------------------------------------------
!> @brief internal microstructure state for all plasticity constitutive models
!----------------------------------------------------------------------------------------------------
submodule(constitutive) constitutive_plastic
interface
module subroutine plastic_none_init
end subroutine plastic_none_init
module subroutine plastic_isotropic_init
end subroutine plastic_isotropic_init
module subroutine plastic_phenopowerlaw_init
end subroutine plastic_phenopowerlaw_init
module subroutine plastic_kinehardening_init
end subroutine plastic_kinehardening_init
module subroutine plastic_dislotwin_init
end subroutine plastic_dislotwin_init
module subroutine plastic_disloUCLA_init
end subroutine plastic_disloUCLA_init
module subroutine plastic_nonlocal_init
end subroutine plastic_nonlocal_init
module subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
end subroutine plastic_isotropic_LpAndItsTangent
pure module subroutine plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
end subroutine plastic_phenopowerlaw_LpAndItsTangent
pure module subroutine plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
integer, intent(in) :: &
instance, &
of
end subroutine plastic_kinehardening_LpAndItsTangent
module subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of)
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
real(pReal), intent(in) :: &
T
integer, intent(in) :: &
instance, &
of
end subroutine plastic_dislotwin_LpAndItsTangent
pure module subroutine plastic_disloUCLA_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of)
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
real(pReal), intent(in) :: &
T
integer, intent(in) :: &
instance, &
of
end subroutine plastic_disloUCLA_LpAndItsTangent
module subroutine plastic_nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
Mp,Temperature,instance,of,ip,el)
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(3,3,3,3), intent(out) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real(pReal), dimension(3,3), intent(in) :: &
Mp !< Mandel stress
real(pReal), intent(in) :: &
Temperature
integer, intent(in) :: &
instance, &
of, &
ip, & !< current integration point
el !< current element number
end subroutine plastic_nonlocal_LpAndItsTangent
module subroutine plastic_dislotwin_dependentState(T,instance,of)
integer, intent(in) :: &
instance, &
of
real(pReal), intent(in) :: &
T
end subroutine plastic_dislotwin_dependentState
module subroutine plastic_disloUCLA_dependentState(instance,of)
integer, intent(in) :: &
instance, &
of
end subroutine plastic_disloUCLA_dependentState
module subroutine plastic_nonlocal_dependentState(F, Fp, instance, of, ip, el)
real(pReal), dimension(3,3), intent(in) :: &
F, & !< deformation gradient
Fp !< plastic deformation gradient
integer, intent(in) :: &
instance, &
of, &
ip, & !< current integration point
el !< current element number
end subroutine plastic_nonlocal_dependentState
module subroutine plastic_isotropic_results(instance,group)
integer, intent(in) :: instance
character(len=*), intent(in) :: group
end subroutine plastic_isotropic_results
module subroutine plastic_phenopowerlaw_results(instance,group)
integer, intent(in) :: instance
character(len=*), intent(in) :: group
end subroutine plastic_phenopowerlaw_results
module subroutine plastic_kinehardening_results(instance,group)
integer, intent(in) :: instance
character(len=*), intent(in) :: group
end subroutine plastic_kinehardening_results
module subroutine plastic_dislotwin_results(instance,group)
integer, intent(in) :: instance
character(len=*), intent(in) :: group
end subroutine plastic_dislotwin_results
module subroutine plastic_disloUCLA_results(instance,group)
integer, intent(in) :: instance
character(len=*), intent(in) :: group
end subroutine plastic_disloUCLA_results
module subroutine plastic_nonlocal_results(instance,group)
integer, intent(in) :: instance
character(len=*), intent(in) :: group
end subroutine plastic_nonlocal_results
end interface
contains
!--------------------------------------------------------------------------------------------------
!> @brief Initialize constitutive models for plasticity
!--------------------------------------------------------------------------------------------------
module subroutine plastic_init
if (any(phase_plasticity == PLASTICITY_NONE_ID)) call plastic_none_init
if (any(phase_plasticity == PLASTICITY_ISOTROPIC_ID)) call plastic_isotropic_init
if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call plastic_phenopowerlaw_init
if (any(phase_plasticity == PLASTICITY_KINEHARDENING_ID)) call plastic_kinehardening_init
if (any(phase_plasticity == PLASTICITY_DISLOTWIN_ID)) call plastic_dislotwin_init
if (any(phase_plasticity == PLASTICITY_DISLOUCLA_ID)) call plastic_disloucla_init
if (any(phase_plasticity == PLASTICITY_NONLOCAL_ID)) then
call plastic_nonlocal_init
else
call geometry_plastic_nonlocal_disable
endif
end subroutine plastic_init
!--------------------------------------------------------------------------------------------------
!> @brief calls microstructure function of the different plasticity constitutive models
!--------------------------------------------------------------------------------------------------
module subroutine constitutive_plastic_dependentState(F, Fp, ipc, ip, el)
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(3,3) :: &
F, & !< elastic deformation gradient
Fp !< plastic deformation gradient
integer :: &
ho, & !< homogenization
tme, & !< thermal member position
instance, of
ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el)
of = material_phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phaseAt(ipc,el))
plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el)))
case (PLASTICITY_DISLOTWIN_ID) plasticityType
call plastic_dislotwin_dependentState(temperature(ho)%p(tme),instance,of)
case (PLASTICITY_DISLOUCLA_ID) plasticityType
call plastic_disloUCLA_dependentState(instance,of)
case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_dependentState (F,Fp,instance,of,ip,el)
end select plasticityType
end subroutine constitutive_plastic_dependentState
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient
! ToDo: Discuss whether it makes sense if crystallite handles the configuration conversion, i.e.
! Mp in, dLp_dMp out
!--------------------------------------------------------------------------------------------------
module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
S, Fi, ipc, ip, el)
integer, intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(3,3) :: &
S, & !< 2nd Piola-Kirchhoff stress
Fi !< intermediate deformation gradient
real(pReal), intent(out), dimension(3,3) :: &
Lp !< plastic velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLp_dS, &
dLp_dFi !< derivative of Lp with respect to Fi
real(pReal), dimension(3,3,3,3) :: &
dLp_dMp !< derivative of Lp with respect to Mandel stress
real(pReal), dimension(3,3) :: &
Mp !< Mandel stress work conjugate with Lp
integer :: &
ho, & !< homogenization
tme !< thermal member position
integer :: &
i, j, instance, of
ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el)
Mp = matmul(matmul(transpose(Fi),Fi),S)
of = material_phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phaseAt(ipc,el))
plasticityType: select case (phase_plasticity(material_phaseAt(ipc,el)))
case (PLASTICITY_NONE_ID) plasticityType
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
case (PLASTICITY_ISOTROPIC_ID) plasticityType
call plastic_isotropic_LpAndItsTangent (Lp,dLp_dMp,Mp,instance,of)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
call plastic_phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
case (PLASTICITY_KINEHARDENING_ID) plasticityType
call plastic_kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,instance,of)
case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMp,Mp, temperature(ho)%p(tme),instance,of,ip,el)
case (PLASTICITY_DISLOTWIN_ID) plasticityType
call plastic_dislotwin_LpAndItsTangent (Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,of)
case (PLASTICITY_DISLOUCLA_ID) plasticityType
call plastic_disloucla_LpAndItsTangent (Lp,dLp_dMp,Mp,temperature(ho)%p(tme),instance,of)
end select plasticityType
do i=1,3; do j=1,3
dLp_dFi(i,j,1:3,1:3) = matmul(matmul(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + &
matmul(matmul(Fi,dLp_dMp(i,j,1:3,1:3)),S)
dLp_dS(i,j,1:3,1:3) = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi)
enddo; enddo
end subroutine constitutive_plastic_LpAndItsTangents
!--------------------------------------------------------------------------------------------
!> @brief writes plasticity constitutive results to HDF5 output file
!--------------------------------------------------------------------------------------------
module subroutine plastic_results
integer :: p
character(len=pStringLen) :: group
plasticityLoop: do p=1,size(config_name_phase)
group = trim('current/constituent')//'/'//trim(config_name_phase(p))
call results_closeGroup(results_addGroup(group))
group = trim(group)//'/plastic'
call results_closeGroup(results_addGroup(group))
select case(phase_plasticity(p))
case(PLASTICITY_ISOTROPIC_ID)
call plastic_isotropic_results(phase_plasticityInstance(p),group)
case(PLASTICITY_PHENOPOWERLAW_ID)
call plastic_phenopowerlaw_results(phase_plasticityInstance(p),group)
case(PLASTICITY_KINEHARDENING_ID)
call plastic_kinehardening_results(phase_plasticityInstance(p),group)
case(PLASTICITY_DISLOTWIN_ID)
call plastic_dislotwin_results(phase_plasticityInstance(p),group)
case(PLASTICITY_DISLOUCLA_ID)
call plastic_disloUCLA_results(phase_plasticityInstance(p),group)
case(PLASTICITY_NONLOCAL_ID)
call plastic_nonlocal_results(phase_plasticityInstance(p),group)
end select
enddo plasticityLoop
end subroutine plastic_results
end submodule constitutive_plastic

2
src/constitutive_plastic_disloUCLA.f90
View File

@ -5,7 +5,7 @@
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief crystal plasticity model for bcc metals, especially Tungsten !> @brief crystal plasticity model for bcc metals, especially Tungsten
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
submodule(constitutive) plastic_disloUCLA submodule(constitutive:constitutive_plastic) plastic_disloUCLA
real(pReal), parameter :: & real(pReal), parameter :: &
kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin

2
src/constitutive_plastic_dislotwin.f90
View File

@ -7,7 +7,7 @@
!> @brief material subroutine incoprorating dislocation and twinning physics !> @brief material subroutine incoprorating dislocation and twinning physics
!> @details to be done !> @details to be done
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
submodule(constitutive) plastic_dislotwin submodule(constitutive:constitutive_plastic) plastic_dislotwin
real(pReal), parameter :: & real(pReal), parameter :: &
kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin kB = 1.38e-23_pReal !< Boltzmann constant in J/Kelvin

2
src/constitutive_plastic_isotropic.f90
View File

@ -7,7 +7,7 @@
!! resolving the stress on the slip systems. Will give the response of phenopowerlaw for an !! resolving the stress on the slip systems. Will give the response of phenopowerlaw for an
!! untextured polycrystal !! untextured polycrystal
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
submodule(constitutive) plastic_isotropic submodule(constitutive:constitutive_plastic) plastic_isotropic
type :: tParameters type :: tParameters
real(pReal) :: & real(pReal) :: &

2
src/constitutive_plastic_kinehardening.f90
View File

@ -5,7 +5,7 @@
!> @brief Phenomenological crystal plasticity using a power law formulation for the shear rates !> @brief Phenomenological crystal plasticity using a power law formulation for the shear rates
!! and a Voce-type kinematic hardening rule !! and a Voce-type kinematic hardening rule
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
submodule(constitutive) plastic_kinehardening submodule(constitutive:constitutive_plastic) plastic_kinehardening
type :: tParameters type :: tParameters
real(pReal) :: & real(pReal) :: &

2
src/constitutive_plastic_none.f90
View File

@ -4,7 +4,7 @@
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Dummy plasticity for purely elastic material !> @brief Dummy plasticity for purely elastic material
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
submodule(constitutive) plastic_none submodule(constitutive:constitutive_plastic) plastic_none
contains contains

2
src/constitutive_plastic_nonlocal.f90
View File

@ -4,7 +4,7 @@
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine for plasticity including dislocation flux !> @brief material subroutine for plasticity including dislocation flux
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
submodule(constitutive) plastic_nonlocal submodule(constitutive:constitutive_plastic) plastic_nonlocal
use geometry_plastic_nonlocal, only: & use geometry_plastic_nonlocal, only: &
nIPneighbors => geometry_plastic_nonlocal_nIPneighbors, & nIPneighbors => geometry_plastic_nonlocal_nIPneighbors, &
IPneighborhood => geometry_plastic_nonlocal_IPneighborhood, & IPneighborhood => geometry_plastic_nonlocal_IPneighborhood, &

2
src/constitutive_plastic_phenopowerlaw.f90
View File

@ -4,7 +4,7 @@
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief phenomenological crystal plasticity formulation using a powerlaw fitting !> @brief phenomenological crystal plasticity formulation using a powerlaw fitting
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
submodule(constitutive) plastic_phenopowerlaw submodule(constitutive:constitutive_plastic) plastic_phenopowerlaw
type :: tParameters type :: tParameters
real(pReal) :: & real(pReal) :: &

116
src/constitutive_thermal.f90 Normal file
View File

@ -0,0 +1,116 @@
!----------------------------------------------------------------------------------------------------
!> @brief internal microstructure state for all thermal sources and kinematics constitutive models
!----------------------------------------------------------------------------------------------------
submodule(constitutive) constitutive_thermal
interface
module subroutine source_thermal_dissipation_init
end subroutine source_thermal_dissipation_init
module subroutine source_thermal_externalheat_init
end subroutine source_thermal_externalheat_init
module subroutine kinematics_thermal_expansion_init
end subroutine kinematics_thermal_expansion_init
module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase)
integer, intent(in) :: &
phase !< phase ID of element
real(pReal), intent(in), dimension(3,3) :: &
Tstar !< 2nd Piola Kirchoff stress tensor for a given element
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
module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phase, of)
integer, intent(in) :: &
phase, &
of
real(pReal), intent(out) :: &
TDot, &
dTDot_dT
end subroutine source_thermal_externalheat_getRateAndItsTangent
end interface
contains
!----------------------------------------------------------------------------------------------
!< @brief initializes thermal sources and kinematics mechanism
!----------------------------------------------------------------------------------------------
module subroutine thermal_init
! initialize source mechanisms
if (any(phase_source == SOURCE_thermal_dissipation_ID)) call source_thermal_dissipation_init
if (any(phase_source == SOURCE_thermal_externalheat_ID)) call source_thermal_externalheat_init
!--------------------------------------------------------------------------------------------------
!initialize kinematic mechanisms
if (any(phase_kinematics == KINEMATICS_thermal_expansion_ID)) call kinematics_thermal_expansion_init
end subroutine thermal_init
!----------------------------------------------------------------------------------------------
!< @brief calculates thermal dissipation rate
!----------------------------------------------------------------------------------------------
module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, S, Lp, ip, el)
integer, intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), intent(in) :: &
T
real(pReal), intent(in), dimension(:,:,:,:,:) :: &
S, & !< current 2nd Piola Kirchoff stress
Lp !< plastic velocity gradient
real(pReal), intent(inout) :: &
TDot, &
dTDot_dT
real(pReal) :: &
my_Tdot, &
my_dTdot_dT
integer :: &
phase, &
homog, &
instance, &
grain, &
source, &
constituent
homog = material_homogenizationAt(el)
instance = thermal_typeInstance(homog)
do grain = 1, homogenization_Ngrains(homog)
phase = material_phaseAt(grain,el)
constituent = material_phasememberAt(grain,ip,el)
do source = 1, phase_Nsources(phase)
select case(phase_source(source,phase))
case (SOURCE_thermal_dissipation_ID)
call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, &
S(1:3,1:3,grain,ip,el), &
Lp(1:3,1:3,grain,ip,el), &
phase)
case (SOURCE_thermal_externalheat_ID)
call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, &
phase, constituent)
case default
my_Tdot = 0.0_pReal
my_dTdot_dT = 0.0_pReal
end select
Tdot = Tdot + my_Tdot
dTdot_dT = dTdot_dT + my_dTdot_dT
enddo
enddo
end subroutine constitutive_thermal_getRateAndItsTangents
end submodule constitutive_thermal

2
src/crystallite.f90
View File

@ -641,7 +641,7 @@ subroutine crystallite_orientations
if (plasticState(material_phaseAt(1,e))%nonlocal) then if (plasticState(material_phaseAt(1,e))%nonlocal) then
do i = FEsolving_execIP(1),FEsolving_execIP(2) do i = FEsolving_execIP(1),FEsolving_execIP(2)
call plastic_nonlocal_updateCompatibility(crystallite_orientation, & call plastic_nonlocal_updateCompatibility(crystallite_orientation, &
phase_plasticityInstance(material_phaseAt(i,e)),i,e) phase_plasticityInstance(material_phaseAt(1,e)),i,e)
enddo enddo
endif endif
enddo enddo

38
src/damage_local.f90
View File

@ -9,10 +9,7 @@ module damage_local
use config use config
use numerics use numerics
use YAML_types use YAML_types
use source_damage_isoBrittle use constitutive
use source_damage_isoDuctile
use source_damage_anisoBrittle
use source_damage_anisoDuctile
use results use results
implicit none implicit none
@ -128,42 +125,13 @@ subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el
el !< element number el !< element number
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
phi phi
integer :: &
phase, &
grain, &
source, &
constituent
real(pReal) :: & real(pReal) :: &
phiDot, dPhiDot_dPhi, localphiDot, dLocalphiDot_dPhi phiDot, dPhiDot_dPhi
phiDot = 0.0_pReal phiDot = 0.0_pReal
dPhiDot_dPhi = 0.0_pReal dPhiDot_dPhi = 0.0_pReal
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el))
phase = material_phaseAt(grain,el)
constituent = material_phasememberAt(grain,ip,el)
do source = 1, phase_Nsources(phase)
select case(phase_source(source,phase))
case (SOURCE_damage_isoBrittle_ID)
call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case (SOURCE_damage_isoDuctile_ID) call constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el)
call source_damage_isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case (SOURCE_damage_anisoBrittle_ID)
call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case (SOURCE_damage_anisoDuctile_ID)
call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case default
localphiDot = 0.0_pReal
dLocalphiDot_dPhi = 0.0_pReal
end select
phiDot = phiDot + localphiDot
dPhiDot_dPhi = dPhiDot_dPhi + dLocalphiDot_dPhi
enddo
enddo
phiDot = phiDot/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) phiDot = phiDot/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal)
dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal)

39
src/damage_nonlocal.f90
View File

@ -10,10 +10,7 @@ module damage_nonlocal
use YAML_types use YAML_types
use crystallite use crystallite
use lattice use lattice
use source_damage_isoBrittle use constitutive
use source_damage_isoDuctile
use source_damage_anisoBrittle
use source_damage_anisoDuctile
use results use results
implicit none implicit none
@ -97,43 +94,13 @@ subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip,
el !< element number el !< element number
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
phi phi
integer :: &
phase, &
grain, &
source, &
constituent
real(pReal) :: & real(pReal) :: &
phiDot, dPhiDot_dPhi, localphiDot, dLocalphiDot_dPhi phiDot, dPhiDot_dPhi
phiDot = 0.0_pReal phiDot = 0.0_pReal
dPhiDot_dPhi = 0.0_pReal dPhiDot_dPhi = 0.0_pReal
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el))
phase = material_phaseAt(grain,el)
constituent = material_phasememberAt(grain,ip,el)
do source = 1, phase_Nsources(phase)
select case(phase_source(source,phase))
case (SOURCE_damage_isoBrittle_ID)
call source_damage_isobrittle_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case (SOURCE_damage_isoDuctile_ID)
call source_damage_isoductile_getRateAndItsTangent (localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case (SOURCE_damage_anisoBrittle_ID)
call source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case (SOURCE_damage_anisoDuctile_ID)
call source_damage_anisoductile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
case default
localphiDot = 0.0_pReal
dLocalphiDot_dPhi = 0.0_pReal
end select
phiDot = phiDot + localphiDot
dPhiDot_dPhi = dPhiDot_dPhi + dLocalphiDot_dPhi
enddo
enddo
call constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el)
phiDot = phiDot/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) phiDot = phiDot/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal)
dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal)

20
src/kinematics_cleavage_opening.f90
View File

@ -4,16 +4,7 @@
!> @brief material subroutine incorporating kinematics resulting from opening of cleavage planes !> @brief material subroutine incorporating kinematics resulting from opening of cleavage planes
!> @details to be done !> @details to be done
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module kinematics_cleavage_opening submodule(constitutive:constitutive_damage) kinematics_cleavage_opening
use prec
use IO
use config
use math
use lattice
use material
implicit none
private
integer, dimension(:), allocatable :: kinematics_cleavage_opening_instance integer, dimension(:), allocatable :: kinematics_cleavage_opening_instance
@ -31,9 +22,6 @@ module kinematics_cleavage_opening
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance)
public :: &
kinematics_cleavage_opening_init, &
kinematics_cleavage_opening_LiAndItsTangent
contains contains
@ -42,7 +30,7 @@ contains
!> @brief module initialization !> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks !> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine kinematics_cleavage_opening_init module subroutine kinematics_cleavage_opening_init
integer :: Ninstance,p integer :: Ninstance,p
integer, dimension(:), allocatable :: N_cl !< active number of cleavage systems per family integer, dimension(:), allocatable :: N_cl !< active number of cleavage systems per family
@ -95,7 +83,7 @@ end subroutine kinematics_cleavage_opening_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient !> @brief contains the constitutive equation for calculating the velocity gradient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
ipc, & !< grain number ipc, & !< grain number
@ -158,4 +146,4 @@ subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, i
end subroutine kinematics_cleavage_opening_LiAndItsTangent end subroutine kinematics_cleavage_opening_LiAndItsTangent
end module kinematics_cleavage_opening end submodule kinematics_cleavage_opening

20
src/kinematics_slipplane_opening.f90
View File

@ -4,16 +4,7 @@
!> @brief material subroutine incorporating kinematics resulting from opening of slip planes !> @brief material subroutine incorporating kinematics resulting from opening of slip planes
!> @details to be done !> @details to be done
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module kinematics_slipplane_opening submodule(constitutive:constitutive_damage) kinematics_slipplane_opening
use prec
use config
use IO
use math
use lattice
use material
implicit none
private
integer, dimension(:), allocatable :: kinematics_slipplane_opening_instance integer, dimension(:), allocatable :: kinematics_slipplane_opening_instance
@ -33,9 +24,6 @@ module kinematics_slipplane_opening
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance)
public :: &
kinematics_slipplane_opening_init, &
kinematics_slipplane_opening_LiAndItsTangent
contains contains
@ -44,7 +32,7 @@ contains
!> @brief module initialization !> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks !> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine kinematics_slipplane_opening_init module subroutine kinematics_slipplane_opening_init
integer :: Ninstance,p,i integer :: Ninstance,p,i
character(len=pStringLen) :: extmsg = '' character(len=pStringLen) :: extmsg = ''
@ -107,7 +95,7 @@ end subroutine kinematics_slipplane_opening_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient !> @brief contains the constitutive equation for calculating the velocity gradient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el) module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
ipc, & !< grain number ipc, & !< grain number
@ -183,4 +171,4 @@ subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ipc,
end subroutine kinematics_slipplane_opening_LiAndItsTangent end subroutine kinematics_slipplane_opening_LiAndItsTangent
end module kinematics_slipplane_opening end submodule kinematics_slipplane_opening

47
src/kinematics_thermal_expansion.f90
View File

@ -3,16 +3,7 @@
!> @brief material subroutine incorporating kinematics resulting from thermal expansion !> @brief material subroutine incorporating kinematics resulting from thermal expansion
!> @details to be done !> @details to be done
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module kinematics_thermal_expansion submodule(constitutive:constitutive_thermal) kinematics_thermal_expansion
use prec
use IO
use config
use math
use lattice
use material
implicit none
private
integer, dimension(:), allocatable :: kinematics_thermal_expansion_instance integer, dimension(:), allocatable :: kinematics_thermal_expansion_instance
@ -25,10 +16,6 @@ module kinematics_thermal_expansion
type(tParameters), dimension(:), allocatable :: param type(tParameters), dimension(:), allocatable :: param
public :: &
kinematics_thermal_expansion_init, &
kinematics_thermal_expansion_initialStrain, &
kinematics_thermal_expansion_LiAndItsTangent
contains contains
@ -37,7 +24,7 @@ contains
!> @brief module initialization !> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks !> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine kinematics_thermal_expansion_init module subroutine kinematics_thermal_expansion_init
integer :: Ninstance,p,i integer :: Ninstance,p,i
real(pReal), dimension(:), allocatable :: temp real(pReal), dimension(:), allocatable :: temp
@ -79,22 +66,22 @@ end subroutine kinematics_thermal_expansion_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief report initial thermal strain based on current temperature deviation from reference !> @brief report initial thermal strain based on current temperature deviation from reference
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function kinematics_thermal_expansion_initialStrain(homog,phase,offset) pure module function kinematics_thermal_expansion_initialStrain(homog,phase,offset) result(initialStrain)
integer, intent(in) :: & integer, intent(in) :: &
phase, & phase, &
homog, & homog, &
offset offset
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
kinematics_thermal_expansion_initialStrain !< initial thermal strain (should be small strain, though) initialStrain !< initial thermal strain (should be small strain, though)
associate(prm => param(kinematics_thermal_expansion_instance(phase))) associate(prm => param(kinematics_thermal_expansion_instance(phase)))
kinematics_thermal_expansion_initialStrain = & initialStrain = &
(temperature(homog)%p(offset) - prm%T_ref)**1 / 1. * prm%expansion(1:3,1:3,1) + & ! constant coefficient (temperature(homog)%p(offset) - prm%T_ref)**1 / 1. * prm%expansion(1:3,1:3,1) + & ! constant coefficient
(temperature(homog)%p(offset) - prm%T_ref)**2 / 2. * prm%expansion(1:3,1:3,2) + & ! linear coefficient (temperature(homog)%p(offset) - prm%T_ref)**2 / 2. * prm%expansion(1:3,1:3,2) + & ! linear coefficient
(temperature(homog)%p(offset) - prm%T_ref)**3 / 3. * prm%expansion(1:3,1:3,3) ! quadratic coefficient (temperature(homog)%p(offset) - prm%T_ref)**3 / 3. * prm%expansion(1:3,1:3,3) ! quadratic coefficient
end associate end associate
end function kinematics_thermal_expansion_initialStrain end function kinematics_thermal_expansion_initialStrain
@ -102,7 +89,7 @@ end function kinematics_thermal_expansion_initialStrain
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief constitutive equation for calculating the velocity gradient !> @brief constitutive equation for calculating the velocity gradient
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, el) module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
ipc, & !< grain number ipc, & !< grain number
@ -140,4 +127,4 @@ subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, ipc, ip,
end subroutine kinematics_thermal_expansion_LiAndItsTangent end subroutine kinematics_thermal_expansion_LiAndItsTangent
end module kinematics_thermal_expansion end submodule kinematics_thermal_expansion

88
src/marc/discretization_marc.f90
View File

@ -119,7 +119,7 @@ subroutine discretization_marc_init
unscaledNormals = IPareaNormal(elem,nElems,connectivity_cell,node0_cell) unscaledNormals = IPareaNormal(elem,nElems,connectivity_cell,node0_cell)
call geometry_plastic_nonlocal_setIParea(norm2(unscaledNormals,1)) call geometry_plastic_nonlocal_setIParea(norm2(unscaledNormals,1))
call geometry_plastic_nonlocal_setIPareaNormal(unscaledNormals/spread(norm2(unscaledNormals,1),1,3)) call geometry_plastic_nonlocal_setIPareaNormal(unscaledNormals/spread(norm2(unscaledNormals,1),1,3))
!call geometry_plastic_nonlocal_setIPneighborhood ToDo: Support nonlocal call geometry_plastic_nonlocal_setIPneighborhood(IPneighborhood(elem,connectivity_cell))
call geometry_plastic_nonlocal_results call geometry_plastic_nonlocal_results
end subroutine discretization_marc_init end subroutine discretization_marc_init
@ -733,10 +733,10 @@ end subroutine inputRead_microstructureAndHomogenization
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Determine cell connectivity and definition of cell nodes !> @brief Calculates cell node coordinates from element node coordinates
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine buildCells(connectivity_cell,cellNodeDefinition, & pure subroutine buildCells(connectivity_cell,cellNodeDefinition, &
elem,connectivity_elem) elem,connectivity_elem)
type(tCellNodeDefinition), dimension(:), intent(out) :: cellNodeDefinition ! definition of cell nodes for increasing number of parents type(tCellNodeDefinition), dimension(:), intent(out) :: cellNodeDefinition ! definition of cell nodes for increasing number of parents
integer, dimension(:,:,:),intent(out) :: connectivity_cell integer, dimension(:,:,:),intent(out) :: connectivity_cell
@ -747,8 +747,7 @@ subroutine buildCells(connectivity_cell,cellNodeDefinition, &
integer,dimension(:), allocatable :: candidates_local integer,dimension(:), allocatable :: candidates_local
integer,dimension(:,:), allocatable :: parentsAndWeights,candidates_global integer,dimension(:,:), allocatable :: parentsAndWeights,candidates_global
integer :: e,n,c,p,s,i,m,j,& integer :: e, n, c, p, s,i,m,j,nParentNodes,nCellNode,Nelem,candidateID
nParentNodes,nCellNode,Nelem,candidateID
Nelem = size(connectivity_elem,2) Nelem = size(connectivity_elem,2)
@ -762,7 +761,9 @@ subroutine buildCells(connectivity_cell,cellNodeDefinition, &
do e = 1, Nelem do e = 1, Nelem
do c = 1, elem%NcellNodes do c = 1, elem%NcellNodes
realNode: if (count(elem%cellNodeParentNodeWeights(:,c) /= 0) == 1) then realNode: if (count(elem%cellNodeParentNodeWeights(:,c) /= 0) == 1) then
where(connectivity_cell(:,:,e) == -c) connectivity_cell(:,:,e) = connectivity_elem(c,e) where(connectivity_cell(:,:,e) == -c)
connectivity_cell(:,:,e) = connectivity_elem(c,e)
end where
endif realNode endif realNode
enddo enddo
enddo enddo
@ -889,10 +890,10 @@ end subroutine buildCells
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Calculate cell node coordinates from element node coordinates !> @brief Calculates cell node coordinates from element node coordinates
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine buildCellNodes(node_cell, & pure subroutine buildCellNodes(node_cell, &
definition,node_elem) definition,node_elem)
real(pReal), dimension(:,:), intent(out) :: node_cell !< cell node coordinates real(pReal), dimension(:,:), intent(out) :: node_cell !< cell node coordinates
type(tCellNodeDefinition), dimension(:), intent(in) :: definition !< cell node definition (weights and parents) type(tCellNodeDefinition), dimension(:), intent(in) :: definition !< cell node definition (weights and parents)
@ -919,10 +920,10 @@ end subroutine buildCellNodes
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Calculate IP coordinates as center of cell !> @brief Calculates IP coordinates as center of cell
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine buildIPcoordinates(IPcoordinates, & pure subroutine buildIPcoordinates(IPcoordinates, &
connectivity_cell,node_cell) connectivity_cell,node_cell)
real(pReal), dimension(:,:), intent(out):: IPcoordinates !< cell-center/IP coordinates real(pReal), dimension(:,:), intent(out):: IPcoordinates !< cell-center/IP coordinates
integer, dimension(:,:), intent(in) :: connectivity_cell !< connectivity for each cell integer, dimension(:,:), intent(in) :: connectivity_cell !< connectivity for each cell
@ -947,7 +948,7 @@ end subroutine buildIPcoordinates
!> @details The IP volume is calculated differently depending on the cell type. !> @details The IP volume is calculated differently depending on the cell type.
!> 2D cells assume an element depth of 1.0 !> 2D cells assume an element depth of 1.0
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
function IPvolume(elem,node,connectivity) pure function IPvolume(elem,node,connectivity)
type(tElement), intent(in) :: elem type(tElement), intent(in) :: elem
real(pReal), dimension(:,:), intent(in) :: node real(pReal), dimension(:,:), intent(in) :: node
@ -1005,7 +1006,7 @@ end function IPvolume
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculation of IP interface areas !> @brief calculation of IP interface areas
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function IPareaNormal(elem,nElem,connectivity,node) pure function IPareaNormal(elem,nElem,connectivity,node)
type(tElement), intent(in) :: elem type(tElement), intent(in) :: elem
integer, intent(in) :: nElem integer, intent(in) :: nElem
@ -1051,6 +1052,63 @@ function IPareaNormal(elem,nElem,connectivity,node)
end function IPareaNormal end function IPareaNormal
!--------------------------------------------------------------------------------------------------
!> @brief IP neighborhood
!--------------------------------------------------------------------------------------------------
function IPneighborhood(elem,connectivity)
type(tElement), intent(in) :: elem ! definition of the element in use
integer, dimension(:,:,:), intent(in) :: connectivity ! cell connectivity
integer, dimension(3,size(elem%cellFace,2), &
size(connectivity,2),size(connectivity,3)) :: IPneighborhood ! neighboring IPs as [element ID, IP ID, face ID]
integer, dimension(size(elem%cellFace,1)+3,&
size(elem%cellFace,2)*size(connectivity,2)*size(connectivity,3)) :: face
integer, dimension(size(connectivity,1)) :: myConnectivity
integer, dimension(size(elem%cellFace,1)) :: face_unordered
integer :: e,i,f,n,c,s
c = 0
do e = 1, size(connectivity,3)
do i = 1, size(connectivity,2)
myConnectivity = connectivity(:,i,e)
do f = 1, size(elem%cellFace,2)
c = c + 1
face_unordered = myConnectivity(elem%cellFace(:,f))
do n = 1, size(face_unordered)
face(n,c) = minval(face_unordered)
face_unordered(minloc(face_unordered)) = huge(face_unordered)
enddo
face(n:n+3,c) = [e,i,f]
enddo
enddo; enddo
!--------------------------------------------------------------------------------------------------
! sort face definitions
call math_sort(face,sortDim=1)
do c=2, size(face,1)-4
s = 1
e = 1
do while (e < size(face,2))
e = e + 1
if(any(face(:c,s) /= face(:c,e))) then
if(e-1/=s) call math_sort(face(:,s:e-1),sortDim=c)
s = e
endif
enddo
enddo
IPneighborhood = 0
do c=1, size(face,2) - 1
if(all(face(:n-1,c) == face(:n-1,c+1))) then
IPneighborhood(:,face(n+2,c+1),face(n+1,c+1),face(n+0,c+1)) = face(n:n+3,c+0)
IPneighborhood(:,face(n+2,c+0),face(n+1,c+0),face(n+0,c+0)) = face(n:n+3,c+1)
endif
enddo
end function IPneighborhood
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return integer list corresponding to items in consecutive lines. !> @brief return integer list corresponding to items in consecutive lines.
!! First integer in array is counter !! First integer in array is counter

18
src/math.f90
View File

@ -126,11 +126,12 @@ end subroutine math_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Quicksort algorithm for two-dimensional integer arrays !> @brief Sorting of two-dimensional integer arrays
! Sorting is done with respect to array(sort,:) and keeps array(/=sort,:) linked to it. !> @details Based on quicksort.
! default: sort=1 ! Sorting is done with respect to array(sortDim,:) and keeps array(/=sortDim,:) linked to it.
! Default: sortDim=1
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
recursive subroutine math_sort(a, istart, iend, sortDim) pure recursive subroutine math_sort(a, istart, iend, sortDim)
integer, dimension(:,:), intent(inout) :: a integer, dimension(:,:), intent(inout) :: a
integer, intent(in),optional :: istart,iend, sortDim integer, intent(in),optional :: istart,iend, sortDim
@ -155,7 +156,7 @@ recursive subroutine math_sort(a, istart, iend, sortDim)
endif endif
if (s < e) then if (s < e) then
ipivot = qsort_partition(a,s, e, d) call qsort_partition(a,ipivot, s,e, d)
call math_sort(a, s, ipivot-1, d) call math_sort(a, s, ipivot-1, d)
call math_sort(a, ipivot+1, e, d) call math_sort(a, ipivot+1, e, d)
endif endif
@ -166,9 +167,10 @@ recursive subroutine math_sort(a, istart, iend, sortDim)
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
!> @brief Partitioning required for quicksort !> @brief Partitioning required for quicksort
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
integer function qsort_partition(a, istart, iend, sort) pure subroutine qsort_partition(a,p, istart, iend, sort)
integer, dimension(:,:), intent(inout) :: a integer, dimension(:,:), intent(inout) :: a
integer, intent(out) :: p ! Pivot element
integer, intent(in) :: istart,iend,sort integer, intent(in) :: istart,iend,sort
integer, dimension(size(a,1)) :: tmp integer, dimension(size(a,1)) :: tmp
integer :: i,j integer :: i,j
@ -186,7 +188,7 @@ recursive subroutine math_sort(a, istart, iend, sortDim)
tmp = a(:,istart) tmp = a(:,istart)
a(:,istart) = a(:,j) a(:,istart) = a(:,j)
a(:,j) = tmp a(:,j) = tmp
qsort_partition = j p = j
return return
else cross ! exchange values else cross ! exchange values
tmp = a(:,i) tmp = a(:,i)
@ -195,7 +197,7 @@ recursive subroutine math_sort(a, istart, iend, sortDim)
endif cross endif cross
enddo enddo
end function qsort_partition end subroutine qsort_partition
end subroutine math_sort end subroutine math_sort

29
src/source_damage_anisoBrittle.f90
View File

@ -4,18 +4,7 @@
!> @brief material subroutine incorporating anisotropic brittle damage source mechanism !> @brief material subroutine incorporating anisotropic brittle damage source mechanism
!> @details to be done !> @details to be done
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module source_damage_anisoBrittle submodule (constitutive:constitutive_damage) source_damage_anisoBrittle
use prec
use IO
use math
use discretization
use material
use config
use lattice
use results
implicit none
private
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
source_damage_anisoBrittle_offset, & !< which source is my current source mechanism? source_damage_anisoBrittle_offset, & !< which source is my current source mechanism?
@ -39,12 +28,6 @@ module source_damage_anisoBrittle
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance)
public :: &
source_damage_anisoBrittle_init, &
source_damage_anisoBrittle_dotState, &
source_damage_anisobrittle_getRateAndItsTangent, &
source_damage_anisoBrittle_results
contains contains
@ -52,7 +35,7 @@ contains
!> @brief module initialization !> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks !> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine source_damage_anisoBrittle_init module subroutine source_damage_anisoBrittle_init
integer :: Ninstance,sourceOffset,NipcMyPhase,p integer :: Ninstance,sourceOffset,NipcMyPhase,p
integer, dimension(:), allocatable :: N_cl integer, dimension(:), allocatable :: N_cl
@ -123,7 +106,7 @@ end subroutine source_damage_anisoBrittle_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state !> @brief calculates derived quantities from state
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el) module subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
@ -172,7 +155,7 @@ end subroutine source_damage_anisoBrittle_dotState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns local part of nonlocal damage driving force !> @brief returns local part of nonlocal damage driving force
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) module subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
integer, intent(in) :: & integer, intent(in) :: &
phase, & phase, &
@ -199,7 +182,7 @@ end subroutine source_damage_anisoBrittle_getRateAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file !> @brief writes results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine source_damage_anisoBrittle_results(phase,group) module subroutine source_damage_anisoBrittle_results(phase,group)
integer, intent(in) :: phase integer, intent(in) :: phase
character(len=*), intent(in) :: group character(len=*), intent(in) :: group
@ -218,4 +201,4 @@ subroutine source_damage_anisoBrittle_results(phase,group)
end subroutine source_damage_anisoBrittle_results end subroutine source_damage_anisoBrittle_results
end module source_damage_anisoBrittle end submodule source_damage_anisoBrittle

33
src/source_damage_anisoDuctile.f90
View File

@ -4,23 +4,13 @@
!> @brief material subroutine incorporating anisotropic ductile damage source mechanism !> @brief material subroutine incorporating anisotropic ductile damage source mechanism
!> @details to be done !> @details to be done
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module source_damage_anisoDuctile submodule(constitutive:constitutive_damage) source_damage_anisoDuctile
use prec
use IO
use math
use discretization
use material
use config
use results
implicit none
private
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
source_damage_anisoDuctile_offset, & !< which source is my current damage mechanism? source_damage_anisoDuctile_offset, & !< which source is my current damage mechanism?
source_damage_anisoDuctile_instance !< instance of damage source mechanism source_damage_anisoDuctile_instance !< instance of damage source mechanism
type, private :: tParameters !< container type for internal constitutive parameters type :: tParameters !< container type for internal constitutive parameters
real(pReal) :: & real(pReal) :: &
n n
real(pReal), dimension(:), allocatable :: & real(pReal), dimension(:), allocatable :: &
@ -29,14 +19,7 @@ module source_damage_anisoDuctile
output output
end type tParameters end type tParameters
type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance)
public :: &
source_damage_anisoDuctile_init, &
source_damage_anisoDuctile_dotState, &
source_damage_anisoDuctile_getRateAndItsTangent, &
source_damage_anisoDuctile_results
contains contains
@ -45,7 +28,7 @@ contains
!> @brief module initialization !> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks !> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine source_damage_anisoDuctile_init module subroutine source_damage_anisoDuctile_init
integer :: Ninstance,sourceOffset,NipcMyPhase,p integer :: Ninstance,sourceOffset,NipcMyPhase,p
integer, dimension(:), allocatable :: N_sl integer, dimension(:), allocatable :: N_sl
@ -105,7 +88,7 @@ end subroutine source_damage_anisoDuctile_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state !> @brief calculates derived quantities from state
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine source_damage_anisoDuctile_dotState(ipc, ip, el) module subroutine source_damage_anisoDuctile_dotState(ipc, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
@ -136,7 +119,7 @@ end subroutine source_damage_anisoDuctile_dotState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns local part of nonlocal damage driving force !> @brief returns local part of nonlocal damage driving force
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) module subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
integer, intent(in) :: & integer, intent(in) :: &
phase, & phase, &
@ -163,7 +146,7 @@ end subroutine source_damage_anisoDuctile_getRateAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file !> @brief writes results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine source_damage_anisoDuctile_results(phase,group) module subroutine source_damage_anisoDuctile_results(phase,group)
integer, intent(in) :: phase integer, intent(in) :: phase
character(len=*), intent(in) :: group character(len=*), intent(in) :: group
@ -182,4 +165,4 @@ subroutine source_damage_anisoDuctile_results(phase,group)
end subroutine source_damage_anisoDuctile_results end subroutine source_damage_anisoDuctile_results
end module source_damage_anisoDuctile end submodule source_damage_anisoDuctile

31
src/source_damage_isoBrittle.f90
View File

@ -4,23 +4,13 @@
!> @brief material subroutine incoprorating isotropic brittle damage source mechanism !> @brief material subroutine incoprorating isotropic brittle damage source mechanism
!> @details to be done !> @details to be done
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module source_damage_isoBrittle submodule(constitutive:constitutive_damage) source_damage_isoBrittle
use prec
use IO
use math
use discretization
use material
use config
use results
implicit none
private
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
source_damage_isoBrittle_offset, & source_damage_isoBrittle_offset, &
source_damage_isoBrittle_instance source_damage_isoBrittle_instance
type, private :: tParameters !< container type for internal constitutive parameters type :: tParameters !< container type for internal constitutive parameters
real(pReal) :: & real(pReal) :: &
critStrainEnergy, & critStrainEnergy, &
N N
@ -30,13 +20,6 @@ module source_damage_isoBrittle
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance)
public :: &
source_damage_isoBrittle_init, &
source_damage_isoBrittle_deltaState, &
source_damage_isoBrittle_getRateAndItsTangent, &
source_damage_isoBrittle_results
contains contains
@ -44,7 +27,7 @@ contains
!> @brief module initialization !> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks !> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine source_damage_isoBrittle_init module subroutine source_damage_isoBrittle_init
integer :: Ninstance,sourceOffset,NipcMyPhase,p integer :: Ninstance,sourceOffset,NipcMyPhase,p
character(len=pStringLen) :: extmsg = '' character(len=pStringLen) :: extmsg = ''
@ -99,7 +82,7 @@ end subroutine source_damage_isoBrittle_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state !> @brief calculates derived quantities from state
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el) module subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
@ -145,7 +128,7 @@ end subroutine source_damage_isoBrittle_deltaState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns local part of nonlocal damage driving force !> @brief returns local part of nonlocal damage driving force
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) module subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
integer, intent(in) :: & integer, intent(in) :: &
phase, & phase, &
@ -174,7 +157,7 @@ end subroutine source_damage_isoBrittle_getRateAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file !> @brief writes results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine source_damage_isoBrittle_results(phase,group) module subroutine source_damage_isoBrittle_results(phase,group)
integer, intent(in) :: phase integer, intent(in) :: phase
character(len=*), intent(in) :: group character(len=*), intent(in) :: group
@ -193,4 +176,4 @@ subroutine source_damage_isoBrittle_results(phase,group)
end subroutine source_damage_isoBrittle_results end subroutine source_damage_isoBrittle_results
end module source_damage_isoBrittle end submodule source_damage_isoBrittle

31
src/source_damage_isoDuctile.f90
View File

@ -4,22 +4,13 @@
!> @brief material subroutine incorporating isotropic ductile damage source mechanism !> @brief material subroutine incorporating isotropic ductile damage source mechanism
!> @details to be done !> @details to be done
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module source_damage_isoDuctile submodule (constitutive:constitutive_damage) source_damage_isoDuctile
use prec
use IO
use discretization
use material
use config
use results
implicit none
private
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
source_damage_isoDuctile_offset, & !< which source is my current damage mechanism? source_damage_isoDuctile_offset, & !< which source is my current damage mechanism?
source_damage_isoDuctile_instance !< instance of damage source mechanism source_damage_isoDuctile_instance !< instance of damage source mechanism
type, private :: tParameters !< container type for internal constitutive parameters type:: tParameters !< container type for internal constitutive parameters
real(pReal) :: & real(pReal) :: &
critPlasticStrain, & critPlasticStrain, &
N N
@ -27,15 +18,9 @@ module source_damage_isoDuctile
output output
end type tParameters end type tParameters
type(tParameters), dimension(:), allocatable, private :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance)
public :: &
source_damage_isoDuctile_init, &
source_damage_isoDuctile_dotState, &
source_damage_isoDuctile_getRateAndItsTangent, &
source_damage_isoDuctile_Results
contains contains
@ -43,7 +28,7 @@ contains
!> @brief module initialization !> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks !> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine source_damage_isoDuctile_init module subroutine source_damage_isoDuctile_init
integer :: Ninstance,sourceOffset,NipcMyPhase,p integer :: Ninstance,sourceOffset,NipcMyPhase,p
character(len=pStringLen) :: extmsg = '' character(len=pStringLen) :: extmsg = ''
@ -98,7 +83,7 @@ end subroutine source_damage_isoDuctile_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state !> @brief calculates derived quantities from state
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine source_damage_isoDuctile_dotState(ipc, ip, el) module subroutine source_damage_isoDuctile_dotState(ipc, ip, el)
integer, intent(in) :: & integer, intent(in) :: &
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
@ -129,7 +114,7 @@ end subroutine source_damage_isoDuctile_dotState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns local part of nonlocal damage driving force !> @brief returns local part of nonlocal damage driving force
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent) module subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, phase, constituent)
integer, intent(in) :: & integer, intent(in) :: &
phase, & phase, &
@ -156,7 +141,7 @@ end subroutine source_damage_isoDuctile_getRateAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file !> @brief writes results to HDF5 output file
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine source_damage_isoDuctile_results(phase,group) module subroutine source_damage_isoDuctile_results(phase,group)
integer, intent(in) :: phase integer, intent(in) :: phase
character(len=*), intent(in) :: group character(len=*), intent(in) :: group
@ -175,4 +160,4 @@ subroutine source_damage_isoDuctile_results(phase,group)
end subroutine source_damage_isoDuctile_results end subroutine source_damage_isoDuctile_results
end module source_damage_isoDuctile end submodule source_damage_isoDuctile

19
src/source_thermal_dissipation.f90
View File

@ -4,14 +4,7 @@
!> @brief material subroutine for thermal source due to plastic dissipation !> @brief material subroutine for thermal source due to plastic dissipation
!> @details to be done !> @details to be done
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module source_thermal_dissipation submodule(constitutive:constitutive_thermal) source_thermal_dissipation
use prec
use discretization
use material
use config
implicit none
private
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
source_thermal_dissipation_offset, & !< which source is my current thermal dissipation mechanism? source_thermal_dissipation_offset, & !< which source is my current thermal dissipation mechanism?
@ -25,10 +18,6 @@ module source_thermal_dissipation
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance)
public :: &
source_thermal_dissipation_init, &
source_thermal_dissipation_getRateAndItsTangent
contains contains
@ -36,7 +25,7 @@ contains
!> @brief module initialization !> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks !> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine source_thermal_dissipation_init module subroutine source_thermal_dissipation_init
integer :: Ninstance,sourceOffset,NipcMyPhase,p integer :: Ninstance,sourceOffset,NipcMyPhase,p
@ -76,7 +65,7 @@ end subroutine source_thermal_dissipation_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Ninstances dissipation rate !> @brief Ninstances dissipation rate
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase) module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase)
integer, intent(in) :: & integer, intent(in) :: &
phase phase
@ -96,4 +85,4 @@ subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar
end subroutine source_thermal_dissipation_getRateAndItsTangent end subroutine source_thermal_dissipation_getRateAndItsTangent
end module source_thermal_dissipation end submodule source_thermal_dissipation

21
src/source_thermal_externalheat.f90
View File

@ -4,14 +4,8 @@
!> @author Philip Eisenlohr, Michigan State University !> @author Philip Eisenlohr, Michigan State University
!> @brief material subroutine for variable heat source !> @brief material subroutine for variable heat source
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module source_thermal_externalheat submodule(constitutive:constitutive_thermal) source_thermal_externalheat
use prec
use discretization
use material
use config
implicit none
private
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
source_thermal_externalheat_offset, & !< which source is my current thermal dissipation mechanism? source_thermal_externalheat_offset, & !< which source is my current thermal dissipation mechanism?
@ -28,11 +22,6 @@ module source_thermal_externalheat
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance)
public :: &
source_thermal_externalheat_init, &
source_thermal_externalheat_dotState, &
source_thermal_externalheat_getRateAndItsTangent
contains contains
@ -40,7 +29,7 @@ contains
!> @brief module initialization !> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks !> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine source_thermal_externalheat_init module subroutine source_thermal_externalheat_init
integer :: Ninstance,sourceOffset,NipcMyPhase,p integer :: Ninstance,sourceOffset,NipcMyPhase,p
@ -84,7 +73,7 @@ end subroutine source_thermal_externalheat_init
!> @brief rate of change of state !> @brief rate of change of state
!> @details state only contains current time to linearly interpolate given heat powers !> @details state only contains current time to linearly interpolate given heat powers
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine source_thermal_externalheat_dotState(phase, of) module subroutine source_thermal_externalheat_dotState(phase, of)
integer, intent(in) :: & integer, intent(in) :: &
phase, & phase, &
@ -103,7 +92,7 @@ end subroutine source_thermal_externalheat_dotState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns local heat generation rate !> @brief returns local heat generation rate
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phase, of) module subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phase, of)
integer, intent(in) :: & integer, intent(in) :: &
phase, & phase, &
@ -135,4 +124,4 @@ subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, phas
end subroutine source_thermal_externalheat_getRateAndItsTangent end subroutine source_thermal_externalheat_getRateAndItsTangent
end module source_thermal_externalheat end submodule source_thermal_externalheat

40
src/thermal_adiabatic.f90
View File

@ -8,8 +8,7 @@ module thermal_adiabatic
use numerics use numerics
use material use material
use results use results
use source_thermal_dissipation use constitutive
use source_thermal_externalheat
use crystallite use crystallite
use lattice use lattice
@ -125,45 +124,14 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
T T
real(pReal), intent(out) :: & real(pReal), intent(out) :: &
Tdot, dTdot_dT Tdot, dTdot_dT
real(pReal) :: &
my_Tdot, my_dTdot_dT
integer :: & integer :: &
phase, & homog
homog, &
instance, &
grain, &
source, &
constituent
homog = material_homogenizationAt(el)
instance = thermal_typeInstance(homog)
Tdot = 0.0_pReal Tdot = 0.0_pReal
dTdot_dT = 0.0_pReal dTdot_dT = 0.0_pReal
do grain = 1, homogenization_Ngrains(homog)
phase = material_phaseAt(grain,el)
constituent = material_phasememberAt(grain,ip,el)
do source = 1, phase_Nsources(phase)
select case(phase_source(source,phase))
case (SOURCE_thermal_dissipation_ID)
call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, &
crystallite_S(1:3,1:3,grain,ip,el), &
crystallite_Lp(1:3,1:3,grain,ip,el), &
phase)
case (SOURCE_thermal_externalheat_ID) homog = material_homogenizationAt(el)
call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, crystallite_S, crystallite_Lp, ip, el)
phase, constituent)
case default
my_Tdot = 0.0_pReal
my_dTdot_dT = 0.0_pReal
end select
Tdot = Tdot + my_Tdot
dTdot_dT = dTdot_dT + my_dTdot_dT
enddo
enddo
Tdot = Tdot/real(homogenization_Ngrains(homog),pReal) Tdot = Tdot/real(homogenization_Ngrains(homog),pReal)
dTdot_dT = dTdot_dT/real(homogenization_Ngrains(homog),pReal) dTdot_dT = dTdot_dT/real(homogenization_Ngrains(homog),pReal)

41
src/thermal_conduction.f90
View File

@ -9,8 +9,7 @@ module thermal_conduction
use lattice use lattice
use results use results
use crystallite use crystallite
use source_thermal_dissipation use constitutive
use source_thermal_externalheat
implicit none implicit none
private private
@ -84,46 +83,14 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
T T
real(pReal), intent(out) :: & real(pReal), intent(out) :: &
Tdot, dTdot_dT Tdot, dTdot_dT
real(pReal) :: &
my_Tdot, my_dTdot_dT
integer :: & integer :: &
phase, & homog
homog, &
offset, &
instance, &
grain, &
source, &
constituent
homog = material_homogenizationAt(el)
offset = material_homogenizationMemberAt(ip,el)
instance = thermal_typeInstance(homog)
Tdot = 0.0_pReal Tdot = 0.0_pReal
dTdot_dT = 0.0_pReal dTdot_dT = 0.0_pReal
do grain = 1, homogenization_Ngrains(homog)
phase = material_phaseAt(grain,el)
constituent = material_phasememberAt(grain,ip,el)
do source = 1, phase_Nsources(phase)
select case(phase_source(source,phase))
case (SOURCE_thermal_dissipation_ID)
call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, &
crystallite_S(1:3,1:3,grain,ip,el), &
crystallite_Lp(1:3,1:3,grain,ip,el), &
phase)
case (SOURCE_thermal_externalheat_ID) homog = material_homogenizationAt(el)
call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, & call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, crystallite_S,crystallite_Lp ,ip, el)
phase, constituent)
case default
my_Tdot = 0.0_pReal
my_dTdot_dT = 0.0_pReal
end select
Tdot = Tdot + my_Tdot
dTdot_dT = dTdot_dT + my_dTdot_dT
enddo
enddo
Tdot = Tdot/real(homogenization_Ngrains(homog),pReal) Tdot = Tdot/real(homogenization_Ngrains(homog),pReal)
dTdot_dT = dTdot_dT/real(homogenization_Ngrains(homog),pReal) dTdot_dT = dTdot_dT/real(homogenization_Ngrains(homog),pReal)