2016-01-09 01:15:20 +05:30
!--------------------------------------------------------------------------------------------------
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine for isotropic (ISOTROPIC) plasticity
!> @details Isotropic (ISOTROPIC) Plasticity which resembles the phenopowerlaw plasticity without
!! resolving the stress on the slip systems. Will give the response of phenopowerlaw for an
!! untextured polycrystal
!--------------------------------------------------------------------------------------------------
module plastic_isotropic
use prec , only : &
pReal , &
2017-05-04 04:02:44 +05:30
pInt
2016-01-09 01:15:20 +05:30
implicit none
private
2016-01-22 13:43:05 +05:30
integer ( pInt ) , dimension ( : ) , allocatable , public , protected :: &
plastic_isotropic_sizePostResults !< cumulative size of post results
integer ( pInt ) , dimension ( : , : ) , allocatable , target , public :: &
plastic_isotropic_sizePostResult !< size of each post result output
character ( len = 64 ) , dimension ( : , : ) , allocatable , target , public :: &
plastic_isotropic_output !< name of each post result output
integer ( pInt ) , dimension ( : ) , allocatable , target , public :: &
plastic_isotropic_Noutput !< number of outputs per instance
2016-01-09 01:15:20 +05:30
enum , bind ( c )
enumerator :: undefined_ID , &
flowstress_ID , &
strainrate_ID
end enum
2016-01-22 06:38:36 +05:30
type , private :: tParameters !< container type for internal constitutive parameters
integer ( kind ( undefined_ID ) ) , allocatable , dimension ( : ) :: &
outputID
real ( pReal ) :: &
2017-05-04 04:02:44 +05:30
fTaylor , &
tau0 , &
gdot0 , &
n , &
h0 , &
2016-04-10 20:22:43 +05:30
h0_slopeLnRate = 0.0_pReal , &
2017-05-04 04:02:44 +05:30
tausat , &
a , &
2016-04-10 20:22:43 +05:30
aTolFlowstress = 1.0_pReal , &
aTolShear = 1.0e-6_pReal , &
tausat_SinhFitA = 0.0_pReal , &
tausat_SinhFitB = 0.0_pReal , &
tausat_SinhFitC = 0.0_pReal , &
tausat_SinhFitD = 0.0_pReal
2016-01-22 06:38:36 +05:30
logical :: &
2016-03-22 20:10:21 +05:30
dilatation = . false .
2016-01-22 06:38:36 +05:30
end type
2018-05-25 03:26:09 +05:30
type ( tParameters ) , dimension ( : ) , allocatable , target , private :: param !< containers of constitutive parameters (len Ninstance)
2016-01-22 06:38:36 +05:30
type , private :: tIsotropicState !< internal state aliases
real ( pReal ) , pointer , dimension ( : ) :: & ! scalars along NipcMyInstance
2016-01-09 01:15:20 +05:30
flowstress , &
2016-01-22 06:38:36 +05:30
accumulatedShear
end type
2017-09-19 05:12:27 +05:30
2016-01-22 06:38:36 +05:30
type , private :: tIsotropicAbsTol !< internal alias for abs tolerance in state
2017-09-30 04:04:18 +05:30
real ( pReal ) , pointer :: & ! scalars
2016-01-22 06:38:36 +05:30
flowstress , &
accumulatedShear
end type
2017-09-19 05:12:27 +05:30
2016-01-22 06:38:36 +05:30
type ( tIsotropicState ) , allocatable , dimension ( : ) , private :: & !< state aliases per instance
state , &
state0 , &
dotState
2017-09-19 05:12:27 +05:30
2016-01-22 06:38:36 +05:30
type ( tIsotropicAbsTol ) , allocatable , dimension ( : ) , private :: & !< state aliases per instance
stateAbsTol
2016-01-09 01:15:20 +05:30
public :: &
plastic_isotropic_init , &
plastic_isotropic_LpAndItsTangent , &
plastic_isotropic_LiAndItsTangent , &
plastic_isotropic_dotState , &
plastic_isotropic_postResults
contains
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
!> @details reads in material parameters, allocates arrays, and does sanity checks
!--------------------------------------------------------------------------------------------------
subroutine plastic_isotropic_init ( fileUnit )
2018-02-02 17:06:09 +05:30
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
2017-10-05 20:05:34 +05:30
use , intrinsic :: iso_fortran_env , only : &
compiler_version , &
compiler_options
#endif
2016-01-09 01:15:20 +05:30
use debug , only : &
debug_level , &
debug_constitutive , &
debug_levelBasic
use numerics , only : &
numerics_integrator
use math , only : &
math_Mandel3333to66 , &
math_Voigt66to3333
use IO , only : &
IO_read , &
IO_lc , &
IO_getTag , &
IO_isBlank , &
IO_stringPos , &
IO_stringValue , &
IO_floatValue , &
IO_error , &
IO_timeStamp , &
IO_EOF
use material , only : &
phase_plasticity , &
phase_plasticityInstance , &
phase_Noutput , &
PLASTICITY_ISOTROPIC_label , &
PLASTICITY_ISOTROPIC_ID , &
material_phase , &
plasticState , &
MATERIAL_partPhase
use lattice
implicit none
integer ( pInt ) , intent ( in ) :: fileUnit
2018-05-25 03:54:58 +05:30
type ( tParameters ) , pointer :: p
2016-01-09 01:15:20 +05:30
integer ( pInt ) , allocatable , dimension ( : ) :: chunkPos
integer ( pInt ) :: &
o , &
phase , &
instance , &
2016-01-22 06:38:36 +05:30
maxNinstance , &
2016-01-09 01:15:20 +05:30
mySize , &
sizeDotState , &
sizeState , &
sizeDeltaState
character ( len = 65536 ) :: &
2016-01-22 06:38:36 +05:30
tag = '' , &
line = '' , &
extmsg = ''
2016-04-25 02:15:33 +05:30
character ( len = 64 ) :: &
outputtag = ''
2017-09-19 05:12:27 +05:30
integer ( pInt ) :: NipcMyPhase
2016-01-09 01:15:20 +05:30
2016-07-25 23:42:00 +05:30
write ( 6 , '(/,a)' ) ' <<<+- constitutive_' / / PLASTICITY_ISOTROPIC_label / / ' init -+>>>'
2018-04-22 13:37:49 +05:30
write ( 6 , '(/,a)' ) ' Ma et al., Computational Materials Science, 109:323– 329, 2015'
2018-04-18 18:31:03 +05:30
write ( 6 , '(/,a)' ) ' https://doi.org/10.1016/j.commatsci.2015.07.041'
2016-07-25 23:42:00 +05:30
write ( 6 , '(a15,a)' ) ' Current time: ' , IO_timeStamp ( )
2016-01-09 01:15:20 +05:30
#include "compilation_info.f90"
maxNinstance = int ( count ( phase_plasticity == PLASTICITY_ISOTROPIC_ID ) , pInt )
if ( maxNinstance == 0_pInt ) return
if ( iand ( debug_level ( debug_constitutive ) , debug_levelBasic ) / = 0_pInt ) &
write ( 6 , '(a16,1x,i5,/)' ) '# instances:' , maxNinstance
2016-01-22 13:43:05 +05:30
allocate ( plastic_isotropic_sizePostResults ( maxNinstance ) , source = 0_pInt )
allocate ( plastic_isotropic_sizePostResult ( maxval ( phase_Noutput ) , maxNinstance ) , source = 0_pInt )
allocate ( plastic_isotropic_output ( maxval ( phase_Noutput ) , maxNinstance ) )
plastic_isotropic_output = ''
allocate ( plastic_isotropic_Noutput ( maxNinstance ) , source = 0_pInt )
2016-01-22 06:38:36 +05:30
allocate ( param ( maxNinstance ) ) ! one container of parameters per instance
2016-03-08 03:42:37 +05:30
2016-01-09 01:15:20 +05:30
rewind ( fileUnit )
phase = 0_pInt
do while ( trim ( line ) / = IO_EOF . and . IO_lc ( IO_getTag ( line , '<' , '>' ) ) / = material_partPhase ) ! wind forward to <phase>
line = IO_read ( fileUnit )
enddo
parsingFile : do while ( trim ( line ) / = IO_EOF ) ! read through sections of phase part
line = IO_read ( fileUnit )
if ( IO_isBlank ( line ) ) cycle ! skip empty lines
if ( IO_getTag ( line , '<' , '>' ) / = '' ) then ! stop at next part
line = IO_read ( fileUnit , . true . ) ! reset IO_read
exit
endif
if ( IO_getTag ( line , '[' , ']' ) / = '' ) then ! next section
phase = phase + 1_pInt ! advance section counter
if ( phase_plasticity ( phase ) == PLASTICITY_ISOTROPIC_ID ) then
2018-05-25 03:26:09 +05:30
p = > param ( phase_plasticityInstance ( phase ) ) ! shorthand pointer to parameter object of my constitutive law
allocate ( p % outputID ( phase_Noutput ( phase ) ) ) ! allocate space for IDs of every requested output
2016-01-09 01:15:20 +05:30
endif
cycle ! skip to next line
endif
2016-01-22 06:38:36 +05:30
if ( phase > 0_pInt ) then ; if ( phase_plasticity ( phase ) == PLASTICITY_ISOTROPIC_ID ) then ! one of my phases. Do not short-circuit here (.and. between if-statements), it's not safe in Fortran
2016-01-09 01:15:20 +05:30
instance = phase_plasticityInstance ( phase ) ! which instance of my plasticity is present phase
2018-05-25 03:26:09 +05:30
p = > param ( instance )
2016-01-09 01:15:20 +05:30
chunkPos = IO_stringPos ( line )
2016-01-22 06:38:36 +05:30
tag = IO_lc ( IO_stringValue ( line , chunkPos , 1_pInt ) ) ! extract key
2016-01-09 01:15:20 +05:30
select case ( tag )
case ( '(output)' )
2016-01-22 06:38:36 +05:30
outputtag = IO_lc ( IO_stringValue ( line , chunkPos , 2_pInt ) )
select case ( outputtag )
2016-01-09 01:15:20 +05:30
case ( 'flowstress' )
2016-01-22 13:43:05 +05:30
plastic_isotropic_Noutput ( instance ) = plastic_isotropic_Noutput ( instance ) + 1_pInt
2018-05-25 03:26:09 +05:30
p % outputID ( plastic_isotropic_Noutput ( instance ) ) = flowstress_ID
2016-01-22 13:43:05 +05:30
plastic_isotropic_output ( plastic_isotropic_Noutput ( instance ) , instance ) = outputtag
2016-01-09 01:15:20 +05:30
case ( 'strainrate' )
2016-01-22 13:43:05 +05:30
plastic_isotropic_Noutput ( instance ) = plastic_isotropic_Noutput ( instance ) + 1_pInt
2018-05-25 03:26:09 +05:30
p % outputID ( plastic_isotropic_Noutput ( instance ) ) = strainrate_ID
2016-01-22 13:43:05 +05:30
plastic_isotropic_output ( plastic_isotropic_Noutput ( instance ) , instance ) = outputtag
2016-01-09 01:15:20 +05:30
end select
2016-01-22 06:38:36 +05:30
2016-01-09 01:15:20 +05:30
case ( '/dilatation/' )
2018-05-25 03:26:09 +05:30
p % dilatation = . true .
2016-01-22 06:38:36 +05:30
2016-01-09 01:15:20 +05:30
case ( 'tau0' )
2018-05-25 03:26:09 +05:30
p % tau0 = IO_floatValue ( line , chunkPos , 2_pInt )
2016-01-22 06:38:36 +05:30
2016-01-09 01:15:20 +05:30
case ( 'gdot0' )
2018-05-25 03:26:09 +05:30
p % gdot0 = IO_floatValue ( line , chunkPos , 2_pInt )
2016-01-22 06:38:36 +05:30
2016-01-09 01:15:20 +05:30
case ( 'n' )
2018-05-25 03:26:09 +05:30
p % n = IO_floatValue ( line , chunkPos , 2_pInt )
2016-01-22 06:38:36 +05:30
2016-01-09 01:15:20 +05:30
case ( 'h0' )
2018-05-25 03:26:09 +05:30
p % h0 = IO_floatValue ( line , chunkPos , 2_pInt )
2016-01-22 06:38:36 +05:30
2016-01-09 01:15:20 +05:30
case ( 'h0_slope' , 'slopelnrate' )
2018-05-25 03:26:09 +05:30
p % h0_slopeLnRate = IO_floatValue ( line , chunkPos , 2_pInt )
2016-01-22 06:38:36 +05:30
2016-01-09 01:15:20 +05:30
case ( 'tausat' )
2018-05-25 03:26:09 +05:30
p % tausat = IO_floatValue ( line , chunkPos , 2_pInt )
2016-01-22 06:38:36 +05:30
2016-01-09 01:15:20 +05:30
case ( 'tausat_sinhfita' )
2018-05-25 03:26:09 +05:30
p % tausat_SinhFitA = IO_floatValue ( line , chunkPos , 2_pInt )
2016-01-22 06:38:36 +05:30
2016-01-09 01:15:20 +05:30
case ( 'tausat_sinhfitb' )
2018-05-25 03:26:09 +05:30
p % tausat_SinhFitB = IO_floatValue ( line , chunkPos , 2_pInt )
2016-01-22 06:38:36 +05:30
2016-01-09 01:15:20 +05:30
case ( 'tausat_sinhfitc' )
2018-05-25 03:26:09 +05:30
p % tausat_SinhFitC = IO_floatValue ( line , chunkPos , 2_pInt )
2016-01-22 06:38:36 +05:30
2016-01-09 01:15:20 +05:30
case ( 'tausat_sinhfitd' )
2018-05-25 03:26:09 +05:30
p % tausat_SinhFitD = IO_floatValue ( line , chunkPos , 2_pInt )
2016-01-22 06:38:36 +05:30
2016-01-09 01:15:20 +05:30
case ( 'a' , 'w0' )
2018-05-25 03:26:09 +05:30
p % a = IO_floatValue ( line , chunkPos , 2_pInt )
2016-01-22 06:38:36 +05:30
2016-01-09 01:15:20 +05:30
case ( 'taylorfactor' )
2018-05-25 03:26:09 +05:30
p % fTaylor = IO_floatValue ( line , chunkPos , 2_pInt )
2016-01-22 06:38:36 +05:30
case ( 'atol_flowstress' )
2018-05-25 03:26:09 +05:30
p % aTolFlowstress = IO_floatValue ( line , chunkPos , 2_pInt )
2016-01-22 06:38:36 +05:30
2016-01-09 01:15:20 +05:30
case ( 'atol_shear' )
2018-05-25 03:26:09 +05:30
p % aTolShear = IO_floatValue ( line , chunkPos , 2_pInt )
2016-01-09 01:15:20 +05:30
case default
end select
endif ; endif
enddo parsingFile
2016-01-22 06:38:36 +05:30
allocate ( state ( maxNinstance ) ) ! internal state aliases
allocate ( state0 ( maxNinstance ) )
allocate ( dotState ( maxNinstance ) )
allocate ( stateAbsTol ( maxNinstance ) )
initializeInstances : do phase = 1_pInt , size ( phase_plasticity ) ! loop over every plasticity
myPhase : if ( phase_plasticity ( phase ) == PLASTICITY_isotropic_ID ) then ! isolate instances of own constitutive description
NipcMyPhase = count ( material_phase == phase ) ! number of own material points (including point components ipc)
2016-01-09 01:15:20 +05:30
instance = phase_plasticityInstance ( phase )
2018-05-25 03:26:09 +05:30
p = > param ( instance )
2017-05-01 02:40:31 +05:30
extmsg = ''
2016-01-09 01:15:20 +05:30
!--------------------------------------------------------------------------------------------------
! sanity checks
2018-05-25 03:26:09 +05:30
if ( p % aTolShear < = 0.0_pReal ) p % aTolShear = 1.0e-6_pReal ! default absolute tolerance 1e-6
if ( p % tau0 < 0.0_pReal ) extmsg = trim ( extmsg ) / / ' tau0'
if ( p % gdot0 < = 0.0_pReal ) extmsg = trim ( extmsg ) / / ' gdot0'
if ( p % n < = 0.0_pReal ) extmsg = trim ( extmsg ) / / ' n'
if ( p % tausat < = 0.0_pReal ) extmsg = trim ( extmsg ) / / ' tausat'
if ( p % a < = 0.0_pReal ) extmsg = trim ( extmsg ) / / ' a'
if ( p % fTaylor < = 0.0_pReal ) extmsg = trim ( extmsg ) / / ' taylorfactor'
if ( p % aTolFlowstress < = 0.0_pReal ) extmsg = trim ( extmsg ) / / ' atol_flowstress'
2017-05-01 02:40:31 +05:30
if ( extmsg / = '' ) then
extmsg = trim ( extmsg ) / / ' (' / / PLASTICITY_ISOTROPIC_label / / ')' ! prepare error message identifier
call IO_error ( 211_pInt , ip = instance , ext_msg = extmsg )
endif
2016-01-09 01:15:20 +05:30
!--------------------------------------------------------------------------------------------------
! Determine size of postResults array
2016-01-22 13:43:05 +05:30
outputsLoop : do o = 1_pInt , plastic_isotropic_Noutput ( instance )
2018-05-25 03:26:09 +05:30
select case ( p % outputID ( o ) )
2016-01-09 01:15:20 +05:30
case ( flowstress_ID , strainrate_ID )
mySize = 1_pInt
case default
end select
outputFound : if ( mySize > 0_pInt ) then
plastic_isotropic_sizePostResult ( o , instance ) = mySize
plastic_isotropic_sizePostResults ( instance ) = &
plastic_isotropic_sizePostResults ( instance ) + mySize
endif outputFound
enddo outputsLoop
!--------------------------------------------------------------------------------------------------
! allocate state arrays
2017-09-19 05:12:27 +05:30
sizeDotState = 2_pInt ! flowstress, accumulated_shear
2016-01-22 06:38:36 +05:30
sizeDeltaState = 0_pInt ! no sudden jumps in state
2017-09-19 05:12:27 +05:30
sizeState = sizeDotState + sizeDeltaState
2016-01-09 01:15:20 +05:30
plasticState ( phase ) % sizeState = sizeState
plasticState ( phase ) % sizeDotState = sizeDotState
plasticState ( phase ) % sizeDeltaState = sizeDeltaState
plasticState ( phase ) % sizePostResults = plastic_isotropic_sizePostResults ( instance )
plasticState ( phase ) % nSlip = 1
plasticState ( phase ) % nTwin = 0
plasticState ( phase ) % nTrans = 0
allocate ( plasticState ( phase ) % aTolState ( sizeState ) )
2016-01-22 06:38:36 +05:30
allocate ( plasticState ( phase ) % state0 ( sizeState , NipcMyPhase ) , source = 0.0_pReal )
allocate ( plasticState ( phase ) % partionedState0 ( sizeState , NipcMyPhase ) , source = 0.0_pReal )
allocate ( plasticState ( phase ) % subState0 ( sizeState , NipcMyPhase ) , source = 0.0_pReal )
allocate ( plasticState ( phase ) % state ( sizeState , NipcMyPhase ) , source = 0.0_pReal )
allocate ( plasticState ( phase ) % dotState ( sizeDotState , NipcMyPhase ) , source = 0.0_pReal )
allocate ( plasticState ( phase ) % deltaState ( sizeDeltaState , NipcMyPhase ) , source = 0.0_pReal )
2016-01-09 01:15:20 +05:30
if ( any ( numerics_integrator == 1_pInt ) ) then
2016-01-22 06:38:36 +05:30
allocate ( plasticState ( phase ) % previousDotState ( sizeDotState , NipcMyPhase ) , source = 0.0_pReal )
allocate ( plasticState ( phase ) % previousDotState2 ( sizeDotState , NipcMyPhase ) , source = 0.0_pReal )
2016-01-09 01:15:20 +05:30
endif
if ( any ( numerics_integrator == 4_pInt ) ) &
2016-01-22 06:38:36 +05:30
allocate ( plasticState ( phase ) % RK4dotState ( sizeDotState , NipcMyPhase ) , source = 0.0_pReal )
2016-01-09 01:15:20 +05:30
if ( any ( numerics_integrator == 5_pInt ) ) &
2016-01-22 06:38:36 +05:30
allocate ( plasticState ( phase ) % RKCK45dotState ( 6 , sizeDotState , NipcMyPhase ) , source = 0.0_pReal )
!--------------------------------------------------------------------------------------------------
! globally required state aliases
plasticState ( phase ) % slipRate = > plasticState ( phase ) % dotState ( 2 : 2 , 1 : NipcMyPhase )
plasticState ( phase ) % accumulatedSlip = > plasticState ( phase ) % state ( 2 : 2 , 1 : NipcMyPhase )
!--------------------------------------------------------------------------------------------------
! locally defined state aliases
state ( instance ) % flowstress = > plasticState ( phase ) % state ( 1 , 1 : NipcMyPhase )
state0 ( instance ) % flowstress = > plasticState ( phase ) % state0 ( 1 , 1 : NipcMyPhase )
dotState ( instance ) % flowstress = > plasticState ( phase ) % dotState ( 1 , 1 : NipcMyPhase )
stateAbsTol ( instance ) % flowstress = > plasticState ( phase ) % aTolState ( 1 )
state ( instance ) % accumulatedShear = > plasticState ( phase ) % state ( 2 , 1 : NipcMyPhase )
state0 ( instance ) % accumulatedShear = > plasticState ( phase ) % state0 ( 2 , 1 : NipcMyPhase )
dotState ( instance ) % accumulatedShear = > plasticState ( phase ) % dotState ( 2 , 1 : NipcMyPhase )
stateAbsTol ( instance ) % accumulatedShear = > plasticState ( phase ) % aTolState ( 2 )
!--------------------------------------------------------------------------------------------------
! init state
2018-05-25 03:26:09 +05:30
state0 ( instance ) % flowstress = p % tau0
2016-01-22 06:38:36 +05:30
state0 ( instance ) % accumulatedShear = 0.0_pReal
!--------------------------------------------------------------------------------------------------
! init absolute state tolerances
2018-05-25 03:26:09 +05:30
stateAbsTol ( instance ) % flowstress = p % aTolFlowstress
stateAbsTol ( instance ) % accumulatedShear = p % aTolShear
2016-01-22 06:38:36 +05:30
2016-01-09 01:15:20 +05:30
endif myPhase
enddo initializeInstances
end subroutine plastic_isotropic_init
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
subroutine plastic_isotropic_LpAndItsTangent ( Lp , dLp_dTstar99 , Tstar_v , ipc , ip , el )
use debug , only : &
debug_level , &
debug_constitutive , &
debug_levelBasic , &
debug_levelExtensive , &
debug_levelSelective , &
debug_e , &
debug_i , &
debug_g
use math , only : &
math_mul6x6 , &
math_Mandel6to33 , &
math_Plain3333to99 , &
math_deviatoric33 , &
math_mul33xx33 , &
math_transpose33
use material , only : &
2016-04-25 02:15:33 +05:30
phasememberAt , &
2016-01-09 01:15:20 +05:30
material_phase , &
phase_plasticityInstance
implicit none
real ( pReal ) , dimension ( 3 , 3 ) , intent ( out ) :: &
Lp !< plastic velocity gradient
real ( pReal ) , dimension ( 9 , 9 ) , intent ( out ) :: &
dLp_dTstar99 !< 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
integer ( pInt ) , intent ( in ) :: &
ipc , & !< component-ID of integration point
ip , & !< integration point
el !< element
2018-05-25 04:01:32 +05:30
type ( tParameters ) , pointer :: p
2018-05-25 03:26:09 +05:30
2016-01-09 01:15:20 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: &
Tstar_dev_33 !< deviatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: &
dLp_dTstar_3333 !< derivative of Lp with respect to Tstar as 4th order tensor
real ( pReal ) :: &
gamma_dot , & !< strainrate
norm_Tstar_dev , & !< euclidean norm of Tstar_dev
squarenorm_Tstar_dev !< square of the euclidean norm of Tstar_dev
integer ( pInt ) :: &
2016-01-22 06:38:36 +05:30
instance , of , &
2016-01-09 01:15:20 +05:30
k , l , m , n
2016-01-22 06:38:36 +05:30
of = phasememberAt ( ipc , ip , el ) ! phasememberAt should be tackled by material and be renamed to material_phasemember
2016-04-25 02:15:33 +05:30
instance = phase_plasticityInstance ( material_phase ( ipc , ip , el ) )
2018-05-25 03:26:09 +05:30
p = > param ( instance )
2016-01-09 01:15:20 +05:30
Tstar_dev_33 = math_deviatoric33 ( math_Mandel6to33 ( Tstar_v ) ) ! deviatoric part of 2nd Piola-Kirchhoff stress
squarenorm_Tstar_dev = math_mul33xx33 ( Tstar_dev_33 , Tstar_dev_33 )
norm_Tstar_dev = sqrt ( squarenorm_Tstar_dev )
if ( norm_Tstar_dev < = 0.0_pReal ) then ! Tstar == 0 --> both Lp and dLp_dTstar are zero
Lp = 0.0_pReal
dLp_dTstar99 = 0.0_pReal
else
2018-05-25 03:26:09 +05:30
gamma_dot = p % gdot0 &
* ( sqrt ( 1.5_pReal ) * norm_Tstar_dev / p % fTaylor / state ( instance ) % flowstress ( of ) ) &
** p % n
2016-01-09 01:15:20 +05:30
2018-05-25 03:26:09 +05:30
Lp = Tstar_dev_33 / norm_Tstar_dev * gamma_dot / p % fTaylor
2016-01-09 01:15:20 +05:30
if ( iand ( debug_level ( debug_constitutive ) , debug_levelExtensive ) / = 0_pInt &
. and . ( ( el == debug_e . and . ip == debug_i . and . ipc == debug_g ) &
. or . . not . iand ( debug_level ( debug_constitutive ) , debug_levelSelective ) / = 0_pInt ) ) then
write ( 6 , '(a,i8,1x,i2,1x,i3)' ) '<< CONST isotropic >> at el ip g ' , el , ip , ipc
write ( 6 , '(/,a,/,3(12x,3(f12.4,1x)/))' ) '<< CONST isotropic >> Tstar (dev) / MPa' , &
math_transpose33 ( Tstar_dev_33 ( 1 : 3 , 1 : 3 ) ) * 1.0e-6_pReal
write ( 6 , '(/,a,/,f12.5)' ) '<< CONST isotropic >> norm Tstar / MPa' , norm_Tstar_dev * 1.0e-6_pReal
write ( 6 , '(/,a,/,f12.5)' ) '<< CONST isotropic >> gdot' , gamma_dot
end if
!--------------------------------------------------------------------------------------------------
! Calculation of the tangent of Lp
forall ( k = 1_pInt : 3_pInt , l = 1_pInt : 3_pInt , m = 1_pInt : 3_pInt , n = 1_pInt : 3_pInt ) &
2018-05-25 03:26:09 +05:30
dLp_dTstar_3333 ( k , l , m , n ) = ( p % n - 1.0_pReal ) * &
2016-01-09 01:15:20 +05:30
Tstar_dev_33 ( k , l ) * Tstar_dev_33 ( m , n ) / squarenorm_Tstar_dev
forall ( k = 1_pInt : 3_pInt , l = 1_pInt : 3_pInt ) &
dLp_dTstar_3333 ( k , l , k , l ) = dLp_dTstar_3333 ( k , l , k , l ) + 1.0_pReal
forall ( k = 1_pInt : 3_pInt , m = 1_pInt : 3_pInt ) &
dLp_dTstar_3333 ( k , k , m , m ) = dLp_dTstar_3333 ( k , k , m , m ) - 1.0_pReal / 3.0_pReal
2018-05-25 03:26:09 +05:30
dLp_dTstar99 = math_Plain3333to99 ( gamma_dot / p % fTaylor * &
2016-01-09 01:15:20 +05:30
dLp_dTstar_3333 / norm_Tstar_dev )
end if
end subroutine plastic_isotropic_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
subroutine plastic_isotropic_LiAndItsTangent ( Li , dLi_dTstar_3333 , Tstar_v , ipc , ip , el )
use math , only : &
math_mul6x6 , &
math_Mandel6to33 , &
math_Plain3333to99 , &
math_spherical33 , &
math_mul33xx33
use material , only : &
2016-03-23 22:47:47 +05:30
phasememberAt , &
2016-01-09 01:15:20 +05:30
material_phase , &
phase_plasticityInstance
implicit none
real ( pReal ) , dimension ( 3 , 3 ) , intent ( out ) :: &
Li !< plastic velocity gradient
2016-04-10 20:22:43 +05:30
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) , intent ( out ) :: &
dLi_dTstar_3333 !< derivative of Li with respect to Tstar as 4th order tensor
2016-01-09 01:15:20 +05:30
real ( pReal ) , dimension ( 6 ) , intent ( in ) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
integer ( pInt ) , intent ( in ) :: &
ipc , & !< component-ID of integration point
ip , & !< integration point
el !< element
2018-05-25 04:01:32 +05:30
type ( tParameters ) , pointer :: p
2018-05-25 03:26:09 +05:30
2016-01-09 01:15:20 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: &
Tstar_sph_33 !< sphiatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor
2018-05-25 03:26:09 +05:30
real ( pReal ) :: &
2016-01-09 01:15:20 +05:30
gamma_dot , & !< strainrate
norm_Tstar_sph , & !< euclidean norm of Tstar_sph
squarenorm_Tstar_sph !< square of the euclidean norm of Tstar_sph
integer ( pInt ) :: &
2016-01-22 06:38:36 +05:30
instance , of , &
2016-01-09 01:15:20 +05:30
k , l , m , n
2016-01-22 06:38:36 +05:30
of = phasememberAt ( ipc , ip , el ) ! phasememberAt should be tackled by material and be renamed to material_phasemember
2016-03-23 22:47:47 +05:30
instance = phase_plasticityInstance ( material_phase ( ipc , ip , el ) )
2018-05-25 03:26:09 +05:30
p = > param ( instance )
2016-01-09 01:15:20 +05:30
Tstar_sph_33 = math_spherical33 ( math_Mandel6to33 ( Tstar_v ) ) ! spherical part of 2nd Piola-Kirchhoff stress
squarenorm_Tstar_sph = math_mul33xx33 ( Tstar_sph_33 , Tstar_sph_33 )
norm_Tstar_sph = sqrt ( squarenorm_Tstar_sph )
2018-05-25 03:26:09 +05:30
if ( p % dilatation . and . norm_Tstar_sph > 0.0_pReal ) then ! Tstar == 0 or J2 plascitiy --> both Li and dLi_dTstar are zero
gamma_dot = p % gdot0 &
* ( sqrt ( 1.5_pReal ) * norm_Tstar_sph / p % fTaylor / state ( instance ) % flowstress ( of ) ) &
** p % n
2016-04-22 15:35:06 +05:30
2018-05-25 03:26:09 +05:30
Li = Tstar_sph_33 / norm_Tstar_sph * gamma_dot / p % fTaylor
2016-04-22 15:35:06 +05:30
!--------------------------------------------------------------------------------------------------
! Calculation of the tangent of Li
forall ( k = 1_pInt : 3_pInt , l = 1_pInt : 3_pInt , m = 1_pInt : 3_pInt , n = 1_pInt : 3_pInt ) &
2018-05-25 03:26:09 +05:30
dLi_dTstar_3333 ( k , l , m , n ) = ( p % n - 1.0_pReal ) * &
2016-04-22 15:35:06 +05:30
Tstar_sph_33 ( k , l ) * Tstar_sph_33 ( m , n ) / squarenorm_Tstar_sph
forall ( k = 1_pInt : 3_pInt , l = 1_pInt : 3_pInt ) &
dLi_dTstar_3333 ( k , l , k , l ) = dLi_dTstar_3333 ( k , l , k , l ) + 1.0_pReal
2018-05-25 03:26:09 +05:30
dLi_dTstar_3333 = gamma_dot / p % fTaylor * &
2016-04-22 15:35:06 +05:30
dLi_dTstar_3333 / norm_Tstar_sph
2016-04-10 20:22:43 +05:30
else
Li = 0.0_pReal
dLi_dTstar_3333 = 0.0_pReal
2016-01-09 09:11:56 +05:30
endif
2018-05-25 04:01:32 +05:30
end subroutine plastic_isotropic_LiAndItsTangent
2016-01-09 01:15:20 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
subroutine plastic_isotropic_dotState ( Tstar_v , ipc , ip , el )
2016-05-29 14:15:03 +05:30
use prec , only : &
2016-10-29 13:09:08 +05:30
dEq0
2016-01-09 01:15:20 +05:30
use math , only : &
math_mul6x6
use material , only : &
2016-03-23 22:47:47 +05:30
phasememberAt , &
2016-01-09 01:15:20 +05:30
material_phase , &
phase_plasticityInstance
implicit none
real ( pReal ) , dimension ( 6 ) , intent ( in ) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
integer ( pInt ) , intent ( in ) :: &
ipc , & !< component-ID of integration point
ip , & !< integration point
el !< element
2018-05-25 04:01:32 +05:30
type ( tParameters ) , pointer :: p
2016-01-09 01:15:20 +05:30
real ( pReal ) , dimension ( 6 ) :: &
Tstar_dev_v !< deviatoric 2nd Piola Kirchhoff stress tensor in Mandel notation
real ( pReal ) :: &
gamma_dot , & !< strainrate
hardening , & !< hardening coefficient
2016-01-22 06:38:36 +05:30
saturation , & !< saturation flowstress
2016-01-09 01:15:20 +05:30
norm_Tstar_v !< euclidean norm of Tstar_dev
integer ( pInt ) :: &
instance , & !< instance of my instance (unique number of my constitutive model)
2016-01-22 06:38:36 +05:30
of !< shortcut notation for offset position in state array
2016-01-09 01:15:20 +05:30
2016-01-22 06:38:36 +05:30
of = phasememberAt ( ipc , ip , el ) ! phasememberAt should be tackled by material and be renamed to material_phasemember
2016-03-23 22:47:47 +05:30
instance = phase_plasticityInstance ( material_phase ( ipc , ip , el ) )
2018-05-25 03:26:09 +05:30
p = > param ( instance )
2016-01-09 01:15:20 +05:30
!--------------------------------------------------------------------------------------------------
! norm of (deviatoric) 2nd Piola-Kirchhoff stress
2018-05-25 03:26:09 +05:30
if ( p % dilatation ) then
2016-01-09 01:15:20 +05:30
norm_Tstar_v = sqrt ( math_mul6x6 ( Tstar_v , Tstar_v ) )
else
Tstar_dev_v ( 1 : 3 ) = Tstar_v ( 1 : 3 ) - sum ( Tstar_v ( 1 : 3 ) ) / 3.0_pReal
Tstar_dev_v ( 4 : 6 ) = Tstar_v ( 4 : 6 )
norm_Tstar_v = sqrt ( math_mul6x6 ( Tstar_dev_v , Tstar_dev_v ) )
end if
!--------------------------------------------------------------------------------------------------
! strain rate
2018-05-25 03:26:09 +05:30
gamma_dot = p % gdot0 * ( sqrt ( 1.5_pReal ) * norm_Tstar_v &
2016-01-09 01:15:20 +05:30
/ & !-----------------------------------------------------------------------------------
2018-05-25 03:26:09 +05:30
( p % fTaylor * state ( instance ) % flowstress ( of ) ) ) ** p % n
2016-01-09 01:15:20 +05:30
!--------------------------------------------------------------------------------------------------
! hardening coefficient
if ( abs ( gamma_dot ) > 1e-12_pReal ) then
2018-05-25 03:26:09 +05:30
if ( dEq0 ( p % tausat_SinhFitA ) ) then
saturation = p % tausat
2016-01-09 01:15:20 +05:30
else
2018-05-25 03:26:09 +05:30
saturation = p % tausat &
+ asinh ( ( gamma_dot / p % tausat_SinhFitA &
) ** ( 1.0_pReal / p % tausat_SinhFitD ) &
) ** ( 1.0_pReal / p % tausat_SinhFitC ) &
/ ( p % tausat_SinhFitB &
* ( gamma_dot / p % gdot0 ) ** ( 1.0_pReal / p % n ) &
2017-05-10 11:10:26 +05:30
)
2016-01-09 01:15:20 +05:30
endif
2018-05-25 03:26:09 +05:30
hardening = ( p % h0 + p % h0_slopeLnRate * log ( gamma_dot ) ) &
* abs ( 1.0_pReal - state ( instance ) % flowstress ( of ) / saturation ) ** p % a &
2016-01-22 06:38:36 +05:30
* sign ( 1.0_pReal , 1.0_pReal - state ( instance ) % flowstress ( of ) / saturation )
2016-01-09 01:15:20 +05:30
else
hardening = 0.0_pReal
endif
2016-01-22 06:38:36 +05:30
dotState ( instance ) % flowstress ( of ) = hardening * gamma_dot
dotState ( instance ) % accumulatedShear ( of ) = gamma_dot
2016-01-09 01:15:20 +05:30
end subroutine plastic_isotropic_dotState
!--------------------------------------------------------------------------------------------------
!> @brief return array of constitutive results
!--------------------------------------------------------------------------------------------------
function plastic_isotropic_postResults ( Tstar_v , ipc , ip , el )
use math , only : &
math_mul6x6
use material , only : &
material_phase , &
2016-03-23 22:47:47 +05:30
phasememberAt , &
2016-01-09 01:15:20 +05:30
phase_plasticityInstance
implicit none
real ( pReal ) , dimension ( 6 ) , intent ( in ) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
integer ( pInt ) , intent ( in ) :: &
ipc , & !< component-ID of integration point
ip , & !< integration point
el !< element
2018-05-25 03:26:09 +05:30
2018-05-25 04:01:32 +05:30
type ( tParameters ) , pointer :: p
2018-05-25 03:26:09 +05:30
2016-01-09 01:15:20 +05:30
real ( pReal ) , dimension ( plastic_isotropic_sizePostResults ( phase_plasticityInstance ( material_phase ( ipc , ip , el ) ) ) ) :: &
plastic_isotropic_postResults
real ( pReal ) , dimension ( 6 ) :: &
Tstar_dev_v !< deviatoric 2nd Piola Kirchhoff stress tensor in Mandel notation
real ( pReal ) :: &
norm_Tstar_v ! euclidean norm of Tstar_dev
integer ( pInt ) :: &
instance , & !< instance of my instance (unique number of my constitutive model)
of , & !< shortcut notation for offset position in state array
c , &
o
2016-01-22 06:38:36 +05:30
of = phasememberAt ( ipc , ip , el ) ! phasememberAt should be tackled by material and be renamed to material_phasemember
2016-03-23 22:47:47 +05:30
instance = phase_plasticityInstance ( material_phase ( ipc , ip , el ) )
2018-05-25 03:26:09 +05:30
p = > param ( instance )
2016-01-09 01:15:20 +05:30
!--------------------------------------------------------------------------------------------------
! norm of (deviatoric) 2nd Piola-Kirchhoff stress
2018-05-25 03:26:09 +05:30
if ( p % dilatation ) then
2016-01-09 01:15:20 +05:30
norm_Tstar_v = sqrt ( math_mul6x6 ( Tstar_v , Tstar_v ) )
else
Tstar_dev_v ( 1 : 3 ) = Tstar_v ( 1 : 3 ) - sum ( Tstar_v ( 1 : 3 ) ) / 3.0_pReal
Tstar_dev_v ( 4 : 6 ) = Tstar_v ( 4 : 6 )
norm_Tstar_v = sqrt ( math_mul6x6 ( Tstar_dev_v , Tstar_dev_v ) )
end if
c = 0_pInt
plastic_isotropic_postResults = 0.0_pReal
2016-01-22 13:43:05 +05:30
outputsLoop : do o = 1_pInt , plastic_isotropic_Noutput ( instance )
2018-05-25 03:26:09 +05:30
select case ( p % outputID ( o ) )
2016-01-09 01:15:20 +05:30
case ( flowstress_ID )
2016-01-22 06:38:36 +05:30
plastic_isotropic_postResults ( c + 1_pInt ) = state ( instance ) % flowstress ( of )
2016-01-09 01:15:20 +05:30
c = c + 1_pInt
case ( strainrate_ID )
plastic_isotropic_postResults ( c + 1_pInt ) = &
2018-05-25 03:26:09 +05:30
p % gdot0 * ( sqrt ( 1.5_pReal ) * norm_Tstar_v &
2016-01-09 01:15:20 +05:30
/ & !----------------------------------------------------------------------------------
2018-05-25 03:26:09 +05:30
( p % fTaylor * state ( instance ) % flowstress ( of ) ) ) ** p % n
2016-01-09 01:15:20 +05:30
c = c + 1_pInt
end select
enddo outputsLoop
end function plastic_isotropic_postResults
end module plastic_isotropic