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

quite some changes: 

# non-greedy memory allocation
# generation of outputConstitutive to allow for script-based T16 extraction
# exchange of phenomenological by more general phenopowerlaw
# lattice is based on slip and twin families which can be treated as individual entities (switched on/off, separate hardening, etc.)
# nicer debugging output
# changed some error/warning codes
# plus potentially some minor additional brushes here and there
This commit is contained in:
Philip Eisenlohr 2009-07-22 16:07:19 +00:00
parent 290410b3fc
commit f337847f35
15 changed files with 2050 additions and 1782 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

44
code/CPFEM.f90
View File

@ -36,16 +36,16 @@ subroutine CPFEM_init()
! initialize stress and jacobian to zero ! initialize stress and jacobian to zero
allocate(CPFEM_cs(6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_cs = 0.0_pReal allocate(CPFEM_cs(6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_cs = 0.0_pReal
allocate(CPFEM_dcsdE(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dcsde = 0.0_pReal allocate(CPFEM_dcsdE(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dcsdE = 0.0_pReal
allocate(CPFEM_dcsdE_knownGood(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dcsde_knownGood = 0.0_pReal allocate(CPFEM_dcsdE_knownGood(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dcsdE_knownGood = 0.0_pReal
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,*) write(6,*)
write(6,*) '<<<+- cpfem init -+>>>' write(6,*) '<<<+- cpfem init -+>>>'
write(6,*) write(6,*)
write(6,'(a32,x,6(i5,x))') 'CPFEM_cs: ', shape(CPFEM_cs) write(6,'(a32,x,6(i5,x))') 'CPFEM_cs: ', shape(CPFEM_cs)
write(6,'(a32,x,6(i5,x))') 'CPFEM_dcsde: ', shape(CPFEM_dcsde) write(6,'(a32,x,6(i5,x))') 'CPFEM_dcsdE: ', shape(CPFEM_dcsdE)
write(6,'(a32,x,6(i5,x))') 'CPFEM_dcsde_knownGood: ', shape(CPFEM_dcsde_knownGood) write(6,'(a32,x,6(i5,x))') 'CPFEM_dcsdE_knownGood: ', shape(CPFEM_dcsdE_knownGood)
write(6,*) write(6,*)
write(6,*) 'parallelExecution: ', parallelExecution write(6,*) 'parallelExecution: ', parallelExecution
write(6,*) 'symmetricSolver: ', symmetricSolver write(6,*) 'symmetricSolver: ', symmetricSolver
@ -61,7 +61,7 @@ endsubroutine
!*** call the actual material model *** !*** call the actual material model ***
!*********************************************************************** !***********************************************************************
subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchyStress, jacobian, ngens) subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchyStress, jacobian, ngens)
! note: cauchyStress = Cauchy stress cs(6) and jacobian = Consistent tangent dcs/de ! note: cauchyStress = Cauchy stress cs(6) and jacobian = Consistent tangent dcs/dE
!*** variables and functions from other modules ***! !*** variables and functions from other modules ***!
use prec, only: pReal, & use prec, only: pReal, &
@ -90,12 +90,10 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
use mesh, only: mesh_init, & use mesh, only: mesh_init, &
mesh_FEasCP, & mesh_FEasCP, &
mesh_NcpElems, & mesh_NcpElems, &
mesh_maxNips, & mesh_maxNips
mesh_element
use lattice, only: lattice_init use lattice, only: lattice_init
use material, only: material_init, & use material, only: material_init, &
homogenization_maxNgrains, & homogenization_maxNgrains
homogenization_Ngrains
use constitutive, only: constitutive_init,& use constitutive, only: constitutive_init,&
constitutive_state0,constitutive_state constitutive_state0,constitutive_state
use crystallite, only: crystallite_init, & use crystallite, only: crystallite_init, &
@ -144,8 +142,6 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
real(pReal), dimension (3,3,3,3) :: H, & real(pReal), dimension (3,3,3,3) :: H, &
H_sym H_sym
integer(pInt) cp_en, & ! crystal plasticity element number integer(pInt) cp_en, & ! crystal plasticity element number
e, &
g, &
i, & i, &
j, & j, &
k, & k, &
@ -183,6 +179,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
cp_en = mesh_FEasCP('elem',element) cp_en = mesh_FEasCP('elem',element)
if (cp_en == 1 .and. IP == 1) then if (cp_en == 1 .and. IP == 1) then
write(6,*)
write(6,*) '#####################################' write(6,*) '#####################################'
write(6,'(a10,1x,f8.4,1x,a10,1x,i4,1x,a10,1x,i3,1x,a10,1x,i2,x,a10,1x,i2)') & write(6,'(a10,1x,f8.4,1x,a10,1x,i4,1x,a10,1x,i3,1x,a10,1x,i2,x,a10,1x,i2)') &
'theTime',theTime,'theInc',theInc,'theCycle',theCycle,'theLovl',theLovl,& 'theTime',theTime,'theInc',theInc,'theCycle',theCycle,'theLovl',theLovl,&
@ -200,18 +197,15 @@ 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
do e = 1,mesh_NcpElems forall ( i = 1:homogenization_maxNgrains, &
do g = 1,homogenization_Ngrains(mesh_element(3,e)) j = 1:mesh_maxNips, &
do i = 1,mesh_maxNips k = 1:mesh_NcpElems ) &
constitutive_state0(g,i,e)%p = constitutive_state(g,i,e)%p ! microstructure of crystallites constitutive_state0(i,j,k)%p = constitutive_state(i,j,k)%p ! microstructure of crystallites
enddo write(6,'(a10,/,4(3(f20.8,x),/))') 'aged state',constitutive_state(1,1,1)%p
enddo do k = 1,mesh_NcpElems
enddo do j = 1,mesh_maxNips
write(6,'(a10,/,4(3(f10.3,x),/))') 'aged state',constitutive_state(1,1,1)%p/1e6 if (homogenization_sizeState(j,k) > 0_pInt) &
do e = 1,mesh_NcpElems homogenization_state0(j,k)%p = homogenization_state(j,k)%p ! internal state of homogenization scheme
do i = 1,mesh_maxNips
if (homogenization_sizeState(i,e) > 0_pInt) &
homogenization_state0(i,e)%p = homogenization_state(i,e)%p ! internal state of homogenization scheme
enddo enddo
enddo enddo
endif endif
@ -269,7 +263,6 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
! --+>> COLLECTION OF FEM DATA AND RETURN OF ODD STRESS AND JACOBIAN <<+-- ! --+>> COLLECTION OF FEM DATA AND RETURN OF ODD STRESS AND JACOBIAN <<+--
case (3) case (3)
if (IP==1.AND.cp_en==1) write(6,*) 'Temp from CPFEM', Temperature
materialpoint_Temperature(IP,cp_en) = Temperature materialpoint_Temperature(IP,cp_en) = Temperature
materialpoint_F0(:,:,IP,cp_en) = ffn materialpoint_F0(:,:,IP,cp_en) = ffn
materialpoint_F(:,:,IP,cp_en) = ffn1 materialpoint_F(:,:,IP,cp_en) = ffn1
@ -294,14 +287,13 @@ subroutine CPFEM_general(mode, ffn, ffn1, Temperature, dt, element, IP, cauchySt
! return the local stress and the jacobian from storage ! return the local stress and the jacobian from storage
cauchyStress(1:ngens) = CPFEM_cs(1:ngens,IP,cp_en) cauchyStress(1:ngens) = CPFEM_cs(1:ngens,IP,cp_en)
jacobian(1:ngens,1:ngens) = CPFEM_dcsdE(1:ngens,1:ngens,IP,cp_en) jacobian(1:ngens,1:ngens) = CPFEM_dcsdE(1:ngens,1:ngens,IP,cp_en)
if (IP == 1 .and. cp_en == 1) write(6,'(a,/,6(6(f10.3,x)/))') 'jacobian/GPa at ip 1 el 1',jacobian/1e9
! return temperature ! return temperature
if (theInc > 0_pInt) Temperature = materialpoint_Temperature(IP,cp_en) ! homogenized result except for potentially non-isothermal starting condition. if (theInc > 0_pInt) Temperature = materialpoint_Temperature(IP,cp_en) ! homogenized result except for potentially non-isothermal starting condition.
return return
end subroutine end subroutine
END MODULE CPFEM END MODULE CPFEM

58
code/IO.f90
View File

@ -58,8 +58,9 @@
use prec, only: pReal, pInt use prec, only: pReal, pInt
implicit none implicit none
integer(pInt), intent(in) :: unit
integer(pInt) extPos
character(256) outName character(256) outName
integer(pInt) unit, extPos
character(3) ext character(3) ext
inquire(6, name=outName) ! determine outputfileName inquire(6, name=outName) ! determine outputfileName
@ -78,6 +79,29 @@
endfunction endfunction
!********************************************************************
! open (write) file related to current job
! but with different extension to given unit
!********************************************************************
logical function IO_open_jobFile(unit,newExt)
use prec, only: pReal, pInt
implicit none
integer(pInt), intent(in) :: unit
character(*), intent(in) :: newExt
character(256) outName
inquire(6, name=outName) ! determine outputfileName
open(unit,status='replace',err=100,file=outName(1:len_trim(outName)-3)//newExt)
IO_open_jobFile = .true.
return
100 IO_open_jobFile = .false.
return
endfunction
!******************************************************************** !********************************************************************
! hybrid IA repetition counter ! hybrid IA repetition counter
!******************************************************************** !********************************************************************
@ -771,6 +795,8 @@ endfunction
select case (ID) select case (ID)
case (0) case (0)
msg = 'Unable to open input file' msg = 'Unable to open input file'
case (50)
msg = 'Error writing constitutive output description'
case (100) case (100)
msg = 'Error reading from configuration file' msg = 'Error reading from configuration file'
case (105) case (105)
@ -793,22 +819,16 @@ endfunction
msg = 'Unknown constitution specified' msg = 'Unknown constitution specified'
case (201) case (201)
msg = 'Unknown homogenization specified' msg = 'Unknown homogenization specified'
case (202)
msg = 'Number of slip systems too small'
case (203)
msg = 'Negative initial slip resistance'
case (204)
msg = 'Non-positive reference shear rate'
case (205) case (205)
msg = 'Unknown lattice structure encountered'
case (210)
msg = 'Negative initial resistance'
case (211)
msg = 'Non-positive reference shear rate'
case (212)
msg = 'Non-positive stress exponent' msg = 'Non-positive stress exponent'
case (206) case (213)
msg = 'Non-positive initial hardening slope'
case (207)
msg = 'Non-positive saturation stress' msg = 'Non-positive saturation stress'
case (208)
msg = 'Non-positive w0'
case (209)
msg = 'Negative latent hardening ratio'
case (220) case (220)
msg = 'Negative initial dislocation density' msg = 'Negative initial dislocation density'
case (221) case (221)
@ -818,9 +838,13 @@ endfunction
case (223) case (223)
msg = 'Negative self diffusion energy' msg = 'Negative self diffusion energy'
case (224) case (224)
msg = 'Negative diffusion constant' msg = 'Non-positive diffusion prefactor'
case (225)
msg = 'No slip systems specified'
case (240) case (240)
msg = 'Non-positive Taylor factor' msg = 'Non-positive Taylor factor'
case (241)
msg = 'Non-positive hardening exponent'
case (260) case (260)
msg = 'Non-positive relevant strain' msg = 'Non-positive relevant strain'
case (261) case (261)
@ -853,8 +877,6 @@ endfunction
msg = 'Non-positive relative maximum value (upper bound) for GIA residual' msg = 'Non-positive relative maximum value (upper bound) for GIA residual'
case (275) case (275)
msg = 'Limit for GIA iteration too small' msg = 'Limit for GIA iteration too small'
case (276)
msg = 'Non-positive relative tolerance for temperature'
case (300) case (300)
msg = 'This material can only be used with elements with three direct stress components' msg = 'This material can only be used with elements with three direct stress components'
case (500) case (500)
@ -866,7 +888,7 @@ endfunction
case (700) case (700)
msg = 'Singular matrix in stress iteration' msg = 'Singular matrix in stress iteration'
case (800) case (800)
msg = 'GIA requires 8 grains per IP (bonehead, you!)' msg = 'RGC requires 8 grains per IP (bonehead, you!) -- but now outdated'
case default case default
msg = 'Unknown error number...' msg = 'Unknown error number...'
end select end select

109
code/constitutive.f90
View File

@ -39,24 +39,64 @@ subroutine constitutive_init()
!* Module initialization * !* Module initialization *
!************************************** !**************************************
use prec, only: pReal,pInt use prec, only: pReal,pInt
use IO, only: IO_error, IO_open_file use debug, only: debugger
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
use constitutive_phenomenological
use constitutive_j2 use constitutive_j2
use constitutive_phenopowerlaw
use constitutive_dislobased use constitutive_dislobased
integer(pInt), parameter :: fileunit = 200 integer(pInt), parameter :: fileunit = 200
integer(pInt) e,i,g,myInstance,myNgrains integer(pInt) e,i,g,p,myInstance,myNgrains
integer(pInt), dimension(:,:), pointer :: thisSize
character(64), dimension(:,:), pointer :: thisOutput
logical knownConstitution
if(.not. IO_open_file(fileunit,material_configFile)) call IO_error (100) ! corrupt config file if(.not. IO_open_file(fileunit,material_configFile)) call IO_error (100) ! corrupt config file
call constitutive_phenomenological_init(fileunit) ! parse all phases of this constitution call constitutive_j2_init(fileunit) ! parse all phases of this constitution
call constitutive_j2_init(fileunit) call constitutive_phenopowerlaw_init(fileunit)
call constitutive_dislobased_init(fileunit) call constitutive_dislobased_init(fileunit)
close(fileunit) close(fileunit)
! write description file for constitutive phase output
if(.not. IO_open_jobFile(fileunit,'outputConstitutive')) call IO_error (50) ! problems in writing file
do p = 1,material_Nphase
i = phase_constitutionInstance(p) ! which instance of a constitution is present phase
knownConstitution = .true. ! assume valid
select case(phase_constitution(p)) ! split per constitiution
case (constitutive_j2_label)
thisOutput => constitutive_j2_output
thisSize => constitutive_j2_sizePostResult
case (constitutive_phenopowerlaw_label)
thisOutput => constitutive_phenopowerlaw_output
thisSize => constitutive_phenopowerlaw_sizePostResult
case (constitutive_dislobased_label)
thisOutput => constitutive_dislobased_output
thisSize => constitutive_dislobased_sizePostResult
case default
knownConstitution = .false.
end select
write(fileunit,*)
write(fileunit,'(a)') '['//trim(phase_name(p))//']'
write(fileunit,*)
if (knownConstitution) then
write(fileunit,'(a)') '#'//char(9)//'constitution'//char(9)//trim(phase_constitution(p))
do e = 1,phase_Noutput(p)
write(fileunit,'(a,i4)') trim(thisOutput(e,i))//char(9),thisSize(e,i)
enddo
endif
enddo
close(fileunit)
! allocate memory for state management
allocate(constitutive_state0(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(constitutive_state0(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(constitutive_partionedState0(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(constitutive_partionedState0(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(constitutive_subState0(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems)) allocate(constitutive_subState0(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems))
@ -69,17 +109,9 @@ subroutine constitutive_init()
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = 1,FE_Nips(mesh_element(2,e)) ! loop over IPs do i = 1,FE_Nips(mesh_element(2,e)) ! loop over IPs
do g = 1,myNgrains ! loop over grains do g = 1,myNgrains ! loop over grains
debugger = (e == 1 .and. i == 1 .and. g == 1)
myInstance = phase_constitutionInstance(material_phase(g,i,e)) myInstance = phase_constitutionInstance(material_phase(g,i,e))
select case(phase_constitution(material_phase(g,i,e))) select case(phase_constitution(material_phase(g,i,e)))
case (constitutive_phenomenological_label)
allocate(constitutive_state0(g,i,e)%p(constitutive_phenomenological_sizeState(myInstance)))
allocate(constitutive_partionedState0(g,i,e)%p(constitutive_phenomenological_sizeState(myInstance)))
allocate(constitutive_subState0(g,i,e)%p(constitutive_phenomenological_sizeState(myInstance)))
allocate(constitutive_state(g,i,e)%p(constitutive_phenomenological_sizeState(myInstance)))
constitutive_state0(g,i,e)%p = constitutive_phenomenological_stateInit(myInstance)
constitutive_sizeDotState(g,i,e) = constitutive_phenomenological_sizeDotState(myInstance)
constitutive_sizeState(g,i,e) = constitutive_phenomenological_sizeState(myInstance)
constitutive_sizePostResults(g,i,e) = constitutive_phenomenological_sizePostResults(myInstance)
case (constitutive_j2_label) case (constitutive_j2_label)
allocate(constitutive_state0(g,i,e)%p(constitutive_j2_sizeState(myInstance))) allocate(constitutive_state0(g,i,e)%p(constitutive_j2_sizeState(myInstance)))
allocate(constitutive_partionedState0(g,i,e)%p(constitutive_j2_sizeState(myInstance))) allocate(constitutive_partionedState0(g,i,e)%p(constitutive_j2_sizeState(myInstance)))
@ -89,6 +121,15 @@ subroutine constitutive_init()
constitutive_sizeDotState(g,i,e) = constitutive_j2_sizeDotState(myInstance) constitutive_sizeDotState(g,i,e) = constitutive_j2_sizeDotState(myInstance)
constitutive_sizeState(g,i,e) = constitutive_j2_sizeState(myInstance) constitutive_sizeState(g,i,e) = constitutive_j2_sizeState(myInstance)
constitutive_sizePostResults(g,i,e) = constitutive_j2_sizePostResults(myInstance) constitutive_sizePostResults(g,i,e) = constitutive_j2_sizePostResults(myInstance)
case (constitutive_phenopowerlaw_label)
allocate(constitutive_state0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance)))
allocate(constitutive_partionedState0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance)))
allocate(constitutive_subState0(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance)))
allocate(constitutive_state(g,i,e)%p(constitutive_phenopowerlaw_sizeState(myInstance)))
constitutive_state0(g,i,e)%p = constitutive_phenopowerlaw_stateInit(myInstance)
constitutive_sizeDotState(g,i,e) = constitutive_phenopowerlaw_sizeDotState(myInstance)
constitutive_sizeState(g,i,e) = constitutive_phenopowerlaw_sizeState(myInstance)
constitutive_sizePostResults(g,i,e) = constitutive_phenopowerlaw_sizePostResults(myInstance)
case (constitutive_dislobased_label) case (constitutive_dislobased_label)
allocate(constitutive_state0(g,i,e)%p(constitutive_dislobased_sizeState(myInstance))) allocate(constitutive_state0(g,i,e)%p(constitutive_dislobased_sizeState(myInstance)))
allocate(constitutive_partionedState0(g,i,e)%p(constitutive_dislobased_sizeState(myInstance))) allocate(constitutive_partionedState0(g,i,e)%p(constitutive_dislobased_sizeState(myInstance)))
@ -139,8 +180,8 @@ function constitutive_homogenizedC(ipc,ip,el)
!********************************************************************* !*********************************************************************
use prec, only: pReal,pInt use prec, only: pReal,pInt
use material, only: phase_constitution,material_phase use material, only: phase_constitution,material_phase
use constitutive_phenomenological
use constitutive_j2 use constitutive_j2
use constitutive_phenopowerlaw
use constitutive_dislobased use constitutive_dislobased
implicit none implicit none
@ -149,13 +190,12 @@ function constitutive_homogenizedC(ipc,ip,el)
real(pReal), dimension(6,6) :: constitutive_homogenizedC real(pReal), dimension(6,6) :: constitutive_homogenizedC
select case (phase_constitution(material_phase(ipc,ip,el))) select case (phase_constitution(material_phase(ipc,ip,el)))
case (constitutive_phenomenological_label)
constitutive_homogenizedC = constitutive_phenomenological_homogenizedC(constitutive_state,ipc,ip,el)
case (constitutive_j2_label) case (constitutive_j2_label)
constitutive_homogenizedC = constitutive_j2_homogenizedC(constitutive_state,ipc,ip,el) constitutive_homogenizedC = constitutive_j2_homogenizedC(constitutive_state,ipc,ip,el)
case (constitutive_phenopowerlaw_label)
constitutive_homogenizedC = constitutive_phenopowerlaw_homogenizedC(constitutive_state,ipc,ip,el)
case (constitutive_dislobased_label) case (constitutive_dislobased_label)
constitutive_homogenizedC = constitutive_dislobased_homogenizedC(constitutive_state,ipc,ip,el) constitutive_homogenizedC = constitutive_dislobased_homogenizedC(constitutive_state,ipc,ip,el)
end select end select
return return
@ -174,8 +214,8 @@ subroutine constitutive_microstructure(Temperature,ipc,ip,el)
!********************************************************************* !*********************************************************************
use prec, only: pReal,pInt use prec, only: pReal,pInt
use material, only: phase_constitution,material_phase use material, only: phase_constitution,material_phase
use constitutive_phenomenological
use constitutive_j2 use constitutive_j2
use constitutive_phenopowerlaw
use constitutive_dislobased use constitutive_dislobased
implicit none implicit none
@ -184,13 +224,12 @@ integer(pInt) ipc,ip,el
real(pReal) Temperature real(pReal) Temperature
select case (phase_constitution(material_phase(ipc,ip,el))) select case (phase_constitution(material_phase(ipc,ip,el)))
case (constitutive_phenomenological_label)
call constitutive_phenomenological_microstructure(Temperature,constitutive_state,ipc,ip,el)
case (constitutive_j2_label) case (constitutive_j2_label)
call constitutive_j2_microstructure(Temperature,constitutive_state,ipc,ip,el) call constitutive_j2_microstructure(Temperature,constitutive_state,ipc,ip,el)
case (constitutive_phenopowerlaw_label)
call constitutive_phenopowerlaw_microstructure(Temperature,constitutive_state,ipc,ip,el)
case (constitutive_dislobased_label) case (constitutive_dislobased_label)
call constitutive_dislobased_microstructure(Temperature,constitutive_state,ipc,ip,el) call constitutive_dislobased_microstructure(Temperature,constitutive_state,ipc,ip,el)
end select end select
end subroutine end subroutine
@ -211,8 +250,8 @@ subroutine constitutive_LpAndItsTangent(Lp,dLp_dTstar, Tstar_v,Temperature,ipc,i
!********************************************************************* !*********************************************************************
use prec, only: pReal,pInt use prec, only: pReal,pInt
use material, only: phase_constitution,material_phase use material, only: phase_constitution,material_phase
use constitutive_phenomenological
use constitutive_j2 use constitutive_j2
use constitutive_phenopowerlaw
use constitutive_dislobased use constitutive_dislobased
implicit none implicit none
@ -224,13 +263,12 @@ subroutine constitutive_LpAndItsTangent(Lp,dLp_dTstar, Tstar_v,Temperature,ipc,i
real(pReal), dimension(9,9) :: dLp_dTstar real(pReal), dimension(9,9) :: dLp_dTstar
select case (phase_constitution(material_phase(ipc,ip,el))) select case (phase_constitution(material_phase(ipc,ip,el)))
case (constitutive_phenomenological_label)
call constitutive_phenomenological_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_j2_label) case (constitutive_j2_label)
call constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el) call constitutive_j2_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_phenopowerlaw_label)
call constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_dislobased_label) case (constitutive_dislobased_label)
call constitutive_dislobased_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el) call constitutive_dislobased_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,constitutive_state,ipc,ip,el)
end select end select
return return
@ -252,8 +290,8 @@ function constitutive_dotState(Tstar_v,Temperature,ipc,ip,el)
!********************************************************************* !*********************************************************************
use prec, only: pReal,pInt use prec, only: pReal,pInt
use material, only: phase_constitution,material_phase use material, only: phase_constitution,material_phase
use constitutive_phenomenological
use constitutive_j2 use constitutive_j2
use constitutive_phenopowerlaw
use constitutive_dislobased use constitutive_dislobased
implicit none implicit none
@ -264,13 +302,12 @@ function constitutive_dotState(Tstar_v,Temperature,ipc,ip,el)
real(pReal), dimension(constitutive_sizeDotState(ipc,ip,el)) :: constitutive_dotState real(pReal), dimension(constitutive_sizeDotState(ipc,ip,el)) :: constitutive_dotState
select case (phase_constitution(material_phase(ipc,ip,el))) select case (phase_constitution(material_phase(ipc,ip,el)))
case (constitutive_phenomenological_label)
constitutive_dotState = constitutive_phenomenological_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_j2_label) case (constitutive_j2_label)
constitutive_dotState = constitutive_j2_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el) constitutive_dotState = constitutive_j2_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_phenopowerlaw_label)
constitutive_dotState = constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_dislobased_label) case (constitutive_dislobased_label)
constitutive_dotState = constitutive_dislobased_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el) constitutive_dotState = constitutive_dislobased_dotState(Tstar_v,Temperature,constitutive_state,ipc,ip,el)
end select end select
return return
end function end function
@ -291,8 +328,8 @@ function constitutive_dotTemperature(Tstar_v,Temperature,ipc,ip,el)
!********************************************************************* !*********************************************************************
use prec, only: pReal,pInt use prec, only: pReal,pInt
use material, only: phase_constitution,material_phase use material, only: phase_constitution,material_phase
use constitutive_phenomenological
use constitutive_j2 use constitutive_j2
use constitutive_phenopowerlaw
use constitutive_dislobased use constitutive_dislobased
implicit none implicit none
@ -303,10 +340,10 @@ function constitutive_dotTemperature(Tstar_v,Temperature,ipc,ip,el)
real(pReal) constitutive_dotTemperature real(pReal) constitutive_dotTemperature
select case (phase_constitution(material_phase(ipc,ip,el))) select case (phase_constitution(material_phase(ipc,ip,el)))
case (constitutive_phenomenological_label)
constitutive_dotTemperature = constitutive_phenomenological_dotTemperature(Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_j2_label) case (constitutive_j2_label)
constitutive_dotTemperature = constitutive_j2_dotTemperature(Tstar_v,Temperature,constitutive_state,ipc,ip,el) constitutive_dotTemperature = constitutive_j2_dotTemperature(Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_phenopowerlaw_label)
constitutive_dotTemperature = constitutive_phenopowerlaw_dotTemperature(Tstar_v,Temperature,constitutive_state,ipc,ip,el)
case (constitutive_dislobased_label) case (constitutive_dislobased_label)
constitutive_dotTemperature = constitutive_dislobased_dotTemperature(Tstar_v,Temperature,constitutive_state,ipc,ip,el) constitutive_dotTemperature = constitutive_dislobased_dotTemperature(Tstar_v,Temperature,constitutive_state,ipc,ip,el)
@ -327,8 +364,8 @@ pure function constitutive_postResults(Tstar_v,Temperature,dt,ipc,ip,el)
!********************************************************************* !*********************************************************************
use prec, only: pReal,pInt use prec, only: pReal,pInt
use material, only: phase_constitution,material_phase use material, only: phase_constitution,material_phase
use constitutive_phenomenological
use constitutive_j2 use constitutive_j2
use constitutive_phenopowerlaw
use constitutive_dislobased use constitutive_dislobased
implicit none implicit none
@ -340,10 +377,10 @@ pure function constitutive_postResults(Tstar_v,Temperature,dt,ipc,ip,el)
constitutive_postResults = 0.0_pReal constitutive_postResults = 0.0_pReal
select case (phase_constitution(material_phase(ipc,ip,el))) select case (phase_constitution(material_phase(ipc,ip,el)))
case (constitutive_phenomenological_label)
constitutive_postResults = constitutive_phenomenological_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el)
case (constitutive_j2_label) case (constitutive_j2_label)
constitutive_postResults = constitutive_j2_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el) constitutive_postResults = constitutive_j2_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el)
case (constitutive_phenopowerlaw_label)
constitutive_postResults = constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el)
case (constitutive_dislobased_label) case (constitutive_dislobased_label)
constitutive_postResults = constitutive_dislobased_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el) constitutive_postResults = constitutive_dislobased_postResults(Tstar_v,Temperature,dt,constitutive_state,ipc,ip,el)

