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

added comments (doxygen conform) to phenopowerlaw, added warning on specifying h0_sliptwin as it has no effect

This commit is contained in:
Martin Diehl 2012-10-11 14:49:12 +00:00
parent c783aa1265
commit 324dfda5a2
3 changed files with 255 additions and 327 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
code/IO.f90
View File

@ -1518,6 +1518,8 @@ subroutine IO_warning(warning_ID,e,i,g,ext_msg)
msg = 'Found Spectral solver parameter ' msg = 'Found Spectral solver parameter '
case (41_pInt) case (41_pInt)
msg = 'Found PETSc solver parameter ' msg = 'Found PETSc solver parameter '
case (42_pInt)
msg = 'parameter has no effect '
case (47_pInt) case (47_pInt)
msg = 'No valid parameter for FFTW given, using FFTW_PATIENT' msg = 'No valid parameter for FFTW given, using FFTW_PATIENT'
case (101_pInt) case (101_pInt)

2
code/config/material.config
View File

@ -164,7 +164,7 @@ twin_c 0
twin_d 0 twin_d 0
twin_e 0 twin_e 0
h0_slipslip 75e6 h0_slipslip 75e6
h0_sliptwin 0 #h0_sliptwin 0 no effect
h0_twinslip 0 h0_twinslip 0
h0_twintwin 0 h0_twintwin 0
interaction_slipslip 1 1 1.4 1.4 1.4 1.4 interaction_slipslip 1 1 1.4 1.4 1.4 1.4

342
code/constitutive_phenopowerlaw.f90
View File

