2009-10-12 21:56:16 +05:30
!* $Id$
2009-10-06 20:49:19 +05:30
!************************************
!* Module: CONSTITUTIVE *
!************************************
MODULE constitutive_dislotwin
2009-10-10 01:58:23 +05:30
!* Include other modules
use prec , only : pReal , pInt
implicit none
2009-10-06 20:49:19 +05:30
2009-10-10 01:58:23 +05:30
!* Lists of states and physical parameters
character ( len = * ) , parameter :: constitutive_dislotwin_label = 'dislotwin'
character ( len = 18 ) , dimension ( 2 ) , parameter :: constitutive_dislotwin_listBasicSlipStates = ( / 'rhoEdge ' , &
2009-10-12 21:56:16 +05:30
'rhoEdgeDip' / )
2010-08-04 05:17:00 +05:30
character ( len = 18 ) , dimension ( 1 ) , parameter :: constitutive_dislotwin_listBasicTwinStates = ( / 'twinFraction' / )
2009-10-10 01:58:23 +05:30
character ( len = 18 ) , dimension ( 4 ) , parameter :: constitutive_dislotwin_listDependentSlipStates = ( / 'invLambdaSlip ' , &
2009-10-12 21:56:16 +05:30
'invLambdaSlipTwin' , &
'meanFreePathSlip ' , &
'tauSlipThreshold ' / )
2009-10-10 01:58:23 +05:30
character ( len = 18 ) , dimension ( 4 ) , parameter :: constitutive_dislotwin_listDependentTwinStates = ( / 'invLambdaTwin ' , &
2010-08-04 05:17:00 +05:30
'meanFreePathTwin' , &
2009-10-12 21:56:16 +05:30
'tauTwinThreshold' , &
'twinVolume ' / )
2009-10-10 01:58:23 +05:30
real ( pReal ) , parameter :: kB = 1.38e-23_pReal ! Boltzmann constant in J/Kelvin
!* Definition of global variables
integer ( pInt ) , dimension ( : ) , allocatable :: constitutive_dislotwin_sizeDotState , & ! number of dotStates
constitutive_dislotwin_sizeState , & ! total number of microstructural state variables
constitutive_dislotwin_sizePostResults ! cumulative size of post results
integer ( pInt ) , dimension ( : , : ) , allocatable , target :: constitutive_dislotwin_sizePostResult ! size of each post result output
2010-08-17 19:53:55 +05:30
character ( len = 64 ) , dimension ( : , : ) , allocatable , target :: constitutive_dislotwin_output ! name of each post result output
2009-10-10 01:58:23 +05:30
character ( len = 32 ) , dimension ( : ) , allocatable :: constitutive_dislotwin_structureName ! name of the lattice structure
integer ( pInt ) , dimension ( : ) , allocatable :: constitutive_dislotwin_structure , & ! number representing the kind of lattice structure
constitutive_dislotwin_totalNslip , & ! total number of active slip systems for each instance
constitutive_dislotwin_totalNtwin ! total number of active twin systems for each instance
integer ( pInt ) , dimension ( : , : ) , allocatable :: constitutive_dislotwin_Nslip , & ! number of active slip systems for each family and instance
constitutive_dislotwin_Ntwin , & ! number of active twin systems for each family and instance
constitutive_dislotwin_slipFamily , & ! lookup table relating active slip system to slip family for each instance
constitutive_dislotwin_twinFamily , & ! lookup table relating active twin system to twin family for each instance
constitutive_dislotwin_slipSystemLattice , & ! lookup table relating active slip system index to lattice slip system index for each instance
constitutive_dislotwin_twinSystemLattice ! lookup table relating active twin system index to lattice twin system index for each instance
real ( pReal ) , dimension ( : ) , allocatable :: constitutive_dislotwin_CoverA , & ! c/a ratio for hex type lattice
constitutive_dislotwin_C11 , & ! C11 element in elasticity matrix
constitutive_dislotwin_C12 , & ! C12 element in elasticity matrix
constitutive_dislotwin_C13 , & ! C13 element in elasticity matrix
constitutive_dislotwin_C33 , & ! C33 element in elasticity matrix
constitutive_dislotwin_C44 , & ! C44 element in elasticity matrix
constitutive_dislotwin_Gmod , & ! shear modulus
constitutive_dislotwin_CAtomicVolume , & ! atomic volume in Bugers vector unit
constitutive_dislotwin_D0 , & ! prefactor for self-diffusion coefficient
constitutive_dislotwin_Qsd , & ! activation energy for dislocation climb
constitutive_dislotwin_GrainSize , & ! grain size
constitutive_dislotwin_p , & ! p-exponent in glide velocity
constitutive_dislotwin_q , & ! q-exponent in glide velocity
constitutive_dislotwin_MaxTwinFraction , & ! maximum allowed total twin volume fraction
constitutive_dislotwin_r , & ! r-exponent in twin nucleation rate
2010-08-17 19:53:55 +05:30
constitutive_dislotwin_CEdgeDipMinDistance , & !
constitutive_dislotwin_Cmfptwin , & !
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_Cthresholdtwin , & !
2010-08-17 19:53:55 +05:30
constitutive_dislotwin_SolidSolutionStrength , & ! Strength due to elements in solid solution
constitutive_dislotwin_L0 , & ! Length of twin nuclei in Burgers vectors
2010-10-26 18:46:37 +05:30
constitutive_dislotwin_aTolRho ! absolute tolerance for integration of dislocation density
2009-10-10 01:58:23 +05:30
real ( pReal ) , dimension ( : , : , : ) , allocatable :: constitutive_dislotwin_Cslip_66 ! elasticity matrix in Mandel notation for each instance
real ( pReal ) , dimension ( : , : , : , : ) , allocatable :: constitutive_dislotwin_Ctwin_66 ! twin elasticity matrix in Mandel notation for each instance
real ( pReal ) , dimension ( : , : , : , : , : ) , allocatable :: constitutive_dislotwin_Cslip_3333 ! elasticity matrix for each instance
real ( pReal ) , dimension ( : , : , : , : , : , : ) , allocatable :: constitutive_dislotwin_Ctwin_3333 ! twin elasticity matrix for each instance
real ( pReal ) , dimension ( : , : ) , allocatable :: constitutive_dislotwin_rhoEdge0 , & ! initial edge dislocation density per slip system for each family and instance
constitutive_dislotwin_rhoEdgeDip0 , & ! initial edge dipole density per slip system for each family and instance
constitutive_dislotwin_burgersPerSlipFamily , & ! absolute length of burgers vector [m] for each slip family and instance
constitutive_dislotwin_burgersPerSlipSystem , & ! absolute length of burgers vector [m] for each slip system and instance
constitutive_dislotwin_burgersPerTwinFamily , & ! absolute length of burgers vector [m] for each twin family and instance
constitutive_dislotwin_burgersPerTwinSystem , & ! absolute length of burgers vector [m] for each twin system and instance
constitutive_dislotwin_QedgePerSlipFamily , & ! activation energy for glide [J] for each slip family and instance
constitutive_dislotwin_QedgePerSlipSystem , & ! activation energy for glide [J] for each slip system and instance
constitutive_dislotwin_v0PerSlipFamily , & ! dislocation velocity prefactor [m/s] for each family and instance
constitutive_dislotwin_v0PerSlipSystem , & ! dislocation velocity prefactor [m/s] for each slip system and instance
constitutive_dislotwin_Ndot0PerTwinFamily , & ! twin nucleation rate [1/m³s] for each twin family and instance
constitutive_dislotwin_Ndot0PerTwinSystem , & ! twin nucleation rate [1/m³s] for each twin system and instance
constitutive_dislotwin_twinsizePerTwinFamily , & ! twin thickness [m] for each twin family and instance
constitutive_dislotwin_twinsizePerTwinSystem , & ! twin thickness [m] for each twin system and instance
constitutive_dislotwin_CLambdaSlipPerSlipFamily , & ! Adj. parameter for distance between 2 forest dislocations for each slip family and instance
2009-10-12 21:56:16 +05:30
constitutive_dislotwin_CLambdaSlipPerSlipSystem , & ! Adj. parameter for distance between 2 forest dislocations for each slip system and instance
constitutive_dislotwin_interactionSlipSlip , & ! coefficients for slip-slip interaction for each interaction type and instance
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_interactionSlipTwin , & ! coefficients for slip-twin interaction for each interaction type and instance
constitutive_dislotwin_interactionTwinSlip , & ! coefficients for twin-slip interaction for each interaction type and instance
constitutive_dislotwin_interactionTwinTwin ! coefficients for twin-twin interaction for each interaction type and instance
real ( pReal ) , dimension ( : , : , : ) , allocatable :: constitutive_dislotwin_interactionMatrixSlipSlip , & ! interaction matrix of the different slip systems for each instance
constitutive_dislotwin_interactionMatrixSlipTwin , & ! interaction matrix of slip systems with twin systems for each instance
constitutive_dislotwin_interactionMatrixTwinSlip , & ! interaction matrix of twin systems with slip systems for each instance
2010-08-17 19:53:55 +05:30
constitutive_dislotwin_interactionMatrixTwinTwin , & ! interaction matrix of the different twin systems for each instance
constitutive_dislotwin_forestProjectionEdge ! matrix of forest projections of edge dislocations for each instance
2009-10-06 20:49:19 +05:30
CONTAINS
!****************************************
2009-10-10 01:58:23 +05:30
!* - constitutive_dislotwin_init
!* - constitutive_dislotwin_stateInit
!* - constitutive_dislotwin_relevantState
!* - constitutive_dislotwin_homogenizedC
!* - constitutive_dislotwin_microstructure
!* - constitutive_dislotwin_LpAndItsTangent
!* - consistutive_dislotwin_dotState
!* - constitutive_dislotwin_dotTemperature
!* - consistutive_dislotwin_postResults
2009-10-06 20:49:19 +05:30
!****************************************
subroutine constitutive_dislotwin_init ( file )
!**************************************
!* Module initialization *
!**************************************
2009-10-10 01:58:23 +05:30
use prec , only : pInt , pReal
use math , only : math_Mandel3333to66 , math_Voigt66to3333 , math_mul3x3
use IO
use material
use lattice
!* Input variables
integer ( pInt ) , intent ( in ) :: file
!* Local variables
integer ( pInt ) , parameter :: maxNchunks = 21
integer ( pInt ) , dimension ( 1 + 2 * maxNchunks ) :: positions
integer ( pInt ) section , maxNinstance , f , i , j , k , l , m , n , o , p , q , r , s , s1 , s2 , t1 , t2 , ns , nt , output , mySize , myStructure , maxTotalNslip , maxTotalNtwin
character ( len = 64 ) tag
character ( len = 1024 ) line
!write(6,*)
!write(6,'(a20,a20,a12)') '<<<+- constitutive_',constitutive_dislotwin_label,' init -+>>>'
2009-10-12 21:56:16 +05:30
!write(6,*) '$Id$'
2009-10-10 01:58:23 +05:30
!write(6,*)
maxNinstance = count ( phase_constitution == constitutive_dislotwin_label )
if ( maxNinstance == 0 ) return
2009-10-12 21:56:16 +05:30
2009-10-10 01:58:23 +05:30
!* Space allocation for global variables
2010-08-17 19:53:55 +05:30
allocate ( constitutive_dislotwin_sizeDotState ( maxNinstance ) )
allocate ( constitutive_dislotwin_sizeState ( maxNinstance ) )
allocate ( constitutive_dislotwin_sizePostResults ( maxNinstance ) )
2009-10-10 01:58:23 +05:30
allocate ( constitutive_dislotwin_sizePostResult ( maxval ( phase_Noutput ) , maxNinstance ) )
2009-10-12 21:56:16 +05:30
allocate ( constitutive_dislotwin_output ( maxval ( phase_Noutput ) , maxNinstance ) )
constitutive_dislotwin_sizeDotState = 0_pInt
constitutive_dislotwin_sizeState = 0_pInt
constitutive_dislotwin_sizePostResults = 0_pInt
constitutive_dislotwin_sizePostResult = 0_pInt
constitutive_dislotwin_output = ''
2009-10-10 01:58:23 +05:30
2010-08-17 19:53:55 +05:30
allocate ( constitutive_dislotwin_structureName ( maxNinstance ) )
2009-10-10 01:58:23 +05:30
allocate ( constitutive_dislotwin_structure ( maxNinstance ) )
2009-10-12 21:56:16 +05:30
allocate ( constitutive_dislotwin_Nslip ( lattice_maxNslipFamily , maxNinstance ) )
2010-08-17 19:53:55 +05:30
allocate ( constitutive_dislotwin_Ntwin ( lattice_maxNtwinFamily , maxNinstance ) )
2009-10-12 21:56:16 +05:30
allocate ( constitutive_dislotwin_slipFamily ( lattice_maxNslip , maxNinstance ) )
2009-10-10 01:58:23 +05:30
allocate ( constitutive_dislotwin_twinFamily ( lattice_maxNtwin , maxNinstance ) )
allocate ( constitutive_dislotwin_slipSystemLattice ( lattice_maxNslip , maxNinstance ) )
2010-08-17 19:53:55 +05:30
allocate ( constitutive_dislotwin_twinSystemLattice ( lattice_maxNtwin , maxNinstance ) )
allocate ( constitutive_dislotwin_totalNslip ( maxNinstance ) )
allocate ( constitutive_dislotwin_totalNtwin ( maxNinstance ) )
2009-10-12 21:56:16 +05:30
constitutive_dislotwin_structureName = ''
constitutive_dislotwin_structure = 0_pInt
constitutive_dislotwin_Nslip = 0_pInt
constitutive_dislotwin_Ntwin = 0_pInt
constitutive_dislotwin_slipFamily = 0_pInt
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_twinFamily = 0_pInt
constitutive_dislotwin_slipSystemLattice = 0.0_pReal
2009-10-12 21:56:16 +05:30
constitutive_dislotwin_twinSystemLattice = 0.0_pReal
constitutive_dislotwin_totalNslip = 0_pInt
constitutive_dislotwin_totalNtwin = 0_pInt
allocate ( constitutive_dislotwin_CoverA ( maxNinstance ) )
allocate ( constitutive_dislotwin_C11 ( maxNinstance ) )
allocate ( constitutive_dislotwin_C12 ( maxNinstance ) )
allocate ( constitutive_dislotwin_C13 ( maxNinstance ) )
allocate ( constitutive_dislotwin_C33 ( maxNinstance ) )
allocate ( constitutive_dislotwin_C44 ( maxNinstance ) )
allocate ( constitutive_dislotwin_Gmod ( maxNinstance ) )
allocate ( constitutive_dislotwin_CAtomicVolume ( maxNinstance ) )
2009-10-10 01:58:23 +05:30
allocate ( constitutive_dislotwin_D0 ( maxNinstance ) )
allocate ( constitutive_dislotwin_Qsd ( maxNinstance ) )
allocate ( constitutive_dislotwin_GrainSize ( maxNinstance ) )
allocate ( constitutive_dislotwin_p ( maxNinstance ) )
allocate ( constitutive_dislotwin_q ( maxNinstance ) )
allocate ( constitutive_dislotwin_MaxTwinFraction ( maxNinstance ) )
allocate ( constitutive_dislotwin_r ( maxNinstance ) )
2009-10-12 21:56:16 +05:30
allocate ( constitutive_dislotwin_CEdgeDipMinDistance ( maxNinstance ) )
allocate ( constitutive_dislotwin_Cmfptwin ( maxNinstance ) )
2009-10-10 01:58:23 +05:30
allocate ( constitutive_dislotwin_Cthresholdtwin ( maxNinstance ) )
2010-08-17 19:53:55 +05:30
allocate ( constitutive_dislotwin_SolidSolutionStrength ( maxNinstance ) )
allocate ( constitutive_dislotwin_L0 ( maxNinstance ) )
2010-10-26 18:46:37 +05:30
allocate ( constitutive_dislotwin_aTolRho ( maxNinstance ) )
2009-10-12 21:56:16 +05:30
allocate ( constitutive_dislotwin_Cslip_66 ( 6 , 6 , maxNinstance ) )
allocate ( constitutive_dislotwin_Cslip_3333 ( 3 , 3 , 3 , 3 , maxNinstance ) )
2010-10-26 18:46:37 +05:30
constitutive_dislotwin_CoverA = 0.0_pReal
constitutive_dislotwin_C11 = 0.0_pReal
constitutive_dislotwin_C12 = 0.0_pReal
constitutive_dislotwin_C13 = 0.0_pReal
constitutive_dislotwin_C33 = 0.0_pReal
constitutive_dislotwin_C44 = 0.0_pReal
constitutive_dislotwin_Gmod = 0.0_pReal
constitutive_dislotwin_CAtomicVolume = 0.0_pReal
constitutive_dislotwin_D0 = 0.0_pReal
constitutive_dislotwin_Qsd = 0.0_pReal
constitutive_dislotwin_GrainSize = 0.0_pReal
constitutive_dislotwin_p = 0.0_pReal
constitutive_dislotwin_q = 0.0_pReal
constitutive_dislotwin_MaxTwinFraction = 0.0_pReal
constitutive_dislotwin_r = 0.0_pReal
constitutive_dislotwin_CEdgeDipMinDistance = 0.0_pReal
constitutive_dislotwin_Cmfptwin = 0.0_pReal
constitutive_dislotwin_Cthresholdtwin = 0.0_pReal
constitutive_dislotwin_SolidSolutionStrength = 0.0_pReal
constitutive_dislotwin_L0 = 0.0_pReal
constitutive_dislotwin_aTolRho = 0.0_pReal
constitutive_dislotwin_Cslip_66 = 0.0_pReal
constitutive_dislotwin_Cslip_3333 = 0.0_pReal
2009-10-10 01:58:23 +05:30
allocate ( constitutive_dislotwin_rhoEdge0 ( lattice_maxNslipFamily , maxNinstance ) )
2010-08-17 19:53:55 +05:30
allocate ( constitutive_dislotwin_rhoEdgeDip0 ( lattice_maxNslipFamily , maxNinstance ) )
2009-10-10 01:58:23 +05:30
allocate ( constitutive_dislotwin_burgersPerSlipFamily ( lattice_maxNslipFamily , maxNinstance ) )
allocate ( constitutive_dislotwin_burgersPerTwinFamily ( lattice_maxNtwinFamily , maxNinstance ) )
allocate ( constitutive_dislotwin_QedgePerSlipFamily ( lattice_maxNslipFamily , maxNinstance ) )
allocate ( constitutive_dislotwin_v0PerSlipFamily ( lattice_maxNslipFamily , maxNinstance ) )
allocate ( constitutive_dislotwin_Ndot0PerTwinFamily ( lattice_maxNtwinFamily , maxNinstance ) )
allocate ( constitutive_dislotwin_twinsizePerTwinFamily ( lattice_maxNtwinFamily , maxNinstance ) )
allocate ( constitutive_dislotwin_CLambdaSlipPerSlipFamily ( lattice_maxNslipFamily , maxNinstance ) )
2009-10-12 21:56:16 +05:30
constitutive_dislotwin_rhoEdge0 = 0.0_pReal
constitutive_dislotwin_rhoEdgeDip0 = 0.0_pReal
constitutive_dislotwin_burgersPerSlipFamily = 0.0_pReal
constitutive_dislotwin_burgersPerTwinFamily = 0.0_pReal
constitutive_dislotwin_QedgePerSlipFamily = 0.0_pReal
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_v0PerSlipFamily = 0.0_pReal
constitutive_dislotwin_Ndot0PerTwinFamily = 0.0_pReal
2009-10-12 21:56:16 +05:30
constitutive_dislotwin_twinsizePerTwinFamily = 0.0_pReal
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_CLambdaSlipPerSlipFamily = 0.0_pReal
2010-08-17 19:53:55 +05:30
allocate ( constitutive_dislotwin_interactionSlipSlip ( lattice_maxNinteraction , maxNinstance ) )
allocate ( constitutive_dislotwin_interactionSlipTwin ( lattice_maxNinteraction , maxNinstance ) )
allocate ( constitutive_dislotwin_interactionTwinSlip ( lattice_maxNinteraction , maxNinstance ) )
allocate ( constitutive_dislotwin_interactionTwinTwin ( lattice_maxNinteraction , maxNinstance ) )
2010-10-26 18:46:37 +05:30
constitutive_dislotwin_interactionSlipSlip = 0.0_pReal
constitutive_dislotwin_interactionSlipTwin = 0.0_pReal
constitutive_dislotwin_interactionTwinSlip = 0.0_pReal
constitutive_dislotwin_interactionTwinTwin = 0.0_pReal
2009-10-12 21:56:16 +05:30
2009-10-10 01:58:23 +05:30
!* Readout data from material.config file
rewind ( file )
line = ''
section = 0
2010-08-17 19:53:55 +05:30
2009-10-10 01:58:23 +05:30
do while ( IO_lc ( IO_getTag ( line , '<' , '>' ) ) / = 'phase' ) ! wind forward to <phase>
2009-10-06 20:49:19 +05:30
read ( file , '(a1024)' , END = 100 ) line
2009-10-10 01:58:23 +05:30
enddo
2009-10-06 20:49:19 +05:30
2009-10-10 01:58:23 +05:30
do ! read thru sections of phase part
2009-10-06 20:49:19 +05:30
read ( file , '(a1024)' , END = 100 ) line
if ( IO_isBlank ( line ) ) cycle ! skip empty lines
if ( IO_getTag ( line , '<' , '>' ) / = '' ) exit ! stop at next part
if ( IO_getTag ( line , '[' , ']' ) / = '' ) then ! next section
section = section + 1
output = 0 ! reset output counter
endif
if ( section > 0 . and . phase_constitution ( section ) == constitutive_dislotwin_label ) then ! one of my sections
i = phase_constitutionInstance ( section ) ! which instance of my constitution is present phase
positions = IO_stringPos ( line , maxNchunks )
tag = IO_lc ( IO_stringValue ( line , positions , 1 ) ) ! extract key
select case ( tag )
case ( '(output)' )
output = output + 1
constitutive_dislotwin_output ( output , i ) = IO_lc ( IO_stringValue ( line , positions , 2 ) )
case ( 'lattice_structure' )
constitutive_dislotwin_structureName ( i ) = IO_lc ( IO_stringValue ( line , positions , 2 ) )
case ( 'covera_ratio' )
constitutive_dislotwin_CoverA ( i ) = IO_floatValue ( line , positions , 2 )
case ( 'c11' )
constitutive_dislotwin_C11 ( i ) = IO_floatValue ( line , positions , 2 )
case ( 'c12' )
constitutive_dislotwin_C12 ( i ) = IO_floatValue ( line , positions , 2 )
case ( 'c13' )
constitutive_dislotwin_C13 ( i ) = IO_floatValue ( line , positions , 2 )
case ( 'c33' )
constitutive_dislotwin_C33 ( i ) = IO_floatValue ( line , positions , 2 )
case ( 'c44' )
constitutive_dislotwin_C44 ( i ) = IO_floatValue ( line , positions , 2 )
case ( 'nslip' )
2009-10-12 21:56:16 +05:30
forall ( j = 1 : lattice_maxNslipFamily ) &
constitutive_dislotwin_Nslip ( j , i ) = IO_intValue ( line , positions , 1 + j )
2009-10-06 20:49:19 +05:30
case ( 'ntwin' )
2009-10-12 21:56:16 +05:30
forall ( j = 1 : lattice_maxNtwinFamily ) &
constitutive_dislotwin_Ntwin ( j , i ) = IO_intValue ( line , positions , 1 + j )
case ( 'rhoedge0' )
forall ( j = 1 : lattice_maxNslipFamily ) &
constitutive_dislotwin_rhoEdge0 ( j , i ) = IO_floatValue ( line , positions , 1 + j )
case ( 'rhoedgedip0' )
forall ( j = 1 : lattice_maxNslipFamily ) &
constitutive_dislotwin_rhoEdgeDip0 ( j , i ) = IO_floatValue ( line , positions , 1 + j )
2009-10-10 01:58:23 +05:30
case ( 'slipburgers' )
2009-10-12 21:56:16 +05:30
forall ( j = 1 : lattice_maxNslipFamily ) &
constitutive_dislotwin_burgersPerSlipFamily ( j , i ) = IO_floatValue ( line , positions , 1 + j )
case ( 'twinburgers' )
forall ( j = 1 : lattice_maxNtwinFamily ) &
constitutive_dislotwin_burgersPerTwinFamily ( j , i ) = IO_floatValue ( line , positions , 1 + j )
2009-10-06 20:49:19 +05:30
case ( 'qedge' )
2009-10-12 21:56:16 +05:30
forall ( j = 1 : lattice_maxNslipFamily ) &
constitutive_dislotwin_QedgePerSlipFamily ( j , i ) = IO_floatValue ( line , positions , 1 + j )
case ( 'v0' )
forall ( j = 1 : lattice_maxNslipFamily ) &
constitutive_dislotwin_v0PerSlipFamily ( j , i ) = IO_floatValue ( line , positions , 1 + j )
case ( 'ndot0' )
forall ( j = 1 : lattice_maxNtwinFamily ) &
constitutive_dislotwin_Ndot0PerTwinFamily ( j , i ) = IO_floatValue ( line , positions , 1 + j )
case ( 'twinsize' )
forall ( j = 1 : lattice_maxNtwinFamily ) &
constitutive_dislotwin_twinsizePerTwinFamily ( j , i ) = IO_floatValue ( line , positions , 1 + j )
2009-10-10 01:58:23 +05:30
case ( 'clambdaslip' )
2009-10-12 21:56:16 +05:30
forall ( j = 1 : lattice_maxNslipFamily ) &
constitutive_dislotwin_CLambdaSlipPerSlipFamily ( j , i ) = IO_floatValue ( line , positions , 1 + j )
case ( 'grainsize' )
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_GrainSize ( i ) = IO_floatValue ( line , positions , 2 )
case ( 'maxtwinfraction' )
constitutive_dislotwin_MaxTwinFraction ( i ) = IO_floatValue ( line , positions , 2 )
2009-10-12 21:56:16 +05:30
case ( 'pexponent' )
constitutive_dislotwin_p ( i ) = IO_floatValue ( line , positions , 2 )
case ( 'qexponent' )
constitutive_dislotwin_q ( i ) = IO_floatValue ( line , positions , 2 )
case ( 'rexponent' )
constitutive_dislotwin_r ( i ) = IO_floatValue ( line , positions , 2 )
case ( 'd0' )
constitutive_dislotwin_D0 ( i ) = IO_floatValue ( line , positions , 2 )
case ( 'qsd' )
constitutive_dislotwin_Qsd ( i ) = IO_floatValue ( line , positions , 2 )
2010-10-26 18:46:37 +05:30
case ( 'atol_rho' )
constitutive_dislotwin_aTolRho ( i ) = IO_floatValue ( line , positions , 2 )
2010-08-17 19:53:55 +05:30
case ( 'cmfptwin' )
2009-10-06 20:49:19 +05:30
constitutive_dislotwin_Cmfptwin ( i ) = IO_floatValue ( line , positions , 2 )
2010-08-17 19:53:55 +05:30
case ( 'cthresholdtwin' )
2009-10-06 20:49:19 +05:30
constitutive_dislotwin_Cthresholdtwin ( i ) = IO_floatValue ( line , positions , 2 )
2010-10-15 18:29:00 +05:30
case ( 'solidsolutionstrength' )
2010-08-17 19:53:55 +05:30
constitutive_dislotwin_SolidSolutionStrength ( i ) = IO_floatValue ( line , positions , 2 )
2010-10-26 18:46:37 +05:30
case ( 'l0' )
2010-08-17 19:53:55 +05:30
constitutive_dislotwin_L0 ( i ) = IO_floatValue ( line , positions , 2 )
case ( 'cedgedipmindistance' )
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_CEdgeDipMinDistance ( i ) = IO_floatValue ( line , positions , 2 )
2010-08-17 19:53:55 +05:30
case ( 'catomicvolume' )
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_CAtomicVolume ( i ) = IO_floatValue ( line , positions , 2 )
case ( 'interactionslipslip' )
2009-10-06 20:49:19 +05:30
forall ( j = 1 : lattice_maxNinteraction ) &
2009-10-12 21:56:16 +05:30
constitutive_dislotwin_interactionSlipSlip ( j , i ) = IO_floatValue ( line , positions , 1 + j )
2009-10-10 01:58:23 +05:30
case ( 'interactionsliptwin' )
2009-10-06 20:49:19 +05:30
forall ( j = 1 : lattice_maxNinteraction ) &
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_interactionSlipTwin ( j , i ) = IO_floatValue ( line , positions , 1 + j )
case ( 'interactiontwinslip' )
2009-10-06 20:49:19 +05:30
forall ( j = 1 : lattice_maxNinteraction ) &
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_interactionTwinSlip ( j , i ) = IO_floatValue ( line , positions , 1 + j )
case ( 'interactiontwintwin' )
2009-10-06 20:49:19 +05:30
forall ( j = 1 : lattice_maxNinteraction ) &
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_interactionTwinTwin ( j , i ) = IO_floatValue ( line , positions , 1 + j )
2009-10-06 20:49:19 +05:30
end select
endif
2009-10-10 01:58:23 +05:30
enddo
2010-08-17 19:53:55 +05:30
2009-10-06 20:49:19 +05:30
100 do i = 1 , maxNinstance
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_structure ( i ) = &
lattice_initializeStructure ( constitutive_dislotwin_structureName ( i ) , constitutive_dislotwin_CoverA ( i ) )
myStructure = constitutive_dislotwin_structure ( i )
!* Sanity checks
if ( myStructure < 1 . or . myStructure > 3 ) call IO_error ( 205 )
if ( sum ( constitutive_dislotwin_Nslip ( : , i ) ) < = 0_pInt ) call IO_error ( 225 )
if ( sum ( constitutive_dislotwin_Ntwin ( : , i ) ) < 0_pInt ) call IO_error ( 225 ) !***
do f = 1 , lattice_maxNslipFamily
2010-08-17 19:53:55 +05:30
if ( constitutive_dislotwin_Nslip ( f , i ) > 0_pInt ) then
2009-10-12 21:56:16 +05:30
if ( constitutive_dislotwin_rhoEdge0 ( f , i ) < 0.0_pReal ) call IO_error ( 220 )
if ( constitutive_dislotwin_rhoEdgeDip0 ( f , i ) < 0.0_pReal ) call IO_error ( 220 )
if ( constitutive_dislotwin_burgersPerSlipFamily ( f , i ) < = 0.0_pReal ) call IO_error ( 221 )
if ( constitutive_dislotwin_v0PerSlipFamily ( f , i ) < = 0.0_pReal ) call IO_error ( 226 )
2009-10-10 01:58:23 +05:30
endif
enddo
do f = 1 , lattice_maxNtwinFamily
2010-08-17 19:53:55 +05:30
if ( constitutive_dislotwin_Nslip ( f , i ) > 0_pInt ) then
2009-10-12 21:56:16 +05:30
if ( constitutive_dislotwin_burgersPerTwinFamily ( f , i ) < = 0.0_pReal ) call IO_error ( 221 ) !***
if ( constitutive_dislotwin_Ndot0PerTwinFamily ( f , i ) < 0.0_pReal ) call IO_error ( 226 ) !***
2009-10-10 01:58:23 +05:30
endif
enddo
! if (any(constitutive_dislotwin_interactionSlipSlip(1:maxval(lattice_interactionSlipSlip(:,:,myStructure)),i) < 1.0_pReal)) call IO_error(229)
if ( constitutive_dislotwin_CAtomicVolume ( i ) < = 0.0_pReal ) call IO_error ( 230 )
if ( constitutive_dislotwin_D0 ( i ) < = 0.0_pReal ) call IO_error ( 231 )
if ( constitutive_dislotwin_Qsd ( i ) < = 0.0_pReal ) call IO_error ( 232 )
2010-10-26 18:46:37 +05:30
if ( constitutive_dislotwin_aTolRho ( i ) < = 0.0_pReal ) call IO_error ( 233 )
2010-08-17 19:53:55 +05:30
2009-10-10 01:58:23 +05:30
!* Determine total number of active slip or twin systems
constitutive_dislotwin_Nslip ( : , i ) = min ( lattice_NslipSystem ( : , myStructure ) , constitutive_dislotwin_Nslip ( : , i ) )
constitutive_dislotwin_Ntwin ( : , i ) = min ( lattice_NtwinSystem ( : , myStructure ) , constitutive_dislotwin_Ntwin ( : , i ) )
2009-10-06 20:49:19 +05:30
constitutive_dislotwin_totalNslip ( i ) = sum ( constitutive_dislotwin_Nslip ( : , i ) )
constitutive_dislotwin_totalNtwin ( i ) = sum ( constitutive_dislotwin_Ntwin ( : , i ) )
2010-08-17 19:53:55 +05:30
enddo
2009-10-10 01:58:23 +05:30
!* Allocation of variables whose size depends on the total number of active slip systems
maxTotalNslip = maxval ( constitutive_dislotwin_totalNslip )
2010-08-17 19:53:55 +05:30
maxTotalNtwin = maxval ( constitutive_dislotwin_totalNtwin )
2009-10-12 21:56:16 +05:30
2009-10-10 01:58:23 +05:30
allocate ( constitutive_dislotwin_burgersPerSlipSystem ( maxTotalNslip , maxNinstance ) )
allocate ( constitutive_dislotwin_burgersPerTwinSystem ( maxTotalNtwin , maxNinstance ) )
allocate ( constitutive_dislotwin_QedgePerSlipSystem ( maxTotalNslip , maxNinstance ) )
allocate ( constitutive_dislotwin_v0PerSlipSystem ( maxTotalNslip , maxNinstance ) )
allocate ( constitutive_dislotwin_Ndot0PerTwinSystem ( maxTotalNtwin , maxNinstance ) )
allocate ( constitutive_dislotwin_twinsizePerTwinSystem ( maxTotalNtwin , maxNinstance ) )
allocate ( constitutive_dislotwin_CLambdaSlipPerSlipSystem ( maxTotalNslip , maxNinstance ) )
constitutive_dislotwin_burgersPerSlipSystem = 0.0_pReal
constitutive_dislotwin_burgersPerTwinSystem = 0.0_pReal
constitutive_dislotwin_QedgePerSlipSystem = 0.0_pReal
constitutive_dislotwin_v0PerSlipSystem = 0.0_pReal
constitutive_dislotwin_Ndot0PerTwinSystem = 0.0_pReal
constitutive_dislotwin_twinsizePerTwinSystem = 0.0_pReal
constitutive_dislotwin_CLambdaSlipPerSlipSystem = 0.0_pReal
allocate ( constitutive_dislotwin_interactionMatrixSlipSlip ( maxTotalNslip , maxTotalNslip , maxNinstance ) )
allocate ( constitutive_dislotwin_interactionMatrixSlipTwin ( maxTotalNslip , maxTotalNtwin , maxNinstance ) )
allocate ( constitutive_dislotwin_interactionMatrixTwinSlip ( maxTotalNtwin , maxTotalNslip , maxNinstance ) )
allocate ( constitutive_dislotwin_interactionMatrixTwinTwin ( maxTotalNtwin , maxTotalNtwin , maxNinstance ) )
allocate ( constitutive_dislotwin_forestProjectionEdge ( maxTotalNslip , maxTotalNslip , maxNinstance ) )
constitutive_dislotwin_interactionMatrixSlipSlip = 0.0_pReal
constitutive_dislotwin_interactionMatrixSlipTwin = 0.0_pReal
constitutive_dislotwin_interactionMatrixTwinSlip = 0.0_pReal
constitutive_dislotwin_interactionMatrixTwinTwin = 0.0_pReal
constitutive_dislotwin_forestProjectionEdge = 0.0_pReal
allocate ( constitutive_dislotwin_Ctwin_66 ( 6 , 6 , maxTotalNtwin , maxNinstance ) )
allocate ( constitutive_dislotwin_Ctwin_3333 ( 3 , 3 , 3 , 3 , maxTotalNtwin , maxNinstance ) )
constitutive_dislotwin_Ctwin_66 = 0.0_pReal
constitutive_dislotwin_Ctwin_3333 = 0.0_pReal
2010-08-17 19:53:55 +05:30
do i = 1 , maxNinstance
2009-10-10 01:58:23 +05:30
myStructure = constitutive_dislotwin_structure ( i )
2009-10-06 20:49:19 +05:30
2009-10-10 01:58:23 +05:30
!* Inverse lookup of my slip system family
l = 0_pInt
do f = 1 , lattice_maxNslipFamily
do k = 1 , constitutive_dislotwin_Nslip ( f , i )
l = l + 1
constitutive_dislotwin_slipFamily ( l , i ) = f
constitutive_dislotwin_slipSystemLattice ( l , i ) = sum ( lattice_NslipSystem ( 1 : f - 1 , myStructure ) ) + k
enddo ; enddo
2010-08-17 19:53:55 +05:30
2009-10-10 01:58:23 +05:30
!* Inverse lookup of my twin system family
l = 0_pInt
do f = 1 , lattice_maxNtwinFamily
do k = 1 , constitutive_dislotwin_Ntwin ( f , i )
l = l + 1
constitutive_dislotwin_twinFamily ( l , i ) = f
constitutive_dislotwin_twinSystemLattice ( l , i ) = sum ( lattice_NtwinSystem ( 1 : f - 1 , myStructure ) ) + k
enddo ; enddo
2010-08-17 19:53:55 +05:30
!* Determine size of state array
2009-10-10 01:58:23 +05:30
ns = constitutive_dislotwin_totalNslip ( i )
nt = constitutive_dislotwin_totalNtwin ( i )
constitutive_dislotwin_sizeDotState ( i ) = &
size ( constitutive_dislotwin_listBasicSlipStates ) * ns + size ( constitutive_dislotwin_listBasicTwinStates ) * nt
constitutive_dislotwin_sizeState ( i ) = &
constitutive_dislotwin_sizeDotState ( i ) + &
size ( constitutive_dislotwin_listDependentSlipStates ) * ns + size ( constitutive_dislotwin_listDependentTwinStates ) * nt
2010-08-17 19:53:55 +05:30
!* Determine size of postResults array
2009-10-10 01:58:23 +05:30
do o = 1 , maxval ( phase_Noutput )
select case ( constitutive_dislotwin_output ( o , i ) )
2009-10-12 21:56:16 +05:30
case ( 'edge_density' , &
'dipole_density' , &
'shear_rate_slip' , &
'mfp_slip' , &
'resolved_stress_slip' , &
2009-11-06 23:16:12 +05:30
'threshold_stress_slip' , &
'edge_dipole_distance' , &
2010-08-04 05:17:00 +05:30
'stress_exponent' &
2009-10-12 21:56:16 +05:30
)
mySize = constitutive_dislotwin_totalNslip ( i )
case ( 'twin_fraction' , &
'shear_rate_twin' , &
'mfp_twin' , &
'resolved_stress_twin' , &
'threshold_stress_twin' &
)
mySize = constitutive_dislotwin_totalNtwin ( i )
case default
mySize = 0_pInt
2009-10-10 01:58:23 +05:30
end select
2009-10-06 20:49:19 +05:30
2010-08-17 19:53:55 +05:30
if ( mySize > 0_pInt ) then ! any meaningful output found
2009-10-12 21:56:16 +05:30
constitutive_dislotwin_sizePostResult ( o , i ) = mySize
constitutive_dislotwin_sizePostResults ( i ) = constitutive_dislotwin_sizePostResults ( i ) + mySize
endif
2009-10-06 20:49:19 +05:30
enddo
2009-10-10 01:58:23 +05:30
!* Elasticity matrix and shear modulus according to material.config
select case ( myStructure )
2009-10-06 20:49:19 +05:30
case ( 1 : 2 ) ! cubic(s)
forall ( k = 1 : 3 )
forall ( j = 1 : 3 ) &
constitutive_dislotwin_Cslip_66 ( k , j , i ) = constitutive_dislotwin_C12 ( i )
constitutive_dislotwin_Cslip_66 ( k , k , i ) = constitutive_dislotwin_C11 ( i )
constitutive_dislotwin_Cslip_66 ( k + 3 , k + 3 , i ) = constitutive_dislotwin_C44 ( i )
end forall
case ( 3 : ) ! all hex
constitutive_dislotwin_Cslip_66 ( 1 , 1 , i ) = constitutive_dislotwin_C11 ( i )
constitutive_dislotwin_Cslip_66 ( 2 , 2 , i ) = constitutive_dislotwin_C11 ( i )
constitutive_dislotwin_Cslip_66 ( 3 , 3 , i ) = constitutive_dislotwin_C33 ( i )
constitutive_dislotwin_Cslip_66 ( 1 , 2 , i ) = constitutive_dislotwin_C12 ( i )
constitutive_dislotwin_Cslip_66 ( 2 , 1 , i ) = constitutive_dislotwin_C12 ( i )
constitutive_dislotwin_Cslip_66 ( 1 , 3 , i ) = constitutive_dislotwin_C13 ( i )
constitutive_dislotwin_Cslip_66 ( 3 , 1 , i ) = constitutive_dislotwin_C13 ( i )
constitutive_dislotwin_Cslip_66 ( 2 , 3 , i ) = constitutive_dislotwin_C13 ( i )
constitutive_dislotwin_Cslip_66 ( 3 , 2 , i ) = constitutive_dislotwin_C13 ( i )
constitutive_dislotwin_Cslip_66 ( 4 , 4 , i ) = constitutive_dislotwin_C44 ( i )
constitutive_dislotwin_Cslip_66 ( 5 , 5 , i ) = constitutive_dislotwin_C44 ( i )
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_Cslip_66 ( 6 , 6 , i ) = 0.5_pReal * ( constitutive_dislotwin_C11 ( i ) - constitutive_dislotwin_C12 ( i ) )
2009-10-06 20:49:19 +05:30
end select
constitutive_dislotwin_Cslip_66 ( : , : , i ) = math_Mandel3333to66 ( math_Voigt66to3333 ( constitutive_dislotwin_Cslip_66 ( : , : , i ) ) )
constitutive_dislotwin_Cslip_3333 ( : , : , : , : , i ) = math_Voigt66to3333 ( constitutive_dislotwin_Cslip_66 ( : , : , i ) )
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_Gmod ( i ) = &
0.2_pReal * ( constitutive_dislotwin_C11 ( i ) - constitutive_dislotwin_C12 ( i ) ) + 0.3_pReal * constitutive_dislotwin_C44 ( i )
2010-08-17 19:53:55 +05:30
2009-10-06 20:49:19 +05:30
!* Construction of the twin elasticity matrices
do j = 1 , lattice_maxNtwinFamily
2010-08-17 19:53:55 +05:30
do k = 1 , constitutive_dislotwin_Ntwin ( j , i )
2009-10-06 20:49:19 +05:30
do l = 1 , 3 ; do m = 1 , 3 ; do n = 1 , 3 ; do o = 1 , 3 ; do p = 1 , 3 ; do q = 1 , 3 ; do r = 1 , 3 ; do s = 1 , 3
2009-10-12 21:56:16 +05:30
constitutive_dislotwin_Ctwin_3333 ( l , m , n , o , sum ( constitutive_dislotwin_Nslip ( 1 : j - 1 , i ) ) + k , i ) = &
constitutive_dislotwin_Ctwin_3333 ( l , m , n , o , sum ( constitutive_dislotwin_Nslip ( 1 : j - 1 , i ) ) + k , i ) + &
constitutive_dislotwin_Cslip_3333 ( p , q , r , s , i ) * &
lattice_Qtwin ( l , p , sum ( lattice_NslipSystem ( 1 : j - 1 , myStructure ) ) + k , myStructure ) * &
lattice_Qtwin ( m , q , sum ( lattice_NslipSystem ( 1 : j - 1 , myStructure ) ) + k , myStructure ) * &
lattice_Qtwin ( n , r , sum ( lattice_NslipSystem ( 1 : j - 1 , myStructure ) ) + k , myStructure ) * &
lattice_Qtwin ( o , s , sum ( lattice_NslipSystem ( 1 : j - 1 , myStructure ) ) + k , myStructure )
enddo ; enddo ; enddo ; enddo ; enddo ; enddo ; enddo ; enddo
2009-10-06 20:49:19 +05:30
constitutive_dislotwin_Ctwin_66 ( : , : , k , i ) = math_Mandel3333to66 ( constitutive_dislotwin_Ctwin_3333 ( : , : , : , : , k , i ) )
2009-10-12 21:56:16 +05:30
enddo
2009-10-06 20:49:19 +05:30
enddo
2010-08-17 19:53:55 +05:30
!* Burgers vector, dislocation velocity prefactor, mean free path prefactor and minimum dipole distance for each slip system
do s = 1 , constitutive_dislotwin_totalNslip ( i )
f = constitutive_dislotwin_slipFamily ( s , i )
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_burgersPerSlipSystem ( s , i ) = constitutive_dislotwin_burgersPerSlipFamily ( f , i )
constitutive_dislotwin_QedgePerSlipSystem ( s , i ) = constitutive_dislotwin_QedgePerSlipFamily ( f , i )
constitutive_dislotwin_v0PerSlipSystem ( s , i ) = constitutive_dislotwin_v0PerSlipFamily ( f , i )
2009-10-12 21:56:16 +05:30
constitutive_dislotwin_CLambdaSlipPerSlipSystem ( s , i ) = constitutive_dislotwin_CLambdaSlipPerSlipFamily ( f , i )
2010-08-17 19:53:55 +05:30
enddo
!* Burgers vector, nucleation rate prefactor and twin size for each twin system
do s = 1 , constitutive_dislotwin_totalNtwin ( i )
f = constitutive_dislotwin_twinFamily ( s , i )
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_burgersPerTwinSystem ( s , i ) = constitutive_dislotwin_burgersPerTwinFamily ( f , i )
constitutive_dislotwin_Ndot0PerTwinSystem ( s , i ) = constitutive_dislotwin_Ndot0PerTwinFamily ( f , i )
constitutive_dislotwin_twinsizePerTwinSystem ( s , i ) = constitutive_dislotwin_twinsizePerTwinFamily ( f , i )
2010-08-17 19:53:55 +05:30
enddo
2009-10-10 01:58:23 +05:30
!* Construction of interaction matrices
do s1 = 1 , constitutive_dislotwin_totalNslip ( i )
2010-08-17 19:53:55 +05:30
do s2 = 1 , constitutive_dislotwin_totalNslip ( i )
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_interactionMatrixSlipSlip ( s1 , s2 , i ) = &
constitutive_dislotwin_interactionSlipSlip ( lattice_interactionSlipSlip ( constitutive_dislotwin_slipSystemLattice ( s1 , i ) , &
constitutive_dislotwin_slipSystemLattice ( s2 , i ) , &
2009-10-12 21:56:16 +05:30
myStructure ) , i )
2009-10-10 01:58:23 +05:30
enddo ; enddo
2010-08-17 19:53:55 +05:30
2009-10-10 01:58:23 +05:30
do s1 = 1 , constitutive_dislotwin_totalNslip ( i )
2010-08-17 19:53:55 +05:30
do t2 = 1 , constitutive_dislotwin_totalNtwin ( i )
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_interactionMatrixSlipTwin ( s1 , t2 , i ) = &
constitutive_dislotwin_interactionSlipTwin ( lattice_interactionSlipTwin ( constitutive_dislotwin_slipSystemLattice ( s1 , i ) , &
constitutive_dislotwin_twinSystemLattice ( t2 , i ) , &
2010-08-17 19:53:55 +05:30
myStructure ) , i )
2009-10-10 01:58:23 +05:30
enddo ; enddo
2010-08-17 19:53:55 +05:30
2009-10-10 01:58:23 +05:30
do t1 = 1 , constitutive_dislotwin_totalNtwin ( i )
2010-08-17 19:53:55 +05:30
do s2 = 1 , constitutive_dislotwin_totalNslip ( i )
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_interactionMatrixTwinSlip ( t1 , s2 , i ) = &
constitutive_dislotwin_interactionTwinSlip ( lattice_interactionTwinSlip ( constitutive_dislotwin_twinSystemLattice ( t1 , i ) , &
constitutive_dislotwin_slipSystemLattice ( s2 , i ) , &
2010-08-17 19:53:55 +05:30
myStructure ) , i )
2009-10-10 01:58:23 +05:30
enddo ; enddo
2009-10-06 20:49:19 +05:30
2009-10-10 01:58:23 +05:30
do t1 = 1 , constitutive_dislotwin_totalNtwin ( i )
2010-08-17 19:53:55 +05:30
do t2 = 1 , constitutive_dislotwin_totalNtwin ( i )
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_interactionMatrixTwinTwin ( t1 , t2 , i ) = &
constitutive_dislotwin_interactionTwinTwin ( lattice_interactionTwinTwin ( constitutive_dislotwin_twinSystemLattice ( t1 , i ) , &
constitutive_dislotwin_twinSystemLattice ( t2 , i ) , &
2010-08-17 19:53:55 +05:30
myStructure ) , i )
2009-10-10 01:58:23 +05:30
enddo ; enddo
2010-08-17 19:53:55 +05:30
!* Calculation of forest projections for edge dislocations
2009-10-10 01:58:23 +05:30
do s1 = 1 , constitutive_dislotwin_totalNslip ( i )
2010-08-17 19:53:55 +05:30
do s2 = 1 , constitutive_dislotwin_totalNslip ( i )
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_forestProjectionEdge ( s1 , s2 , i ) = &
abs ( math_mul3x3 ( lattice_sn ( : , constitutive_dislotwin_slipSystemLattice ( s1 , i ) , myStructure ) , &
2010-08-17 19:53:55 +05:30
lattice_st ( : , constitutive_dislotwin_slipSystemLattice ( s2 , i ) , myStructure ) ) )
2009-10-10 01:58:23 +05:30
enddo ; enddo
2010-08-17 19:53:55 +05:30
2009-10-10 01:58:23 +05:30
enddo
2009-10-06 20:49:19 +05:30
2009-10-10 01:58:23 +05:30
return
2009-10-06 20:49:19 +05:30
end subroutine
function constitutive_dislotwin_stateInit ( myInstance )
!*********************************************************************
!* initial microstructural state *
!*********************************************************************
2009-10-12 21:56:16 +05:30
use prec , only : pReal , pInt
2009-10-10 01:58:23 +05:30
use math , only : pi
use lattice , only : lattice_maxNslipFamily , lattice_maxNtwinFamily
implicit none
!* Input-Output variables
integer ( pInt ) :: myInstance
real ( pReal ) , dimension ( constitutive_dislotwin_sizeState ( myInstance ) ) :: constitutive_dislotwin_stateInit
!* Local variables
integer ( pInt ) s0 , s1 , s , t , f , ns , nt
2009-10-12 21:56:16 +05:30
real ( pReal ) , dimension ( constitutive_dislotwin_totalNslip ( myInstance ) ) :: rhoEdge0 , &
rhoEdgeDip0 , &
invLambdaSlip0 , &
MeanFreePathSlip0 , &
tauSlipThreshold0
real ( pReal ) , dimension ( constitutive_dislotwin_totalNtwin ( myInstance ) ) :: MeanFreePathTwin0 , TwinVolume0
2009-10-10 01:58:23 +05:30
2009-10-12 21:56:16 +05:30
ns = constitutive_dislotwin_totalNslip ( myInstance )
2009-10-10 01:58:23 +05:30
nt = constitutive_dislotwin_totalNtwin ( myInstance )
constitutive_dislotwin_stateInit = 0.0_pReal
2009-10-12 21:56:16 +05:30
!* Initialize basic slip state variables
s1 = 0_pInt
do f = 1 , lattice_maxNslipFamily
s0 = s1 + 1_pInt
2010-08-17 19:53:55 +05:30
s1 = s0 + constitutive_dislotwin_Nslip ( f , myInstance ) - 1_pInt
2009-10-12 21:56:16 +05:30
do s = s0 , s1
rhoEdge0 ( s ) = constitutive_dislotwin_rhoEdge0 ( f , myInstance )
rhoEdgeDip0 ( s ) = constitutive_dislotwin_rhoEdgeDip0 ( f , myInstance )
2010-08-17 19:53:55 +05:30
enddo
2009-10-12 21:56:16 +05:30
enddo
constitutive_dislotwin_stateInit ( 1 : ns ) = rhoEdge0
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_stateInit ( ns + 1 : 2 * ns ) = rhoEdgeDip0
!* Initialize dependent slip microstructural variables
forall ( s = 1 : ns ) &
invLambdaSlip0 ( s ) = sqrt ( dot_product ( ( rhoEdge0 + rhoEdgeDip0 ) , constitutive_dislotwin_forestProjectionEdge ( 1 : ns , s , myInstance ) ) ) / &
2010-08-17 19:53:55 +05:30
constitutive_dislotwin_CLambdaSlipPerSlipSystem ( s , myInstance )
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_stateInit ( 2 * ns + nt + 1 : 3 * ns + nt ) = invLambdaSlip0
forall ( s = 1 : ns ) &
MeanFreePathSlip0 ( s ) = &
2010-08-17 19:53:55 +05:30
constitutive_dislotwin_GrainSize ( myInstance ) / ( 1.0_pReal + invLambdaSlip0 ( s ) * constitutive_dislotwin_GrainSize ( myInstance ) )
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_stateInit ( 4 * ns + 2 * nt + 1 : 5 * ns + 2 * nt ) = MeanFreePathSlip0
forall ( s = 1 : ns ) &
2010-08-17 19:53:55 +05:30
tauSlipThreshold0 ( s ) = constitutive_dislotwin_SolidSolutionStrength ( myInstance ) + &
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_Gmod ( myInstance ) * constitutive_dislotwin_burgersPerSlipSystem ( s , myInstance ) * &
2009-10-13 20:35:21 +05:30
sqrt ( dot_product ( ( rhoEdge0 + rhoEdgeDip0 ) , constitutive_dislotwin_interactionMatrixSlipSlip ( 1 : ns , s , myInstance ) ) )
2009-10-12 21:56:16 +05:30
constitutive_dislotwin_stateInit ( 5 * ns + 3 * nt + 1 : 6 * ns + 3 * nt ) = tauSlipThreshold0
2009-10-10 01:58:23 +05:30
!* Initialize dependent twin microstructural variables
forall ( t = 1 : nt ) &
MeanFreePathTwin0 ( t ) = constitutive_dislotwin_GrainSize ( myInstance )
constitutive_dislotwin_stateInit ( 5 * ns + 2 * nt + 1 : 5 * ns + 3 * nt ) = MeanFreePathTwin0
forall ( t = 1 : nt ) &
2010-08-17 19:53:55 +05:30
TwinVolume0 ( t ) = &
2009-10-10 01:58:23 +05:30
( pi / 6.0_pReal ) * constitutive_dislotwin_twinsizePerTwinSystem ( t , myInstance ) * MeanFreePathTwin0 ( t ) ** ( 2.0_pReal )
constitutive_dislotwin_stateInit ( 6 * ns + 4 * nt + 1 : 6 * ns + 5 * nt ) = TwinVolume0
2009-10-12 21:56:16 +05:30
!write(6,*) '#STATEINIT#'
!write(6,*)
!write(6,'(a,/,4(3(f30.20,x)/))') 'RhoEdge',rhoEdge0
!write(6,'(a,/,4(3(f30.20,x)/))') 'RhoEdgedip',rhoEdgeDip0
!write(6,'(a,/,4(3(f30.20,x)/))') 'invLambdaSlip',invLambdaSlip0
!write(6,'(a,/,4(3(f30.20,x)/))') 'MeanFreePathSlip',MeanFreePathSlip0
!write(6,'(a,/,4(3(f30.20,x)/))') 'tauSlipThreshold', tauSlipThreshold0
!write(6,'(a,/,4(3(f30.20,x)/))') 'MeanFreePathTwin', MeanFreePathTwin0
!write(6,'(a,/,4(3(f30.20,x)/))') 'TwinVolume', TwinVolume0
2009-10-10 01:58:23 +05:30
return
2009-10-06 20:49:19 +05:30
end function
2010-10-26 18:46:37 +05:30
pure function constitutive_dislotwin_aTolState ( myInstance )
2009-10-06 20:49:19 +05:30
!*********************************************************************
2010-10-26 18:46:37 +05:30
!* absolute state tolerance *
2009-10-06 20:49:19 +05:30
!*********************************************************************
2009-10-10 01:58:23 +05:30
use prec , only : pReal , pInt
2009-10-06 20:49:19 +05:30
implicit none
2009-10-10 01:58:23 +05:30
!* Input-Output variables
integer ( pInt ) , intent ( in ) :: myInstance
2010-10-26 18:46:37 +05:30
real ( pReal ) , dimension ( constitutive_dislotwin_sizeState ( myInstance ) ) :: constitutive_dislotwin_aTolState
2009-10-06 20:49:19 +05:30
2010-10-26 18:46:37 +05:30
constitutive_dislotwin_aTolState = constitutive_dislotwin_aTolRho ( myInstance )
2009-10-06 20:49:19 +05:30
2009-10-10 01:58:23 +05:30
return
2009-10-06 20:49:19 +05:30
endfunction
2009-10-10 01:58:23 +05:30
pure function constitutive_dislotwin_homogenizedC ( state , g , ip , el )
2009-10-06 20:49:19 +05:30
!*********************************************************************
!* calculates homogenized elacticity matrix *
!* - state : microstructure quantities *
2009-10-10 01:58:23 +05:30
!* - g : component-ID of current integration point *
2009-10-06 20:49:19 +05:30
!* - ip : current integration point *
!* - el : current element *
!*********************************************************************
2009-10-10 01:58:23 +05:30
use prec , only : pReal , pInt , p_vec
use mesh , only : mesh_NcpElems , mesh_maxNips
use material , only : homogenization_maxNgrains , material_phase , phase_constitutionInstance
implicit none
!* Input-Output variables
integer ( pInt ) , intent ( in ) :: g , ip , el
type ( p_vec ) , dimension ( homogenization_maxNgrains , mesh_maxNips , mesh_NcpElems ) , intent ( in ) :: state
real ( pReal ) , dimension ( 6 , 6 ) :: constitutive_dislotwin_homogenizedC
2010-08-17 19:53:55 +05:30
!* Local variables
2009-10-10 01:58:23 +05:30
integer ( pInt ) myInstance , ns , nt , i
real ( pReal ) sumf
2010-08-17 19:53:55 +05:30
2009-10-10 01:58:23 +05:30
!* Shortened notation
myInstance = phase_constitutionInstance ( material_phase ( g , ip , el ) )
ns = constitutive_dislotwin_totalNslip ( myInstance )
nt = constitutive_dislotwin_totalNtwin ( myInstance )
2009-10-12 21:56:16 +05:30
!* Total twin volume fraction
sumf = sum ( state ( g , ip , el ) % p ( ( 2 * ns + 1 ) : ( 2 * ns + nt ) ) ) ! safe for nt == 0
2009-10-10 01:58:23 +05:30
!* Homogenized elasticity matrix
constitutive_dislotwin_homogenizedC = ( 1.0_pReal - sumf ) * constitutive_dislotwin_Cslip_66 ( : , : , myInstance )
do i = 1 , nt
2009-10-12 21:56:16 +05:30
constitutive_dislotwin_homogenizedC = &
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_homogenizedC + state ( g , ip , el ) % p ( 2 * ns + i ) * constitutive_dislotwin_Ctwin_66 ( : , : , i , myInstance )
2010-08-17 19:53:55 +05:30
enddo
2009-10-10 01:58:23 +05:30
return
2009-10-06 20:49:19 +05:30
end function
2009-10-10 01:58:23 +05:30
subroutine constitutive_dislotwin_microstructure ( Temperature , state , g , ip , el )
2009-10-06 20:49:19 +05:30
!*********************************************************************
!* calculates quantities characterizing the microstructure *
!* - Temperature : temperature *
!* - state : microstructure quantities *
!* - ipc : component-ID of current integration point *
!* - ip : current integration point *
!* - el : current element *
!*********************************************************************
2009-10-10 01:58:23 +05:30
use prec , only : pReal , pInt , p_vec
use math , only : pi
use mesh , only : mesh_NcpElems , mesh_maxNips
use material , only : homogenization_maxNgrains , material_phase , phase_constitutionInstance
use lattice , only : lattice_interactionSlipTwin , lattice_interactionTwinTwin
!use debug, only: debugger
implicit none
!* Input-Output variables
integer ( pInt ) , intent ( in ) :: g , ip , el
real ( pReal ) , intent ( in ) :: Temperature
type ( p_vec ) , dimension ( homogenization_maxNgrains , mesh_maxNips , mesh_NcpElems ) , intent ( inout ) :: state
!* Local variables
2010-08-17 19:53:55 +05:30
integer ( pInt ) myInstance , myStructure , ns , nt , s , t
2009-10-10 01:58:23 +05:30
real ( pReal ) sumf , sfe
real ( pReal ) , dimension ( constitutive_dislotwin_totalNtwin ( phase_constitutionInstance ( material_phase ( g , ip , el ) ) ) ) :: fOverStacksize
2010-08-17 19:53:55 +05:30
2009-10-10 01:58:23 +05:30
!* Shortened notation
myInstance = phase_constitutionInstance ( material_phase ( g , ip , el ) )
myStructure = constitutive_dislotwin_structure ( myInstance )
ns = constitutive_dislotwin_totalNslip ( myInstance )
nt = constitutive_dislotwin_totalNtwin ( myInstance )
2009-10-12 21:56:16 +05:30
!* State: 1 : ns rho_edge
2009-10-10 01:58:23 +05:30
!* State: ns+1 : 2*ns rho_dipole
!* State: 2*ns+1 : 2*ns+nt f
!* State: 2*ns+nt+1 : 3*ns+nt 1/lambda_slip
!* State: 3*ns+nt+1 : 4*ns+nt 1/lambda_sliptwin
!* State: 4*ns+nt+1 : 4*ns+2*nt 1/lambda_twin
!* State: 4*ns+2*nt+1 : 5*ns+2*nt mfp_slip
!* State: 5*ns+2*nt+1 : 5*ns+3*nt mfp_twin
!* State: 5*ns+3*nt+1 : 6*ns+3*nt threshold_stress_slip
!* State: 6*ns+3*nt+1 : 6*ns+4*nt threshold_stress_twin
2009-10-12 21:56:16 +05:30
!* State: 6*ns+4*nt+1 : 6*ns+5*nt twin volume
2009-10-10 01:58:23 +05:30
2009-10-12 21:56:16 +05:30
!* Total twin volume fraction
sumf = sum ( state ( g , ip , el ) % p ( ( 2 * ns + 1 ) : ( 2 * ns + nt ) ) ) ! safe for nt == 0
!* Stacking fault energy
sfe = 0.0002_pReal * Temperature - 0.0396_pReal
2009-10-10 01:58:23 +05:30
!* rescaled twin volume fraction for topology
forall ( t = 1 : nt ) &
2009-10-12 21:56:16 +05:30
fOverStacksize ( t ) = &
state ( g , ip , el ) % p ( 2 * ns + t ) / constitutive_dislotwin_twinsizePerTwinSystem ( t , myInstance )
2010-08-17 19:53:55 +05:30
2009-10-10 01:58:23 +05:30
!* 1/mean free distance between 2 forest dislocations seen by a moving dislocation
forall ( s = 1 : ns ) &
2009-10-12 21:56:16 +05:30
state ( g , ip , el ) % p ( 2 * ns + nt + s ) = &
sqrt ( dot_product ( ( state ( g , ip , el ) % p ( 1 : ns ) + state ( g , ip , el ) % p ( ns + 1 : 2 * ns ) ) , &
constitutive_dislotwin_forestProjectionEdge ( 1 : ns , s , myInstance ) ) ) / &
2010-08-17 19:53:55 +05:30
constitutive_dislotwin_CLambdaSlipPerSlipSystem ( s , myInstance )
2009-10-10 01:58:23 +05:30
!* 1/mean free distance between 2 twin stacks from different systems seen by a moving dislocation
!$OMP CRITICAL (evilmatmul)
state ( g , ip , el ) % p ( ( 3 * ns + nt + 1 ) : ( 4 * ns + nt ) ) = 0.0_pReal
2009-10-12 21:56:16 +05:30
if ( nt > 0_pInt ) &
state ( g , ip , el ) % p ( ( 3 * ns + nt + 1 ) : ( 4 * ns + nt ) ) = &
matmul ( constitutive_dislotwin_interactionMatrixSlipTwin ( 1 : ns , 1 : nt , myInstance ) , fOverStacksize ( 1 : nt ) ) / ( 1.0_pReal - sumf )
2009-10-10 01:58:23 +05:30
!$OMP END CRITICAL (evilmatmul)
!* 1/mean free distance between 2 twin stacks from different systems seen by a growing twin
!$OMP CRITICAL (evilmatmul)
2009-10-12 21:56:16 +05:30
if ( nt > 0_pInt ) &
state ( g , ip , el ) % p ( ( 4 * ns + nt + 1 ) : ( 4 * ns + 2 * nt ) ) = &
matmul ( constitutive_dislotwin_interactionMatrixTwinTwin ( 1 : nt , 1 : nt , myInstance ) , fOverStacksize ( 1 : nt ) ) / ( 1.0_pReal - sumf )
2009-10-10 01:58:23 +05:30
!$OMP END CRITICAL (evilmatmul)
!* mean free path between 2 obstacles seen by a moving dislocation
do s = 1 , ns
if ( nt > 0_pInt ) then
2009-10-12 21:56:16 +05:30
state ( g , ip , el ) % p ( 4 * ns + 2 * nt + s ) = &
constitutive_dislotwin_GrainSize ( myInstance ) / ( 1.0_pReal + constitutive_dislotwin_GrainSize ( myInstance ) * &
2010-08-17 19:53:55 +05:30
( state ( g , ip , el ) % p ( 2 * ns + nt + s ) + state ( g , ip , el ) % p ( 3 * ns + nt + s ) ) )
2009-10-10 01:58:23 +05:30
else
2009-10-12 21:56:16 +05:30
state ( g , ip , el ) % p ( 4 * ns + s ) = &
constitutive_dislotwin_GrainSize ( myInstance ) / &
( 1.0_pReal + constitutive_dislotwin_GrainSize ( myInstance ) * ( state ( g , ip , el ) % p ( 2 * ns + s ) ) )
2009-10-10 01:58:23 +05:30
endif
enddo
!* mean free path between 2 obstacles seen by a growing twin
forall ( t = 1 : nt ) &
2009-10-12 21:56:16 +05:30
state ( g , ip , el ) % p ( 5 * ns + 2 * nt + t ) = &
( constitutive_dislotwin_Cmfptwin ( myInstance ) * constitutive_dislotwin_GrainSize ( myInstance ) ) / &
2010-08-17 19:53:55 +05:30
( 1.0_pReal + constitutive_dislotwin_GrainSize ( myInstance ) * state ( g , ip , el ) % p ( 4 * ns + nt + t ) )
2009-10-10 01:58:23 +05:30
!* threshold stress for dislocation motion
forall ( s = 1 : ns ) &
2010-08-17 19:53:55 +05:30
state ( g , ip , el ) % p ( 5 * ns + 3 * nt + s ) = constitutive_dislotwin_SolidSolutionStrength ( myInstance ) + &
2009-10-12 21:56:16 +05:30
constitutive_dislotwin_Gmod ( myInstance ) * constitutive_dislotwin_burgersPerSlipSystem ( s , myInstance ) * &
2009-10-13 20:35:21 +05:30
sqrt ( dot_product ( ( state ( g , ip , el ) % p ( 1 : ns ) + state ( g , ip , el ) % p ( ns + 1 : 2 * ns ) ) , &
2010-08-04 05:17:00 +05:30
constitutive_dislotwin_interactionMatrixSlipSlip ( 1 : ns , s , myInstance ) ) )
2009-10-12 21:56:16 +05:30
2009-10-10 01:58:23 +05:30
!* threshold stress for growing twin
forall ( t = 1 : nt ) &
2009-10-12 21:56:16 +05:30
state ( g , ip , el ) % p ( 6 * ns + 3 * nt + t ) = &
constitutive_dislotwin_Cthresholdtwin ( myInstance ) * &
( sfe / ( 3.0_pReal * constitutive_dislotwin_burgersPerTwinSystem ( t , myInstance ) ) + &
2010-03-24 13:47:27 +05:30
3.0_pReal * constitutive_dislotwin_burgersPerTwinSystem ( t , myInstance ) * constitutive_dislotwin_Gmod ( myInstance ) / &
2010-08-17 19:53:55 +05:30
( constitutive_dislotwin_L0 ( myInstance ) * constitutive_dislotwin_burgersPerSlipSystem ( t , myInstance ) ) )
2009-10-10 01:58:23 +05:30
!* final twin volume after growth
forall ( t = 1 : nt ) &
2009-10-12 21:56:16 +05:30
state ( g , ip , el ) % p ( 6 * ns + 4 * nt + t ) = &
( pi / 6.0_pReal ) * constitutive_dislotwin_twinsizePerTwinSystem ( t , myInstance ) * state ( g , ip , el ) % p ( 5 * ns + 2 * nt + t ) ** ( 2.0_pReal )
!if ((ip==1).and.(el==1)) then
! write(6,*) '#MICROSTRUCTURE#'
! write(6,*)
! write(6,'(a,/,4(3(f10.4,x)/))') 'rhoEdge',state(g,ip,el)%p(1:ns)/1e9
! write(6,'(a,/,4(3(f10.4,x)/))') 'rhoEdgeDip',state(g,ip,el)%p(ns+1:2*ns)/1e9
! write(6,'(a,/,4(3(f10.4,x)/))') 'Fraction',state(g,ip,el)%p(2*ns+1:2*ns+nt)
!endif
2009-10-10 01:58:23 +05:30
return
2009-10-06 20:49:19 +05:30
end subroutine
2009-10-10 01:58:23 +05:30
subroutine constitutive_dislotwin_LpAndItsTangent ( Lp , dLp_dTstar , Tstar_v , Temperature , state , g , ip , el )
2009-10-06 20:49:19 +05:30
!*********************************************************************
!* calculates plastic velocity gradient and its tangent *
!* INPUT: *
!* - Temperature : temperature *
!* - state : microstructure quantities *
!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) *
!* - ipc : component-ID at current integration point *
!* - ip : current integration point *
!* - el : current element *
!* OUTPUT: *
!* - Lp : plastic velocity gradient *
!* - dLp_dTstar : derivative of Lp (4th-rank tensor) *
!*********************************************************************
2009-10-10 01:58:23 +05:30
use prec , only : pReal , pInt , p_vec
use math , only : math_Plain3333to99
use mesh , only : mesh_NcpElems , mesh_maxNips
use material , only : homogenization_maxNgrains , material_phase , phase_constitutionInstance
use lattice , only : lattice_Sslip , lattice_Sslip_v , lattice_Stwin , lattice_Stwin_v , lattice_maxNslipFamily , lattice_maxNtwinFamily , &
2009-10-12 21:56:16 +05:30
lattice_NslipSystem , lattice_NtwinSystem , lattice_shearTwin
2009-10-10 01:58:23 +05:30
implicit none
!* Input-Output variables
integer ( pInt ) , intent ( in ) :: g , ip , el
real ( pReal ) , intent ( in ) :: Temperature
real ( pReal ) , dimension ( 6 ) , intent ( in ) :: Tstar_v
type ( p_vec ) , dimension ( homogenization_maxNgrains , mesh_maxNips , mesh_NcpElems ) , intent ( inout ) :: state
real ( pReal ) , dimension ( 3 , 3 ) , intent ( out ) :: Lp
real ( pReal ) , dimension ( 9 , 9 ) , intent ( out ) :: dLp_dTstar
!* Local variables
integer ( pInt ) myInstance , myStructure , ns , nt , f , i , j , k , l , m , n , index_myFamily
real ( pReal ) sumf , StressRatio_p , StressRatio_pminus1 , StressRatio_r , BoltzmannRatio , DotGamma0
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: dLp_dTstar3333
real ( pReal ) , dimension ( constitutive_dislotwin_totalNslip ( phase_constitutionInstance ( material_phase ( g , ip , el ) ) ) ) :: &
2009-10-06 20:49:19 +05:30
gdot_slip , dgdot_dtauslip , tau_slip
2009-10-10 01:58:23 +05:30
real ( pReal ) , dimension ( constitutive_dislotwin_totalNtwin ( phase_constitutionInstance ( material_phase ( g , ip , el ) ) ) ) :: &
2009-10-06 20:49:19 +05:30
gdot_twin , dgdot_dtautwin , tau_twin
2009-10-10 01:58:23 +05:30
!* Shortened notation
myInstance = phase_constitutionInstance ( material_phase ( g , ip , el ) )
2010-08-17 19:53:55 +05:30
myStructure = constitutive_dislotwin_structure ( myInstance )
2009-10-10 01:58:23 +05:30
ns = constitutive_dislotwin_totalNslip ( myInstance )
nt = constitutive_dislotwin_totalNtwin ( myInstance )
2009-10-12 21:56:16 +05:30
!* Total twin volume fraction
sumf = sum ( state ( g , ip , el ) % p ( ( 2 * ns + 1 ) : ( 2 * ns + nt ) ) ) ! safe for nt == 0
2009-10-10 01:58:23 +05:30
Lp = 0.0_pReal
dLp_dTstar3333 = 0.0_pReal
dLp_dTstar = 0.0_pReal
!* Dislocation glide part
gdot_slip = 0.0_pReal
dgdot_dtauslip = 0.0_pReal
j = 0_pInt
do f = 1 , lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum ( lattice_NslipSystem ( 1 : f - 1 , myStructure ) ) ! at which index starts my family
do i = 1 , constitutive_dislotwin_Nslip ( f , myInstance ) ! process each (active) slip system in family
j = j + 1_pInt
2009-10-12 21:56:16 +05:30
!* Calculation of Lp
!* Resolved shear stress on slip system
2010-08-17 19:53:55 +05:30
tau_slip ( j ) = dot_product ( Tstar_v , lattice_Sslip_v ( : , index_myFamily + i , myStructure ) )
2009-10-12 21:56:16 +05:30
!* Stress ratios
StressRatio_p = ( abs ( tau_slip ( j ) ) / state ( g , ip , el ) % p ( 5 * ns + 3 * nt + j ) ) ** constitutive_dislotwin_p ( myInstance )
StressRatio_pminus1 = ( abs ( tau_slip ( j ) ) / state ( g , ip , el ) % p ( 5 * ns + 3 * nt + j ) ) ** ( constitutive_dislotwin_p ( myInstance ) - 1.0_pReal )
!* Boltzmann ratio
BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem ( f , myInstance ) / ( kB * Temperature )
!* Initial shear rates
DotGamma0 = &
state ( g , ip , el ) % p ( j ) * constitutive_dislotwin_burgersPerSlipSystem ( f , myInstance ) * &
constitutive_dislotwin_v0PerSlipSystem ( f , myInstance )
!* Shear rates due to slip
2010-08-04 05:17:00 +05:30
gdot_slip ( j ) = DotGamma0 * exp ( - BoltzmannRatio * ( 1 - StressRatio_p ) ** constitutive_dislotwin_q ( myInstance ) ) * &
sign ( 1.0_pReal , tau_slip ( j ) )
2009-10-12 21:56:16 +05:30
!* Derivatives of shear rates
dgdot_dtauslip ( j ) = &
2009-11-06 23:16:12 +05:30
( ( abs ( gdot_slip ( j ) ) * BoltzmannRatio * &
2009-10-12 21:56:16 +05:30
constitutive_dislotwin_p ( myInstance ) * constitutive_dislotwin_q ( myInstance ) ) / state ( g , ip , el ) % p ( 5 * ns + 3 * nt + j ) ) * &
StressRatio_pminus1 * ( 1 - StressRatio_p ) ** ( constitutive_dislotwin_q ( myInstance ) - 1.0_pReal )
2009-10-10 01:58:23 +05:30
!* Plastic velocity gradient for dislocation glide
Lp = Lp + ( 1.0_pReal - sumf ) * gdot_slip ( j ) * lattice_Sslip ( : , : , index_myFamily + i , myStructure )
2009-10-12 21:56:16 +05:30
!* Calculation of the tangent of Lp
2009-10-10 01:58:23 +05:30
forall ( k = 1 : 3 , l = 1 : 3 , m = 1 : 3 , n = 1 : 3 ) &
2009-10-12 21:56:16 +05:30
dLp_dTstar3333 ( k , l , m , n ) = &
dLp_dTstar3333 ( k , l , m , n ) + dgdot_dtauslip ( j ) * &
lattice_Sslip ( k , l , index_myFamily + i , myStructure ) * &
2010-08-17 19:53:55 +05:30
lattice_Sslip ( m , n , index_myFamily + i , myStructure )
2009-10-10 01:58:23 +05:30
enddo
enddo
!* Mechanical twinning part
gdot_twin = 0.0_pReal
dgdot_dtautwin = 0.0_pReal
j = 0_pInt
do f = 1 , lattice_maxNtwinFamily ! loop over all slip families
index_myFamily = sum ( lattice_NtwinSystem ( 1 : f - 1 , myStructure ) ) ! at which index starts my family
do i = 1 , constitutive_dislotwin_Ntwin ( f , myInstance ) ! process each (active) slip system in family
j = j + 1_pInt
2009-10-12 21:56:16 +05:30
!* Calculation of Lp
!* Resolved shear stress on twin system
2010-08-17 19:53:55 +05:30
tau_twin ( j ) = dot_product ( Tstar_v , lattice_Stwin_v ( : , index_myFamily + i , myStructure ) )
2010-08-04 05:17:00 +05:30
!* Stress ratios
2010-08-17 19:53:55 +05:30
StressRatio_r = ( state ( g , ip , el ) % p ( 6 * ns + 3 * nt + j ) / tau_twin ( j ) ) ** constitutive_dislotwin_r ( myInstance )
2010-08-04 05:17:00 +05:30
!* Shear rates and their derivatives due to twin
2010-08-17 19:53:55 +05:30
if ( tau_twin ( j ) > 0.0_pReal ) then
2009-10-12 21:56:16 +05:30
gdot_twin ( j ) = &
( constitutive_dislotwin_MaxTwinFraction ( myInstance ) - sumf ) * lattice_shearTwin ( index_myFamily + i , myStructure ) * &
2010-08-17 19:53:55 +05:30
state ( g , ip , el ) % p ( 6 * ns + 4 * nt + j ) * constitutive_dislotwin_Ndot0PerTwinSystem ( f , myInstance ) * exp ( - StressRatio_r )
2009-10-12 21:56:16 +05:30
dgdot_dtautwin ( j ) = ( ( gdot_twin ( j ) * constitutive_dislotwin_r ( myInstance ) ) / tau_twin ( j ) ) * StressRatio_r
endif
2010-08-17 19:53:55 +05:30
!* Plastic velocity gradient for mechanical twinning
2009-10-10 01:58:23 +05:30
Lp = Lp + gdot_twin ( j ) * lattice_Stwin ( : , : , index_myFamily + i , myStructure )
!* Calculation of the tangent of Lp
forall ( k = 1 : 3 , l = 1 : 3 , m = 1 : 3 , n = 1 : 3 ) &
2009-10-12 21:56:16 +05:30
dLp_dTstar3333 ( k , l , m , n ) = &
dLp_dTstar3333 ( k , l , m , n ) + dgdot_dtautwin ( j ) * &
lattice_Stwin ( k , l , index_myFamily + i , myStructure ) * &
lattice_Stwin ( m , n , index_myFamily + i , myStructure )
2009-10-10 01:58:23 +05:30
enddo
enddo
2009-10-12 21:56:16 +05:30
dLp_dTstar = math_Plain3333to99 ( dLp_dTstar3333 )
!if ((ip==1).and.(el==1)) then
2009-10-13 20:35:21 +05:30
! write(6,*) '#LP/TANGENT#'
! write(6,*)
! write(6,*) 'Tstar_v', Tstar_v
! write(6,*) 'tau_slip', tau_slip
! write(6,'(a10,/,4(3(e20.8,x),/))') 'state',state(1,1,1)%p
! write(6,'(a,/,3(3(f10.4,x)/))') 'Lp',Lp
! write(6,'(a,/,9(9(f10.4,x)/))') 'dLp_dTstar',dLp_dTstar
2009-10-12 21:56:16 +05:30
!endif
2009-10-10 01:58:23 +05:30
return
2009-10-06 20:49:19 +05:30
end subroutine
2009-10-10 01:58:23 +05:30
function constitutive_dislotwin_dotState ( Tstar_v , Temperature , state , g , ip , el )
2009-10-06 20:49:19 +05:30
!*********************************************************************
!* rate of change of microstructure *
!* INPUT: *
!* - Temperature : temperature *
!* - state : microstructure quantities *
!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) *
!* - ipc : component-ID at current integration point *
!* - ip : current integration point *
!* - el : current element *
!* OUTPUT: *
!* - constitutive_dotState : evolution of state variable *
!*********************************************************************
2009-10-12 21:56:16 +05:30
use prec , only : pReal , pInt , p_vec
2009-10-10 01:58:23 +05:30
use math , only : pi
2010-08-04 05:17:00 +05:30
use mesh , only : mesh_NcpElems , mesh_maxNips
use material , only : homogenization_maxNgrains , material_phase , phase_constitutionInstance
use lattice , only : lattice_Sslip , lattice_Sslip_v , &
lattice_Stwin , lattice_Stwin_v , &
lattice_maxNslipFamily , lattice_maxNtwinFamily , &
lattice_NslipSystem , lattice_NtwinSystem , &
lattice_shearTwin
2009-10-10 01:58:23 +05:30
implicit none
2009-10-06 20:49:19 +05:30
2009-10-10 01:58:23 +05:30
!* Input-Output variables
integer ( pInt ) , intent ( in ) :: g , ip , el
real ( pReal ) , intent ( in ) :: Temperature
2009-10-12 21:56:16 +05:30
real ( pReal ) , dimension ( 6 ) , intent ( in ) :: Tstar_v
2009-10-10 01:58:23 +05:30
type ( p_vec ) , dimension ( homogenization_maxNgrains , mesh_maxNips , mesh_NcpElems ) , intent ( in ) :: state
real ( pReal ) , dimension ( constitutive_dislotwin_sizeDotState ( phase_constitutionInstance ( material_phase ( g , ip , el ) ) ) ) :: &
constitutive_dislotwin_dotState
!* Local variables
2010-08-17 19:53:55 +05:30
integer ( pInt ) MyInstance , MyStructure , ns , nt , f , i , j , index_myFamily
2009-10-12 21:56:16 +05:30
real ( pReal ) sumf , StressRatio_p , StressRatio_pminus1 , BoltzmannRatio , DotGamma0 , &
EdgeDipMinDistance , AtomicVolume , VacancyDiffusion , StressRatio_r
2009-10-10 01:58:23 +05:30
real ( pReal ) , dimension ( constitutive_dislotwin_totalNslip ( phase_constitutionInstance ( material_phase ( g , ip , el ) ) ) ) :: &
2009-10-12 21:56:16 +05:30
gdot_slip , tau_slip , DotRhoMultiplication , EdgeDipDistance , DotRhoEdgeEdgeAnnihilation , DotRhoEdgeDipAnnihilation , &
2009-10-10 01:58:23 +05:30
ClimbVelocity , DotRhoEdgeDipClimb , DotRhoDipFormation
2010-08-04 05:17:00 +05:30
real ( pReal ) , dimension ( constitutive_dislotwin_totalNtwin ( phase_constitutionInstance ( material_phase ( g , ip , el ) ) ) ) :: &
2010-08-17 19:53:55 +05:30
tau_twin
2009-10-10 01:58:23 +05:30
!* Shortened notation
myInstance = phase_constitutionInstance ( material_phase ( g , ip , el ) )
2010-08-17 19:53:55 +05:30
MyStructure = constitutive_dislotwin_structure ( myInstance )
2009-10-10 01:58:23 +05:30
ns = constitutive_dislotwin_totalNslip ( myInstance )
nt = constitutive_dislotwin_totalNtwin ( myInstance )
2009-10-12 21:56:16 +05:30
!* Total twin volume fraction
sumf = sum ( state ( g , ip , el ) % p ( ( 2 * ns + 1 ) : ( 2 * ns + nt ) ) ) ! safe for nt == 0
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_dotState = 0.0_pReal
!* Dislocation density evolution
gdot_slip = 0.0_pReal
j = 0_pInt
do f = 1 , lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum ( lattice_NslipSystem ( 1 : f - 1 , MyStructure ) ) ! at which index starts my family
do i = 1 , constitutive_dislotwin_Nslip ( f , myInstance ) ! process each (active) slip system in family
j = j + 1_pInt
2009-10-12 21:56:16 +05:30
!* Resolved shear stress on slip system
2010-08-17 19:53:55 +05:30
tau_slip ( j ) = dot_product ( Tstar_v , lattice_Sslip_v ( : , index_myFamily + i , myStructure ) )
2009-10-12 21:56:16 +05:30
!* Stress ratios
StressRatio_p = ( abs ( tau_slip ( j ) ) / state ( g , ip , el ) % p ( 5 * ns + 3 * nt + j ) ) ** constitutive_dislotwin_p ( myInstance )
StressRatio_pminus1 = ( abs ( tau_slip ( j ) ) / state ( g , ip , el ) % p ( 5 * ns + 3 * nt + j ) ) ** ( constitutive_dislotwin_p ( myInstance ) - 1.0_pReal )
!* Boltzmann ratio
BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem ( f , myInstance ) / ( kB * Temperature )
!* Initial shear rates
DotGamma0 = &
state ( g , ip , el ) % p ( j ) * constitutive_dislotwin_burgersPerSlipSystem ( f , myInstance ) * &
2010-08-17 19:53:55 +05:30
constitutive_dislotwin_v0PerSlipSystem ( f , myInstance )
2009-10-12 21:56:16 +05:30
!* Shear rates due to slip
2010-08-04 05:17:00 +05:30
gdot_slip ( j ) = DotGamma0 * exp ( - BoltzmannRatio * ( 1 - StressRatio_p ) ** constitutive_dislotwin_q ( myInstance ) ) * &
sign ( 1.0_pReal , tau_slip ( j ) )
2009-10-12 21:56:16 +05:30
2009-10-10 01:58:23 +05:30
!* Multiplication
2009-10-12 21:56:16 +05:30
DotRhoMultiplication ( j ) = abs ( gdot_slip ( j ) ) / &
( constitutive_dislotwin_burgersPerSlipSystem ( f , myInstance ) * state ( g , ip , el ) % p ( 4 * ns + 2 * nt + j ) )
!* Dipole formation
2009-11-06 23:16:12 +05:30
EdgeDipMinDistance = &
constitutive_dislotwin_CEdgeDipMinDistance ( myInstance ) * constitutive_dislotwin_burgersPerSlipSystem ( f , myInstance )
2009-10-12 21:56:16 +05:30
if ( tau_slip ( j ) == 0.0_pReal ) then
2009-10-13 20:35:21 +05:30
DotRhoDipFormation ( j ) = 0.0_pReal
2010-08-17 19:53:55 +05:30
else
2009-10-12 21:56:16 +05:30
EdgeDipDistance ( j ) = &
( 3.0_pReal * constitutive_dislotwin_Gmod ( myInstance ) * constitutive_dislotwin_burgersPerSlipSystem ( f , myInstance ) ) / &
( 1 6.0_pReal * pi * abs ( tau_slip ( j ) ) )
2010-08-04 05:17:00 +05:30
if ( EdgeDipDistance ( j ) > state ( g , ip , el ) % p ( 4 * ns + 2 * nt + j ) ) EdgeDipDistance ( j ) = state ( g , ip , el ) % p ( 4 * ns + 2 * nt + j )
2010-08-17 19:53:55 +05:30
if ( EdgeDipDistance ( j ) < EdgeDipMinDistance ) EdgeDipDistance ( j ) = EdgeDipMinDistance
2009-10-12 21:56:16 +05:30
DotRhoDipFormation ( j ) = &
( ( 2.0_pReal * EdgeDipDistance ( j ) ) / constitutive_dislotwin_burgersPerSlipSystem ( f , myInstance ) ) * &
state ( g , ip , el ) % p ( j ) * abs ( gdot_slip ( j ) )
endif
2010-08-17 19:53:55 +05:30
2009-10-12 21:56:16 +05:30
!* Spontaneous annihilation of 2 single edge dislocations
DotRhoEdgeEdgeAnnihilation ( j ) = &
( ( 2.0_pReal * EdgeDipMinDistance ) / constitutive_dislotwin_burgersPerSlipSystem ( f , myInstance ) ) * &
state ( g , ip , el ) % p ( j ) * abs ( gdot_slip ( j ) )
2010-08-17 19:53:55 +05:30
2009-10-12 21:56:16 +05:30
!* Spontaneous annihilation of a single edge dislocation with a dipole constituent
DotRhoEdgeDipAnnihilation ( j ) = &
( ( 2.0_pReal * EdgeDipMinDistance ) / constitutive_dislotwin_burgersPerSlipSystem ( f , myInstance ) ) * &
state ( g , ip , el ) % p ( ns + j ) * abs ( gdot_slip ( j ) )
!* Dislocation dipole climb
AtomicVolume = &
constitutive_dislotwin_CAtomicVolume ( myInstance ) * constitutive_dislotwin_burgersPerSlipSystem ( f , myInstance ) ** ( 3.0_pReal )
VacancyDiffusion = &
constitutive_dislotwin_D0 ( myInstance ) * exp ( - constitutive_dislotwin_Qsd ( myInstance ) / ( kB * Temperature ) )
2009-10-13 20:35:21 +05:30
if ( tau_slip ( j ) == 0.0_pReal ) then
2010-08-04 05:17:00 +05:30
DotRhoEdgeDipClimb ( j ) = 0.0_pReal
else
ClimbVelocity ( j ) = &
2009-10-13 20:35:21 +05:30
( ( 3.0_pReal * constitutive_dislotwin_Gmod ( myInstance ) * VacancyDiffusion * AtomicVolume ) / ( 2.0_pReal * pi * kB * Temperature ) ) * &
( 1 / ( EdgeDipDistance ( j ) + EdgeDipMinDistance ) )
2010-08-04 05:17:00 +05:30
DotRhoEdgeDipClimb ( j ) = &
2009-11-06 23:16:12 +05:30
( 4.0_pReal * ClimbVelocity ( j ) * state ( g , ip , el ) % p ( ns + j ) ) / ( EdgeDipDistance ( j ) - EdgeDipMinDistance )
2010-08-04 05:17:00 +05:30
endif
2009-10-12 21:56:16 +05:30
!* Edge dislocation density rate of change
constitutive_dislotwin_dotState ( j ) = &
DotRhoMultiplication ( j ) - DotRhoDipFormation ( j ) - DotRhoEdgeEdgeAnnihilation ( j )
!* Edge dislocation dipole density rate of change
constitutive_dislotwin_dotState ( ns + j ) = &
DotRhoDipFormation ( j ) - DotRhoEdgeDipAnnihilation ( j ) - DotRhoEdgeDipClimb ( j )
2009-10-10 01:58:23 +05:30
enddo
enddo
!* Twin volume fraction evolution
j = 0_pInt
do f = 1 , lattice_maxNtwinFamily ! loop over all twin families
index_myFamily = sum ( lattice_NtwinSystem ( 1 : f - 1 , MyStructure ) ) ! at which index starts my family
do i = 1 , constitutive_dislotwin_Ntwin ( f , myInstance ) ! process each (active) twin system in family
j = j + 1_pInt
2009-10-12 21:56:16 +05:30
!* Resolved shear stress on twin system
tau_twin ( j ) = dot_product ( Tstar_v , lattice_Stwin_v ( : , index_myFamily + i , myStructure ) )
!* Stress ratios
StressRatio_r = ( state ( g , ip , el ) % p ( 6 * ns + 3 * nt + j ) / tau_twin ( j ) ) ** constitutive_dislotwin_r ( myInstance )
2010-08-17 19:53:55 +05:30
2009-10-12 21:56:16 +05:30
!* Shear rates and their derivatives due to twin
if ( tau_twin ( j ) > 0.0_pReal ) then
2009-10-13 20:35:21 +05:30
constitutive_dislotwin_dotState ( 2 * ns + j ) = &
2009-10-12 21:56:16 +05:30
( constitutive_dislotwin_MaxTwinFraction ( myInstance ) - sumf ) * &
2010-08-17 19:53:55 +05:30
state ( g , ip , el ) % p ( 6 * ns + 4 * nt + j ) * constitutive_dislotwin_Ndot0PerTwinSystem ( f , myInstance ) * exp ( - StressRatio_r )
2009-10-12 21:56:16 +05:30
endif
2009-10-10 01:58:23 +05:30
enddo
enddo
2009-10-12 21:56:16 +05:30
2009-10-13 20:35:21 +05:30
!write(6,*) '#DOTSTATE#'
!write(6,*)
!write(6,'(a,/,4(3(f30.20,x)/))') 'tau slip',tau_slip
!write(6,'(a,/,4(3(f30.20,x)/))') 'gamma slip',gdot_slip
!write(6,'(a,/,4(3(f30.20,x)/))') 'RhoEdge',state(g,ip,el)%p(1:ns)
!write(6,'(a,/,4(3(f30.20,x)/))') 'Threshold Slip', state(g,ip,el)%p(5*ns+3*nt+1:6*ns+3*nt)
!write(6,'(a,/,4(3(f30.20,x)/))') 'Multiplication',DotRhoMultiplication
!write(6,'(a,/,4(3(f30.20,x)/))') 'DipFormation',DotRhoDipFormation
!write(6,'(a,/,4(3(f30.20,x)/))') 'SingleSingle',DotRhoEdgeEdgeAnnihilation
!write(6,'(a,/,4(3(f30.20,x)/))') 'SingleDipole',DotRhoEdgeDipAnnihilation
2010-08-17 19:53:55 +05:30
!write(6,'(a,/,4(3(f30.20,x)/))') 'DipClimb',DotRhoEdgeDipClimb
2009-10-12 21:56:16 +05:30
2009-10-10 01:58:23 +05:30
return
2009-10-06 20:49:19 +05:30
end function
2009-10-10 01:58:23 +05:30
pure function constitutive_dislotwin_dotTemperature ( Tstar_v , Temperature , state , g , ip , el )
2009-10-06 20:49:19 +05:30
!*********************************************************************
!* rate of change of microstructure *
!* INPUT: *
!* - Temperature : temperature *
!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) *
!* - ipc : component-ID at current integration point *
!* - ip : current integration point *
!* - el : current element *
!* OUTPUT: *
!* - constitutive_dotTemperature : evolution of Temperature *
!*********************************************************************
2009-10-10 01:58:23 +05:30
use prec , only : pReal , pInt , p_vec
use mesh , only : mesh_NcpElems , mesh_maxNips
use material , only : homogenization_maxNgrains
implicit none
2009-10-06 20:49:19 +05:30
2009-10-10 01:58:23 +05:30
!* Input-Output variables
integer ( pInt ) , intent ( in ) :: g , ip , el
real ( pReal ) , intent ( in ) :: Temperature
real ( pReal ) , dimension ( 6 ) , intent ( in ) :: Tstar_v
type ( p_vec ) , dimension ( homogenization_maxNgrains , mesh_maxNips , mesh_NcpElems ) , intent ( in ) :: state
real ( pReal ) constitutive_dislotwin_dotTemperature
constitutive_dislotwin_dotTemperature = 0.0_pReal
2010-08-17 19:53:55 +05:30
2009-10-10 01:58:23 +05:30
return
2009-10-06 20:49:19 +05:30
end function
2009-10-10 01:58:23 +05:30
pure function constitutive_dislotwin_postResults ( Tstar_v , Temperature , dt , state , g , ip , el )
2009-10-06 20:49:19 +05:30
!*********************************************************************
!* return array of constitutive results *
!* INPUT: *
!* - Temperature : temperature *
!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) *
!* - dt : current time increment *
!* - ipc : component-ID at current integration point *
!* - ip : current integration point *
!* - el : current element *
!*********************************************************************
2009-10-10 01:58:23 +05:30
use prec , only : pReal , pInt , p_vec
2009-11-06 23:16:12 +05:30
use math , only : pi
2009-10-10 01:58:23 +05:30
use mesh , only : mesh_NcpElems , mesh_maxNips
use material , only : homogenization_maxNgrains , material_phase , phase_constitutionInstance , phase_Noutput
use lattice , only : lattice_Sslip_v , lattice_Stwin_v , lattice_maxNslipFamily , lattice_maxNtwinFamily , &
2010-08-17 19:53:55 +05:30
lattice_NslipSystem , lattice_NtwinSystem , lattice_shearTwin
2009-10-10 01:58:23 +05:30
implicit none
2009-10-06 20:49:19 +05:30
!* Definition of variables
2009-10-10 01:58:23 +05:30
integer ( pInt ) , intent ( in ) :: g , ip , el
real ( pReal ) , intent ( in ) :: dt , Temperature
real ( pReal ) , dimension ( 6 ) , intent ( in ) :: Tstar_v
type ( p_vec ) , dimension ( homogenization_maxNgrains , mesh_maxNips , mesh_NcpElems ) , intent ( in ) :: state
integer ( pInt ) myInstance , myStructure , ns , nt , f , o , i , c , j , index_myFamily
2009-11-06 23:16:12 +05:30
real ( pReal ) sumf , tau , StressRatio_p , StressRatio_pminus1 , BoltzmannRatio , DotGamma0 , StressRatio_r , gdot_slip , dgdot_dtauslip
2009-10-10 01:58:23 +05:30
real ( pReal ) , dimension ( constitutive_dislotwin_sizePostResults ( phase_constitutionInstance ( material_phase ( g , ip , el ) ) ) ) :: &
constitutive_dislotwin_postResults
!* Shortened notation
myInstance = phase_constitutionInstance ( material_phase ( g , ip , el ) )
2010-08-17 19:53:55 +05:30
myStructure = constitutive_dislotwin_structure ( myInstance )
2009-10-10 01:58:23 +05:30
ns = constitutive_dislotwin_totalNslip ( myInstance )
nt = constitutive_dislotwin_totalNtwin ( myInstance )
!* Total twin volume fraction
sumf = sum ( state ( g , ip , el ) % p ( ( 2 * ns + 1 ) : ( 2 * ns + nt ) ) ) ! safe for nt == 0
2010-08-17 19:53:55 +05:30
!* Required output
2009-10-10 01:58:23 +05:30
c = 0_pInt
constitutive_dislotwin_postResults = 0.0_pReal
do o = 1 , phase_Noutput ( material_phase ( g , ip , el ) )
select case ( constitutive_dislotwin_output ( o , myInstance ) )
case ( 'edge_density' )
constitutive_dislotwin_postResults ( c + 1 : c + ns ) = state ( g , ip , el ) % p ( 1 : ns )
2009-10-06 20:49:19 +05:30
c = c + ns
2009-10-12 21:56:16 +05:30
case ( 'dipole_density' )
constitutive_dislotwin_postResults ( c + 1 : c + ns ) = state ( g , ip , el ) % p ( ns + 1 : 2 * ns )
c = c + ns
2009-10-10 01:58:23 +05:30
case ( 'shear_rate_slip' )
2009-10-06 20:49:19 +05:30
j = 0_pInt
2009-10-10 01:58:23 +05:30
do f = 1 , lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum ( lattice_NslipSystem ( 1 : f - 1 , myStructure ) ) ! at which index starts my family
do i = 1 , constitutive_dislotwin_Nslip ( f , myInstance ) ! process each (active) slip system in family
2009-10-06 20:49:19 +05:30
j = j + 1_pInt
2009-10-12 21:56:16 +05:30
!* Resolved shear stress on slip system
2010-08-17 19:53:55 +05:30
tau = dot_product ( Tstar_v , lattice_Sslip_v ( : , index_myFamily + i , myStructure ) )
2009-10-12 21:56:16 +05:30
!* Stress ratios
StressRatio_p = ( abs ( tau ) / state ( g , ip , el ) % p ( 5 * ns + 3 * nt + j ) ) ** constitutive_dislotwin_p ( myInstance )
StressRatio_pminus1 = ( abs ( tau ) / state ( g , ip , el ) % p ( 5 * ns + 3 * nt + j ) ) ** ( constitutive_dislotwin_p ( myInstance ) - 1.0_pReal )
!* Boltzmann ratio
BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem ( f , myInstance ) / ( kB * Temperature )
!* Initial shear rates
DotGamma0 = &
state ( g , ip , el ) % p ( j ) * constitutive_dislotwin_burgersPerSlipSystem ( f , myInstance ) * &
constitutive_dislotwin_v0PerSlipSystem ( f , myInstance )
2010-08-17 19:53:55 +05:30
2009-10-12 21:56:16 +05:30
!* Shear rates due to slip
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_postResults ( c + j ) = &
2009-10-12 21:56:16 +05:30
DotGamma0 * exp ( - BoltzmannRatio * ( 1 - StressRatio_p ) ** constitutive_dislotwin_q ( myInstance ) ) * sign ( 1.0_pReal , tau )
enddo ; enddo
2009-10-06 20:49:19 +05:30
c = c + ns
case ( 'mfp_slip' )
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_postResults ( c + 1 : c + ns ) = state ( g , ip , el ) % p ( ( 4 * ns + 2 * nt + 1 ) : ( 5 * ns + 2 * nt ) )
2009-10-06 20:49:19 +05:30
c = c + ns
2009-10-10 01:58:23 +05:30
case ( 'resolved_stress_slip' )
2009-10-06 20:49:19 +05:30
j = 0_pInt
2009-10-10 01:58:23 +05:30
do f = 1 , lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum ( lattice_NslipSystem ( 1 : f - 1 , myStructure ) ) ! at which index starts my family
do i = 1 , constitutive_dislotwin_Nslip ( f , myInstance ) ! process each (active) slip system in family
2009-10-06 20:49:19 +05:30
j = j + 1_pInt
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_postResults ( c + j ) = dot_product ( Tstar_v , lattice_Sslip_v ( : , index_myFamily + i , myStructure ) )
2009-10-12 21:56:16 +05:30
enddo ; enddo
2009-10-06 20:49:19 +05:30
c = c + ns
2009-10-10 01:58:23 +05:30
case ( 'threshold_stress_slip' )
constitutive_dislotwin_postResults ( c + 1 : c + ns ) = state ( g , ip , el ) % p ( ( 5 * ns + 3 * nt + 1 ) : ( 6 * ns + 3 * nt ) )
2009-11-06 23:16:12 +05:30
c = c + ns
2010-08-04 05:17:00 +05:30
case ( 'edge_dipole_distance' )
2009-11-06 23:16:12 +05:30
j = 0_pInt
do f = 1 , lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum ( lattice_NslipSystem ( 1 : f - 1 , myStructure ) ) ! at which index starts my family
do i = 1 , constitutive_dislotwin_Nslip ( f , myInstance ) ! process each (active) slip system in family
j = j + 1_pInt
2010-08-17 19:53:55 +05:30
constitutive_dislotwin_postResults ( c + j ) = &
2010-08-04 05:17:00 +05:30
( 3.0_pReal * constitutive_dislotwin_Gmod ( myInstance ) * constitutive_dislotwin_burgersPerSlipSystem ( f , myInstance ) ) / &
( 1 6.0_pReal * pi * abs ( dot_product ( Tstar_v , lattice_Sslip_v ( : , index_myFamily + i , myStructure ) ) ) )
constitutive_dislotwin_postResults ( c + j ) = min ( constitutive_dislotwin_postResults ( c + j ) , state ( g , ip , el ) % p ( 4 * ns + 2 * nt + j ) )
! constitutive_dislotwin_postResults(c+j) = max(constitutive_dislotwin_postResults(c+j),state(g,ip,el)%p(4*ns+2*nt+j))
2009-11-06 23:16:12 +05:30
enddo ; enddo
2009-10-06 20:49:19 +05:30
c = c + ns
2009-10-10 01:58:23 +05:30
case ( 'twin_fraction' )
constitutive_dislotwin_postResults ( c + 1 : c + nt ) = state ( g , ip , el ) % p ( ( 2 * ns + 1 ) : ( 2 * ns + nt ) )
2009-10-06 20:49:19 +05:30
c = c + nt
2009-10-10 01:58:23 +05:30
case ( 'shear_rate_twin' )
2010-08-17 19:53:55 +05:30
if ( nt > 0_pInt ) then
2009-10-12 21:56:16 +05:30
j = 0_pInt
do f = 1 , lattice_maxNtwinFamily ! loop over all twin families
index_myFamily = sum ( lattice_NtwinSystem ( 1 : f - 1 , myStructure ) ) ! at which index starts my family
do i = 1 , constitutive_dislotwin_Ntwin ( f , myInstance ) ! process each (active) twin system in family
2009-10-06 20:49:19 +05:30
j = j + 1_pInt
2009-10-10 01:58:23 +05:30
2009-10-12 21:56:16 +05:30
!* Resolved shear stress on twin system
tau = dot_product ( Tstar_v , lattice_Stwin_v ( : , index_myFamily + i , myStructure ) )
!* Stress ratios
StressRatio_r = ( state ( g , ip , el ) % p ( 6 * ns + 3 * nt + j ) / tau ) ** constitutive_dislotwin_r ( myInstance )
2009-10-06 20:49:19 +05:30
2009-10-12 21:56:16 +05:30
!* Shear rates and their derivatives due to twin
if ( tau > 0.0_pReal ) then
constitutive_dislotwin_postResults ( c + j ) = &
( constitutive_dislotwin_MaxTwinFraction ( myInstance ) - sumf ) * &
state ( g , ip , el ) % p ( 6 * ns + 4 * nt + j ) * constitutive_dislotwin_Ndot0PerTwinSystem ( f , myInstance ) * exp ( - StressRatio_r )
endif
enddo ; enddo
endif
c = c + nt
2009-10-06 20:49:19 +05:30
case ( 'mfp_twin' )
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_postResults ( c + 1 : c + nt ) = state ( g , ip , el ) % p ( ( 5 * ns + 2 * nt + 1 ) : ( 5 * ns + 3 * nt ) )
2009-10-06 20:49:19 +05:30
c = c + nt
2009-10-10 01:58:23 +05:30
case ( 'resolved_stress_twin' )
2009-10-12 21:56:16 +05:30
if ( nt > 0_pInt ) then
j = 0_pInt
do f = 1 , lattice_maxNtwinFamily ! loop over all slip families
2010-08-04 05:17:00 +05:30
index_myFamily = sum ( lattice_NtwinSystem ( 1 : f - 1 , myStructure ) ) ! at which index starts my family
do i = 1 , constitutive_dislotwin_Ntwin ( f , myInstance ) ! process each (active) slip system in family
2009-10-06 20:49:19 +05:30
j = j + 1_pInt
2009-10-10 01:58:23 +05:30
constitutive_dislotwin_postResults ( c + j ) = dot_product ( Tstar_v , lattice_Stwin_v ( : , index_myFamily + i , myStructure ) )
2009-10-12 21:56:16 +05:30
enddo ; enddo
endif
2009-10-06 20:49:19 +05:30
c = c + nt
2009-10-10 01:58:23 +05:30
case ( 'threshold_stress_twin' )
constitutive_dislotwin_postResults ( c + 1 : c + nt ) = state ( g , ip , el ) % p ( ( 6 * ns + 3 * nt + 1 ) : ( 6 * ns + 4 * nt ) )
2009-10-06 20:49:19 +05:30
c = c + nt
2009-11-06 23:16:12 +05:30
case ( 'stress_exponent' )
j = 0_pInt
do f = 1 , lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum ( lattice_NslipSystem ( 1 : f - 1 , myStructure ) ) ! at which index starts my family
do i = 1 , constitutive_dislotwin_Nslip ( f , myInstance ) ! process each (active) slip system in family
j = j + 1_pInt
!* Resolved shear stress on slip system
2010-08-17 19:53:55 +05:30
tau = dot_product ( Tstar_v , lattice_Sslip_v ( : , index_myFamily + i , myStructure ) )
2009-11-06 23:16:12 +05:30
!* Stress ratios
StressRatio_p = ( abs ( tau ) / state ( g , ip , el ) % p ( 5 * ns + 3 * nt + j ) ) ** constitutive_dislotwin_p ( myInstance )
StressRatio_pminus1 = ( abs ( tau ) / state ( g , ip , el ) % p ( 5 * ns + 3 * nt + j ) ) ** ( constitutive_dislotwin_p ( myInstance ) - 1.0_pReal )
!* Boltzmann ratio
BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem ( f , myInstance ) / ( kB * Temperature )
!* Initial shear rates
DotGamma0 = &
state ( g , ip , el ) % p ( j ) * constitutive_dislotwin_burgersPerSlipSystem ( f , myInstance ) * &
constitutive_dislotwin_v0PerSlipSystem ( f , myInstance )
2010-08-17 19:53:55 +05:30
2009-11-06 23:16:12 +05:30
!* Shear rates due to slip
gdot_slip = &
DotGamma0 * exp ( - BoltzmannRatio * ( 1 - StressRatio_p ) ** constitutive_dislotwin_q ( myInstance ) ) * sign ( 1.0_pReal , tau )
!* Derivatives of shear rates
dgdot_dtauslip = &
( ( abs ( gdot_slip ) * BoltzmannRatio * &
constitutive_dislotwin_p ( myInstance ) * constitutive_dislotwin_q ( myInstance ) ) / state ( g , ip , el ) % p ( 5 * ns + 3 * nt + j ) ) * &
StressRatio_pminus1 * ( 1 - StressRatio_p ) ** ( constitutive_dislotwin_q ( myInstance ) - 1.0_pReal )
!* Stress exponent
2010-08-04 05:17:00 +05:30
if ( gdot_slip == 0.0_pReal ) then
constitutive_dislotwin_postResults ( c + j ) = 0.0_pReal
else
constitutive_dislotwin_postResults ( c + j ) = ( tau / gdot_slip ) * dgdot_dtauslip
endif
2009-11-06 23:16:12 +05:30
enddo ; enddo
c = c + ns
2009-10-06 20:49:19 +05:30
end select
2009-10-10 01:58:23 +05:30
enddo
return
2009-10-06 20:49:19 +05:30
end function
2009-10-12 21:56:16 +05:30
END MODULE