799
code/constitutive_dislobased.f90
View File

File diff suppressed because it is too large Load Diff

108
code/constitutive_j2.f90
View File

@ -18,7 +18,7 @@
! gdot0 0.001 ! gdot0 0.001
! n 20 ! n 20
! h0 75e6 ! h0 75e6
! tau_sat 63e6 ! tausat 63e6
! w0 2.25 ! w0 2.25
MODULE constitutive_j2 MODULE constitutive_j2
@ -32,7 +32,8 @@ MODULE constitutive_j2
integer(pInt), dimension(:), allocatable :: constitutive_j2_sizeDotState, & integer(pInt), dimension(:), allocatable :: constitutive_j2_sizeDotState, &
constitutive_j2_sizeState, & constitutive_j2_sizeState, &
constitutive_j2_sizePostResults constitutive_j2_sizePostResults
character(len=64), dimension(:,:), allocatable :: constitutive_j2_output integer(pInt), dimension(:,:), allocatable,target :: constitutive_j2_sizePostResult ! size of each post result output
character(len=64), dimension(:,:), allocatable,target :: constitutive_j2_output ! name of each post result output
real(pReal), dimension(:), allocatable :: constitutive_j2_C11 real(pReal), dimension(:), allocatable :: constitutive_j2_C11
real(pReal), dimension(:), allocatable :: constitutive_j2_C12 real(pReal), dimension(:), allocatable :: constitutive_j2_C12
real(pReal), dimension(:,:,:), allocatable :: constitutive_j2_Cslip_66 real(pReal), dimension(:,:,:), allocatable :: constitutive_j2_Cslip_66
@ -42,7 +43,7 @@ MODULE constitutive_j2
real(pReal), dimension(:), allocatable :: constitutive_j2_gdot0 real(pReal), dimension(:), allocatable :: constitutive_j2_gdot0
real(pReal), dimension(:), allocatable :: constitutive_j2_n real(pReal), dimension(:), allocatable :: constitutive_j2_n
real(pReal), dimension(:), allocatable :: constitutive_j2_h0 real(pReal), dimension(:), allocatable :: constitutive_j2_h0
real(pReal), dimension(:), allocatable :: constitutive_j2_tau_sat real(pReal), dimension(:), allocatable :: constitutive_j2_tausat
real(pReal), dimension(:), allocatable :: constitutive_j2_w0 real(pReal), dimension(:), allocatable :: constitutive_j2_w0
@ -69,16 +70,22 @@ subroutine constitutive_j2_init(file)
integer(pInt), intent(in) :: file integer(pInt), intent(in) :: file
integer(pInt), parameter :: maxNchunks = 7 integer(pInt), parameter :: maxNchunks = 7
integer(pInt), dimension(1+2*maxNchunks) :: positions integer(pInt), dimension(1+2*maxNchunks) :: positions
integer(pInt) section, maxNinstance, i,j,k,l, output integer(pInt) section, maxNinstance, i,j,k,l, output, mySize
character(len=64) tag character(len=64) tag
character(len=1024) line character(len=1024) line
write(6,*)
write(6,'(a20,a20,a12)') '<<<+- constitutive_',constitutive_j2_label,' init -+>>>'
write(6,*)
maxNinstance = count(phase_constitution == constitutive_j2_label) maxNinstance = count(phase_constitution == constitutive_j2_label)
if (maxNinstance == 0) return if (maxNinstance == 0) return
allocate(constitutive_j2_sizeDotState(maxNinstance)) ; constitutive_j2_sizeDotState = 0_pInt allocate(constitutive_j2_sizeDotState(maxNinstance)) ; constitutive_j2_sizeDotState = 0_pInt
allocate(constitutive_j2_sizeState(maxNinstance)) ; constitutive_j2_sizeState = 0_pInt allocate(constitutive_j2_sizeState(maxNinstance)) ; constitutive_j2_sizeState = 0_pInt
allocate(constitutive_j2_sizePostResults(maxNinstance)); constitutive_j2_sizePostResults = 0_pInt allocate(constitutive_j2_sizePostResults(maxNinstance)); constitutive_j2_sizePostResults = 0_pInt
allocate(constitutive_j2_sizePostResult(maxval(phase_Noutput), &
maxNinstance)) ; constitutive_j2_sizePostResult = 0_pInt
allocate(constitutive_j2_output(maxval(phase_Noutput), & allocate(constitutive_j2_output(maxval(phase_Noutput), &
maxNinstance)) ; constitutive_j2_output = '' maxNinstance)) ; constitutive_j2_output = ''
allocate(constitutive_j2_C11(maxNinstance)) ; constitutive_j2_C11 = 0.0_pReal allocate(constitutive_j2_C11(maxNinstance)) ; constitutive_j2_C11 = 0.0_pReal
@ -89,7 +96,7 @@ subroutine constitutive_j2_init(file)
allocate(constitutive_j2_gdot0(maxNinstance)) ; constitutive_j2_gdot0 = 0.0_pReal allocate(constitutive_j2_gdot0(maxNinstance)) ; constitutive_j2_gdot0 = 0.0_pReal
allocate(constitutive_j2_n(maxNinstance)) ; constitutive_j2_n = 0.0_pReal allocate(constitutive_j2_n(maxNinstance)) ; constitutive_j2_n = 0.0_pReal
allocate(constitutive_j2_h0(maxNinstance)) ; constitutive_j2_h0 = 0.0_pReal allocate(constitutive_j2_h0(maxNinstance)) ; constitutive_j2_h0 = 0.0_pReal
allocate(constitutive_j2_tau_sat(maxNinstance)) ; constitutive_j2_tau_sat = 0.0_pReal allocate(constitutive_j2_tausat(maxNinstance)) ; constitutive_j2_tausat = 0.0_pReal
allocate(constitutive_j2_w0(maxNinstance)) ; constitutive_j2_w0 = 0.0_pReal allocate(constitutive_j2_w0(maxNinstance)) ; constitutive_j2_w0 = 0.0_pReal
rewind(file) rewind(file)
@ -128,8 +135,8 @@ subroutine constitutive_j2_init(file)
constitutive_j2_n(i) = IO_floatValue(line,positions,2) constitutive_j2_n(i) = IO_floatValue(line,positions,2)
case ('h0') case ('h0')
constitutive_j2_h0(i) = IO_floatValue(line,positions,2) constitutive_j2_h0(i) = IO_floatValue(line,positions,2)
case ('tau_sat') case ('tausat')
constitutive_j2_tau_sat(i) = IO_floatValue(line,positions,2) constitutive_j2_tausat(i) = IO_floatValue(line,positions,2)
case ('w0') case ('w0')
constitutive_j2_w0(i) = IO_floatValue(line,positions,2) constitutive_j2_w0(i) = IO_floatValue(line,positions,2)
case ('taylorfactor') case ('taylorfactor')
@ -139,35 +146,40 @@ subroutine constitutive_j2_init(file)
enddo enddo
100 do i = 1,maxNinstance ! sanity checks 100 do i = 1,maxNinstance ! sanity checks
if (constitutive_j2_tau0(i) < 0.0_pReal) call IO_error(203) if (constitutive_j2_tau0(i) < 0.0_pReal) call IO_error(210)
if (constitutive_j2_gdot0(i) <= 0.0_pReal) call IO_error(204) if (constitutive_j2_gdot0(i) <= 0.0_pReal) call IO_error(211)
if (constitutive_j2_n(i) <= 0.0_pReal) call IO_error(205) if (constitutive_j2_n(i) <= 0.0_pReal) call IO_error(212)
if (constitutive_j2_h0(i) <= 0.0_pReal) call IO_error(206) if (constitutive_j2_tausat(i) <= 0.0_pReal) call IO_error(213)
if (constitutive_j2_tau_sat(i) <= 0.0_pReal) call IO_error(207) if (constitutive_j2_w0(i) <= 0.0_pReal) call IO_error(241)
if (constitutive_j2_w0(i) <= 0.0_pReal) call IO_error(208)
if (constitutive_j2_fTaylor(i) <= 0.0_pReal) call IO_error(240) if (constitutive_j2_fTaylor(i) <= 0.0_pReal) call IO_error(240)
enddo enddo
do i = 1,maxNinstance do i = 1,maxNinstance
do j = 1,maxval(phase_Noutput)
select case(constitutive_j2_output(j,i))
case('flowstress')
mySize = 1_pInt
case('strainrate')
mySize = 1_pInt
case default
mySize = 0_pInt
end select
if (mySize > 0_pInt) then ! any meaningful output found
constitutive_j2_sizePostResult(j,i) = mySize
constitutive_j2_sizePostResults(i) = &
constitutive_j2_sizePostResults(i) + mySize
endif
enddo
constitutive_j2_sizeDotState(i) = 1 constitutive_j2_sizeDotState(i) = 1
constitutive_j2_sizeState(i) = 1 constitutive_j2_sizeState(i) = 1
do j = 1,maxval(phase_Noutput)
select case(constitutive_j2_output(j,i))
case('flowstress')
constitutive_j2_sizePostResults(i) = &
constitutive_j2_sizePostResults(i) + 1
case('strainrate')
constitutive_j2_sizePostResults(i) = &
constitutive_j2_sizePostResults(i) + 1
end select
enddo
forall(k=1:3) forall(k=1:3)
forall(j=1:3) & forall(j=1:3) &
constitutive_j2_Cslip_66(k,j,i) = constitutive_j2_C12(i) constitutive_j2_Cslip_66(k,j,i) = constitutive_j2_C12(i)
constitutive_j2_Cslip_66(k,k,i) = constitutive_j2_C11(i) constitutive_j2_Cslip_66(k,k,i) = constitutive_j2_C11(i)
constitutive_j2_Cslip_66(k+3,k+3,i) = 0.5_pReal*(constitutive_j2_C11(i)-constitutive_j2_C12(i)) constitutive_j2_Cslip_66(k+3,k+3,i) = 0.5_pReal*(constitutive_j2_C11(i)-constitutive_j2_C12(i))
end forall end forall
constitutive_j2_Cslip_66(:,:,i) = & constitutive_j2_Cslip_66(:,:,i) = &
math_Mandel3333to66(math_Voigt66to3333(constitutive_j2_Cslip_66(:,:,i))) math_Mandel3333to66(math_Voigt66to3333(constitutive_j2_Cslip_66(:,:,i)))
@ -184,23 +196,12 @@ endsubroutine
!********************************************************************* !*********************************************************************
pure function constitutive_j2_stateInit(myInstance) pure function constitutive_j2_stateInit(myInstance)
!*** variables and functions from other modules ***!
use prec, only: pReal,pInt use prec, only: pReal,pInt
implicit none implicit none
!*** input variables ***!
integer(pInt), intent(in) :: myInstance integer(pInt), intent(in) :: myInstance
!*** output variables ***!
real(pReal), dimension(1) :: constitutive_j2_stateInit real(pReal), dimension(1) :: constitutive_j2_stateInit
!*** local variables ***!
!*** global variables ***!
! constitutive_j2_tau0
constitutive_j2_stateInit = constitutive_j2_tau0(myInstance) constitutive_j2_stateInit = constitutive_j2_tau0(myInstance)
return return
@ -221,7 +222,6 @@ function constitutive_j2_homogenizedC(state,ipc,ip,el)
use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance
implicit none implicit none
!* Definition of variables
integer(pInt), intent(in) :: ipc,ip,el integer(pInt), intent(in) :: ipc,ip,el
integer(pInt) matID integer(pInt) matID
real(pReal), dimension(6,6) :: constitutive_j2_homogenizedC real(pReal), dimension(6,6) :: constitutive_j2_homogenizedC
@ -306,12 +306,6 @@ pure subroutine constitutive_j2_LpAndItsTangent(Lp, dLp_dTstar_99, Tstar_dev_v,
m, & m, &
n n
!*** global variables ***!
! constitutive_j2_gdot0
! constitutive_j2_fTaylor
! constitutive_j2_n
matID = phase_constitutionInstance(material_phase(g,ip,el)) matID = phase_constitutionInstance(material_phase(g,ip,el))
! convert Tstar to matrix and calculate euclidean norm ! convert Tstar to matrix and calculate euclidean norm
@ -384,32 +378,24 @@ pure function constitutive_j2_dotState(Tstar_v, Temperature, state, g, ip, el)
norm_Tstar_dev ! euclidean norm of Tstar_dev norm_Tstar_dev ! euclidean norm of Tstar_dev
integer(pInt) matID integer(pInt) matID
!*** global variables ***!
! constitutive_j2_gdot0
! constitutive_j2_fTaylor
! constitutive_j2_n
! constitutive_j2_h0
! constitutive_j2_tau_sat
! constitutive_j2_w0
matID = phase_constitutionInstance(material_phase(g,ip,el)) matID = phase_constitutionInstance(material_phase(g,ip,el))
! caclulate deviatoric part of 2nd Piola-Kirchhoff stress ! deviatoric part of 2nd Piola-Kirchhoff stress
Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3_pReal 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) Tstar_dev_v(4:6) = Tstar_v(4:6)
norm_Tstar_dev = dsqrt(math_mul6x6(Tstar_dev_v,Tstar_dev_v)) norm_Tstar_dev = dsqrt(math_mul6x6(Tstar_dev_v,Tstar_dev_v))
! Calculation of gamma_dot ! gamma_dot
gamma_dot = constitutive_j2_gdot0(matID) * ( dsqrt(1.5_pReal) * norm_Tstar_dev & gamma_dot = constitutive_j2_gdot0(matID) * ( dsqrt(1.5_pReal) * norm_Tstar_dev &
/ &!--------------------------------------------------- / &!---------------------------------------------------
(constitutive_j2_fTaylor(matID) * state(g,ip,el)%p(1)) ) ** constitutive_j2_n(matID) (constitutive_j2_fTaylor(matID) * state(g,ip,el)%p(1)) ) ** constitutive_j2_n(matID)
! calculate hardening coefficient ! hardening coefficient
hardening = constitutive_j2_h0(matID) * & hardening = constitutive_j2_h0(matID) * &
( 1.0_pReal - state(g,ip,el)%p(1) / constitutive_j2_tau_sat(matID) ) ** constitutive_j2_w0(matID) ( 1.0_pReal - state(g,ip,el)%p(1) / constitutive_j2_tausat(matID) ) ** constitutive_j2_w0(matID)
! finally calculate dotState ! dotState
constitutive_j2_dotState = hardening * gamma_dot constitutive_j2_dotState = hardening * gamma_dot
return return
@ -418,7 +404,7 @@ endfunction
!**************************************************************** !****************************************************************
!* calculates the rate of change of microstructure * !* calculates the rate of change of temperature *
!**************************************************************** !****************************************************************
pure function constitutive_j2_dotTemperature(Tstar_v, Temperature, state, g, ip, el) pure function constitutive_j2_dotTemperature(Tstar_v, Temperature, state, g, ip, el)
@ -438,7 +424,7 @@ pure function constitutive_j2_dotTemperature(Tstar_v, Temperature, state, g, ip,
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state ! state of the current microstructure type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state ! state of the current microstructure
!*** output variables ***! !*** output variables ***!
real(pReal) constitutive_j2_dotTemperature ! evolution of state variable real(pReal) constitutive_j2_dotTemperature ! rate of change of temperature
! calculate dotTemperature ! calculate dotTemperature
constitutive_j2_dotTemperature = 0.0_pReal constitutive_j2_dotTemperature = 0.0_pReal

519
code/constitutive_phenomenological.f90
View File

@ -1,519 +0,0 @@
!*****************************************************
!* Module: CONSTITUTIVE_PHENOMENOLOGICAL *
!*****************************************************
!* contains: *
!* - constitutive equations *
!* - parameters definition *
!*****************************************************
! [Alu]
! constitution phenomenological
! (output) slipresistance
! (output) rateofshear
! lattice_structure 1
! Nslip 12
!
! c11 106.75e9
! c12 60.41e9
! c44 28.34e9
!
! s0_slip 31e6
! gdot0_slip 0.001
! n_slip 20
! h0 75e6
! s_sat 63e6
! w0 2.25
! latent_ratio 1.4
MODULE constitutive_phenomenological
!*** Include other modules ***
use prec, only: pReal,pInt
implicit none
character (len=*), parameter :: constitutive_phenomenological_label = 'phenomenological'
integer(pInt), dimension(:), allocatable :: constitutive_phenomenological_sizeDotState, &
constitutive_phenomenological_sizeState, &
constitutive_phenomenological_sizePostResults
character(len=64), dimension(:,:), allocatable :: constitutive_phenomenological_output
character(len=32), dimension(:), allocatable :: constitutive_phenomenological_structureName
integer(pInt), dimension(:), allocatable :: constitutive_phenomenological_structure
integer(pInt), dimension(:), allocatable :: constitutive_phenomenological_Nslip
real(pReal), dimension(:), allocatable :: constitutive_phenomenological_CoverA
real(pReal), dimension(:), allocatable :: constitutive_phenomenological_C11
real(pReal), dimension(:), allocatable :: constitutive_phenomenological_C12
real(pReal), dimension(:), allocatable :: constitutive_phenomenological_C13
real(pReal), dimension(:), allocatable :: constitutive_phenomenological_C33
real(pReal), dimension(:), allocatable :: constitutive_phenomenological_C44
real(pReal), dimension(:,:,:), allocatable :: constitutive_phenomenological_Cslip_66
!* Visco-plastic constitutive_phenomenological parameters
real(pReal), dimension(:), allocatable :: constitutive_phenomenological_s0_slip
real(pReal), dimension(:), allocatable :: constitutive_phenomenological_gdot0_slip
real(pReal), dimension(:), allocatable :: constitutive_phenomenological_n_slip
real(pReal), dimension(:), allocatable :: constitutive_phenomenological_h0
real(pReal), dimension(:), allocatable :: constitutive_phenomenological_s_sat
real(pReal), dimension(:), allocatable :: constitutive_phenomenological_w0
real(pReal), dimension(:), allocatable :: constitutive_phenomenological_latent
real(pReal), dimension(:,:,:), allocatable :: constitutive_phenomenological_HardeningMatrix
CONTAINS
!****************************************
!* - constitutive_init
!* - constitutive_stateInit
!* - constitutive_homogenizedC
!* - constitutive_microstructure
!* - constitutive_LpAndItsTangent
!* - consistutive_dotState
!* - consistutive_postResults
!****************************************
subroutine constitutive_phenomenological_init(file)
!**************************************
!* Module initialization *
!**************************************
use prec, only: pInt, pReal
use math, only: math_Mandel3333to66, math_Voigt66to3333
use IO
use material
use lattice, only: lattice_initializeStructure
integer(pInt), intent(in) :: file
integer(pInt), parameter :: maxNchunks = 7
integer(pInt), dimension(1+2*maxNchunks) :: positions
integer(pInt) section, maxNinstance, i,j,k, output
character(len=64) tag
character(len=1024) line
maxNinstance = count(phase_constitution == constitutive_phenomenological_label)
if (maxNinstance == 0) return
allocate(constitutive_phenomenological_sizeDotState(maxNinstance)) ; constitutive_phenomenological_sizeDotState = 0_pInt
allocate(constitutive_phenomenological_sizeState(maxNinstance)) ; constitutive_phenomenological_sizeState = 0_pInt
allocate(constitutive_phenomenological_sizePostResults(maxNinstance)); constitutive_phenomenological_sizePostResults = 0_pInt
allocate(constitutive_phenomenological_output(maxval(phase_Noutput), &
maxNinstance)) ; constitutive_phenomenological_output = ''
allocate(constitutive_phenomenological_structureName(maxNinstance)) ; constitutive_phenomenological_structureName = ''
allocate(constitutive_phenomenological_structure(maxNinstance)) ; constitutive_phenomenological_structure = 0_pInt
allocate(constitutive_phenomenological_Nslip(maxNinstance)) ; constitutive_phenomenological_Nslip = 0_pInt
allocate(constitutive_phenomenological_CoverA(maxNinstance)) ; constitutive_phenomenological_CoverA = 0.0_pReal
allocate(constitutive_phenomenological_C11(maxNinstance)) ; constitutive_phenomenological_C11 = 0.0_pReal
allocate(constitutive_phenomenological_C12(maxNinstance)) ; constitutive_phenomenological_C12 = 0.0_pReal
allocate(constitutive_phenomenological_C13(maxNinstance)) ; constitutive_phenomenological_C13 = 0.0_pReal
allocate(constitutive_phenomenological_C33(maxNinstance)) ; constitutive_phenomenological_C33 = 0.0_pReal
allocate(constitutive_phenomenological_C44(maxNinstance)) ; constitutive_phenomenological_C44 = 0.0_pReal
allocate(constitutive_phenomenological_Cslip_66(6,6,maxNinstance)) ; constitutive_phenomenological_Cslip_66 = 0.0_pReal
allocate(constitutive_phenomenological_s0_slip(maxNinstance)) ; constitutive_phenomenological_s0_slip = 0.0_pReal
allocate(constitutive_phenomenological_gdot0_slip(maxNinstance)) ; constitutive_phenomenological_gdot0_slip = 0.0_pReal
allocate(constitutive_phenomenological_n_slip(maxNinstance)) ; constitutive_phenomenological_n_slip = 0.0_pReal
allocate(constitutive_phenomenological_h0(maxNinstance)) ; constitutive_phenomenological_h0 = 0.0_pReal
allocate(constitutive_phenomenological_s_sat(maxNinstance)) ; constitutive_phenomenological_s_sat = 0.0_pReal
allocate(constitutive_phenomenological_w0(maxNinstance)) ; constitutive_phenomenological_w0 = 0.0_pReal
allocate(constitutive_phenomenological_latent(maxNinstance)) ; constitutive_phenomenological_latent = 0.0_pReal
rewind(file)
line = ''
section = 0
do while (IO_lc(IO_getTag(line,'<','>')) /= 'phase') ! wind forward to <phase>
read(file,'(a1024)',END=100) line
enddo
do ! read thru sections of phase part
read(file,'(a1024)',END=100) line
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1
output = 0 ! reset output counter
endif
if (section > 0 .and. phase_constitution(section) == constitutive_phenomenological_label) then ! one of my sections
i = phase_constitutionInstance(section) ! which instance of my constitution is present phase
positions = IO_stringPos(line,maxNchunks)
tag = IO_lc(IO_stringValue(line,positions,1)) ! extract key
select case(tag)
case ('(output)')
output = output + 1
constitutive_phenomenological_output(output,i) = IO_lc(IO_stringValue(line,positions,2))
case ('lattice_structure')
constitutive_phenomenological_structureName(i) = IO_lc(IO_stringValue(line,positions,2))
case ('nslip')
constitutive_phenomenological_Nslip(i) = IO_intValue(line,positions,2)
case ('covera_ratio')
constitutive_phenomenological_CoverA(i) = IO_floatValue(line,positions,2)
case ('c11')
constitutive_phenomenological_C11(i) = IO_floatValue(line,positions,2)
case ('c12')
constitutive_phenomenological_C12(i) = IO_floatValue(line,positions,2)
case ('c13')
constitutive_phenomenological_C13(i) = IO_floatValue(line,positions,2)
case ('c33')
constitutive_phenomenological_C33(i) = IO_floatValue(line,positions,2)
case ('c44')
constitutive_phenomenological_C44(i) = IO_floatValue(line,positions,2)
case ('s0_slip')
constitutive_phenomenological_s0_slip(i) = IO_floatValue(line,positions,2)
case ('gdot0_slip')
constitutive_phenomenological_gdot0_slip(i) = IO_floatValue(line,positions,2)
case ('n_slip')
constitutive_phenomenological_n_slip(i) = IO_floatValue(line,positions,2)
case ('h0')
constitutive_phenomenological_h0(i) = IO_floatValue(line,positions,2)
case ('s_sat')
constitutive_phenomenological_s_sat(i) = IO_floatValue(line,positions,2)
case ('w0')
constitutive_phenomenological_w0(i) = IO_floatValue(line,positions,2)
case ('latent_ratio')
constitutive_phenomenological_latent(i) = IO_floatValue(line,positions,2)
end select
endif
enddo
100 do i = 1,maxNinstance
constitutive_phenomenological_structure(i) = lattice_initializeStructure(constitutive_phenomenological_structureName(i), &
constitutive_phenomenological_CoverA(i)) ! sanity checks
if (constitutive_phenomenological_structure(i) < 1 .or. &
constitutive_phenomenological_structure(i) > 3) call IO_error(201)
if (constitutive_phenomenological_Nslip(i) < 1) call IO_error(202)
if (constitutive_phenomenological_s0_slip(i) < 0.0_pReal) call IO_error(203)
if (constitutive_phenomenological_gdot0_slip(i) <= 0.0_pReal) call IO_error(204)
if (constitutive_phenomenological_n_slip(i) <= 0.0_pReal) call IO_error(205)
if (constitutive_phenomenological_h0(i) <= 0.0_pReal) call IO_error(206)
if (constitutive_phenomenological_s_sat(i) <= 0.0_pReal) call IO_error(207)
if (constitutive_phenomenological_w0(i) <= 0.0_pReal) call IO_error(208)
if (constitutive_phenomenological_latent(i) <= 0.0_pReal) call IO_error(209)
enddo
allocate(constitutive_phenomenological_hardeningMatrix(maxval(constitutive_phenomenological_Nslip),&
maxval(constitutive_phenomenological_Nslip),&
maxNinstance))
do i = 1,maxNinstance
constitutive_phenomenological_sizeDotState(i) = constitutive_phenomenological_Nslip(i)
constitutive_phenomenological_sizeState(i) = constitutive_phenomenological_Nslip(i)
do j = 1,maxval(phase_Noutput)
select case(constitutive_phenomenological_output(j,i))
case('slipresistance')
constitutive_phenomenological_sizePostResults(i) = &
constitutive_phenomenological_sizePostResults(i) + constitutive_phenomenological_Nslip(i)
case('rateofshear')
constitutive_phenomenological_sizePostResults(i) = &
constitutive_phenomenological_sizePostResults(i) + constitutive_phenomenological_Nslip(i)
end select
enddo
select case (constitutive_phenomenological_structure(i))
case(1:2) ! cubic(s)
forall(k=1:3)
forall(j=1:3) &
constitutive_phenomenological_Cslip_66(k,j,i) = constitutive_phenomenological_C12(i)
constitutive_phenomenological_Cslip_66(k,k,i) = constitutive_phenomenological_C11(i)
constitutive_phenomenological_Cslip_66(k+3,k+3,i) = constitutive_phenomenological_C44(i)
end forall
case(3) ! hcp
constitutive_phenomenological_Cslip_66(1,1,i) = constitutive_phenomenological_C11(i)
constitutive_phenomenological_Cslip_66(2,2,i) = constitutive_phenomenological_C11(i)
constitutive_phenomenological_Cslip_66(3,3,i) = constitutive_phenomenological_C33(i)
constitutive_phenomenological_Cslip_66(1,2,i) = constitutive_phenomenological_C12(i)
constitutive_phenomenological_Cslip_66(2,1,i) = constitutive_phenomenological_C12(i)
constitutive_phenomenological_Cslip_66(1,3,i) = constitutive_phenomenological_C13(i)
constitutive_phenomenological_Cslip_66(3,1,i) = constitutive_phenomenological_C13(i)
constitutive_phenomenological_Cslip_66(2,3,i) = constitutive_phenomenological_C13(i)
constitutive_phenomenological_Cslip_66(3,2,i) = constitutive_phenomenological_C13(i)
constitutive_phenomenological_Cslip_66(4,4,i) = constitutive_phenomenological_C44(i)
constitutive_phenomenological_Cslip_66(5,5,i) = constitutive_phenomenological_C44(i)
constitutive_phenomenological_Cslip_66(6,6,i) = 0.5_pReal*(constitutive_phenomenological_C11(i)- &
constitutive_phenomenological_C12(i))
end select
constitutive_phenomenological_Cslip_66(:,:,i) = &
math_Mandel3333to66(math_Voigt66to3333(constitutive_phenomenological_Cslip_66(:,:,i)))
constitutive_phenomenological_hardeningMatrix(:,:,i) = constitutive_phenomenological_latent(i)
forall (j = 1:constitutive_phenomenological_Nslip(i)) &
constitutive_phenomenological_hardeningMatrix(j,j,i) = 1.0_pReal
enddo
return
end subroutine
function constitutive_phenomenological_stateInit(myInstance)
!*********************************************************************
!* initial microstructural state *
!*********************************************************************
use prec, only: pReal,pInt
implicit none
!* Definition of variables
integer(pInt), intent(in) :: myInstance
real(pReal), dimension(constitutive_phenomenological_Nslip(myInstance)) :: constitutive_phenomenological_stateInit
constitutive_phenomenological_stateInit = constitutive_phenomenological_s0_slip(myInstance)
return
end function
function constitutive_phenomenological_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: pReal,pInt,p_vec
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance
implicit none
!* Definition of variables
integer(pInt), intent(in) :: ipc,ip,el
integer(pInt) matID
real(pReal), dimension(6,6) :: constitutive_phenomenological_homogenizedC
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state
matID = phase_constitutionInstance(material_phase(ipc,ip,el))
constitutive_phenomenological_homogenizedC = constitutive_phenomenological_Cslip_66(:,:,matID)
return
end function
subroutine constitutive_phenomenological_microstructure(Temperature,state,ipc,ip,el)
!*********************************************************************
!* calculate derived quantities from state (not used here) *
!* INPUT: *
!* - Tp : temperature *
!* - ipc : component-ID of current integration point *
!* - ip : current integration point *
!* - el : current element *
!*********************************************************************
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance
implicit none
!* Definition of variables
integer(pInt) ipc,ip,el, matID
real(pReal) Temperature
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state
matID = phase_constitutionInstance(material_phase(ipc,ip,el))
end subroutine
subroutine constitutive_phenomenological_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,state,ipc,ip,el)
!*********************************************************************
!* plastic velocity gradient and its tangent *
!* INPUT: *
!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) *
!* - ipc : component-ID at current integration point *
!* - ip : current integration point *
!* - el : current element *
!* OUTPUT: *
!* - Lp : plastic velocity gradient *
!* - dLp_dTstar : derivative of Lp (4th-rank tensor) *
!*********************************************************************
use prec, only: pReal,pInt,p_vec
use math, only: math_Plain3333to99
use lattice, only: lattice_Sslip,lattice_Sslip_v
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance
implicit none
!* Definition of variables
integer(pInt) ipc,ip,el
integer(pInt) matID,i,k,l,m,n
real(pReal) Temperature
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state
real(pReal), dimension(6) :: Tstar_v
real(pReal), dimension(3,3) :: Lp
real(pReal), dimension(3,3,3,3) :: dLp_dTstar3333
real(pReal), dimension(9,9) :: dLp_dTstar
real(pReal), dimension(constitutive_phenomenological_Nslip(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: &
gdot_slip,dgdot_dtauslip,tau_slip
matID = phase_constitutionInstance(material_phase(ipc,ip,el))
!* Calculation of Lp
Lp = 0.0_pReal
do i = 1,constitutive_phenomenological_Nslip(matID)
tau_slip(i) = dot_product(Tstar_v,lattice_Sslip_v(:,i,constitutive_phenomenological_structure(matID)))
gdot_slip(i) = constitutive_phenomenological_gdot0_slip(matID)*(abs(tau_slip(i))/state(ipc,ip,el)%p(i))**&
constitutive_phenomenological_n_slip(matID)*sign(1.0_pReal,tau_slip(i))
Lp = Lp + gdot_slip(i)*lattice_Sslip(:,:,i,constitutive_phenomenological_structure(matID))
enddo
!* Calculation of the tangent of Lp
dLp_dTstar3333 = 0.0_pReal
dLp_dTstar = 0.0_pReal
do i = 1,constitutive_phenomenological_Nslip(matID)
dgdot_dtauslip(i) = constitutive_phenomenological_gdot0_slip(matID)*(abs(tau_slip(i))/state(ipc,ip,el)%p(i))**&
(constitutive_phenomenological_n_slip(matID)-1.0_pReal)*&
constitutive_phenomenological_n_slip(matID)/state(ipc,ip,el)%p(i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
dgdot_dtauslip(i)*lattice_Sslip(k,l,i,constitutive_phenomenological_structure(matID))* &
lattice_Sslip(m,n,i,constitutive_phenomenological_structure(matID))
enddo
dLp_dTstar = math_Plain3333to99(dLp_dTstar3333)
return
end subroutine
function constitutive_phenomenological_dotState(Tstar_v,Temperature,state,ipc,ip,el)
!*********************************************************************
!* rate of change of microstructure *
!* INPUT: *
!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) *
!* - ipc : component-ID at current integration point *
!* - ip : current integration point *
!* - el : current element *
!* OUTPUT: *
!* - constitutive_dotState : evolution of state variable *
!*********************************************************************
use prec, only: pReal,pInt,p_vec
use lattice, only: lattice_Sslip_v
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance
implicit none
!* Definition of variables
integer(pInt) ipc,ip,el
integer(pInt) matID,i,n
real(pReal) Temperature,tau_slip,gdot_slip
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state
real(pReal), dimension(6) :: Tstar_v
real(pReal), dimension(constitutive_phenomenological_Nslip(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: &
constitutive_phenomenological_dotState,self_hardening
matID = phase_constitutionInstance(material_phase(ipc,ip,el))
n = constitutive_phenomenological_Nslip(matID)
!* Self-Hardening of each system
do i = 1,n
tau_slip = dot_product(Tstar_v,lattice_Sslip_v(:,i,constitutive_phenomenological_structure(matID)))
gdot_slip = constitutive_phenomenological_gdot0_slip(matID)*(abs(tau_slip)/state(ipc,ip,el)%p(i))**&
constitutive_phenomenological_n_slip(matID)*sign(1.0_pReal,tau_slip)
self_hardening(i) = constitutive_phenomenological_h0(matID)*(1.0_pReal-state(ipc,ip,el)%p(i)/&
constitutive_phenomenological_s_sat(matID))**constitutive_phenomenological_w0(matID)*abs(gdot_slip)
enddo
!$OMP CRITICAL (evilmatmul)
constitutive_phenomenological_dotState = matmul(constitutive_phenomenological_hardeningMatrix(1:n,1:n,matID),self_hardening)
!$OMP END CRITICAL (evilmatmul)
return
end function
function constitutive_phenomenological_dotTemperature(Tstar_v,Temperature,state,ipc,ip,el)
!*********************************************************************
!* rate of change of microstructure *
!* INPUT: *
!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) *
!* - ipc : component-ID at current integration point *
!* - ip : current integration point *
!* - el : current element *
!* OUTPUT: *
!* - constitutive_dotTemperature : evolution of temperature *
!*********************************************************************
use prec, only: pReal,pInt,p_vec
use lattice, only: lattice_Sslip_v
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance
implicit none
!* Definition of variables
integer(pInt) ipc,ip,el
integer(pInt) matID,i,n
real(pReal) Temperature
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state
real(pReal), dimension(6) :: Tstar_v
real(pReal) constitutive_phenomenological_dotTemperature
constitutive_phenomenological_dotTemperature = 0.0_pReal
return
end function
pure function constitutive_phenomenological_postResults(Tstar_v,Temperature,dt,state,ipc,ip,el)
!*********************************************************************
!* return array of constitutive results *
!* INPUT: *
!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) *
!* - dt : current time increment *
!* - ipc : component-ID at current integration point *
!* - ip : current integration point *
!* - el : current element *
!*********************************************************************
use prec, only: pReal,pInt,p_vec
use lattice, only: lattice_Sslip_v
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance,phase_Noutput
implicit none
!* Definition of variables
integer(pInt), intent(in) :: ipc,ip,el
real(pReal), intent(in) :: dt,Temperature
real(pReal), dimension(6), intent(in) :: Tstar_v
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state
integer(pInt) matID,o,i,c,n
real(pReal) tau_slip
real(pReal), dimension(constitutive_phenomenological_sizePostResults(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: &
constitutive_phenomenological_postResults
matID = phase_constitutionInstance(material_phase(ipc,ip,el))
n = constitutive_phenomenological_Nslip(matID)
c = 0_pInt
constitutive_phenomenological_postResults = 0.0_pReal
do o = 1,phase_Noutput(material_phase(ipc,ip,el))
select case(constitutive_phenomenological_output(o,matID))
case ('slipresistance')
constitutive_phenomenological_postResults(c+1:c+n) = state(ipc,ip,el)%p(1:n)
c = c + n
case ('rateofshear')
do i = 1,n
tau_slip = dot_product(Tstar_v,lattice_Sslip_v(:,i,constitutive_phenomenological_structure(matID)))
constitutive_phenomenological_postResults(c+i) = sign(1.0_pReal,tau_slip)*constitutive_phenomenological_gdot0_slip(matID)*&
(abs(tau_slip)/state(ipc,ip,el)%p(i))**&
constitutive_phenomenological_n_slip(matID)
enddo
c = c + n
end select
enddo
return
end function
END MODULE

912
code/constitutive_phenopowerlaw.f90 Normal file
View File

@ -0,0 +1,912 @@
!*****************************************************
!* Module: CONSTITUTIVE_PHENOPOWERLAW *
!*****************************************************
!* contains: *
!* - constitutive equations *
!* - parameters definition *
!*****************************************************
! [Alu]
! constitution phenopowerlaw
! (output) resistance_slip
! (output) shearrate_slip
! (output) resolvedstress_slip
! (output) totalshear
! (output) resistance_twin
! (output) shearrate_twin
! (output) resolvedstress_twin
! (output) totalvolfrac
! lattice_structure hex
! covera_ratio 1.57
! Nslip 3 3 6 12 # per family
! Ntwin 6 6 6 6 # per family
!
! c11 106.75e9
! c12 60.41e9
! c44 28.34e9
!
! gdot0_slip 0.001
! n_slip 20
! tau0_slip 31e6 31e6 60e6 123e6 # per family
! tausat_slip 63e6 90e6 200e6 400e6 # per family
! gdot0_twin 0.001
! n_twin 20
! tau0_twin 31e6 31e6 60e6 123e6 # per family
! s_pr 100e6 # push-up factor for slip saturation due to twinning
! twin_b 2
! twin_c 25
! twin_d 6
! twin_e 9
! h0_slipslip 75e6
! h0_sliptwin 75e6
! h0_twinslip 75e6
! h0_twintwin 75e6
! interaction_slipslip 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
! interaction_sliptwin 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
! interaction_twinslip 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16
! interaction_twintwin 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
MODULE constitutive_phenopowerlaw
!*** Include other modules ***
use prec, only: pReal,pInt
implicit none
character (len=*), parameter :: constitutive_phenopowerlaw_label = 'phenopowerlaw'
integer(pInt), dimension(:), allocatable :: constitutive_phenopowerlaw_sizeDotState, &
constitutive_phenopowerlaw_sizeState, &
constitutive_phenopowerlaw_sizePostResults ! cumulative size of post results
integer(pInt), dimension(:,:), allocatable,target :: constitutive_phenopowerlaw_sizePostResult ! size of each post result output
character(len=64), dimension(:,:), allocatable,target :: constitutive_phenopowerlaw_output ! name of each post result output
character(len=32), dimension(:), allocatable :: constitutive_phenopowerlaw_structureName
integer(pInt), dimension(:), allocatable :: constitutive_phenopowerlaw_structure
integer(pInt), dimension(:,:), allocatable :: constitutive_phenopowerlaw_Nslip ! active number of slip systems per family
integer(pInt), dimension(:,:), allocatable :: constitutive_phenopowerlaw_Ntwin ! active number of twin systems per family
integer(pInt), dimension(:), allocatable :: constitutive_phenopowerlaw_totalNslip ! no. of slip system used in simulation
integer(pInt), dimension(:), allocatable :: constitutive_phenopowerlaw_totalNtwin ! no. of twin system used in simulation
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_CoverA
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_C11
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_C12
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_C13
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_C33
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_C44
real(pReal), dimension(:,:,:), allocatable :: constitutive_phenopowerlaw_Cslip_66
!* Visco-plastic constitutive_phenomenological parameters
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_gdot0_slip
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_n_slip
real(pReal), dimension(:,:), allocatable :: constitutive_phenopowerlaw_tau0_slip
real(pReal), dimension(:,:), allocatable :: constitutive_phenopowerlaw_tausat_slip
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_gdot0_twin
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_n_twin
real(pReal), dimension(:,:), allocatable :: constitutive_phenopowerlaw_tau0_twin
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_spr
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_twinB
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_twinC
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_twinD
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_twinE
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_h0_slipslip
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_h0_sliptwin
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_h0_twinslip
real(pReal), dimension(:), allocatable :: constitutive_phenopowerlaw_h0_twintwin
real(pReal), dimension(:,:), allocatable :: constitutive_phenopowerlaw_interaction_slipslip
real(pReal), dimension(:,:), allocatable :: constitutive_phenopowerlaw_interaction_sliptwin
real(pReal), dimension(:,:), allocatable :: constitutive_phenopowerlaw_interaction_twinslip
real(pReal), dimension(:,:), allocatable :: constitutive_phenopowerlaw_interaction_twintwin
real(pReal), dimension(:,:,:), allocatable :: constitutive_phenopowerlaw_hardeningMatrix_slipslip
real(pReal), dimension(:,:,:), allocatable :: constitutive_phenopowerlaw_hardeningMatrix_sliptwin
real(pReal), dimension(:,:,:), allocatable :: constitutive_phenopowerlaw_hardeningMatrix_twinslip
real(pReal), dimension(:,:,:), allocatable :: constitutive_phenopowerlaw_hardeningMatrix_twintwin
CONTAINS
!****************************************
!* - constitutive_init
!* - constitutive_stateInit
!* - constitutive_homogenizedC
!* - constitutive_microstructure
!* - constitutive_LpAndItsTangent
!* - consistutive_dotState
!* - consistutive_postResults
!****************************************
subroutine constitutive_phenopowerlaw_init(file)
!**************************************
!* Module initialization *
!**************************************
use prec, only: pInt, pReal
use math, only: math_Mandel3333to66, math_Voigt66to3333
use IO
use material
use lattice, only: lattice_initializeStructure, lattice_maxNslipFamily, lattice_maxNtwinFamily, &
lattice_maxNinteraction, lattice_NslipSystem, lattice_NtwinSystem, &
lattice_interactionSlipSlip,lattice_interactionSlipTwin,lattice_interactionTwinSlip,lattice_interactionTwinTwin
integer(pInt), intent(in) :: file
integer(pInt), parameter :: maxNchunks = 21
integer(pInt), dimension(1+2*maxNchunks) :: positions
integer(pInt) section, maxNinstance, i,j,k,l,m, output, mySize
character(len=64) tag
character(len=1024) line
write(6,*)
write(6,'(a20,a20,a12)') '<<<+- constitutive_',constitutive_phenopowerlaw_label,' init -+>>>'
write(6,*)
maxNinstance = count(phase_constitution == constitutive_phenopowerlaw_label)
if (maxNinstance == 0) return
allocate(constitutive_phenopowerlaw_sizeDotState(maxNinstance)) ; constitutive_phenopowerlaw_sizeDotState = 0_pInt
allocate(constitutive_phenopowerlaw_sizeState(maxNinstance)) ; constitutive_phenopowerlaw_sizeState = 0_pInt
allocate(constitutive_phenopowerlaw_sizePostResults(maxNinstance)); constitutive_phenopowerlaw_sizePostResults = 0_pInt
allocate(constitutive_phenopowerlaw_sizePostResult(maxval(phase_Noutput), &
maxNinstance)) ; constitutive_phenopowerlaw_sizePostResult = 0_pInt
allocate(constitutive_phenopowerlaw_output(maxval(phase_Noutput), &
maxNinstance)) ; constitutive_phenopowerlaw_output = ''
allocate(constitutive_phenopowerlaw_structureName(maxNinstance)) ; constitutive_phenopowerlaw_structureName = ''
allocate(constitutive_phenopowerlaw_structure(maxNinstance)) ; constitutive_phenopowerlaw_structure = 0_pInt
allocate(constitutive_phenopowerlaw_Nslip(lattice_maxNslipFamily,&
maxNinstance)) ; constitutive_phenopowerlaw_Nslip = 0_pInt
allocate(constitutive_phenopowerlaw_Ntwin(lattice_maxNtwinFamily,&
maxNinstance)) ; constitutive_phenopowerlaw_Ntwin = 0_pInt
allocate(constitutive_phenopowerlaw_totalNslip(maxNinstance)) ; constitutive_phenopowerlaw_totalNslip = 0_pInt !no. of slip system used in simulation (YJ.RO)
allocate(constitutive_phenopowerlaw_totalNtwin(maxNinstance)) ; constitutive_phenopowerlaw_totalNtwin = 0_pInt !no. of twin system used in simulation (YJ.RO)
allocate(constitutive_phenopowerlaw_CoverA(maxNinstance)) ; constitutive_phenopowerlaw_CoverA = 0.0_pReal
allocate(constitutive_phenopowerlaw_C11(maxNinstance)) ; constitutive_phenopowerlaw_C11 = 0.0_pReal
allocate(constitutive_phenopowerlaw_C12(maxNinstance)) ; constitutive_phenopowerlaw_C12 = 0.0_pReal
allocate(constitutive_phenopowerlaw_C13(maxNinstance)) ; constitutive_phenopowerlaw_C13 = 0.0_pReal
allocate(constitutive_phenopowerlaw_C33(maxNinstance)) ; constitutive_phenopowerlaw_C33 = 0.0_pReal
allocate(constitutive_phenopowerlaw_C44(maxNinstance)) ; constitutive_phenopowerlaw_C44 = 0.0_pReal
allocate(constitutive_phenopowerlaw_Cslip_66(6,6,maxNinstance)) ; constitutive_phenopowerlaw_Cslip_66 = 0.0_pReal
allocate(constitutive_phenopowerlaw_gdot0_slip(maxNinstance)) ; constitutive_phenopowerlaw_gdot0_slip = 0.0_pReal
allocate(constitutive_phenopowerlaw_n_slip(maxNinstance)) ; constitutive_phenopowerlaw_n_slip = 0.0_pReal
allocate(constitutive_phenopowerlaw_tau0_slip(lattice_maxNslipFamily,&
maxNinstance)) ; constitutive_phenopowerlaw_tau0_slip = 0.0_pReal
allocate(constitutive_phenopowerlaw_tausat_slip(lattice_maxNslipFamily,&
maxNinstance)) ; constitutive_phenopowerlaw_tausat_slip = 0.0_pReal
allocate(constitutive_phenopowerlaw_gdot0_twin(maxNinstance)) ; constitutive_phenopowerlaw_gdot0_twin = 0.0_pReal
allocate(constitutive_phenopowerlaw_n_twin(maxNinstance)) ; constitutive_phenopowerlaw_n_twin = 0.0_pReal
allocate(constitutive_phenopowerlaw_tau0_twin(lattice_maxNtwinFamily,&
maxNinstance)) ; constitutive_phenopowerlaw_tau0_twin = 0.0_pReal
allocate(constitutive_phenopowerlaw_spr(maxNinstance)) ; constitutive_phenopowerlaw_spr = 0.0_pReal
allocate(constitutive_phenopowerlaw_twinB(maxNinstance)) ; constitutive_phenopowerlaw_twinB = 0.0_pReal
allocate(constitutive_phenopowerlaw_twinC(maxNinstance)) ; constitutive_phenopowerlaw_twinC = 0.0_pReal
allocate(constitutive_phenopowerlaw_twinD(maxNinstance)) ; constitutive_phenopowerlaw_twinD = 0.0_pReal
allocate(constitutive_phenopowerlaw_twinE(maxNinstance)) ; constitutive_phenopowerlaw_twinE = 0.0_pReal
allocate(constitutive_phenopowerlaw_h0_slipslip(maxNinstance)) ; constitutive_phenopowerlaw_h0_slipslip = 0.0_pReal
allocate(constitutive_phenopowerlaw_h0_sliptwin(maxNinstance)) ; constitutive_phenopowerlaw_h0_sliptwin = 0.0_pReal
allocate(constitutive_phenopowerlaw_h0_twinslip(maxNinstance)) ; constitutive_phenopowerlaw_h0_twinslip = 0.0_pReal
allocate(constitutive_phenopowerlaw_h0_twintwin(maxNinstance)) ; constitutive_phenopowerlaw_h0_twintwin = 0.0_pReal
allocate(constitutive_phenopowerlaw_interaction_slipslip(lattice_maxNinteraction,&
maxNinstance)) ; constitutive_phenopowerlaw_interaction_slipslip = 0.0_pReal
allocate(constitutive_phenopowerlaw_interaction_sliptwin(lattice_maxNinteraction,&
maxNinstance)) ; constitutive_phenopowerlaw_interaction_sliptwin = 0.0_pReal
allocate(constitutive_phenopowerlaw_interaction_twinslip(lattice_maxNinteraction,&
maxNinstance)) ; constitutive_phenopowerlaw_interaction_twinslip = 0.0_pReal
allocate(constitutive_phenopowerlaw_interaction_twintwin(lattice_maxNinteraction,&
maxNinstance)) ; constitutive_phenopowerlaw_interaction_twintwin = 0.0_pReal
rewind(file)
line = ''
section = 0
do while (IO_lc(IO_getTag(line,'<','>')) /= 'phase') ! wind forward to <phase>
read(file,'(a1024)',END=100) line
enddo
do ! read thru sections of phase part
read(file,'(a1024)',END=100) line
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') exit ! stop at next part
if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1
output = 0 ! reset output counter
endif
if (section > 0 .and. phase_constitution(section) == constitutive_phenopowerlaw_label) then ! one of my sections
i = phase_constitutionInstance(section) ! which instance of my constitution is present phase
positions = IO_stringPos(line,maxNchunks)
tag = IO_lc(IO_stringValue(line,positions,1)) ! extract key
select case(tag)
case ('(output)')
output = output + 1
constitutive_phenopowerlaw_output(output,i) = IO_lc(IO_stringValue(line,positions,2))
case ('lattice_structure')
constitutive_phenopowerlaw_structureName(i) = IO_lc(IO_stringValue(line,positions,2))
case ('covera_ratio')
constitutive_phenopowerlaw_CoverA(i) = IO_floatValue(line,positions,2)
case ('c11')
constitutive_phenopowerlaw_C11(i) = IO_floatValue(line,positions,2)
case ('c12')
constitutive_phenopowerlaw_C12(i) = IO_floatValue(line,positions,2)
case ('c13')
constitutive_phenopowerlaw_C13(i) = IO_floatValue(line,positions,2)
case ('c33')
constitutive_phenopowerlaw_C33(i) = IO_floatValue(line,positions,2)
case ('c44')
constitutive_phenopowerlaw_C44(i) = IO_floatValue(line,positions,2)
case ('nslip')
forall (j = 1:lattice_maxNslipFamily) constitutive_phenopowerlaw_Nslip(j,i) = IO_intValue(line,positions,1+j)
case ('gdot0_slip')
constitutive_phenopowerlaw_gdot0_slip(i) = IO_floatValue(line,positions,2)
case ('n_slip')
constitutive_phenopowerlaw_n_slip(i) = IO_floatValue(line,positions,2)
case ('tau0_slip')
forall (j = 1:lattice_maxNslipFamily) constitutive_phenopowerlaw_tau0_slip(j,i) = IO_floatValue(line,positions,1+j)
case ('tausat_slip')
forall (j = 1:lattice_maxNslipFamily) constitutive_phenopowerlaw_tausat_slip(j,i) = IO_floatValue(line,positions,1+j)
case ('ntwin')
forall (j = 1:lattice_maxNtwinFamily) constitutive_phenopowerlaw_Ntwin(j,i) = IO_intValue(line,positions,1+j)
case ('gdot0_twin')
constitutive_phenopowerlaw_gdot0_twin(i) = IO_floatValue(line,positions,2)
case ('n_twin')
constitutive_phenopowerlaw_n_twin(i) = IO_floatValue(line,positions,2)
case ('tau0_twin')
forall (j = 1:lattice_maxNtwinFamily) constitutive_phenopowerlaw_tau0_twin(j,i) = IO_floatValue(line,positions,1+j)
case ('s_pr')
constitutive_phenopowerlaw_spr(i) = IO_floatValue(line,positions,2)
case ('twin_b')
constitutive_phenopowerlaw_twinB(i) = IO_floatValue(line,positions,2)
case ('twin_c')
constitutive_phenopowerlaw_twinC(i) = IO_floatValue(line,positions,2)
case ('twin_d')
constitutive_phenopowerlaw_twinD(i) = IO_floatValue(line,positions,2)
case ('twin_e')
constitutive_phenopowerlaw_twinE(i) = IO_floatValue(line,positions,2)
case ('h0_slipslip')
constitutive_phenopowerlaw_h0_slipslip(i) = IO_floatValue(line,positions,2)
case ('h0_sliptwin')
constitutive_phenopowerlaw_h0_sliptwin(i) = IO_floatValue(line,positions,2)
case ('h0_twinslip')
constitutive_phenopowerlaw_h0_twinslip(i) = IO_floatValue(line,positions,2)
case ('h0_twintwin')
constitutive_phenopowerlaw_h0_twintwin(i) = IO_floatValue(line,positions,2)
case ('interaction_slipslip')
forall (j = 1:lattice_maxNinteraction) &
constitutive_phenopowerlaw_interaction_slipslip(j,i) = IO_floatValue(line,positions,1+j)
case ('interaction_sliptwin')
forall (j = 1:lattice_maxNinteraction) &
constitutive_phenopowerlaw_interaction_sliptwin(j,i) = IO_floatValue(line,positions,1+j)
case ('interaction_twinslip')
forall (j = 1:lattice_maxNinteraction) &
constitutive_phenopowerlaw_interaction_twinslip(j,i) = IO_floatValue(line,positions,1+j)
case ('interaction_twintwin')
forall (j = 1:lattice_maxNinteraction) &
constitutive_phenopowerlaw_interaction_twintwin(j,i) = IO_floatValue(line,positions,1+j)
end select
endif
enddo
100 do i = 1,maxNinstance
constitutive_phenopowerlaw_structure(i) = lattice_initializeStructure(constitutive_phenopowerlaw_structureName(i), & ! get structure
constitutive_phenopowerlaw_CoverA(i))
constitutive_phenopowerlaw_Nslip(:,i) = min(lattice_NslipSystem(:,constitutive_phenopowerlaw_structure(i)),& ! limit active slip systems per family to max
constitutive_phenopowerlaw_Nslip(:,i))
constitutive_phenopowerlaw_Ntwin(:,i) = min(lattice_NtwinSystem(:,constitutive_phenopowerlaw_structure(i)),& ! limit active twin systems per family to max
constitutive_phenopowerlaw_Ntwin(:,i))
constitutive_phenopowerlaw_totalNslip(i) = sum(constitutive_phenopowerlaw_Nslip(:,i)) ! how many slip systems altogether
constitutive_phenopowerlaw_totalNtwin(i) = sum(constitutive_phenopowerlaw_Ntwin(:,i)) ! how many twin systems altogether
if (constitutive_phenopowerlaw_structure(i) < 1 .or. & ! sanity checks
constitutive_phenopowerlaw_structure(i) > 3) call IO_error(205)
if (any(constitutive_phenopowerlaw_tau0_slip(:,i) < 0.0_pReal .and. &
constitutive_phenopowerlaw_Nslip(:,i) > 0)) call IO_error(210)
if (constitutive_phenopowerlaw_gdot0_slip(i) <= 0.0_pReal) call IO_error(211)
if (constitutive_phenopowerlaw_n_slip(i) <= 0.0_pReal) call IO_error(212)
if (any(constitutive_phenopowerlaw_tausat_slip(:,i) <= 0.0_pReal .and. &
constitutive_phenopowerlaw_Nslip(:,i) > 0)) call IO_error(213)
if (any(constitutive_phenopowerlaw_tau0_twin(:,i) < 0.0_pReal .and. &
constitutive_phenopowerlaw_Ntwin(:,i) > 0)) call IO_error(210)
if (constitutive_phenopowerlaw_gdot0_twin(i) <= 0.0_pReal) call IO_error(211)
if (constitutive_phenopowerlaw_n_twin(i) <= 0.0_pReal) call IO_error(212)
enddo
allocate(constitutive_phenopowerlaw_hardeningMatrix_slipslip(maxval(constitutive_phenopowerlaw_totalNslip),&
maxval(constitutive_phenopowerlaw_totalNslip),&
maxNinstance))
allocate(constitutive_phenopowerlaw_hardeningMatrix_sliptwin(maxval(constitutive_phenopowerlaw_totalNslip),&
maxval(constitutive_phenopowerlaw_totalNtwin),&
maxNinstance))
allocate(constitutive_phenopowerlaw_hardeningMatrix_twinslip(maxval(constitutive_phenopowerlaw_totalNtwin),&
maxval(constitutive_phenopowerlaw_totalNslip),&
maxNinstance))
allocate(constitutive_phenopowerlaw_hardeningMatrix_twintwin(maxval(constitutive_phenopowerlaw_totalNtwin),&
maxval(constitutive_phenopowerlaw_totalNtwin),&
maxNinstance))
do i = 1,maxNinstance
do j = 1,maxval(phase_Noutput)
select case(constitutive_phenopowerlaw_output(j,i))
case('resistance_slip')
mySize = constitutive_phenopowerlaw_totalNslip(i)
case('shearrate_slip')
mySize = constitutive_phenopowerlaw_totalNslip(i)
case('resolvedstress_slip')
mySize = constitutive_phenopowerlaw_totalNslip(i)
case('totalshear')
mySize = 1_pInt
case('resistance_twin')
mySize = constitutive_phenopowerlaw_totalNtwin(i)
case('shearrate_twin')
mySize = constitutive_phenopowerlaw_totalNtwin(i)
case('resolvedstress_twin')
mySize = constitutive_phenopowerlaw_totalNtwin(i)
case('totalvolfrac')
mySize = 1_pInt
case default
mySize = 0_pInt
end select
if (mySize > 0_pInt) then ! any meaningful output found
constitutive_phenopowerlaw_sizePostResult(j,i) = mySize
constitutive_phenopowerlaw_sizePostResults(i) = &
constitutive_phenopowerlaw_sizePostResults(i) + mySize
endif
enddo
constitutive_phenopowerlaw_sizeDotState(i) = constitutive_phenopowerlaw_totalNslip(i)+ &
constitutive_phenopowerlaw_totalNtwin(i)+ 2 ! s_slip, s_twin, sum(gamma), sum(f)
constitutive_phenopowerlaw_sizeState(i) = constitutive_phenopowerlaw_totalNslip(i)+ &
constitutive_phenopowerlaw_totalNtwin(i)+ 2 ! s_slip, s_twin, sum(gamma), sum(f)
select case (constitutive_phenopowerlaw_structure(i)) ! assign elasticity tensor
case(1:2) ! cubic(s)
forall(k=1:3)
forall(j=1:3) &
constitutive_phenopowerlaw_Cslip_66(k,j,i) = constitutive_phenopowerlaw_C12(i)
constitutive_phenopowerlaw_Cslip_66(k,k,i) = constitutive_phenopowerlaw_C11(i)
constitutive_phenopowerlaw_Cslip_66(k+3,k+3,i) = constitutive_phenopowerlaw_C44(i)
end forall
case(3) ! hex
constitutive_phenopowerlaw_Cslip_66(1,1,i) = constitutive_phenopowerlaw_C11(i)
constitutive_phenopowerlaw_Cslip_66(2,2,i) = constitutive_phenopowerlaw_C11(i)
constitutive_phenopowerlaw_Cslip_66(3,3,i) = constitutive_phenopowerlaw_C33(i)
constitutive_phenopowerlaw_Cslip_66(1,2,i) = constitutive_phenopowerlaw_C12(i)
constitutive_phenopowerlaw_Cslip_66(2,1,i) = constitutive_phenopowerlaw_C12(i)
constitutive_phenopowerlaw_Cslip_66(1,3,i) = constitutive_phenopowerlaw_C13(i)
constitutive_phenopowerlaw_Cslip_66(3,1,i) = constitutive_phenopowerlaw_C13(i)
constitutive_phenopowerlaw_Cslip_66(2,3,i) = constitutive_phenopowerlaw_C13(i)
constitutive_phenopowerlaw_Cslip_66(3,2,i) = constitutive_phenopowerlaw_C13(i)
constitutive_phenopowerlaw_Cslip_66(4,4,i) = constitutive_phenopowerlaw_C44(i)
constitutive_phenopowerlaw_Cslip_66(5,5,i) = constitutive_phenopowerlaw_C44(i)
constitutive_phenopowerlaw_Cslip_66(6,6,i) = 0.5_pReal*(constitutive_phenopowerlaw_C11(i)- &
constitutive_phenopowerlaw_C12(i))
end select
constitutive_phenopowerlaw_Cslip_66(:,:,i) = &
math_Mandel3333to66(math_Voigt66to3333(constitutive_phenopowerlaw_Cslip_66(:,:,i)))
do j = 1,lattice_maxNslipFamily
do k = 1,constitutive_phenopowerlaw_Nslip(j,i)
do l = 1,lattice_maxNslipFamily
do m = 1,constitutive_phenopowerlaw_Nslip(l,i)
constitutive_phenopowerlaw_hardeningMatrix_slipslip(sum(constitutive_phenopowerlaw_Nslip(1:j-1,i))+k,&
sum(constitutive_phenopowerlaw_Nslip(1:l-1,i))+m, i) = &
constitutive_phenopowerlaw_interaction_slipslip(lattice_interactionSlipSlip( &
sum(lattice_NslipSystem(1:j-1,constitutive_phenopowerlaw_structure(i)))+k, &
sum(lattice_NslipSystem(1:l-1,constitutive_phenopowerlaw_structure(i)))+m, &
constitutive_phenopowerlaw_structure(i)), i )
enddo; enddo; enddo; enddo
do j = 1,lattice_maxNslipFamily
do k = 1,constitutive_phenopowerlaw_Nslip(j,i)
do l = 1,lattice_maxNtwinFamily
do m = 1,constitutive_phenopowerlaw_Ntwin(l,i)
constitutive_phenopowerlaw_hardeningMatrix_sliptwin(sum(constitutive_phenopowerlaw_Nslip(1:j-1,i))+k,&
sum(constitutive_phenopowerlaw_Ntwin(1:l-1,i))+m, i) = &
constitutive_phenopowerlaw_interaction_sliptwin(lattice_interactionSlipTwin( &
sum(lattice_NslipSystem(1:j-1,constitutive_phenopowerlaw_structure(i)))+k, &
sum(lattice_NtwinSystem(1:l-1,constitutive_phenopowerlaw_structure(i)))+m, &
constitutive_phenopowerlaw_structure(i)) ,i )
enddo; enddo; enddo; enddo
do j = 1,lattice_maxNtwinFamily
do k = 1,constitutive_phenopowerlaw_Ntwin(j,i)
do l = 1,lattice_maxNslipFamily
do m = 1,constitutive_phenopowerlaw_Nslip(l,i)
constitutive_phenopowerlaw_hardeningMatrix_twinslip(sum(constitutive_phenopowerlaw_Ntwin(1:j-1,i))+k,&
sum(constitutive_phenopowerlaw_Nslip(1:l-1,i))+m, i) = &
constitutive_phenopowerlaw_interaction_twinslip(lattice_interactionTwinSlip( &
sum(lattice_NtwinSystem(1:j-1,constitutive_phenopowerlaw_structure(i)))+k, &
sum(lattice_NslipSystem(1:l-1,constitutive_phenopowerlaw_structure(i)))+m, &
constitutive_phenopowerlaw_structure(i)), i )
enddo; enddo; enddo; enddo
do j = 1,lattice_maxNtwinFamily
do k = 1,constitutive_phenopowerlaw_Ntwin(j,i)
do l = 1,lattice_maxNtwinFamily
do m = 1,constitutive_phenopowerlaw_Ntwin(l,i)
constitutive_phenopowerlaw_hardeningMatrix_twintwin(sum(constitutive_phenopowerlaw_Ntwin(1:j-1,i))+k,&
sum(constitutive_phenopowerlaw_Ntwin(1:l-1,i))+m, i) = &
constitutive_phenopowerlaw_interaction_twintwin(lattice_interactionTwinTwin( &
sum(lattice_NtwinSystem(1:j-1,constitutive_phenopowerlaw_structure(i)))+k, &
sum(lattice_NtwinSystem(1:l-1,constitutive_phenopowerlaw_structure(i)))+m, &
constitutive_phenopowerlaw_structure(i)), i )
enddo; enddo; enddo; enddo
enddo
return
end subroutine
function constitutive_phenopowerlaw_stateInit(myInstance)
!*********************************************************************
!* initial microstructural state *
!*********************************************************************
use prec, only: pReal,pInt
use debug, only: debugger
use lattice, only: lattice_maxNslipFamily, lattice_maxNtwinFamily
implicit none
!* Definition of variables
integer(pInt), intent(in) :: myInstance
integer(pInt) i
real(pReal), dimension(constitutive_phenopowerlaw_totalNslip(myInstance) + &
constitutive_phenopowerlaw_totalNtwin(myInstance) + 2) :: constitutive_phenopowerlaw_stateInit
constitutive_phenopowerlaw_stateInit = 0.0_pReal
do i = 1,lattice_maxNslipFamily
constitutive_phenopowerlaw_stateInit(1+&
sum(constitutive_phenopowerlaw_Nslip(1:i-1,myInstance)) : &
sum(constitutive_phenopowerlaw_Nslip(1:i ,myInstance))) = &
constitutive_phenopowerlaw_tau0_slip(i,myInstance)
enddo
do i = 1,lattice_maxNtwinFamily
constitutive_phenopowerlaw_stateInit(1+sum(constitutive_phenopowerlaw_Nslip(:,myInstance))+&
sum(constitutive_phenopowerlaw_Ntwin(1:i-1,myInstance)) : &
sum(constitutive_phenopowerlaw_Nslip(:,myInstance))+&
sum(constitutive_phenopowerlaw_Ntwin(1:i ,myInstance))) = &
constitutive_phenopowerlaw_tau0_twin(i,myInstance)
enddo
return
end function
function constitutive_phenopowerlaw_homogenizedC(state,ipc,ip,el)
!*********************************************************************
!* homogenized elacticity matrix *
!* INPUT: *
!* - state : state variables *
!* - ipc : component-ID of current integration point *
!* - ip : current integration point *
!* - el : current element *
!*********************************************************************
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance
implicit none
!* Definition of variables
integer(pInt), intent(in) :: ipc,ip,el
integer(pInt) matID
real(pReal), dimension(6,6) :: constitutive_phenopowerlaw_homogenizedC
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state
matID = phase_constitutionInstance(material_phase(ipc,ip,el))
constitutive_phenopowerlaw_homogenizedC = constitutive_phenopowerlaw_Cslip_66(:,:,matID)
return
end function
subroutine constitutive_phenopowerlaw_microstructure(Temperature,state,ipc,ip,el)
!*********************************************************************
!* calculate derived quantities from state (not used here) *
!* INPUT: *
!* - Tp : temperature *
!* - ipc : component-ID of current integration point *
!* - ip : current integration point *
!* - el : current element *
!*********************************************************************
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance
implicit none
!* Definition of variables
integer(pInt) ipc,ip,el, matID,nSlip,nTwin
real(pReal) Temperature
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state
matID = phase_constitutionInstance(material_phase(ipc,ip,el))
end subroutine
subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperature,state,ipc,ip,el)
!*********************************************************************
!* plastic velocity gradient and its tangent *
!* INPUT: *
!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) *
!* - ipc : component-ID at current integration point *
!* - ip : current integration point *
!* - el : current element *
!* OUTPUT: *
!* - Lp : plastic velocity gradient *
!* - dLp_dTstar : derivative of Lp (4th-rank tensor) *
!*********************************************************************
use prec, only: pReal,pInt,p_vec
use debug, only: debugger
use math, only: math_Plain3333to99
use lattice, only: lattice_Sslip,lattice_Sslip_v,lattice_Stwin,lattice_Stwin_v, lattice_maxNslipFamily, lattice_maxNtwinFamily, &
lattice_NslipSystem,lattice_NtwinSystem
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance
implicit none
!* Definition of variables
integer(pInt) ipc,ip,el
integer(pInt) matID,nSlip,nTwin,f,i,j,k,l,m,n, structID,index_Gamma,index_F,index_myFamily
real(pReal) Temperature
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state
real(pReal), dimension(6) :: Tstar_v
real(pReal), dimension(3,3) :: Lp
real(pReal), dimension(3,3,3,3) :: dLp_dTstar3333
real(pReal), dimension(9,9) :: dLp_dTstar
real(pReal), dimension(constitutive_phenopowerlaw_totalNslip(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: &
gdot_slip,dgdot_dtauslip,tau_slip
real(pReal), dimension(constitutive_phenopowerlaw_totalNtwin(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: &
gdot_twin,dgdot_dtautwin,tau_twin
matID = phase_constitutionInstance(material_phase(ipc,ip,el))
structID = constitutive_phenopowerlaw_structure(matID)
nSlip = constitutive_phenopowerlaw_totalNslip(matID)
nTwin = constitutive_phenopowerlaw_totalNtwin(matID)
index_Gamma = nSlip + nTwin + 1
index_F = nSlip + nTwin + 2
Lp = 0.0_pReal
dLp_dTstar3333 = 0.0_pReal
dLp_dTstar = 0.0_pReal
j = 0_pInt
do f = 1,lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum(lattice_NslipSystem(1:f-1,structID)) ! at which index starts my family
do i = 1,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family
j = j+1_pInt
!* Calculation of Lp
tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,structID))
gdot_slip(j) = constitutive_phenopowerlaw_gdot0_slip(matID)*(abs(tau_slip(j))/state(ipc,ip,el)%p(j))**&
constitutive_phenopowerlaw_n_slip(matID)*sign(1.0_pReal,tau_slip(j))
Lp = Lp + (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F
gdot_slip(j)*lattice_Sslip(:,:,index_myFamily+i,structID)
!* Calculation of the tangent of Lp
dgdot_dtauslip(j) = gdot_slip(j)*constitutive_phenopowerlaw_n_slip(matID)/tau_slip(j)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
dgdot_dtauslip(j)*lattice_Sslip(k,l,index_myFamily+i,structID)* &
lattice_Sslip(m,n,index_myFamily+i,structID)
enddo
enddo
j = 0_pInt
do f = 1,lattice_maxNtwinFamily ! loop over all twin families
index_myFamily = sum(lattice_NtwinSystem(1:f-1,structID)) ! at which index starts my family
do i = 1,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family
j = j+1_pInt
!* Calculation of Lp
tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,structID))
gdot_twin(j) = constitutive_phenopowerlaw_gdot0_twin(matID)*&
(abs(tau_twin(j))/state(ipc,ip,el)%p(nSlip+j))**&
constitutive_phenopowerlaw_n_twin(matID)*max(0.0_pReal,sign(1.0_pReal,tau_twin(j)))
Lp = Lp + gdot_twin(j)*lattice_Stwin(:,:,index_myFamily+i,structID)
!* Calculation of the tangent of Lp
dgdot_dtautwin(j) = gdot_twin(j)*constitutive_phenopowerlaw_n_twin(matID)/tau_twin(j)
forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
dgdot_dtautwin(j)*lattice_Stwin(k,l,index_myFamily+i,structID)* &
lattice_Stwin(m,n,index_myFamily+i,structID)
enddo
enddo
dLp_dTstar = math_Plain3333to99(dLp_dTstar3333)
return
end subroutine
function constitutive_phenopowerlaw_dotState(Tstar_v,Temperature,state,ipc,ip,el)
!*********************************************************************
!* rate of change of microstructure *
!* INPUT: *
!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) *
!* - ipc : component-ID at current integration point *
!* - ip : current integration point *
!* - el : current element *
!* OUTPUT: *
!* - constitutive_dotState : evolution of state variable *
!*********************************************************************
use prec, only: pReal,pInt,p_vec
use debug, only: debugger
use lattice, only: lattice_Sslip,lattice_Sslip_v,lattice_Stwin,lattice_Stwin_v, lattice_maxNslipFamily, lattice_maxNtwinFamily, &
lattice_NslipSystem,lattice_NtwinSystem,lattice_shearTwin
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance
implicit none
!* Definition of variables
integer(pInt) ipc,ip,el
integer(pInt) matID,nSlip,nTwin,f,i,j,k, structID,index_Gamma,index_F,index_myFamily
real(pReal) Temperature,c_slipslip,c_sliptwin,c_twinslip,c_twintwin, ssat
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: state
real(pReal), dimension(6) :: Tstar_v
real(pReal), dimension(constitutive_phenopowerlaw_totalNslip(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: &
gdot_slip,tau_slip,h_slipslip,h_sliptwin,N_slipslip,N_twinslip
real(pReal), dimension(constitutive_phenopowerlaw_totalNtwin(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: &
gdot_twin,tau_twin,h_twinslip,h_twintwin,N_sliptwin,N_twintwin
real(pReal), dimension(constitutive_phenopowerlaw_sizeDotState(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: &
constitutive_phenopowerlaw_dotState
matID = phase_constitutionInstance(material_phase(ipc,ip,el))
structID = constitutive_phenopowerlaw_structure(matID)
nSlip = constitutive_phenopowerlaw_totalNslip(matID)
nTwin = constitutive_phenopowerlaw_totalNtwin(matID)
index_Gamma = nSlip + nTwin + 1
index_F = nSlip + nTwin + 2
constitutive_phenopowerlaw_dotState = 0.0_pReal
!-- system-independent (nonlinear) prefactors to M_xx matrices
c_slipslip = constitutive_phenopowerlaw_h0_slipslip(matID)*&
(1.0_pReal + &
constitutive_phenopowerlaw_twinC(matID)*state(ipc,ip,el)%p(index_F)**constitutive_phenopowerlaw_twinB(matID))
c_sliptwin = 0.0_pReal
c_twinslip = constitutive_phenopowerlaw_h0_twinslip(matID)*&
state(ipc,ip,el)%p(index_Gamma)**constitutive_phenopowerlaw_twinE(matID)
c_twintwin = constitutive_phenopowerlaw_h0_twintwin(matID)*&
state(ipc,ip,el)%p(index_F)**constitutive_phenopowerlaw_twinD(matID)
!-- add system-dependent part and calculate dot gammas
j = 0_pInt
do f = 1,lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum(lattice_NslipSystem(1:f-1,structID)) ! at which index starts my family
ssat = constitutive_phenopowerlaw_tausat_slip(f,matID) + &
constitutive_phenopowerlaw_spr(matID)*dsqrt(state(ipc,ip,el)%p(index_F))
do i = 1,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family
j = j+1_pInt
h_slipslip(j) = c_slipslip*(1.0_pReal-state(ipc,ip,el)%p(j)/ssat) ! system-dependent prefactor for slip--slip interaction
h_sliptwin(j) = c_sliptwin ! no system-dependent part
!* Calculation of dot gamma
tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,structID))
gdot_slip(j) = constitutive_phenopowerlaw_gdot0_slip(matID)*(abs(tau_slip(j))/state(ipc,ip,el)%p(j))**&
constitutive_phenopowerlaw_n_slip(matID)*sign(1.0_pReal,tau_slip(j))
enddo
enddo
j = 0_pInt
do f = 1,lattice_maxNtwinFamily ! loop over all twin families
index_myFamily = sum(lattice_NtwinSystem(1:f-1,structID)) ! at which index starts my family
do i = 1,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family
j = j+1_pInt
h_twinslip(j) = c_twinslip ! no system-dependent parts
h_twintwin(j) = c_twintwin
!* Calculation of dot vol frac
tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,structID))
gdot_twin(j) = constitutive_phenopowerlaw_gdot0_twin(matID)*&
(abs(tau_twin(j))/state(ipc,ip,el)%p(nSlip+j))**&
constitutive_phenopowerlaw_n_twin(matID)*max(0.0_pReal,sign(1.0_pReal,tau_twin(j)))
enddo
enddo
!-- calculate the overall hardening based on above
j = 0_pInt
do f = 1,lattice_maxNslipFamily ! loop over all slip families
do i = 1,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family
j = j+1_pInt
forall (k=1:nSlip) N_slipslip(k) = constitutive_phenopowerlaw_hardeningMatrix_slipslip(j,k,matID) * &
abs(gdot_slip(k)) ! dot gamma_slip
forall (k=1:nTwin) N_sliptwin(k) = constitutive_phenopowerlaw_hardeningMatrix_sliptwin(j,k,matID) * &
gdot_twin(k) ! dot gamma_twin
constitutive_phenopowerlaw_dotState(j) = h_slipslip(j)*sum(N_slipslip) + & ! evolution of slip resistance j
h_sliptwin(j)*sum(N_sliptwin)
constitutive_phenopowerlaw_dotState(index_Gamma) = constitutive_phenopowerlaw_dotState(index_Gamma) + &
abs(gdot_slip(j))
enddo
enddo
j = 0_pInt
do f = 1,lattice_maxNtwinFamily ! loop over all twin families
index_myFamily = sum(lattice_NtwinSystem(1:f-1,structID)) ! at which index starts my family
do i = 1,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family
j = j+1_pInt
forall (k=1:nSlip) N_twinslip(k) = constitutive_phenopowerlaw_hardeningMatrix_twinslip(j,k,matID) * &
gdot_slip(k) ! dot gamma_slip
forall (k=1:nTwin) N_twintwin(k) = constitutive_phenopowerlaw_hardeningMatrix_twintwin(j,k,matID) * &
gdot_twin(k) ! dot gamma_twin
constitutive_phenopowerlaw_dotState(j+nSlip) = h_twinslip(j)*sum(N_twinslip) + & ! evolution of twin resistance j
h_twintwin(j)*sum(N_twintwin)
constitutive_phenopowerlaw_dotState(index_F) = constitutive_phenopowerlaw_dotState(index_F) + &
gdot_twin(j)/lattice_shearTwin(index_myFamily+i,structID)
enddo
enddo
return
end function
!****************************************************************
!* calculates the rate of change of temperature *
!****************************************************************
pure function constitutive_phenopowerlaw_dotTemperature(Tstar_v,Temperature,state,ipc,ip,el)
!*** variables and functions from other modules ***!
use prec, only: pReal,pInt,p_vec
use lattice, only: lattice_Sslip_v
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance
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):: ipc, & ! 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_phenopowerlaw_dotTemperature ! rate of change of temparature
! calculate dotTemperature
constitutive_phenopowerlaw_dotTemperature = 0.0_pReal
return
endfunction
pure function constitutive_phenopowerlaw_postResults(Tstar_v,Temperature,dt,state,ipc,ip,el)
!*********************************************************************
!* return array of constitutive results *
!* INPUT: *
!* - Tstar_v : 2nd Piola Kirchhoff stress tensor (Mandel) *
!* - dt : current time increment *
!* - ipc : component-ID at current integration point *
!* - ip : current integration point *
!* - el : current element *
!*********************************************************************
use prec, only: pReal,pInt,p_vec
use lattice, only: lattice_Sslip_v,lattice_Stwin_v, lattice_maxNslipFamily, lattice_maxNtwinFamily, &
lattice_NslipSystem,lattice_NtwinSystem
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance,phase_Noutput
implicit none
!* Definition of variables
integer(pInt), intent(in) :: ipc,ip,el
real(pReal), intent(in) :: dt,Temperature
real(pReal), dimension(6), intent(in) :: Tstar_v
type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state
integer(pInt) matID,o,f,i,c,nSlip,nTwin,j, structID,index_Gamma,index_F,index_myFamily
real(pReal) tau
real(pReal), dimension(constitutive_phenopowerlaw_sizePostResults(phase_constitutionInstance(material_phase(ipc,ip,el)))) :: &
constitutive_phenopowerlaw_postResults
matID = phase_constitutionInstance(material_phase(ipc,ip,el))
structID = constitutive_phenopowerlaw_structure(matID)
nSlip = constitutive_phenopowerlaw_totalNslip(matID)
nTwin = constitutive_phenopowerlaw_totalNtwin(matID)
index_Gamma = nSlip + nTwin + 1
index_F = nSlip + nTwin + 2
constitutive_phenopowerlaw_postResults = 0.0_pReal
c = 0_pInt
do o = 1,phase_Noutput(material_phase(ipc,ip,el))
select case(constitutive_phenopowerlaw_output(o,matID))
case ('resistance_slip')
constitutive_phenopowerlaw_postResults(c+1:c+nSlip) = state(ipc,ip,el)%p(1:nSlip)
c = c + nSlip
case ('shearrate_slip')
j = 0_pInt
do f = 1,lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum(lattice_NslipSystem(1:f-1,structID)) ! at which index starts my family
do i = 1,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family
j = j + 1_pInt
tau = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,structID))
constitutive_phenopowerlaw_postResults(c+j) = constitutive_phenopowerlaw_gdot0_slip(matID)*&
(abs(tau)/state(ipc,ip,el)%p(j))**&
constitutive_phenopowerlaw_n_slip(matID)*sign(1.0_pReal,tau)
enddo; enddo
c = c + nSlip
case ('resolvedstress_slip')
j = 0_pInt
do f = 1,lattice_maxNslipFamily ! loop over all slip families
index_myFamily = sum(lattice_NslipSystem(1:f-1,structID)) ! at which index starts my family
do i = 1,constitutive_phenopowerlaw_Nslip(f,matID) ! process each (active) slip system in family
j = j + 1_pInt
constitutive_phenopowerlaw_postResults(c+j) = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,structID))
enddo; enddo
c = c + nSlip
case ('totalshear')
constitutive_phenopowerlaw_postResults(c+1) = state(ipc,ip,el)%p(index_Gamma)
c = c + 1
case ('resistance_twin')
constitutive_phenopowerlaw_postResults(c+1:c+nTwin) = state(ipc,ip,el)%p(1+nSlip:nTwin+nSlip)
c = c + nTwin
case ('shearrate_twin')
j = 0_pInt
do f = 1,lattice_maxNtwinFamily ! loop over all twin families
index_myFamily = sum(lattice_NtwinSystem(1:f-1,structID)) ! at which index starts my family
do i = 1,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family
j = j + 1_pInt
tau = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,structID))
constitutive_phenopowerlaw_postResults(c+j) = constitutive_phenopowerlaw_gdot0_twin(matID)*&
(abs(tau)/state(ipc,ip,el)%p(j+nSlip))**&
constitutive_phenopowerlaw_n_twin(matID)*max(0.0_pReal,sign(1.0_pReal,tau))
enddo; enddo
c = c + nTwin
case ('resolvedstress_twin')
j = 0_pInt
do f = 1,lattice_maxNtwinFamily ! loop over all twin families
index_myFamily = sum(lattice_NtwinSystem(1:f-1,structID)) ! at which index starts my family
do i = 1,constitutive_phenopowerlaw_Ntwin(f,matID) ! process each (active) twin system in family
j = j + 1_pInt
constitutive_phenopowerlaw_postResults(c+j) = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,structID))
enddo; enddo
c = c + nTwin
case ('totalvolfrac')
constitutive_phenopowerlaw_postResults(c+1) = state(ipc,ip,el)%p(index_F)
c = c + 1
end select
enddo
return
end function
END MODULE

