major edit.
combined routines into CPFEM_general, CPFEM_stressIP (incl cutback scheme), CPFEM_stressCrystallite, and CPFEM_timeIntegration (Newton scheme) error management now based on text strings
This commit is contained in:
Philip Eisenlohr
parent
9da4b4cf24
commit
ff6642ea8c
810
trunk/CPFEM.f90
810
trunk/CPFEM.f90
|
|
@ -1,7 +1,7 @@
|
||||||
! last modified 29.03.07
|
|
||||||
! ---------------------------
|
!##############################################################
|
||||||
MODULE CPFEM
|
MODULE CPFEM
|
||||||
! ---------------------------
|
!##############################################################
|
||||||
! *** CPFEM engine ***
|
! *** CPFEM engine ***
|
||||||
!
|
!
|
||||||
use prec, only: pReal,pInt
|
use prec, only: pReal,pInt
|
||||||
|
|
@ -28,62 +28,11 @@
|
||||||
|
|
||||||
CONTAINS
|
CONTAINS
|
||||||
|
|
||||||
!***********************************************************************
|
!*********************************************************
|
||||||
!*** This routine checks for initialization, variables update and ***
|
!*** allocate the arrays defined in module CPFEM ***
|
||||||
!*** calls the actual material model ***
|
!*** and initialize them ***
|
||||||
!***********************************************************************
|
!*********************************************************
|
||||||
subroutine cpfem_general(ffn, ffn1, CPFEM_inc, CPFEM_subinc, CPFEM_cn, CPFEM_dt, cp_en, CPFEM_in)
|
SUBROUTINE CPFEM_init()
|
||||||
!
|
|
||||||
use prec, only: pReal,pInt
|
|
||||||
! use CPFEM, only: CPFEM_ffn_all, CPFEM_ffn1_all, CPFEM_inc_old
|
|
||||||
! use IO, only: IO_init
|
|
||||||
use constitutive, only: constitutive_state_old, constitutive_state_new
|
|
||||||
implicit none
|
|
||||||
!
|
|
||||||
real(pReal) ffn(3,3), ffn1(3,3), CPFEM_dt
|
|
||||||
integer(pInt) CPFEM_inc, CPFEM_subinc, CPFEM_cn, cp_en, CPFEM_in
|
|
||||||
!
|
|
||||||
! initialization step
|
|
||||||
if (CPFEM_first_call) then
|
|
||||||
! three dimensional stress state ?
|
|
||||||
! call IO_init()
|
|
||||||
call mesh_init()
|
|
||||||
call constitutive_init()
|
|
||||||
call math_init()
|
|
||||||
call CPFEM_init()
|
|
||||||
CPFEM_first_call = .false.
|
|
||||||
endif
|
|
||||||
! not a new increment
|
|
||||||
if (CPFEM_inc==CPFEM_inc_old) then
|
|
||||||
! case of a new subincrement:update starting with subinc 2
|
|
||||||
if (CPFEM_subinc > CPFEM_subinc_old) then
|
|
||||||
CPFEM_sigma_old = CPFEM_sigma_new
|
|
||||||
CPFEM_Fp_old = CPFEM_Fp_new
|
|
||||||
constitutive_state_old = constitutive_state_new
|
|
||||||
CPFEM_subinc_old = CPFEM_subinc
|
|
||||||
endif
|
|
||||||
! case of a new increment
|
|
||||||
else
|
|
||||||
CPFEM_sigma_old = CPFEM_sigma_new
|
|
||||||
CPFEM_Fp_old = CPFEM_Fp_new
|
|
||||||
constitutive_state_old = constitutive_state_new
|
|
||||||
CPFEM_inc_old = CPFEM_inc
|
|
||||||
CPFEM_subinc_old = 1_pInt
|
|
||||||
endif
|
|
||||||
!
|
|
||||||
! get cp element number for fe element number
|
|
||||||
CPFEM_ffn_all(:,:,CPFEM_in, cp_en) = ffn
|
|
||||||
CPFEM_ffn1_all(:,:,CPFEM_in, cp_en) = ffn1
|
|
||||||
call CPFEM_general_material(CPFEM_cn, CPFEM_dt, cp_en, CPFEM_in)
|
|
||||||
return
|
|
||||||
end subroutine
|
|
||||||
|
|
||||||
|
|
||||||
!***********************************************************************
|
|
||||||
!*** This routine allocates the arrays defined in module CPFEM ***
|
|
||||||
!*** and initializes them ***
|
|
||||||
!***********************************************************************
|
|
||||||
subroutine CPFEM_init()
|
|
||||||
!
|
!
|
||||||
use prec, only: pReal,pInt
|
use prec, only: pReal,pInt
|
||||||
! use math, only: math_I3
|
! use math, only: math_I3
|
||||||
|
|
@ -93,34 +42,25 @@
|
||||||
implicit none
|
implicit none
|
||||||
!
|
!
|
||||||
! *** mpie.marc parameters ***
|
! *** mpie.marc parameters ***
|
||||||
allocate(CPFEM_ffn_all(3,3,mesh_maxNips,mesh_NcpElems))
|
allocate(CPFEM_ffn_all (3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_ffn_all = 0.0_pReal
|
||||||
allocate(CPFEM_ffn1_all(3,3,mesh_maxNips,mesh_NcpElems))
|
allocate(CPFEM_ffn1_all (3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_ffn1_all = 0.0_pReal
|
||||||
allocate(CPFEM_stress_all(6,mesh_maxNips,mesh_NcpElems))
|
allocate(CPFEM_stress_all( 6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_stress_all = 0.0_pReal
|
||||||
allocate(CPFEM_jacobi_all(6,6,mesh_maxNips,mesh_NcpElems))
|
allocate(CPFEM_jacobi_all(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_jacobi_all = 0.0_pReal
|
||||||
CPFEM_ffn_all = 0.0_pReal
|
|
||||||
CPFEM_ffn1_all = 0.0_pReal
|
|
||||||
CPFEM_stress_all = 0.0_pReal
|
|
||||||
CPFEM_jacobi_all = 0.0_pReal
|
|
||||||
!
|
!
|
||||||
! *** User defined results !!! MISSING incorporate consti_Nresults ***
|
! *** User defined results !!! MISSING incorporate consti_Nresults ***
|
||||||
allocate(CPFEM_results(CPFEM_Nresults+constitutive_maxNresults,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems))
|
allocate(CPFEM_results(CPFEM_Nresults+constitutive_maxNresults,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems))
|
||||||
CPFEM_results = 0.0_pReal
|
CPFEM_results = 0.0_pReal
|
||||||
!
|
!
|
||||||
! *** Second Piola-Kirchoff stress tensor at (t=t0) and (t=t1) ***
|
! *** Second Piola-Kirchoff stress tensor at (t=t0) and (t=t1) ***
|
||||||
allocate(CPFEM_sigma_old(6,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems))
|
allocate(CPFEM_sigma_old(6,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_sigma_old = 0.0_pReal
|
||||||
allocate(CPFEM_sigma_new(6,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems))
|
allocate(CPFEM_sigma_new(6,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_sigma_new = 0.0_pReal
|
||||||
CPFEM_sigma_old = 0.0_pReal
|
|
||||||
CPFEM_sigma_new = 0.0_pReal
|
|
||||||
!
|
!
|
||||||
! *** Plastic deformation gradient at (t=t0) and (t=t1) ***
|
! *** Plastic deformation gradient at (t=t0) and (t=t1) ***
|
||||||
allocate(CPFEM_Fp_old(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems))
|
allocate(CPFEM_Fp_old(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Fp_old = 0.0_pReal
|
||||||
allocate(CPFEM_Fp_new(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems))
|
allocate(CPFEM_Fp_new(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Fp_new = 0.0_pReal
|
||||||
CPFEM_Fp_old = 0.0_pReal
|
|
||||||
CPFEM_Fp_new = 0.0_pReal
|
|
||||||
!
|
!
|
||||||
! *** Old jacobian (consistent tangent) ***
|
! *** Old jacobian (consistent tangent) ***
|
||||||
allocate(CPFEM_jaco_old(6,6,mesh_maxNips,mesh_NcpElems))
|
allocate(CPFEM_jaco_old(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_jaco_old = 0.0_pReal
|
||||||
CPFEM_jaco_old = 0.0_pReal
|
|
||||||
!
|
!
|
||||||
! *** Output to MARC output file ***
|
! *** Output to MARC output file ***
|
||||||
write(6,*)
|
write(6,*)
|
||||||
|
|
@ -138,468 +78,372 @@
|
||||||
write(6,*)
|
write(6,*)
|
||||||
call flush(6)
|
call flush(6)
|
||||||
return
|
return
|
||||||
end subroutine
|
|
||||||
|
END SUBROUTINE
|
||||||
!
|
!
|
||||||
!
|
!
|
||||||
subroutine CPFEM_general_material(&
|
!***********************************************************************
|
||||||
|
!*** perform initialization at first call, update variables and ***
|
||||||
|
!*** call the actual material model ***
|
||||||
|
!***********************************************************************
|
||||||
|
SUBROUTINE CPFEM_general(ffn, ffn1, CPFEM_inc, CPFEM_subinc, CPFEM_cn, CPFEM_dt, cp_en, CPFEM_in)
|
||||||
|
!
|
||||||
|
use prec, only: pReal,pInt
|
||||||
|
use constitutive, only: constitutive_state_old, constitutive_state_new
|
||||||
|
implicit none
|
||||||
|
!
|
||||||
|
real(pReal) ffn(3,3), ffn1(3,3), CPFEM_dt
|
||||||
|
integer(pInt) CPFEM_inc, CPFEM_subinc, CPFEM_cn, cp_en, CPFEM_in
|
||||||
|
!
|
||||||
|
! initialization step
|
||||||
|
if (CPFEM_first_call) then
|
||||||
|
! three dimensional stress state ?
|
||||||
|
call math_init()
|
||||||
|
call mesh_init()
|
||||||
|
call constitutive_init()
|
||||||
|
call CPFEM_init()
|
||||||
|
CPFEM_first_call = .false.
|
||||||
|
endif
|
||||||
|
if (CPFEM_inc==CPFEM_inc_old) then ! not a new increment
|
||||||
|
! case of a new subincrement:update starting with subinc 2
|
||||||
|
if (CPFEM_subinc > CPFEM_subinc_old) then
|
||||||
|
CPFEM_sigma_old = CPFEM_sigma_new
|
||||||
|
CPFEM_Fp_old = CPFEM_Fp_new
|
||||||
|
constitutive_state_old = constitutive_state_new
|
||||||
|
CPFEM_subinc_old = CPFEM_subinc
|
||||||
|
endif
|
||||||
|
else ! new increment
|
||||||
|
CPFEM_sigma_old = CPFEM_sigma_new
|
||||||
|
CPFEM_Fp_old = CPFEM_Fp_new
|
||||||
|
constitutive_state_old = constitutive_state_new
|
||||||
|
CPFEM_inc_old = CPFEM_inc
|
||||||
|
CPFEM_subinc_old = 1_pInt
|
||||||
|
endif
|
||||||
|
!
|
||||||
|
CPFEM_ffn_all(:,:,CPFEM_in, cp_en) = ffn
|
||||||
|
CPFEM_ffn1_all(:,:,CPFEM_in, cp_en) = ffn1
|
||||||
|
call CPFEM_stressIP(CPFEM_cn, CPFEM_dt, cp_en, CPFEM_in)
|
||||||
|
return
|
||||||
|
|
||||||
|
END SUBROUTINE
|
||||||
|
|
||||||
|
|
||||||
|
!**********************************************************
|
||||||
|
!*** calculate the material behaviour at IP level ***
|
||||||
|
!**********************************************************
|
||||||
|
SUBROUTINE CPFEM_stressIP(&
|
||||||
CPFEM_cn,& ! Cycle number
|
CPFEM_cn,& ! Cycle number
|
||||||
CPFEM_dt,& ! Time increment (dt)
|
CPFEM_dt,& ! Time increment (dt)
|
||||||
cp_en,& ! Element number
|
cp_en,& ! Element number
|
||||||
CPFEM_in) ! Integration point number
|
CPFEM_in) ! Integration point number
|
||||||
!***********************************************************************
|
|
||||||
!*** This routine calculates the material behaviour ***
|
use prec, only: pReal,pInt,ijaco,nCutback
|
||||||
!***********************************************************************
|
use IO, only: IO_error
|
||||||
use prec, only: pReal,pInt, ijaco
|
use mesh, only: mesh_element
|
||||||
! use IO, only: IO_error
|
|
||||||
use math
|
|
||||||
use mesh
|
|
||||||
use constitutive
|
use constitutive
|
||||||
!
|
!
|
||||||
implicit none
|
implicit none
|
||||||
!
|
|
||||||
! *** Definition of variables ***
|
|
||||||
! *** Subroutine parameters ***
|
|
||||||
real(pReal) CPFEM_dt
|
|
||||||
integer(pInt) CPFEM_cn, cp_en ,CPFEM_in
|
|
||||||
! *** Local variables ***
|
|
||||||
real(pReal) vf, cs(6), cd(6,6), CPFEM_d(6,6), CPFEM_s(6)
|
|
||||||
integer(pInt) jpara,nori, iori, ising, icut, iconv, CPFEM_en
|
|
||||||
!
|
|
||||||
! *** Flag for recalculation of jacobian ***
|
|
||||||
jpara = 1_pInt
|
|
||||||
! get number of grains from cp element number and integration point number
|
|
||||||
nori = constitutive_Ngrains(CPFEM_in,cp_en) !ÄÄÄ
|
|
||||||
!
|
|
||||||
CPFEM_en = mesh_element(1,cp_en) ! remap back to FE id
|
|
||||||
!
|
|
||||||
CPFEM_s=0
|
|
||||||
CPFEM_d=0
|
|
||||||
!
|
|
||||||
! *** Loop over all the components ***
|
|
||||||
do iori=1,nori
|
|
||||||
!
|
|
||||||
! *** Initialization of the matrices for t=t0 ***
|
|
||||||
! data from constitutive?
|
|
||||||
vf = constitutive_matVolFrac(iori,CPFEM_in,cp_en)*constitutive_texVolFrac(iori,CPFEM_in,cp_en) !ÄÄÄ
|
|
||||||
|
|
||||||
! *** Calculation of the solution at t=t1 ***
|
integer(pInt), parameter :: i_now = 1_pInt,i_then = 2_pInt
|
||||||
! QUESTION use the mod() as flag parameter in the call ??
|
character(len=128) msg
|
||||||
if (mod(CPFEM_cn,ijaco)==0) then !ÄÄÄ
|
integer(pInt) CPFEM_cn,cp_en,CPFEM_in,grain,i
|
||||||
call CPFEM_stress(cs, cd, CPFEM_dt,cp_en,CPFEM_in, iori, ising, icut, iconv, 1_pInt)
|
logical updateJaco
|
||||||
! *** Evaluation of ising ***
|
real(pReal) CPFEM_dt,dt,t,volfrac
|
||||||
! *** ising=2 => singular matrix in jacobi calculation ***
|
real(pReal), dimension(6) :: cs,Tstar_v
|
||||||
! *** => use old jacobi ***
|
real(pReal), dimension(6,6) :: cd,cd_IP
|
||||||
if (ising==2) jpara=0
|
real(pReal), dimension(3,3) :: deltaFg
|
||||||
! *** Calculation of the consistent tangent ***
|
real(pReal), dimension(3) :: Euler
|
||||||
CPFEM_d=CPFEM_d+vf*cd
|
real(pReal), dimension(3,3,2) :: Fg,Fp
|
||||||
else
|
real(pReal), dimension(constitutive_Nstatevars(grain,CPFEM_in,cp_en),2) :: state
|
||||||
call CPFEM_stress(cs, cd, CPFEM_dt,cp_en,CPFEM_in, iori, ising, icut, iconv, 0_pInt)
|
|
||||||
jpara=0
|
|
||||||
endif
|
|
||||||
! *** Cases of unsuccessful calculations ***
|
|
||||||
! *** Evaluation of ising ***
|
|
||||||
! *** ising!=0 => singular matrix ***
|
|
||||||
if (ising==1) then
|
|
||||||
write(6,*) 'Singular matrix!'
|
|
||||||
write(6,*) 'Integration point: ',CPFEM_in
|
|
||||||
write(6,*) 'Element: ',CPFEM_en
|
|
||||||
call IO_error(700)
|
|
||||||
! CPFEM_timefactor=1.e5_pReal
|
|
||||||
return
|
|
||||||
endif
|
|
||||||
! *** Evaluation of icut ***
|
|
||||||
! *** icut!=0 => too many cutbacks ***
|
|
||||||
if (icut==1) then
|
|
||||||
write(6,*) 'Too many cutbacks'
|
|
||||||
write(6,*) 'Integration point: ',CPFEM_in
|
|
||||||
write(6,*) 'Element: ',CPFEM_en
|
|
||||||
call IO_error(600)
|
|
||||||
! CPFEM_timefactor=1.e5_pReal
|
|
||||||
return
|
|
||||||
endif
|
|
||||||
! *** Evaluation of iconv ***
|
|
||||||
! *** iconv!=0 => no convergence ***
|
|
||||||
if (iconv==1) then
|
|
||||||
write(6,*) 'Inner loop did not converge!'
|
|
||||||
write(6,*) 'Integration point: ',CPFEM_in
|
|
||||||
write(6,*) 'Element: ',CPFEM_en
|
|
||||||
call IO_error(600)
|
|
||||||
! CPFEM_timefactor=1.e5_pReal
|
|
||||||
return
|
|
||||||
else if (iconv==2) then
|
|
||||||
write(6,*) 'Outer loop did not converge!'
|
|
||||||
write(6,*) 'Integration point: ',CPFEM_in
|
|
||||||
write(6,*) 'Element: ',CPFEM_en
|
|
||||||
call IO_error(600)
|
|
||||||
! CPFEM_timefactor=1.e5_pReal
|
|
||||||
return
|
|
||||||
endif
|
|
||||||
! *** Evaluation of the average Cauchy stress ***
|
|
||||||
CPFEM_s=CPFEM_s+vf*cs
|
|
||||||
enddo
|
|
||||||
! *** End of the loop over the components ***
|
|
||||||
|
|
||||||
! *************************************
|
updateJaco = (mod(CPFEM_cn,ijaco)==0) ! update consistent tangent every ijaco'th iteration
|
||||||
! *** End of the CP-FEM Calculation ***
|
|
||||||
! *************************************
|
CPFEM_stress_all(:,CPFEM_in,cp_en) = 0.0_pReal ! average Cauchy stress
|
||||||
! *** Store the new stress ***
|
cd_IP = 0.0_pReal ! average consistent tangent
|
||||||
CPFEM_stress_all(:,CPFEM_in,cp_en)=CPFEM_s
|
|
||||||
! *** Store the new jacobian ***
|
! -------------- grain loop -----------------
|
||||||
if (jpara/=0) CPFEM_jaco_old(:,:,CPFEM_in,cp_en)=CPFEM_d
|
do grain = 1,constitutive_Ngrains(CPFEM_in,cp_en)
|
||||||
return
|
! -------------------------------------------
|
||||||
end subroutine
|
|
||||||
!
|
i = 0_pInt ! cutback counter
|
||||||
!
|
state(:,i_now) = constitutive_state_old(:,grain,CPFEM_in,cp_en)
|
||||||
subroutine CPFEM_stress(&
|
Fg(:,:,i_now) = CPFEM_ffn_all(:,:,CPFEM_in,cp_en)
|
||||||
cs,& ! stress vector
|
Fp(:,:,i_now) = CPFEM_Fp_old(:,:,grain,CPFEM_in,cp_en)
|
||||||
cd,& ! Jacoby matrix
|
|
||||||
CPFEM_dt,& ! Time increment (dt)
|
deltaFg = CPFEM_ffn1_all(:,:,CPFEM_in,cp_en)-CPFEM_ffn_all(:,:,CPFEM_in,cp_en)
|
||||||
cp_en,& ! Element number
|
dt = CPFEM_dt
|
||||||
CPFEM_in,& ! Integration point number
|
|
||||||
iori,& ! number of orintation
|
Tstar_v = 0.0_pReal ! fully elastic initial guess
|
||||||
ising,& ! flag for singular matrix
|
Fg(:,:,i_then) = Fg(:,:,i_now)
|
||||||
icut,& ! flag for too many cut backs
|
state(:,i_then) = 0.0_pReal ! state_old as initial guess
|
||||||
iconv,& ! flag for non convergence
|
t = 0.0_pReal
|
||||||
isjaco) ! flag whether to calculate Jacoby matrix
|
|
||||||
!********************************************************************
|
! ------- crystallite integration -----------
|
||||||
! This routine calculates the stress for a single component
|
do
|
||||||
! and manages the independent time incrmentation
|
! -------------------------------------------
|
||||||
!********************************************************************
|
if (t+dt < CPFEM_dt) then ! intermediate solution
|
||||||
use prec, only: pReal,pInt, ncut
|
t = t+dt ! next time inc
|
||||||
use constitutive, only: constitutive_Nstatevars, constitutive_state_old, constitutive_state_new, constitutive_Nresults,&
|
Fg(:,:,i_then) = Fg(:,:,i_then)+deltaFg ! corresponding Fg
|
||||||
constitutive_results
|
else ! full step solution
|
||||||
implicit none
|
t = CPFEM_dt ! final time
|
||||||
!
|
Fg(:,:,i_then) = CPFEM_ffn_all(:,:,CPFEM_in,cp_en) ! final Fg
|
||||||
! *** Definition of variables ***
|
|
||||||
! *** Subroutine parameters ***
|
|
||||||
real(pReal) cs(6), cd(6,6), CPFEM_dt
|
|
||||||
integer(pInt) cp_en ,CPFEM_in, iori, ising, icut, iconv, isjaco
|
|
||||||
! *** Local variables ***
|
|
||||||
real(pReal) Fp_old(3,3), Fp_new(3,3), state_old(constitutive_Nstatevars(iori, CPFEM_in, cp_en))
|
|
||||||
real(pReal) state_new(constitutive_Nstatevars(iori, CPFEM_in, cp_en)), Tstar_v(6), CPFEM_ffn(3,3), CPFEM_ffn1(3,3)
|
|
||||||
real(pReal) Tstar_v_h(6), state_new_h(constitutive_Nstatevars(iori, CPFEM_in, cp_en)), phi1, PHI, phi2, dt_i
|
|
||||||
real(pReal) delta_Fg(3,3), Fg_i(3,3), state_new_i(constitutive_Nstatevars(iori, CPFEM_in, cp_en)), time
|
|
||||||
integer(pInt) jcut
|
|
||||||
!
|
|
||||||
icut=0
|
|
||||||
!
|
|
||||||
! *** Initialization of the matrices for t=t0 ***
|
|
||||||
Fp_old = CPFEM_Fp_old(:,:,iori,CPFEM_in,cp_en)
|
|
||||||
Fp_new = 0.0_pReal
|
|
||||||
state_old = constitutive_state_old(:,iori,CPFEM_in,cp_en)
|
|
||||||
state_new = state_old
|
|
||||||
Tstar_v = CPFEM_sigma_old(:,iori,CPFEM_in,cp_en)
|
|
||||||
CPFEM_ffn = CPFEM_ffn_all(:,:,CPFEM_in,cp_en)
|
|
||||||
CPFEM_ffn1 = CPFEM_ffn1_all(:,:,CPFEM_in,cp_en)
|
|
||||||
!
|
|
||||||
! *** First attempt to calculate Tstar and tauc with initial timestep ***
|
|
||||||
! save copies of Tstar_v and state_new
|
|
||||||
Tstar_v_h = Tstar_v
|
|
||||||
state_new_h = state_new
|
|
||||||
call CPFEM_stress_int(cs, cd, CPFEM_dt, cp_en,CPFEM_in, iori, ising, iconv, isjaco, phi1, PHI, phi2,&
|
|
||||||
CPFEM_ffn, CPFEM_ffn1,Fp_old,Fp_new,state_old, state_new, Tstar_v)
|
|
||||||
if ((iconv==0).AND.(ising==0)) then
|
|
||||||
! *** Update the differents matrices for t=t1 ***
|
|
||||||
CPFEM_Fp_new(:,:,iori,CPFEM_in,cp_en) = Fp_new
|
|
||||||
constitutive_state_new(:,iori,CPFEM_in,cp_en) = state_new
|
|
||||||
CPFEM_sigma_new(:,iori,CPFEM_in,cp_en) = Tstar_v
|
|
||||||
! *** Update the results plotted in MENTAT ***
|
|
||||||
CPFEM_results(1,iori,CPFEM_in,cp_en) = phi1
|
|
||||||
CPFEM_results(2,iori,CPFEM_in,cp_en) = PHI
|
|
||||||
CPFEM_results(3,iori,CPFEM_in,cp_en) = phi2
|
|
||||||
CPFEM_results(4:3+constitutive_Nresults(iori,CPFEM_in,cp_en),iori,CPFEM_in,cp_en)=&
|
|
||||||
constitutive_results(1:constitutive_Nresults(iori,CPFEM_in,cp_en),iori,CPFEM_in,cp_en)!ÄÄÄÄ
|
|
||||||
return
|
|
||||||
endif
|
endif
|
||||||
!
|
|
||||||
! *** Calculation of stress and resistences with a cut timestep ***
|
call CPFEM_stressCrystallite(msg,cs,cd,Tstar_v,Fp(:,:,i_then),state(:,i_then),Euler,&
|
||||||
! *** when first try did not converge ***
|
dt,cp_en,CPFEM_in,grain,updateJaco .and. t==CPFEM_dt,&
|
||||||
jcut=1_pInt
|
Fg(:,:,i_now),Fg(:,:,i_then),Fp(:,:,i_now),state(:,i_now))
|
||||||
dt_i=0.5_pReal*CPFEM_dt
|
if (msg == 'ok') then ! solution converged
|
||||||
delta_Fg=0.5_pReal*(CPFEM_ffn1-CPFEM_ffn)
|
if (t == CPFEM_dt) exit ! reached final "then"
|
||||||
Fg_i=CPFEM_ffn+delta_Fg
|
else ! solution not found
|
||||||
Tstar_v=Tstar_v_h
|
i = i+1_pInt ! inc cutback counter
|
||||||
state_new_i=state_new_h
|
if (i > nCutback) then ! limit exceeded?
|
||||||
! *** Start time ***
|
write(6,*) 'cutback limit --> '//msg
|
||||||
time=dt_i
|
write(6,*) 'Grain: ',grain
|
||||||
do while (time<=CPFEM_dt)
|
write(6,*) 'Integration point: ',CPFEM_in
|
||||||
call CPFEM_stress_int(cs, cd, time, cp_en,CPFEM_in, iori, ising, iconv, isjaco, phi1, PHI, phi2,&
|
write(6,*) 'Element: ',mesh_element(1,cp_en)
|
||||||
CPFEM_ffn, Fg_i,Fp_old,Fp_new,state_old, state_new_i, Tstar_v)
|
call IO_error(600)
|
||||||
if ((iconv==0).AND.(ising==0)) then
|
return ! byebye
|
||||||
time=time+dt_i
|
|
||||||
Fg_i=Fg_i+delta_Fg
|
|
||||||
Tstar_v_h=Tstar_v
|
|
||||||
state_new_h=state_new_i
|
|
||||||
else
|
else
|
||||||
jcut=jcut+1_pInt
|
t = t-dt ! rewind time
|
||||||
if (jcut>ncut) then
|
Fg(:,:,i_then) = Fg(:,:,i_then)-deltaFg ! rewind Fg
|
||||||
icut=1_pInt
|
dt = 0.5_pReal*dt ! cut time-step in half
|
||||||
return
|
deltaFg = 0.5_pReal*deltaFg ! cut Fg-step in half
|
||||||
endif
|
endif
|
||||||
dt_i=0.5_pReal*dt_i
|
|
||||||
time=time-dt_i
|
|
||||||
delta_Fg=0.5_pReal*delta_Fg
|
|
||||||
Fg_i=Fg_i-delta_Fg
|
|
||||||
Tstar_v=Tstar_v_h
|
|
||||||
state_new_i=state_new_h
|
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo ! crystallite integration (cutback loop)
|
||||||
!
|
|
||||||
! *** Final calculation of stress and resistences with full timestep ***
|
! ---- update crystallite matrices at t = t1 ----
|
||||||
state_new=state_new_i
|
CPFEM_Fp_new(:,:,grain,CPFEM_in,cp_en) = Fp(:,:,i_then)
|
||||||
call CPFEM_stress_int(cs, cd, CPFEM_dt, cp_en,CPFEM_in, iori, ising, iconv, isjaco, phi1, PHI, phi2,&
|
constitutive_state_new(:,grain,CPFEM_in,cp_en) = state(:,i_then)
|
||||||
CPFEM_ffn, CPFEM_ffn1,Fp_old,Fp_new,state_old, state_new, Tstar_v)
|
CPFEM_sigma_new(:,grain,CPFEM_in,cp_en) = Tstar_v
|
||||||
! *** Update the differents matrices for t=t1 ***
|
! ---- update results plotted in MENTAT ----
|
||||||
CPFEM_Fp_new(:,:,iori,CPFEM_in,cp_en) = Fp_new
|
CPFEM_results(1:3,grain,CPFEM_in,cp_en) = Euler
|
||||||
constitutive_state_new(:,iori,CPFEM_in,cp_en) = state_new
|
CPFEM_results(4:3+constitutive_Nresults(grain,CPFEM_in,cp_en),grain,CPFEM_in,cp_en) = &
|
||||||
CPFEM_sigma_new(:,iori,CPFEM_in,cp_en) = Tstar_v
|
constitutive_results(1:constitutive_Nresults(grain,CPFEM_in,cp_en),grain,CPFEM_in,cp_en)!ÄÄÄÄ
|
||||||
! *** Update the results plotted in MENTAT ***
|
|
||||||
CPFEM_results(1,iori,CPFEM_in,cp_en) = phi1
|
! ---- contribute to IP result ----
|
||||||
CPFEM_results(2,iori,CPFEM_in,cp_en) = PHI
|
volfrac = constitutive_matVolFrac(grain,CPFEM_in,cp_en)*constitutive_texVolFrac(grain,CPFEM_in,cp_en)
|
||||||
CPFEM_results(3,iori,CPFEM_in,cp_en) = phi2
|
CPFEM_stress_all(:,CPFEM_in,cp_en) = CPFEM_stress_all(:,CPFEM_in,cp_en)+volfrac*cs ! average Cauchy stress
|
||||||
|
if (updateJaco) cd_IP = cd_IP+volfrac*cd ! average consistent tangent
|
||||||
|
|
||||||
|
enddo ! grain loop
|
||||||
|
|
||||||
return
|
return
|
||||||
end subroutine
|
END SUBROUTINE
|
||||||
!
|
|
||||||
!
|
|
||||||
subroutine CPFEM_stress_int(&
|
|
||||||
cs,& ! Cauchy stress vector
|
|
||||||
dcs_de,& ! Consistent tangent
|
|
||||||
dt,& ! Time increment
|
|
||||||
cp_en,& ! Element number
|
|
||||||
CPFEM_in,& ! Integration point number
|
|
||||||
iori,& ! number of orintation
|
|
||||||
ising,& ! flag for singular matrix
|
|
||||||
iconv,& ! flag for non convergence
|
|
||||||
isjaco,& ! flag whether to calculate Jacoby matrix
|
|
||||||
phi1,& ! Euler angle
|
|
||||||
PHI,& ! Euler angle
|
|
||||||
phi2,& ! Euler angle
|
|
||||||
Fg_old,& ! Old global deformation gradient
|
|
||||||
Fg_new,& ! New global deformation gradient
|
|
||||||
Fp_old,& ! Old plastic deformation gradient
|
|
||||||
Fp_new,& ! New plastic deformation gradient
|
|
||||||
state_old,& ! Old state variable array
|
|
||||||
state_new,& ! New state variable array
|
|
||||||
Tstar_v) ! Second Piola-Kirschoff stress tensor
|
|
||||||
!********************************************************************
|
!********************************************************************
|
||||||
! This routine calculates the stress for a single component
|
! Calculates the stress for a single component
|
||||||
! it is based on the paper by Kalidindi et al.:
|
! it is based on the paper by Kalidindi et al.:
|
||||||
! J. Mech. Phys, Solids Vol. 40, No. 3, pp. 537-569, 1992
|
! J. Mech. Phys, Solids Vol. 40, No. 3, pp. 537-569, 1992
|
||||||
! it is modified to use anisotropic elasticity matrix
|
! it is modified to use anisotropic elasticity matrix
|
||||||
!********************************************************************
|
!********************************************************************
|
||||||
|
subroutine CPFEM_stressCrystallite(&
|
||||||
|
msg,& ! return message
|
||||||
|
cs,& ! Cauchy stress vector
|
||||||
|
dcs_de,& ! consistent tangent
|
||||||
|
Tstar_v,& ! second Piola-Kirchoff stress tensor
|
||||||
|
Fp_new,& ! new plastic deformation gradient
|
||||||
|
state_new,& ! new state variable array
|
||||||
|
Euler,& ! Euler angles
|
||||||
|
!
|
||||||
|
dt,& ! time increment
|
||||||
|
cp_en,& ! element number
|
||||||
|
CPFEM_in,& ! integration point number
|
||||||
|
grain,& ! grain number
|
||||||
|
updateJaco,& ! boolean to calculate Jacobi matrix
|
||||||
|
Fg_old,& ! old global deformation gradient
|
||||||
|
Fg_new,& ! new global deformation gradient
|
||||||
|
Fp_old,& ! old plastic deformation gradient
|
||||||
|
state_old) ! old state variable array
|
||||||
|
|
||||||
use prec, only: pReal,pInt,pert_e
|
use prec, only: pReal,pInt,pert_e
|
||||||
use constitutive, only: constitutive_Nstatevars
|
use constitutive, only: constitutive_Nstatevars
|
||||||
use math, only: math_Mandel6to33
|
use math, only: math_Mandel6to33, mapMandel,math_pDecomposition,math_RtoEuler
|
||||||
implicit none
|
implicit none
|
||||||
!
|
|
||||||
! *** Definition of variables ***
|
|
||||||
! *** Subroutine parameters ***
|
|
||||||
integer(pInt) cp_en, CPFEM_in, iori, ising, iconv, isjaco
|
|
||||||
real(pReal) cs(6), dcs_de(6,6), dt, phi1, PHI, phi2, Fg_old(3,3), Fg_new(3,3)
|
|
||||||
real(pReal) Fp_old(3,3), Fp_new(3,3), state_old(constitutive_Nstatevars(iori, CPFEM_in, cp_en))
|
|
||||||
real(pReal) state_new(constitutive_Nstatevars(iori, CPFEM_in, cp_en)), Tstar_v(6)
|
|
||||||
! *** Local variables ***
|
|
||||||
integer(pInt) ic
|
|
||||||
real(pReal) Fe(3,3), R(3,3), U(3,3), Fg_pert(3,3), sgm2(6)
|
|
||||||
real(pReal) state2(constitutive_Nstatevars(iori, CPFEM_in, cp_en)), Fp2(3,3), cs1(6),E_pert(3,3)
|
|
||||||
! *** Error treatment ***
|
|
||||||
iconv = 0
|
|
||||||
ising = 0
|
|
||||||
|
|
||||||
! *********************************************
|
character(len=*) msg
|
||||||
! *** Calculation of the new Cauchy stress ***
|
logical updateJaco,error
|
||||||
! *********************************************
|
integer(pInt) cp_en,CPFEM_in,grain,i
|
||||||
|
real(pReal) dt
|
||||||
|
real(pReal), dimension(3) :: Euler
|
||||||
|
real(pReal), dimension(3,3) :: Fg_old,Fg_new,Fg_pert,Fp_old,Fp_new,Fp_pert,Fe_new,Fe_pert,R,U,E_pert
|
||||||
|
real(pReal), dimension(6) :: cs,Tstar_v,Tstar_v_pert
|
||||||
|
real(pReal), dimension(6,6) :: dcs_de
|
||||||
|
real(pReal), dimension(constitutive_Nstatevars(grain, CPFEM_in, cp_en)) :: state_old,state_new,state_pert
|
||||||
|
|
||||||
! *** Call Newton-Raphson method ***
|
call CPFEM_timeIntegration(msg,Fp_new,Fe_new,Tstar_v,state_new, &
|
||||||
call NEWTON_RAPHSON(dt,cp_en,CPFEM_in,iori,Fg_new,Fp_old,Fp_new,Fe,state_old,state_new,Tstar_v,cs,iconv,ising)
|
dt,cp_en,CPFEM_in,grain,Fg_new,Fp_old,state_old)
|
||||||
!
|
if (msg /= 'ok') return
|
||||||
! *** Calculation of the new orientation ***
|
cs = CPFEM_CauchyStress(Tstar_v,Fe_new) ! Cauchy stress
|
||||||
call math_pDecomposition(Fe,U,R,ising)
|
|
||||||
if (ising==1) then
|
call math_pDecomposition(Fe_new,U,R,error) ! polar decomposition
|
||||||
|
if (error) then
|
||||||
|
msg = 'polar decomposition'
|
||||||
return
|
return
|
||||||
endif
|
endif
|
||||||
call math_RtoEuler(transpose(R),phi1,PHI,phi2)
|
Euler = math_RtoEuler(transpose(R)) ! orientation
|
||||||
!
|
|
||||||
! *** Choice of the calculation of the consistent tangent ***
|
if (updateJaco) then
|
||||||
if (isjaco==0) return
|
|
||||||
!
|
|
||||||
! *********************************************
|
|
||||||
! *** Calculation of the consistent tangent ***
|
|
||||||
! *********************************************
|
|
||||||
!
|
|
||||||
! *** Calculation of the consistent tangent with perturbation ***
|
! *** Calculation of the consistent tangent with perturbation ***
|
||||||
! *** Perturbation on the component of Fg ***
|
! *** Perturbation on the component of Fg ***
|
||||||
do ic=1,6
|
do i = 1,6
|
||||||
!
|
E_pert = 0.0_pReal
|
||||||
! *** Method of small perturbation
|
E_pert(mapMandel(1,i),mapMandel(2,i)) = E_pert(mapMandel(1,i),mapMandel(2,i)) + pert_e/2.0_pReal
|
||||||
! Missing direct matrix perturbation
|
E_pert(mapMandel(2,i),mapMandel(1,i)) = E_pert(mapMandel(2,i),mapMandel(1,i)) + pert_e/2.0_pReal
|
||||||
E_pert=0
|
|
||||||
if(ic<=3) then
|
|
||||||
E_pert(ic,ic) = pert_e
|
|
||||||
else if(ic==4) then
|
|
||||||
E_pert(1,2) = pert_e/2
|
|
||||||
E_pert(2,1) = pert_e/2
|
|
||||||
else if(ic==5) then
|
|
||||||
E_pert(2,3) = pert_e/2
|
|
||||||
E_pert(3,2) = pert_e/2
|
|
||||||
else if(ic==6) then
|
|
||||||
E_pert(1,3) = pert_e/2
|
|
||||||
E_pert(3,1) = pert_e/2
|
|
||||||
end if
|
|
||||||
Fg_pert=Fg_new+matmul(E_pert, Fg_old)
|
|
||||||
sgm2=Tstar_v
|
|
||||||
state2=state_new
|
|
||||||
|
|
||||||
! *** Calculation of the perturbated Cauchy stress ***
|
Fg_pert = Fg_new+matmul(E_pert,Fg_old) ! perturbated Fg
|
||||||
call NEWTON_RAPHSON(dt,cp_en,CPFEM_in,iori,Fg_pert,Fp_old,Fp2,Fe,state_old,state2,sgm2,cs1,iconv,ising)
|
Tstar_v_pert = Tstar_v ! initial guess at center
|
||||||
!
|
state_pert = state_new ! initial guess at center
|
||||||
! *** Consistent tangent *** as cs is Mandel dcs_de(:,4:6) is too large by sqrt(2)
|
|
||||||
dcs_de(:,ic)=(cs1-cs)/pert_e
|
call CPFEM_timeIntegration(msg,Fp_pert,Fe_pert,Tstar_v_pert,state_pert, &
|
||||||
enddo
|
dt,cp_en,CPFEM_in,grain,Fg_pert,Fp_old,state_old)
|
||||||
!
|
if (msg /= 'ok') then
|
||||||
|
msg = 'consistent tangent --> '//msg
|
||||||
return
|
return
|
||||||
end subroutine
|
endif
|
||||||
!
|
! *** MISSING:Consistent tangent, (perturbated) Cauchy stress is Mandel hence dcs_de(:,4:6) is too large by sqrt(2)
|
||||||
!
|
dcs_de(:,i) = (CPFEM_CauchyStress(Tstar_v_pert,Fe_pert)-cs)/pert_e
|
||||||
subroutine NEWTON_RAPHSON(&
|
enddo
|
||||||
dt,&
|
endif
|
||||||
cp_en,& ! Element number
|
return
|
||||||
CPFEM_in,& ! Integration point number
|
|
||||||
iori,& ! number of orientation
|
END SUBROUTINE
|
||||||
Fg_new,&
|
|
||||||
Fp_old,&
|
|
||||||
Fp_new,&
|
|
||||||
Fe,&
|
|
||||||
state_old,&
|
|
||||||
state_new,&
|
|
||||||
Tstar_v,&
|
|
||||||
cs,&
|
|
||||||
iconv,&
|
|
||||||
ising)
|
|
||||||
!***********************************************************************
|
!***********************************************************************
|
||||||
!*** NEWTON-RAPHSON Calculation ***
|
!*** fully-implicit two-level time integration ***
|
||||||
!***********************************************************************
|
!***********************************************************************
|
||||||
use prec, only: pReal,pInt, nouter, tol_outer, ninner, tol_inner, crite
|
SUBROUTINE CPFEM_timeIntegration(&
|
||||||
use constitutive, only: constitutive_Nstatevars, constitutive_HomogenizedC, constitutive_dotState
|
msg,& ! return message
|
||||||
|
Fp_new,& ! new plastic deformation gradient
|
||||||
|
Fe_new,& ! new "elastic" deformation gradient
|
||||||
|
Tstar_v,& ! 2nd PK stress (taken as initial guess if /= 0)
|
||||||
|
state_new,& ! current microstructure at end of time inc (taken as guess if /= 0)
|
||||||
|
!
|
||||||
|
dt,& ! time increment
|
||||||
|
cp_en,& ! element number
|
||||||
|
CPFEM_in,& ! integration point number
|
||||||
|
grain,& ! grain number
|
||||||
|
Fg_new,& ! new total def gradient
|
||||||
|
Fp_old,& ! former plastic def gradient
|
||||||
|
state_old) ! former microstructure
|
||||||
|
|
||||||
|
use prec, only: pReal,pInt, nState,tol_State,nStress,tol_Stress, crite, nReg
|
||||||
|
use constitutive, only: constitutive_Nstatevars,&
|
||||||
|
constitutive_homogenizedC,constitutive_dotState,constitutive_LpAndItsTangent
|
||||||
use math
|
use math
|
||||||
implicit none
|
implicit none
|
||||||
! *** Definition of variables ***
|
|
||||||
! *** Subroutine parameters ***
|
character(len=*) msg
|
||||||
integer(pInt) cp_en, CPFEM_in, iori, iconv, ising
|
integer(pInt) cp_en, CPFEM_in, grain
|
||||||
real(pReal) dt,Fg_new(3,3),Fp_old(3,3),Fp_new(3,3), Fe(3,3)
|
integer(pInt) iState,iStress,dummy, i,j,k,l,m
|
||||||
real(pReal) state_old(constitutive_Nstatevars(iori, CPFEM_in, cp_en)), state_new(constitutive_Nstatevars(iori, CPFEM_in, cp_en))
|
real(pReal) dt,det
|
||||||
real(pReal) Tstar_v(6), cs(6)
|
real(pReal), dimension(6) :: Tstar_v,dTstar_v,Rstress
|
||||||
! *** Local variables ***
|
real(pReal), dimension(6,6) :: C_66,Jacobi,invJacobi,help2
|
||||||
real(pReal) invFp_old(3,3), det, A(3,3), C_66(6,6), Lp(3,3), dLp(3,3,3,3)
|
real(pReal), dimension(3,3) :: Fg_new,Fp_old,Fp_new,Fe_new,invFp_old,invFp_new,Lp,A,B,AB
|
||||||
real(pReal) I3tLp(3,3), help(3,3), help1(3,3,3,3), Tstar0_v(6), R1(6)
|
real(pReal), dimension(3,3,3,3) :: dLp, LTL
|
||||||
real(pReal) dstate(constitutive_Nstatevars(iori, CPFEM_in, cp_en)), R2(constitutive_Nstatevars(iori, CPFEM_in, cp_en))
|
real(pReal), dimension(constitutive_Nstatevars(grain, CPFEM_in, cp_en)) :: state_old,state_new,dstate,Rstate,RstateS
|
||||||
real(pReal) R2s(constitutive_Nstatevars(iori, CPFEM_in, cp_en)), invFp_new(3,3)
|
logical failed
|
||||||
real(pReal) Jacobi(6,6), invJacobi(6,6), dTstar_v(6), help2(6,6)
|
|
||||||
integer(pInt) iouter, iinner , dummy, err, i, j, k, l, m
|
msg = 'ok' ! error-free so far
|
||||||
!
|
|
||||||
! *** Error treatment ***
|
call math_invert3x3(Fp_old,invFp_old,det,failed) ! inversion of Fp
|
||||||
iconv = 0
|
if (failed) then
|
||||||
ising = 0
|
msg = 'inversion Fp_old'
|
||||||
!
|
|
||||||
! initialize new state
|
|
||||||
state_new=state_old
|
|
||||||
! *** Calculation of Fp_old(-1) ***
|
|
||||||
call invert3x3(Fp_old, invFp_old, det, err) !ÄÄÄ
|
|
||||||
if (err==1_pInt) then
|
|
||||||
ising=1
|
|
||||||
return
|
return
|
||||||
endif
|
endif
|
||||||
!
|
|
||||||
! *** Calculation of A and T*0 (see Kalidindi) ***
|
C_66 = constitutive_HomogenizedC(grain, CPFEM_in, cp_en)
|
||||||
A = matmul(Fg_new,invFp_old) ! actually Fe
|
A = matmul(Fg_new,invFp_old) ! actually Fe
|
||||||
A = matmul(transpose(A), A)
|
A = matmul(transpose(A), A)
|
||||||
C_66 = constitutive_HomogenizedC(iori, CPFEM_in, cp_en) !ÄÄÄ
|
|
||||||
Tstar_v = 0.5_pReal*matmul(C_66, math_Mandel33to6(A-math_I3)) ! fully elastic guess ADDED 1/2
|
! former state guessed, if none specified
|
||||||
|
if (all(state_new == 0.0_pReal)) state_new = state_old
|
||||||
|
RstateS = state_new
|
||||||
|
iState = 0_pInt
|
||||||
|
! fully elastic guess (Lp = 0), if none specified
|
||||||
|
if (all(Tstar_v == 0.0_pReal)) Tstar_v = 0.5_pReal*matmul(C_66,math_Mandel33to6(A-math_I3))
|
||||||
! QUESTION follow former plastic slope to guess better?
|
! QUESTION follow former plastic slope to guess better?
|
||||||
!
|
Rstress = Tstar_v
|
||||||
! *** Second level of iterative procedure: Resistences ***
|
iStress = 0_pInt
|
||||||
do iouter=1,nouter
|
|
||||||
! *** First level of iterative procedure: Stresses ***
|
state: do ! outer iteration: state
|
||||||
do iinner=1,ninner
|
iState = iState+1
|
||||||
!
|
if (iState > nState) then
|
||||||
! *** Calculation of gdot_slip ***
|
msg = 'limit state iteration'
|
||||||
call constitutive_LpAndItsTangent(Tstar_v, iori, CPFEM_in, cp_en, Lp, dLp)
|
return
|
||||||
I3tLp = math_I3-dt*Lp
|
endif
|
||||||
help=matmul(transpose(I3tLp),matmul(A, I3tLp))
|
stress: do ! inner iteration: stress
|
||||||
Tstar0_v = 0.5_pReal * matmul(C_66, math_Mandel33to6(help-math_I3))
|
iStress = iStress+1
|
||||||
R1=Tstar_v-Tstar0_v
|
if (iStress > nStress) then ! too many loops required
|
||||||
if (maxval(abs(R1/maxval(abs(Tstar_v)))) < tol_inner) goto 100
|
msg = 'limit stress iteration'
|
||||||
!
|
return
|
||||||
! *** Jacobi Calculation: dRes/dTstar ***
|
endif
|
||||||
help=matmul(A, I3tLp)
|
call constitutive_LpAndItsTangent(Lp,dLp, Tstar_v,state_new,grain,CPFEM_in,cp_en)
|
||||||
help1=0.0_pReal
|
B = math_I3-dt*Lp
|
||||||
|
Rstress = Tstar_v - 0.5_pReal*matmul(C_66,math_Mandel33to6(matmul(transpose(B),matmul(A,B))-math_I3))
|
||||||
|
if (maxval(abs(Rstress/maxval(abs(Tstar_v)))) < tol_Stress) exit stress
|
||||||
|
|
||||||
|
! update stress guess using inverse of dRes/dTstar (Newton--Raphson)
|
||||||
|
AB = matmul(A,B)
|
||||||
|
LTL = 0.0_pReal
|
||||||
do i=1,3
|
do i=1,3
|
||||||
do j=1,3
|
do j=1,3
|
||||||
do k=1,3
|
do k=1,3
|
||||||
do l=1,3
|
do l=1,3
|
||||||
do m=1,3
|
do m=1,3
|
||||||
help1(i,j,k,l)=help1(i,j,k,l)+help(i,m)*dLp(m,j,k,l)+help(j,m)*dLp(m,i,l,k)
|
! LTL(i,j,k,l) = LTL(i,j,k,l) + AB(i,m)*dLp(m,j,k,l) + AB(j,m)*dLp(m,i,l,k) ! old
|
||||||
|
LTL(i,j,k,l) = LTL(i,j,k,l) + dLp(j,i,k,m)*AB(m,l) + AB(m,i)*dLp(m,j,k,l) ! new (and correct??)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
help2=math_Mandel3333to66(help1)
|
|
||||||
Jacobi= 0.5_pReal*matmul(C_66, help2) + math_identity2nd(6)
|
Jacobi = math_identity2nd(6) + 0.5_pReal*dt*matmul(C_66,math_Mandel3333to66(LTL))
|
||||||
call math_invert6x6(Jacobi, invJacobi, dummy, err) !ÄÄÄ
|
j = 0_pInt ; failed = .true.
|
||||||
if (err==1_pInt) then
|
do while (failed .and. j <= nReg)
|
||||||
forall (i=1:6) Jacobi(i,i)=1.05d0*maxval(Jacobi(i,:)) ! regularization
|
call math_invert6x6(Jacobi,invJacobi,dummy,failed)
|
||||||
call math_invert6x6(Jacobi, invJacobi, dummy, err)
|
forall (i=1:6) Jacobi(i,i) = 1.05_pReal*maxval(Jacobi(i,:)) ! regularization
|
||||||
if (err==1_pInt) then ! sorry, can't help here!!
|
j = j+1
|
||||||
ising=1
|
enddo
|
||||||
|
if (failed) then
|
||||||
|
msg = 'regularization Jacobi'
|
||||||
return
|
return
|
||||||
endif
|
endif
|
||||||
endif
|
|
||||||
dTstar_v=matmul(invJacobi,R1) ! correction to Tstar
|
|
||||||
|
|
||||||
! *** Correction (see Kalidindi) ***
|
dTstar_v = matmul(invJacobi,Rstress) ! correction to Tstar
|
||||||
forall(i=1:6, abs(dTstar_v(i)) > crite*maxval(abs(Tstar_v))) &
|
forall(i=1:6, abs(dTstar_v(i)) > crite*maxval(abs(Tstar_v))) &
|
||||||
dTstar_v(i) = sign(crite*maxval(abs(Tstar_v)),dTstar_v(i))
|
dTstar_v(i) = sign(crite*maxval(abs(Tstar_v)),dTstar_v(i)) ! cap to maximum correction
|
||||||
|
|
||||||
Tstar_v = Tstar_v-dTstar_v
|
Tstar_v = Tstar_v-dTstar_v
|
||||||
!
|
|
||||||
enddo
|
|
||||||
iconv=1
|
|
||||||
return
|
|
||||||
! *** End of the first level of iterative procedure ***
|
|
||||||
|
|
||||||
100 dstate=dt*constitutive_dotState(Tstar_v, iori, CPFEM_in, cp_en)
|
enddo stress
|
||||||
! *** Arrays of residuals ***
|
|
||||||
R2=state_new-state_old-dstate
|
dstate = dt*constitutive_dotState(Tstar_v,state_new,grain,CPFEM_in,cp_en) ! evolution of microstructure
|
||||||
R2s=0.0_pReal
|
Rstate = state_new - (state_old+dstate)
|
||||||
forall(i=1:constitutive_Nstatevars(iori, CPFEM_in, cp_en), state_new(i)/=0.0_pReal) R2s(i)=R2(i)/state_new(i)
|
RstateS = 0.0_pReal
|
||||||
if (maxval(abs(R2s)) < tol_outer) goto 200
|
forall (i=1:constitutive_Nstatevars(grain,CPFEM_in,cp_en), state_new(i)/=0.0_pReal) &
|
||||||
|
RstateS(i) = Rstress(i)/state_new(i)
|
||||||
|
if (maxval(abs(RstateS)) < tol_State) exit state
|
||||||
state_new = state_old+dstate
|
state_new = state_old+dstate
|
||||||
enddo
|
|
||||||
iconv=2
|
|
||||||
return
|
|
||||||
! *** End of the second level of iterative procedure ***
|
|
||||||
|
|
||||||
! *** Calculation of Fp(t+dt) (see Kalidindi) ***
|
enddo state
|
||||||
200 invFp_new=matmul(Fp_old, I3tLp)
|
|
||||||
call math_invert3x3(invFp_new, Fp_new, det, err) !ÄÄÄ
|
invFp_new = matmul(invFp_old,B)
|
||||||
if (err==1_pInt) then
|
call math_invert3x3(invFp_new,Fp_new,det,failed)
|
||||||
ising=1
|
if (failed) then
|
||||||
|
msg = 'inversion Fp_new'
|
||||||
return
|
return
|
||||||
endif
|
endif
|
||||||
Fp_new=Fp_new/math_det3x3(Fp_new)**(1.0_pReal/3.0_pReal)
|
Fp_new = Fp_new*det**(1.0_pReal/3.0_pReal) ! det = det(InvFp_new) !!
|
||||||
!
|
Fe_new = matmul(Fg_new,invFp_new)
|
||||||
! *** Calculation of F*(t+dt) (see Kalidindi) ***
|
|
||||||
Fe=matmul(Fg_new,invFp_new)
|
|
||||||
!
|
|
||||||
! *** Calculation of the Cauchy stress ***
|
|
||||||
! QUESTION seems to need Tstar, not Estar..??
|
|
||||||
cs = CPFEM_cauchy_stress(Tstar_v,Fe)
|
|
||||||
!
|
|
||||||
return
|
return
|
||||||
end subroutine
|
END SUBROUTINE
|
||||||
!
|
|
||||||
function CPFEM_cauchy_stress(PK_v, Fe)
|
|
||||||
|
FUNCTION CPFEM_CauchyStress(PK_v,Fe)
|
||||||
!***********************************************************************
|
!***********************************************************************
|
||||||
!*** Cauchy stress calculation ***
|
!*** Cauchy stress calculation ***
|
||||||
!***********************************************************************
|
!***********************************************************************
|
||||||
|
|
@ -607,8 +451,12 @@
|
||||||
use math, only: math_Mandel33to6,math_Mandel6to33,math_det3x3
|
use math, only: math_Mandel33to6,math_Mandel6to33,math_det3x3
|
||||||
implicit none
|
implicit none
|
||||||
! *** Subroutine parameters ***
|
! *** Subroutine parameters ***
|
||||||
real(pReal) PK_v(6), Fe(3,3), CPFEM_cauchy_stress(6)
|
real(pReal) PK_v(6), Fe(3,3), CPFEM_CauchyStress(6)
|
||||||
|
|
||||||
|
CPFEM_CauchyStress = math_Mandel33to6(matmul(matmul(Fe,math_Mandel6to33(PK_v)),transpose(Fe))/math_det3x3(Fe))
|
||||||
|
return
|
||||||
|
END FUNCTION
|
||||||
|
|
||||||
|
|
||||||
|
END MODULE
|
||||||
|
|
||||||
CPFEM_cauchy_stress = math_Mandel33to6(matmul(matmul(Fe,math_Mandel6to33(PK_v)),transpose(Fe))/math_det3x3(Fe))
|
|
||||||
end function
|
|
||||||
end module
|
|
||||||
Loading…
Reference in New Issue