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

more documentation and parameters capitalization unified and using ipc in all modules (sometimes called gr) 

nonlocal: only missing line continuation in string fixed
This commit is contained in:
Martin Diehl 2013-07-01 06:10:42 +00:00
parent 89cea68bc5
commit 40ace5c666
4 changed files with 564 additions and 441 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

82
code/constitutive_j2.f90
View File

@ -19,9 +19,9 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! $Id$ ! $Id$
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Isotropic (J2) Plasticity !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine for isotropic (J2) plasticity
!> @details Isotropic (J2) Plasticity which resembles the phenopowerlaw plasticity without !> @details Isotropic (J2) Plasticity which resembles the phenopowerlaw plasticity without
!! 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
@ -47,12 +47,12 @@ module constitutive_j2
character(len=64), dimension(:,:), allocatable, target, public :: & character(len=64), dimension(:,:), allocatable, target, public :: &
constitutive_j2_output !< name of each post result output constitutive_j2_output !< name of each post result output
integer(pInt), dimension(:), allocatable, private :: &
constitutive_j2_Noutput !< ??
character(len=32), dimension(:), allocatable, private :: & character(len=32), dimension(:), allocatable, private :: &
constitutive_j2_structureName constitutive_j2_structureName
integer(pInt), dimension(:), allocatable, private :: &
constitutive_j2_Noutput !< ??
real(pReal), dimension(:), allocatable, private :: & real(pReal), dimension(:), allocatable, private :: &
constitutive_j2_fTaylor, & !< Taylor factor constitutive_j2_fTaylor, & !< Taylor factor
constitutive_j2_tau0, & !< initial plastic stress constitutive_j2_tau0, & !< initial plastic stress
@ -72,7 +72,6 @@ module constitutive_j2
constitutive_j2_tausat_SinhFitC, & !< fitting parameter for normalized strain rate vs. stress function constitutive_j2_tausat_SinhFitC, & !< fitting parameter for normalized strain rate vs. stress function
constitutive_j2_tausat_SinhFitD !< fitting parameter for normalized strain rate vs. stress function constitutive_j2_tausat_SinhFitD !< fitting parameter for normalized strain rate vs. stress function
real(pReal), dimension(:,:,:), allocatable, private :: & real(pReal), dimension(:,:,:), allocatable, private :: &
constitutive_j2_Cslip_66 constitutive_j2_Cslip_66
@ -93,6 +92,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief module initialization !> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_j2_init(myFile) subroutine constitutive_j2_init(myFile)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
@ -100,6 +100,7 @@ subroutine constitutive_j2_init(myFile)
math_Mandel3333to66, & math_Mandel3333to66, &
math_Voigt66to3333 math_Voigt66to3333
use IO, only: & use IO, only: &
IO_read, &
IO_lc, & IO_lc, &
IO_getTag, & IO_getTag, &
IO_isBlank, & IO_isBlank, &
@ -107,8 +108,7 @@ subroutine constitutive_j2_init(myFile)
IO_stringValue, & IO_stringValue, &
IO_floatValue, & IO_floatValue, &
IO_error, & IO_error, &
IO_timeStamp, & IO_timeStamp
IO_read
use material use material
use debug, only: & use debug, only: &
debug_level, & debug_level, &
@ -133,12 +133,11 @@ subroutine constitutive_j2_init(myFile)
write(6,'(a15,a)') ' Current time: ',IO_timeStamp() write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
maxNinstance = int(count(phase_plasticity == constitutive_j2_label),pInt) maxNinstance = int(count(phase_plasticity == CONSTITUTIVE_J2_label),pInt)
if (maxNinstance == 0_pInt) return if (maxNinstance == 0_pInt) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
endif
allocate(constitutive_j2_sizeDotState(maxNinstance)) allocate(constitutive_j2_sizeDotState(maxNinstance))
constitutive_j2_sizeDotState = 0_pInt constitutive_j2_sizeDotState = 0_pInt
@ -184,7 +183,6 @@ subroutine constitutive_j2_init(myFile)
constitutive_j2_tausat_SinhFitD = 0.0_pReal constitutive_j2_tausat_SinhFitD = 0.0_pReal
rewind(myFile) rewind(myFile)
do while (trim(line) /= '#EOF#' .and. IO_lc(IO_getTag(line,'<','>')) /= 'phase') ! wind forward to <phase> do while (trim(line) /= '#EOF#' .and. IO_lc(IO_getTag(line,'<','>')) /= 'phase') ! wind forward to <phase>
line = IO_read(myFile) line = IO_read(myFile)
enddo enddo
@ -195,12 +193,12 @@ subroutine constitutive_j2_init(myFile)
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') then ! next section if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1_pInt ! advance section counter section = section + 1_pInt ! advance section counter
cycle cycle ! skip to next line
endif endif
if (section > 0_pInt ) then ! do not short-circuit here (.and. with next if-statement). It's not safe in Fortran if (section > 0_pInt ) then ! do not short-circuit here (.and. with next if-statement). It's not safe in Fortran
if (phase_plasticity(section) == CONSTITUTIVE_J2_label) then ! one of my sections if (phase_plasticity(section) == CONSTITUTIVE_J2_label) then ! one of my sections
i = phase_plasticityInstance(section) ! which instance of my plasticity is present phase i = phase_plasticityInstance(section) ! which instance of my plasticity is present phase
positions = IO_stringPos(line,maxNchunks) positions = IO_stringPos(line,MAXNCHUNKS)
tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key
select case(tag) select case(tag)
case ('plasticity','elasticity') case ('plasticity','elasticity')
@ -256,7 +254,7 @@ subroutine constitutive_j2_init(myFile)
case ('atol_resistance') case ('atol_resistance')
constitutive_j2_aTolResistance(i) = IO_floatValue(line,positions,2_pInt) constitutive_j2_aTolResistance(i) = IO_floatValue(line,positions,2_pInt)
case default case default
call IO_error(210_pInt,ext_msg=trim(tag)//' ('//constitutive_j2_label//')') call IO_error(210_pInt,ext_msg=trim(tag)//' ('//CONSTITUTIVE_J2_label//')')
end select end select
endif endif
endif endif
@ -280,8 +278,8 @@ subroutine constitutive_j2_init(myFile)
//CONSTITUTIVE_J2_label//')') //CONSTITUTIVE_J2_label//')')
enddo sanityChecks enddo sanityChecks
do i = 1_pInt,maxNinstance instancesLoop: do i = 1_pInt,maxNinstance
do o = 1_pInt,constitutive_j2_Noutput(i) outputsLoop: do o = 1_pInt,constitutive_j2_Noutput(i)
select case(constitutive_j2_output(o,i)) select case(constitutive_j2_output(o,i))
case('flowstress') case('flowstress')
mySize = 1_pInt mySize = 1_pInt
@ -296,7 +294,7 @@ subroutine constitutive_j2_init(myFile)
constitutive_j2_sizePostResults(i) = & constitutive_j2_sizePostResults(i) = &
constitutive_j2_sizePostResults(i) + mySize constitutive_j2_sizePostResults(i) + mySize
endif endif
enddo enddo outputsLoop
constitutive_j2_sizeDotState(i) = 1_pInt constitutive_j2_sizeDotState(i) = 1_pInt
constitutive_j2_sizeState(i) = 1_pInt constitutive_j2_sizeState(i) = 1_pInt
@ -306,14 +304,14 @@ subroutine constitutive_j2_init(myFile)
constitutive_j2_Cslip_66(1:6,1:6,i) = & constitutive_j2_Cslip_66(1:6,1:6,i) = &
math_Mandel3333to66(math_Voigt66to3333(constitutive_j2_Cslip_66(1:6,1:6,i))) ! todo what is going on here? math_Mandel3333to66(math_Voigt66to3333(constitutive_j2_Cslip_66(1:6,1:6,i))) ! todo what is going on here?
enddo enddo instancesLoop
end subroutine constitutive_j2_init end subroutine constitutive_j2_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief initial microstructural state !> @brief sets the initial microstructural state for a given instance of this plasticity
!> @detail initial microstructural state is set to the value specified by tau0 !> @details initial microstructural state is set to the value specified by tau0
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function constitutive_j2_stateInit(myInstance) pure function constitutive_j2_stateInit(myInstance)
@ -327,7 +325,7 @@ end function constitutive_j2_stateInit
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief relevant state values for the current instance of this plasticity !> @brief sets the relevant state values for a given instance of this plasticity
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function constitutive_j2_aTolState(myInstance) pure function constitutive_j2_aTolState(myInstance)
@ -343,20 +341,22 @@ end function constitutive_j2_aTolState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief homogenized elasticity matrix !> @brief returns the homogenized elasticity matrix
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function constitutive_j2_homogenizedC(state,ipc,ip,el) pure function constitutive_j2_homogenizedC(state,ipc,ip,el)
use prec, only: & use prec, only: &
p_vec p_vec
use mesh, only: & use mesh, only: &
mesh_NcpElems,mesh_maxNips mesh_NcpElems, &
mesh_maxNips
use material, only: & use material, only: &
homogenization_maxNgrains,& homogenization_maxNgrains,&
material_phase, & material_phase, &
phase_plasticityInstance phase_plasticityInstance
implicit none implicit none
real(pReal), dimension(6,6) :: constitutive_j2_homogenizedC real(pReal), dimension(6,6) :: &
constitutive_j2_homogenizedC
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
@ -371,7 +371,8 @@ end function constitutive_j2_homogenizedC
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculate derived quantities from state (not used here) !> @brief calculates derived quantities from state
!> @details dummy subroutine, does nothing
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine constitutive_j2_microstructure(temperature,state,ipc,ip,el) pure subroutine constitutive_j2_microstructure(temperature,state,ipc,ip,el)
use prec, only: & use prec, only: &
@ -380,9 +381,7 @@ pure subroutine constitutive_j2_microstructure(temperature,state,ipc,ip,el)
mesh_NcpElems, & mesh_NcpElems, &
mesh_maxNips mesh_maxNips
use material, only: & use material, only: &
homogenization_maxNgrains, & homogenization_maxNgrains
material_phase, &
phase_plasticityInstance
implicit none implicit none
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
@ -419,6 +418,11 @@ pure subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar_99,Tstar_v,&
phase_plasticityInstance phase_plasticityInstance
implicit none implicit none
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(9,9), intent(out) :: &
dLp_dTstar_99 !< derivative of Lp with respect to 2nd Piola Kirchhoff stress
real(pReal), dimension(6), intent(in) :: & real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
@ -430,11 +434,6 @@ pure subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar_99,Tstar_v,&
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
state !< microstructure state state !< microstructure state
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(9,9), intent(out) :: &
dLp_dTstar_99 !< derivative of Lp with respect to 2nd Piola Kirchhoff stress
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
Tstar_dev_33 !< deviatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor Tstar_dev_33 !< deviatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor
real(pReal), dimension(3,3,3,3) :: & real(pReal), dimension(3,3,3,3) :: &
@ -479,7 +478,7 @@ end subroutine constitutive_j2_LpAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure !> @brief calculates the rate of change of microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function constitutive_j2_dotState(Tstar_v,Temperature,state,ipc,ip, el) pure function constitutive_j2_dotState(Tstar_v,temperature,state,ipc,ip,el)
use prec, only: & use prec, only: &
p_vec p_vec
use math, only: & use math, only: &
@ -498,7 +497,7 @@ pure function constitutive_j2_dotState(Tstar_v,Temperature,state,ipc,ip, el)
real(pReal), dimension(6), intent(in):: & real(pReal), dimension(6), intent(in):: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
Temperature !< temperature at integration point temperature !< temperature at integration point
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point ipc, & !< component-ID of integration point
ip, & !< integration point ip, & !< integration point
@ -561,13 +560,12 @@ end function constitutive_j2_dotState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief (instantaneous) incremental change of microstructure (dummy function) !> @brief (instantaneous) incremental change of microstructure
!> @details dummy function, returns 0.0
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function constitutive_j2_deltaState(Tstar_v,temperature,state,ipc,ip,el) pure function constitutive_j2_deltaState(Tstar_v,temperature,state,ipc,ip,el)
use prec, only: & use prec, only: &
p_vec p_vec
use math, only: &
math_mul6x6
use mesh, only: & use mesh, only: &
mesh_NcpElems, & mesh_NcpElems, &
mesh_maxNips mesh_maxNips
@ -597,7 +595,8 @@ end function constitutive_j2_deltaState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of temperature (dummy function) !> @brief calculates the rate of change of temperature
!> @details dummy function, returns 0.0
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
real(pReal) pure function constitutive_j2_dotTemperature(Tstar_v,temperature,state,ipc,ip,el) real(pReal) pure function constitutive_j2_dotTemperature(Tstar_v,temperature,state,ipc,ip,el)
use prec, only: & use prec, only: &
@ -653,7 +652,8 @@ pure function constitutive_j2_postResults(Tstar_v,temperature,dt,state,ipc,ip,el
ip, & !< integration point ip, & !< integration point
el !< element el !< element
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
state state !< microstructure state
real(pReal), dimension(constitutive_j2_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & real(pReal), dimension(constitutive_j2_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
constitutive_j2_postResults constitutive_j2_postResults

213
code/constitutive_none.f90
View File

@ -30,19 +30,20 @@ module constitutive_none
implicit none implicit none
private private
character (len=*), parameter, public :: constitutive_none_LABEL = 'none' character (len=*), parameter, public :: &
CONSTITUTIVE_NONE_label = 'none' !< label for this constitutive model
integer(pInt), dimension(:), allocatable, public :: & integer(pInt), dimension(:), allocatable, public :: &
constitutive_none_sizeDotState, & constitutive_none_sizeDotState, &
constitutive_none_sizeState, & constitutive_none_sizeState, &
constitutive_none_sizePostResults constitutive_none_sizePostResults
integer(pInt), dimension(:,:), allocatable, target, public :: &
constitutive_none_sizePostResult !< size of each post result output
character(len=32), dimension(:), allocatable, private :: & character(len=32), dimension(:), allocatable, private :: &
constitutive_none_structureName constitutive_none_structureName
integer(pInt), dimension(:,:), allocatable, target, public :: &
constitutive_none_sizePostResult ! size of each post result output
real(pReal), dimension(:,:,:), allocatable, private :: & real(pReal), dimension(:,:,:), allocatable, private :: &
constitutive_none_Cslip_66 constitutive_none_Cslip_66
@ -62,7 +63,8 @@ module constitutive_none
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief reads in material parameters and allocates arrays !> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_none_init(myFile) subroutine constitutive_none_init(myFile)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
@ -84,23 +86,26 @@ subroutine constitutive_none_init(myFile)
debug_level, & debug_level, &
debug_constitutive, & debug_constitutive, &
debug_levelBasic debug_levelBasic
use lattice, only: lattice_symmetrizeC66 use lattice, only: &
lattice_symmetrizeC66
implicit none implicit none
integer(pInt), intent(in) :: myFile integer(pInt), intent(in) :: myFile
integer(pInt), parameter :: MAXNCHUNKS = 7_pInt integer(pInt), parameter :: MAXNCHUNKS = 7_pInt
integer(pInt), dimension(1_pInt+2_pInt*MAXNCHUNKS) :: positions integer(pInt), dimension(1_pInt+2_pInt*MAXNCHUNKS) :: positions
integer(pInt) :: section = 0_pInt, maxNinstance, i integer(pInt) :: section = 0_pInt, maxNinstance, i
character(len=65536) :: tag character(len=65536) :: &
character(len=65536) :: line = '' ! to start initialized tag = '', &
line = '' ! to start initialized
write(6,'(/,a)') ' <<<+- constitutive_'//trim(constitutive_none_LABEL)//' init -+>>>' write(6,'(/,a)') ' <<<+- constitutive_'//trim(CONSTITUTIVE_NONE_label)//' init -+>>>'
write(6,'(a)') ' $Id$' write(6,'(a)') ' $Id$'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp() write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
maxNinstance = int(count(phase_plasticity == constitutive_none_label),pInt) maxNinstance = int(count(phase_plasticity == CONSTITUTIVE_NONE_label),pInt)
if (maxNinstance == 0_pInt) return if (maxNinstance == 0_pInt) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) & if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
@ -123,7 +128,7 @@ subroutine constitutive_none_init(myFile)
line = IO_read(myFile) line = IO_read(myFile)
enddo enddo
do while (trim(line) /= '#EOF#') ! read thru sections of phase part do while (trim(line) /= '#EOF#') ! read through sections of phase part
line = IO_read(myFile) line = IO_read(myFile)
if (IO_isBlank(line)) cycle ! skip empty lines if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
@ -131,8 +136,8 @@ subroutine constitutive_none_init(myFile)
section = section + 1_pInt ! advance section counter section = section + 1_pInt ! advance section counter
cycle cycle
endif endif
if (section > 0_pInt ) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran if (section > 0_pInt ) then ! do not short-circuit here (.and. with next if-statement). It's not safe in Fortran
if ( phase_plasticity(section) == constitutive_none_LABEL) then ! one of my sections if (phase_plasticity(section) == CONSTITUTIVE_NONE_label) then ! one of my sections
i = phase_plasticityInstance(section) ! which instance of my plasticity is present phase i = phase_plasticityInstance(section) ! which instance of my plasticity is present phase
positions = IO_stringPos(line,MAXNCHUNKS) positions = IO_stringPos(line,MAXNCHUNKS)
tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key
@ -160,7 +165,7 @@ subroutine constitutive_none_init(myFile)
case ('c66') case ('c66')
constitutive_none_Cslip_66(6,6,i) = IO_floatValue(line,positions,2_pInt) constitutive_none_Cslip_66(6,6,i) = IO_floatValue(line,positions,2_pInt)
case default case default
call IO_error(210_pInt,ext_msg=trim(tag)//' ('//constitutive_none_label//')') call IO_error(210_pInt,ext_msg=trim(tag)//' ('//CONSTITUTIVE_NONE_label//')')
end select end select
endif endif
endif endif
@ -170,7 +175,7 @@ subroutine constitutive_none_init(myFile)
if (constitutive_none_structureName(i) == '') call IO_error(205_pInt,e=i) if (constitutive_none_structureName(i) == '') call IO_error(205_pInt,e=i)
enddo enddo
do i = 1_pInt,maxNinstance instancesLoop: do i = 1_pInt,maxNinstance
constitutive_none_sizeDotState(i) = 1_pInt constitutive_none_sizeDotState(i) = 1_pInt
constitutive_none_sizeState(i) = 1_pInt constitutive_none_sizeState(i) = 1_pInt
@ -179,20 +184,20 @@ subroutine constitutive_none_init(myFile)
constitutive_none_Cslip_66(:,:,i) = & constitutive_none_Cslip_66(:,:,i) = &
math_Mandel3333to66(math_Voigt66to3333(constitutive_none_Cslip_66(:,:,i))) math_Mandel3333to66(math_Voigt66to3333(constitutive_none_Cslip_66(:,:,i)))
enddo enddo instancesLoop
end subroutine constitutive_none_init end subroutine constitutive_none_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief sets the initial microstructural state !> @brief sets the initial microstructural state for a given instance of this plasticity
!> @details dummy function, returns 0.0
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function constitutive_none_stateInit(myInstance) pure function constitutive_none_stateInit(myInstance)
implicit none implicit none
integer(pInt), intent(in) :: myInstance
real(pReal), dimension(1) :: constitutive_none_stateInit real(pReal), dimension(1) :: constitutive_none_stateInit
integer(pInt), intent(in) :: myInstance !< number specifying the instance of the plasticity
constitutive_none_stateInit = 0.0_pReal constitutive_none_stateInit = 0.0_pReal
@ -200,22 +205,24 @@ end function constitutive_none_stateInit
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief relevant microstructural state (ensures convergence as state is always 0.0) !> @brief sets the relevant state values for a given instance of this plasticity
!> @details ensures convergence as state is always 0.0
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function constitutive_none_aTolState(myInstance) pure function constitutive_none_aTolState(myInstance)
implicit none implicit none
integer(pInt), intent(in) :: myInstance !< number specifying the current instance of the plasticity integer(pInt), intent(in) :: myInstance !< number specifying the instance of the plasticity
real(pReal), dimension(constitutive_none_sizeState(myInstance)) :: & real(pReal), dimension(constitutive_none_sizeState(myInstance)) :: &
constitutive_none_aTolState !< relevant state values for the current instance of this plasticity constitutive_none_aTolState
constitutive_none_aTolState = 1.0_pReal constitutive_none_aTolState = 1.0_pReal
end function constitutive_none_aTolState end function constitutive_none_aTolState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief homogenized elacticity matrix !> @brief returns the homogenized elasticity matrix
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function constitutive_none_homogenizedC(state,ipc,ip,el) pure function constitutive_none_homogenizedC(state,ipc,ip,el)
use prec, only: & use prec, only: &
@ -229,51 +236,53 @@ pure function constitutive_none_homogenizedC(state,ipc,ip,el)
phase_plasticityInstance phase_plasticityInstance
implicit none implicit none
real(pReal), dimension(6,6) :: constitutive_none_homogenizedC real(pReal), dimension(6,6) :: &
constitutive_none_homogenizedC
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
ipc, & !< component-ID of current integration point ipc, & !< component-ID of integration point
ip, & !< current integration point ip, & !< integration point
el !< current element el !< element
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
integer(pInt) :: matID state !< microstructure state
matID = phase_plasticityInstance(material_phase(ipc,ip,el)) constitutive_none_homogenizedC = constitutive_none_Cslip_66(1:6,1:6,&
constitutive_none_homogenizedC = constitutive_none_Cslip_66(1:6,1:6,matID) phase_plasticityInstance(material_phase(ipc,ip,el)))
end function constitutive_none_homogenizedC end function constitutive_none_homogenizedC
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state (not used here) !> @brief calculates derived quantities from state
!> @details dummy subroutine, does nothing
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine constitutive_none_microstructure(Temperature,state,ipc,ip,el) pure subroutine constitutive_none_microstructure(temperature,state,ipc,ip,el)
use prec, only: & use prec, only: &
p_vec p_vec
use mesh, only: & use mesh, only: &
mesh_NcpElems, & mesh_NcpElems, &
mesh_maxNips mesh_maxNips
use material, only: & use material, only: &
homogenization_maxNgrains, & homogenization_maxNgrains
material_phase, &
phase_plasticityInstance
implicit none implicit none
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
ipc, & !< component-ID of current integration point ipc, & !< component-ID of integration point
ip, & !< current integration point ip, & !< integration point
el !< current element el !< element
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
Temperature !< temperature temperature !< temperature at IP
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
state !< microstructure state
end subroutine constitutive_none_microstructure end subroutine constitutive_none_microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent !> @brief calculates plastic velocity gradient and its tangent
!> @details dummy function, returns 0.0 and Identity
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine constitutive_none_LpAndItsTangent(Lp, dLp_dTstar_99, Tstar_dev_v, Temperature, & pure subroutine constitutive_none_LpAndItsTangent(Lp,dLp_dTstar_99,Tstar_dev_v, &
state, gr, ip, el) temperature, state, ipc, ip, el)
use prec, only: & use prec, only: &
p_vec p_vec
use math, only: & use math, only: &
@ -287,31 +296,35 @@ pure subroutine constitutive_none_LpAndItsTangent(Lp, dLp_dTstar_99, Tstar_dev_v
phase_plasticityInstance phase_plasticityInstance
implicit none implicit none
real(pReal), dimension(6), intent(in) :: Tstar_dev_v !< deviatoric part of the 2nd Piola Kirchhoff stress tensor in Mandel notation real(pReal), dimension(3,3), intent(out) :: &
real(pReal), intent(in) :: Temperature Lp !< plastic velocity gradient
real(pReal), dimension(9,9), intent(out) :: &
dLp_dTstar_99 !< derivative of Lp with respect to 2nd Piola Kirchhoff stress
real(pReal), dimension(6), intent(in) :: &
Tstar_dev_v !< deviatoric part of 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: &
temperature !< temperature at IP
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
gr, & !< grain number ipc, & !< component-ID of integration point
ip, & !< integration point number ip, & !< integration point
el !< element number el !< element
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state !< state of the current microstructure type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
state !< microstructure state
real(pReal), dimension(3,3), intent(out) :: Lp !< plastic velocity gradient Lp = 0.0_pReal ! set Lp to zero
real(pReal), dimension(9,9), intent(out) :: dLp_dTstar_99 !< derivative of Lp with respect to Tstar (9x9 matrix) dLp_dTstar_99 = math_identity2nd(9) ! set dLp_dTstar to Identity
Lp = 0.0_pReal !< set Lp to zero
dLp_dTstar_99 = math_identity2nd(9) !< set dLp_dTstar to Identity
end subroutine constitutive_none_LpAndItsTangent end subroutine constitutive_none_LpAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure !> @brief calculates the rate of change of microstructure
!> @details dummy function, returns 0.0
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function constitutive_none_dotState(Tstar_v, Temperature, state, gr, ip, el) pure function constitutive_none_dotState(Tstar_v,temperature,state,ipc,ip,el)
use prec, only: & use prec, only: &
p_vec p_vec
use math, only: &
math_identity2nd
use mesh, only: & use mesh, only: &
mesh_NcpElems, & mesh_NcpElems, &
mesh_maxNips mesh_maxNips
@ -321,15 +334,18 @@ pure function constitutive_none_dotState(Tstar_v, Temperature, state, gr, ip, el
phase_plasticityInstance phase_plasticityInstance
implicit none implicit none
real(pReal), dimension(1) :: constitutive_none_dotState real(pReal), dimension(1) :: &
constitutive_none_dotState
real(pReal), dimension(6), intent(in) :: Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation real(pReal), dimension(6), intent(in):: &
real(pReal), intent(in) :: Temperature Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: &
temperature !< temperature at integration point
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
gr, & !< grain number ipc, & !< component-ID of integration point
ip, & !< integration point number ip, & !< integration point
el !< element number el !< element
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state !< state of the current microstructure type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
state !< microstructure state
constitutive_none_dotState = 0.0_pReal constitutive_none_dotState = 0.0_pReal
@ -338,12 +354,11 @@ end function constitutive_none_dotState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief (instantaneous) incremental change of microstructure !> @brief (instantaneous) incremental change of microstructure
!> @details dummy function, returns 0.0
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function constitutive_none_deltaState(Tstar_v, Temperature, state, gr, ip, el) function constitutive_none_deltaState(Tstar_v,temperature,state,ipc,ip,el)
use prec, only: & use prec, only: &
p_vec p_vec
use math, only: &
math_identity2nd
use mesh, only: & use mesh, only: &
mesh_NcpElems, & mesh_NcpElems, &
mesh_maxNips mesh_maxNips
@ -353,16 +368,19 @@ function constitutive_none_deltaState(Tstar_v, Temperature, state, gr, ip, el)
phase_plasticityInstance phase_plasticityInstance
implicit none implicit none
real(pReal), dimension(6), intent(in) :: Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation real(pReal), dimension(6), intent(in):: &
real(pReal), intent(in) :: Temperature Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: &
Temperature !< temperature at integration point
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
gr, & !< grain number ipc, & !< component-ID of integration point
ip, & !< integration point number ip, & !< integration point
el !< element number el !< element
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state !< state of the current microstructure type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
state !< microstructure state
real(pReal), dimension(constitutive_none_sizeDotState(phase_plasticityInstance(& real(pReal), dimension(constitutive_none_sizeDotState(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
material_phase(gr,ip,el)))) :: constitutive_none_deltaState constitutive_none_deltaState
constitutive_none_deltaState = 0.0_pReal constitutive_none_deltaState = 0.0_pReal
@ -372,12 +390,11 @@ end function constitutive_none_deltaState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of temperature !> @brief calculates the rate of change of temperature
!> @details dummy function, returns 0.0
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure real(pReal) function constitutive_none_dotTemperature(Tstar_v, Temperature, state, gr, ip, el) real(pReal) pure function constitutive_none_dotTemperature(Tstar_v,temperature,state,ipc,ip,el)
use prec, only: & use prec, only: &
p_vec p_vec
use math, only: &
math_identity2nd
use mesh, only: & use mesh, only: &
mesh_NcpElems, & mesh_NcpElems, &
mesh_maxNips mesh_maxNips
@ -385,13 +402,16 @@ pure real(pReal) function constitutive_none_dotTemperature(Tstar_v, Temperature,
homogenization_maxNgrains homogenization_maxNgrains
implicit none implicit none
real(pReal), dimension(6), intent(in) :: Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation real(pReal), dimension(6), intent(in) :: &
real(pReal), intent(in) :: Temperature Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: &
temperature !< temperature at integration point
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
gr, & !< grain number ipc, & !< component-ID of integration point
ip, & !< integration point number ip, & !< integration point
el !< element number el !< element
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state !< state of the current microstructure type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
state !< microstructure state
constitutive_none_dotTemperature = 0.0_pReal constitutive_none_dotTemperature = 0.0_pReal
@ -400,12 +420,11 @@ end function constitutive_none_dotTemperature
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return array of constitutive results !> @brief return array of constitutive results
!> @details dummy function, returns 0.0
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function constitutive_none_postResults(Tstar_v, Temperature, dt, state, gr, ip, el) pure function constitutive_none_postResults(Tstar_v,temperature,dt,state,ipc,ip,el)
use prec, only: & use prec, only: &
p_vec p_vec
use math, only: &
math_mul6x6
use mesh, only: & use mesh, only: &
mesh_NcpElems, & mesh_NcpElems, &
mesh_maxNips mesh_maxNips
@ -416,18 +435,20 @@ pure function constitutive_none_postResults(Tstar_v, Temperature, dt, state, gr,
phase_Noutput phase_Noutput
implicit none implicit none
real(pReal), dimension(6), intent(in) :: Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
Temperature, & temperature, & !< temperature at integration point
dt !< current time increment dt
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
gr, & !< grain number ipc, & !< component-ID of integration point
ip, & !< integration point number ip, & !< integration point
el !< element number el !< element
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state !< state of the current microstructure type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
state !< microstructure state
real(pReal), dimension(constitutive_none_sizePostResults(phase_plasticityInstance(& real(pReal), dimension(constitutive_none_sizePostResults(phase_plasticityInstance(&
material_phase(gr,ip,el)))) :: constitutive_none_postResults material_phase(ipc,ip,el)))) :: constitutive_none_postResults
constitutive_none_postResults = 0.0_pReal constitutive_none_postResults = 0.0_pReal

4
code/constitutive_nonlocal.f90
View File

@ -937,8 +937,8 @@ do i = 1,maxNinstance
case('dislocationstress') case('dislocationstress')
mySize = 6_pInt mySize = 6_pInt
case default case default
call IO_error(212_pInt,ext_msg=constitutive_nonlocal_output(o,i)//'& call IO_error(212_pInt,ext_msg=constitutive_nonlocal_output(o,i)//&
('//CONSTITUTIVE_NONLOCAL_LABEL//')') '('//CONSTITUTIVE_NONLOCAL_LABEL//')')
end select end select
if (mySize > 0_pInt) then ! any meaningful output found if (mySize > 0_pInt) then ! any meaningful output found

558
code/constitutive_phenopowerlaw.f90
View File

@ -21,16 +21,18 @@
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @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 phenomenological crystal plasticity formulation using a powerlaw fitting !> @brief material subroutine for phenomenological crystal plasticity formulation using a powerlaw
!! fitting
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module constitutive_phenopowerlaw module constitutive_phenopowerlaw
use prec, only: &
use prec, only: pReal,pInt pReal,&
pInt
implicit none implicit none
private private
character (len=*), parameter, public :: & character (len=*), parameter, public :: &
constitutive_phenopowerlaw_LABEL = 'phenopowerlaw' CONSTITUTIVE_PHENOPOWERLAW_label = 'phenopowerlaw'
integer(pInt), dimension(:), allocatable, public :: & integer(pInt), dimension(:), allocatable, public :: &
constitutive_phenopowerlaw_sizeDotState, & constitutive_phenopowerlaw_sizeDotState, &
@ -38,37 +40,31 @@ module constitutive_phenopowerlaw
constitutive_phenopowerlaw_sizePostResults, & !< cumulative size of post results constitutive_phenopowerlaw_sizePostResults, & !< cumulative size of post results
constitutive_phenopowerlaw_structure constitutive_phenopowerlaw_structure
integer(pInt), dimension(:), allocatable, private :: &
constitutive_phenopowerlaw_Noutput, & !< number of outputs per instance of this constitution
constitutive_phenopowerlaw_totalNslip, & !< no. of slip system used in simulation
constitutive_phenopowerlaw_totalNtwin !< no. of twin system used in simulation
integer(pInt), dimension(:,:), allocatable, target, public :: & integer(pInt), dimension(:,:), allocatable, target, public :: &
constitutive_phenopowerlaw_sizePostResult !< size of each post result output constitutive_phenopowerlaw_sizePostResult !< size of each post result output
integer(pInt), dimension(:,:), allocatable, private :: &
constitutive_phenopowerlaw_Nslip, & !< active number of slip systems per family (input parameter, per family)
constitutive_phenopowerlaw_Ntwin !< active number of twin systems per family (input parameter, per family)
character(len=64), dimension(:,:), allocatable, target, public :: & character(len=64), dimension(:,:), allocatable, target, public :: &
constitutive_phenopowerlaw_output !< name of each post result output constitutive_phenopowerlaw_output !< name of each post result output
character(len=32), dimension(:), allocatable, public :: & character(len=32), dimension(:), allocatable, public :: &
constitutive_phenopowerlaw_structureName constitutive_phenopowerlaw_structureName
integer(pInt), dimension(:), allocatable, private :: &
constitutive_phenopowerlaw_Noutput, & !< number of outputs per instance of this constitution
constitutive_phenopowerlaw_totalNslip, & !< no. of slip system used in simulation
constitutive_phenopowerlaw_totalNtwin !< no. of twin system used in simulation
integer(pInt), dimension(:,:), allocatable, private :: &
constitutive_phenopowerlaw_Nslip, & !< active number of slip systems per family (input parameter, per family)
constitutive_phenopowerlaw_Ntwin !< active number of twin systems per family (input parameter, per family)
real(pReal), dimension(:), allocatable, private :: & real(pReal), dimension(:), allocatable, private :: &
constitutive_phenopowerlaw_CoverA, & !< c/a of the crystal (input parameter) constitutive_phenopowerlaw_CoverA, & !< c/a of the crystal (input parameter)
constitutive_phenopowerlaw_gdot0_slip, & !< reference shear strain rate for slip (input parameter) constitutive_phenopowerlaw_gdot0_slip, & !< reference shear strain rate for slip (input parameter)
constitutive_phenopowerlaw_gdot0_twin, & !< reference shear strain rate for twin (input parameter) constitutive_phenopowerlaw_gdot0_twin, & !< reference shear strain rate for twin (input parameter)
constitutive_phenopowerlaw_n_slip, & !< stress exponent for slip (input parameter) constitutive_phenopowerlaw_n_slip, & !< stress exponent for slip (input parameter)
constitutive_phenopowerlaw_n_twin !< stress exponent for twin (input parameter) constitutive_phenopowerlaw_n_twin, & !< stress exponent for twin (input parameter)
real(pReal), dimension(:,:), allocatable, private :: &
constitutive_phenopowerlaw_tau0_slip, & !< initial critical shear stress for slip (input parameter, per family)
constitutive_phenopowerlaw_tau0_twin, & !< initial critical shear stress for twin (input parameter, per family)
constitutive_phenopowerlaw_tausat_slip !< maximum critical shear stress for slip (input parameter, per family)
real(pReal), dimension(:), allocatable, private :: &
constitutive_phenopowerlaw_spr, & !< push-up factor for slip saturation due to twinning constitutive_phenopowerlaw_spr, & !< push-up factor for slip saturation due to twinning
constitutive_phenopowerlaw_twinB, & constitutive_phenopowerlaw_twinB, &
constitutive_phenopowerlaw_twinC, & constitutive_phenopowerlaw_twinC, &
@ -84,6 +80,11 @@ module constitutive_phenopowerlaw
constitutive_phenopowerlaw_aTolTwinfrac constitutive_phenopowerlaw_aTolTwinfrac
real(pReal), dimension(:,:), allocatable, private :: & real(pReal), dimension(:,:), allocatable, private :: &
constitutive_phenopowerlaw_tau0_slip, & !< initial critical shear stress for slip (input parameter, per family)
constitutive_phenopowerlaw_tau0_twin, & !< initial critical shear stress for twin (input parameter, per family)
constitutive_phenopowerlaw_tausat_slip, & !< maximum critical shear stress for slip (input parameter, per family)
constitutive_phenopowerlaw_nonSchmidCoeff, &
constitutive_phenopowerlaw_interaction_SlipSlip, & !< interaction factors slip - slip (input parameter) constitutive_phenopowerlaw_interaction_SlipSlip, & !< interaction factors slip - slip (input parameter)
constitutive_phenopowerlaw_interaction_SlipTwin, & !< interaction factors slip - twin (input parameter) constitutive_phenopowerlaw_interaction_SlipTwin, & !< interaction factors slip - twin (input parameter)
constitutive_phenopowerlaw_interaction_TwinSlip, & !< interaction factors twin - slip (input parameter) constitutive_phenopowerlaw_interaction_TwinSlip, & !< interaction factors twin - slip (input parameter)
@ -96,61 +97,65 @@ module constitutive_phenopowerlaw
constitutive_phenopowerlaw_hardeningMatrix_TwinTwin, & constitutive_phenopowerlaw_hardeningMatrix_TwinTwin, &
constitutive_phenopowerlaw_Cslip_66 constitutive_phenopowerlaw_Cslip_66
real(pReal), dimension(:,:), allocatable, private :: &
constitutive_phenopowerlaw_nonSchmidCoeff
public :: & public :: &
constitutive_phenopowerlaw_init, & constitutive_phenopowerlaw_init, &
constitutive_phenopowerlaw_homogenizedC, & constitutive_phenopowerlaw_stateInit, &
constitutive_phenopowerlaw_aTolState, & constitutive_phenopowerlaw_aTolState, &
constitutive_phenopowerlaw_homogenizedC, &
constitutive_phenopowerlaw_microstructure, &
constitutive_phenopowerlaw_LpAndItsTangent, &
constitutive_phenopowerlaw_dotState, & constitutive_phenopowerlaw_dotState, &
constitutive_phenopowerlaw_deltaState, & constitutive_phenopowerlaw_deltaState, &
constitutive_phenopowerlaw_dotTemperature, & constitutive_phenopowerlaw_dotTemperature, &
constitutive_phenopowerlaw_microstructure, & constitutive_phenopowerlaw_postResults
constitutive_phenopowerlaw_LpAndItsTangent, &
constitutive_phenopowerlaw_postResults, &
constitutive_phenopowerlaw_stateInit
contains contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief reading in parameters from material config and doing consistency checks !> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_phenopowerlaw_init(myFile) subroutine constitutive_phenopowerlaw_init(myFile)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use math, only: math_Mandel3333to66, & use math, only: &
math_Mandel3333to66, &
math_Voigt66to3333 math_Voigt66to3333
use IO use IO
use material use material
use debug, only: debug_level,& use debug, only: &
debug_level, &
debug_constitutive,& debug_constitutive,&
debug_levelBasic debug_levelBasic
use lattice use lattice
implicit none implicit none
integer(pInt), intent(in) :: myFile integer(pInt), intent(in) :: myFile
integer(pInt), parameter :: MAXNCHUNKS = lattice_maxNinteraction + 1_pInt integer(pInt), parameter :: MAXNCHUNKS = lattice_maxNinteraction + 1_pInt
integer(pInt), dimension(1+2*MAXNCHUNKS) :: positions integer(pInt), dimension(1_pInt+2_pInt*MAXNCHUNKS) :: positions
integer(pInt), dimension(6) :: configNchunks integer(pInt), dimension(6) :: configNchunks
integer(pInt) :: section = 0_pInt, maxNinstance, i,j,k, f,o, & integer(pInt) :: &
maxNinstance, &
i,j,k, f,o, &
Nchunks_SlipSlip, Nchunks_SlipTwin, Nchunks_TwinSlip, Nchunks_TwinTwin, & Nchunks_SlipSlip, Nchunks_SlipTwin, Nchunks_TwinSlip, Nchunks_TwinTwin, &
Nchunks_SlipFamilies, Nchunks_TwinFamilies, & Nchunks_SlipFamilies, Nchunks_TwinFamilies, &
mySize=0_pInt, myStructure, index_myFamily, index_otherFamily myStructure, index_myFamily, index_otherFamily, &
character(len=65536) :: tag mySize=0_pInt, section = 0_pInt
character(len=65536) :: line = '' ! to start initialized character(len=65536) :: &
tag = '', &
line = '' ! to start initialized
write(6,'(/,a)') ' <<<+- constitutive_'//trim(constitutive_phenopowerlaw_LABEL)//' init -+>>>' write(6,'(/,a)') ' <<<+- constitutive_'//trim(CONSTITUTIVE_PHENOPOWERLAW_label)//' init -+>>>'
write(6,'(a)') ' $Id$' write(6,'(a)') ' $Id$'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp() write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90" #include "compilation_info.f90"
maxNinstance = int(count(phase_plasticity == constitutive_phenopowerlaw_label),pInt) maxNinstance = int(count(phase_plasticity == CONSTITUTIVE_PHENOPOWERLAW_label),pInt)
if (maxNinstance == 0) return if (maxNinstance == 0_pInt) return
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5)') '# instances:',maxNinstance write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
write(6,*)
endif
Nchunks_SlipFamilies = lattice_maxNslipFamily Nchunks_SlipFamilies = lattice_maxNslipFamily
Nchunks_TwinFamilies = lattice_maxNtwinFamily Nchunks_TwinFamilies = lattice_maxNtwinFamily
@ -239,12 +244,10 @@ subroutine constitutive_phenopowerlaw_init(myFile)
constitutive_phenopowerlaw_nonSchmidCoeff = 0.0_pReal constitutive_phenopowerlaw_nonSchmidCoeff = 0.0_pReal
rewind(myFile) rewind(myFile)
do while (trim(line) /= '#EOF#' .and. IO_lc(IO_getTag(line,'<','>')) /= 'phase') ! wind forward to <phase> do while (trim(line) /= '#EOF#' .and. IO_lc(IO_getTag(line,'<','>')) /= 'phase') ! wind forward to <phase>
line = IO_read(myFile) line = IO_read(myFile)
enddo enddo
do while (trim(line) /= '#EOF#') ! read through sections of phase part
do while (trim(line) /= '#EOF#') ! read thru sections of phase part
line = IO_read(myFile) line = IO_read(myFile)
if (IO_isBlank(line)) cycle ! skip empty lines if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
@ -252,8 +255,8 @@ subroutine constitutive_phenopowerlaw_init(myFile)
section = section + 1_pInt ! advance section counter section = section + 1_pInt ! advance section counter
cycle ! skip to next line cycle ! skip to next line
endif endif
if (section > 0_pInt ) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran if (section > 0_pInt ) then ! do not short-circuit here (.and. with next if-statement). It's not safe in Fortran
if (phase_plasticity(section) == constitutive_phenopowerlaw_LABEL) then ! one of my sections if (phase_plasticity(section) == CONSTITUTIVE_PHENOPOWERLAW_label) then ! one of my sections
i = phase_plasticityInstance(section) ! which instance of my plasticity is present phase i = phase_plasticityInstance(section) ! which instance of my plasticity is present phase
positions = IO_stringPos(line,MAXNCHUNKS) positions = IO_stringPos(line,MAXNCHUNKS)
tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key
@ -337,7 +340,7 @@ subroutine constitutive_phenopowerlaw_init(myFile)
constitutive_phenopowerlaw_h0_SlipSlip(i) = IO_floatValue(line,positions,2_pInt) constitutive_phenopowerlaw_h0_SlipSlip(i) = IO_floatValue(line,positions,2_pInt)
case ('h0_sliptwin') case ('h0_sliptwin')
constitutive_phenopowerlaw_h0_SlipTwin(i) = IO_floatValue(line,positions,2_pInt) constitutive_phenopowerlaw_h0_SlipTwin(i) = IO_floatValue(line,positions,2_pInt)
call IO_warning(42_pInt,ext_msg=trim(tag)//' ('//constitutive_phenopowerlaw_label//')') call IO_warning(42_pInt,ext_msg=trim(tag)//' ('//CONSTITUTIVE_PHENOPOWERLAW_label//')')
case ('h0_twinslip') case ('h0_twinslip')
constitutive_phenopowerlaw_h0_TwinSlip(i) = IO_floatValue(line,positions,2_pInt) constitutive_phenopowerlaw_h0_TwinSlip(i) = IO_floatValue(line,positions,2_pInt)
case ('h0_twintwin') case ('h0_twintwin')
@ -369,14 +372,13 @@ subroutine constitutive_phenopowerlaw_init(myFile)
constitutive_phenopowerlaw_nonSchmidCoeff(j,i) = IO_floatValue(line,positions,1_pInt+j) constitutive_phenopowerlaw_nonSchmidCoeff(j,i) = IO_floatValue(line,positions,1_pInt+j)
enddo enddo
case default case default
call IO_error(210_pInt,ext_msg=trim(tag)//' ('//constitutive_phenopowerlaw_label//')') call IO_error(210_pInt,ext_msg=trim(tag)//' ('//CONSTITUTIVE_PHENOPOWERLAW_label//')')
end select end select
endif endif
endif endif
enddo enddo
do i = 1_pInt,maxNinstance sanityChecks: do i = 1_pInt,maxNinstance
constitutive_phenopowerlaw_structure(i) = lattice_initializeStructure(constitutive_phenopowerlaw_structureName(i), & ! get structure constitutive_phenopowerlaw_structure(i) = lattice_initializeStructure(constitutive_phenopowerlaw_structureName(i), & ! get structure
constitutive_phenopowerlaw_CoverA(i)) constitutive_phenopowerlaw_CoverA(i))
constitutive_phenopowerlaw_Nslip(1:lattice_maxNslipFamily,i) = & constitutive_phenopowerlaw_Nslip(1:lattice_maxNslipFamily,i) = &
@ -391,26 +393,26 @@ subroutine constitutive_phenopowerlaw_init(myFile)
if (constitutive_phenopowerlaw_structure(i) < 1 ) call IO_error(205_pInt,e=i) if (constitutive_phenopowerlaw_structure(i) < 1 ) call IO_error(205_pInt,e=i)
if (any(constitutive_phenopowerlaw_tau0_slip(:,i) < 0.0_pReal .and. & if (any(constitutive_phenopowerlaw_tau0_slip(:,i) < 0.0_pReal .and. &
constitutive_phenopowerlaw_Nslip(:,i) > 0)) call IO_error(211_pInt,e=i,ext_msg='tau0_slip (' & constitutive_phenopowerlaw_Nslip(:,i) > 0)) call IO_error(211_pInt,e=i,ext_msg='tau0_slip (' &
//constitutive_phenopowerlaw_label//')') //CONSTITUTIVE_PHENOPOWERLAW_label//')')
if (constitutive_phenopowerlaw_gdot0_slip(i) <= 0.0_pReal) call IO_error(211_pInt,e=i,ext_msg='gdot0_slip (' & if (constitutive_phenopowerlaw_gdot0_slip(i) <= 0.0_pReal) call IO_error(211_pInt,e=i,ext_msg='gdot0_slip (' &
//constitutive_phenopowerlaw_label//')') //CONSTITUTIVE_PHENOPOWERLAW_label//')')
if (constitutive_phenopowerlaw_n_slip(i) <= 0.0_pReal) call IO_error(211_pInt,e=i,ext_msg='n_slip (' & if (constitutive_phenopowerlaw_n_slip(i) <= 0.0_pReal) call IO_error(211_pInt,e=i,ext_msg='n_slip (' &
//constitutive_phenopowerlaw_label//')') //CONSTITUTIVE_PHENOPOWERLAW_label//')')
if (any(constitutive_phenopowerlaw_tausat_slip(:,i) <= 0.0_pReal .and. & if (any(constitutive_phenopowerlaw_tausat_slip(:,i) <= 0.0_pReal .and. &
constitutive_phenopowerlaw_Nslip(:,i) > 0)) call IO_error(211_pInt,e=i,ext_msg='tausat_slip (' & constitutive_phenopowerlaw_Nslip(:,i) > 0)) call IO_error(211_pInt,e=i,ext_msg='tausat_slip (' &
//constitutive_phenopowerlaw_label//')') //CONSTITUTIVE_PHENOPOWERLAW_label//')')
if (any(constitutive_phenopowerlaw_a_slip(i) == 0.0_pReal .and. & if (any(constitutive_phenopowerlaw_a_slip(i) == 0.0_pReal .and. &
constitutive_phenopowerlaw_Nslip(:,i) > 0)) call IO_error(211_pInt,e=i,ext_msg='a_slip (' & constitutive_phenopowerlaw_Nslip(:,i) > 0)) call IO_error(211_pInt,e=i,ext_msg='a_slip (' &
//constitutive_phenopowerlaw_label//')') //CONSTITUTIVE_PHENOPOWERLAW_label//')')
if (any(constitutive_phenopowerlaw_tau0_twin(:,i) < 0.0_pReal .and. & if (any(constitutive_phenopowerlaw_tau0_twin(:,i) < 0.0_pReal .and. &
constitutive_phenopowerlaw_Ntwin(:,i) > 0)) call IO_error(211_pInt,e=i,ext_msg='tau0_twin (' & constitutive_phenopowerlaw_Ntwin(:,i) > 0)) call IO_error(211_pInt,e=i,ext_msg='tau0_twin (' &
//constitutive_phenopowerlaw_label//')') //CONSTITUTIVE_PHENOPOWERLAW_label//')')
if ( constitutive_phenopowerlaw_gdot0_twin(i) <= 0.0_pReal .and. & if ( constitutive_phenopowerlaw_gdot0_twin(i) <= 0.0_pReal .and. &
any(constitutive_phenopowerlaw_Ntwin(:,i) > 0)) call IO_error(211_pInt,e=i,ext_msg='gdot0_twin (' & any(constitutive_phenopowerlaw_Ntwin(:,i) > 0)) call IO_error(211_pInt,e=i,ext_msg='gdot0_twin (' &
//constitutive_phenopowerlaw_label//')') //CONSTITUTIVE_PHENOPOWERLAW_label//')')
if ( constitutive_phenopowerlaw_n_twin(i) <= 0.0_pReal .and. & if ( constitutive_phenopowerlaw_n_twin(i) <= 0.0_pReal .and. &
any(constitutive_phenopowerlaw_Ntwin(:,i) > 0)) call IO_error(211_pInt,e=i,ext_msg='n_twin (' & any(constitutive_phenopowerlaw_Ntwin(:,i) > 0)) call IO_error(211_pInt,e=i,ext_msg='n_twin (' &
//constitutive_phenopowerlaw_label//')') //CONSTITUTIVE_PHENOPOWERLAW_label//')')
if (constitutive_phenopowerlaw_aTolResistance(i) <= 0.0_pReal) & if (constitutive_phenopowerlaw_aTolResistance(i) <= 0.0_pReal) &
constitutive_phenopowerlaw_aTolResistance(i) = 1.0_pReal ! default absolute tolerance 1 Pa constitutive_phenopowerlaw_aTolResistance(i) = 1.0_pReal ! default absolute tolerance 1 Pa
if (constitutive_phenopowerlaw_aTolShear(i) <= 0.0_pReal) & if (constitutive_phenopowerlaw_aTolShear(i) <= 0.0_pReal) &
@ -418,7 +420,7 @@ subroutine constitutive_phenopowerlaw_init(myFile)
if (constitutive_phenopowerlaw_aTolTwinfrac(i) <= 0.0_pReal) & if (constitutive_phenopowerlaw_aTolTwinfrac(i) <= 0.0_pReal) &
constitutive_phenopowerlaw_aTolTwinfrac(i) = 1.0e-6_pReal ! default absolute tolerance 1e-6 constitutive_phenopowerlaw_aTolTwinfrac(i) = 1.0e-6_pReal ! default absolute tolerance 1e-6
enddo enddo sanityChecks
allocate(constitutive_phenopowerlaw_hardeningMatrix_SlipSlip(maxval(constitutive_phenopowerlaw_totalNslip),& ! slip resistance from slip activity allocate(constitutive_phenopowerlaw_hardeningMatrix_SlipSlip(maxval(constitutive_phenopowerlaw_totalNslip),& ! slip resistance from slip activity
maxval(constitutive_phenopowerlaw_totalNslip),& maxval(constitutive_phenopowerlaw_totalNslip),&
@ -437,8 +439,8 @@ subroutine constitutive_phenopowerlaw_init(myFile)
constitutive_phenopowerlaw_hardeningMatrix_TwinSlip = 0.0_pReal constitutive_phenopowerlaw_hardeningMatrix_TwinSlip = 0.0_pReal
constitutive_phenopowerlaw_hardeningMatrix_TwinTwin = 0.0_pReal constitutive_phenopowerlaw_hardeningMatrix_TwinTwin = 0.0_pReal
do i = 1_pInt,maxNinstance instancesLoop: do i = 1_pInt,maxNinstance
do o = 1_pInt,constitutive_phenopowerlaw_Noutput(i) outputsLoop: do o = 1_pInt,constitutive_phenopowerlaw_Noutput(i)
select case(constitutive_phenopowerlaw_output(o,i)) select case(constitutive_phenopowerlaw_output(o,i))
case('resistance_slip', & case('resistance_slip', &
'shearrate_slip', & 'shearrate_slip', &
@ -457,7 +459,7 @@ subroutine constitutive_phenopowerlaw_init(myFile)
) )
mySize = 1_pInt mySize = 1_pInt
case default case default
call IO_error(212_pInt,ext_msg=constitutive_phenopowerlaw_output(o,i)//' ('//constitutive_phenopowerlaw_label//')') call IO_error(212_pInt,ext_msg=constitutive_phenopowerlaw_output(o,i)//' ('//CONSTITUTIVE_PHENOPOWERLAW_label//')')
end select end select
if (mySize > 0_pInt) then ! any meaningful output found if (mySize > 0_pInt) then ! any meaningful output found
@ -465,7 +467,7 @@ subroutine constitutive_phenopowerlaw_init(myFile)
constitutive_phenopowerlaw_sizePostResults(i) = & constitutive_phenopowerlaw_sizePostResults(i) = &
constitutive_phenopowerlaw_sizePostResults(i) + mySize constitutive_phenopowerlaw_sizePostResults(i) + mySize
endif endif
enddo ! outputs enddo outputsLoop
constitutive_phenopowerlaw_sizeDotState(i) = constitutive_phenopowerlaw_totalNslip(i)+ & constitutive_phenopowerlaw_sizeDotState(i) = constitutive_phenopowerlaw_totalNslip(i)+ &
constitutive_phenopowerlaw_totalNtwin(i)+ & constitutive_phenopowerlaw_totalNtwin(i)+ &
@ -533,23 +535,26 @@ subroutine constitutive_phenopowerlaw_init(myFile)
enddo; enddo enddo; enddo
! report to out file... enddo instancesLoop
enddo
end subroutine constitutive_phenopowerlaw_init end subroutine constitutive_phenopowerlaw_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief initial microstructural state !> @brief sets the initial microstructural state for a given instance of this plasticity
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function constitutive_phenopowerlaw_stateInit(myInstance) pure function constitutive_phenopowerlaw_stateInit(myInstance)
use lattice, only: lattice_maxNslipFamily, lattice_maxNtwinFamily use lattice, only: &
lattice_maxNslipFamily, &
lattice_maxNtwinFamily
implicit none implicit none
integer(pInt), intent(in) :: myInstance integer(pInt), intent(in) :: &
integer(pInt) :: i myInstance !< number specifying the instance of the plasticity
real(pReal), dimension(constitutive_phenopowerlaw_sizeDotState(myInstance)) :: constitutive_phenopowerlaw_stateInit real(pReal), dimension(constitutive_phenopowerlaw_sizeDotState(myInstance)) :: &
constitutive_phenopowerlaw_stateInit
integer(pInt) :: &
i
constitutive_phenopowerlaw_stateInit = 0.0_pReal constitutive_phenopowerlaw_stateInit = 0.0_pReal
@ -572,14 +577,15 @@ end function constitutive_phenopowerlaw_stateInit
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief absolute state tolerance !> @brief sets the relevant state values for a given instance of this plasticity
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function constitutive_phenopowerlaw_aTolState(myInstance) pure function constitutive_phenopowerlaw_aTolState(myInstance)
implicit none
integer(pInt), intent(in) :: myInstance ! number specifying the current instance of the plasticity
real(pReal), dimension(constitutive_phenopowerlaw_sizeState(myInstance)) :: &
constitutive_phenopowerlaw_aTolState ! relevant state values for the current instance of this plasticity
implicit none
integer(pInt), intent(in) :: myInstance !< number specifying the instance of the plasticity
real(pReal), dimension(constitutive_phenopowerlaw_sizeState(myInstance)) :: &
constitutive_phenopowerlaw_aTolState
constitutive_phenopowerlaw_aTolState(1:constitutive_phenopowerlaw_totalNslip(myInstance)+ & constitutive_phenopowerlaw_aTolState(1:constitutive_phenopowerlaw_totalNslip(myInstance)+ &
constitutive_phenopowerlaw_totalNtwin(myInstance)) = & constitutive_phenopowerlaw_totalNtwin(myInstance)) = &
@ -600,33 +606,9 @@ end function constitutive_phenopowerlaw_aTolState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief homogenized elacticity matrix !> @brief returns the homogenized elasticity matrix
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function constitutive_phenopowerlaw_homogenizedC(state,ipc,ip,el) pure function constitutive_phenopowerlaw_homogenizedC(state,ipc,ip,el)
use prec, only: p_vec
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase, phase_plasticityInstance
implicit none
integer(pInt), intent(in) :: &
ipc, & !component-ID of current integration point
ip, & !current integration point
el !current element
integer(pInt) matID
real(pReal), dimension(6,6) :: constitutive_phenopowerlaw_homogenizedC
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
state ! state variables
matID = phase_plasticityInstance(material_phase(ipc,ip,el))
constitutive_phenopowerlaw_homogenizedC = constitutive_phenopowerlaw_Cslip_66(1:6,1:6,matID)
end function constitutive_phenopowerlaw_homogenizedC
!--------------------------------------------------------------------------------------------------
!> @brief calculate derived quantities from state (dummy subroutine, not used here)
!--------------------------------------------------------------------------------------------------
pure subroutine constitutive_phenopowerlaw_microstructure(Temperature,state,ipc,ip,el)
use prec, only: & use prec, only: &
p_vec p_vec
use mesh, only: & use mesh, only: &
@ -638,43 +620,101 @@ pure subroutine constitutive_phenopowerlaw_microstructure(Temperature,state,ipc,
phase_plasticityInstance phase_plasticityInstance
implicit none implicit none
real(pReal), dimension(6,6) :: &
constitutive_phenopowerlaw_homogenizedC
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
ipc, & !component-ID of current integration point ipc, & !< component-ID of integration point
ip, & !current integration point ip, & !< integration point
el !current element el !< element
integer(pInt) :: matID type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
real(pReal), intent(in) :: Temperature ! temperature state !< microstructure state
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state
matID = phase_plasticityInstance(material_phase(ipc,ip,el)) constitutive_phenopowerlaw_homogenizedC = constitutive_phenopowerlaw_Cslip_66(1:6,1:6,&
phase_plasticityInstance(material_phase(ipc,ip,el)))
end function constitutive_phenopowerlaw_homogenizedC
!--------------------------------------------------------------------------------------------------
!> @brief calculates derived quantities from state
!> @details dummy subroutine, does nothing
!--------------------------------------------------------------------------------------------------
pure subroutine constitutive_phenopowerlaw_microstructure(temperature,state,ipc,ip,el)
use prec, only: &
p_vec
use mesh, only: &
mesh_NcpElems, &
mesh_maxNips
use material, only: &
homogenization_maxNgrains
implicit none
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
real(pReal), intent(in) :: &
temperature !< temperature at IP
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
state !< microstructure state
end subroutine constitutive_phenopowerlaw_microstructure end subroutine constitutive_phenopowerlaw_microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief plastic velocity gradient and its tangent !> @brief calculates plastic velocity gradient and its tangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,state,ipc,ip,el) pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar_99,Tstar_v,&
use prec, only: p_vec temperature,state,ipc,ip,el)
use math, only: math_Plain3333to99,math_Mandel6to33 use prec, only: &
use lattice, only: lattice_Sslip,lattice_Sslip_v,lattice_Stwin,lattice_Stwin_v, lattice_maxNslipFamily, lattice_maxNtwinFamily, & p_vec
lattice_NslipSystem,lattice_NtwinSystem,NnonSchmid use math, only: &
use mesh, only: mesh_NcpElems,mesh_maxNips math_Plain3333to99, &
use material, only: homogenization_maxNgrains,material_phase, phase_plasticityInstance math_Mandel6to33
use lattice, only: &
lattice_Sslip, &
lattice_Sslip_v, &
lattice_Stwin, &
lattice_Stwin_v, &
lattice_maxNslipFamily, &
lattice_maxNtwinFamily, &
lattice_NslipSystem, &
lattice_NtwinSystem, &
NnonSchmid
use mesh, only: &
mesh_NcpElems, &
mesh_maxNips
use material, only: &
homogenization_maxNgrains, &
material_phase, &
phase_plasticityInstance
implicit none implicit none
real(pReal), dimension(3,3), intent(out) :: &
Lp !< plastic velocity gradient
real(pReal), dimension(9,9), intent(out) :: &
dLp_dTstar_99 !< derivative of Lp with respect to 2nd Piola Kirchhoff stress
real(pReal), dimension(6), intent(in) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: &
temperature !< temperature at IP
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
ipc, & ! component-ID at current integration point ipc, & !< component-ID of integration point
ip, & ! current integration point ip, & !< integration point
el ! current element el !< element
integer(pInt) matID,nSlip,nTwin,f,i,j,k,l,m,n, structID,index_Gamma,index_F,index_myFamily type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
real(pReal) Temperature state !< microstructure state
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state
real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola Kirchhoff stress tensor (Mandel) integer(pInt) :: &
real(pReal), dimension(3,3), intent(out) :: Lp ! plastic velocity gradient matID, &
real(pReal), dimension(3,3,3,3) :: dLp_dTstar3333 ! derivative of Lp (4th-rank tensor) nSlip, &
real(pReal), dimension(9,9), intent(out) :: dLp_dTstar nTwin,structID,index_Gamma,index_F,index_myFamily, &
real(pReal), dimension(3,3,2) :: nonSchmid_tensor f,i,j,k,l,m,n
real(pReal), dimension(3,3,3,3) :: &
dLp_dTstar3333 !< derivative of Lp with respect to Tstar as 4th order tensor
real(pReal), dimension(3,3,2) :: &
nonSchmid_tensor
real(pReal), dimension(constitutive_phenopowerlaw_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & real(pReal), dimension(constitutive_phenopowerlaw_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg,tau_slip_pos,tau_slip_neg gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg,tau_slip_pos,tau_slip_neg
real(pReal), dimension(constitutive_phenopowerlaw_totalNtwin(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & real(pReal), dimension(constitutive_phenopowerlaw_totalNtwin(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
@ -691,10 +731,10 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temp
Lp = 0.0_pReal Lp = 0.0_pReal
dLp_dTstar3333 = 0.0_pReal dLp_dTstar3333 = 0.0_pReal
dLp_dTstar = 0.0_pReal dLp_dTstar_99 = 0.0_pReal
j = 0_pInt j = 0_pInt
do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families slipFamiliesLoop: do f = 1_pInt,lattice_maxNslipFamily
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,structID)) ! at which index starts my family index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,structID)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family
j = j+1_pInt j = j+1_pInt
@ -742,10 +782,10 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temp
nonSchmid_tensor(m,n,2) nonSchmid_tensor(m,n,2)
endif endif
enddo enddo
enddo enddo slipFamiliesLoop
j = 0_pInt j = 0_pInt
do f = 1_pInt,lattice_maxNtwinFamily ! loop over all twin families twinFamiliesLoop: do f = 1_pInt,lattice_maxNtwinFamily
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,structID)) ! at which index starts my family index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,structID)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family
j = j+1_pInt j = j+1_pInt
@ -769,39 +809,60 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temp
lattice_Stwin(m,n,index_myFamily+i,structID) lattice_Stwin(m,n,index_myFamily+i,structID)
endif endif
enddo enddo
enddo enddo twinFamiliesLoop
dLp_dTstar = math_Plain3333to99(dLp_dTstar3333) dLp_dTstar_99 = math_Plain3333to99(dLp_dTstar3333)
end subroutine constitutive_phenopowerlaw_LpAndItsTangent end subroutine constitutive_phenopowerlaw_LpAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief of change of microstructure, evolution of state variable !> @brief calculates the rate of change of microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el) function constitutive_phenopowerlaw_dotState(Tstar_v,temperature,state,ipc,ip,el)
use prec, only: p_vec use prec, only: &
use lattice, only: lattice_Sslip_v, lattice_Stwin_v, lattice_maxNslipFamily, lattice_maxNtwinFamily, & p_vec
lattice_NslipSystem,lattice_NtwinSystem,lattice_shearTwin,NnonSchmid use lattice, only: &
use mesh, only: mesh_NcpElems,mesh_maxNips lattice_Sslip_v, &
use material, only: homogenization_maxNgrains,material_phase, phase_plasticityInstance lattice_Stwin_v, &
lattice_maxNslipFamily, &
lattice_maxNtwinFamily, &
lattice_NslipSystem, &
lattice_NtwinSystem, &
lattice_shearTwin, &
NnonSchmid
use mesh, only: &
mesh_NcpElems,&
mesh_maxNips
use material, only: &
homogenization_maxNgrains, &
material_phase, &
phase_plasticityInstance
implicit none implicit none
real(pReal), dimension(6), intent(in):: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: &
temperature !< temperature at integration point
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
ipc, & !< component-ID at current integration point ipc, & !< component-ID of integration point
ip, & !< current integration point ip, & !< integration point
el !< current element el !< element
integer(pInt) matID,nSlip,nTwin,f,i,j,k,structID, & type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
state !< microstructure state
real(pReal), dimension(constitutive_phenopowerlaw_sizeDotState(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
constitutive_phenopowerlaw_dotState
integer(pInt) :: matID,nSlip,nTwin,f,i,j,k,structID, &
index_Gamma,index_F,offset_accshear_slip,offset_accshear_twin,index_myFamily index_Gamma,index_F,offset_accshear_slip,offset_accshear_twin,index_myFamily
real(pReal) Temperature,c_SlipSlip,c_SlipTwin,c_TwinSlip,c_TwinTwin, ssat_offset real(pReal) :: c_SlipSlip,c_SlipTwin,c_TwinSlip,c_TwinTwin, ssat_offset
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state
real(pReal), dimension(6), intent(in) :: Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel)
real(pReal), dimension(constitutive_phenopowerlaw_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & real(pReal), dimension(constitutive_phenopowerlaw_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
gdot_slip,tau_slip_pos,tau_slip_neg,left_SlipSlip,left_SlipTwin,right_SlipSlip,right_TwinSlip gdot_slip,tau_slip_pos,tau_slip_neg,left_SlipSlip,left_SlipTwin,right_SlipSlip,right_TwinSlip
real(pReal), dimension(constitutive_phenopowerlaw_totalNtwin(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & real(pReal), dimension(constitutive_phenopowerlaw_totalNtwin(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
gdot_twin,tau_twin,left_TwinSlip,left_TwinTwin,right_SlipTwin,right_TwinTwin gdot_twin,tau_twin,left_TwinSlip,left_TwinTwin,right_SlipTwin,right_TwinTwin
real(pReal), dimension(constitutive_phenopowerlaw_sizeDotState(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
constitutive_phenopowerlaw_dotState
matID = phase_plasticityInstance(material_phase(ipc,ip,el)) matID = phase_plasticityInstance(material_phase(ipc,ip,el))
structID = constitutive_phenopowerlaw_structure(matID) structID = constitutive_phenopowerlaw_structure(matID)
@ -819,19 +880,19 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices ! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices
c_SlipSlip = constitutive_phenopowerlaw_h0_SlipSlip(matID)*& c_SlipSlip = constitutive_phenopowerlaw_h0_SlipSlip(matID)*&
(1.0_pReal + & (1.0_pReal + constitutive_phenopowerlaw_twinC(matID)*state(ipc,ip,el)%p(index_F)**&
constitutive_phenopowerlaw_twinC(matID)*state(ipc,ip,el)%p(index_F)**constitutive_phenopowerlaw_twinB(matID)) constitutive_phenopowerlaw_twinB(matID))
c_SlipTwin = 0.0_pReal c_SlipTwin = 0.0_pReal
c_TwinSlip = constitutive_phenopowerlaw_h0_TwinSlip(matID)*& c_TwinSlip = constitutive_phenopowerlaw_h0_TwinSlip(matID)*&
state(ipc,ip,el)%p(index_Gamma)**constitutive_phenopowerlaw_twinE(matID) state(ipc,ip,el)%p(index_Gamma)**constitutive_phenopowerlaw_twinE(matID)
c_TwinTwin = constitutive_phenopowerlaw_h0_TwinTwin(matID)*& c_TwinTwin = constitutive_phenopowerlaw_h0_TwinTwin(matID)*&
state(ipc,ip,el)%p(index_F)**constitutive_phenopowerlaw_twinD(matID) state(ipc,ip,el)%p(index_F)**constitutive_phenopowerlaw_twinD(matID)
!-- calculate left and right vectors and calculate dot gammas !--------------------------------------------------------------------------------------------------
! calculate left and right vectors and calculate dot gammas
ssat_offset = constitutive_phenopowerlaw_spr(matID)*sqrt(state(ipc,ip,el)%p(index_F)) ssat_offset = constitutive_phenopowerlaw_spr(matID)*sqrt(state(ipc,ip,el)%p(index_F))
j = 0_pInt j = 0_pInt
do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families slipFamiliesLoop1: do f = 1_pInt,lattice_maxNslipFamily
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,structID)) ! at which index starts my family index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,structID)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family
j = j+1_pInt j = j+1_pInt
@ -859,10 +920,10 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el
+(abs(tau_slip_neg(j))/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(matID))& +(abs(tau_slip_neg(j))/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(matID))&
*sign(1.0_pReal,tau_slip_pos(j)) *sign(1.0_pReal,tau_slip_pos(j))
enddo enddo
enddo enddo slipFamiliesLoop1
j = 0_pInt j = 0_pInt
do f = 1_pInt,lattice_maxNtwinFamily ! loop over all twin families twinFamiliesLoop1: do f = 1_pInt,lattice_maxNtwinFamily
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,structID)) ! at which index starts my family index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,structID)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family
j = j+1_pInt j = j+1_pInt
@ -871,20 +932,20 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el
right_SlipTwin(j) = 1.0_pReal ! no system-dependent right part right_SlipTwin(j) = 1.0_pReal ! no system-dependent right part
right_TwinTwin(j) = 1.0_pReal ! no system-dependent right part right_TwinTwin(j) = 1.0_pReal ! no system-dependent right part
!* Calculation of dot vol frac !--------------------------------------------------------------------------------------------------
! Calculation of dot vol frac
tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,structID)) tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,structID))
gdot_twin(j) = (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F gdot_twin(j) = (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F
constitutive_phenopowerlaw_gdot0_twin(matID)*& constitutive_phenopowerlaw_gdot0_twin(matID)*&
(abs(tau_twin(j))/state(ipc,ip,el)%p(nSlip+j))**& (abs(tau_twin(j))/state(ipc,ip,el)%p(nSlip+j))**&
constitutive_phenopowerlaw_n_twin(matID)*max(0.0_pReal,sign(1.0_pReal,tau_twin(j))) constitutive_phenopowerlaw_n_twin(matID)*max(0.0_pReal,sign(1.0_pReal,tau_twin(j)))
enddo enddo
enddo enddo twinFamiliesLoop1
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate the overall hardening based on above ! calculate the overall hardening based on above
j = 0_pInt j = 0_pInt
do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families slipFamiliesLoop2: do f = 1_pInt,lattice_maxNslipFamily
do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family
j = j+1_pInt j = j+1_pInt
constitutive_phenopowerlaw_dotState(j) = & ! evolution of slip resistance j constitutive_phenopowerlaw_dotState(j) = & ! evolution of slip resistance j
@ -898,10 +959,10 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el
abs(gdot_slip(j)) abs(gdot_slip(j))
constitutive_phenopowerlaw_dotState(offset_accshear_slip+j) = abs(gdot_slip(j)) constitutive_phenopowerlaw_dotState(offset_accshear_slip+j) = abs(gdot_slip(j))
enddo enddo
enddo enddo slipFamiliesLoop2
j = 0_pInt j = 0_pInt
do f = 1_pInt,lattice_maxNtwinFamily ! loop over all twin families twinFamiliesLoop2: do f = 1_pInt,lattice_maxNtwinFamily
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,structID)) ! at which index starts my family index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,structID)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family
j = j+1_pInt j = j+1_pInt
@ -916,38 +977,40 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el
gdot_twin(j)/lattice_shearTwin(index_myFamily+i,structID) gdot_twin(j)/lattice_shearTwin(index_myFamily+i,structID)
constitutive_phenopowerlaw_dotState(offset_accshear_twin+j) = abs(gdot_twin(j)) constitutive_phenopowerlaw_dotState(offset_accshear_twin+j) = abs(gdot_twin(j))
enddo enddo
enddo enddo twinFamiliesLoop2
end function constitutive_phenopowerlaw_dotState end function constitutive_phenopowerlaw_dotState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief (instantaneous) incremental change of microstructure !> @brief (instantaneous) incremental change of microstructure
!> @details dummy function, returns 0.0
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function constitutive_phenopowerlaw_deltaState(Tstar_v, Temperature, state, g,ip,el) function constitutive_phenopowerlaw_deltaState(Tstar_v,temperature,state,ipc,ip,el)
use prec, only: &
use prec, only: pReal, &
pInt, &
p_vec p_vec
use mesh, only: mesh_NcpElems, & use mesh, only: &
mesh_NcpElems, &
mesh_maxNips mesh_maxNips
use material, only: homogenization_maxNgrains, & use material, only: &
homogenization_maxNgrains, &
material_phase, & material_phase, &
phase_plasticityInstance phase_plasticityInstance
implicit none implicit none
real(pReal), dimension(6), intent(in):: &
integer(pInt), intent(in) :: g, & ! current grain number Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
ip, & ! current integration point real(pReal), intent(in) :: &
el ! current element number Temperature !< temperature at integration point
real(pReal), intent(in) :: Temperature ! temperature integer(pInt), intent(in) :: &
real(pReal), dimension(6), intent(in) :: Tstar_v ! current 2nd Piola-Kirchhoff stress in Mandel notation ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
state ! current microstructural state state !< microstructure state
real(pReal), dimension(constitutive_phenopowerlaw_sizeDotState(phase_plasticityInstance(material_phase(g,ip,el)))) :: &
constitutive_phenopowerlaw_deltaState ! change of state variables / microstructure
real(pReal), dimension(constitutive_phenopowerlaw_sizeDotState(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
constitutive_phenopowerlaw_deltaState
constitutive_phenopowerlaw_deltaState = 0.0_pReal constitutive_phenopowerlaw_deltaState = 0.0_pReal
@ -955,22 +1018,29 @@ end function constitutive_phenopowerlaw_deltaState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of temperature (dummy function) !> @brief calculates the rate of change of temperature
!> @details dummy function, returns 0.0
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function constitutive_phenopowerlaw_dotTemperature(Tstar_v,Temperature,state,ipc,ip,el) real(pReal) pure function constitutive_phenopowerlaw_dotTemperature(Tstar_v,temperature,state,ipc,ip,el)
use prec, only: pReal,pInt,p_vec use prec, only: &
use mesh, only: mesh_NcpElems, mesh_maxNips p_vec
use material, only: homogenization_maxNgrains use mesh, only: &
mesh_NcpElems, &
mesh_maxNips
use material, only: &
homogenization_maxNgrains
implicit none implicit none
real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola Kirchhoff stress tensor in Mandel notation real(pReal), dimension(6), intent(in) :: &
real(pReal), intent(in) :: Temperature Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
integer(pInt), intent(in):: ipc, & ! grain number real(pReal), intent(in) :: &
ip, & ! integration point number temperature !< temperature at integration point
el ! element number integer(pInt), intent(in) :: &
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state ! state of the current microstructure ipc, & !< component-ID of integration point
ip, & !< integration point
real(pReal) constitutive_phenopowerlaw_dotTemperature ! rate of change of temparature el !< element
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
state
constitutive_phenopowerlaw_dotTemperature = 0.0_pReal constitutive_phenopowerlaw_dotTemperature = 0.0_pReal
@ -980,29 +1050,53 @@ end function constitutive_phenopowerlaw_dotTemperature
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief return array of constitutive results !> @brief return array of constitutive results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,state,ipc,ip,el) pure function constitutive_phenopowerlaw_postResults(Tstar_v,temperature,dt,state,ipc,ip,el)
use prec, only: pReal,pInt,p_vec use prec, only: &
use lattice, only: lattice_Sslip_v,lattice_Stwin_v, lattice_maxNslipFamily, lattice_maxNtwinFamily, & p_vec
lattice_NslipSystem,lattice_NtwinSystem,NnonSchmid use mesh, only: &
use mesh, only: mesh_NcpElems,mesh_maxNips mesh_NcpElems, &
use material, only: homogenization_maxNgrains,material_phase,phase_plasticityInstance,phase_Noutput mesh_maxNips
use material, only: &
homogenization_maxNgrains, &
material_phase, &
phase_plasticityInstance, &
phase_Noutput
use lattice, only: &
lattice_Sslip_v, &
lattice_Stwin_v, &
lattice_maxNslipFamily, &
lattice_maxNtwinFamily, &
lattice_NslipSystem, &
lattice_NtwinSystem,NnonSchmid
use mesh, only: &
mesh_NcpElems, &
mesh_maxNips
implicit none implicit none
integer(pInt), intent(in) :: & real(pReal), dimension(6), intent(in) :: &
ipc, & !component-ID at current integration point Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
ip, & !current integration point
el !current element
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
dt, & !current time increment temperature, & !< temperature at integration point
Temperature dt
real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola Kirchhoff stress tensor (Mandel) integer(pInt), intent(in) :: &
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state ipc, & !< component-ID of integration point
integer(pInt) matID,o,f,i,c,nSlip,nTwin,j,k,structID, & ip, & !< integration point
index_Gamma,index_F,index_accshear_slip,index_accshear_twin,index_myFamily el !< element
real(pReal) tau_slip_pos,tau_slip_neg,tau type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
state !< microstructure state
real(pReal), dimension(constitutive_phenopowerlaw_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & real(pReal), dimension(constitutive_phenopowerlaw_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
constitutive_phenopowerlaw_postResults constitutive_phenopowerlaw_postResults
integer(pInt) :: &
matID,structID, &
nSlip,nTwin, &
o,f,i,c,j,k, &
index_Gamma,index_F,index_accshear_slip,index_accshear_twin,index_myFamily
real(pReal) :: &
tau_slip_pos,tau_slip_neg,tau
matID = phase_plasticityInstance(material_phase(ipc,ip,el)) matID = phase_plasticityInstance(material_phase(ipc,ip,el))
structID = constitutive_phenopowerlaw_structure(matID) structID = constitutive_phenopowerlaw_structure(matID)
@ -1017,7 +1111,7 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,stat
constitutive_phenopowerlaw_postResults = 0.0_pReal constitutive_phenopowerlaw_postResults = 0.0_pReal
c = 0_pInt c = 0_pInt
do o = 1_pInt,phase_Noutput(material_phase(ipc,ip,el)) outputsLoop: do o = 1_pInt,phase_Noutput(material_phase(ipc,ip,el))
select case(constitutive_phenopowerlaw_output(o,matID)) select case(constitutive_phenopowerlaw_output(o,matID))
case ('resistance_slip') case ('resistance_slip')
constitutive_phenopowerlaw_postResults(c+1_pInt:c+nSlip) = state(ipc,ip,el)%p(1:nSlip) constitutive_phenopowerlaw_postResults(c+1_pInt:c+nSlip) = state(ipc,ip,el)%p(1:nSlip)
@ -1030,7 +1124,7 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,stat
case ('shearrate_slip') case ('shearrate_slip')
j = 0_pInt j = 0_pInt
do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families slipFamiliesLoop1: do f = 1_pInt,lattice_maxNslipFamily
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,structID)) ! at which index starts my family index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,structID)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family
j = j + 1_pInt j = j + 1_pInt
@ -1046,35 +1140,40 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,stat
((abs(tau_slip_pos)/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(matID) & ((abs(tau_slip_pos)/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(matID) &
+(abs(tau_slip_neg)/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(matID))& +(abs(tau_slip_neg)/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(matID))&
*sign(1.0_pReal,tau_slip_pos) *sign(1.0_pReal,tau_slip_pos)
enddo; enddo enddo
enddo slipFamiliesLoop1
c = c + nSlip c = c + nSlip
case ('resolvedstress_slip') case ('resolvedstress_slip')
j = 0_pInt j = 0_pInt
do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families slipFamiliesLoop2: do f = 1_pInt,lattice_maxNslipFamily
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,structID)) ! at which index starts my family index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,structID)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family
j = j + 1_pInt j = j + 1_pInt
constitutive_phenopowerlaw_postResults(c+j) = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,structID)) constitutive_phenopowerlaw_postResults(c+j) = &
enddo; enddo dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,structID))
enddo
enddo slipFamiliesLoop2
c = c + nSlip c = c + nSlip
case ('totalshear') case ('totalshear')
constitutive_phenopowerlaw_postResults(c+1_pInt) = state(ipc,ip,el)%p(index_Gamma) constitutive_phenopowerlaw_postResults(c+1_pInt) = &
state(ipc,ip,el)%p(index_Gamma)
c = c + 1_pInt c = c + 1_pInt
case ('resistance_twin') case ('resistance_twin')
constitutive_phenopowerlaw_postResults(c+1_pInt:c+nTwin) = state(ipc,ip,el)%p(1_pInt+nSlip:nTwin+nSlip) constitutive_phenopowerlaw_postResults(c+1_pInt:c+nTwin) = &
state(ipc,ip,el)%p(1_pInt+nSlip:nTwin+nSlip)
c = c + nTwin c = c + nTwin
case ('accumulatedshear_twin') case ('accumulatedshear_twin')
constitutive_phenopowerlaw_postResults(c+1_pInt:c+nTwin) = state(ipc,ip,el)%p(index_accshear_twin:& constitutive_phenopowerlaw_postResults(c+1_pInt:c+nTwin) = &
index_accshear_twin+nTwin) state(ipc,ip,el)%p(index_accshear_twin:index_accshear_twin+nTwin)
c = c + nTwin c = c + nTwin
case ('shearrate_twin') case ('shearrate_twin')
j = 0_pInt j = 0_pInt
do f = 1_pInt,lattice_maxNtwinFamily ! loop over all twin families twinFamiliesLoop1: do f = 1_pInt,lattice_maxNtwinFamily
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,structID)) ! at which index starts my family index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,structID)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family
j = j + 1_pInt j = j + 1_pInt
@ -1083,17 +1182,20 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,stat
constitutive_phenopowerlaw_gdot0_twin(matID)*& constitutive_phenopowerlaw_gdot0_twin(matID)*&
(abs(tau)/state(ipc,ip,el)%p(j+nSlip))**& (abs(tau)/state(ipc,ip,el)%p(j+nSlip))**&
constitutive_phenopowerlaw_n_twin(matID)*max(0.0_pReal,sign(1.0_pReal,tau)) constitutive_phenopowerlaw_n_twin(matID)*max(0.0_pReal,sign(1.0_pReal,tau))
enddo; enddo enddo
enddo twinFamiliesLoop1
c = c + nTwin c = c + nTwin
case ('resolvedstress_twin') case ('resolvedstress_twin')
j = 0_pInt j = 0_pInt
do f = 1_pInt,lattice_maxNtwinFamily ! loop over all twin families twinFamiliesLoop2: do f = 1_pInt,lattice_maxNtwinFamily
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,structID)) ! at which index starts my family index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,structID)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family
j = j + 1_pInt j = j + 1_pInt
constitutive_phenopowerlaw_postResults(c+j) = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,structID)) constitutive_phenopowerlaw_postResults(c+j) = &
enddo; enddo dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,structID))
enddo
enddo twinFamiliesLoop2
c = c + nTwin c = c + nTwin
case ('totalvolfrac') case ('totalvolfrac')
@ -1101,7 +1203,7 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,stat
c = c + 1_pInt c = c + 1_pInt
end select end select
enddo enddo outputsLoop
end function constitutive_phenopowerlaw_postResults end function constitutive_phenopowerlaw_postResults