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

* now able to choose method for state integration (integrator and integratorStiffness in numerics.config); standard is "Fixed Point Iteration", which is basically the same as before; others available are "Explicit Euler", "AdaptiveEuler", "Classical Runge-Kutta" and "Runge-Kutta Cash-Karp" 

* now remembering stiffness similar to how we do it for Lp etc.; avoids undefined stiffness values for nonconverged stiffness calculation

* non-local stuff:
   * changed non-local kinetics (Gilman2002)
   * enforce zero shearrate for overall carrrier density below relevant density
   * enforce zero density for those states that become negative and were below relevant density before
   * dislocation velocity is not limited by V^(1/3) / dt anymore
This commit is contained in:
Christoph Kords 2010-10-01 12:18:49 +00:00
parent 9927cd7adb
commit fce7590c17
12 changed files with 1899 additions and 731 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
code/CPFEM.f90
View File

@ -118,6 +118,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
crystallite_Fp, & crystallite_Fp, &
crystallite_Lp0, & crystallite_Lp0, &
crystallite_Lp, & crystallite_Lp, &
crystallite_dPdF0, &
crystallite_dPdF, &
crystallite_Tstar0_v, & crystallite_Tstar0_v, &
crystallite_Tstar_v crystallite_Tstar_v
use homogenization, only: homogenization_init, & use homogenization, only: homogenization_init, &
@ -242,6 +244,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
crystallite_F0 = crystallite_partionedF ! crystallite deformation (_subF is perturbed...) crystallite_F0 = crystallite_partionedF ! crystallite deformation (_subF is perturbed...)
crystallite_Fp0 = crystallite_Fp ! crystallite plastic deformation crystallite_Fp0 = crystallite_Fp ! crystallite plastic deformation
crystallite_Lp0 = crystallite_Lp ! crystallite plastic velocity crystallite_Lp0 = crystallite_Lp ! crystallite plastic velocity
crystallite_dPdF0 = crystallite_dPdF ! crystallite stiffness
crystallite_Tstar0_v = crystallite_Tstar_v ! crystallite 2nd Piola Kirchhoff stress crystallite_Tstar0_v = crystallite_Tstar_v ! crystallite 2nd Piola Kirchhoff stress
forall ( i = 1:homogenization_maxNgrains, & forall ( i = 1:homogenization_maxNgrains, &
j = 1:mesh_maxNips, & j = 1:mesh_maxNips, &

2
code/IO.f90
View File

@ -1189,6 +1189,8 @@ endfunction
case (294) case (294)
msg = 'Non-positive tolerance for deformation gradient' msg = 'Non-positive tolerance for deformation gradient'
case (298)
msg = 'Chosen integration method does not exist'
case (299) case (299)
msg = 'Chosen perturbation method does not exist' msg = 'Chosen perturbation method does not exist'

111
code/constitutive.f90
View File

@ -19,10 +19,12 @@ MODULE constitutive
constitutive_state, & ! pointer array to current microstructure (end of converged time step) constitutive_state, & ! pointer array to current microstructure (end of converged time step)
constitutive_state_backup, & ! pointer array to backed up microstructure (end of converged time step) constitutive_state_backup, & ! pointer array to backed up microstructure (end of converged time step)
constitutive_dotState, & ! pointer array to evolution of current microstructure constitutive_dotState, & ! pointer array to evolution of current microstructure
constitutive_dotState_backup, & ! pointer array to backed up evolution of current microstructure
constitutive_previousDotState,& ! pointer array to previous evolution of current microstructure constitutive_previousDotState,& ! pointer array to previous evolution of current microstructure
constitutive_previousDotState2,&! pointer array to 2nd previous evolution of current microstructure constitutive_previousDotState2,&! pointer array to 2nd previous evolution of current microstructure
constitutive_dotState_backup, & ! pointer array to backed up evolution of current microstructure
constitutive_RK4dotState, & ! pointer array to evolution of microstructure defined by classical Runge-Kutta method
constitutive_relevantState ! relevant state values constitutive_relevantState ! relevant state values
type(p_vec), dimension(:,:,:,:), allocatable :: constitutive_RKCK45dotState ! pointer array to evolution of microstructure used by Cash-Karp Runge-Kutta method
integer(pInt), dimension(:,:,:), allocatable :: constitutive_sizeDotState, & ! size of dotState array integer(pInt), dimension(:,:,:), allocatable :: constitutive_sizeDotState, & ! size of dotState array
constitutive_sizeState, & ! size of state array per grain constitutive_sizeState, & ! size of state array per grain
constitutive_sizePostResults ! size of postResults array per grain constitutive_sizePostResults ! size of postResults array per grain
@ -48,6 +50,7 @@ subroutine constitutive_init()
!************************************** !**************************************
use prec, only: pReal,pInt use prec, only: pReal,pInt
use debug, only: debugger, selectiveDebugger, debug_e, debug_i, debug_g use debug, only: debugger, selectiveDebugger, debug_e, debug_i, debug_g
use numerics, only: integrator, integratorStiffness
use IO, only: IO_error, IO_open_file, IO_open_jobFile use IO, only: IO_error, IO_open_file, IO_open_jobFile
use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips use mesh, only: mesh_maxNips,mesh_NcpElems,mesh_element,FE_Nips
use material use material
@ -58,7 +61,7 @@ subroutine constitutive_init()
use constitutive_nonlocal use constitutive_nonlocal
integer(pInt), parameter :: fileunit = 200 integer(pInt), parameter :: fileunit = 200
integer(pInt) e,i,g,p,myInstance,myNgrains integer(pInt) e,i,g,p,s,myInstance,myNgrains
integer(pInt), dimension(:,:), pointer :: thisSize integer(pInt), dimension(:,:), pointer :: thisSize
character(len=64), dimension(:,:), pointer :: thisOutput character(len=64), dimension(:,:), pointer :: thisOutput
logical knownConstitution logical knownConstitution
@ -71,7 +74,6 @@ subroutine constitutive_init()
call constitutive_dislotwin_init(fileunit) call constitutive_dislotwin_init(fileunit)
call constitutive_nonlocal_init(fileunit) call constitutive_nonlocal_init(fileunit)
close(fileunit) close(fileunit)
! write description file for constitutive phase output ! write description file for constitutive phase output
@ -123,12 +125,18 @@ subroutine constitutive_init()
allocate(constitutive_state_backup(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(constitutive_state_backup(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(constitutive_dotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(constitutive_dotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(constitutive_dotState_backup(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(constitutive_dotState_backup(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(constitutive_previousDotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(constitutive_previousDotState2(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(constitutive_relevantState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(constitutive_relevantState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(constitutive_sizeDotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_sizeDotState = 0_pInt allocate(constitutive_sizeDotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_sizeDotState = 0_pInt
allocate(constitutive_sizeState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_sizeState = 0_pInt allocate(constitutive_sizeState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; constitutive_sizeState = 0_pInt
allocate(constitutive_sizePostResults(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)); constitutive_sizePostResults = 0_pInt allocate(constitutive_sizePostResults(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)); constitutive_sizePostResults = 0_pInt
if (integrator == 1 .or. integratorStiffness == 1) then
allocate(constitutive_previousDotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(constitutive_previousDotState2(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
endif
if (integrator == 3 .or. integratorStiffness == 3) &
allocate(constitutive_RK4dotState(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
if (integrator == 4 .or. integratorStiffness == 4) &
allocate(constitutive_RKCK45dotState(6,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
do e = 1,mesh_NcpElems ! loop over elements do e = 1,mesh_NcpElems ! loop over elements
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
@ -147,8 +155,17 @@ subroutine constitutive_init()
allocate(constitutive_relevantState(g,i,e)%p(constitutive_j2_sizeState(myInstance))) allocate(constitutive_relevantState(g,i,e)%p(constitutive_j2_sizeState(myInstance)))
allocate(constitutive_dotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) allocate(constitutive_dotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
if (integrator == 1 .or. integratorStiffness == 1) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) allocate(constitutive_previousDotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_j2_sizeDotState(myInstance))) allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
endif
if (integrator == 3 .or. integratorStiffness == 3) &
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
if (integrator == 4 .or. integratorStiffness == 4) then
do s = 1,6
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_j2_sizeDotState(myInstance)))
enddo
endif
constitutive_state0(g,i,e)%p = constitutive_j2_stateInit(myInstance) constitutive_state0(g,i,e)%p = constitutive_j2_stateInit(myInstance)
constitutive_relevantState(g,i,e)%p = constitutive_j2_relevantState(myInstance) constitutive_relevantState(g,i,e)%p = constitutive_j2_relevantState(myInstance)
constitutive_sizeState(g,i,e) = constitutive_j2_sizeState(myInstance) constitutive_sizeState(g,i,e) = constitutive_j2_sizeState(myInstance)
@ -164,8 +181,17 @@ subroutine constitutive_init()
allocate(constitutive_relevantState(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance))) allocate(constitutive_relevantState(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance)))
allocate(constitutive_dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) allocate(constitutive_dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
if (integrator == 1 .or. integratorStiffness == 1) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) allocate(constitutive_previousDotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance))) allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
endif
if (integrator == 3 .or. integratorStiffness == 3) &
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
if (integrator == 4 .or. integratorStiffness == 4) then
do s = 1,6
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_phenopowerlaw_sizeDotState(myInstance)))
enddo
endif
constitutive_state0(g,i,e)%p = constitutive_phenopowerlaw_stateInit(myInstance) constitutive_state0(g,i,e)%p = constitutive_phenopowerlaw_stateInit(myInstance)
constitutive_relevantState(g,i,e)%p = constitutive_phenopowerlaw_relevantState(myInstance) constitutive_relevantState(g,i,e)%p = constitutive_phenopowerlaw_relevantState(myInstance)
constitutive_sizeState(g,i,e) = constitutive_phenopowerlaw_sizeState(myInstance) constitutive_sizeState(g,i,e) = constitutive_phenopowerlaw_sizeState(myInstance)
@ -181,8 +207,17 @@ subroutine constitutive_init()
allocate(constitutive_relevantState(g,i,e)%p(constitutive_titanmod_sizeState(myInstance))) allocate(constitutive_relevantState(g,i,e)%p(constitutive_titanmod_sizeState(myInstance)))
allocate(constitutive_dotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) allocate(constitutive_dotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
if (integrator == 1 .or. integratorStiffness == 1) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) allocate(constitutive_previousDotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance))) allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
endif
if (integrator == 3 .or. integratorStiffness == 3) &
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
if (integrator == 4 .or. integratorStiffness == 4) then
do s = 1,6
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_titanmod_sizeDotState(myInstance)))
enddo
endif
constitutive_state0(g,i,e)%p = constitutive_titanmod_stateInit(myInstance) constitutive_state0(g,i,e)%p = constitutive_titanmod_stateInit(myInstance)
constitutive_relevantState(g,i,e)%p = constitutive_titanmod_relevantState(myInstance) constitutive_relevantState(g,i,e)%p = constitutive_titanmod_relevantState(myInstance)
constitutive_sizeState(g,i,e) = constitutive_titanmod_sizeState(myInstance) constitutive_sizeState(g,i,e) = constitutive_titanmod_sizeState(myInstance)
@ -198,8 +233,17 @@ subroutine constitutive_init()
allocate(constitutive_relevantState(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance))) allocate(constitutive_relevantState(g,i,e)%p(constitutive_dislotwin_sizeState(myInstance)))
allocate(constitutive_dotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) allocate(constitutive_dotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
if (integrator == 1 .or. integratorStiffness == 1) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) allocate(constitutive_previousDotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance))) allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
endif
if (integrator == 3 .or. integratorStiffness == 3) &
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
if (integrator == 4 .or. integratorStiffness == 4) then
do s = 1,6
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_dislotwin_sizeDotState(myInstance)))
enddo
endif
constitutive_state0(g,i,e)%p = constitutive_dislotwin_stateInit(myInstance) constitutive_state0(g,i,e)%p = constitutive_dislotwin_stateInit(myInstance)
constitutive_relevantState(g,i,e)%p = constitutive_dislotwin_relevantState(myInstance) constitutive_relevantState(g,i,e)%p = constitutive_dislotwin_relevantState(myInstance)
constitutive_sizeState(g,i,e) = constitutive_dislotwin_sizeState(myInstance) constitutive_sizeState(g,i,e) = constitutive_dislotwin_sizeState(myInstance)
@ -215,8 +259,17 @@ subroutine constitutive_init()
allocate(constitutive_relevantState(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance))) allocate(constitutive_relevantState(g,i,e)%p(constitutive_nonlocal_sizeState(myInstance)))
allocate(constitutive_dotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) allocate(constitutive_dotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) allocate(constitutive_dotState_backup(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
if (integrator == 1 .or. integratorStiffness == 1) then
allocate(constitutive_previousDotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) allocate(constitutive_previousDotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance))) allocate(constitutive_previousDotState2(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
endif
if (integrator == 3 .or. integratorStiffness == 3) &
allocate(constitutive_RK4dotState(g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
if (integrator == 4 .or. integratorStiffness == 4) then
do s = 1,6
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_nonlocal_sizeDotState(myInstance)))
enddo
endif
constitutive_state0(g,i,e)%p = constitutive_nonlocal_stateInit(myInstance) constitutive_state0(g,i,e)%p = constitutive_nonlocal_stateInit(myInstance)
constitutive_relevantState(g,i,e)%p = constitutive_nonlocal_relevantState(myInstance) constitutive_relevantState(g,i,e)%p = constitutive_nonlocal_relevantState(myInstance)
constitutive_sizeState(g,i,e) = constitutive_nonlocal_sizeState(myInstance) constitutive_sizeState(g,i,e) = constitutive_nonlocal_sizeState(myInstance)
@ -247,7 +300,6 @@ subroutine constitutive_init()
write(6,'(a32,x,7(i5,x))') 'constitutive_state: ', shape(constitutive_state) write(6,'(a32,x,7(i5,x))') 'constitutive_state: ', shape(constitutive_state)
write(6,'(a32,x,7(i5,x))') 'constitutive_relevantState: ', shape(constitutive_relevantState) write(6,'(a32,x,7(i5,x))') 'constitutive_relevantState: ', shape(constitutive_relevantState)
write(6,'(a32,x,7(i5,x))') 'constitutive_dotState: ', shape(constitutive_dotState) write(6,'(a32,x,7(i5,x))') 'constitutive_dotState: ', shape(constitutive_dotState)
write(6,'(a32,x,7(i5,x))') 'constitutive_previousDotState:', shape(constitutive_previousDotState)
write(6,'(a32,x,7(i5,x))') 'constitutive_sizeState: ', shape(constitutive_sizeState) write(6,'(a32,x,7(i5,x))') 'constitutive_sizeState: ', shape(constitutive_sizeState)
write(6,'(a32,x,7(i5,x))') 'constitutive_sizeDotState: ', shape(constitutive_sizeDotState) write(6,'(a32,x,7(i5,x))') 'constitutive_sizeDotState: ', shape(constitutive_sizeDotState)
write(6,'(a32,x,7(i5,x))') 'constitutive_sizePostResults: ', shape(constitutive_sizePostResults) write(6,'(a32,x,7(i5,x))') 'constitutive_sizePostResults: ', shape(constitutive_sizePostResults)
@ -448,7 +500,8 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, Temperature, ip
call constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el) call constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_nonlocal_label) case (constitutive_nonlocal_label)
call constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, Temperature, constitutive_state, ipc, ip, el) call constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, Temperature, constitutive_state, &
constitutive_relevantState, ipc, ip, el)
end select end select
@ -478,11 +531,16 @@ use mesh, only: mesh_NcpElems, &
use material, only: phase_constitution, & use material, only: phase_constitution, &
material_phase, & material_phase, &
homogenization_maxNgrains homogenization_maxNgrains
use constitutive_j2 use constitutive_j2, only: constitutive_j2_dotState, &
use constitutive_phenopowerlaw constitutive_j2_label
use constitutive_titanmod use constitutive_phenopowerlaw, only: constitutive_phenopowerlaw_dotState, &
use constitutive_dislotwin constitutive_phenopowerlaw_label
use constitutive_nonlocal use constitutive_titanmod, only: constitutive_titanmod_dotState, &
constitutive_titanmod_label
use constitutive_dislotwin, only: constitutive_dislotwin_dotState, &
constitutive_dislotwin_label
use constitutive_nonlocal, only: constitutive_nonlocal_dotState, &
constitutive_nonlocal_label
implicit none implicit none
@ -549,12 +607,19 @@ function constitutive_dotTemperature(Tstar_v,Temperature,ipc,ip,el)
!* - constitutive_dotTemperature : evolution of temperature * !* - constitutive_dotTemperature : evolution of temperature *
!********************************************************************* !*********************************************************************
use prec, only: pReal,pInt use prec, only: pReal,pInt
use debug, only: debug_cumDotTemperatureCalls, &
debug_cumDotTemperatureTicks
use material, only: phase_constitution,material_phase use material, only: phase_constitution,material_phase
use constitutive_j2 use constitutive_j2, only: constitutive_j2_dotTemperature, &
use constitutive_phenopowerlaw constitutive_j2_label
use constitutive_titanmod use constitutive_phenopowerlaw, only: constitutive_phenopowerlaw_dotTemperature, &
use constitutive_dislotwin constitutive_phenopowerlaw_label
use constitutive_nonlocal use constitutive_titanmod, only: constitutive_titanmod_dotTemperature, &
constitutive_titanmod_label
use constitutive_dislotwin, only: constitutive_dislotwin_dotTemperature, &
constitutive_dislotwin_label
use constitutive_nonlocal, only: constitutive_nonlocal_dotTemperature, &
constitutive_nonlocal_label
implicit none implicit none
!* Definition of variables !* Definition of variables
@ -562,6 +627,12 @@ function constitutive_dotTemperature(Tstar_v,Temperature,ipc,ip,el)
real(pReal) Temperature real(pReal) Temperature
real(pReal) constitutive_dotTemperature real(pReal) constitutive_dotTemperature
real(pReal), dimension(6) :: Tstar_v real(pReal), dimension(6) :: Tstar_v
integer(pLongInt) tick, &
tock, &
tickrate, &
maxticks
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
select case (phase_constitution(material_phase(ipc,ip,el))) select case (phase_constitution(material_phase(ipc,ip,el)))
@ -581,6 +652,12 @@ function constitutive_dotTemperature(Tstar_v,Temperature,ipc,ip,el)
constitutive_dotTemperature = constitutive_nonlocal_dotTemperature(Tstar_v,Temperature,constitutive_state,ipc,ip,el) constitutive_dotTemperature = constitutive_nonlocal_dotTemperature(Tstar_v,Temperature,constitutive_state,ipc,ip,el)
end select end select
call system_clock(count=tock,count_rate=tickrate,count_max=maxticks)
debug_cumDotTemperatureCalls = debug_cumDotTemperatureCalls + 1_pInt
debug_cumDotTemperatureTicks = debug_cumDotTemperatureTicks + tock-tick
if (tock < tick) debug_cumDotTemperatureTicks = debug_cumDotTemperatureTicks + maxticks
return return
endfunction endfunction

78
code/constitutive_nonlocal.f90
View File

@ -257,7 +257,6 @@ constitutive_nonlocal_interactionSlipSlip = 0.0_pReal
constitutive_nonlocal_dLowerEdgePerSlipFamily = 0.0_pReal constitutive_nonlocal_dLowerEdgePerSlipFamily = 0.0_pReal
constitutive_nonlocal_dLowerScrewPerSlipFamily = 0.0_pReal constitutive_nonlocal_dLowerScrewPerSlipFamily = 0.0_pReal
!*** readout data from material.config file !*** readout data from material.config file
rewind(file) rewind(file)
@ -973,7 +972,7 @@ endsubroutine
!********************************************************************* !*********************************************************************
!* calculates kinetics * !* calculates kinetics *
!********************************************************************* !*********************************************************************
subroutine constitutive_nonlocal_kinetics(Tstar_v, Temperature, state, g, ip, el) subroutine constitutive_nonlocal_kinetics(Tstar_v, Temperature, state, g, ip, el, dv_dtau)
use prec, only: pReal, & use prec, only: pReal, &
pInt, & pInt, &
@ -1003,6 +1002,8 @@ type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), in
real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola-Kirchhoff stress in Mandel notation real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola-Kirchhoff stress in Mandel notation
!*** output variables !*** output variables
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))), &
intent(out), optional :: dv_dtau ! velocity derivative with respect to resolved shear stress
!*** local variables !*** local variables
integer(pInt) myInstance, & ! current instance of this constitution integer(pInt) myInstance, & ! current instance of this constitution
@ -1013,31 +1014,42 @@ integer(pInt) myInstance, & ! curren
real(pReal), dimension(6) :: Tdislocation_v ! dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress real(pReal), dimension(6) :: Tdislocation_v ! dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: & real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: &
tauThreshold, & ! threshold shear stress tauThreshold, & ! threshold shear stress
tau ! resolved shear stress tau, & ! resolved shear stress
rhoForest ! forest dislocation density
real(pReal) boltzmannProbability
myInstance = phase_constitutionInstance(material_phase(g,ip,el)) myInstance = phase_constitutionInstance(material_phase(g,ip,el))
myStructure = constitutive_nonlocal_structure(myInstance) myStructure = constitutive_nonlocal_structure(myInstance)
ns = constitutive_nonlocal_totalNslip(myInstance) ns = constitutive_nonlocal_totalNslip(myInstance)
rhoForest = state(g,ip,el)%p(10*ns+1:11*ns)
tauThreshold = state(g,ip,el)%p(11*ns+1:12*ns) tauThreshold = state(g,ip,el)%p(11*ns+1:12*ns)
Tdislocation_v = state(g,ip,el)%p(12*ns+1:12*ns+6) Tdislocation_v = state(g,ip,el)%p(12*ns+1:12*ns+6)
tau = 0.0_pReal tau = 0.0_pReal
constitutive_nonlocal_v(:,:,g,ip,el) = 0.0_pReal constitutive_nonlocal_v(:,:,g,ip,el) = 0.0_pReal
if (present(dv_dtau)) dv_dtau = 0.0_pReal
if ( Temperature > 0.0_pReal ) then
do s = 1,ns do s = 1,ns
if ((tauThreshold(s) > 0.0_pReal) .and. (Temperature > 0.0_pReal)) then
tau(s) = math_mul6x6( Tstar_v + Tdislocation_v, & tau(s) = math_mul6x6( Tstar_v + Tdislocation_v, &
lattice_Sslip_v(:,constitutive_nonlocal_slipSystemLattice(s,myInstance),myStructure) ) lattice_Sslip_v(:,constitutive_nonlocal_slipSystemLattice(s,myInstance),myStructure) )
if ( abs(tau(s)) > 0.0_pReal ) then
boltzmannProbability = dexp( -constitutive_nonlocal_Q0(myInstance) * dsqrt(rhoForest(s)) / ( abs(tau(s)) * kB * Temperature) )
constitutive_nonlocal_v(s,:,g,ip,el) = constitutive_nonlocal_v0PerSlipSystem(s,myInstance) & constitutive_nonlocal_v(s,:,g,ip,el) = constitutive_nonlocal_v0PerSlipSystem(s,myInstance) &
* exp( - constitutive_nonlocal_Q0(myInstance) / ( kB * Temperature) * (1.0_pReal - (abs(tau(s))/tauThreshold(s)) ) ) & / ( constitutive_nonlocal_burgersPerSlipSystem(s,myInstance) * dsqrt(rhoForest(s)) ) &
* sign(1.0_pReal,tau(s)) * boltzmannProbability * sign(1.0_pReal,tau(s))
if (present(dv_dtau)) &
dv_dtau(s) = constitutive_nonlocal_Q0(myInstance) * constitutive_nonlocal_v0PerSlipSystem(s,myInstance) &
/ ( constitutive_nonlocal_burgersPerSlipSystem(s,myInstance) * tau(s)**2.0_pReal * kB * Temperature ) &
* boltzmannProbability
endif endif
enddo enddo
endif
if (verboseDebugger .and. selectiveDebugger) then if (verboseDebugger .and. selectiveDebugger) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
@ -1058,7 +1070,7 @@ endsubroutine
!********************************************************************* !*********************************************************************
!* calculates plastic velocity gradient and its tangent * !* calculates plastic velocity gradient and its tangent *
!********************************************************************* !*********************************************************************
subroutine constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar99, Tstar_v, Temperature, state, g, ip, el) subroutine constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar99, Tstar_v, Temperature, state, relevantState, g, ip, el)
use prec, only: pReal, & use prec, only: pReal, &
pInt, & pInt, &
@ -1085,7 +1097,8 @@ integer(pInt), intent(in) :: g, & ! curren
el ! current element number el ! current element number
real(pReal), intent(in) :: Temperature ! temperature real(pReal), intent(in) :: Temperature ! temperature
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 ! microstructural state state, & ! microstructural state
relevantState ! relevant microstructural state
real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola-Kirchhoff stress in Mandel notation real(pReal), dimension(6), intent(in) :: Tstar_v ! 2nd Piola-Kirchhoff stress in Mandel notation
!*** output variables !*** output variables
@ -1111,7 +1124,9 @@ real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstan
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: & real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: &
tauThreshold, & ! threshold shear stress tauThreshold, & ! threshold shear stress
gdotTotal, & ! shear rate gdotTotal, & ! shear rate
dgdotTotal_dtau ! derivative of the shear rate with respect to the shear stress dv_dtau, & ! velocity derivative with respect to the shear stress
dgdotTotal_dtau, & ! derivative of the shear rate with respect to the shear stress
rhoForest ! forest dislocation density
!*** initialize local variables !*** initialize local variables
@ -1131,17 +1146,18 @@ forall (t = 1:8) &
forall (s = 1:ns, t = 5:8, rhoSgl(s,t) * constitutive_nonlocal_v(s,t-4,g,ip,el) < 0.0_pReal) & ! contribution of used rho for changing sign of v forall (s = 1:ns, t = 5:8, rhoSgl(s,t) * constitutive_nonlocal_v(s,t-4,g,ip,el) < 0.0_pReal) & ! contribution of used rho for changing sign of v
rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t)) rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t))
rhoForest = state(g,ip,el)%p(10*ns+1:11*ns)
tauThreshold = state(g,ip,el)%p(11*ns+1:12*ns) tauThreshold = state(g,ip,el)%p(11*ns+1:12*ns)
call constitutive_nonlocal_kinetics(Tstar_v, Temperature, state, g, ip, el) ! update dislocation velocity call constitutive_nonlocal_kinetics(Tstar_v, Temperature, state, g, ip, el, dv_dtau) ! update dislocation velocity
!*** Calculation of gdot and its tangent !*** Calculation of gdot and its tangent
forall (t = 1:4 ) & forall (s = 1:ns, t = 1:4, rhoSgl(s,t) > relevantState(g,ip,el)%p((t-1)*ns+s)) & ! no shear rate contribution for densities below relevant state
gdot(:,t) = rhoSgl(:,t) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) * constitutive_nonlocal_v(:,t,g,ip,el) gdot(s,t) = rhoSgl(s,t) * constitutive_nonlocal_burgersPerSlipSystem(s,myInstance) * constitutive_nonlocal_v(s,t,g,ip,el)
gdotTotal = sum(gdot,2) gdotTotal = sum(gdot,2)
dgdotTotal_dtau = abs(gdotTotal) * constitutive_nonlocal_Q0(myInstance) / ( kB * Temperature * tauThreshold ) dgdotTotal_dtau = sum(rhoSgl,2) * constitutive_nonlocal_burgersPerSlipSystem(:,myInstance) * dv_dtau
!*** Calculation of Lp and its tangent !*** Calculation of Lp and its tangent
@ -1217,6 +1233,7 @@ use lattice, only: lattice_Sslip, &
lattice_st, & lattice_st, &
lattice_maxNslipFamily, & lattice_maxNslipFamily, &
lattice_NslipSystem lattice_NslipSystem
use FEsolving, only:theInc
implicit none implicit none
!*** input variables !*** input variables
@ -1352,16 +1369,16 @@ if (timestep <= 0.0_pReal) then
return return
endif endif
if (any(constitutive_nonlocal_v(:,:,g,ip,el)*timestep > mesh_ipVolume(ip,el)**(1.0_pReal/3.0_pReal))) then ! if timestep is too large,... !if (any(constitutive_nonlocal_v(:,:,g,ip,el)*timestep > mesh_ipVolume(ip,el)**(1.0_pReal/3.0_pReal))) then ! if timestep is too large,...
dotState(1,ip,el)%p(1:10*ns) = NaN ! ...assign NaN and enforce a cutback ! dotState(1,ip,el)%p(1:10*ns) = NaN ! ...assign NaN and enforce a cutback
if (verboseDebugger) then ! if (verboseDebugger) then
!$OMP CRITICAL (write2out) ! !$OMP CRITICAL (write2out)
write(6,*) 'exceeded maximum allowed dislocation velocity at ',g,ip,el ! write(6,*) 'exceeded maximum allowed dislocation velocity at ',g,ip,el
write(6,*) ! write(6,*)
!$OMPEND CRITICAL (write2out) ! !$OMPEND CRITICAL (write2out)
endif ! endif
return ! return
endif !endif
!**************************************************************************** !****************************************************************************
@ -1604,19 +1621,22 @@ forall (t = 1:10) &
+ rhoDotRemobilization(:,t) & + rhoDotRemobilization(:,t) &
+ rhoDotMultiplication(:,t) & + rhoDotMultiplication(:,t) &
+ rhoDotThermalAnnihilation(:,t) + rhoDotThermalAnnihilation(:,t)
! + rhoDotDipole2SingleStressChange(:,t) & ! + rhoDotDipole2SingleStressChange(:,t)
! + rhoDotSingle2DipoleStressChange(:,t) ! + rhoDotSingle2DipoleStressChange(:,t)
dotState(g,ip,el)%p(1:10*ns) = dotState(g,ip,el)%p(1:10*ns) + reshape(rhoDot,(/10*ns/)) dotState(g,ip,el)%p(1:10*ns) = reshape(rhoDot,(/10*ns/))
do i = 1,4*ns do i = 1,10*ns
if (i > 4*ns .and. i <= 8*ns) & ! skip immobile densities
continue
if (previousState(g,ip,el)%p(i) + dotState(g,ip,el)%p(i)*timestep < 0.0_pReal) then ! if single mobile densities become negative... if (previousState(g,ip,el)%p(i) + dotState(g,ip,el)%p(i)*timestep < 0.0_pReal) then ! if single mobile densities become negative...
if (previousState(g,ip,el)%p(i) < relevantState(g,ip,el)%p(i)) then ! ... and density is already below relevance... if (previousState(g,ip,el)%p(i) < relevantState(g,ip,el)%p(i)) then ! ... and density is already below relevance...
dotState(g,ip,el)%p(i) = 0.0_pReal ! ... set dotState to zero dotState(g,ip,el)%p(i) = - previousState(g,ip,el)%p(i) / timestep ! ... set dotState to zero
else ! ... otherwise... else ! ... otherwise...
if (verboseDebugger) then if (verboseDebugger) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) 'negative dislocation density at ',g,ip,el write(6,*) 'negative dislocation density: enforced cutback at ',g,ip,el,i
write(6,'(a12,x,e15.8)') 'dotState was',dotState(g,ip,el)%p(i)
write(6,*) write(6,*)
!$OMPEND CRITICAL (write2out) !$OMPEND CRITICAL (write2out)
endif endif
@ -1671,7 +1691,7 @@ disorientationAxisAngle = math_QuaternionToAxisAngle(disorientation)
disorientationAxis = disorientationAxisAngle(1:3) disorientationAxis = disorientationAxisAngle(1:3)
disorientationAngle = disorientationAxisAngle(4) disorientationAngle = disorientationAxisAngle(4)
if (disorientationAngle < 3.0_pReal) then if (disorientationAngle < 5.0_pReal) then
constitutive_nonlocal_transmissivity = 1.0_pReal constitutive_nonlocal_transmissivity = 1.0_pReal
else else
constitutive_nonlocal_transmissivity = 0.5_pReal constitutive_nonlocal_transmissivity = 0.5_pReal