149
code/crystallite.f90
View File

@ -55,7 +55,6 @@ logical, dimension (:,:,:), allocatable :: crystallite_localConstit
crystallite_onTrack, & ! flag to indicate ongoing calculation crystallite_onTrack, & ! flag to indicate ongoing calculation
crystallite_converged ! convergence flag crystallite_converged ! convergence flag
CONTAINS CONTAINS
!******************************************************************** !********************************************************************
@ -128,7 +127,7 @@ subroutine crystallite_init(Temperature)
allocate(crystallite_subdt(gMax,iMax,eMax)); crystallite_subdt = 0.0_pReal allocate(crystallite_subdt(gMax,iMax,eMax)); crystallite_subdt = 0.0_pReal
allocate(crystallite_subFrac(gMax,iMax,eMax)); crystallite_subFrac = 0.0_pReal allocate(crystallite_subFrac(gMax,iMax,eMax)); crystallite_subFrac = 0.0_pReal
allocate(crystallite_subStep(gMax,iMax,eMax)); crystallite_subStep = 0.0_pReal allocate(crystallite_subStep(gMax,iMax,eMax)); crystallite_subStep = 0.0_pReal
allocate(crystallite_localConstitution(gMax,iMax,eMax)); allocate(crystallite_localConstitution(gMax,iMax,eMax)); crystallite_localConstitution = .true.
allocate(crystallite_requested(gMax,iMax,eMax)); crystallite_requested = .false. allocate(crystallite_requested(gMax,iMax,eMax)); crystallite_requested = .false.
allocate(crystallite_onTrack(gMax,iMax,eMax)); crystallite_onTrack = .false. allocate(crystallite_onTrack(gMax,iMax,eMax)); crystallite_onTrack = .false.
allocate(crystallite_converged(gMax,iMax,eMax)); crystallite_converged = .true. allocate(crystallite_converged(gMax,iMax,eMax)); crystallite_converged = .true.
@ -138,7 +137,7 @@ subroutine crystallite_init(Temperature)
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element
do g = 1,myNgrains do g = 1,myNgrains
crystallite_partionedTemperature0(g,i,e) = Temperature ! isothermal assumption crystallite_partionedTemperature0(g,i,e) = Temperature ! isothermal assumption
crystallite_Fp0(:,:,g,i,e) = math_EulerToR(material_EulerAngles(:,g,i,e)) ! plastic def gradient reflects init orientation crystallite_Fp0(:,:,g,i,e) = math_EulerToR(material_EulerAngles(:,g,i,e)) ! plastic def gradient reflects init orientation
crystallite_F0(:,:,g,i,e) = math_I3 crystallite_F0(:,:,g,i,e) = math_I3
crystallite_partionedFp0(:,:,g,i,e) = crystallite_Fp0(:,:,g,i,e) crystallite_partionedFp0(:,:,g,i,e) = crystallite_Fp0(:,:,g,i,e)
@ -146,8 +145,7 @@ subroutine crystallite_init(Temperature)
crystallite_partionedF(:,:,g,i,e) = crystallite_F0(:,:,g,i,e) crystallite_partionedF(:,:,g,i,e) = crystallite_F0(:,:,g,i,e)
crystallite_requested(g,i,e) = .true. crystallite_requested(g,i,e) = .true.
crystallite_localConstitution(g,i,e) = phase_localConstitution(material_phase(g,i,e)) crystallite_localConstitution(g,i,e) = phase_localConstitution(material_phase(g,i,e))
enddo
enddo
enddo enddo
enddo enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
@ -274,16 +272,17 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
converged ! flag indicating if iteration converged converged ! flag indicating if iteration converged
! ------ initialize to starting condition ------ ! ------ initialize to starting condition ------
write (6,*) !$OMP CRITICAL (write2out)
write (6,*) 'Crystallite request from Materialpoint' ! write (6,*)
write (6,'(a,/,(f12.7,x))') 'crystallite_partionedTemp0 of 1 1 1' ,crystallite_partionedTemperature0(1,1,1) ! write (6,*) 'Crystallite request from Materialpoint'
write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedF0 of 1 1 1' ,crystallite_partionedF0(1:3,:,1,1,1) ! write (6,'(a,/,(f12.7,x))') 'crystallite_partionedTemp0 of 1 1 1' ,crystallite_partionedTemperature0(1,1,1)
write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedFp0 of 1 1 1' ,crystallite_partionedFp0(1:3,:,1,1,1) ! write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedF0 of 1 1 1' ,crystallite_partionedF0(1:3,:,1,1,1)
write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedF of 1 1 1' ,crystallite_partionedF(1:3,:,1,1,1) ! write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedFp0 of 1 1 1' ,crystallite_partionedFp0(1:3,:,1,1,1)
write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedLp0 of 1 1 1' ,crystallite_partionedLp0(1:3,:,1,1,1) ! write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedF of 1 1 1' ,crystallite_partionedF(1:3,:,1,1,1)
! write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedLp0 of 1 1 1' ,crystallite_partionedLp0(1:3,:,1,1,1)
!$OMPEND CRITICAL (write2out)
!$OMP PARALLEL DO !$OMP PARALLEL DO
@ -329,7 +328,16 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains do g = 1,myNgrains
debugger = (e == 1 .and. i == 1 .and. g == 1)
if (crystallite_converged(g,i,e)) then if (crystallite_converged(g,i,e)) then
if (debugger) then
!$OMP CRITICAL (write2out)
write(6,'(a21,f10.8,a33,f10.8,a35)') 'winding forward from ', &
crystallite_subFrac(g,i,e),' to current crystallite_subfrac ', &
crystallite_subFrac(g,i,e)+crystallite_subStep(g,i,e),' in crystallite_stressAndItsTangent'
write(6,*)
!$OMPEND CRITICAL (write2out)
endif
crystallite_subFrac(g,i,e) = crystallite_subFrac(g,i,e) + crystallite_subStep(g,i,e) crystallite_subFrac(g,i,e) = crystallite_subFrac(g,i,e) + crystallite_subStep(g,i,e)
crystallite_subStep(g,i,e) = min(1.0_pReal-crystallite_subFrac(g,i,e), 2.0_pReal * crystallite_subStep(g,i,e)) crystallite_subStep(g,i,e) = min(1.0_pReal-crystallite_subFrac(g,i,e), 2.0_pReal * crystallite_subStep(g,i,e))
if (crystallite_subStep(g,i,e) > subStepMin) then if (crystallite_subStep(g,i,e) > subStepMin) then
@ -340,14 +348,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
constitutive_subState0(g,i,e)%p = constitutive_state(g,i,e)%p ! ...microstructure constitutive_subState0(g,i,e)%p = constitutive_state(g,i,e)%p ! ...microstructure
crystallite_subTstar0_v(:,g,i,e) = crystallite_Tstar_v(:,g,i,e) ! ...2nd PK stress crystallite_subTstar0_v(:,g,i,e) = crystallite_Tstar_v(:,g,i,e) ! ...2nd PK stress
endif endif
if (debugger) then
!$OMP CRITICAL (write2out)
write(6,'(a21,f6.4,a28,f6.4,a35)') 'winding forward from ', &
crystallite_subFrac(g,i,e)-crystallite_subStep(g,i,e),' to new crystallite_subfrac ', &
crystallite_subFrac(g,i,e),' in crystallite_stressAndItsTangent'
write(6,*)
!$OMPEND CRITICAL (write2out)
endif
else else
crystallite_subStep(g,i,e) = 0.5_pReal * crystallite_subStep(g,i,e) ! cut step in half and restore... crystallite_subStep(g,i,e) = 0.5_pReal * crystallite_subStep(g,i,e) ! cut step in half and restore...
crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) ! ...temperature crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) ! ...temperature
@ -357,8 +357,8 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
crystallite_Tstar_v(:,g,i,e) = crystallite_subTstar0_v(:,g,i,e) ! ...2nd PK stress crystallite_Tstar_v(:,g,i,e) = crystallite_subTstar0_v(:,g,i,e) ! ...2nd PK stress
if (debugger) then if (debugger) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write(6,'(a78,f6.4)') 'cutback step in crystallite_stressAndItsTangent with new crystallite_subStep: ',& write(6,'(a78,f10.8)') 'cutback step in crystallite_stressAndItsTangent with new crystallite_subStep: ',&
crystallite_subStep(g,i,e) crystallite_subStep(g,i,e)
write(6,*) write(6,*)
!$OMPEND CRITICAL (write2out) !$OMPEND CRITICAL (write2out)
endif endif
@ -389,8 +389,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
! based on constitutive_subState0 ! based on constitutive_subState0
! results in constitutive_state ! results in constitutive_state
if (debugger) write(6,*) 'state integration started'
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
@ -431,6 +429,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains do g = 1,myNgrains
debugger = (e == 1 .and. i == 1 .and. g == 1)
if ( crystallite_requested(g,i,e) & if ( crystallite_requested(g,i,e) &
.and. crystallite_onTrack(g,i,e) & .and. crystallite_onTrack(g,i,e) &
.and. .not. crystallite_converged(g,i,e) ) & ! all undone crystallites .and. .not. crystallite_converged(g,i,e) ) & ! all undone crystallites
@ -451,11 +450,13 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains do g = 1,myNgrains
debugger = (e == 1 .and. i == 1 .and. g == 1)
if ( crystallite_requested(g,i,e) & if ( crystallite_requested(g,i,e) &
.and. crystallite_onTrack(g,i,e) & .and. crystallite_onTrack(g,i,e) &
.and. .not. crystallite_converged(g,i,e)) then ! all undone crystallites .and. .not. crystallite_converged(g,i,e)) then ! all undone crystallites
crystallite_converged(g,i,e) = (crystallite_updateState(g,i,e).AND.& crystallite_converged(g,i,e) = (crystallite_updateState(g,i,e) .and. &
crystallite_updateTemperature(g,i,e)) crystallite_updateTemperature(g,i,e))
if (debugger) write (6,*) g,i,e,'converged after updState',crystallite_converged(g,i,e)
if (crystallite_converged(g,i,e)) then if (crystallite_converged(g,i,e)) then
!$OMP CRITICAL (distributionState) !$OMP CRITICAL (distributionState)
debug_StateLoopDistribution(NiterationState) = debug_StateLoopDistribution(NiterationState) + 1 debug_StateLoopDistribution(NiterationState) = debug_StateLoopDistribution(NiterationState) + 1
@ -473,8 +474,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
enddo ! crystallite convergence loop enddo ! crystallite convergence loop
if (debugger) write(6,*) 'state integration converged'
enddo ! cutback loop enddo ! cutback loop
@ -504,13 +503,12 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
! --+>> stiffness calculation <<+-- ! --+>> stiffness calculation <<+--
if(updateJaco) then ! Jacobian required if(updateJaco) then ! Jacobian required
if (debugger) write (6,*) 'Stiffness calculation started'
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains do g = 1,myNgrains
! debugger = (g == 1 .and. i == 1 .and. e == 1)
if (crystallite_converged(g,i,e)) then ! grain converged in above iteration if (crystallite_converged(g,i,e)) then ! grain converged in above iteration
mySizeState = constitutive_sizeState(g,i,e) ! number of state variables for this grain mySizeState = constitutive_sizeState(g,i,e) ! number of state variables for this grain
myState(1:mySizeState) = constitutive_state(g,i,e)%p ! remember unperturbed, converged state... myState(1:mySizeState) = constitutive_state(g,i,e)%p ! remember unperturbed, converged state...
@ -522,40 +520,45 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
myTstar_v = crystallite_Tstar_v(:,g,i,e) myTstar_v = crystallite_Tstar_v(:,g,i,e)
myP = crystallite_P(:,:,g,i,e) myP = crystallite_P(:,:,g,i,e)
if (debugger) then if (debugger) then
!$OMP CRITICAL (write2out)
write (6,*) '#############' write (6,*) '#############'
write (6,*) 'central solution' write (6,*) 'central solution'
write (6,*) '#############' write (6,*) '#############'
write (6,'(a,/,3(3(f12.4,x)/))') ' P of 1 1 1',myP(1:3,:)/1e6 write (6,'(a,/,3(3(f12.4,x)/))') ' P of 1 1 1',myP(1:3,:)/1e6
write (6,'(a,/,3(3(f12.8,x)/))') ' Fp of 1 1 1',myFp(1:3,:) write (6,'(a,/,3(3(f12.8,x)/))') ' Fp of 1 1 1',myFp(1:3,:)
write (6,'(a,/,3(3(f12.8,x)/))') ' Lp of 1 1 1',myLp(1:3,:) write (6,'(a,/,3(3(f12.8,x)/))') ' Lp of 1 1 1',myLp(1:3,:)
write (6,'(a,/,f12.4)') 'state of 1 1 1',myState/1e6 write (6,'(a,/,16(6(e12.4,x)/),2(f12.4,x))') 'state of 1 1 1',myState/1e6
!$OMPEND CRITICAL (write2out)
endif endif
do k = 1,3 ! perturbation... do k = 1,3 ! perturbation...
do l = 1,3 ! ...components do l = 1,3 ! ...components
crystallite_subF(:,:,g,i,e) = myF ! initialize perturbed F to match converged crystallite_subF(:,:,g,i,e) = myF ! initialize perturbed F to match converged
crystallite_subF(k,l,g,i,e) = crystallite_subF(k,l,g,i,e) + pert_Fg ! perturb single component crystallite_subF(k,l,g,i,e) = crystallite_subF(k,l,g,i,e) + pert_Fg ! perturb single component
if (debugger) then if (debugger) then
!$OMP CRITICAL (write2out)
write (6,*) '=============' write (6,*) '============='
write (6,'(i1,x,i1)') k,l write (6,'(i1,x,i1)') k,l
write (6,*) '=============' write (6,*) '============='
write (6,'(a,/,3(3(f12.6,x)/))') 'pertF of 1 1 1',crystallite_subF(1:3,:,g,i,e) write (6,'(a,/,3(3(f12.6,x)/))') 'pertF of 1 1 1',crystallite_subF(1:3,:,g,i,e)
!$OMPEND CRITICAL (write2out)
endif endif
onTrack = .true. onTrack = .true.
converged = .false. converged = .false.
NiterationState = 0_pInt NiterationState = 0_pInt
do while(.not. converged .and. onTrack .and. NiterationState < nState) ! keep cycling until done (potentially non-converged) do while(.not. converged .and. onTrack .and. NiterationState < nState) ! keep cycling until done (potentially non-converged)
NiterationState = NiterationState + 1_pInt NiterationState = NiterationState + 1_pInt
if (debugger) write (6,'(a4,x,i6)') 'loop',NiterationState
onTrack = crystallite_integrateStress(g,i,e) ! stress of perturbed situation (overwrites _P,_Tstar_v,_Fp,_Lp,_Fe) onTrack = crystallite_integrateStress(g,i,e) ! stress of perturbed situation (overwrites _P,_Tstar_v,_Fp,_Lp,_Fe)
if (onTrack) converged = crystallite_updateState(g,i,e).AND.& ! update state if (onTrack) converged = crystallite_updateState(g,i,e).AND.& ! update state
crystallite_updateTemperature(g,i,e) ! update temperature crystallite_updateTemperature(g,i,e) ! update temperature
if (debugger) then if (debugger) then
!$OMP CRITICAL (write2out)
write (6,*) '-------------' write (6,*) '-------------'
write (6,'(l,x,l)') onTrack,converged write (6,'(l,x,l)') onTrack,converged
write (6,'(a,/,3(3(f12.4,x)/))') 'pertP of 1 1 1',crystallite_P(1:3,:,g,i,e)/1e6 write (6,'(a,/,3(3(f12.4,x)/))') 'pertP of 1 1 1',crystallite_P(1:3,:,g,i,e)/1e6
write (6,'(a,/,3(3(f12.4,x)/))') 'DP of 1 1 1',(crystallite_P(1:3,:,g,i,e)-myP(1:3,:))/1e6 write (6,'(a,/,3(3(f12.4,x)/))') 'DP of 1 1 1',(crystallite_P(1:3,:,g,i,e)-myP(1:3,:))/1e6
write (6,'(a,/,f12.4)') 'state of 1 1 1',constitutive_state(g,i,e)%p/1e6 write (6,'(a,/,16(6(e12.4,x)/),/,2(f12.4,x))') 'state of 1 1 1',constitutive_state(g,i,e)%p/1e6
write (6,'(a,/,f12.4)') 'Dstate of 1 1 1',(constitutive_state(g,i,e)%p-myState)/1e6 write (6,'(a,/,16(6(e12.4,x)/),/,2(f12.4,x))') 'Dstate of 1 1 1',(constitutive_state(g,i,e)%p-myState)/1e6
!$OMPEND CRITICAL (write2out)
endif endif
enddo enddo
if (converged) & ! converged state warrants stiffness update if (converged) & ! converged state warrants stiffness update
@ -573,7 +576,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
!$OMPEND CRITICAL (out) !$OMPEND CRITICAL (out)
enddo enddo
enddo enddo
if (debugger) write (6,'(a,/,9(9(f12.4,x)/))') 'dPdF/GPa',crystallite_dPdF(:,:,:,:,1,1,1)/1e9
else ! grain did not converged else ! grain did not converged
crystallite_dPdF(:,:,:,:,g,i,e) = crystallite_fallbackdPdF(:,:,:,:,g,i,e) ! use fallback crystallite_dPdF(:,:,:,:,g,i,e) = crystallite_fallbackdPdF(:,:,:,:,g,i,e) ! use fallback
@ -583,8 +585,6 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
enddo enddo
!$OMPEND PARALLEL DO !$OMPEND PARALLEL DO
if (debugger) write (6,*) 'Stiffness calculation finished'
endif endif
endsubroutine endsubroutine
@ -630,11 +630,6 @@ endsubroutine
tickrate, & tickrate, &
maxticks maxticks
!*** global variables ***!
! crystallite_Tstar_v
! crystallite_subdt
! crystallite_Temperature
mySize = constitutive_sizeDotState(g,i,e) mySize = constitutive_sizeDotState(g,i,e)
! calculate the residuum ! calculate the residuum
@ -649,19 +644,19 @@ endsubroutine
! if NaN occured then return without changing the state ! if NaN occured then return without changing the state
if (any(residuum/=residuum)) then if (any(residuum/=residuum)) then
crystallite_updateState = .false. ! indicate state update failed crystallite_updateState = .false. ! indicate state update failed
if (debugger) write(6,*) '::: updateState encountered NaN' !$OMP CRITICAL (write2out)
write(6,*) '::: updateState encountered NaN',e,i,g
!$OMPEND CRITICAL (write2out)
return return
endif endif
! update the microstructure ! update the microstructure
constitutive_state(g,i,e)%p(1:mySize) = constitutive_state(g,i,e)%p(1:mySize) - residuum constitutive_state(g,i,e)%p(1:mySize) = constitutive_state(g,i,e)%p(1:mySize) - residuum
! setting flag to true if state is below relative Tolerance, otherwise set it to false ! setting flag to true if state is below relative tolerance, otherwise set it to false
crystallite_updateState = maxval(abs(residuum/constitutive_state(g,i,e)%p(1:mySize)), & crystallite_updateState = maxval(abs(residuum/constitutive_state(g,i,e)%p(1:mySize)), &
constitutive_state(g,i,e)%p(1:mySize) /= 0.0_pReal) < rTol_crystalliteState constitutive_state(g,i,e)%p(1:mySize) /= 0.0_pReal) < rTol_crystalliteState
if (debugger) write(6,'(a,/,f12.4)') 'updated state: ', constitutive_state(g,i,e)%p(1)
return return
endfunction endfunction
@ -714,17 +709,21 @@ endsubroutine
! if NaN occured then return without changing the state ! if NaN occured then return without changing the state
if (residuum/=residuum) then if (residuum/=residuum) then
crystallite_updateTemperature = .false. ! indicate update failed crystallite_updateTemperature = .false. ! indicate update failed
if (debugger) write(6,*) '::: updateTemperature encountered NaN' !$OMP CRITICAL (write2out)
write(6,*) '::: updateTemperature encountered NaN',e,i,g
!$OMPEND CRITICAL (write2out)
return return
endif endif
! update the microstructure ! update the microstructure
crystallite_Temperature(g,i,e) = crystallite_Temperature(g,i,e) - residuum crystallite_Temperature(g,i,e) = crystallite_Temperature(g,i,e) - residuum
! setting flag to true if state is below relative Tolerance, otherwise set it to false ! setting flag to true if residuum is below relative tolerance (or zero Kelvin), otherwise set it to false
crystallite_updateTemperature = abs(residuum/crystallite_Temperature(g,i,e)) < rTol_crystalliteTemperature if (crystallite_Temperature(g,i,e) /= 0.0_pReal) then
crystallite_updateTemperature = abs(residuum/crystallite_Temperature(g,i,e)) < rTol_crystalliteTemperature
if (debugger) write(6,'(a,/,f12.4)') 'updated temperature: ', crystallite_Temperature(g,i,e) else
crystallite_updateTemperature = .true.
endif
return return
@ -824,16 +823,6 @@ endsubroutine
tickrate, & tickrate, &
maxticks maxticks
!*** global variables ***!
! crystallite_subF
! crystallite_subFp0
! crystallite_Tstar_v
! crystallite_Lp
! crystallite_subdt
! crystallite_Temperature
if (debugger) write(6,*) '::: integrateStress started'
! be pessimistic ! be pessimistic
crystallite_integrateStress = .false. crystallite_integrateStress = .false.
@ -847,7 +836,9 @@ endsubroutine
! inversion of Fp_current... ! inversion of Fp_current...
invFp_current = math_inv3x3(Fp_current) invFp_current = math_inv3x3(Fp_current)
if (all(invFp_current == 0.0_pReal)) then ! ... failed? if (all(invFp_current == 0.0_pReal)) then ! ... failed?
if (debugger) write(6,*) '::: integrateStress failed on invFp_current inversion' !$OMP CRITICAL (write2out)
write(6,*) '::: integrateStress failed on invFp_current inversion',e,i,g
!$OMPEND CRITICAL (write2out)
return return
endif endif
@ -858,6 +849,7 @@ endsubroutine
! get elasticity tensor ! get elasticity tensor
C_66 = constitutive_homogenizedC(g,i,e) C_66 = constitutive_homogenizedC(g,i,e)
! if (debugger) write(6,'(a,/,6(6(f10.4,x)/))') 'elasticity',C_66(1:6,:)/1e9
C = math_Mandel66to3333(C_66) C = math_Mandel66to3333(C_66)
! start LpLoop with no acceleration ! start LpLoop with no acceleration
@ -872,10 +864,7 @@ LpLoop: do
NiterationStress = NiterationStress + 1 NiterationStress = NiterationStress + 1
! too many loops required ? ! too many loops required ?
if (NiterationStress > nStress) then if (NiterationStress > nStress) return
if (debugger) write(6,*) '::: integrateStress exceeded nStress loopcount'
return
endif
B = math_I3 - crystallite_subdt(g,i,e)*Lpguess B = math_I3 - crystallite_subdt(g,i,e)*Lpguess
BT = transpose(B) BT = transpose(B)
@ -910,7 +899,9 @@ LpLoop: do
! NaN occured at regular speed? ! NaN occured at regular speed?
if (any(residuum/=residuum) .and. leapfrog == 1.0) then if (any(residuum/=residuum) .and. leapfrog == 1.0) then
if (debugger) write(6,*) '::: integrateStress encountered NaN at iteration', NiterationStress !$OMP CRITICAL (write2out)
write(6,*) '::: integrateStress encountered NaN at iteration', NiterationStress,'at',e,i,g
!$OMPEND CRITICAL (write2out)
return return
! something went wrong at accelerated speed? ! something went wrong at accelerated speed?
@ -941,7 +932,11 @@ LpLoop: do
invdRdLp = 0.0_pReal invdRdLp = 0.0_pReal
call math_invert(9,dRdLp,invdRdLp,dummy,error) ! invert dR/dLp --> dLp/dR call math_invert(9,dRdLp,invdRdLp,dummy,error) ! invert dR/dLp --> dLp/dR
if (error) then if (error) then
if (debugger) write(6,*) '::: integrateStress failed on dR/dLp inversion at iteration', NiterationStress if (debugger) then
!$OMP CRITICAL (write2out)
write(6,*) '::: integrateStress failed on dR/dLp inversion at iteration', NiterationStress
!$OMPEND CRITICAL (write2out)
endif
return return
endif endif
endif endif
@ -965,7 +960,11 @@ LpLoop: do
invFp_new = invFp_new/math_det3x3(invFp_new)**(1.0_pReal/3.0_pReal) ! regularize by det invFp_new = invFp_new/math_det3x3(invFp_new)**(1.0_pReal/3.0_pReal) ! regularize by det
call math_invert3x3(invFp_new,Fp_new,det,error) call math_invert3x3(invFp_new,Fp_new,det,error)
if (error) then if (error) then
if (debugger) write(6,*) '::: integrateStress failed on invFp_new inversion at iteration', NiterationStress if (debugger) then
!$OMP CRITICAL (write2out)
write(6,*) '::: integrateStress failed on invFp_new inversion at iteration', NiterationStress
!$OMPEND CRITICAL (write2out)
endif
return return
endif endif
Fe_new = math_mul33x33(Fg_new,invFp_new) ! calc resulting Fe Fe_new = math_mul33x33(Fg_new,invFp_new) ! calc resulting Fe
@ -985,10 +984,12 @@ LpLoop: do
! set return flag to true ! set return flag to true
crystallite_integrateStress = .true. crystallite_integrateStress = .true.
if (debugger) then if (debugger) then
!$OMP CRITICAL (write2out)
write(6,*) '::: integrateStress converged at iteration', NiterationStress write(6,*) '::: integrateStress converged at iteration', NiterationStress
write(6,*) write(6,*)
write(6,'(a,/,3(3(e15.7,x)/))') 'P of 1 1 1',crystallite_P(:,:,1,1,1) write(6,'(a,/,3(3(f12.7,x)/))') 'P / MPa',crystallite_P(:,:,g,i,e)/1e6
write(6,'(a,/,3(3(f12.7,x)/))') 'Lp of 1 1 1',crystallite_Lp(:,:,1,1,1) write(6,'(a,/,3(3(f12.7,x)/))') 'Lp',crystallite_Lp(:,:,g,i,e)
!$OMP CRITICAL (write2out)
endif endif
!$OMP CRITICAL (distributionStress) !$OMP CRITICAL (distributionStress)
@ -1044,10 +1045,6 @@ function crystallite_postResults(&
R R
logical error logical error
!*** global variables ***!
! crystallite_Nresults
! crystallite_Fe
if (crystallite_Nresults >= 2) then if (crystallite_Nresults >= 2) then
crystallite_postResults(1) = material_phase(g,i,e) crystallite_postResults(1) = material_phase(g,i,e)
crystallite_postResults(2) = material_volume(g,i,e) crystallite_postResults(2) = material_volume(g,i,e)

14
code/debug.f90
View File

@ -81,7 +81,7 @@ endsubroutine
if (debug_cumLpCalls > 0_pInt) then if (debug_cumLpCalls > 0_pInt) then
call system_clock(count_rate=tickrate) call system_clock(count_rate=tickrate)
write(6,'(a33,x,f12.6)') 'avg CPU time/microsecs per call :',dble(debug_cumLpTicks)/tickrate/1.0e-6_pReal/debug_cumLpCalls write(6,'(a33,x,f12.6)') 'avg CPU time/microsecs per call :',dble(debug_cumLpTicks)/tickrate/1.0e-6_pReal/debug_cumLpCalls
write(6,'(a33,x,i12)') 'total CPU ticks :',debug_cumLpTicks write(6,'(a33,x,f12.3)') 'total CPU time/s :',dble(debug_cumLpTicks)/tickrate
endif endif
write(6,*) write(6,*)
write(6,'(a33,x,i12)') 'total calls to dotState :',debug_cumDotStateCalls write(6,'(a33,x,i12)') 'total calls to dotState :',debug_cumDotStateCalls
@ -89,7 +89,7 @@ endsubroutine
call system_clock(count_rate=tickrate) call system_clock(count_rate=tickrate)
write(6,'(a33,x,f12.6)') 'avg CPU time/microsecs per call :',& write(6,'(a33,x,f12.6)') 'avg CPU time/microsecs per call :',&
dble(debug_cumDotStateTicks)/tickrate/1.0e-6_pReal/debug_cumDotStateCalls dble(debug_cumDotStateTicks)/tickrate/1.0e-6_pReal/debug_cumDotStateCalls
write(6,'(a33,x,i12)') 'total CPU ticks :',debug_cumDotStateTicks write(6,'(a33,x,f12.3)') 'total CPU time/s :',dble(debug_cumDotStateTicks)/tickrate
endif endif
write(6,*) write(6,*)
write(6,'(a33,x,i12)') 'total calls to dotTemperature :',debug_cumDotTemperatureCalls write(6,'(a33,x,i12)') 'total calls to dotTemperature :',debug_cumDotTemperatureCalls
@ -97,7 +97,7 @@ endsubroutine
call system_clock(count_rate=tickrate) call system_clock(count_rate=tickrate)
write(6,'(a33,x,f12.6)') 'avg CPU time/microsecs per call :',& write(6,'(a33,x,f12.6)') 'avg CPU time/microsecs per call :',&
dble(debug_cumDotTemperatureTicks)/tickrate/1.0e-6_pReal/debug_cumDotTemperatureCalls dble(debug_cumDotTemperatureTicks)/tickrate/1.0e-6_pReal/debug_cumDotTemperatureCalls
write(6,'(a33,x,i12)') 'total CPU ticks :',debug_cumDotTemperatureTicks write(6,'(a33,x,f12.3)') 'total CPU time/s :',dble(debug_cumDotTemperatureTicks)/tickrate
endif endif
integral = 0_pInt integral = 0_pInt
@ -109,7 +109,7 @@ endsubroutine
write(6,'(i25,i10)') i,debug_StressLoopDistribution(i) write(6,'(i25,i10)') i,debug_StressLoopDistribution(i)
endif endif
enddo enddo
write(6,'(a15,i10,i10)') ' total',sum(debug_StressLoopDistribution),integral write(6,'(a15,i10,i10)') ' total',integral,sum(debug_StressLoopDistribution)
integral = 0_pInt integral = 0_pInt
write(6,*) write(6,*)
@ -120,7 +120,7 @@ endsubroutine
write(6,'(i25,i10)') i,debug_StateLoopDistribution(i) write(6,'(i25,i10)') i,debug_StateLoopDistribution(i)
endif endif
enddo enddo
write(6,'(a15,i10,i10)') ' total',sum(debug_StateLoopDistribution),integral write(6,'(a15,i10,i10)') ' total',integral,sum(debug_StateLoopDistribution)
integral = 0_pInt integral = 0_pInt
write(6,*) write(6,*)
@ -131,7 +131,7 @@ endsubroutine
write(6,'(i25,i10)') i,debug_StiffnessStateLoopDistribution(i) write(6,'(i25,i10)') i,debug_StiffnessStateLoopDistribution(i)
endif endif
enddo enddo
write(6,'(a15,i10,i10)') ' total',sum(debug_StiffnessStateLoopDistribution),integral write(6,'(a15,i10,i10)') ' total',integral,sum(debug_StiffnessStateLoopDistribution)
integral = 0_pInt integral = 0_pInt
write(6,*) write(6,*)
@ -142,7 +142,7 @@ endsubroutine
write(6,'(i25,i10)') i,debug_CrystalliteLoopDistribution(i) write(6,'(i25,i10)') i,debug_CrystalliteLoopDistribution(i)
endif endif
enddo enddo
write(6,'(a15,i10,i10)') ' total',sum(debug_CrystalliteLoopDistribution),integral write(6,'(a15,i10,i10)') ' total',integral,sum(debug_CrystalliteLoopDistribution)
write(6,*) write(6,*)
endsubroutine endsubroutine

7
code/homogenization.f90
View File

@ -119,7 +119,7 @@ subroutine homogenization_init(Temperature)
homogenization_maxSizeState = maxval(homogenization_sizeState) homogenization_maxSizeState = maxval(homogenization_sizeState)
homogenization_maxSizePostResults = maxval(homogenization_sizePostResults) homogenization_maxSizePostResults = maxval(homogenization_sizePostResults)
allocate(materialpoint_results( 1+homogenization_maxSizePostResults + & allocate(materialpoint_results( 1+ 1+homogenization_maxSizePostResults + & ! grain count, homogSize, homogResult
homogenization_maxNgrains*(1+crystallite_Nresults+constitutive_maxSizePostResults), mesh_maxNips,mesh_NcpElems)) homogenization_maxNgrains*(1+crystallite_Nresults+constitutive_maxSizePostResults), mesh_maxNips,mesh_NcpElems))
@ -372,7 +372,7 @@ subroutine materialpoint_stressAndItsTangent(&
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
write (6,*) 'Material Point finished' write (6,*) 'Material Point finished'
write (6,'(a,/,3(3(f12.7,x)/))') 'Lp of 1 8 1',crystallite_Lp(1:3,:,1,8,1) write (6,'(a,/,3(3(f12.7,x)/))') 'Lp of 1 1 1',crystallite_Lp(1:3,:,1,1,1)
! how to deal with stiffness? ! how to deal with stiffness?
return return
@ -401,9 +401,10 @@ subroutine materialpoint_postResults(dt)
myNgrains = homogenization_Ngrains(mesh_element(3,e)) myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
c = 0_pInt c = 0_pInt
materialpoint_results(c+1,i,e) = myNgrains; c = c+1_pInt ! tell number of grains at materialpoint
d = homogenization_sizePostResults(i,e) d = homogenization_sizePostResults(i,e)
materialpoint_results(c+1,i,e) = d; c = c+1_pInt ! tell size of homogenization results materialpoint_results(c+1,i,e) = d; c = c+1_pInt ! tell size of homogenization results
if (d > 0_pInt) then ! any homogenization results to mention? if (d > 0_pInt) then ! any homogenization results to mention?
materialpoint_results(c+1:c+d,i,e) = & ! tell homogenization results materialpoint_results(c+1:c+d,i,e) = & ! tell homogenization results
homogenization_postResults(i,e); c = c+d homogenization_postResults(i,e); c = c+d
endif endif

380
code/lattice.f90
View File

@ -20,41 +20,52 @@ implicit none
integer(pInt) lattice_Nhexagonal, & ! # of hexagonal lattice structure (from tag CoverA_ratio) integer(pInt) lattice_Nhexagonal, & ! # of hexagonal lattice structure (from tag CoverA_ratio)
lattice_Nstructure ! # of lattice structures (1: fcc,2: bcc,3+: hexagonal) lattice_Nstructure ! # of lattice structures (1: fcc,2: bcc,3+: hexagonal)
integer(pInt), parameter :: lattice_maxNslipFamily = 4 ! max # of slip system families over lattice structures
integer(pInt), parameter :: lattice_maxNtwinFamily = 4 ! max # of twin system families over lattice structures
integer(pInt), parameter :: lattice_maxNslip = 48 ! max # of slip systems over lattice structures integer(pInt), parameter :: lattice_maxNslip = 48 ! max # of slip systems over lattice structures
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), pointer, dimension(:,:) :: interactionSlipSlip, & integer(pInt), pointer, dimension(:,:) :: interactionSlipSlip, &
interactionSlipTwin, & interactionSlipTwin, &
interactionTwinTwin, & interactionTwinSlip, &
interactionTwinSlip interactionTwinTwin
! Schmid matrices, normal, shear direction and nxd of slip systems ! Schmid matrices, normal, shear direction and d x n of slip systems
real(pReal), allocatable, dimension(:,:,:,:) :: lattice_Sslip real(pReal), allocatable, dimension(:,:,:,:) :: lattice_Sslip
real(pReal), allocatable, dimension(:,:,:) :: lattice_Sslip_v real(pReal), allocatable, dimension(:,:,:) :: lattice_Sslip_v
real(pReal), allocatable, dimension(:,:,:) :: lattice_sn real(pReal), allocatable, dimension(:,:,:) :: lattice_sn, &
real(pReal), allocatable, dimension(:,:,:) :: lattice_sd lattice_sd, &
real(pReal), allocatable, dimension(:,:,:) :: lattice_st lattice_st
! Rotation and Schmid matrices, normal, shear direction and nxd of twin systems ! rotation and Schmid matrices, normal, shear direction and d x n of twin systems
real(pReal), allocatable, dimension(:,:,:,:) :: lattice_Qtwin real(pReal), allocatable, dimension(:,:,:,:) :: lattice_Qtwin
real(pReal), allocatable, dimension(:,:,:,:) :: lattice_Stwin real(pReal), allocatable, dimension(:,:,:,:) :: lattice_Stwin
real(pReal), allocatable, dimension(:,:,:) :: lattice_Stwin_v real(pReal), allocatable, dimension(:,:,:) :: lattice_Stwin_v
real(pReal), allocatable, dimension(:,:,:) :: lattice_tn real(pReal), allocatable, dimension(:,:,:) :: lattice_tn, &
real(pReal), allocatable, dimension(:,:,:) :: lattice_td lattice_td, &
real(pReal), allocatable, dimension(:,:,:) :: lattice_tt lattice_tt
! characteristic twin shear
real(pReal), allocatable, dimension(:,:) :: lattice_shearTwin real(pReal), allocatable, dimension(:,:) :: lattice_shearTwin
integer(pInt), allocatable, dimension(:,:,:) :: lattice_interactionSlipSlip ! number of slip and twin systems in each family
integer(pInt), allocatable, dimension(:,:,:) :: lattice_interactionSlipTwin integer(pInt), allocatable, dimension(:,:) :: lattice_NslipSystem, &
integer(pInt), allocatable, dimension(:,:,:) :: lattice_interactionTwinTwin lattice_NtwinSystem
integer(pInt), allocatable, dimension(:,:,:) :: lattice_interactionTwinSlip
! interaction type of slip and twin systems among each other
integer(pInt), allocatable, dimension(:,:,:) :: lattice_interactionSlipSlip, &
lattice_interactionSlipTwin, &
lattice_interactionTwinSlip, &
lattice_interactionTwinTwin
!============================== fcc (1) ================================= !============================== fcc (1) =================================
integer(pInt), parameter :: lattice_fcc_Nslip = 12_pInt integer(pInt), parameter, dimension(lattice_maxNslipFamily) :: lattice_fcc_NslipSystem = (/12, 0, 0, 0/)
integer(pInt), parameter :: lattice_fcc_Ntwin = 12_pInt integer(pInt), parameter, dimension(lattice_maxNtwinFamily) :: lattice_fcc_NtwinSystem = (/12, 0, 0, 0/)
integer(pInt), parameter :: lattice_fcc_Nslip = 12 ! sum(lattice_fcc_NslipSystem)
integer(pInt), parameter :: lattice_fcc_Ntwin = 12 ! sum(lattice_fcc_NtwinSystem)
integer(pInt) :: lattice_fcc_Nstructure = 0_pInt integer(pInt) :: lattice_fcc_Nstructure = 0_pInt
real(pReal), dimension(3+3,lattice_fcc_Nslip), parameter :: lattice_fcc_systemSlip = & real(pReal), dimension(3+3,lattice_fcc_Nslip), parameter :: lattice_fcc_systemSlip = &
@ -161,8 +172,10 @@ integer(pInt), allocatable, dimension(:,:,:) :: lattice_interactionTwinSlip
!============================== bcc (2) ================================= !============================== bcc (2) =================================
integer(pInt), parameter :: lattice_bcc_Nslip = 48_pInt integer(pInt), parameter, dimension(lattice_maxNslipFamily) :: lattice_bcc_NslipSystem = (/12,12,24, 0/)
integer(pInt), parameter :: lattice_bcc_Ntwin = 12_pInt integer(pInt), parameter, dimension(lattice_maxNtwinFamily) :: lattice_bcc_NtwinSystem = (/12, 0, 0, 0/)
integer(pInt), parameter :: lattice_bcc_Nslip = 48 ! sum(lattice_bcc_NslipSystem)
integer(pInt), parameter :: lattice_bcc_Ntwin = 12 ! sum(lattice_bcc_NtwinSystem)
integer(pInt) :: lattice_bcc_Nstructure = 0_pInt integer(pInt) :: lattice_bcc_Nstructure = 0_pInt
real(pReal), dimension(3+3,lattice_bcc_Nslip), parameter :: lattice_bcc_systemSlip = & real(pReal), dimension(3+3,lattice_bcc_Nslip), parameter :: lattice_bcc_systemSlip = &
@ -352,8 +365,10 @@ integer(pInt), allocatable, dimension(:,:,:) :: lattice_interactionTwinSlip
!============================== hex (3+) ================================= !============================== hex (3+) =================================
integer(pInt), parameter :: lattice_hex_Nslip = 24_pInt integer(pInt), parameter, dimension(lattice_maxNslipFamily) :: lattice_hex_NslipSystem = (/ 3, 3, 6,12/)
integer(pInt), parameter :: lattice_hex_Ntwin = 24_pInt integer(pInt), parameter, dimension(lattice_maxNtwinFamily) :: lattice_hex_NtwinSystem = (/ 6, 6, 6, 6/)
integer(pInt), parameter :: lattice_hex_Nslip = 24 ! sum(lattice_hex_NslipSystem)
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 YJ.Ro and Philip !* sorted by YJ.Ro and Philip
@ -391,29 +406,29 @@ integer(pInt), allocatable, dimension(:,:,:) :: lattice_interactionTwinSlip
real(pReal), dimension(4+4,lattice_hex_Ntwin), parameter :: lattice_hex_systemTwin = & real(pReal), dimension(4+4,lattice_hex_Ntwin), parameter :: lattice_hex_systemTwin = &
reshape((/& reshape((/&
0, 1, -1, 1, 0, -1, 1, 2, & ! <1011>{1012} Twin: shear 0.169 -1.26 tensile 0, 1, -1, 1, 0, -1, 1, 2, & ! <1011>{1012} Twin: shear 0.169 -1.26 compression
-1, 1, 0, 1, 1, -1, 0, 2, & -1, 1, 0, 1, 1, -1, 0, 2, &
-1, 0, 1, 1, 1, 0, -1, 2, & -1, 0, 1, 1, 1, 0, -1, 2, &
0, -1, 1, 1, 0, 1, -1, 2, & 0, -1, 1, 1, 0, 1, -1, 2, &
1, -1, 0, 1, -1, 1, 0, 2, & 1, -1, 0, 1, -1, 1, 0, 2, &
1, 0, -1, 1, -1, 0, 1, 2, & 1, 0, -1, 1, -1, 0, 1, 2, &
2, -1, -1, -3, 2, -1, -1, 2, & ! <211-2>{2112} Twin: shear 0.224 1.19 compressive 2, -1, -1, -3, 2, -1, -1, 2, & ! <211-2>{2112} Twin: shear 0.224 1.19 tension
1, 1, -2, -3, 1, 1, -2, 2, & 1, 1, -2, -3, 1, 1, -2, 2, &
-1, 2, -1, -3, -1, 2, -1, 2, & -1, 2, -1, -3, -1, 2, -1, 2, &
-2, 1, 1, -3, -2, 1, 1, 2, & -2, 1, 1, -3, -2, 1, 1, 2, &
-1, -1, 2, -3, -1, -1, 2, 2, & -1, -1, 2, -3, -1, -1, 2, 2, &
1, -2, 1, -3, 1, -2, 1, 2, & 1, -2, 1, -3, 1, -2, 1, 2, &
-2, 1, 1, 6, 2, -1, -1, 1, & ! <211-6>{2111} Twin: shear 0.628 -0.39 tensile -2, 1, 1, 6, 2, -1, -1, 1, & ! <211-6>{2111} Twin: shear 0.628 -0.39 compression
-1, -1, 2, 6, 1, 1, -2, 1, & -1, -1, 2, 6, 1, 1, -2, 1, &
1, -2, 1, 6, -1, 2, -1, 1, & 1, -2, 1, 6, -1, 2, -1, 1, &
2, -1, -1, 6, -2, 1, 1, 1, & 2, -1, -1, 6, -2, 1, 1, 1, &
1, 1, -2, 6, -1, -1, 2, 1, & 1, 1, -2, 6, -1, -1, 2, 1, &
-1, 2, -1, 6, 1, -2, 1, 1, & -1, 2, -1, 6, 1, -2, 1, 1, &
0, -1, 1, -2, 0, -1, 1, 1, & ! <101-2>{1011} Twin: shear 0.103 1.09 compressive 1, 0, -1, -2, 1, 0, -1, 1, & ! <101-2>{1011} Twin: shear 0.103 1.09 tension
1, -1, 0, -2, 1, -1, 0, 1, &
1, 0, -1, -2, 1, 0, -1, 1, &
0, 1, -1, -2, 0, 1, -1, 1, &
-1, 0, 1, -2, -1, 0, 1, 1, & -1, 0, 1, -2, -1, 0, 1, 1, &
0, 1, -1, -2, 0, 1, -1, 1, &
0, -1, 1, -2, 0, -1, 1, 1, &
1, -1, 0, -2, 1, -1, 0, 1, &
-1, 1, 0, -2, -1, 1, 0, 1 & -1, 1, 0, -2, -1, 1, 0, 1 &
/),(/4+4,lattice_hex_Ntwin/)) !* Sort? Numbering of twin system follows Prof. Tom Bieler's scheme (to be consistent with his work); but numbering in data was restarted from 1 & /),(/4+4,lattice_hex_Ntwin/)) !* Sort? Numbering of twin system follows Prof. Tom Bieler's scheme (to be consistent with his work); but numbering in data was restarted from 1 &
@ -446,131 +461,124 @@ integer(pInt), allocatable, dimension(:,:,:) :: lattice_interactionTwinSlip
/),(/lattice_hex_Ntwin/)) /),(/lattice_hex_Ntwin/))
!* four different interaction type matrix !* four different interaction type matrix
!* 1. slip-slip interaction !* 1. slip-slip interaction - 20 types
!- wihtin the same slip familiy, self or latent hardening; 1 (self), 2 ~ 5 (latent for each slip; bas, pri, pyr_a, pyr_c+a in sequence) !* 2. slip-twin interaction - 16 types
!- interaction between different slip families; 6 ~ 11 !* 3. twin-twin interaction - 20 types
!* 2. slip-twin interaction: current scheme -- indirect effect of twin volume !* 4. twin-slip interaction - 16 types
!- all 1
!* 3. twin-twin interaction
!- wihtin the same slip familiy, self or latent hardening; 1 (self), 2 ~ 5
!- interaction between different slip families; 6 ~ 11
!* 4. twin-slip interaction
!- T1 interacting slip; 1 ~ 4, C1 interacting slip; 5 ~ 8
!- T2 interacting slip; 9 ~ 12, C1 interacting slip; 13 ~ 16
integer(pInt), target, dimension(lattice_hex_Nslip,lattice_hex_Nslip) :: lattice_hex_interactionSlipSlip = & integer(pInt), target, dimension(lattice_hex_Nslip,lattice_hex_Nslip) :: lattice_hex_interactionSlipSlip = &
reshape((/& reshape((/&
1,2,2,6,6,6,7,7,7,7,7,7,8,8,8,8,8,8,8,8,8,8,8,8, & 1, 5, 5, 9, 9, 9,12,12,12,12,12,12,14,14,14,14,14,14,14,14,14,14,14,14, &
2,1,2,6,6,6,7,7,7,7,7,7,8,8,8,8,8,8,8,8,8,8,8,8, & 5, 1, 5, 9, 9, 9,12,12,12,12,12,12,14,14,14,14,14,14,14,14,14,14,14,14, &
2,2,1,6,6,6,7,7,7,7,7,7,8,8,8,8,8,8,8,8,8,8,8,8, & 5, 5, 1, 9, 9, 9,12,12,12,12,12,12,14,14,14,14,14,14,14,14,14,14,14,14, &
6,6,6,1,3,3,9,9,9,9,9,9,10,10,10,10,10,10,10,10,10,10,10,10, & 15,15,15, 2, 6, 6,10,10,10,10,10,10,13,13,13,13,13,13,13,13,13,13,13,13, &
6,6,6,3,1,3,9,9,9,9,9,9,10,10,10,10,10,10,10,10,10,10,10,10, & 15,15,15, 6, 2, 6,10,10,10,10,10,10,13,13,13,13,13,13,13,13,13,13,13,13, &
6,6,6,3,3,1,9,9,9,9,9,9,10,10,10,10,10,10,10,10,10,10,10,10, & 15,15,15, 6, 6, 2,10,10,10,10,10,10,13,13,13,13,13,13,13,13,13,13,13,13, &
7,7,7,9,9,9,1,4,4,4,4,4,11,11,11,11,11,11,11,11,11,11,11,11, & 18,18,18,16,16,16, 3, 7, 7, 7, 7, 7,11,11,11,11,11,11,11,11,11,11,11,11, &
7,7,7,9,9,9,4,1,4,4,4,4,11,11,11,11,11,11,11,11,11,11,11,11, & 18,18,18,16,16,16, 7, 3, 7, 7, 7, 7,11,11,11,11,11,11,11,11,11,11,11,11, &
7,7,7,9,9,9,4,4,1,4,4,4,11,11,11,11,11,11,11,11,11,11,11,11, & 18,18,18,16,16,16, 7, 7, 3, 7, 7, 7,11,11,11,11,11,11,11,11,11,11,11,11, &
7,7,7,9,9,9,4,4,4,1,4,4,11,11,11,11,11,11,11,11,11,11,11,11, & 18,18,18,16,16,16, 7, 7, 7, 3, 7, 7,11,11,11,11,11,11,11,11,11,11,11,11, &
7,7,7,9,9,9,4,4,4,4,1,4,11,11,11,11,11,11,11,11,11,11,11,11, & 18,18,18,16,16,16, 7, 7, 7, 7, 3, 7,11,11,11,11,11,11,11,11,11,11,11,11, &
7,7,7,9,9,9,4,4,4,4,4,1,11,11,11,11,11,11,11,11,11,11,11,11, & 18,18,18,16,16,16, 7, 7, 7, 7, 7, 3,11,11,11,11,11,11,11,11,11,11,11,11, &
8,8,8,10,10,10,11,11,11,11,11,11,1,5,1,5,5,5,5,5,5,5,5,5, & 20,20,20,19,19,19,17,17,17,17,17,17, 4, 8, 4, 8, 8, 8, 8, 8, 8, 8, 8, 8, &
8,8,8,10,10,10,11,11,11,11,11,11,5,1,5,1,5,5,5,5,5,5,5,5, & 20,20,20,19,19,19,17,17,17,17,17,17, 8, 4, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, &
8,8,8,10,10,10,11,11,11,11,11,11,1,5,1,5,5,5,5,5,5,5,5,5, & 20,20,20,19,19,19,17,17,17,17,17,17, 4, 8, 4, 8, 8, 8, 8, 8, 8, 8, 8, 8, &
8,8,8,10,10,10,11,11,11,11,11,11,5,1,5,1,5,5,5,5,5,5,5,5, & 20,20,20,19,19,19,17,17,17,17,17,17, 8, 4, 8, 4, 8, 8, 8, 8, 8, 8, 8, 8, &
8,8,8,10,10,10,11,11,11,11,11,11,5,5,5,5,1,5,5,5,1,5,5,5, & 20,20,20,19,19,19,17,17,17,17,17,17, 8, 8, 8, 8, 4, 8, 8, 8, 4, 8, 8, 8, &
8,8,8,10,10,10,11,11,11,11,11,11,5,5,5,5,5,1,5,1,5,5,5,5, & 20,20,20,19,19,19,17,17,17,17,17,17, 8, 8, 8, 8, 8, 4, 8, 4, 8, 8, 8, 8, &
8,8,8,10,10,10,11,11,11,11,11,11,5,5,5,5,5,5,1,5,5,5,1,5, & 20,20,20,19,19,19,17,17,17,17,17,17, 8, 8, 8, 8, 8, 8, 4, 8, 8, 8, 4, 8, &
8,8,8,10,10,10,11,11,11,11,11,11,5,5,5,5,5,1,5,1,5,5,5,5, & 20,20,20,19,19,19,17,17,17,17,17,17, 8, 8, 8, 8, 8, 4, 8, 4, 8, 8, 8, 8, &
8,8,8,10,10,10,11,11,11,11,11,11,5,5,5,5,1,5,5,5,1,5,5,5, & 20,20,20,19,19,19,17,17,17,17,17,17, 8, 8, 8, 8, 4, 8, 8, 8, 4, 8, 8, 8, &
8,8,8,10,10,10,11,11,11,11,11,11,5,5,5,5,5,5,5,5,5,1,5,1, & 20,20,20,19,19,19,17,17,17,17,17,17, 8, 8, 8, 8, 8, 8, 8, 8, 8, 4, 8, 4, &
8,8,8,10,10,10,11,11,11,11,11,11,5,5,5,5,5,5,1,5,5,5,1,5, & 20,20,20,19,19,19,17,17,17,17,17,17, 8, 8, 8, 8, 8, 8, 4, 8, 8, 8, 4, 8, &
8,8,8,10,10,10,11,11,11,11,11,11,5,5,5,5,5,5,5,5,5,1,5,1 & 20,20,20,19,19,19,17,17,17,17,17,17, 8, 8, 8, 8, 8, 8, 8, 8, 8, 4, 8, 4 &
/),(/lattice_hex_Nslip,lattice_hex_Nslip/)) /),(/lattice_hex_Nslip,lattice_hex_Nslip/))
!* isotropic interaction at the moment !* isotropic interaction at the moment
integer(pInt), target, dimension(lattice_hex_Nslip,lattice_hex_Ntwin) :: lattice_hex_interactionSlipTwin = & integer(pInt), target, dimension(lattice_hex_Nslip,lattice_hex_Ntwin) :: lattice_hex_interactionSlipTwin = &
reshape((/& reshape((/&
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 5, 5, 5, 5, 5, 5, 6, 6, 6, 6, 6, 6, 7, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 8, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 9, 9, 9, 9, 9, 9,10,10,10,10,10,10,11,11,11,11,11,11,12,12,12,12,12,12, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 9, 9, 9, 9, 9, 9,10,10,10,10,10,10,11,11,11,11,11,11,12,12,12,12,12,12, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 9, 9, 9, 9, 9, 9,10,10,10,10,10,10,11,11,11,11,11,11,12,12,12,12,12,12, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 9, 9, 9, 9, 9, 9,10,10,10,10,10,10,11,11,11,11,11,11,12,12,12,12,12,12, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 9, 9, 9, 9, 9, 9,10,10,10,10,10,10,11,11,11,11,11,11,12,12,12,12,12,12, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 9, 9, 9, 9, 9, 9,10,10,10,10,10,10,11,11,11,11,11,11,12,12,12,12,12,12, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 13,13,13,13,13,13,14,14,14,14,14,14,15,15,15,15,15,15,16,16,16,16,16,16, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 13,13,13,13,13,13,14,14,14,14,14,14,15,15,15,15,15,15,16,16,16,16,16,16, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 13,13,13,13,13,13,14,14,14,14,14,14,15,15,15,15,15,15,16,16,16,16,16,16, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 13,13,13,13,13,13,14,14,14,14,14,14,15,15,15,15,15,15,16,16,16,16,16,16, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 13,13,13,13,13,13,14,14,14,14,14,14,15,15,15,15,15,15,16,16,16,16,16,16, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 13,13,13,13,13,13,14,14,14,14,14,14,15,15,15,15,15,15,16,16,16,16,16,16, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 13,13,13,13,13,13,14,14,14,14,14,14,15,15,15,15,15,15,16,16,16,16,16,16, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 13,13,13,13,13,13,14,14,14,14,14,14,15,15,15,15,15,15,16,16,16,16,16,16, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 13,13,13,13,13,13,14,14,14,14,14,14,15,15,15,15,15,15,16,16,16,16,16,16, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 13,13,13,13,13,13,14,14,14,14,14,14,15,15,15,15,15,15,16,16,16,16,16,16, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 13,13,13,13,13,13,14,14,14,14,14,14,15,15,15,15,15,15,16,16,16,16,16,16, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1 & 13,13,13,13,13,13,14,14,14,14,14,14,15,15,15,15,15,15,16,16,16,16,16,16 &
/),(/lattice_hex_Nslip,lattice_hex_Ntwin/)) /),(/lattice_hex_Nslip,lattice_hex_Ntwin/))
integer(pInt), target, dimension(lattice_hex_Ntwin,lattice_hex_Ntwin) :: lattice_hex_interactionTwinTwin = & integer(pInt), target, dimension(lattice_hex_Ntwin,lattice_hex_Ntwin) :: lattice_hex_interactionTwinTwin = &
reshape((/& reshape((/&
1,2,2,2,2,2,6,6,6,6,6,6,7,7,7,7,7,7,8,8,8,8,8,8, & 1, 5, 5, 5, 5, 5, 9, 9, 9, 9, 9, 9,12,12,12,12,12,12,14,14,14,14,14,14, &
2,1,2,2,2,2,6,6,6,6,6,6,7,7,7,7,7,7,8,8,8,8,8,8, & 5, 1, 5, 5, 5, 5, 9, 9, 9, 9, 9, 9,12,12,12,12,12,12,14,14,14,14,14,14, &
2,2,1,2,2,2,6,6,6,6,6,6,7,7,7,7,7,7,8,8,8,8,8,8, & 5, 5, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9,12,12,12,12,12,12,14,14,14,14,14,14, &
2,2,2,1,2,2,6,6,6,6,6,6,7,7,7,7,7,7,8,8,8,8,8,8, & 5, 5, 5, 1, 5, 5, 9, 9, 9, 9, 9, 9,12,12,12,12,12,12,14,14,14,14,14,14, &
2,2,2,2,1,2,6,6,6,6,6,6,7,7,7,7,7,7,8,8,8,8,8,8, & 5, 5, 5, 5, 1, 5, 9, 9, 9, 9, 9, 9,12,12,12,12,12,12,14,14,14,14,14,14, &
2,2,2,2,2,1,6,6,6,6,6,6,7,7,7,7,7,7,8,8,8,8,8,8, & 5, 5, 5, 5, 5, 1, 9, 9, 9, 9, 9, 9,12,12,12,12,12,12,14,14,14,14,14,14, &
6,6,6,6,6,6,1,3,3,3,3,3,9,9,9,9,9,9,10,10,10,10,10,10, & 15,15,15,15,15,15, 2, 6, 6, 6, 6, 6,10,10,10,10,10,10,13,13,13,13,13,13, &
6,6,6,6,6,6,3,1,3,3,3,3,9,9,9,9,9,9,10,10,10,10,10,10, & 15,15,15,15,15,15, 6, 2, 6, 6, 6, 6,10,10,10,10,10,10,13,13,13,13,13,13, &
6,6,6,6,6,6,3,3,1,3,3,3,9,9,9,9,9,9,10,10,10,10,10,10, & 15,15,15,15,15,15, 6, 6, 2, 6, 6, 6,10,10,10,10,10,10,13,13,13,13,13,13, &
6,6,6,6,6,6,3,3,3,1,3,3,9,9,9,9,9,9,10,10,10,10,10,10, & 15,15,15,15,15,15, 6, 6, 6, 2, 6, 6,10,10,10,10,10,10,13,13,13,13,13,13, &
6,6,6,6,6,6,3,3,3,3,1,3,9,9,9,9,9,9,10,10,10,10,10,10, & 15,15,15,15,15,15, 6, 6, 6, 6, 2, 6,10,10,10,10,10,10,13,13,13,13,13,13, &
6,6,6,6,6,6,3,3,3,3,3,1,9,9,9,9,9,9,10,10,10,10,10,10, & 15,15,15,15,15,15, 6, 6, 6, 6, 6, 2,10,10,10,10,10,10,13,13,13,13,13,13, &
7,7,7,7,7,7,9,9,9,9,9,9,1,4,4,4,4,4,11,11,11,11,11,11, & 18,18,18,18,18,18,16,16,16,16,16,16, 3, 7, 7, 7, 7, 7,11,11,11,11,11,11, &
7,7,7,7,7,7,9,9,9,9,9,9,4,1,4,4,4,4,11,11,11,11,11,11, & 18,18,18,18,18,18,16,16,16,16,16,16, 7, 3, 7, 7, 7, 7,11,11,11,11,11,11, &
7,7,7,7,7,7,9,9,9,9,9,9,4,4,1,4,4,4,11,11,11,11,11,11, & 18,18,18,18,18,18,16,16,16,16,16,16, 7, 7, 3, 7, 7, 7,11,11,11,11,11,11, &
7,7,7,7,7,7,9,9,9,9,9,9,4,4,4,1,4,4,11,11,11,11,11,11, & 18,18,18,18,18,18,16,16,16,16,16,16, 7, 7, 7, 3, 7, 7,11,11,11,11,11,11, &
7,7,7,7,7,7,9,9,9,9,9,9,4,4,4,4,1,4,11,11,11,11,11,11, & 18,18,18,18,18,18,16,16,16,16,16,16, 7, 7, 7, 7, 3, 7,11,11,11,11,11,11, &
7,7,7,7,7,7,9,9,9,9,9,9,4,4,4,4,4,1,11,11,11,11,11,11, & 18,18,18,18,18,18,16,16,16,16,16,16, 7, 7, 7, 7, 7, 3,11,11,11,11,11,11, &
8,8,8,8,8,8,10,10,10,10,10,10,11,11,11,11,11,11,1,5,5,5,5,5, & 20,20,20,20,20,20,19,19,19,19,19,19,17,17,17,17,17,17, 4, 8, 8, 8, 8, 8, &
8,8,8,8,8,8,10,10,10,10,10,10,11,11,11,11,11,11,5,1,5,5,5,5, & 20,20,20,20,20,20,19,19,19,19,19,19,17,17,17,17,17,17, 8, 4, 8, 8, 8, 8, &
8,8,8,8,8,8,10,10,10,10,10,10,11,11,11,11,11,11,5,5,1,5,5,5, & 20,20,20,20,20,20,19,19,19,19,19,19,17,17,17,17,17,17, 8, 8, 4, 8, 8, 8, &
8,8,8,8,8,8,10,10,10,10,10,10,11,11,11,11,11,11,5,5,5,1,5,5, & 20,20,20,20,20,20,19,19,19,19,19,19,17,17,17,17,17,17, 8, 8, 8, 4, 8, 8, &
8,8,8,8,8,8,10,10,10,10,10,10,11,11,11,11,11,11,5,5,5,5,1,5, & 20,20,20,20,20,20,19,19,19,19,19,19,17,17,17,17,17,17, 8, 8, 8, 8, 4, 8, &
8,8,8,8,8,8,10,10,10,10,10,10,11,11,11,11,11,11,5,5,5,5,5,1 & 20,20,20,20,20,20,19,19,19,19,19,19,17,17,17,17,17,17, 8, 8, 8, 8, 8, 4 &
/),(/lattice_hex_Ntwin,lattice_hex_Ntwin/)) /),(/lattice_hex_Ntwin,lattice_hex_Ntwin/))
!* isotropic interaction at the moment !* isotropic interaction at the moment
integer(pInt), target, dimension(lattice_hex_Ntwin,lattice_hex_Nslip) :: lattice_hex_interactionTwinSlip = & integer(pInt), target, dimension(lattice_hex_Ntwin,lattice_hex_Nslip) :: lattice_hex_interactionTwinSlip = &
reshape((/& reshape((/&
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 1, 1, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9,13,13,13,13,13,13,13,13,13,13,13,13, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 1, 1, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9,13,13,13,13,13,13,13,13,13,13,13,13, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 1, 1, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9,13,13,13,13,13,13,13,13,13,13,13,13, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 1, 1, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9,13,13,13,13,13,13,13,13,13,13,13,13, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 1, 1, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9,13,13,13,13,13,13,13,13,13,13,13,13, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 1, 1, 1, 5, 5, 5, 9, 9, 9, 9, 9, 9,13,13,13,13,13,13,13,13,13,13,13,13, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 2, 2, 2, 6, 6, 6,10,10,10,10,10,10,14,14,14,14,14,14,14,14,14,14,14,14, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 2, 2, 2, 6, 6, 6,10,10,10,10,10,10,14,14,14,14,14,14,14,14,14,14,14,14, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 2, 2, 2, 6, 6, 6,10,10,10,10,10,10,14,14,14,14,14,14,14,14,14,14,14,14, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 2, 2, 2, 6, 6, 6,10,10,10,10,10,10,14,14,14,14,14,14,14,14,14,14,14,14, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 2, 2, 2, 6, 6, 6,10,10,10,10,10,10,14,14,14,14,14,14,14,14,14,14,14,14, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 2, 2, 2, 6, 6, 6,10,10,10,10,10,10,14,14,14,14,14,14,14,14,14,14,14,14, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 3, 3, 3, 7, 7, 7,11,11,11,11,11,11,15,15,15,15,15,15,15,15,15,15,15,15, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 3, 3, 3, 7, 7, 7,11,11,11,11,11,11,15,15,15,15,15,15,15,15,15,15,15,15, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 3, 3, 3, 7, 7, 7,11,11,11,11,11,11,15,15,15,15,15,15,15,15,15,15,15,15, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 3, 3, 3, 7, 7, 7,11,11,11,11,11,11,15,15,15,15,15,15,15,15,15,15,15,15, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 3, 3, 3, 7, 7, 7,11,11,11,11,11,11,15,15,15,15,15,15,15,15,15,15,15,15, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 3, 3, 3, 7, 7, 7,11,11,11,11,11,11,15,15,15,15,15,15,15,15,15,15,15,15, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 3, 3, 3, 7, 7, 7,11,11,11,11,11,11,15,15,15,15,15,15,15,15,15,15,15,15, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 4, 4, 4, 8, 8, 8,12,12,12,12,12,12,16,16,16,16,16,16,16,16,16,16,16,16, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 4, 4, 4, 8, 8, 8,12,12,12,12,12,12,16,16,16,16,16,16,16,16,16,16,16,16, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 4, 4, 4, 8, 8, 8,12,12,12,12,12,12,16,16,16,16,16,16,16,16,16,16,16,16, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1, & 4, 4, 4, 8, 8, 8,12,12,12,12,12,12,16,16,16,16,16,16,16,16,16,16,16,16, &
1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1 & 4, 4, 4, 8, 8, 8,12,12,12,12,12,12,16,16,16,16,16,16,16,16,16,16,16,16 &
/),(/lattice_hex_Ntwin,lattice_hex_Nslip/)) /),(/lattice_hex_Ntwin,lattice_hex_Nslip/))
@ -585,42 +593,46 @@ subroutine lattice_init()
!************************************** !**************************************
!* Module initialization * !* Module initialization *
!************************************** !**************************************
use IO, only: IO_open_file,IO_countSections,IO_countTagInPart,IO_error use IO, only: IO_open_file,IO_countSections,IO_countTagInPart,IO_error
use material, only: material_configfile,material_partPhase use material, only: material_configfile,material_partPhase
implicit none implicit none
integer(pInt), parameter :: fileunit = 200 integer(pInt), parameter :: fileunit = 200
integer(pInt) i,Nsections integer(pInt) i,Nsections
write(6,*) write(6,*)
write(6,*) '<<<+- lattice init -+>>>' write(6,*) '<<<+- lattice init -+>>>'
write(6,*) write(6,*)
if(.not. IO_open_file(fileunit,material_configFile)) call IO_error (100) ! corrupt config file if(.not. IO_open_file(fileunit,material_configFile)) call IO_error (100) ! corrupt config file
Nsections = IO_countSections(fileunit,material_partPhase) Nsections = IO_countSections(fileunit,material_partPhase)
lattice_Nstructure = 2_pInt + sum(IO_countTagInPart(fileunit,material_partPhase,'covera_ratio',Nsections)) ! fcc + bcc + all hex lattice_Nstructure = 2_pInt + sum(IO_countTagInPart(fileunit,material_partPhase,'covera_ratio',Nsections)) ! fcc + bcc + all hex
close(fileunit) close(fileunit)
allocate(lattice_Sslip(3,3,lattice_maxNslip,lattice_Nstructure)); lattice_Sslip = 0.0_pReal allocate(lattice_Sslip(3,3,lattice_maxNslip,lattice_Nstructure)); lattice_Sslip = 0.0_pReal
allocate(lattice_Sslip_v(6,lattice_maxNslip,lattice_Nstructure)); lattice_Sslip_v = 0.0_pReal allocate(lattice_Sslip_v(6,lattice_maxNslip,lattice_Nstructure)); lattice_Sslip_v = 0.0_pReal
allocate(lattice_sd(3,lattice_maxNslip,lattice_Nstructure)); lattice_sd = 0.0_pReal allocate(lattice_sd(3,lattice_maxNslip,lattice_Nstructure)); lattice_sd = 0.0_pReal
allocate(lattice_st(3,lattice_maxNslip,lattice_Nstructure)); lattice_st = 0.0_pReal allocate(lattice_st(3,lattice_maxNslip,lattice_Nstructure)); lattice_st = 0.0_pReal
allocate(lattice_sn(3,lattice_maxNslip,lattice_Nstructure)); lattice_sn = 0.0_pReal allocate(lattice_sn(3,lattice_maxNslip,lattice_Nstructure)); lattice_sn = 0.0_pReal
allocate(lattice_Qtwin(3,3,lattice_maxNtwin,lattice_Nstructure)); lattice_Qtwin = 0.0_pReal allocate(lattice_Qtwin(3,3,lattice_maxNtwin,lattice_Nstructure)); lattice_Qtwin = 0.0_pReal
allocate(lattice_Stwin(3,3,lattice_maxNtwin,lattice_Nstructure)); lattice_Stwin = 0.0_pReal allocate(lattice_Stwin(3,3,lattice_maxNtwin,lattice_Nstructure)); lattice_Stwin = 0.0_pReal
allocate(lattice_Stwin_v(6,lattice_maxNtwin,lattice_Nstructure)); lattice_Stwin_v = 0.0_pReal allocate(lattice_Stwin_v(6,lattice_maxNtwin,lattice_Nstructure)); lattice_Stwin_v = 0.0_pReal
allocate(lattice_td(3,lattice_maxNtwin,lattice_Nstructure)); lattice_td = 0.0_pReal allocate(lattice_td(3,lattice_maxNtwin,lattice_Nstructure)); lattice_td = 0.0_pReal
allocate(lattice_tt(3,lattice_maxNtwin,lattice_Nstructure)); lattice_tt = 0.0_pReal allocate(lattice_tt(3,lattice_maxNtwin,lattice_Nstructure)); lattice_tt = 0.0_pReal
allocate(lattice_tn(3,lattice_maxNtwin,lattice_Nstructure)); lattice_tn = 0.0_pReal allocate(lattice_tn(3,lattice_maxNtwin,lattice_Nstructure)); lattice_tn = 0.0_pReal
allocate(lattice_shearTwin(lattice_maxNtwin,lattice_Nstructure)); lattice_shearTwin = 0.0_pReal
allocate(lattice_interactionSlipSlip(lattice_maxNslip,lattice_maxNslip,lattice_Nstructure)); lattice_interactionSlipSlip = 0_pInt allocate(lattice_shearTwin(lattice_maxNtwin,lattice_Nstructure)); lattice_shearTwin = 0.0_pReal
allocate(lattice_interactionSlipTwin(lattice_maxNslip,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionSlipTwin = 0_pInt
allocate(lattice_interactionTwinTwin(lattice_maxNtwin,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionTwinTwin = 0_pInt
allocate(lattice_interactionTwinSlip(lattice_maxNtwin,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionTwinSlip = 0_pInt
endsubroutine allocate(lattice_NslipSystem(lattice_maxNslipFamily,lattice_Nstructure)); lattice_NslipSystem = 0.0_pReal
allocate(lattice_NtwinSystem(lattice_maxNtwinFamily,lattice_Nstructure)); lattice_NtwinSystem = 0.0_pReal
allocate(lattice_interactionSlipSlip(lattice_maxNslip,lattice_maxNslip,lattice_Nstructure)); lattice_interactionSlipSlip = 0_pInt
allocate(lattice_interactionSlipTwin(lattice_maxNslip,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionSlipTwin = 0_pInt
allocate(lattice_interactionTwinTwin(lattice_maxNtwin,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionTwinTwin = 0_pInt
allocate(lattice_interactionTwinSlip(lattice_maxNtwin,lattice_maxNtwin,lattice_Nstructure)); lattice_interactionTwinSlip = 0_pInt
end subroutine
function lattice_initializeStructure(struct,CoverA) function lattice_initializeStructure(struct,CoverA)
@ -643,14 +655,20 @@ function lattice_initializeStructure(struct,CoverA)
real(pReal), dimension(lattice_maxNtwin) :: ts = 0.0_pReal real(pReal), dimension(lattice_maxNtwin) :: ts = 0.0_pReal
real(pReal), dimension(3) :: hex_d = 0.0_pReal, & real(pReal), dimension(3) :: hex_d = 0.0_pReal, &
hex_n = 0.0_pReal hex_n = 0.0_pReal
integer(pInt), dimension(lattice_maxNslipFamily) :: myNslipSystem = 0_pInt
integer(pInt), dimension(lattice_maxNtwinFamily) :: myNtwinSystem = 0_pInt
integer(pInt) :: i,myNslip,myNtwin,myStructure = 0_pInt integer(pInt) :: i,myNslip,myNtwin,myStructure = 0_pInt
logical :: processMe = .false. logical :: processMe
integer(pInt) lattice_initializeStructure integer(pInt) lattice_initializeStructure
processMe = .false.
select case(struct(1:3)) ! check first three chars of structure name select case(struct(1:3)) ! check first three chars of structure name
case ('fcc') case ('fcc')
myStructure = 1_pInt myStructure = 1_pInt
myNslipSystem = lattice_fcc_NslipSystem
myNtwinSystem = lattice_fcc_NtwinSystem
myNslip = lattice_fcc_Nslip myNslip = lattice_fcc_Nslip
myNtwin = lattice_fcc_Ntwin myNtwin = lattice_fcc_Ntwin
lattice_fcc_Nstructure = lattice_fcc_Nstructure + 1_pInt lattice_fcc_Nstructure = lattice_fcc_Nstructure + 1_pInt
@ -675,6 +693,8 @@ function lattice_initializeStructure(struct,CoverA)
case ('bcc') case ('bcc')
myStructure = 2_pInt myStructure = 2_pInt
myNslipSystem = lattice_bcc_NslipSystem
myNtwinSystem = lattice_bcc_NtwinSystem
myNslip = lattice_bcc_Nslip myNslip = lattice_bcc_Nslip
myNtwin = lattice_bcc_Ntwin myNtwin = lattice_bcc_Ntwin
lattice_bcc_Nstructure = lattice_bcc_Nstructure + 1_pInt lattice_bcc_Nstructure = lattice_bcc_Nstructure + 1_pInt
@ -701,29 +721,31 @@ function lattice_initializeStructure(struct,CoverA)
if (CoverA > 0.0_pReal) then if (CoverA > 0.0_pReal) then
lattice_hex_Nstructure = lattice_hex_Nstructure + 1_pInt lattice_hex_Nstructure = lattice_hex_Nstructure + 1_pInt
myStructure = 2_pInt + lattice_hex_Nstructure myStructure = 2_pInt + lattice_hex_Nstructure
myNslipSystem = lattice_hex_NslipSystem
myNtwinSystem = lattice_hex_NtwinSystem
myNslip = lattice_hex_Nslip myNslip = lattice_hex_Nslip
myNtwin = lattice_hex_Ntwin myNtwin = lattice_hex_Ntwin
processMe = .true. processMe = .true.
! converting from 4 axes coordinate system (a1=a2=a3=c) to ortho-hexgonal system (a, b, c) ! converting from 4 axes coordinate system (a1=a2=a3=c) to ortho-hexgonal system (a, b, c)
do i = 1,myNslip do i = 1,myNslip
hex_d(1) = lattice_hex_systemSlip(1,i) ! direction [uvtw]->[3u/2 (u+2v)*sqrt(3)/2 w*(c/a)] hex_d(1) = lattice_hex_systemSlip(1,i)*1.5_pReal ! direction [uvtw]->[3u/2 (u+2v)*sqrt(3)/2 w*(c/a)]
hex_d(2) = (lattice_hex_systemSlip(1,i)+2.0_pReal*lattice_hex_systemSlip(2,i))/dsqrt(3.0_pReal) hex_d(2) = (lattice_hex_systemSlip(1,i)+2.0_pReal*lattice_hex_systemSlip(2,i))*(0.5_pReal*dsqrt(3.0_pReal))
hex_d(3) = lattice_hex_systemSlip(4,i)/CoverA hex_d(3) = lattice_hex_systemSlip(4,i)*CoverA
hex_n(1) = lattice_hex_systemSlip(5,i)*1.5_pReal ! plane (hkil)->(h (h+2k)/sqrt(3) l/(c/a)) hex_n(1) = lattice_hex_systemSlip(5,i) ! plane (hkil)->(h (h+2k)/sqrt(3) l/(c/a))
hex_n(2) = (lattice_hex_systemSlip(5,i)+2.0_pReal*lattice_hex_systemSlip(6,i))*(0.5_pReal*dsqrt(3.0_pReal)) hex_n(2) = (lattice_hex_systemSlip(5,i)+2.0_pReal*lattice_hex_systemSlip(6,i))/dsqrt(3.0_pReal)
hex_n(3) = lattice_hex_systemSlip(8,i)*CoverA hex_n(3) = lattice_hex_systemSlip(8,i)/CoverA
sd(:,i) = hex_d/dsqrt(math_mul3x3(hex_d,hex_d)) sd(:,i) = hex_d/dsqrt(math_mul3x3(hex_d,hex_d))
sn(:,i) = hex_n/dsqrt(math_mul3x3(hex_n,hex_n)) sn(:,i) = hex_n/dsqrt(math_mul3x3(hex_n,hex_n))
st(:,i) = math_vectorproduct(sd(:,i),sn(:,i)) st(:,i) = math_vectorproduct(sd(:,i),sn(:,i))
enddo enddo
do i = 1,myNtwin do i = 1,myNtwin
hex_d(1) = lattice_hex_systemTwin(1,i) hex_d(1) = lattice_hex_systemTwin(1,i)*1.5_pReal
hex_d(2) = (lattice_hex_systemTwin(1,i)+2.0_pReal*lattice_hex_systemTwin(2,i))/dsqrt(3.0_pReal) hex_d(2) = (lattice_hex_systemTwin(1,i)+2.0_pReal*lattice_hex_systemTwin(2,i))*(0.5_pReal*dsqrt(3.0_pReal))
hex_d(3) = lattice_hex_systemTwin(4,i)/CoverA hex_d(3) = lattice_hex_systemTwin(4,i)*CoverA
hex_n(1) = lattice_hex_systemTwin(5,i)*1.5_pReal hex_n(1) = lattice_hex_systemTwin(4,i)
hex_n(2) = (lattice_hex_systemTwin(5,i)+2.0_pReal*lattice_hex_systemTwin(6,i))*(0.5_pReal*dsqrt(3.0_pReal)) hex_n(2) = (lattice_hex_systemTwin(4,i)+2.0_pReal*lattice_hex_systemTwin(6,i))/dsqrt(3.0_pReal)
hex_n(3) = lattice_hex_systemTwin(8,i)*CoverA hex_n(3) = lattice_hex_systemTwin(8,i)/CoverA
td(:,i) = hex_d/dsqrt(math_mul3x3(hex_d,hex_d)) td(:,i) = hex_d/dsqrt(math_mul3x3(hex_d,hex_d))
tn(:,i) = hex_n/dsqrt(math_mul3x3(hex_n,hex_n)) tn(:,i) = hex_n/dsqrt(math_mul3x3(hex_n,hex_n))
@ -754,15 +776,17 @@ function lattice_initializeStructure(struct,CoverA)
lattice_Qtwin(:,:,i,myStructure) = math_RodrigToR(tn(:,i),180.0_pReal*inRad) lattice_Qtwin(:,:,i,myStructure) = math_RodrigToR(tn(:,i),180.0_pReal*inRad)
lattice_shearTwin(i,myStructure) = ts(i) lattice_shearTwin(i,myStructure) = ts(i)
enddo enddo
lattice_NslipSystem(1:lattice_maxNslipFamily,myStructure) = myNslipSystem ! number of slip systems in each family
lattice_NtwinSystem(1:lattice_maxNtwinFamily,myStructure) = myNtwinSystem ! number of twin systems in each family
lattice_interactionSlipSlip(1:myNslip,1:myNslip,myStructure) = interactionSlipSlip(1:myNslip,1:myNslip) lattice_interactionSlipSlip(1:myNslip,1:myNslip,myStructure) = interactionSlipSlip(1:myNslip,1:myNslip)
lattice_interactionSlipTwin(1:myNslip,1:myNtwin,myStructure) = interactionSlipTwin(1:myNslip,1:myNtwin) lattice_interactionSlipTwin(1:myNslip,1:myNtwin,myStructure) = interactionSlipTwin(1:myNslip,1:myNtwin)
lattice_interactionTwinTwin(1:myNtwin,1:myNtwin,myStructure) = interactionTwinTwin(1:myNtwin,1:myNtwin)
lattice_interactionTwinSlip(1:myNtwin,1:myNslip,myStructure) = interactionTwinSlip(1:myNtwin,1:myNslip) lattice_interactionTwinSlip(1:myNtwin,1:myNslip,myStructure) = interactionTwinSlip(1:myNtwin,1:myNslip)
lattice_interactionTwinTwin(1:myNtwin,1:myNtwin,myStructure) = interactionTwinTwin(1:myNtwin,1:myNtwin)
endif endif
lattice_initializeStructure = myStructure lattice_initializeStructure = myStructure
endfunction end function
END MODULE END MODULE

106
code/material.config
View File

@ -61,58 +61,70 @@ Ngrains 100
<phase> <phase>
##################### #####################
[Aluminum] # below given format will not work. need to select one constitution block from it. [Aluminum_J2isotropic]
constitution j2
(output) flowstress
(output) strainrate
c11 110.9e9
c12 58.34e9
taylorfactor 3
tau0 31e6
gdot0 0.001
n 20
h0 75e6
tausat 63e6
w0 1
[Aluminum_phenopowerlaw]
# slip only
constitution phenopowerlaw
(output) resistance_slip
(output) shearrate_slip
(output) resolvedstress_slip
(output) totalshear
(output) resistance_twin
(output) shearrate_twin
(output) resolvedstress_twin
(output) totalvolfrac
lattice_structure fcc
Nslip 12 0 0 0 # per family
Ntwin 0 0 0 0 # per family
c11 106.75e9 c11 106.75e9
c12 60.41e9 c12 60.41e9
c44 28.34e9 c44 28.34e9
lattice_structure fcc
Nslip 12
Ntwin 12
constitution phenomenological
(output) slipResistance
(output) rateOfShear
### phenomenological constitutive parameters ###
s0_slip 31e6
gdot0_slip 0.001 gdot0_slip 0.001
n_slip 20 n_slip 20
h0 75e6 s0_slip 31e6 0 0 0 # per family
s_sat 63e6 ssat_slip 63e6 0 0 0 # per family
w0 2.25 gdot0_twin 0.001
latent_ratio 1.4 n_twin 20
s0_twin 31e6 0 0 0 # per family
s_pr 0 # push-up factor for slip saturation due to twinning
twin_b 0
twin_c 0
twin_d 0
twin_e 0
h0_slipslip 75e6
h0_sliptwin 0
h0_twinslip 0
h0_twintwin 0
interaction_slipslip 1 1 1.4 1.4 1.4 1.4
interaction_sliptwin 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
interaction_twinslip 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
interaction_twintwin 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
#constitution j2
#c11 110.9e9
#c12 58.34e9
#(output) flowstress
#(output) strainrate
#taylorfactor 3
#s0 31e6
#gdot0 0.001
#n 20
#h0 75e6
#s_sat 63e6
#w0 2.25
[Copper]
[Ferrite]
[Martensite]
[TWIP steel FeMnC] [TWIP steel FeMnC]
C11 175.0e9 # elastic constants in Pa
C12 115.0e9
C44 135.0e9
lattice_structure fcc
Nslip 12
Ntwin 12
#constitution phenomenological
#(output) slipResistance
#(output) rateOfShear
constitution dislobased constitution dislobased
(output) state_slip (output) state_slip
(output) rateofshear_slip (output) rateofshear_slip
(output) mfp_slip (output) mfp_slip
@ -122,14 +134,12 @@ constitution dislobased
(output) mfp_twin (output) mfp_twin
(output) thresholdstress_twin (output) thresholdstress_twin
### phenomenological constitutive parameters ### C11 175.0e9 # elastic constants in Pa
#s0_slip 85.0e6 # initial slip resistance C12 115.0e9
#gdot0_slip 0.001 # reference shear rate C44 135.0e9
#n_slip 100.0 # stress exponent lattice_structure fcc
#h0 355.0e6 # initial hardening slope Nslip 12
#s_sat 265.0e6 # saturation stress Ntwin 12
#w0 1.0 # exponent
#latent_ratio 1.4 # latent/self hardening ratio
### dislocation density-based constitutive parameters ### ### dislocation density-based constitutive parameters ###
burgers 2.56e-10 # Burgers vector [m] burgers 2.56e-10 # Burgers vector [m]

38
code/material.f90
View File

@ -28,7 +28,7 @@ integer(pInt) material_Nhomogenization, &
material_Nphase, & ! number of phases material_Nphase, & ! number of phases
material_Ntexture, & ! number of textures material_Ntexture, & ! number of textures
microstructure_maxNconstituents, & ! max number of constituents in any phase microstructure_maxNconstituents, & ! max number of constituents in any phase
homogenization_maxNgrains, & ! max number of grains in any homogenization homogenization_maxNgrains, & ! max number of grains in any USED homogenization
texture_maxNgauss, & ! max number of Gauss components in any texture texture_maxNgauss, & ! max number of Gauss components in any texture
texture_maxNfiber ! max number of Fiber components in any texture texture_maxNfiber ! max number of Fiber components in any texture
character(len=64), dimension(:), allocatable :: homogenization_name, & ! name of each homogenization character(len=64), dimension(:), allocatable :: homogenization_name, & ! name of each homogenization
@ -44,10 +44,12 @@ integer(pInt), dimension(:), allocatable :: homogenization_Ngrains, &
microstructure_Nconstituents, & ! number of constituents in each microstructure microstructure_Nconstituents, & ! number of constituents in each microstructure
phase_constitutionInstance, & ! instance of particular constitution of each phase phase_constitutionInstance, & ! instance of particular constitution of each phase
phase_Noutput, & ! number of '(output)' items per phase phase_Noutput, & ! number of '(output)' items per phase
phase_localConstitution, & ! flag phases with local constitutive law
texture_symmetry, & ! number of symmetric orientations per texture texture_symmetry, & ! number of symmetric orientations per texture
texture_Ngauss, & ! number of Gauss components per texture texture_Ngauss, & ! number of Gauss components per texture
texture_Nfiber ! number of Fiber components per texture texture_Nfiber ! number of Fiber components per texture
logical, dimension(:), allocatable :: homogenization_active, & !
microstructure_active, & !
phase_localConstitution ! flags phases with local constitutive law
integer(pInt), dimension(:,:), allocatable :: microstructure_phase, & ! phase IDs of each microstructure integer(pInt), dimension(:,:), allocatable :: microstructure_phase, & ! phase IDs of each microstructure
microstructure_texture ! texture IDs of each microstructure microstructure_texture ! texture IDs of each microstructure
real(pReal), dimension(:,:), allocatable :: microstructure_fraction ! vol fraction of each constituent in microstructure real(pReal), dimension(:,:), allocatable :: microstructure_fraction ! vol fraction of each constituent in microstructure
@ -71,7 +73,7 @@ subroutine material_init()
!* Definition of variables !* Definition of variables
integer(pInt), parameter :: fileunit = 200 integer(pInt), parameter :: fileunit = 200
integer(pInt) i integer(pInt) i,j
write(6,*) write(6,*)
write(6,*) '<<<+- material init -+>>>' write(6,*) '<<<+- material init -+>>>'
@ -94,14 +96,22 @@ subroutine material_init()
write (6,*) write (6,*)
write (6,*) 'MATERIAL configuration' write (6,*) 'MATERIAL configuration'
write (6,*) write (6,*)
write (6,*) 'Homogenization' write (6,'(a32,x,a16,x,a6)') 'homogenization ','type ','grains'
do i = 1,material_Nhomogenization do i = 1,material_Nhomogenization
write (6,'(a32,x,a16,x,i4)') homogenization_name(i),homogenization_type(i),homogenization_Ngrains(i) write (6,'(x,a32,x,a16,x,i4)') homogenization_name(i),homogenization_type(i),homogenization_Ngrains(i)
enddo enddo
write (6,*) write (6,*)
write (6,*) 'Microstructure' write (6,'(a32,x,a12)') 'microstructure ','constituents'
do i = 1,material_Nmicrostructure do i = 1,material_Nmicrostructure
write (6,'(a32,x,i4)') microstructure_name(i),microstructure_Nconstituents(i) write (6,'(a32,x,i4)') microstructure_name(i),microstructure_Nconstituents(i)
if (microstructure_Nconstituents(i) > 0_pInt) then
do j = 1,microstructure_Nconstituents(i)
write (6,'(a1,x,a32,x,a32,x,f6.4)') '>',phase_name(microstructure_phase(j,i)),&
texture_name(microstructure_texture(j,i)),&
microstructure_fraction(j,i)
enddo
write (6,*)
endif
enddo enddo
call material_populateGrains() call material_populateGrains()
@ -115,6 +125,7 @@ subroutine material_parseHomogenization(file,myPart)
use prec, only: pInt use prec, only: pInt
use IO use IO
use mesh, only: mesh_element
implicit none implicit none
character(len=*), intent(in) :: myPart character(len=*), intent(in) :: myPart
@ -128,16 +139,15 @@ subroutine material_parseHomogenization(file,myPart)
Nsections = IO_countSections(file,myPart) Nsections = IO_countSections(file,myPart)
material_Nhomogenization = Nsections material_Nhomogenization = Nsections
write (6,*) 'homogenization sections found',material_Nhomogenization
allocate(homogenization_name(Nsections)); homogenization_name = '' allocate(homogenization_name(Nsections)); homogenization_name = ''
allocate(homogenization_type(Nsections)); homogenization_type = '' allocate(homogenization_type(Nsections)); homogenization_type = ''
allocate(homogenization_typeInstance(Nsections)); homogenization_typeInstance = 0_pInt allocate(homogenization_typeInstance(Nsections)); homogenization_typeInstance = 0_pInt
allocate(homogenization_Ngrains(Nsections)); homogenization_Ngrains = 0_pInt allocate(homogenization_Ngrains(Nsections)); homogenization_Ngrains = 0_pInt
allocate(homogenization_Noutput(Nsections)); homogenization_Noutput = 0_pInt allocate(homogenization_Noutput(Nsections)); homogenization_Noutput = 0_pInt
allocate(homogenization_active(Nsections)); homogenization_active = .false.
write(6,*) 'scanning for (output)',homogenization_Noutput forall (s = 1:Nsections) homogenization_active(s) = any(mesh_element(3,:) == s) ! current homogenization used in model?
homogenization_Noutput = IO_countTagInPart(file,myPart,'(output)',Nsections) homogenization_Noutput = IO_countTagInPart(file,myPart,'(output)',Nsections)
write(6,*) 'count of (output)',homogenization_Noutput
rewind(file) rewind(file)
line = '' line = ''
@ -171,7 +181,7 @@ subroutine material_parseHomogenization(file,myPart)
endif endif
enddo enddo
100 homogenization_maxNgrains = maxval(homogenization_Ngrains) 100 homogenization_maxNgrains = maxval(homogenization_Ngrains,homogenization_active)
return return
endsubroutine endsubroutine
@ -183,6 +193,7 @@ subroutine material_parseMicrostructure(file,myPart)
use prec, only: pInt use prec, only: pInt
use IO use IO
use mesh, only: mesh_element
implicit none implicit none
character(len=*), intent(in) :: myPart character(len=*), intent(in) :: myPart
@ -197,6 +208,9 @@ subroutine material_parseMicrostructure(file,myPart)
material_Nmicrostructure = Nsections material_Nmicrostructure = Nsections
allocate(microstructure_name(Nsections)); microstructure_name = '' allocate(microstructure_name(Nsections)); microstructure_name = ''
allocate(microstructure_Nconstituents(Nsections)) allocate(microstructure_Nconstituents(Nsections))
allocate(microstructure_active(Nsections))
forall (i = 1:Nsections) microstructure_active(i) = any(mesh_element(4,:) == i) ! current microstructure used in model?
microstructure_Nconstituents = IO_countTagInPart(file,myPart,'(constituent)',Nsections) microstructure_Nconstituents = IO_countTagInPart(file,myPart,'(constituent)',Nsections)
microstructure_maxNconstituents = maxval(microstructure_Nconstituents) microstructure_maxNconstituents = maxval(microstructure_Nconstituents)
@ -452,7 +466,7 @@ subroutine material_populateGrains()
allocate(Ngrains(material_Nhomogenization,material_Nmicrostructure)); Ngrains = 0_pInt allocate(Ngrains(material_Nhomogenization,material_Nmicrostructure)); Ngrains = 0_pInt
! count grains per homog/micro pair ! identify maximum grain count per IP (from element) and find grains per homog/micro pair
do e = 1,mesh_NcpElems do e = 1,mesh_NcpElems
homog = mesh_element(3,e) homog = mesh_element(3,e)
micro = mesh_element(4,e) micro = mesh_element(4,e)
@ -463,7 +477,7 @@ subroutine material_populateGrains()
Ngrains(homog,micro) = Ngrains(homog,micro) + homogenization_Ngrains(homog) * FE_Nips(mesh_element(2,e)) Ngrains(homog,micro) = Ngrains(homog,micro) + homogenization_Ngrains(homog) * FE_Nips(mesh_element(2,e))
enddo enddo
allocate(volumeOfGrain(maxval(Ngrains))) ! reserve memory for maximum case allocate(volumeOfGrain(maxval(Ngrains))) ! reserve memory for maximum case
allocate(phaseOfGrain(maxval(Ngrains))) ! reserve memory for maximum case allocate(phaseOfGrain(maxval(Ngrains))) ! reserve memory for maximum case
allocate(orientationOfGrain(3,maxval(Ngrains))) ! reserve memory for maximum case allocate(orientationOfGrain(3,maxval(Ngrains))) ! reserve memory for maximum case

2
code/mpie_cpfem_marc.f90
View File

@ -43,7 +43,7 @@
include "mesh.f90" ! uses prec, math, IO, FEsolving include "mesh.f90" ! uses prec, math, IO, FEsolving
include "material.f90" ! uses prec, math, IO, mesh include "material.f90" ! uses prec, math, IO, mesh
include "lattice.f90" ! uses prec, math, IO, material include "lattice.f90" ! uses prec, math, IO, material
include "constitutive_phenomenological.f90" ! uses prec, math, IO, lattice, material, debug include "constitutive_phenopowerlaw.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive_j2.f90" ! uses prec, math, IO, lattice, material, debug include "constitutive_j2.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive_dislobased.f90" ! uses prec, math, IO, lattice, material, debug include "constitutive_dislobased.f90" ! uses prec, math, IO, lattice, material, debug
include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug include "constitutive.f90" ! uses prec, IO, math, lattice, mesh, debug

19
code/numerics.f90
View File

@ -57,25 +57,6 @@ subroutine numerics_init()
character(len=64) tag character(len=64) tag
character(len=1024) line character(len=1024) line
!*** global variables ***!
! relevantStrain
! iJacoStiffness
! iJacoLpresiduum
! pert_Fg
! nHomog
! nCryst
! nState
! nStress
! subStepMin
! rTol_crystalliteState
! rTol_crystalliteTemperature
! rTol_crystalliteStress
! aTol_crystalliteStress
! resToler
! resAbsol
! resBound
! NRiterMax
write(6,*) write(6,*)
write(6,*) '<<<+- numerics init -+>>>' write(6,*) '<<<+- numerics init -+>>>'
write(6,*) write(6,*)