2013-06-29 00:28:10 +05:30
!--------------------------------------------------------------------------------------------------
! $Id$
!--------------------------------------------------------------------------------------------------
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
2013-07-01 11:40:42 +05:30
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief material subroutine for isotropic (J2) plasticity
2013-06-29 00:28:10 +05:30
!> @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
!! untextured polycrystal
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
module constitutive_j2
2014-04-15 16:19:24 +05:30
#ifdef HDF
use hdf5 , only : &
HID_T
#endif
2013-06-29 00:28:10 +05:30
use prec , only : &
pReal , &
pInt
2012-03-09 01:55:28 +05:30
2009-03-06 15:43:08 +05:30
implicit none
2012-03-09 01:55:28 +05:30
private
2013-12-12 03:33:09 +05:30
integer ( pInt ) , dimension ( : ) , allocatable , public , protected :: &
2014-05-22 16:07:57 +05:30
#ifndef NEWSTATE
2013-10-08 21:57:26 +05:30
constitutive_j2_sizeDotState , & !< number of dotStates
constitutive_j2_sizeState , & !< total number of microstructural variables
2014-05-22 16:07:57 +05:30
#endif
2013-10-08 21:57:26 +05:30
constitutive_j2_sizePostResults !< cumulative size of post results
2012-03-09 01:55:28 +05:30
2013-12-12 03:33:09 +05:30
integer ( pInt ) , dimension ( : , : ) , allocatable , target , public :: &
2013-06-29 00:28:10 +05:30
constitutive_j2_sizePostResult !< size of each post result output
2012-03-09 01:55:28 +05:30
2013-12-12 03:33:09 +05:30
character ( len = 64 ) , dimension ( : , : ) , allocatable , target , public :: &
2013-06-29 00:28:10 +05:30
constitutive_j2_output !< name of each post result output
2009-03-06 15:43:08 +05:30
2013-12-12 03:33:09 +05:30
integer ( pInt ) , dimension ( : ) , allocatable , private :: &
2013-10-08 21:57:26 +05:30
constitutive_j2_Noutput !< number of outputs per instance
2013-12-12 03:33:09 +05:30
real ( pReal ) , dimension ( : ) , allocatable , private :: &
2013-06-29 00:28:10 +05:30
constitutive_j2_fTaylor , & !< Taylor factor
constitutive_j2_tau0 , & !< initial plastic stress
constitutive_j2_gdot0 , & !< reference velocity
constitutive_j2_n , & !< Visco-plastic parameter
!--------------------------------------------------------------------------------------------------
! h0 as function of h0 = A + B log (gammadot)
2012-03-09 01:55:28 +05:30
constitutive_j2_h0 , &
2012-04-11 19:31:02 +05:30
constitutive_j2_h0_slopeLnRate , &
2013-06-29 00:28:10 +05:30
constitutive_j2_tausat , & !< final plastic stress
2012-03-09 01:55:28 +05:30
constitutive_j2_a , &
2012-04-11 19:31:02 +05:30
constitutive_j2_aTolResistance , &
2013-06-29 00:28:10 +05:30
!--------------------------------------------------------------------------------------------------
! tausat += (asinh((gammadot / SinhFitA)**(1 / SinhFitD)))**(1 / SinhFitC) / (SinhFitB * (gammadot / gammadot0)**(1/n))
constitutive_j2_tausat_SinhFitA , & !< fitting parameter for normalized strain rate vs. stress function
constitutive_j2_tausat_SinhFitB , & !< 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
2013-12-12 03:33:09 +05:30
enum , bind ( c )
enumerator :: undefined_ID , &
flowstress_ID , &
strainrate_ID
end enum
integer ( kind ( undefined_ID ) ) , dimension ( : , : ) , allocatable , private :: &
constitutive_j2_outputID !< ID of each post result output
2014-04-15 16:19:24 +05:30
#ifdef HDF
type constitutive_j2_tOutput
real ( pReal ) , dimension ( : ) , allocatable , private :: &
flowstress , &
strainrate
logical :: flowstressActive = . false . , strainrateActive = . false . ! if we can write the output block wise, this is not needed anymore because we can do an if(allocated(xxx))
end type constitutive_j2_tOutput
type ( constitutive_j2_tOutput ) , allocatable , dimension ( : ) :: constitutive_j2_Output2
integer ( HID_T ) , allocatable , dimension ( : ) :: outID
#endif
2013-06-29 00:28:10 +05:30
public :: &
constitutive_j2_init , &
2014-05-22 16:07:57 +05:30
#ifndef NEWSTATE
2013-06-29 00:28:10 +05:30
constitutive_j2_stateInit , &
constitutive_j2_aTolState , &
2014-05-22 16:07:57 +05:30
#endif
2013-06-29 00:28:10 +05:30
constitutive_j2_LpAndItsTangent , &
constitutive_j2_dotState , &
constitutive_j2_postResults
2012-03-09 01:55:28 +05:30
contains
2013-06-29 00:28:10 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief module initialization
2013-07-01 11:40:42 +05:30
!> @details reads in material parameters, allocates arrays, and does sanity checks
2013-06-29 00:28:10 +05:30
!--------------------------------------------------------------------------------------------------
2013-12-12 03:33:09 +05:30
subroutine constitutive_j2_init ( fileUnit )
2013-06-29 00:28:10 +05:30
use , intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
2014-04-15 16:19:24 +05:30
#ifdef HDF
use hdf5
#endif
2014-05-08 20:25:19 +05:30
use debug , only : &
2013-12-19 14:19:47 +05:30
debug_level , &
debug_constitutive , &
debug_levelBasic
2014-05-08 20:25:19 +05:30
use numerics , only : &
numerics_integrator
2012-04-11 19:31:02 +05:30
use math , only : &
math_Mandel3333to66 , &
math_Voigt66to3333
use IO , only : &
2013-07-01 11:40:42 +05:30
IO_read , &
2012-04-11 19:31:02 +05:30
IO_lc , &
IO_getTag , &
IO_isBlank , &
IO_stringPos , &
IO_stringValue , &
IO_floatValue , &
2013-02-25 22:04:59 +05:30
IO_error , &
2013-12-12 03:33:09 +05:30
IO_timeStamp , &
2014-04-15 16:19:24 +05:30
#ifdef HDF
tempResults , &
HDF5_addGroup , &
HDF5_addScalarDataset , &
#endif
2013-12-12 03:33:09 +05:30
IO_EOF
2013-12-19 14:19:47 +05:30
use material , only : &
homogenization_maxNgrains , &
phase_plasticity , &
phase_plasticityInstance , &
phase_Noutput , &
PLASTICITY_J2_label , &
PLASTICITY_J2_ID , &
2014-05-08 20:25:19 +05:30
material_phase , &
#ifdef NEWSTATE
plasticState , &
#endif
2013-12-19 14:19:47 +05:30
MATERIAL_partPhase
2014-05-08 20:25:19 +05:30
2013-11-27 21:50:27 +05:30
use lattice
2012-03-09 01:55:28 +05:30
implicit none
2013-12-12 03:33:09 +05:30
integer ( pInt ) , intent ( in ) :: fileUnit
2012-03-09 01:55:28 +05:30
2013-06-29 00:28:10 +05:30
integer ( pInt ) , parameter :: MAXNCHUNKS = 7_pInt
2012-03-09 01:55:28 +05:30
2013-06-29 00:28:10 +05:30
integer ( pInt ) , dimension ( 1_pInt + 2_pInt * MAXNCHUNKS ) :: positions
2014-04-15 16:19:24 +05:30
integer ( pInt ) :: phase , maxNinstance , instance , o , mySize , myConstituents
2013-06-29 00:28:10 +05:30
character ( len = 65536 ) :: &
tag = '' , &
2013-12-12 22:39:59 +05:30
line = ''
2014-05-08 20:25:19 +05:30
integer ( pInt ) :: NofMyPhase
2014-04-15 16:19:24 +05:30
#ifdef HDF
character ( len = 5 ) :: &
str1
integer ( HID_T ) :: ID , ID2 , ID4
#endif
2014-05-08 20:25:19 +05:30
#ifdef NEWSTATE
integer ( pInt ) :: sizeDotState , sizeState
#endif
2013-11-27 13:34:05 +05:30
write ( 6 , '(/,a)' ) ' <<<+- constitutive_' / / PLASTICITY_J2_label / / ' init -+>>>'
2013-05-28 23:01:55 +05:30
write ( 6 , '(a)' ) ' $Id$'
write ( 6 , '(a15,a)' ) ' Current time: ' , IO_timeStamp ( )
2012-02-01 00:48:55 +05:30
#include "compilation_info.f90"
2009-03-06 15:43:08 +05:30
2013-11-27 13:34:05 +05:30
maxNinstance = int ( count ( phase_plasticity == PLASTICITY_J2_ID ) , pInt )
2012-02-21 21:30:00 +05:30
if ( maxNinstance == 0_pInt ) return
2009-03-06 15:43:08 +05:30
2014-04-15 16:19:24 +05:30
#ifdef HDF
allocate ( constitutive_j2_Output2 ( maxNinstance ) )
allocate ( outID ( maxNinstance ) )
#endif
2013-07-01 11:40:42 +05:30
if ( iand ( debug_level ( debug_constitutive ) , debug_levelBasic ) / = 0_pInt ) &
2013-06-29 00:28:10 +05:30
write ( 6 , '(a16,1x,i5,/)' ) '# instances:' , maxNinstance
2014-05-22 16:07:57 +05:30
#ifndef NEWSTATE
2013-12-20 14:06:15 +05:30
allocate ( constitutive_j2_sizeDotState ( maxNinstance ) , source = 1_pInt )
allocate ( constitutive_j2_sizeState ( maxNinstance ) , source = 1_pInt )
2014-05-22 16:07:57 +05:30
#endif
2013-12-12 03:33:09 +05:30
allocate ( constitutive_j2_sizePostResults ( maxNinstance ) , source = 0_pInt )
allocate ( constitutive_j2_sizePostResult ( maxval ( phase_Noutput ) , maxNinstance ) , source = 0_pInt )
2012-03-09 01:55:28 +05:30
allocate ( constitutive_j2_output ( maxval ( phase_Noutput ) , maxNinstance ) )
constitutive_j2_output = ''
2013-12-12 03:33:09 +05:30
allocate ( constitutive_j2_outputID ( maxval ( phase_Noutput ) , maxNinstance ) , source = undefined_ID )
allocate ( constitutive_j2_Noutput ( maxNinstance ) , source = 0_pInt )
allocate ( constitutive_j2_fTaylor ( maxNinstance ) , source = 0.0_pReal )
allocate ( constitutive_j2_tau0 ( maxNinstance ) , source = 0.0_pReal )
allocate ( constitutive_j2_gdot0 ( maxNinstance ) , source = 0.0_pReal )
allocate ( constitutive_j2_n ( maxNinstance ) , source = 0.0_pReal )
allocate ( constitutive_j2_h0 ( maxNinstance ) , source = 0.0_pReal )
allocate ( constitutive_j2_h0_slopeLnRate ( maxNinstance ) , source = 0.0_pReal )
allocate ( constitutive_j2_tausat ( maxNinstance ) , source = 0.0_pReal )
allocate ( constitutive_j2_a ( maxNinstance ) , source = 0.0_pReal )
allocate ( constitutive_j2_aTolResistance ( maxNinstance ) , source = 0.0_pReal )
allocate ( constitutive_j2_tausat_SinhFitA ( maxNinstance ) , source = 0.0_pReal )
allocate ( constitutive_j2_tausat_SinhFitB ( maxNinstance ) , source = 0.0_pReal )
allocate ( constitutive_j2_tausat_SinhFitC ( maxNinstance ) , source = 0.0_pReal )
allocate ( constitutive_j2_tausat_SinhFitD ( maxNinstance ) , source = 0.0_pReal )
2009-03-06 15:43:08 +05:30
2013-12-12 03:33:09 +05:30
rewind ( fileUnit )
2014-03-09 02:20:31 +05:30
phase = 0_pInt
2013-12-19 14:19:47 +05:30
do while ( trim ( line ) / = IO_EOF . and . IO_lc ( IO_getTag ( line , '<' , '>' ) ) / = material_partPhase ) ! wind forward to <phase>
2013-12-12 03:33:09 +05:30
line = IO_read ( fileUnit )
2009-03-06 15:43:08 +05:30
enddo
2012-02-14 05:00:59 +05:30
2014-04-15 16:19:24 +05:30
parsingFile : do while ( trim ( line ) / = IO_EOF ) ! read through sections of phase part
2013-12-12 03:33:09 +05:30
line = IO_read ( fileUnit )
2013-06-29 00:28:10 +05:30
if ( IO_isBlank ( line ) ) cycle ! skip empty lines
2013-12-12 03:33:09 +05:30
if ( IO_getTag ( line , '<' , '>' ) / = '' ) then ! stop at next part
line = IO_read ( fileUnit , . true . ) ! reset IO_read
exit
endif
2013-06-29 00:28:10 +05:30
if ( IO_getTag ( line , '[' , ']' ) / = '' ) then ! next section
2014-04-15 16:19:24 +05:30
myConstituents = 0_pInt
2014-03-09 02:20:31 +05:30
phase = phase + 1_pInt ! advance section counter
if ( phase_plasticity ( phase ) == PLASTICITY_J2_ID ) then
instance = phase_plasticityInstance ( phase )
2014-04-15 16:19:24 +05:30
myConstituents = count ( material_phase == phase )
#ifdef HDF
2014-05-08 20:25:19 +05:30
outID ( instance ) = HDF5_addGroup ( str1 , tempResults )
2014-04-15 16:19:24 +05:30
#endif
2013-06-12 01:46:40 +05:30
endif
2014-02-10 20:01:19 +05:30
cycle ! skip to next line
2009-03-06 15:43:08 +05:30
endif
2014-04-15 16:19:24 +05:30
if ( myConstituents > 0_pInt ) then
2014-03-09 02:20:31 +05:30
instance = phase_plasticityInstance ( phase ) ! which instance of my plasticity is present phase
2014-02-10 20:01:19 +05:30
positions = IO_stringPos ( line , MAXNCHUNKS )
tag = IO_lc ( IO_stringValue ( line , positions , 1_pInt ) ) ! extract key
2014-04-15 16:19:24 +05:30
2014-02-10 20:01:19 +05:30
select case ( tag )
case ( '(output)' )
select case ( IO_lc ( IO_stringValue ( line , positions , 2_pInt ) ) )
case ( 'flowstress' )
2014-06-25 04:23:25 +05:30
constitutive_j2_Noutput ( instance ) = constitutive_j2_Noutput ( instance ) + 1_pInt
2014-06-26 19:23:12 +05:30
constitutive_j2_outputID ( constitutive_j2_Noutput ( instance ) , instance ) = flowstress_ID
2014-06-25 04:23:25 +05:30
constitutive_j2_output ( constitutive_j2_Noutput ( instance ) , instance ) = &
IO_lc ( IO_stringValue ( line , positions , 2_pInt ) )
2014-04-15 16:19:24 +05:30
#ifdef HDF
2014-05-08 20:25:19 +05:30
call HDF5_addScalarDataset ( outID ( instance ) , myConstituents , 'flowstress' , 'MPa' )
allocate ( constitutive_j2_Output2 ( instance ) % flowstress ( myConstituents ) )
2014-04-15 16:19:24 +05:30
constitutive_j2_Output2 ( instance ) % flowstressActive = . true .
#endif
2014-02-10 20:01:19 +05:30
case ( 'strainrate' )
2014-06-25 04:23:25 +05:30
constitutive_j2_Noutput ( instance ) = constitutive_j2_Noutput ( instance ) + 1_pInt
2014-06-26 19:23:12 +05:30
constitutive_j2_outputID ( constitutive_j2_Noutput ( instance ) , instance ) = strainrate_ID
2014-06-25 04:23:25 +05:30
constitutive_j2_output ( constitutive_j2_Noutput ( instance ) , instance ) = &
IO_lc ( IO_stringValue ( line , positions , 2_pInt ) )
2014-04-15 16:19:24 +05:30
#ifdef HDF
2014-05-08 20:25:19 +05:30
call HDF5_addScalarDataset ( outID ( instance ) , myConstituents , 'strainrate' , '1/s' )
allocate ( constitutive_j2_Output2 ( instance ) % strainrate ( myConstituents ) )
2014-04-15 16:19:24 +05:30
constitutive_j2_Output2 ( instance ) % strainrateActive = . true .
#endif
2014-02-10 20:01:19 +05:30
case default
2014-06-25 04:23:25 +05:30
2014-02-10 20:01:19 +05:30
end select
case ( 'tau0' )
2014-02-28 15:48:40 +05:30
constitutive_j2_tau0 ( instance ) = IO_floatValue ( line , positions , 2_pInt )
if ( constitutive_j2_tau0 ( instance ) < 0.0_pReal ) &
2014-02-10 20:01:19 +05:30
call IO_error ( 211_pInt , ext_msg = trim ( tag ) / / ' (' / / PLASTICITY_J2_label / / ')' )
case ( 'gdot0' )
2014-02-28 15:48:40 +05:30
constitutive_j2_gdot0 ( instance ) = IO_floatValue ( line , positions , 2_pInt )
if ( constitutive_j2_gdot0 ( instance ) < = 0.0_pReal ) &
2014-02-10 20:01:19 +05:30
call IO_error ( 211_pInt , ext_msg = trim ( tag ) / / ' (' / / PLASTICITY_J2_label / / ')' )
case ( 'n' )
2014-02-28 15:48:40 +05:30
constitutive_j2_n ( instance ) = IO_floatValue ( line , positions , 2_pInt )
if ( constitutive_j2_n ( instance ) < = 0.0_pReal ) &
2014-02-10 20:01:19 +05:30
call IO_error ( 211_pInt , ext_msg = trim ( tag ) / / ' (' / / PLASTICITY_J2_label / / ')' )
case ( 'h0' )
2014-02-28 15:48:40 +05:30
constitutive_j2_h0 ( instance ) = IO_floatValue ( line , positions , 2_pInt )
2014-02-10 20:01:19 +05:30
case ( 'h0_slope' , 'slopelnrate' )
2014-02-28 15:48:40 +05:30
constitutive_j2_h0_slopeLnRate ( instance ) = IO_floatValue ( line , positions , 2_pInt )
2014-02-10 20:01:19 +05:30
case ( 'tausat' )
2014-02-28 15:48:40 +05:30
constitutive_j2_tausat ( instance ) = IO_floatValue ( line , positions , 2_pInt )
if ( constitutive_j2_tausat ( instance ) < = 0.0_pReal ) &
2014-02-10 20:01:19 +05:30
call IO_error ( 211_pInt , ext_msg = trim ( tag ) / / ' (' / / PLASTICITY_J2_label / / ')' )
case ( 'tausat_sinhfita' )
2014-02-28 15:48:40 +05:30
constitutive_j2_tausat_SinhFitA ( instance ) = IO_floatValue ( line , positions , 2_pInt )
2014-02-10 20:01:19 +05:30
case ( 'tausat_sinhfitb' )
2014-02-28 15:48:40 +05:30
constitutive_j2_tausat_SinhFitB ( instance ) = IO_floatValue ( line , positions , 2_pInt )
2014-02-10 20:01:19 +05:30
case ( 'tausat_sinhfitc' )
2014-02-28 15:48:40 +05:30
constitutive_j2_tausat_SinhFitC ( instance ) = IO_floatValue ( line , positions , 2_pInt )
2014-02-10 20:01:19 +05:30
case ( 'tausat_sinhfitd' )
2014-02-28 15:48:40 +05:30
constitutive_j2_tausat_SinhFitD ( instance ) = IO_floatValue ( line , positions , 2_pInt )
2014-02-10 20:01:19 +05:30
case ( 'a' , 'w0' )
2014-02-28 15:48:40 +05:30
constitutive_j2_a ( instance ) = IO_floatValue ( line , positions , 2_pInt )
if ( constitutive_j2_a ( instance ) < = 0.0_pReal ) &
2014-02-10 20:01:19 +05:30
call IO_error ( 211_pInt , ext_msg = trim ( tag ) / / ' (' / / PLASTICITY_J2_label / / ')' )
case ( 'taylorfactor' )
2014-02-28 15:48:40 +05:30
constitutive_j2_fTaylor ( instance ) = IO_floatValue ( line , positions , 2_pInt )
if ( constitutive_j2_fTaylor ( instance ) < = 0.0_pReal ) &
2014-02-10 20:01:19 +05:30
call IO_error ( 211_pInt , ext_msg = trim ( tag ) / / ' (' / / PLASTICITY_J2_label / / ')' )
case ( 'atol_resistance' )
2014-02-28 15:48:40 +05:30
constitutive_j2_aTolResistance ( instance ) = IO_floatValue ( line , positions , 2_pInt )
if ( constitutive_j2_aTolResistance ( instance ) < = 0.0_pReal ) &
2014-02-10 20:01:19 +05:30
call IO_error ( 211_pInt , ext_msg = trim ( tag ) / / ' (' / / PLASTICITY_J2_label / / ')' )
case default
2014-06-25 04:23:25 +05:30
2014-02-10 20:01:19 +05:30
end select
2014-04-15 16:19:24 +05:30
endif
2014-03-09 02:20:31 +05:30
enddo parsingFile
2009-03-06 15:43:08 +05:30
2014-04-15 16:19:24 +05:30
initializeInstances : do phase = 1_pInt , size ( phase_plasticity )
2014-05-08 20:25:19 +05:30
NofMyPhase = count ( material_phase == phase )
if ( phase_plasticity ( phase ) == PLASTICITY_j2_ID . and . NofMyPhase / = 0 ) then
instance = phase_plasticityInstance ( phase )
2014-04-15 16:19:24 +05:30
outputsLoop : do o = 1_pInt , constitutive_j2_Noutput ( instance )
select case ( constitutive_j2_outputID ( o , instance ) )
case ( flowstress_ID , strainrate_ID )
mySize = 1_pInt
case default
end select
2010-10-26 18:46:37 +05:30
2014-04-15 16:19:24 +05:30
if ( mySize > 0_pInt ) then ! any meaningful output found
constitutive_j2_sizePostResult ( o , instance ) = mySize
constitutive_j2_sizePostResults ( instance ) = &
constitutive_j2_sizePostResults ( instance ) + mySize
endif
enddo outputsLoop
2014-05-08 20:25:19 +05:30
#ifdef NEWSTATE
2014-05-22 16:07:57 +05:30
sizeState = 1
2014-05-27 20:16:03 +05:30
plasticState ( phase ) % sizeState = sizeState
2014-06-11 22:22:18 +05:30
sizeDotState = sizeState
plasticState ( phase ) % sizeDotState = sizeDotState
2014-05-08 20:25:19 +05:30
allocate ( plasticState ( phase ) % state0 ( sizeState , NofMyPhase ) , source = constitutive_j2_tau0 ( instance ) )
allocate ( plasticState ( phase ) % partionedState0 ( sizeState , NofMyPhase ) , source = constitutive_j2_tau0 ( instance ) )
allocate ( plasticState ( phase ) % subState0 ( sizeState , NofMyPhase ) , source = 0.0_pReal )
allocate ( plasticState ( phase ) % state ( sizeState , NofMyPhase ) , source = constitutive_j2_tau0 ( instance ) )
allocate ( plasticState ( phase ) % state_backup ( sizeState , NofMyPhase ) , source = 0.0_pReal )
2014-05-22 16:07:57 +05:30
allocate ( plasticState ( phase ) % aTolState ( NofMyPhase ) , source = constitutive_j2_aTolResistance ( instance ) )
2014-05-08 20:25:19 +05:30
allocate ( plasticState ( phase ) % dotState ( sizeDotState , NofMyPhase ) , source = 0.0_pReal )
allocate ( plasticState ( phase ) % dotState_backup ( sizeDotState , NofMyPhase ) , source = 0.0_pReal )
if ( any ( numerics_integrator == 1_pInt ) ) then
allocate ( plasticState ( phase ) % previousDotState ( sizeDotState , NofMyPhase ) , source = 0.0_pReal )
allocate ( plasticState ( phase ) % previousDotState2 ( sizeDotState , NofMyPhase ) , source = 0.0_pReal )
endif
if ( any ( numerics_integrator == 4_pInt ) ) &
allocate ( plasticState ( phase ) % RK4dotState ( sizeDotState , NofMyPhase ) , source = 0.0_pReal )
if ( any ( numerics_integrator == 5_pInt ) ) &
2014-05-15 14:22:16 +05:30
allocate ( plasticState ( phase ) % RKCK45dotState ( 6 , sizeDotState , NofMyPhase ) , source = 0.0_pReal )
2014-05-08 20:25:19 +05:30
#endif
2014-04-15 16:19:24 +05:30
endif
enddo initializeInstances
2009-03-06 15:43:08 +05:30
2012-03-09 01:55:28 +05:30
end subroutine constitutive_j2_init
2009-03-06 15:43:08 +05:30
2014-05-22 16:07:57 +05:30
#ifndef NEWSTATE
2013-06-29 00:28:10 +05:30
!--------------------------------------------------------------------------------------------------
2013-07-01 11:40:42 +05:30
!> @brief sets the initial microstructural state for a given instance of this plasticity
!> @details initial microstructural state is set to the value specified by tau0
2014-05-08 20:25:19 +05:30
! not needed for new state
2013-06-29 00:28:10 +05:30
!--------------------------------------------------------------------------------------------------
2014-02-28 15:48:40 +05:30
pure function constitutive_j2_stateInit ( instance )
2009-06-23 16:09:29 +05:30
2012-03-09 01:55:28 +05:30
implicit none
2013-06-29 00:28:10 +05:30
real ( pReal ) , dimension ( 1 ) :: constitutive_j2_stateInit
2014-02-28 15:48:40 +05:30
integer ( pInt ) , intent ( in ) :: instance !< number specifying the instance of the plasticity
2013-06-29 00:28:10 +05:30
2014-02-28 15:48:40 +05:30
constitutive_j2_stateInit = constitutive_j2_tau0 ( instance )
2009-06-23 16:09:29 +05:30
2014-04-15 16:19:24 +05:30
2012-03-09 01:55:28 +05:30
end function constitutive_j2_stateInit
2009-03-06 15:43:08 +05:30
2013-06-29 00:28:10 +05:30
!--------------------------------------------------------------------------------------------------
2013-07-01 11:40:42 +05:30
!> @brief sets the relevant state values for a given instance of this plasticity
2014-05-08 20:25:19 +05:30
! not needed for new state
2013-06-29 00:28:10 +05:30
!--------------------------------------------------------------------------------------------------
2014-02-28 15:48:40 +05:30
pure function constitutive_j2_aTolState ( instance )
2009-09-18 21:07:14 +05:30
2013-07-12 12:27:15 +05:30
implicit none
2014-03-09 02:20:31 +05:30
real ( pReal ) , dimension ( 1 ) :: constitutive_j2_aTolState
integer ( pInt ) , intent ( in ) :: instance !< number specifying the instance of the plasticity
2009-09-18 21:07:14 +05:30
2014-02-28 15:48:40 +05:30
constitutive_j2_aTolState = constitutive_j2_aTolResistance ( instance )
2009-09-18 21:07:14 +05:30
2012-03-09 01:55:28 +05:30
end function constitutive_j2_aTolState
2009-09-18 21:07:14 +05:30
2014-05-22 16:07:57 +05:30
#endif
2009-09-18 21:07:14 +05:30
2013-06-29 00:28:10 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
2013-10-16 18:34:59 +05:30
pure subroutine constitutive_j2_LpAndItsTangent ( Lp , dLp_dTstar99 , Tstar_v , state , ipc , ip , el )
2013-06-29 00:28:10 +05:30
use prec , only : &
p_vec
use math , only : &
math_mul6x6 , &
math_Mandel6to33 , &
math_Plain3333to99 , &
math_deviatoric33 , &
math_mul33xx33
use mesh , only : &
mesh_NcpElems , &
mesh_maxNips
use material , only : &
homogenization_maxNgrains , &
material_phase , &
phase_plasticityInstance
2009-06-23 16:09:29 +05:30
2013-06-29 00:28:10 +05:30
implicit none
2014-03-13 12:13:49 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( out ) :: &
2013-07-01 11:40:42 +05:30
Lp !< plastic velocity gradient
2014-03-13 12:13:49 +05:30
real ( pReal ) , dimension ( 9 , 9 ) , intent ( out ) :: &
2014-02-10 20:01:19 +05:30
dLp_dTstar99 !< derivative of Lp with respect to 2nd Piola Kirchhoff stress
2013-07-01 11:40:42 +05:30
2014-03-13 12:13:49 +05:30
real ( pReal ) , dimension ( 6 ) , intent ( in ) :: &
2013-06-29 00:28:10 +05:30
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
2014-03-13 12:13:49 +05:30
integer ( pInt ) , intent ( in ) :: &
2013-06-29 00:28:10 +05:30
ipc , & !< component-ID of integration point
ip , & !< integration point
el !< element
2014-05-22 16:07:57 +05:30
#ifdef NEWSTATE
real ( pReal ) , dimension ( 1 ) , intent ( in ) :: &
state
#else
type ( p_vec ) , intent ( in ) :: &
2013-06-29 00:28:10 +05:30
state !< microstructure state
2014-05-22 16:07:57 +05:30
#endif
2013-06-29 00:28:10 +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 ) :: &
2014-02-28 15:48:40 +05:30
instance , &
2013-06-29 00:28:10 +05:30
k , l , m , n
2009-06-23 16:09:29 +05:30
2014-02-28 15:48:40 +05:30
instance = phase_plasticityInstance ( material_phase ( ipc , ip , el ) )
2013-06-29 00:28:10 +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
2013-07-12 12:27:15 +05:30
dLp_dTstar99 = 0.0_pReal
2013-06-29 00:28:10 +05:30
else
2014-05-22 16:07:57 +05:30
#ifdef NEWSTATE
gamma_dot = constitutive_j2_gdot0 ( instance ) &
* ( sqrt ( 1.5_pReal ) * norm_Tstar_dev / constitutive_j2_fTaylor ( instance ) / state ( 1 ) ) &
** constitutive_j2_n ( instance )
#else
2014-02-28 15:48:40 +05:30
gamma_dot = constitutive_j2_gdot0 ( instance ) &
2014-03-13 12:13:49 +05:30
* ( sqrt ( 1.5_pReal ) * norm_Tstar_dev / constitutive_j2_fTaylor ( instance ) / state % p ( 1 ) ) &
2014-02-28 15:48:40 +05:30
** constitutive_j2_n ( instance )
2014-05-22 16:07:57 +05:30
#endif
2014-02-28 15:48:40 +05:30
Lp = Tstar_dev_33 / norm_Tstar_dev * gamma_dot / constitutive_j2_fTaylor ( instance )
2013-06-29 00:28:10 +05:30
!--------------------------------------------------------------------------------------------------
! 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 ) &
2014-02-28 15:48:40 +05:30
dLp_dTstar_3333 ( k , l , m , n ) = ( constitutive_j2_n ( instance ) - 1.0_pReal ) * &
2013-06-29 00:28:10 +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
2014-01-16 14:48:26 +05:30
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
2014-02-28 15:48:40 +05:30
dLp_dTstar99 = math_Plain3333to99 ( gamma_dot / constitutive_j2_fTaylor ( instance ) * &
2013-06-29 00:28:10 +05:30
dLp_dTstar_3333 / norm_Tstar_dev )
end if
2012-03-09 01:55:28 +05:30
end subroutine constitutive_j2_LpAndItsTangent
2009-03-06 15:43:08 +05:30
2013-06-29 00:28:10 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculates the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
2013-10-16 18:34:59 +05:30
pure function constitutive_j2_dotState ( Tstar_v , state , ipc , ip , el )
2013-06-29 00:28:10 +05:30
use prec , only : &
p_vec
use math , only : &
math_mul6x6
use mesh , only : &
mesh_NcpElems , &
mesh_maxNips
use material , only : &
homogenization_maxNgrains , &
material_phase , &
phase_plasticityInstance
2012-05-16 20:13:26 +05:30
2013-06-29 00:28:10 +05:30
implicit none
real ( pReal ) , dimension ( 1 ) :: &
constitutive_j2_dotState
2014-05-22 16:07:57 +05:30
real ( pReal ) :: &
tempState
2014-03-13 12:13:49 +05:30
real ( pReal ) , dimension ( 6 ) , intent ( in ) :: &
2013-06-29 00:28:10 +05:30
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
2014-03-13 12:13:49 +05:30
integer ( pInt ) , intent ( in ) :: &
2013-06-29 00:28:10 +05:30
ipc , & !< component-ID of integration point
ip , & !< integration point
el !< element
2014-05-22 16:07:57 +05:30
#ifdef NEWSTATE
real ( pReal ) , dimension ( 1 ) , intent ( in ) :: &
state
#else
2014-03-13 12:13:49 +05:30
type ( p_vec ) , intent ( in ) :: &
2013-06-29 00:28:10 +05:30
state !< microstructure state
2014-05-22 16:07:57 +05:30
#endif
2013-06-29 00:28:10 +05:30
real ( pReal ) , dimension ( 6 ) :: &
Tstar_dev_v !< deviatoric part of the 2nd Piola Kirchhoff stress tensor in Mandel notation
real ( pReal ) :: &
gamma_dot , & !< strainrate
hardening , & !< hardening coefficient
saturation , & !< saturation resistance
norm_Tstar_dev !< euclidean norm of Tstar_dev
integer ( pInt ) :: &
2014-02-28 15:48:40 +05:30
instance
2014-05-22 16:07:57 +05:30
#ifdef NEWSTATE
tempState = state ( 1 )
#else
tempState = state % p ( 1 )
#endif
2012-05-16 20:13:26 +05:30
2014-02-28 15:48:40 +05:30
instance = phase_plasticityInstance ( material_phase ( ipc , ip , el ) )
2013-06-29 00:28:10 +05:30
!--------------------------------------------------------------------------------------------------
! norm of deviatoric part of 2nd Piola-Kirchhoff stress
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_dev = sqrt ( math_mul6x6 ( Tstar_dev_v , Tstar_dev_v ) )
!--------------------------------------------------------------------------------------------------
! strain rate
2014-02-28 15:48:40 +05:30
gamma_dot = constitutive_j2_gdot0 ( instance ) * ( sqrt ( 1.5_pReal ) * norm_Tstar_dev &
2013-06-29 00:28:10 +05:30
/ & !-----------------------------------------------------------------------------------
2014-05-22 16:07:57 +05:30
( constitutive_j2_fTaylor ( instance ) * tempState ) ) ** constitutive_j2_n ( instance )
2013-06-29 00:28:10 +05:30
!--------------------------------------------------------------------------------------------------
! hardening coefficient
if ( abs ( gamma_dot ) > 1e-12_pReal ) then
2014-02-28 15:48:40 +05:30
if ( constitutive_j2_tausat_SinhFitA ( instance ) == 0.0_pReal ) then
saturation = constitutive_j2_tausat ( instance )
2013-06-29 00:28:10 +05:30
else
2014-02-28 15:48:40 +05:30
saturation = ( constitutive_j2_tausat ( instance ) &
+ ( log ( ( gamma_dot / constitutive_j2_tausat_SinhFitA ( instance ) &
) ** ( 1.0_pReal / constitutive_j2_tausat_SinhFitD ( instance ) ) &
+ sqrt ( ( gamma_dot / constitutive_j2_tausat_SinhFitA ( instance ) &
) ** ( 2.0_pReal / constitutive_j2_tausat_SinhFitD ( instance ) ) &
2013-06-29 00:28:10 +05:30
+ 1.0_pReal ) &
) & ! asinh(K) = ln(K + sqrt(K^2 +1))
2014-02-28 15:48:40 +05:30
) ** ( 1.0_pReal / constitutive_j2_tausat_SinhFitC ( instance ) ) &
/ ( constitutive_j2_tausat_SinhFitB ( instance ) &
* ( gamma_dot / constitutive_j2_gdot0 ( instance ) ) ** ( 1.0_pReal / constitutive_j2_n ( instance ) ) &
2013-06-29 00:28:10 +05:30
) &
)
endif
2014-02-28 15:48:40 +05:30
hardening = ( constitutive_j2_h0 ( instance ) + constitutive_j2_h0_slopeLnRate ( instance ) * log ( gamma_dot ) ) &
2014-05-22 16:07:57 +05:30
* abs ( 1.0_pReal - tempState / saturation ) ** constitutive_j2_a ( instance ) &
* sign ( 1.0_pReal , 1.0_pReal - tempState / saturation )
2013-06-29 00:28:10 +05:30
else
hardening = 0.0_pReal
endif
constitutive_j2_dotState = hardening * gamma_dot
2012-05-16 20:13:26 +05:30
2013-06-29 00:28:10 +05:30
end function constitutive_j2_dotState
2012-05-16 20:13:26 +05:30
2013-06-29 00:28:10 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief return array of constitutive results
!--------------------------------------------------------------------------------------------------
2013-10-16 18:34:59 +05:30
pure function constitutive_j2_postResults ( Tstar_v , state , ipc , ip , el )
2013-06-29 00:28:10 +05:30
use prec , only : &
p_vec
use math , only : &
math_mul6x6
use mesh , only : &
mesh_NcpElems , &
mesh_maxNips
use material , only : &
homogenization_maxNgrains , &
material_phase , &
phase_plasticityInstance , &
phase_Noutput
implicit none
2014-03-13 12:13:49 +05:30
real ( pReal ) , dimension ( 6 ) , intent ( in ) :: &
2013-06-29 00:28:10 +05:30
Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation
2014-03-13 12:13:49 +05:30
integer ( pInt ) , intent ( in ) :: &
2013-06-29 00:28:10 +05:30
ipc , & !< component-ID of integration point
ip , & !< integration point
el !< element
2014-05-22 16:07:57 +05:30
real ( pReal ) :: &
tempState
#ifdef NEWSTATE
real ( pReal ) , dimension ( 1 ) , intent ( in ) :: &
state
#else
type ( p_vec ) , intent ( in ) :: &
2013-07-01 11:40:42 +05:30
state !< microstructure state
2014-05-22 16:07:57 +05:30
#endif
2013-06-29 00:28:10 +05:30
real ( pReal ) , dimension ( constitutive_j2_sizePostResults ( phase_plasticityInstance ( material_phase ( ipc , ip , el ) ) ) ) :: &
constitutive_j2_postResults
2014-05-22 16:07:57 +05:30
2013-06-29 00:28:10 +05:30
real ( pReal ) , dimension ( 6 ) :: &
Tstar_dev_v ! deviatoric part of the 2nd Piola Kirchhoff stress tensor in Mandel notation
real ( pReal ) :: &
norm_Tstar_dev ! euclidean norm of Tstar_dev
integer ( pInt ) :: &
2014-02-28 15:48:40 +05:30
instance , &
2013-06-29 00:28:10 +05:30
o , &
c
2014-05-22 16:07:57 +05:30
#ifdef NEWSTATE
tempState = state ( 1 )
#else
tempState = state % p ( 1 )
#endif
2013-06-29 00:28:10 +05:30
2014-02-28 15:48:40 +05:30
instance = phase_plasticityInstance ( material_phase ( ipc , ip , el ) )
2013-06-29 00:28:10 +05:30
!--------------------------------------------------------------------------------------------------
! calculate deviatoric part of 2nd Piola-Kirchhoff stress and its norm
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_dev = sqrt ( math_mul6x6 ( Tstar_dev_v , Tstar_dev_v ) )
c = 0_pInt
constitutive_j2_postResults = 0.0_pReal
2014-06-25 04:23:25 +05:30
outputsLoop : do o = 1_pInt , constitutive_j2_Noutput ( instance )
2014-02-28 15:48:40 +05:30
select case ( constitutive_j2_outputID ( o , instance ) )
2013-12-12 03:33:09 +05:30
case ( flowstress_ID )
2014-05-22 16:07:57 +05:30
constitutive_j2_postResults ( c + 1_pInt ) = tempState
2013-06-29 00:28:10 +05:30
c = c + 1_pInt
2013-12-12 03:33:09 +05:30
case ( strainrate_ID )
2013-06-29 00:28:10 +05:30
constitutive_j2_postResults ( c + 1_pInt ) = &
2014-02-28 15:48:40 +05:30
constitutive_j2_gdot0 ( instance ) * ( sqrt ( 1.5_pReal ) * norm_Tstar_dev &
2013-06-29 00:28:10 +05:30
/ & !----------------------------------------------------------------------------------
2014-05-22 16:07:57 +05:30
( constitutive_j2_fTaylor ( instance ) * tempState ) ) ** constitutive_j2_n ( instance )
2013-06-29 00:28:10 +05:30
c = c + 1_pInt
end select
enddo outputsLoop
2009-03-06 15:43:08 +05:30
2012-03-09 01:55:28 +05:30
end function constitutive_j2_postResults
2009-03-06 15:43:08 +05:30
2014-04-15 16:19:24 +05:30
2012-03-09 01:55:28 +05:30
end module constitutive_j2