1935
code/crystallite.f90
View File

File diff suppressed because it is too large Load Diff

8
code/debug.config
View File

@ -3,7 +3,7 @@
debug 1 # >0 true to switch on debugging debug 1 # >0 true to switch on debugging
verbose 1 # >0 true to switch on verbose output verbose 1 # >0 true to switch on verbose output
selective 1 # >0 true to switch on e,i,g seelctive debugging selective 1 # >0 true to switch on e,i,g selective debugging
element 123 # selected element for debugging (synonymous: "el", "e") element 1 # selected element for debugging (synonymous: "el", "e")
ip 3 # selected integration point for debugging (synonymous: "integrationpoint", "i") ip 1 # selected integration point for debugging (synonymous: "integrationpoint", "i")
grain 5 # selected grain at ip for debugging (synonymous: "gr", "g") grain 1 # selected grain at ip for debugging (synonymous: "gr", "g")

31
code/debug.f90
View File

@ -9,8 +9,7 @@
integer(pInt), dimension(:,:), allocatable :: debug_StressLoopDistribution integer(pInt), dimension(:,:), allocatable :: debug_StressLoopDistribution
integer(pInt), dimension(:,:), allocatable :: debug_LeapfrogBreakDistribution integer(pInt), dimension(:,:), allocatable :: debug_LeapfrogBreakDistribution
integer(pInt), dimension(:), allocatable :: debug_CrystalliteStateLoopDistribution integer(pInt), dimension(:,:), allocatable :: debug_StateLoopDistribution
integer(pInt), dimension(:), allocatable :: debug_StiffnessStateLoopDistribution
integer(pInt), dimension(:), allocatable :: debug_CrystalliteLoopDistribution integer(pInt), dimension(:), allocatable :: debug_CrystalliteLoopDistribution
integer(pInt), dimension(:), allocatable :: debug_MaterialpointStateLoopDistribution integer(pInt), dimension(:), allocatable :: debug_MaterialpointStateLoopDistribution
integer(pInt), dimension(:), allocatable :: debug_MaterialpointLoopDistribution integer(pInt), dimension(:), allocatable :: debug_MaterialpointLoopDistribution
@ -70,8 +69,7 @@ subroutine debug_init()
allocate(debug_StressLoopDistribution(nStress,2)) ; debug_StressLoopDistribution = 0_pInt allocate(debug_StressLoopDistribution(nStress,2)) ; debug_StressLoopDistribution = 0_pInt
allocate(debug_LeapfrogBreakDistribution(nStress,2)) ; debug_LeapfrogBreakDistribution = 0_pInt allocate(debug_LeapfrogBreakDistribution(nStress,2)) ; debug_LeapfrogBreakDistribution = 0_pInt
allocate(debug_CrystalliteStateLoopDistribution(nState)) ; debug_CrystalliteStateLoopDistribution = 0_pInt allocate(debug_StateLoopDistribution(nState,2)) ; debug_StateLoopDistribution = 0_pInt
allocate(debug_StiffnessStateLoopDistribution(nState)) ; debug_StiffnessStateLoopDistribution = 0_pInt
allocate(debug_CrystalliteLoopDistribution(nCryst+1)) ; debug_CrystalliteLoopDistribution = 0_pInt allocate(debug_CrystalliteLoopDistribution(nCryst+1)) ; debug_CrystalliteLoopDistribution = 0_pInt
allocate(debug_MaterialpointStateLoopDistribution(nMPstate)) ; debug_MaterialpointStateLoopDistribution = 0_pInt allocate(debug_MaterialpointStateLoopDistribution(nMPstate)) ; debug_MaterialpointStateLoopDistribution = 0_pInt
allocate(debug_MaterialpointLoopDistribution(nHomog+1)) ; debug_MaterialpointLoopDistribution = 0_pInt allocate(debug_MaterialpointLoopDistribution(nHomog+1)) ; debug_MaterialpointLoopDistribution = 0_pInt
@ -141,8 +139,7 @@ subroutine debug_reset()
debug_StressLoopDistribution = 0_pInt ! initialize debugging data debug_StressLoopDistribution = 0_pInt ! initialize debugging data
debug_LeapfrogBreakDistribution = 0_pInt debug_LeapfrogBreakDistribution = 0_pInt
debug_CrystalliteStateLoopDistribution = 0_pInt debug_StateLoopDistribution = 0_pInt
debug_StiffnessStateLoopDistribution = 0_pInt
debug_CrystalliteLoopDistribution = 0_pInt debug_CrystalliteLoopDistribution = 0_pInt
debug_MaterialpointStateLoopDistribution = 0_pInt debug_MaterialpointStateLoopDistribution = 0_pInt
debug_MaterialpointLoopDistribution = 0_pInt debug_MaterialpointLoopDistribution = 0_pInt
@ -216,12 +213,14 @@ endsubroutine
write(6,*) write(6,*)
write(6,*) 'distribution_CrystalliteStateLoop :' write(6,*) 'distribution_CrystalliteStateLoop :'
do i=1,nState do i=1,nState
if (debug_CrystalliteStateLoopDistribution(i) /= 0) then if (any(debug_StateLoopDistribution(i,:) /= 0)) then
integral = integral + i*debug_CrystalliteStateLoopDistribution(i) integral = integral + i*debug_StateLoopDistribution(i,1) + i*debug_StateLoopDistribution(i,2)
write(6,'(i25,x,i10)') i,debug_CrystalliteStateLoopDistribution(i) write(6,'(i25,x,i10,12x,i10)') i,debug_StateLoopDistribution(i,1),debug_StateLoopDistribution(i,2)
endif endif
enddo enddo
write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_CrystalliteStateLoopDistribution) write(6,'(a15,i10,x,i10,12x,i10)') ' total',integral,&
sum(debug_StateLoopDistribution(:,1)), &
sum(debug_StateLoopDistribution(:,2))
integral = 0_pInt integral = 0_pInt
write(6,*) write(6,*)
@ -238,18 +237,6 @@ endsubroutine
enddo enddo
write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_CrystalliteLoopDistribution) write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_CrystalliteLoopDistribution)
integral = 0_pInt
write(6,*)
write(6,*) 'distribution_StiffnessStateLoop :'
do i=1,nState
if (debug_StiffnessStateLoopDistribution(i) /= 0) then
integral = integral + i*debug_StiffnessStateLoopDistribution(i)
write(6,'(i25,x,i10)') i,debug_StiffnessStateLoopDistribution(i)
endif
enddo
write(6,'(a15,i10,x,i10)') ' total',integral,sum(debug_StiffnessStateLoopDistribution)
!* Material point loop counter <<<updated 31.07.2009>>>
integral = 0_pInt integral = 0_pInt
write(6,*) write(6,*)
write(6,*) write(6,*)