@ -16,60 +16,13 @@
! You should have received a copy of the GNU General Public License ! You should have received a copy of the GNU General Public License
! along with DAMASK. If not, see <http://www.gnu.org/licenses/>. ! along with DAMASK. If not, see <http://www.gnu.org/licenses/>.
! !
!############################################################## !--------------------------------------------------------------------------------------------------
!* $Id$ ! $Id$
!***************************************************** !--------------------------------------------------------------------------------------------------
!* Module: CONSTITUTIVE_PHENOPOWERLAW * !> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!***************************************************** !> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!* contains: * !> @brief phenomenological crystal plasticity formulation using a powerlaw fitting
!* - constitutive equations * !--------------------------------------------------------------------------------------------------
!* - parameters definition *
!*****************************************************
![Alu]
!plasticity phenopowerlaw
!(output) resistance_slip
!(output) shearrate_slip
!(output) resolvedstress_slip
!(output) totalshear
!(output) resistance_twin
!(output) shearrate_twin
!(output) resolvedstress_twin
!(output) totalvolfrac
!lattice_structure hex
!covera_ratio 1.587
!Nslip 3 3 6 12 # per family
!Ntwin 6 6 6 6 # per family
!
!c11 162.2e9
!c12 91.8e9
!c13 68.8e9
!c33 180.5e9
!c44 46.7e9
!
!gdot0_slip 0.001
!n_slip 50
!tau0_slip 65e6 22e6 52e6 50e6 # per family
!tausat_slip 80e6 180e6 140e6 140e6 # per family
!a_slip 1
!gdot0_twin 0.001
!n_twin 50
!tau0_twin 52e6 52e6 52e6 52e6 # per family
!s_pr 50e6 # push-up stress for slip saturation due to twinning
!twin_b 2
!twin_C 25
!twin_d 0.1
!twin_e 0.1
!h0_slipslip 10e6
!h0_sliptwin 0
!h0_twinslip 625e6
!h0_twintwin 400e6
!interaction_slipslip 5.5 5.5 1.0 52.0 5.5 5.5 1.0 52.0 27.5 0.2 72.8 1.0 72.8 72.8 27.5 1.1 1.4 5.5 7.7 7.7
!interaction_sliptwin 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
!interaction_twinslip 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
!interaction_twintwin 1 1 1 1 1 1 1 1 10 10 10 10 10 10 10 10 10 10 10 10
!relevantResistance 1
module constitutive_phenopowerlaw module constitutive_phenopowerlaw
use prec, only: pReal,pInt use prec, only: pReal,pInt
@ -82,62 +35,63 @@ module constitutive_phenopowerlaw
integer(pInt), dimension(:), allocatable, public :: & integer(pInt), dimension(:), allocatable, public :: &
constitutive_phenopowerlaw_sizeDotState, & constitutive_phenopowerlaw_sizeDotState, &
constitutive_phenopowerlaw_sizeState, & constitutive_phenopowerlaw_sizeState, &
constitutive_phenopowerlaw_sizePostResults, & ! cumulative size of post results constitutive_phenopowerlaw_sizePostResults, & !< cumulative size of post results
constitutive_phenopowerlaw_structure constitutive_phenopowerlaw_structure
integer(pInt), dimension(:), allocatable, private :: & integer(pInt), dimension(:), allocatable, private :: &
constitutive_phenopowerlaw_Noutput, & ! number of outputs per instance of this constitution constitutive_phenopowerlaw_Noutput, & !< number of outputs per instance of this constitution
constitutive_phenopowerlaw_totalNslip, & ! no. of slip system used in simulation constitutive_phenopowerlaw_totalNslip, & !< no. of slip system used in simulation
constitutive_phenopowerlaw_totalNtwin ! no. of twin system used in simulation constitutive_phenopowerlaw_totalNtwin !< no. of twin system used in simulation
integer(pInt), dimension(:,:), allocatable, target, public :: & integer(pInt), dimension(:,:), allocatable, target, public :: &
constitutive_phenopowerlaw_sizePostResult ! size of each post result output constitutive_phenopowerlaw_sizePostResult !< size of each post result output
integer(pInt), dimension(:,:), allocatable, private :: & integer(pInt), dimension(:,:), allocatable, private :: &
constitutive_phenopowerlaw_Nslip, & ! active number of slip systems per family constitutive_phenopowerlaw_Nslip, & !< active number of slip systems per family (input parameter, per family)
constitutive_phenopowerlaw_Ntwin ! active number of twin systems per family constitutive_phenopowerlaw_Ntwin !< active number of twin systems per family (input parameter, per family)
character(len=64), dimension(:,:), allocatable, target, public :: & character(len=64), dimension(:,:), allocatable, target, public :: &
constitutive_phenopowerlaw_output ! name of each post result output constitutive_phenopowerlaw_output !< name of each post result output
character(len=32), dimension(:), allocatable, private :: & character(len=32), dimension(:), allocatable, private :: &
constitutive_phenopowerlaw_structureName constitutive_phenopowerlaw_structureName
real(pReal), dimension(:), allocatable, private :: & real(pReal), dimension(:), allocatable, private :: &
constitutive_phenopowerlaw_CoverA, & constitutive_phenopowerlaw_CoverA, & !< c/a of the crystal (input parameter)
constitutive_phenopowerlaw_C11, & constitutive_phenopowerlaw_C11, & !< component 11 of the stiffness matrix (input parameter)
constitutive_phenopowerlaw_C12, & constitutive_phenopowerlaw_C12, & !< component 12 of the stiffness matrix (input parameter)
constitutive_phenopowerlaw_C13, & constitutive_phenopowerlaw_C13, & !< component 13 of the stiffness matrix (input parameter)
constitutive_phenopowerlaw_C33, & constitutive_phenopowerlaw_C33, & !< component 33 of the stiffness matrix (input parameter)
constitutive_phenopowerlaw_C44, & constitutive_phenopowerlaw_C44, & !< component 44 of the stiffness matrix (input parameter)
constitutive_phenopowerlaw_gdot0_slip, & constitutive_phenopowerlaw_gdot0_slip, & !< reference shear strain rate for slip (input parameter)
constitutive_phenopowerlaw_n_slip, & constitutive_phenopowerlaw_gdot0_twin, & !< reference shear strain rate for twin (input parameter)
constitutive_phenopowerlaw_n_twin, & constitutive_phenopowerlaw_n_slip, & !< (inverse?) of the stress exponent for slip (input parameter)
constitutive_phenopowerlaw_gdot0_twin constitutive_phenopowerlaw_n_twin !< (inverse?) of the stress exponent for twin (input parameter)
real(pReal), dimension(:,:), allocatable, private :: & real(pReal), dimension(:,:), allocatable, private :: &
constitutive_phenopowerlaw_tau0_slip, & constitutive_phenopowerlaw_tau0_slip, & !< initial critical shear stress for slip (input parameter, per family)
constitutive_phenopowerlaw_tausat_slip, & constitutive_phenopowerlaw_tau0_twin, & !< initial critical shear stress for twin (input parameter, per family)
constitutive_phenopowerlaw_tau0_twin constitutive_phenopowerlaw_tausat_slip !< maximum critical shear stress for slip (input parameter, per family)
real(pReal), dimension(:), allocatable, private :: & real(pReal), dimension(:), allocatable, private :: &
constitutive_phenopowerlaw_spr, & constitutive_phenopowerlaw_spr, & !< push-up factor for slip saturation due to twinning
constitutive_phenopowerlaw_twinB, & constitutive_phenopowerlaw_twinB, &
constitutive_phenopowerlaw_twinC, & constitutive_phenopowerlaw_twinC, &
constitutive_phenopowerlaw_twinD, & constitutive_phenopowerlaw_twinD, &
constitutive_phenopowerlaw_twinE, & constitutive_phenopowerlaw_twinE, &
constitutive_phenopowerlaw_h0_slipslip, & constitutive_phenopowerlaw_h0_slipslip, & !< reference hardening slip - slip (input parameter)
constitutive_phenopowerlaw_h0_sliptwin, & constitutive_phenopowerlaw_h0_sliptwin, & !< reference hardening slip - twin (input parameter, no effect at the moment)
constitutive_phenopowerlaw_h0_twinslip, & constitutive_phenopowerlaw_h0_twinslip, & !< reference hardening twin - slip (input parameter)
constitutive_phenopowerlaw_h0_twintwin, & constitutive_phenopowerlaw_h0_twintwin, & !< reference hardening twin - twin (input parameter)
constitutive_phenopowerlaw_a_slip, & constitutive_phenopowerlaw_a_slip, &
constitutive_phenopowerlaw_aTolResistance constitutive_phenopowerlaw_aTolResistance
real(pReal), dimension(:,:), allocatable, private :: & real(pReal), dimension(:,:), allocatable, private :: &
constitutive_phenopowerlaw_interaction_slipslip, & constitutive_phenopowerlaw_interaction_slipslip, & !< interaction factors slip - slip (input parameter)
constitutive_phenopowerlaw_interaction_sliptwin, & constitutive_phenopowerlaw_interaction_sliptwin, & !< interaction factors slip - twin (input parameter)
constitutive_phenopowerlaw_interaction_twinslip, & constitutive_phenopowerlaw_interaction_twinslip, & !< interaction factors twin - slip (input parameter)
constitutive_phenopowerlaw_interaction_twintwin constitutive_phenopowerlaw_interaction_twintwin !< interaction factors twin - twin (input parameter)
real(pReal), dimension(:,:,:), allocatable, private :: & real(pReal), dimension(:,:,:), allocatable, private :: &
constitutive_phenopowerlaw_hardeningMatrix_slipslip, & constitutive_phenopowerlaw_hardeningMatrix_slipslip, &
@ -160,10 +114,10 @@ module constitutive_phenopowerlaw
contains contains
!--------------------------------------------------------------------------------------------------
!> @brief reading in parameters from material config and doing consistency checks
!--------------------------------------------------------------------------------------------------
subroutine constitutive_phenopowerlaw_init(myFile) subroutine constitutive_phenopowerlaw_init(myFile)
!**************************************
!* Module initialization *
!**************************************
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment) use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use math, only: math_Mandel3333to66, & use math, only: math_Mandel3333to66, &
math_Voigt66to3333 math_Voigt66to3333
@ -312,7 +266,8 @@ subroutine constitutive_phenopowerlaw_init(myFile)
cycle cycle
case ('(output)') case ('(output)')
constitutive_phenopowerlaw_Noutput(i) = constitutive_phenopowerlaw_Noutput(i) + 1_pInt constitutive_phenopowerlaw_Noutput(i) = constitutive_phenopowerlaw_Noutput(i) + 1_pInt
constitutive_phenopowerlaw_output(constitutive_phenopowerlaw_Noutput(i),i) = IO_lc(IO_stringValue(line,positions,2_pInt)) constitutive_phenopowerlaw_output(constitutive_phenopowerlaw_Noutput(i),i) = &
IO_lc(IO_stringValue(line,positions,2_pInt))
case ('lattice_structure') case ('lattice_structure')
constitutive_phenopowerlaw_structureName(i) = IO_lc(IO_stringValue(line,positions,2_pInt)) constitutive_phenopowerlaw_structureName(i) = IO_lc(IO_stringValue(line,positions,2_pInt))
case ('covera_ratio') case ('covera_ratio')
@ -366,6 +321,7 @@ subroutine constitutive_phenopowerlaw_init(myFile)
constitutive_phenopowerlaw_h0_slipslip(i) = IO_floatValue(line,positions,2_pInt) constitutive_phenopowerlaw_h0_slipslip(i) = IO_floatValue(line,positions,2_pInt)
case ('h0_sliptwin') case ('h0_sliptwin')
constitutive_phenopowerlaw_h0_sliptwin(i) = IO_floatValue(line,positions,2_pInt) constitutive_phenopowerlaw_h0_sliptwin(i) = IO_floatValue(line,positions,2_pInt)
call IO_warning(42_pInt,ext_msg=trim(tag)//' ('//constitutive_phenopowerlaw_label//')')
case ('h0_twinslip') case ('h0_twinslip')
constitutive_phenopowerlaw_h0_twinslip(i) = IO_floatValue(line,positions,2_pInt) constitutive_phenopowerlaw_h0_twinslip(i) = IO_floatValue(line,positions,2_pInt)
case ('h0_twintwin') case ('h0_twintwin')
@ -509,7 +465,9 @@ subroutine constitutive_phenopowerlaw_init(myFile)
constitutive_phenopowerlaw_Cslip_66(:,:,i) = & constitutive_phenopowerlaw_Cslip_66(:,:,i) = &
math_Mandel3333to66(math_Voigt66to3333(constitutive_phenopowerlaw_Cslip_66(:,:,i))) math_Mandel3333to66(math_Voigt66to3333(constitutive_phenopowerlaw_Cslip_66(:,:,i)))
do f = 1_pInt,lattice_maxNslipFamily ! >>> interaction slip -- X !--------------------------------------------------------------------------------------------------
! interaction slip -- X
do f = 1_pInt,lattice_maxNslipFamily
index_myFamily = sum(constitutive_phenopowerlaw_Nslip(1:f-1_pInt,i)) index_myFamily = sum(constitutive_phenopowerlaw_Nslip(1:f-1_pInt,i))
do j = 1_pInt,constitutive_phenopowerlaw_Nslip(f,i) ! loop over (active) systems in my family (slip) do j = 1_pInt,constitutive_phenopowerlaw_Nslip(f,i) ! loop over (active) systems in my family (slip)
do o = 1_pInt,lattice_maxNslipFamily do o = 1_pInt,lattice_maxNslipFamily
@ -534,7 +492,9 @@ subroutine constitutive_phenopowerlaw_init(myFile)
enddo; enddo enddo; enddo
do f = 1_pInt,lattice_maxNtwinFamily ! >>> interaction twin -- X !--------------------------------------------------------------------------------------------------
! interaction twin -- X
do f = 1_pInt,lattice_maxNtwinFamily
index_myFamily = sum(constitutive_phenopowerlaw_Ntwin(1:f-1_pInt,i)) index_myFamily = sum(constitutive_phenopowerlaw_Ntwin(1:f-1_pInt,i))
do j = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,i) ! loop over (active) systems in my family (twin) do j = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,i) ! loop over (active) systems in my family (twin)
@ -561,18 +521,15 @@ subroutine constitutive_phenopowerlaw_init(myFile)
enddo; enddo enddo; enddo
! report to out file... ! report to out file...
enddo enddo
return
end subroutine constitutive_phenopowerlaw_init end subroutine constitutive_phenopowerlaw_init
!--------------------------------------------------------------------------------------------------
!> @brief initial microstructural state
!--------------------------------------------------------------------------------------------------
function constitutive_phenopowerlaw_stateInit(myInstance) function constitutive_phenopowerlaw_stateInit(myInstance)
!*********************************************************************
!* initial microstructural state *
!*********************************************************************
use lattice, only: lattice_maxNslipFamily, lattice_maxNtwinFamily use lattice, only: lattice_maxNslipFamily, lattice_maxNtwinFamily
implicit none implicit none
@ -596,74 +553,64 @@ function constitutive_phenopowerlaw_stateInit(myInstance)
sum(constitutive_phenopowerlaw_Ntwin(1:i ,myInstance))) = & sum(constitutive_phenopowerlaw_Ntwin(1:i ,myInstance))) = &
constitutive_phenopowerlaw_tau0_twin(i,myInstance) constitutive_phenopowerlaw_tau0_twin(i,myInstance)
enddo enddo
return
end function constitutive_phenopowerlaw_stateInit end function constitutive_phenopowerlaw_stateInit
!********************************************************************* !--------------------------------------------------------------------------------------------------
!* absolute state tolerance * !> @brief absolute state tolerance
!********************************************************************* !--------------------------------------------------------------------------------------------------
pure function constitutive_phenopowerlaw_aTolState(myInstance) pure function constitutive_phenopowerlaw_aTolState(myInstance)
implicit none
implicit none integer(pInt), intent(in) :: myInstance ! number specifying the current instance of the plasticity
!*** input variables real(pReal), dimension(constitutive_phenopowerlaw_sizeState(myInstance)) :: &
integer(pInt), intent(in) :: myInstance ! number specifying the current instance of the plasticity
!*** output variables
real(pReal), dimension(constitutive_phenopowerlaw_sizeState(myInstance)) :: &
constitutive_phenopowerlaw_aTolState ! relevant state values for the current instance of this plasticity constitutive_phenopowerlaw_aTolState ! relevant state values for the current instance of this plasticity
!*** local variables
constitutive_phenopowerlaw_aTolState = constitutive_phenopowerlaw_aTolResistance(myInstance) constitutive_phenopowerlaw_aTolState = constitutive_phenopowerlaw_aTolResistance(myInstance)
end function constitutive_phenopowerlaw_aTolState end function constitutive_phenopowerlaw_aTolState
!--------------------------------------------------------------------------------------------------
!> @brief homogenized elacticity matrix
!--------------------------------------------------------------------------------------------------
function constitutive_phenopowerlaw_homogenizedC(state,ipc,ip,el) function constitutive_phenopowerlaw_homogenizedC(state,ipc,ip,el)
!*********************************************************************
!* homogenized elacticity matrix *
!* INPUT: *
!* - state : state variables *
!* - ipc : component-ID of current integration point *
!* - ip : current integration point *
!* - el : current element *
!*********************************************************************
use prec, only: p_vec use prec, only: p_vec
use mesh, only: mesh_NcpElems,mesh_maxNips use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase, phase_plasticityInstance use material, only: homogenization_maxNgrains,material_phase, phase_plasticityInstance
implicit none implicit none
integer(pInt), intent(in) :: ipc,ip,el integer(pInt), intent(in) :: &
ipc, & !component-ID of current integration point
ip, & !current integration point
el !current element
integer(pInt) matID integer(pInt) matID
real(pReal), dimension(6,6) :: constitutive_phenopowerlaw_homogenizedC real(pReal), dimension(6,6) :: constitutive_phenopowerlaw_homogenizedC
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
state ! state variables
matID = phase_plasticityInstance(material_phase(ipc,ip,el)) matID = phase_plasticityInstance(material_phase(ipc,ip,el))
constitutive_phenopowerlaw_homogenizedC = constitutive_phenopowerlaw_Cslip_66(:,:,matID) constitutive_phenopowerlaw_homogenizedC = constitutive_phenopowerlaw_Cslip_66(:,:,matID)
return
end function constitutive_phenopowerlaw_homogenizedC end function constitutive_phenopowerlaw_homogenizedC
!--------------------------------------------------------------------------------------------------
!> @brief calculate derived quantities from state (dummy subroutine, not used here)
!--------------------------------------------------------------------------------------------------
subroutine constitutive_phenopowerlaw_microstructure(Temperature,state,ipc,ip,el) subroutine constitutive_phenopowerlaw_microstructure(Temperature,state,ipc,ip,el)
!*********************************************************************
!* calculate derived quantities from state (not used here) *
!* INPUT: *
!* - Tp : temperature *
!* - ipc : component-ID of current integration point *
!* - ip : current integration point *
!* - el : current element *
!*********************************************************************
use prec, only: pReal,pInt,p_vec use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_NcpElems,mesh_maxNips use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase, phase_plasticityInstance use material, only: homogenization_maxNgrains,material_phase, phase_plasticityInstance
implicit none implicit none
integer(pInt) ipc,ip,el, matID integer(pInt), intent(in) :: &
real(pReal) Temperature ipc, & !component-ID of current integration point
ip, & !current integration point
el !current element
integer(pInt) :: matID
real(pReal), intent(in) :: Temperature ! temperature
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state
matID = phase_plasticityInstance(material_phase(ipc,ip,el)) matID = phase_plasticityInstance(material_phase(ipc,ip,el))
@ -671,18 +618,10 @@ subroutine constitutive_phenopowerlaw_microstructure(Temperature,state,ipc,ip,el
end subroutine constitutive_phenopowerlaw_microstructure end subroutine constitutive_phenopowerlaw_microstructure
!--------------------------------------------------------------------------------------------------
!> @brief plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,state,ipc,ip,el) subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,state,ipc,ip,el)
!*********************************************************************
!* plastic velocity gradient and its tangent *
!* INPUT: *
!* - 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) *
!*********************************************************************
use prec, only: p_vec use prec, only: p_vec
use math, only: math_Plain3333to99 use math, only: math_Plain3333to99
use lattice, only: lattice_Sslip,lattice_Sslip_v,lattice_Stwin,lattice_Stwin_v, lattice_maxNslipFamily, lattice_maxNtwinFamily, & use lattice, only: lattice_Sslip,lattice_Sslip_v,lattice_Stwin,lattice_Stwin_v, lattice_maxNslipFamily, lattice_maxNtwinFamily, &
@ -691,14 +630,17 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temp
use material, only: homogenization_maxNgrains,material_phase, phase_plasticityInstance use material, only: homogenization_maxNgrains,material_phase, phase_plasticityInstance
implicit none implicit none
integer(pInt) ipc,ip,el integer(pInt), intent(in) :: &
ipc, & ! component-ID at current integration point
ip, & ! current integration point
el ! current element
integer(pInt) matID,nSlip,nTwin,f,i,j,k,l,m,n, structID,index_Gamma,index_F,index_myFamily integer(pInt) matID,nSlip,nTwin,f,i,j,k,l,m,n, structID,index_Gamma,index_F,index_myFamily
real(pReal) Temperature real(pReal) Temperature
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state
real(pReal), dimension(6) :: Tstar_v real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola Kirchhoff stress tensor (Mandel)
real(pReal), dimension(3,3) :: Lp real(pReal), dimension(3,3), intent(out) :: Lp ! plastic velocity gradient
real(pReal), dimension(3,3,3,3) :: dLp_dTstar3333 real(pReal), dimension(3,3,3,3):: dLp_dTstar3333 ! derivative of Lp (4th-rank tensor)
real(pReal), dimension(9,9) :: dLp_dTstar real(pReal), dimension(9,9), intent(out) :: dLp_dTstar
real(pReal), dimension(constitutive_phenopowerlaw_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & real(pReal), dimension(constitutive_phenopowerlaw_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
gdot_slip,dgdot_dtauslip,tau_slip gdot_slip,dgdot_dtauslip,tau_slip
real(pReal), dimension(constitutive_phenopowerlaw_totalNtwin(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & real(pReal), dimension(constitutive_phenopowerlaw_totalNtwin(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
@ -723,16 +665,16 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temp
do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family
j = j+1_pInt j = j+1_pInt
!* Calculation of Lp !--------------------------------------------------------------------------------------------------
! Calculation of Lp
tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(1:6,index_myFamily+i,structID)) tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(1:6,index_myFamily+i,structID))
gdot_slip(j) = constitutive_phenopowerlaw_gdot0_slip(matID)*(abs(tau_slip(j))/state(ipc,ip,el)%p(j))**& gdot_slip(j) = constitutive_phenopowerlaw_gdot0_slip(matID)*(abs(tau_slip(j))/state(ipc,ip,el)%p(j))**&
constitutive_phenopowerlaw_n_slip(matID)*sign(1.0_pReal,tau_slip(j)) constitutive_phenopowerlaw_n_slip(matID)*sign(1.0_pReal,tau_slip(j))
Lp = Lp + (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F Lp = Lp + (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F
gdot_slip(j)*lattice_Sslip(1:3,1:3,index_myFamily+i,structID) gdot_slip(j)*lattice_Sslip(1:3,1:3,index_myFamily+i,structID)
!* Calculation of the tangent of Lp !--------------------------------------------------------------------------------------------------
! Calculation of the tangent of Lp
if (gdot_slip(j) /= 0.0_pReal) then if (gdot_slip(j) /= 0.0_pReal) then
dgdot_dtauslip(j) = gdot_slip(j)*constitutive_phenopowerlaw_n_slip(matID)/tau_slip(j) dgdot_dtauslip(j) = gdot_slip(j)*constitutive_phenopowerlaw_n_slip(matID)/tau_slip(j)
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
@ -749,8 +691,8 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temp
do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family
j = j+1_pInt j = j+1_pInt
!* Calculation of Lp !--------------------------------------------------------------------------------------------------
! Calculation of Lp
tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,structID)) tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,structID))
gdot_twin(j) = (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F gdot_twin(j) = (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F
constitutive_phenopowerlaw_gdot0_twin(matID)*& constitutive_phenopowerlaw_gdot0_twin(matID)*&
@ -758,8 +700,8 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temp
constitutive_phenopowerlaw_n_twin(matID)*max(0.0_pReal,sign(1.0_pReal,tau_twin(j))) constitutive_phenopowerlaw_n_twin(matID)*max(0.0_pReal,sign(1.0_pReal,tau_twin(j)))
Lp = Lp + gdot_twin(j)*lattice_Stwin(1:3,1:3,index_myFamily+i,structID) Lp = Lp + gdot_twin(j)*lattice_Stwin(1:3,1:3,index_myFamily+i,structID)
!* Calculation of the tangent of Lp !--------------------------------------------------------------------------------------------------
! Calculation of the tangent of Lp
if (gdot_twin(j) /= 0.0_pReal) then if (gdot_twin(j) /= 0.0_pReal) then
dgdot_dtautwin(j) = gdot_twin(j)*constitutive_phenopowerlaw_n_twin(matID)/tau_twin(j) dgdot_dtautwin(j) = gdot_twin(j)*constitutive_phenopowerlaw_n_twin(matID)/tau_twin(j)
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
@ -772,21 +714,13 @@ subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temp
dLp_dTstar = math_Plain3333to99(dLp_dTstar3333) dLp_dTstar = math_Plain3333to99(dLp_dTstar3333)
return
end subroutine constitutive_phenopowerlaw_LpAndItsTangent end subroutine constitutive_phenopowerlaw_LpAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief of change of microstructure, evolution of state variable
!--------------------------------------------------------------------------------------------------
function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el) function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el)
!*********************************************************************
!* rate of change of microstructure *
!* INPUT: *
!* - 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 *
!*********************************************************************
use prec, only: p_vec use prec, only: p_vec
use lattice, only: lattice_Sslip_v, lattice_Stwin_v, lattice_maxNslipFamily, lattice_maxNtwinFamily, & use lattice, only: lattice_Sslip_v, lattice_Stwin_v, lattice_maxNslipFamily, lattice_maxNtwinFamily, &
lattice_NslipSystem,lattice_NtwinSystem,lattice_shearTwin lattice_NslipSystem,lattice_NtwinSystem,lattice_shearTwin
@ -794,11 +728,14 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el
use material, only: homogenization_maxNgrains,material_phase, phase_plasticityInstance use material, only: homogenization_maxNgrains,material_phase, phase_plasticityInstance
implicit none implicit none
integer(pInt) ipc,ip,el integer(pInt), intent(in) :: &
ipc, & !< component-ID at current integration point
ip, & !< current integration point
el !< current element
integer(pInt) matID,nSlip,nTwin,f,i,j, structID,index_Gamma,index_F,index_myFamily integer(pInt) matID,nSlip,nTwin,f,i,j, structID,index_Gamma,index_F,index_myFamily
real(pReal) Temperature,c_slipslip,c_sliptwin,c_twinslip,c_twintwin, ssat_offset real(pReal) Temperature,c_slipslip,c_sliptwin,c_twinslip,c_twintwin, ssat_offset
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state
real(pReal), dimension(6) :: Tstar_v real(pReal), dimension(6), intent(in) :: Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel)
real(pReal), dimension(constitutive_phenopowerlaw_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & real(pReal), dimension(constitutive_phenopowerlaw_totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
gdot_slip,tau_slip,h_slipslip,h_sliptwin gdot_slip,tau_slip,h_slipslip,h_sliptwin
real(pReal), dimension(constitutive_phenopowerlaw_totalNtwin(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & real(pReal), dimension(constitutive_phenopowerlaw_totalNtwin(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
@ -817,8 +754,8 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el
constitutive_phenopowerlaw_dotState = 0.0_pReal constitutive_phenopowerlaw_dotState = 0.0_pReal
!-- system-independent (nonlinear) prefactors to M_xx matrices !--------------------------------------------------------------------------------------------------
! system-independent (nonlinear) prefactors to M_xx matrices
c_slipslip = constitutive_phenopowerlaw_h0_slipslip(matID)*& c_slipslip = constitutive_phenopowerlaw_h0_slipslip(matID)*&
(1.0_pReal + & (1.0_pReal + &
constitutive_phenopowerlaw_twinC(matID)*state(ipc,ip,el)%p(index_F)**constitutive_phenopowerlaw_twinB(matID)) constitutive_phenopowerlaw_twinC(matID)*state(ipc,ip,el)%p(index_F)**constitutive_phenopowerlaw_twinB(matID))
@ -828,8 +765,8 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el
c_twintwin = constitutive_phenopowerlaw_h0_twintwin(matID)*& c_twintwin = constitutive_phenopowerlaw_h0_twintwin(matID)*&
state(ipc,ip,el)%p(index_F)**constitutive_phenopowerlaw_twinD(matID) state(ipc,ip,el)%p(index_F)**constitutive_phenopowerlaw_twinD(matID)
!-- add system-dependent part and calculate dot gammas !--------------------------------------------------------------------------------------------------
! add system-dependent part and calculate dot gammas
ssat_offset = constitutive_phenopowerlaw_spr(matID)*sqrt(state(ipc,ip,el)%p(index_F)) ssat_offset = constitutive_phenopowerlaw_spr(matID)*sqrt(state(ipc,ip,el)%p(index_F))
j = 0_pInt j = 0_pInt
do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families
@ -842,8 +779,8 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el
h_sliptwin(j) = c_sliptwin ! no system-dependent part h_sliptwin(j) = c_sliptwin ! no system-dependent part
!* Calculation of dot gamma !--------------------------------------------------------------------------------------------------
! Calculation of dot gamma
tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(1:6,index_myFamily+i,structID)) tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(1:6,index_myFamily+i,structID))
gdot_slip(j) = constitutive_phenopowerlaw_gdot0_slip(matID)*(abs(tau_slip(j))/state(ipc,ip,el)%p(j))**& gdot_slip(j) = constitutive_phenopowerlaw_gdot0_slip(matID)*(abs(tau_slip(j))/state(ipc,ip,el)%p(j))**&
constitutive_phenopowerlaw_n_slip(matID)*sign(1.0_pReal,tau_slip(j)) constitutive_phenopowerlaw_n_slip(matID)*sign(1.0_pReal,tau_slip(j))
@ -858,8 +795,8 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el
h_twinslip(j) = c_twinslip ! no system-dependent parts h_twinslip(j) = c_twinslip ! no system-dependent parts
h_twintwin(j) = c_twintwin h_twintwin(j) = c_twintwin
!* Calculation of dot vol frac !--------------------------------------------------------------------------------------------------
! Calculation of dot vol frac
tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,structID)) tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,structID))
gdot_twin(j) = (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F gdot_twin(j) = (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F
constitutive_phenopowerlaw_gdot0_twin(matID)*& constitutive_phenopowerlaw_gdot0_twin(matID)*&
@ -868,8 +805,8 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el
enddo enddo
enddo enddo
!-- calculate the overall hardening based on above !--------------------------------------------------------------------------------------------------
! calculate the overall hardening based on above
j = 0_pInt j = 0_pInt
do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families do f = 1_pInt,lattice_maxNslipFamily ! loop over all slip families
do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family
@ -898,9 +835,9 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el
end function constitutive_phenopowerlaw_dotState end function constitutive_phenopowerlaw_dotState
!********************************************************************* !--------------------------------------------------------------------------------------------------
!* (instantaneous) incremental change of microstructure * !> @brief (instantaneous) incremental change of microstructure
!********************************************************************* !--------------------------------------------------------------------------------------------------
function constitutive_phenopowerlaw_deltaState(Tstar_v, Temperature, state, g,ip,el) function constitutive_phenopowerlaw_deltaState(Tstar_v, Temperature, state, g,ip,el)
use prec, only: pReal, & use prec, only: pReal, &
@ -914,7 +851,6 @@ use material, only: homogenization_maxNgrains, &
implicit none implicit none
!*** input variables
integer(pInt), intent(in) :: g, & ! current grain number integer(pInt), intent(in) :: g, & ! current grain number
ip, & ! current integration point ip, & ! current integration point
el ! current element number el ! current element number
@ -923,30 +859,24 @@ real(pReal), dimension(6), intent(in) :: Tstar_v ! current
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
state ! current microstructural state state ! current microstructural state
!*** output variables
real(pReal), dimension(constitutive_phenopowerlaw_sizeDotState(phase_plasticityInstance(material_phase(g,ip,el)))) :: & real(pReal), dimension(constitutive_phenopowerlaw_sizeDotState(phase_plasticityInstance(material_phase(g,ip,el)))) :: &
constitutive_phenopowerlaw_deltaState ! change of state variables / microstructure constitutive_phenopowerlaw_deltaState ! change of state variables / microstructure
!*** local variables
constitutive_phenopowerlaw_deltaState = 0.0_pReal constitutive_phenopowerlaw_deltaState = 0.0_pReal
endfunction end function constitutive_phenopowerlaw_deltaState
!**************************************************************** !--------------------------------------------------------------------------------------------------
!* calculates the rate of change of temperature * !> @brief calculates the rate of change of temperature (dummy function)
!**************************************************************** !--------------------------------------------------------------------------------------------------
pure function constitutive_phenopowerlaw_dotTemperature(Tstar_v,Temperature,state,ipc,ip,el) pure function constitutive_phenopowerlaw_dotTemperature(Tstar_v,Temperature,state,ipc,ip,el)
!*** variables and functions from other modules ***!
use prec, only: pReal,pInt,p_vec use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_NcpElems, mesh_maxNips use mesh, only: mesh_NcpElems, mesh_maxNips
use material, only: homogenization_maxNgrains use material, only: homogenization_maxNgrains
implicit none implicit none
!*** input variables ***!
real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola Kirchhoff stress tensor in Mandel notation real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola Kirchhoff stress tensor in Mandel notation
real(pReal), intent(in) :: Temperature real(pReal), intent(in) :: Temperature
integer(pInt), intent(in):: ipc, & ! grain number integer(pInt), intent(in):: ipc, & ! grain number
@ -954,26 +884,17 @@ pure function constitutive_phenopowerlaw_dotTemperature(Tstar_v,Temperature,stat
el ! element number el ! element number
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state ! state of the current microstructure type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state ! state of the current microstructure
!*** output variables ***!
real(pReal) constitutive_phenopowerlaw_dotTemperature ! rate of change of temparature real(pReal) constitutive_phenopowerlaw_dotTemperature ! rate of change of temparature
! calculate dotTemperature
constitutive_phenopowerlaw_dotTemperature = 0.0_pReal constitutive_phenopowerlaw_dotTemperature = 0.0_pReal
end function constitutive_phenopowerlaw_dotTemperature end function constitutive_phenopowerlaw_dotTemperature
!--------------------------------------------------------------------------------------------------
!> @brief return array of constitutive results
!--------------------------------------------------------------------------------------------------
pure function constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,state,ipc,ip,el) pure function constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,state,ipc,ip,el)
!*********************************************************************
!* return array of constitutive results *
!* INPUT: *
!* - 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 *
!*********************************************************************
use prec, only: pReal,pInt,p_vec use prec, only: pReal,pInt,p_vec
use lattice, only: lattice_Sslip_v,lattice_Stwin_v, lattice_maxNslipFamily, lattice_maxNtwinFamily, & use lattice, only: lattice_Sslip_v,lattice_Stwin_v, lattice_maxNslipFamily, lattice_maxNtwinFamily, &
lattice_NslipSystem,lattice_NtwinSystem lattice_NslipSystem,lattice_NtwinSystem
@ -981,9 +902,14 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,stat
use material, only: homogenization_maxNgrains,material_phase,phase_plasticityInstance,phase_Noutput use material, only: homogenization_maxNgrains,material_phase,phase_plasticityInstance,phase_Noutput
implicit none implicit none
integer(pInt), intent(in) :: ipc,ip,el integer(pInt), intent(in) :: &
real(pReal), intent(in) :: dt,Temperature ipc, & !component-ID at current integration point
real(pReal), dimension(6), intent(in) :: Tstar_v ip, & !current integration point
el !current element
real(pReal), intent(in) :: &
dt, & !current time increment
Temperature
real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola Kirchhoff stress tensor (Mandel)
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state
integer(pInt) matID,o,f,i,c,nSlip,nTwin,j, structID,index_Gamma,index_F,index_myFamily integer(pInt) matID,o,f,i,c,nSlip,nTwin,j, structID,index_Gamma,index_F,index_myFamily
real(pReal) tau real(pReal) tau