From ef2c2af47138efdd3c2dec8f721888a26a9a263d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 28 Jun 2013 18:58:10 +0000 Subject: [PATCH] doxygen documentation for J2 --- code/constitutive_j2.f90 | 833 ++++++++++++++++++++------------------- 1 file changed, 417 insertions(+), 416 deletions(-) diff --git a/code/constitutive_j2.f90 b/code/constitutive_j2.f90 index f4c47d90f..923d9fe0e 100644 --- a/code/constitutive_j2.f90 +++ b/code/constitutive_j2.f90 @@ -16,100 +16,90 @@ ! You should have received a copy of the GNU General Public License ! along with DAMASK. If not, see . ! -!############################################################## -!* $Id$ -!***************************************************** -!* Module: CONSTITUTIVE_J2 * -!***************************************************** -!* contains: * -!* - constitutive equations * -!* - parameters definition * -!***************************************************** - -! [Alu] -! plasticity j2 -! (output) flowstress -! (output) strainrate -! c11 110.9e9 # (3 C11 + 2 C12 + 2 C44) / 5 ... with C44 = C11-C12 !! -! c12 58.34e9 # (1 C11 + 4 C12 - 1 C44) / 5 -! taylorfactor 3 -! tau0 31e6 -! gdot0 0.001 -! n 20 -! h0 75e6 -! tausat 63e6 -! a 2.25 - +!-------------------------------------------------------------------------------------------------- +! $Id$ +!-------------------------------------------------------------------------------------------------- +!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH +!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH +!> @brief Isotropic (J2) Plasticity +!> @details Isotropic (J2) Plasticity which resembles the phenopowerlaw plasticity without +!! resolving the stress on the slip systems. Will give the response of phenopowerlaw for an +!! untextured polycrystal +!-------------------------------------------------------------------------------------------------- module constitutive_j2 - - use prec, only: pReal,pInt + use prec, only: & + pReal,& + pInt implicit none private - character (len=*), parameter, public :: constitutive_j2_LABEL = 'j2' + character (len=*), parameter, public :: & + CONSTITUTIVE_J2_label = 'j2' !< label for this constitutive model - integer(pInt), dimension(:), allocatable, public :: & + integer(pInt), dimension(:), allocatable, public :: & constitutive_j2_sizeDotState, & constitutive_j2_sizeState, & constitutive_j2_sizePostResults - integer(pInt), dimension(:,:), allocatable, target, public :: & - constitutive_j2_sizePostResult ! size of each post result output + integer(pInt), dimension(:,:), allocatable, target, public :: & + constitutive_j2_sizePostResult !< size of each post result output - character(len=64), dimension(:,:), allocatable, target, public :: & - constitutive_j2_output ! name of each post result output + character(len=64), dimension(:,:), allocatable, target, public :: & + constitutive_j2_output !< name of each post result output - integer(pInt), dimension(:), allocatable, private :: & - constitutive_j2_Noutput !< name of each post result output + integer(pInt), dimension(:), allocatable, private :: & + constitutive_j2_Noutput !< ?? - character(len=32), dimension(:), allocatable, private :: & + character(len=32), dimension(:), allocatable, private :: & constitutive_j2_structureName - real(pReal), dimension(:), allocatable, private :: & - !* Visco-plastic constitutive_j2 parameters - constitutive_j2_fTaylor, & - constitutive_j2_tau0, & - constitutive_j2_gdot0, & - constitutive_j2_n, & - !* h0 as function of h0 = A + B log (gammadot) + real(pReal), dimension(:), allocatable, private :: & + constitutive_j2_fTaylor, & !< Taylor factor + constitutive_j2_tau0, & !< initial plastic stress + constitutive_j2_gdot0, & !< reference velocity + constitutive_j2_n, & !< Visco-plastic parameter +!-------------------------------------------------------------------------------------------------- +! h0 as function of h0 = A + B log (gammadot) constitutive_j2_h0, & constitutive_j2_h0_slopeLnRate, & - constitutive_j2_tausat, & + constitutive_j2_tausat, & !< final plastic stress constitutive_j2_a, & constitutive_j2_aTolResistance, & - !* Parameters of normalized strain rate vs. stress function: - !* tausat += (asinh((gammadot / SinhFitA)**(1 / SinhFitD)))**(1 / SinhFitC) / (SinhFitB * (gammadot / gammadot0)**(1/n)) - constitutive_j2_tausat_SinhFitA, & - constitutive_j2_tausat_SinhFitB, & - constitutive_j2_tausat_SinhFitC, & - constitutive_j2_tausat_SinhFitD +!-------------------------------------------------------------------------------------------------- +! tausat += (asinh((gammadot / SinhFitA)**(1 / SinhFitD)))**(1 / SinhFitC) / (SinhFitB * (gammadot / gammadot0)**(1/n)) + constitutive_j2_tausat_SinhFitA, & !< fitting parameter for normalized strain rate vs. stress function + constitutive_j2_tausat_SinhFitB, & !< fitting parameter for normalized strain rate vs. stress function + constitutive_j2_tausat_SinhFitC, & !< fitting parameter for normalized strain rate vs. stress function + constitutive_j2_tausat_SinhFitD !< fitting parameter for normalized strain rate vs. stress function - real(pReal), dimension(:,:,:), allocatable, private :: & + + real(pReal), dimension(:,:,:), allocatable, private :: & constitutive_j2_Cslip_66 - public :: constitutive_j2_init, & - constitutive_j2_stateInit, & - constitutive_j2_aTolState, & - constitutive_j2_homogenizedC, & - constitutive_j2_microstructure, & - constitutive_j2_LpAndItsTangent, & - constitutive_j2_dotState, & - constitutive_j2_deltaState, & - constitutive_j2_dotTemperature, & - constitutive_j2_postResults + public :: & + constitutive_j2_init, & + constitutive_j2_stateInit, & + constitutive_j2_aTolState, & + constitutive_j2_homogenizedC, & + constitutive_j2_microstructure, & + constitutive_j2_LpAndItsTangent, & + constitutive_j2_dotState, & + constitutive_j2_deltaState, & + constitutive_j2_dotTemperature, & + constitutive_j2_postResults contains + +!-------------------------------------------------------------------------------------------------- +!> @brief module initialization +!-------------------------------------------------------------------------------------------------- subroutine constitutive_j2_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, & math_Voigt66to3333 use IO, only: & - IO_read, & IO_lc, & IO_getTag, & IO_isBlank, & @@ -117,25 +107,28 @@ subroutine constitutive_j2_init(myFile) IO_stringValue, & IO_floatValue, & IO_error, & - IO_timeStamp + IO_timeStamp, & + IO_read use material use debug, only: & debug_level, & debug_constitutive, & debug_levelBasic - use lattice, only: lattice_symmetrizeC66 + use lattice, only: & + lattice_symmetrizeC66 implicit none integer(pInt), intent(in) :: myFile - integer(pInt), parameter :: maxNchunks = 7_pInt + integer(pInt), parameter :: MAXNCHUNKS = 7_pInt - integer(pInt), dimension(1_pInt+2_pInt*maxNchunks) :: positions + integer(pInt), dimension(1_pInt+2_pInt*MAXNCHUNKS) :: positions integer(pInt) :: section = 0_pInt, maxNinstance, i,o, mySize - character(len=65536) :: tag - character(len=65536) :: line = '' ! to start initialized + character(len=65536) :: & + tag = '', & + line = '' ! to start initialized - write(6,'(/,a)') ' <<<+- constitutive_'//trim(constitutive_j2_LABEL)//' init -+>>>' + write(6,'(/,a)') ' <<<+- constitutive_'//trim(CONSTITUTIVE_J2_label)//' init -+>>>' write(6,'(a)') ' $Id$' write(6,'(a15,a)') ' Current time: ',IO_timeStamp() #include "compilation_info.f90" @@ -144,8 +137,7 @@ subroutine constitutive_j2_init(myFile) if (maxNinstance == 0_pInt) return if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) then - write(6,'(a16,1x,i5)') '# instances:',maxNinstance - write(6,*) + write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance endif allocate(constitutive_j2_sizeDotState(maxNinstance)) @@ -193,29 +185,30 @@ subroutine constitutive_j2_init(myFile) rewind(myFile) - do while (trim(line) /= '#EOF#' .and. IO_lc(IO_getTag(line,'<','>')) /= 'phase') ! wind forward to + do while (trim(line) /= '#EOF#' .and. IO_lc(IO_getTag(line,'<','>')) /= 'phase') ! wind forward to line = IO_read(myFile) enddo - do while (trim(line) /= '#EOF#') ! read thru sections of phase part + do while (trim(line) /= '#EOF#') ! read through sections of phase part line = IO_read(myFile) - 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_pInt ! advance section counter + 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_pInt ! advance section counter cycle endif - if (section > 0_pInt ) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran - if (phase_plasticity(section) == constitutive_j2_LABEL) then ! one of my sections - i = phase_plasticityInstance(section) ! which instance of my plasticity is present phase + if (section > 0_pInt ) then ! do not short-circuit here (.and. with next if-statement). It's not safe in Fortran + if (phase_plasticity(section) == CONSTITUTIVE_J2_label) then ! one of my sections + i = phase_plasticityInstance(section) ! which instance of my plasticity is present phase positions = IO_stringPos(line,maxNchunks) - tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key + tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key select case(tag) case ('plasticity','elasticity') cycle case ('(output)') constitutive_j2_Noutput(i) = constitutive_j2_Noutput(i) + 1_pInt - constitutive_j2_output(constitutive_j2_Noutput(i),i) = IO_lc(IO_stringValue(line,positions,2_pInt)) + constitutive_j2_output(constitutive_j2_Noutput(i),i) = & + IO_lc(IO_stringValue(line,positions,2_pInt)) case ('lattice_structure') constitutive_j2_structureName(i) = IO_lc(IO_stringValue(line,positions,2_pInt)) case ('c11') @@ -269,23 +262,23 @@ subroutine constitutive_j2_init(myFile) endif enddo - do i = 1_pInt,maxNinstance ! sanity checks - if (constitutive_j2_structureName(i) == '') call IO_error(205_pInt,e=i) - if (constitutive_j2_tau0(i) < 0.0_pReal) call IO_error(211_pInt,ext_msg='tau0 (' & - //constitutive_j2_label//')') - if (constitutive_j2_gdot0(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='gdot0 (' & - //constitutive_j2_label//')') - if (constitutive_j2_n(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='n (' & - //constitutive_j2_label//')') - if (constitutive_j2_tausat(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='tausat (' & - //constitutive_j2_label//')') - if (constitutive_j2_a(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='a (' & - //constitutive_j2_label//')') - if (constitutive_j2_fTaylor(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='taylorfactor (' & - //constitutive_j2_label//')') - if (constitutive_j2_aTolResistance(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='aTol_resistance (' & - //constitutive_j2_label//')') - enddo + sanityChecks: do i = 1_pInt,maxNinstance + if (constitutive_j2_structureName(i) == '') call IO_error(205_pInt,e=i) + if (constitutive_j2_tau0(i) < 0.0_pReal) call IO_error(211_pInt,ext_msg='tau0 (' & + //CONSTITUTIVE_J2_label//')') + if (constitutive_j2_gdot0(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='gdot0 (' & + //CONSTITUTIVE_J2_label//')') + if (constitutive_j2_n(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='n (' & + //CONSTITUTIVE_J2_label//')') + if (constitutive_j2_tausat(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='tausat (' & + //CONSTITUTIVE_J2_label//')') + if (constitutive_j2_a(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='a (' & + //CONSTITUTIVE_J2_label//')') + if (constitutive_j2_fTaylor(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='taylorfactor (' & + //CONSTITUTIVE_J2_label//')') + if (constitutive_j2_aTolResistance(i) <= 0.0_pReal) call IO_error(211_pInt,ext_msg='aTol_resistance (' & + //CONSTITUTIVE_J2_label//')') + enddo sanityChecks do i = 1_pInt,maxNinstance do o = 1_pInt,constitutive_j2_Noutput(i) @@ -295,10 +288,10 @@ subroutine constitutive_j2_init(myFile) case('strainrate') mySize = 1_pInt case default - call IO_error(212_pInt,ext_msg=constitutive_j2_output(o,i)//' ('//constitutive_j2_label//')') + call IO_error(212_pInt,ext_msg=constitutive_j2_output(o,i)//' ('//CONSTITUTIVE_J2_label//')') end select - if (mySize > 0_pInt) then ! any meaningful output found + if (mySize > 0_pInt) then ! any meaningful output found constitutive_j2_sizePostResult(o,i) = mySize constitutive_j2_sizePostResults(i) = & constitutive_j2_sizePostResults(i) + mySize @@ -308,387 +301,395 @@ subroutine constitutive_j2_init(myFile) constitutive_j2_sizeDotState(i) = 1_pInt constitutive_j2_sizeState(i) = 1_pInt - constitutive_j2_Cslip_66(:,:,i) = lattice_symmetrizeC66(constitutive_j2_structureName(i),& - constitutive_j2_Cslip_66(:,:,i)) + constitutive_j2_Cslip_66(1:6,1:6,i) = lattice_symmetrizeC66(constitutive_j2_structureName(i),& + constitutive_j2_Cslip_66(1:6,1:6,i)) constitutive_j2_Cslip_66(1:6,1:6,i) = & - math_Mandel3333to66(math_Voigt66to3333(constitutive_j2_Cslip_66(1:6,1:6,i))) + math_Mandel3333to66(math_Voigt66to3333(constitutive_j2_Cslip_66(1:6,1:6,i))) ! todo what is going on here? enddo end subroutine constitutive_j2_init -!********************************************************************* -!* initial microstructural state * -!********************************************************************* +!-------------------------------------------------------------------------------------------------- +!> @brief initial microstructural state +!> @detail initial microstructural state is set to the value specified by tau0 +!-------------------------------------------------------------------------------------------------- pure function constitutive_j2_stateInit(myInstance) implicit none - integer(pInt), intent(in) :: myInstance - real(pReal), dimension(1) :: constitutive_j2_stateInit - + real(pReal), dimension(1) :: constitutive_j2_stateInit + integer(pInt), intent(in) :: myInstance !< number specifying the instance of the plasticity + constitutive_j2_stateInit = constitutive_j2_tau0(myInstance) end function constitutive_j2_stateInit -!********************************************************************* -!* relevant microstructural state * -!********************************************************************* +!-------------------------------------------------------------------------------------------------- +!> @brief relevant state values for the current instance of this plasticity +!-------------------------------------------------------------------------------------------------- pure function constitutive_j2_aTolState(myInstance) implicit none - !*** input variables - integer(pInt), intent(in) :: myInstance ! number specifying the current instance of the plasticity + integer(pInt), intent(in) :: myInstance !< number specifying the instance of the plasticity - !*** output variables real(pReal), dimension(constitutive_j2_sizeState(myInstance)) :: & - constitutive_j2_aTolState ! relevant state values for the current instance of this plasticity + constitutive_j2_aTolState constitutive_j2_aTolState = constitutive_j2_aTolResistance(myInstance) end function constitutive_j2_aTolState +!-------------------------------------------------------------------------------------------------- +!> @brief homogenized elasticity matrix +!-------------------------------------------------------------------------------------------------- pure function constitutive_j2_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 mesh, only: mesh_NcpElems,mesh_maxNips - use material, only: homogenization_maxNgrains,material_phase, phase_plasticityInstance + use prec, only: & + p_vec + use mesh, only: & + mesh_NcpElems,mesh_maxNips + use material, only: & + homogenization_maxNgrains,& + material_phase, & + phase_plasticityInstance implicit none - integer(pInt), intent(in) :: ipc,ip,el - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state - integer(pInt) :: matID real(pReal), dimension(6,6) :: constitutive_j2_homogenizedC + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + state !< microstructure state - matID = phase_plasticityInstance(material_phase(ipc,ip,el)) - constitutive_j2_homogenizedC = constitutive_j2_Cslip_66(1:6,1:6,matID) + constitutive_j2_homogenizedC = constitutive_j2_Cslip_66(1:6,1:6,& + phase_plasticityInstance(material_phase(ipc,ip,el))) end function constitutive_j2_homogenizedC -pure subroutine constitutive_j2_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: p_vec - use mesh, only: mesh_NcpElems,mesh_maxNips - use material, only: homogenization_maxNgrains,material_phase, phase_plasticityInstance +!-------------------------------------------------------------------------------------------------- +!> @brief calculate derived quantities from state (not used here) +!-------------------------------------------------------------------------------------------------- +pure subroutine constitutive_j2_microstructure(temperature,state,ipc,ip,el) + use prec, only: & + p_vec + use mesh, only: & + mesh_NcpElems,& + mesh_maxNips + use material, only: & + homogenization_maxNgrains, & + material_phase, & + phase_plasticityInstance implicit none -!* Definition of variables - integer(pInt), intent(in) :: ipc,ip,el - real(pReal), intent(in) :: Temperature - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state - integer(pInt) :: matID - - matID = phase_plasticityInstance(material_phase(ipc,ip,el)) + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + real(pReal), intent(in) :: & + temperature !< temperature at IP + type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + state !< microstructure state end subroutine constitutive_j2_microstructure -!**************************************************************** -!* calculates plastic velocity gradient and its tangent * -!**************************************************************** -pure subroutine constitutive_j2_LpAndItsTangent(Lp, dLp_dTstar_99, Tstar_v, Temperature, state, g, ip, el) +!-------------------------------------------------------------------------------------------------- +!> @brief calculates plastic velocity gradient and its tangent +!-------------------------------------------------------------------------------------------------- +pure subroutine constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar_99,Tstar_v,& + temperature,state,ipc,ip,el) + use prec, only: & + p_vec + use math, only: & + math_mul6x6, & + math_Mandel6to33, & + math_Plain3333to99, & + math_deviatoric33, & + math_mul33xx33 + use mesh, only: & + mesh_NcpElems, & + mesh_maxNips + use material, only: & + homogenization_maxNgrains, & + material_phase, & + phase_plasticityInstance - !*** variables and functions from other modules ***! - use prec, only: p_vec - use math, only: math_mul6x6, & - math_Mandel6to33, & - math_Plain3333to99, & - math_deviatoric33, & - math_mul33xx33 - use mesh, only: mesh_NcpElems, & - mesh_maxNips - use material, only: homogenization_maxNgrains, & - material_phase, & - phase_plasticityInstance + implicit none + real(pReal), dimension(6), intent(in) :: & + Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation + real(pReal), intent(in) :: & + temperature !< temperature at IP + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + state !< microstructure state - implicit none - !*** input variables ***! - real(pReal), dimension(6), intent(in):: Tstar_v ! 2nd Piola Kirchhoff stress tensor in Mandel notation - real(pReal), intent(in):: Temperature - integer(pInt), intent(in):: g, & ! grain number - ip, & ! integration point number - el ! element number - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in):: state ! state of the current microstructure + real(pReal), dimension(3,3), intent(out) :: & + Lp !< plastic velocity gradient + real(pReal), dimension(9,9), intent(out) :: & + dLp_dTstar_99 !< derivative of Lp with respect to 2nd Piola Kirchhoff stress - !*** output variables ***! - real(pReal), dimension(3,3), intent(out) :: Lp ! plastic velocity gradient - real(pReal), dimension(9,9), intent(out) :: dLp_dTstar_99 ! derivative of Lp with respect to Tstar (9x9 matrix) + real(pReal), dimension(3,3) :: & + Tstar_dev_33 !< deviatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor + real(pReal), dimension(3,3,3,3) :: & + dLp_dTstar_3333 !< derivative of Lp with respect to Tstar as 4th order tensor + real(pReal) :: & + gamma_dot, & !< strainrate + norm_Tstar_dev, & !< euclidean norm of Tstar_dev + squarenorm_Tstar_dev !< square of the euclidean norm of Tstar_dev + integer(pInt) :: & + matID, & + k, l, m, n - !*** local variables ***! - real(pReal), dimension(3,3) :: Tstar_dev_33 ! deviatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor - real(pReal), dimension(3,3,3,3) :: dLp_dTstar_3333 ! derivative of Lp with respect to Tstar as 4th order tensor - real(pReal) gamma_dot, & ! strainrate - norm_Tstar_dev, & ! euclidean norm of Tstar_dev - squarenorm_Tstar_dev ! square of the euclidean norm of Tstar_dev - integer(pInt) matID, & - k, & - l, & - m, & - n - - matID = phase_plasticityInstance(material_phase(g,ip,el)) + matID = phase_plasticityInstance(material_phase(ipc,ip,el)) + Tstar_dev_33 = math_deviatoric33(math_Mandel6to33(Tstar_v)) ! deviatoric part of 2nd Piola-Kirchhoff stress + squarenorm_Tstar_dev = math_mul33xx33(Tstar_dev_33,Tstar_dev_33) + norm_Tstar_dev = sqrt(squarenorm_Tstar_dev) - ! deviatoric part of 2nd Piola-Kirchhoff stress - Tstar_dev_33 = math_deviatoric33(math_Mandel6to33(Tstar_v)) + if (norm_Tstar_dev <= 0.0_pReal) then ! Tstar == 0 --> both Lp and dLp_dTstar are zero + Lp = 0.0_pReal + dLp_dTstar_99 = 0.0_pReal + else + gamma_dot = constitutive_j2_gdot0(matID) * ( sqrt(1.5_pReal) * norm_Tstar_dev & + / &!---------------------------------------------------------------------------------- + (constitutive_j2_fTaylor(matID) * state(ipc,ip,el)%p(1)) ) **constitutive_j2_n(matID) - squarenorm_Tstar_dev = math_mul33xx33(Tstar_dev_33,Tstar_dev_33) - norm_Tstar_dev = sqrt(squarenorm_Tstar_dev) + Lp = Tstar_dev_33/norm_Tstar_dev * gamma_dot/constitutive_j2_fTaylor(matID) - ! Initialization of Lp and dLp_dTstar - Lp = 0.0_pReal - dLp_dTstar_99 = 0.0_pReal - - ! Tstar == 0 --> both Lp and dLp_dTstar are zero - if (norm_Tstar_dev > 0.0_pReal) then - - ! Calculation of gamma_dot - gamma_dot = constitutive_j2_gdot0(matID) * ( sqrt(1.5_pReal) * norm_Tstar_dev & - / &!--------------------------------------------------- - (constitutive_j2_fTaylor(matID) * state(g,ip,el)%p(1)) ) **constitutive_j2_n(matID) - - ! Calculation of Lp - Lp = Tstar_dev_33/norm_Tstar_dev * gamma_dot/constitutive_j2_fTaylor(matID) - - !* Calculation of the tangent of Lp - forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & - dLp_dTstar_3333(k,l,m,n) = (constitutive_j2_n(matID)-1.0_pReal) * Tstar_dev_33(k,l)*Tstar_dev_33(m,n) / squarenorm_Tstar_dev - forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) & - dLp_dTstar_3333(k,l,k,l) = dLp_dTstar_3333(k,l,k,l) + 1.0_pReal - dLp_dTstar_99 = math_Plain3333to99(gamma_dot / constitutive_j2_fTaylor(matID) * dLp_dTstar_3333 / norm_Tstar_dev) - end if +!-------------------------------------------------------------------------------------------------- +! Calculation of the tangent of Lp + forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) & + dLp_dTstar_3333(k,l,m,n) = (constitutive_j2_n(matID)-1.0_pReal) * & + Tstar_dev_33(k,l)*Tstar_dev_33(m,n) / squarenorm_Tstar_dev + forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) & + dLp_dTstar_3333(k,l,k,l) = dLp_dTstar_3333(k,l,k,l) + 1.0_pReal + dLp_dTstar_99 = math_Plain3333to99(gamma_dot / constitutive_j2_fTaylor(matID) * & + dLp_dTstar_3333 / norm_Tstar_dev) + end if end subroutine constitutive_j2_LpAndItsTangent -!**************************************************************** -!* calculates the rate of change of microstructure * -!**************************************************************** -pure function constitutive_j2_dotState(Tstar_v, Temperature, state, g, ip, el) +!-------------------------------------------------------------------------------------------------- +!> @brief calculates the rate of change of microstructure +!-------------------------------------------------------------------------------------------------- +pure function constitutive_j2_dotState(Tstar_v,Temperature,state,ipc,ip, el) + use prec, only: & + p_vec + use math, only: & + math_mul6x6 + use mesh, only: & + mesh_NcpElems, & + mesh_maxNips + use material, only: & + homogenization_maxNgrains, & + material_phase, & + phase_plasticityInstance + + implicit none + real(pReal), dimension(1) :: & + constitutive_j2_dotState + real(pReal), dimension(6), intent(in):: & + Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation + real(pReal), intent(in) :: & + Temperature !< temperature at integration point + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + state !< microstructure state - use prec, only: & - p_vec - use math, only: & - math_mul6x6 - use mesh, only: & - mesh_NcpElems, & - mesh_maxNips - use material, only: & - homogenization_maxNgrains, & - material_phase, & - phase_plasticityInstance - - implicit none - !*** input variables ***! - real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola Kirchhoff stress tensor in Mandel notation - real(pReal), intent(in) :: Temperature - integer(pInt), intent(in):: g, & ! grain number - ip, & ! integration point number - el ! element number - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state ! state of the current microstructure - - !*** output variables ***! - real(pReal), dimension(1) :: constitutive_j2_dotState ! evolution of state variable - - !*** local variables ***! - real(pReal), dimension(6) :: Tstar_dev_v ! deviatoric part of the 2nd Piola Kirchhoff stress tensor in Mandel notation - real(pReal) gamma_dot, & ! strainrate - hardening, & ! hardening coefficient - saturation, & ! saturation resistance - norm_Tstar_dev ! euclidean norm of Tstar_dev - integer(pInt) matID + real(pReal), dimension(6) :: & + Tstar_dev_v !< deviatoric part of the 2nd Piola Kirchhoff stress tensor in Mandel notation + real(pReal) :: & + gamma_dot, & !< strainrate + hardening, & !< hardening coefficient + saturation, & !< saturation resistance + norm_Tstar_dev !< euclidean norm of Tstar_dev + integer(pInt) :: & + matID - matID = phase_plasticityInstance(material_phase(g,ip,el)) + matID = phase_plasticityInstance(material_phase(ipc,ip,el)) +!-------------------------------------------------------------------------------------------------- +! norm of deviatoric part of 2nd Piola-Kirchhoff stress + Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal + Tstar_dev_v(4:6) = Tstar_v(4:6) + norm_Tstar_dev = sqrt(math_mul6x6(Tstar_dev_v,Tstar_dev_v)) - ! deviatoric part of 2nd Piola-Kirchhoff stress - Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal - Tstar_dev_v(4:6) = Tstar_v(4:6) - - norm_Tstar_dev = sqrt(math_mul6x6(Tstar_dev_v,Tstar_dev_v)) - - ! gamma_dot - gamma_dot = constitutive_j2_gdot0(matID) * ( sqrt(1.5_pReal) * norm_Tstar_dev & - / &!--------------------------------------------------- - (constitutive_j2_fTaylor(matID) * state(g,ip,el)%p(1)) ) ** constitutive_j2_n(matID) - - ! hardening coefficient - if (abs(gamma_dot) > 1e-12_pReal) then - if (constitutive_j2_tausat_SinhFitA(matID) == 0.0_pReal) then - saturation = constitutive_j2_tausat(matID) - else - saturation = ( constitutive_j2_tausat(matID) & - + ( log( ( gamma_dot / constitutive_j2_tausat_SinhFitA(matID)& - )**(1.0_pReal / constitutive_j2_tausat_SinhFitD(matID))& - + sqrt( ( gamma_dot / constitutive_j2_tausat_SinhFitA(matID) & - )**(2.0_pReal / constitutive_j2_tausat_SinhFitD(matID)) & - + 1.0_pReal ) & - ) & ! asinh(K) = ln(K + sqrt(K^2 +1)) - )**(1.0_pReal / constitutive_j2_tausat_SinhFitC(matID)) & - / ( constitutive_j2_tausat_SinhFitB(matID) & - * (gamma_dot / constitutive_j2_gdot0(matID))**(1.0_pReal / constitutive_j2_n(matID)) & - ) & - ) - endif - hardening = ( constitutive_j2_h0(matID) + constitutive_j2_h0_slopeLnRate(matID) * log(gamma_dot) ) & - * abs( 1.0_pReal - state(g,ip,el)%p(1)/saturation )**constitutive_j2_a(matID) & - * sign(1.0_pReal, 1.0_pReal - state(g,ip,el)%p(1)/saturation) - else - hardening = 0.0_pReal - endif - - ! dotState - constitutive_j2_dotState = hardening * gamma_dot +!-------------------------------------------------------------------------------------------------- +! strain rate + gamma_dot = constitutive_j2_gdot0(matID) * ( sqrt(1.5_pReal) * norm_Tstar_dev & + / &!----------------------------------------------------------------------------------- + (constitutive_j2_fTaylor(matID) * state(ipc,ip,el)%p(1)) ) ** constitutive_j2_n(matID) + +!-------------------------------------------------------------------------------------------------- +! hardening coefficient + if (abs(gamma_dot) > 1e-12_pReal) then + if (constitutive_j2_tausat_SinhFitA(matID) == 0.0_pReal) then + saturation = constitutive_j2_tausat(matID) + else + saturation = ( constitutive_j2_tausat(matID) & + + ( log( ( gamma_dot / constitutive_j2_tausat_SinhFitA(matID)& + )**(1.0_pReal / constitutive_j2_tausat_SinhFitD(matID))& + + sqrt( ( gamma_dot / constitutive_j2_tausat_SinhFitA(matID) & + )**(2.0_pReal / constitutive_j2_tausat_SinhFitD(matID)) & + + 1.0_pReal ) & + ) & ! asinh(K) = ln(K + sqrt(K^2 +1)) + )**(1.0_pReal / constitutive_j2_tausat_SinhFitC(matID)) & + / ( constitutive_j2_tausat_SinhFitB(matID) & + * (gamma_dot / constitutive_j2_gdot0(matID))**(1.0_pReal / constitutive_j2_n(matID)) & + ) & + ) + endif + hardening = ( constitutive_j2_h0(matID) + constitutive_j2_h0_slopeLnRate(matID) * log(gamma_dot) ) & + * abs( 1.0_pReal - state(ipc,ip,el)%p(1)/saturation )**constitutive_j2_a(matID) & + * sign(1.0_pReal, 1.0_pReal - state(ipc,ip,el)%p(1)/saturation) + else + hardening = 0.0_pReal + endif + + constitutive_j2_dotState = hardening * gamma_dot end function constitutive_j2_dotState - -!********************************************************************* -!* (instantaneous) incremental change of microstructure * -!********************************************************************* -function constitutive_j2_deltaState(Tstar_v, Temperature, state, g,ip,el) - -use prec, only: pReal, & - pInt, & - p_vec -use mesh, only: mesh_NcpElems, & - mesh_maxNips -use material, only: homogenization_maxNgrains, & - material_phase, & - phase_plasticityInstance - -implicit none - -!*** input variables -integer(pInt), intent(in) :: g, & ! current grain number - ip, & ! current integration point - el ! current element number -real(pReal), intent(in) :: Temperature ! temperature -real(pReal), dimension(6), intent(in) :: Tstar_v ! current 2nd Piola-Kirchhoff stress in Mandel notation -type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & - state ! current microstructural state - -!*** output variables -real(pReal), dimension(constitutive_j2_sizeDotState(phase_plasticityInstance(material_phase(g,ip,el)))) :: & - constitutive_j2_deltaState ! change of state variables / microstructure +!-------------------------------------------------------------------------------------------------- +!> @brief (instantaneous) incremental change of microstructure (dummy function) +!-------------------------------------------------------------------------------------------------- +pure function constitutive_j2_deltaState(Tstar_v,temperature,state,ipc,ip,el) + use prec, only: & + p_vec + use math, only: & + math_mul6x6 + use mesh, only: & + mesh_NcpElems, & + mesh_maxNips + use material, only: & + homogenization_maxNgrains, & + material_phase, & + phase_plasticityInstance -!*** local variables + implicit none + real(pReal), dimension(6), intent(in):: & + Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation + real(pReal), intent(in) :: & + Temperature !< temperature at integration point + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + state !< microstructure state + + real(pReal), dimension(constitutive_j2_sizeDotState(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & + constitutive_j2_deltaState + + constitutive_j2_deltaState = 0.0_pReal + +end function constitutive_j2_deltaState -constitutive_j2_deltaState = 0.0_pReal +!-------------------------------------------------------------------------------------------------- +!> @brief calculates the rate of change of temperature (dummy function) +!-------------------------------------------------------------------------------------------------- +real(pReal) pure function constitutive_j2_dotTemperature(Tstar_v,temperature,state,ipc,ip,el) + use prec, only: & + p_vec + use mesh, only: & + mesh_NcpElems, & + mesh_maxNips + use material, only: & + homogenization_maxNgrains + + implicit none + real(pReal), dimension(6), intent(in) :: & + Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation + real(pReal), intent(in) :: & + temperature !< temperature at integration point + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + state !< microstructure state -endfunction - - -!**************************************************************** -!* calculates the rate of change of temperature * -!**************************************************************** -pure function constitutive_j2_dotTemperature(Tstar_v, Temperature, state, g, ip, el) - - !*** variables and functions from other modules ***! - use prec, only: p_vec - use mesh, only: mesh_NcpElems,mesh_maxNips - use material, only: homogenization_maxNgrains - - implicit none - !*** input variables ***! - real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola Kirchhoff stress tensor in Mandel notation - real(pReal), intent(in) :: Temperature - integer(pInt), intent(in):: g, & ! grain number - ip, & ! integration point number - el ! element number - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state ! state of the current microstructure - - !*** output variables ***! - real(pReal) constitutive_j2_dotTemperature ! rate of change of temperature - - ! calculate dotTemperature - constitutive_j2_dotTemperature = 0.0_pReal + constitutive_j2_dotTemperature = 0.0_pReal end function constitutive_j2_dotTemperature -!********************************************************************* -!* return array of constitutive results * -!********************************************************************* -pure function constitutive_j2_postResults(Tstar_v, Temperature, dt, state, g, ip, el) +!-------------------------------------------------------------------------------------------------- +!> @brief return array of constitutive results +!-------------------------------------------------------------------------------------------------- +pure function constitutive_j2_postResults(Tstar_v,temperature,dt,state,ipc,ip,el) + use prec, only: & + p_vec + use math, only: & + math_mul6x6 + use mesh, only: & + mesh_NcpElems, & + mesh_maxNips + use material, only: & + homogenization_maxNgrains, & + material_phase, & + phase_plasticityInstance, & + phase_Noutput -!*** variables and functions from other modules ***! - use prec, only: p_vec - use math, only: math_mul6x6 - use mesh, only: mesh_NcpElems, & - mesh_maxNips - use material, only: homogenization_maxNgrains, & - material_phase, & - phase_plasticityInstance, & - phase_Noutput + implicit none + real(pReal), dimension(6), intent(in) :: & + Tstar_v !< 2nd Piola Kirchhoff stress tensor in Mandel notation + real(pReal), intent(in) :: & + temperature, & !< temperature at integration point + dt + integer(pInt), intent(in) :: & + ipc, & !< component-ID of integration point + ip, & !< integration point + el !< element + type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & + state + real(pReal), dimension(constitutive_j2_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & + constitutive_j2_postResults + + real(pReal), dimension(6) :: & + Tstar_dev_v ! deviatoric part of the 2nd Piola Kirchhoff stress tensor in Mandel notation + real(pReal) :: & + norm_Tstar_dev ! euclidean norm of Tstar_dev + integer(pInt) :: & + matID, & + o, & + c + + matID = phase_plasticityInstance(material_phase(ipc,ip,el)) + +!-------------------------------------------------------------------------------------------------- +! calculate deviatoric part of 2nd Piola-Kirchhoff stress and its norm + Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal + Tstar_dev_v(4:6) = Tstar_v(4:6) + norm_Tstar_dev = sqrt(math_mul6x6(Tstar_dev_v,Tstar_dev_v)) + + c = 0_pInt + constitutive_j2_postResults = 0.0_pReal - implicit none - !*** input variables ***! - real(pReal), dimension(6), intent(in):: Tstar_v ! 2nd Piola Kirchhoff stress tensor in Mandel notation - real(pReal), intent(in):: Temperature, & - dt ! current time increment - integer(pInt), intent(in):: g, & ! grain number - ip, & ! integration point number - el ! element number - type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state ! state of the current microstructure - - !*** output variables ***! - real(pReal), dimension(constitutive_j2_sizePostResults(phase_plasticityInstance(material_phase(g,ip,el)))) :: & - constitutive_j2_postResults - - !*** local variables ***! - real(pReal), dimension(6) :: Tstar_dev_v ! deviatoric part of the 2nd Piola Kirchhoff stress tensor in Mandel notation - real(pReal) norm_Tstar_dev ! euclidean norm of Tstar_dev - integer(pInt) matID, & - o, & - c - - !*** global variables ***! - ! constitutive_j2_gdot0 - ! constitutive_j2_fTaylor - ! constitutive_j2_n - - - matID = phase_plasticityInstance(material_phase(g,ip,el)) - - ! calculate deviatoric part of 2nd Piola-Kirchhoff stress and its norm - Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal - Tstar_dev_v(4:6) = Tstar_v(4:6) - norm_Tstar_dev = sqrt(math_mul6x6(Tstar_dev_v,Tstar_dev_v)) - - c = 0_pInt - constitutive_j2_postResults = 0.0_pReal - - do o = 1_pInt,phase_Noutput(material_phase(g,ip,el)) - select case(constitutive_j2_output(o,matID)) - case ('flowstress') - constitutive_j2_postResults(c+1_pInt) = state(g,ip,el)%p(1) - c = c + 1_pInt - case ('strainrate') - constitutive_j2_postResults(c+1_pInt) = & - constitutive_j2_gdot0(matID) * ( sqrt(1.5_pReal) * norm_Tstar_dev & - / &!--------------------------------------------------- - (constitutive_j2_fTaylor(matID) * state(g,ip,el)%p(1)) ) ** constitutive_j2_n(matID) - c = c + 1_pInt - end select - enddo + outputsLoop: do o = 1_pInt,phase_Noutput(material_phase(ipc,ip,el)) + select case(constitutive_j2_output(o,matID)) + case ('flowstress') + constitutive_j2_postResults(c+1_pInt) = state(ipc,ip,el)%p(1) + c = c + 1_pInt + case ('strainrate') + constitutive_j2_postResults(c+1_pInt) = & + constitutive_j2_gdot0(matID) * ( sqrt(1.5_pReal) * norm_Tstar_dev & + / &!---------------------------------------------------------------------------------- + (constitutive_j2_fTaylor(matID) * state(ipc,ip,el)%p(1)) ) ** constitutive_j2_n(matID) + c = c + 1_pInt + end select + enddo outputsLoop end function constitutive_j2_postResults