7
code/homogenization.f90
View File

@ -229,6 +229,7 @@ subroutine materialpoint_stressAndItsTangent(&
stepIncreaseHomog, & stepIncreaseHomog, &
nHomog, & nHomog, &
nMPstate nMPstate
use math, only: math_det3x3
use FEsolving, only: FEsolving_execElem, & use FEsolving, only: FEsolving_execElem, &
FEsolving_execIP, & FEsolving_execIP, &
terminallyIll terminallyIll
@ -243,6 +244,8 @@ subroutine materialpoint_stressAndItsTangent(&
crystallite_Fp, & crystallite_Fp, &
crystallite_Lp0, & crystallite_Lp0, &
crystallite_Lp, & crystallite_Lp, &
crystallite_dPdF, &
crystallite_dPdF0, &
crystallite_Tstar0_v, & crystallite_Tstar0_v, &
crystallite_Tstar_v, & crystallite_Tstar_v, &
crystallite_partionedTemperature0, & crystallite_partionedTemperature0, &
@ -250,6 +253,7 @@ subroutine materialpoint_stressAndItsTangent(&
crystallite_partionedF, & crystallite_partionedF, &
crystallite_partionedFp0, & crystallite_partionedFp0, &
crystallite_partionedLp0, & crystallite_partionedLp0, &
crystallite_partioneddPdF0, &
crystallite_partionedTstar0_v, & crystallite_partionedTstar0_v, &
crystallite_dt, & crystallite_dt, &
crystallite_requested, & crystallite_requested, &
@ -296,6 +300,7 @@ subroutine materialpoint_stressAndItsTangent(&
crystallite_partionedTemperature0(1:myNgrains,i,e) = materialpoint_Temperature(i,e) ! ...temperatures crystallite_partionedTemperature0(1:myNgrains,i,e) = materialpoint_Temperature(i,e) ! ...temperatures
crystallite_partionedFp0(:,:,1:myNgrains,i,e) = crystallite_Fp0(:,:,1:myNgrains,i,e) ! ...plastic def grads crystallite_partionedFp0(:,:,1:myNgrains,i,e) = crystallite_Fp0(:,:,1:myNgrains,i,e) ! ...plastic def grads
crystallite_partionedLp0(:,:,1:myNgrains,i,e) = crystallite_Lp0(:,:,1:myNgrains,i,e) ! ...plastic velocity grads crystallite_partionedLp0(:,:,1:myNgrains,i,e) = crystallite_Lp0(:,:,1:myNgrains,i,e) ! ...plastic velocity grads
crystallite_partioneddPdF0(:,:,:,:,1:myNgrains,i,e) = crystallite_dPdF0(:,:,:,:,1:myNgrains,i,e) ! ...stiffness
crystallite_partionedF0(:,:,1:myNgrains,i,e) = crystallite_F0(:,:,1:myNgrains,i,e) ! ...def grads crystallite_partionedF0(:,:,1:myNgrains,i,e) = crystallite_F0(:,:,1:myNgrains,i,e) ! ...def grads
crystallite_partionedTstar0_v(:,1:myNgrains,i,e)= crystallite_Tstar0_v(:,1:myNgrains,i,e) ! ...2nd PK stress crystallite_partionedTstar0_v(:,1:myNgrains,i,e)= crystallite_Tstar0_v(:,1:myNgrains,i,e) ! ...2nd PK stress
@ -348,6 +353,7 @@ subroutine materialpoint_stressAndItsTangent(&
crystallite_partionedF0(:,:,1:myNgrains,i,e) = crystallite_partionedF(:,:,1:myNgrains,i,e) ! ...def grads crystallite_partionedF0(:,:,1:myNgrains,i,e) = crystallite_partionedF(:,:,1:myNgrains,i,e) ! ...def grads
crystallite_partionedFp0(:,:,1:myNgrains,i,e) = crystallite_Fp(:,:,1:myNgrains,i,e) ! ...plastic def grads crystallite_partionedFp0(:,:,1:myNgrains,i,e) = crystallite_Fp(:,:,1:myNgrains,i,e) ! ...plastic def grads
crystallite_partionedLp0(:,:,1:myNgrains,i,e) = crystallite_Lp(:,:,1:myNgrains,i,e) ! ...plastic velocity grads crystallite_partionedLp0(:,:,1:myNgrains,i,e) = crystallite_Lp(:,:,1:myNgrains,i,e) ! ...plastic velocity grads
crystallite_partioneddPdF0(:,:,:,:,1:myNgrains,i,e) = crystallite_dPdF(:,:,:,:,1:myNgrains,i,e)! ...stiffness
crystallite_partionedTstar0_v(:,1:myNgrains,i,e) = crystallite_Tstar_v(:,1:myNgrains,i,e) ! ...2nd PK stress crystallite_partionedTstar0_v(:,1:myNgrains,i,e) = crystallite_Tstar_v(:,1:myNgrains,i,e) ! ...2nd PK stress
forall (g = 1:myNgrains) constitutive_partionedState0(g,i,e)%p = constitutive_state(g,i,e)%p ! ...microstructures forall (g = 1:myNgrains) constitutive_partionedState0(g,i,e)%p = constitutive_state(g,i,e)%p ! ...microstructures
if (homogenization_sizeState(i,e) > 0_pInt) & if (homogenization_sizeState(i,e) > 0_pInt) &
@ -382,6 +388,7 @@ subroutine materialpoint_stressAndItsTangent(&
! ...initial def grad unchanged ! ...initial def grad unchanged
crystallite_Fp(:,:,1:myNgrains,i,e) = crystallite_partionedFp0(:,:,1:myNgrains,i,e) ! ...plastic def grads crystallite_Fp(:,:,1:myNgrains,i,e) = crystallite_partionedFp0(:,:,1:myNgrains,i,e) ! ...plastic def grads
crystallite_Lp(:,:,1:myNgrains,i,e) = crystallite_partionedLp0(:,:,1:myNgrains,i,e) ! ...plastic velocity grads crystallite_Lp(:,:,1:myNgrains,i,e) = crystallite_partionedLp0(:,:,1:myNgrains,i,e) ! ...plastic velocity grads
crystallite_dPdF(:,:,:,:,1:myNgrains,i,e) = crystallite_partioneddPdF0(:,:,:,:,1:myNgrains,i,e) ! ...stiffness
crystallite_Tstar_v(:,1:myNgrains,i,e) = crystallite_partionedTstar0_v(:,1:myNgrains,i,e) ! ...2nd PK stress crystallite_Tstar_v(:,1:myNgrains,i,e) = crystallite_partionedTstar0_v(:,1:myNgrains,i,e) ! ...2nd PK stress
forall (g = 1:myNgrains) constitutive_state(g,i,e)%p = constitutive_partionedState0(g,i,e)%p ! ...microstructures forall (g = 1:myNgrains) constitutive_state(g,i,e)%p = constitutive_partionedState0(g,i,e)%p ! ...microstructures
if (homogenization_sizeState(i,e) > 0_pInt) & if (homogenization_sizeState(i,e) > 0_pInt) &

79
code/lattice.f90
View File

@ -26,6 +26,9 @@ integer(pInt), parameter :: lattice_maxNslip = 48 ! max # of sli
integer(pInt), parameter :: lattice_maxNtwin = 24 ! max # of twin systems over lattice structures integer(pInt), parameter :: lattice_maxNtwin = 24 ! max # of twin systems over lattice structures
integer(pInt), parameter :: lattice_maxNinteraction = 20 ! max # of interaction types (in hardening matrix part) integer(pInt), parameter :: lattice_maxNinteraction = 20 ! max # of interaction types (in hardening matrix part)
integer(pInt), parameter, dimension(3) :: lattice_symmetryTypes =(/1, 1, 2/) ! maps crystal structures to symmetry tpyes
integer(pInt), pointer, dimension(:,:) :: interactionSlipSlip, & integer(pInt), pointer, dimension(:,:) :: interactionSlipSlip, &
interactionSlipTwin, & interactionSlipTwin, &
interactionTwinSlip, & interactionTwinSlip, &
@ -429,37 +432,37 @@ integer(pInt), allocatable, dimension(:,:,:) :: lattice_interactionSlipSlip, &
integer(pInt), parameter :: lattice_hex_Ntwin = 24 ! sum(lattice_hex_NtwinSystem) integer(pInt), parameter :: lattice_hex_Ntwin = 24 ! sum(lattice_hex_NtwinSystem)
integer(pInt) :: lattice_hex_Nstructure = 0_pInt integer(pInt) :: lattice_hex_Nstructure = 0_pInt
!* sorted by A. Alankar & P. Eisenlohr !* sorted by YJ.Ro and Philip
real(pReal), dimension(4+4,lattice_hex_Nslip), parameter :: lattice_hex_systemSlip = & real(pReal), dimension(4+4,lattice_hex_Nslip), parameter :: lattice_hex_systemSlip = &
reshape((/& reshape((/&
! Basal systems <1120>{0001} (independent of c/a-ratio, Bravais notation (4 coordinate base)) ! Basal systems <1120>{0001} (independent of c/a-ratio, Bravais notation (4 coordinate base))
2, -1, -1, 0, 0, 0, 0, 1, & !A1 2, -1, -1, 0, 0, 0, 0, 1, &
-1, 2, -1, 0, 0, 0, 0, 1, & !A2 1, 1, -2, 0, 0, 0, 0, 1, &
-1, -1, 2, 0, 0, 0, 0, 1, & !A3 -1, 2, -1, 0, 0, 0, 0, 1, &
! 1st type prismatic systems <1120>{1010} (independent of c/a-ratio) ! 1st type prismatic systems <1120>{1010} (independent of c/a-ratio)
2, -1, -1, 0, 0, 1, -1, 0, & !B1 2, -1, -1, 0, 0, -1, 1, 0, &
-1, 2, -1, 0, -1, 0, 1, 0, & !C2 1, 1, 2, 0, 1, -1, 0, 0, &
-1, -1, 2, 0, 1, -1, 0, 0, & !D3 -1, 2, -1, 0, 1, 0, -1, 0, &
! 1st type 1st order pyramidal systems <1120>{1011} -- plane normals depend on the c/a-ratio ! 1st type 1st order pyramidal systems <1120>{1011}
2, -1, -1, 0, 0, 1, -1, 1, & !E1 2, -1, -1, 0, 0, -1, 1, 1, &
1, 1, -2, 0, -1, 1, 0, 1, & !F(-3) 1, 1, -2, 0, 1, -1, 0, 1, &
-1, 2, -1, 0, -1, 0, 1, 1, & !G2 -1, 2, -1, 0, 1, 0, -1, 1, &
-2, 1, 1, 0, 0, -1, 1, 1, & !H(-1) -2, 1, 1, 0, 0, 1, -1, 1, &
-1, -1, 2, 0, 1, -1, 0, 1, & !I3 -1, -1, 2, 0, -1, 1, 0, 1, &
1, -2, 1, 0, 1, 0, -1, 1, & !J(-2) 1, -2, 1, 0, -1, 0, 1, 1, &
! pyramidal system: c+a slip <2113>{1011} -- plane normals depend on the c/a-ratio ! pyramidal system: c+a slip <2113>{1011} -- plane normals depend on the c/a-ratio
-1, 2, -1, 3, 0, 1, -1, 1, & !E4 2, -1, -1, -3, 1, -1, 0, 1, &
1, 1, -2, 3, 0, 1, -1, 1, & !E5 1, 1, -2, -3, 1, 0, -1, 1, &
-2, 1, 1, 3, -1, 1, 0, 1, & !F6 -1, 2, -1, -3, 1, -1, 0, 1, &
-1, 2, -1, 3, -1, 1, 0, 1, & !F4 -2, 1, 1, -3, 1, 0, -1, 1, &
-1, -1, 2, 3, -1, 0, 1, 1, & !G7 -1, -1, 2, -3, 0, 1, -1, 1, &
-2, 1, 1, 3, -1, 0, 1, 1, & !G6 1, -2, 1, -3, 0, -1, 1, 1, &
1, -2, 1, 3, 0, -1, 1, 1, & !H8 -2, 1, 1, -3, -1, 0, 1, 1, &
-1, -1, 2, 3, 0, -1, 1, 1, & !H7 -1, -1, 2, -3, 0, -1, 1, 1, &
2, -1, -1, 3, 1, -1, 0, 1, & !I9 1, -2, 1, -3, 0, 1, -1, 1, &
1, -2, 1, 3, 1, -1, 0, 1, & !I8 2, -1, -1, -3, -1, 1, 0, 1, &
1, 1, -2, 3, 1, 0, -1, 1, & !J5 1, 1, -2, -3, -1, 0, 1, 1, &
2, -1, -1, 3, 1, 0, -1, 1 & !J9 -1, 2, -1, -3, -1, 1, 0, 1 &
/),(/4+4,lattice_hex_Nslip/)) /),(/4+4,lattice_hex_Nslip/))
real(pReal), dimension(4+4,lattice_hex_Ntwin), parameter :: lattice_hex_systemTwin = & real(pReal), dimension(4+4,lattice_hex_Ntwin), parameter :: lattice_hex_systemTwin = &
@ -658,30 +661,6 @@ CONTAINS
!* - lattice_initializeStructure !* - lattice_initializeStructure
!**************************************** !****************************************
pure function lattice_symmetryType(structID)
!**************************************
!* maps structure to symmetry type *
!* fcc(1) and bcc(2) are cubic(1) *
!* hex(3+) is hexagonal(2) *
!**************************************
implicit none
integer(pInt), intent(in) :: structID
integer(pInt) lattice_symmetryType
select case(structID)
case (1,2)
lattice_symmetryType = 1_pInt
case (3:)
lattice_symmetryType = 2_pInt
case default
lattice_symmetryType = 0_pInt
end select
return
end function
subroutine lattice_init() subroutine lattice_init()
!************************************** !**************************************

6
code/material.config
View File

@ -230,9 +230,9 @@ rhoSglScrewPos0 1e11 0 0 0 # Initial positive screw single
rhoSglScrewNeg0 1e11 0 0 0 # Initial negative screw single dislocation density in m/m**3 rhoSglScrewNeg0 1e11 0 0 0 # Initial negative screw single dislocation density in m/m**3
rhoDipEdge0 1e8 0 0 0 # Initial edge dipole dislocation density in m/m**3 rhoDipEdge0 1e8 0 0 0 # Initial edge dipole dislocation density in m/m**3
rhoDipScrew0 1e8 0 0 0 # Initial screw dipole dislocation density in m/m**3 rhoDipScrew0 1e8 0 0 0 # Initial screw dipole dislocation density in m/m**3
r 1e-5 # cutoff radius for dislocation stress in m r 10e-6 # cutoff radius for dislocation stress in m
v0 1e-4 0 0 0 # prefactor for dislocation velocity v0 3500 0 0 0 # maximum dislocation velocity (velocity of sound)
Q0 3e-19 # activation energy for dislocation glide Q0 5e-19 # activation energy for dislocation glide
dDipMinEdge 1e-9 0 0 0 # minimum distance for stable edge dipoles in m dDipMinEdge 1e-9 0 0 0 # minimum distance for stable edge dipoles in m
dDipMinScrew 1e-9 0 0 0 # minimum distance for stable screw dipoles in m dDipMinScrew 1e-9 0 0 0 # minimum distance for stable screw dipoles in m
lambda0 100 0 0 0 # prefactor for mean free path lambda0 100 0 0 0 # prefactor for mean free path

2
code/numerics.config
View File

@ -7,6 +7,8 @@ iJacoStiffness 1 # frequency of stiffness update
iJacoLpresiduum 1 # frequency of Jacobian update of residuum in Lp iJacoLpresiduum 1 # frequency of Jacobian update of residuum in Lp
pert_Fg 1.0e-7 # deformation gradient perturbation for grain tangent pert_Fg 1.0e-7 # deformation gradient perturbation for grain tangent
pert_method 1 # perturbation method (1 = forward, 2 = backward or 3 = central) pert_method 1 # perturbation method (1 = forward, 2 = backward or 3 = central)
integrator 1 # integration method (1 = Fixed Point Iteration, 2 = Euler, 3 = Adaptive Euler, 4 = classical 4th order Runge-Kutta, 5 = 5th order Runge-Kutta Cash-Karp)
integratorStiffness 1 # integration method used for stiffness (1 = Fixed Point Iteration, 2 = Euler, 3 = Adaptive Euler, 4 = classical 4th order Runge-Kutta, 5 = 5th order Runge-Kutta Cash-Karp)
## crystallite numerical parameters ## ## crystallite numerical parameters ##
nCryst 20 # crystallite loop limit (only for debugging info, loop limit is determined by "subStepMinCryst") nCryst 20 # crystallite loop limit (only for debugging info, loop limit is determined by "subStepMinCryst")

16
code/numerics.f90
View File

@ -14,7 +14,9 @@ integer(pInt) iJacoStiffness, & ! freque
nCryst, & ! crystallite loop limit (only for debugging info, loop limit is determined by "subStepMinCryst") nCryst, & ! crystallite loop limit (only for debugging info, loop limit is determined by "subStepMinCryst")
nState, & ! state loop limit nState, & ! state loop limit
nStress, & ! stress loop limit nStress, & ! stress loop limit
pert_method ! method used in perturbation technique for tangent pert_method, & ! method used in perturbation technique for tangent
integrator, & ! method used for state integration
integratorStiffness ! method used for stiffness state integration
real(pReal) relevantStrain, & ! strain increment considered significant (used by crystallite to determine whether strain inc is considered significant) real(pReal) relevantStrain, & ! strain increment considered significant (used by crystallite to determine whether strain inc is considered significant)
defgradTolerance, & ! deviation of deformation gradient that is still allowed (used by CPFEM to determine outdated ffn1) defgradTolerance, & ! deviation of deformation gradient that is still allowed (used by CPFEM to determine outdated ffn1)
pert_Fg, & ! strain perturbation for FEM Jacobi pert_Fg, & ! strain perturbation for FEM Jacobi
@ -106,6 +108,8 @@ subroutine numerics_init()
rTol_crystalliteTemperature = 1.0e-6_pReal rTol_crystalliteTemperature = 1.0e-6_pReal
rTol_crystalliteStress = 1.0e-6_pReal rTol_crystalliteStress = 1.0e-6_pReal
aTol_crystalliteStress = 1.0e-8_pReal ! residuum is in Lp (hence strain on the order of 1e-8 here) aTol_crystalliteStress = 1.0e-8_pReal ! residuum is in Lp (hence strain on the order of 1e-8 here)
integrator = 1
integratorStiffness = 1
!* RGC parameters: added <<<updated 17.12.2009>>> with moderate setting !* RGC parameters: added <<<updated 17.12.2009>>> with moderate setting
absTol_RGC = 1.0e+4 absTol_RGC = 1.0e+4
@ -181,6 +185,10 @@ subroutine numerics_init()
rTol_crystalliteStress = IO_floatValue(line,positions,2) rTol_crystalliteStress = IO_floatValue(line,positions,2)
case ('atol_crystallitestress') case ('atol_crystallitestress')
aTol_crystalliteStress = IO_floatValue(line,positions,2) aTol_crystalliteStress = IO_floatValue(line,positions,2)
case ('integrator')
integrator = IO_intValue(line,positions,2)
case ('integratorstiffness')
integratorStiffness = IO_intValue(line,positions,2)
!* RGC parameters: added <<<updated 17.12.2009>>> !* RGC parameters: added <<<updated 17.12.2009>>>
case ('atol_rgc') case ('atol_rgc')
@ -242,6 +250,8 @@ subroutine numerics_init()
write(6,'(a24,x,e8.1)') 'rTol_crystalliteTemp: ',rTol_crystalliteTemperature write(6,'(a24,x,e8.1)') 'rTol_crystalliteTemp: ',rTol_crystalliteTemperature
write(6,'(a24,x,e8.1)') 'rTol_crystalliteStress: ',rTol_crystalliteStress write(6,'(a24,x,e8.1)') 'rTol_crystalliteStress: ',rTol_crystalliteStress
write(6,'(a24,x,e8.1)') 'aTol_crystalliteStress: ',aTol_crystalliteStress write(6,'(a24,x,e8.1)') 'aTol_crystalliteStress: ',aTol_crystalliteStress
write(6,'(a24,x,i8)') 'integrator: ',integrator
write(6,'(a24,x,i8)') 'integratorStiffness: ',integratorStiffness
write(6,*) write(6,*)
write(6,'(a24,x,i8)') 'nHomog: ',nHomog write(6,'(a24,x,i8)') 'nHomog: ',nHomog
@ -294,6 +304,10 @@ subroutine numerics_init()
if (rTol_crystalliteTemperature <= 0.0_pReal) call IO_error(276) !! oops !! if (rTol_crystalliteTemperature <= 0.0_pReal) call IO_error(276) !! oops !!
if (rTol_crystalliteStress <= 0.0_pReal) call IO_error(270) if (rTol_crystalliteStress <= 0.0_pReal) call IO_error(270)
if (aTol_crystalliteStress <= 0.0_pReal) call IO_error(271) if (aTol_crystalliteStress <= 0.0_pReal) call IO_error(271)
if (integrator <= 0_pInt .or. integrator >= 6_pInt) &
call IO_error(298)
if (integratorStiffness <= 0_pInt .or. integratorStiffness >= 6_pInt) &
call IO_error(298)
!* RGC parameters: added <<<updated 17.11.2009>>> !* RGC parameters: added <<<updated 17.11.2009>>>
if (absTol_RGC <= 0.0_pReal) call IO_error(272) if (absTol_RGC <= 0.0_pReal) call IO_error(272)