DAMASK_EICMD/trunk/CPFEM.f90

614 lines
23 KiB
Fortran
Raw Normal View History

! last modified 29.03.07
! ---------------------------
MODULE CPFEM
! ---------------------------
! *** CPFEM engine ***
!
use prec, only: pReal,pInt
implicit none
!
! ****************************************************************
! *** General variables for the material behaviour calculation ***
! ****************************************************************
2007-03-26 15:57:34 +05:30
real(pReal), dimension (:,:,:), allocatable :: CPFEM_stress_all
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_jacobi_all
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_ffn_all
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_ffn1_all
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_results
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_ini_ori
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_sigma_old
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_sigma_new
real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_Fp_old
real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_Fp_new
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_jaco_old
integer(pInt) :: CPFEM_inc_old = 0_pInt
integer(pInt) :: CPFEM_subinc_old = 1_pInt
2007-03-26 15:57:34 +05:30
integer(pInt) :: CPFEM_Nresults = 3_pInt
2007-03-22 20:19:42 +05:30
logical :: CPFEM_first_call = .true.
CONTAINS
!***********************************************************************
!*** This routine checks for initialization, variables update and ***
!*** calls the actual material model ***
!***********************************************************************
2007-03-22 20:19:42 +05:30
subroutine cpfem_general(ffn, ffn1, CPFEM_inc, CPFEM_subinc, CPFEM_cn, CPFEM_dt, cp_en, CPFEM_in)
!
use prec, only: pReal,pInt
2007-03-26 15:57:34 +05:30
! 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
2007-03-22 20:19:42 +05:30
integer(pInt) CPFEM_inc, CPFEM_subinc, CPFEM_cn, cp_en, CPFEM_in
!
! initialization step
2007-03-22 20:19:42 +05:30
if (CPFEM_first_call) then
! three dimensional stress state ?
! call IO_init()
2007-03-22 20:19:42 +05:30
call mesh_init()
call constitutive_init()
call math_init()
call CPFEM_init()
2007-03-26 15:57:34 +05:30
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
2007-03-26 15:57:34 +05:30
end subroutine
!***********************************************************************
!*** This routine allocates the arrays defined in module CPFEM ***
!*** and initializes them ***
!***********************************************************************
2007-03-22 20:19:42 +05:30
subroutine CPFEM_init()
!
use prec, only: pReal,pInt
! use math, only: math_I3
use mesh
use constitutive
!
implicit none
!
! *** mpie.marc parameters ***
2007-03-22 20:19:42 +05:30
allocate(CPFEM_ffn_all(3,3,mesh_maxNips,mesh_NcpElems))
allocate(CPFEM_ffn1_all(3,3,mesh_maxNips,mesh_NcpElems))
allocate(CPFEM_stress_all(6,mesh_maxNips,mesh_NcpElems))
allocate(CPFEM_jacobi_all(6,6,mesh_maxNips,mesh_NcpElems))
CPFEM_ffn_all = 0.0_pReal
CPFEM_ffn1_all = 0.0_pReal
CPFEM_stress_all = 0.0_pReal
CPFEM_jacobi_all = 0.0_pReal
!
2007-03-22 20:19:42 +05:30
! *** User defined results !!! MISSING incorporate consti_Nresults ***
allocate(CPFEM_results(CPFEM_Nresults+constitutive_maxNresults,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems))
CPFEM_results = 0.0_pReal
!
! *** 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_new(6,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems))
CPFEM_sigma_old = 0.0_pReal
CPFEM_sigma_new = 0.0_pReal
!
! *** Plastic deformation gradient at (t=t0) and (t=t1) ***
allocate(CPFEM_Fp_old(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems))
allocate(CPFEM_Fp_new(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems))
CPFEM_Fp_old = 0.0_pReal
CPFEM_Fp_new = 0.0_pReal
!
! *** Old jacobian (consistent tangent) ***
2007-03-22 20:19:42 +05:30
allocate(CPFEM_jaco_old(6,6,mesh_maxNips,mesh_NcpElems))
CPFEM_jaco_old = 0.0_pReal
!
! *** Output to MARC output file ***
write(6,*)
write(6,*) 'Arrays allocated:'
write(6,*) 'CPFEM_ffn_all: ', shape(CPFEM_ffn_all)
write(6,*) 'CPFEM_ffn1_all: ', shape(CPFEM_ffn1_all)
write(6,*) 'CPFEM_stress_all: ', shape(CPFEM_stress_all)
write(6,*) 'CPFEM_jacobi_all: ', shape(CPFEM_jacobi_all)
write(6,*) 'CPFEM_results: ', shape(CPFEM_results)
write(6,*) 'CPFEM_sigma_old: ', shape(CPFEM_sigma_old)
write(6,*) 'CPFEM_sigma_new: ', shape(CPFEM_sigma_new)
write(6,*) 'CPFEM_Fp_old: ', shape(CPFEM_Fp_old)
write(6,*) 'CPFEM_Fp_new: ', shape(CPFEM_Fp_new)
write(6,*) 'CPFEM_jaco_old: ', shape(CPFEM_jaco_old)
write(6,*)
call flush(6)
return
2007-03-22 20:19:42 +05:30
end subroutine
!
!
subroutine CPFEM_general_material(&
CPFEM_cn,& ! Cycle number
CPFEM_dt,& ! Time increment (dt)
2007-03-22 20:19:42 +05:30
cp_en,& ! Element number
CPFEM_in) ! Integration point number
!***********************************************************************
!*** This routine calculates the material behaviour ***
!***********************************************************************
use prec, only: pReal,pInt, ijaco
2007-03-26 15:57:34 +05:30
! use IO, only: IO_error
use math
use mesh
use constitutive
!
implicit none
!
! *** Definition of variables ***
2007-03-26 15:57:34 +05:30
! *** Subroutine parameters ***
real(pReal) CPFEM_dt
integer(pInt) CPFEM_cn, cp_en ,CPFEM_in
2007-03-26 15:57:34 +05:30
! *** 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
2007-03-22 20:19:42 +05:30
nori = constitutive_Ngrains(CPFEM_in,cp_en) !<21><><EFBFBD>
!
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) !<21><><EFBFBD>
! *** Calculation of the solution at t=t1 ***
2007-03-22 20:19:42 +05:30
! QUESTION use the mod() as flag parameter in the call ??
if (mod(CPFEM_cn,ijaco)==0) then !<21><><EFBFBD>
2007-03-26 15:57:34 +05:30
call CPFEM_stress(cs, cd, CPFEM_dt,cp_en,CPFEM_in, iori, ising, icut, iconv, 1_pInt)
! *** Evaluation of ising ***
! *** ising=2 => singular matrix in jacobi calculation ***
! *** => use old jacobi ***
if (ising==2) jpara=0
! *** Calculation of the consistent tangent ***
CPFEM_d=CPFEM_d+vf*cd
else
2007-03-26 15:57:34 +05:30
call CPFEM_stress(cs, cd, CPFEM_dt,cp_en,CPFEM_in, iori, ising, icut, iconv, 0_pInt)
jpara=0
endif
! *** Cases of unsuccessful calculations ***
2007-03-22 20:19:42 +05:30
! *** 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)
2007-03-26 15:57:34 +05:30
! 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)
2007-03-26 15:57:34 +05:30
! CPFEM_timefactor=1.e5_pReal
return
endif
! *** Evaluation of iconv ***
! *** iconv!=0 => no convergence ***
if (iconv==1) then
2007-03-22 20:19:42 +05:30
write(6,*) 'Inner loop did not converge!'
write(6,*) 'Integration point: ',CPFEM_in
write(6,*) 'Element: ',CPFEM_en
call IO_error(600)
2007-03-26 15:57:34 +05:30
! CPFEM_timefactor=1.e5_pReal
return
else if (iconv==2) then
2007-03-22 20:19:42 +05:30
write(6,*) 'Outer loop did not converge!'
write(6,*) 'Integration point: ',CPFEM_in
write(6,*) 'Element: ',CPFEM_en
call IO_error(600)
2007-03-26 15:57:34 +05:30
! 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 ***
2007-03-26 15:57:34 +05:30
! *************************************
! *** End of the CP-FEM Calculation ***
! *************************************
2007-03-26 15:57:34 +05:30
! *** Store the new stress ***
CPFEM_stress_all(:,CPFEM_in,cp_en)=CPFEM_s
! *** Store the new jacobian ***
if (jpara/=0) CPFEM_jaco_old(:,:,CPFEM_in,cp_en)=CPFEM_d
return
2007-03-22 20:19:42 +05:30
end subroutine
2007-03-26 15:57:34 +05:30
!
!
subroutine CPFEM_stress(&
cs,& ! stress vector
cd,& ! Jacoby matrix
CPFEM_dt,& ! Time increment (dt)
cp_en,& ! Element number
CPFEM_in,& ! Integration point number
iori,& ! number of orintation
ising,& ! flag for singular matrix
icut,& ! flag for too many cut backs
iconv,& ! flag for non convergence
isjaco) ! flag whether to calculate Jacoby matrix
!********************************************************************
! This routine calculates the stress for a single component
! and manages the independent time incrmentation
!********************************************************************
use prec, only: pReal,pInt, ncut
use constitutive, only: constitutive_Nstatevars, constitutive_state_old, constitutive_state_new, constitutive_Nresults,&
constitutive_results
implicit none
!
! *** Definition of variables ***
2007-03-26 15:57:34 +05:30
! *** 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
2007-03-26 15:57:34 +05:30
!
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)
!
2007-03-26 15:57:34 +05:30
! *** 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
2007-03-26 15:57:34 +05:30
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)!<21><><EFBFBD><EFBFBD>
return
endif
!
! *** Calculation of stress and resistences with a cut timestep ***
! *** when first try did not converge ***
jcut=1_pInt
dt_i=0.5_pReal*CPFEM_dt
delta_Fg=0.5_pReal*(CPFEM_ffn1-CPFEM_ffn)
Fg_i=CPFEM_ffn+delta_Fg
Tstar_v=Tstar_v_h
state_new_i=state_new_h
! *** Start time ***
time=dt_i
do while (time<=CPFEM_dt)
call CPFEM_stress_int(cs, cd, time, cp_en,CPFEM_in, iori, ising, iconv, isjaco, phi1, PHI, phi2,&
CPFEM_ffn, Fg_i,Fp_old,Fp_new,state_old, state_new_i, Tstar_v)
if ((iconv==0).AND.(ising==0)) then
time=time+dt_i
Fg_i=Fg_i+delta_Fg
Tstar_v_h=Tstar_v
state_new_h=state_new_i
else
jcut=jcut+1_pInt
2007-03-26 15:57:34 +05:30
if (jcut>ncut) then
icut=1_pInt
return
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
2007-03-26 15:57:34 +05:30
state_new_i=state_new_h
endif
enddo
!
2007-03-22 20:19:42 +05:30
! *** Final calculation of stress and resistences with full timestep ***
state_new=state_new_i
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)
! *** 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
return
2007-03-22 20:19:42 +05:30
end subroutine
2007-03-26 15:57:34 +05:30
!
!
subroutine CPFEM_stress_int(&
cs,& ! Cauchy stress vector
dcs_de,& ! Consistent tangent
dt,& ! Time increment
2007-03-26 15:57:34 +05:30
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
2007-03-26 15:57:34 +05:30
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
! it is based on the paper by Kalidindi et al.:
! J. Mech. Phys, Solids Vol. 40, No. 3, pp. 537-569, 1992
! it is modified to use anisotropic elasticity matrix
!********************************************************************
use prec, only: pReal,pInt,pert_e
use constitutive, only: constitutive_Nstatevars
use math, only: math_Mandel6to33
implicit none
2007-03-26 15:57:34 +05:30
!
! *** Definition of variables ***
! *** Subroutine parameters ***
integer(pInt) cp_en, CPFEM_in, iori, ising, iconv, isjaco
2007-03-26 15:57:34 +05:30
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)
2007-03-26 15:57:34 +05:30
! *** 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
! *********************************************
! *** Calculation of the new Cauchy stress ***
! *********************************************
! *** Call Newton-Raphson method ***
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)
!
! *** Calculation of the new orientation ***
call math_pDecomposition(Fe,U,R,ising)
if (ising==1) then
return
endif
call math_RtoEuler(transpose(R),phi1,PHI,phi2)
!
! *** Choice of the calculation of the consistent tangent ***
if (isjaco==0) return
!
! *********************************************
! *** Calculation of the consistent tangent ***
! *********************************************
!
! *** Calculation of the consistent tangent with perturbation ***
! *** Perturbation on the component of Fg ***
do ic=1,6
!
! *** Method of small perturbation
! Missing direct matrix perturbation
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 ***
call NEWTON_RAPHSON(dt,cp_en,CPFEM_in,iori,Fg_pert,Fp_old,Fp2,Fe,state_old,state2,sgm2,cs1,iconv,ising)
!
! *** Consistent tangent *** as cs is Mandel dcs_de(:,4:6) is too large by sqrt(2)
dcs_de(:,ic)=(cs1-cs)/pert_e
enddo
!
return
2007-03-22 20:19:42 +05:30
end subroutine
!
!
2007-03-26 15:57:34 +05:30
subroutine NEWTON_RAPHSON(&
dt,&
cp_en,& ! Element number
CPFEM_in,& ! Integration point number
iori,& ! number of orintation
Fg_new,&
Fp_old,&
Fp_new,&
Fe,&
state_old,&
state_new,&
Tstar_v,&
cs,&
iconv,&
ising)
!***********************************************************************
!*** NEWTON-RAPHSON Calculation ***
!***********************************************************************
use prec, only: pReal,pInt, nouter, tol_outer, ninner, tol_inner, crite
use constitutive, only: constitutive_Nstatevars, constitutive_HomogenizedC, constitutive_dotState
2007-03-26 15:57:34 +05:30
use math
implicit none
! *** Definition of variables ***
! *** Subroutine parameters ***
integer(pInt) cp_en, CPFEM_in, iori, iconv, ising
real(pReal) dt,Fg_new(3,3),Fp_old(3,3),Fp_new(3,3), Fe(3,3)
real(pReal) state_old(constitutive_Nstatevars(iori, CPFEM_in, cp_en)), state_new(constitutive_Nstatevars(iori, CPFEM_in, cp_en))
2007-03-26 15:57:34 +05:30
real(pReal) Tstar_v(6), cs(6)
! *** Local variables ***
real(pReal) invFp_old(3,3), det, A(3,3), C_66(6,6), Lp(3,3), dLp(3,3,3,3)
real(pReal) I3tLp(3,3), help(3,3), help1(3,3,3,3), Tstar0_v(6), R1(6)
real(pReal) dstate(constitutive_Nstatevars(iori, CPFEM_in, cp_en)), R2(constitutive_Nstatevars(iori, CPFEM_in, cp_en))
real(pReal) R2s(constitutive_Nstatevars(iori, CPFEM_in, cp_en)), invFp_new(3,3)
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
!
! *** Error treatment ***
iconv = 0
ising = 0
!
2007-03-26 15:57:34 +05:30
! initialize new state
state_new=state_old
! *** Calculation of Fp_old(-1) ***
call invert3x3(Fp_old, invFp_old, det, err) !<21><><EFBFBD>
if (err==1_pInt) then
ising=1
return
endif
!
! *** Calculation of A and T*0 (see Kalidindi) ***
2007-03-28 19:02:25 +05:30
A = matmul(Fg_new,invFp_old) ! actually Fe
2007-03-26 15:57:34 +05:30
A = matmul(transpose(A), A)
C_66=constitutive_HomogenizedC(iori, CPFEM_in, cp_en) !<21><><EFBFBD>
2007-03-28 19:02:25 +05:30
Tstar_v=matmul(C_66, math_Mandel33to6(A-math_I3)) ! fully elastic guess
! QUESTION follow former plastic slope to guess better?
!
! *** Second level of iterative procedure: Resistences ***
do iouter=1,nouter
! *** First level of iterative procedure: Stresses ***
do iinner=1,ninner
!
! *** Calculation of gdot_slip ***
call constitutive_LpAndItsTangent(Tstar_v, iori, CPFEM_in, cp_en, Lp, dLp)
2007-03-26 15:57:34 +05:30
I3tLp = math_I3-dt*Lp
help=matmul(transpose(I3tLp),matmul(A, I3tLp))-math_I3
Tstar0_v = 0.5_pReal * matmul(C_66, math_Mandel33to6(help))
R1=Tstar_v-Tstar0_v
2007-03-28 19:02:25 +05:30
if (maxval(abs(R1/maxval(abs(Tstar_v)))) < tol_inner) goto 100
!
2007-03-28 19:02:25 +05:30
! *** Jacobi Calculation: dRes/dTstar ***
2007-03-26 20:21:01 +05:30
help=matmul(A, I3tLp)
help1=0
2007-03-28 19:42:41 +05:30
do i=1,3
do j=1,3
do k=1,3
do l=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)
enddo
enddo
enddo
enddo
enddo
help2=math_Mandel3333to66(help1)
Jacobi= 0.5_pReal*matmul(C_66, help2) + math_identity2nd(6)
2007-03-26 15:57:34 +05:30
call math_invert6x6(Jacobi, invJacobi, dummy, err) !<21><><EFBFBD>
if (err==1_pInt) then
forall (i=1:6) Jacobi(i,i)=1.05d0*maxval(Jacobi(i,:)) ! regularization
2007-03-26 15:57:34 +05:30
call math_invert6x6(Jacobi, invJacobi, dummy, err)
if (err==1_pInt) then ! sorry, can't help here!!
ising=1
return
endif
endif
dTstar_v=matmul(invJacobi,R1) ! correction to Tstar
! *** Correction (see Kalidindi) ***
2007-03-28 19:02:25 +05:30
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))
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)
! *** Arrays of residuals ***
R2=state_new-state_old-dstate
2007-03-28 19:02:25 +05:30
R2s=0.0_pReal
forall(i=1:constitutive_Nstatevars(iori, CPFEM_in, cp_en), state_new(i)/=0.0_pReal) R2s(i)=R2(i)/state_new(i)
if (maxval(abs(R2s)) < tol_outer) goto 200
state_new=state_old+dstate
enddo
iconv=2
return
! *** End of the second level of iterative procedure ***
! *** Calculation of Fp(t+dt) (see Kalidindi) ***
200 invFp_new=matmul(Fp_old, I3tLp)
2007-03-26 20:21:01 +05:30
call math_invert3x3(invFp_new, Fp_new, det, err) !<21><><EFBFBD>
if (err==1_pInt) then
ising=1
return
endif
2007-03-22 20:19:42 +05:30
Fp_new=Fp_new/math_det3x3(Fp_new)**(1.0_pReal/3.0_pReal)
!
! *** 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
2007-03-26 15:57:34 +05:30
end subroutine
!
function CPFEM_cauchy_stress(PK_v, Fe)
!***********************************************************************
!*** Cauchy stress calculation ***
!***********************************************************************
use prec, only: pReal,pInt
2007-03-28 19:42:41 +05:30
use math, only: math_Mandel33to6,math_Mandel6to33,math_det3x3
implicit none
! *** Subroutine parameters ***
real(pReal) PK_v(6), Fe(3,3), CPFEM_cauchy_stress(6)
2007-03-28 19:02:25 +05:30
CPFEM_cauchy_stress = math_Mandel33to6(matmul(matmul(Fe,math_Mandel6to33(PK_v)),transpose(Fe))/math_det3x3(Fe))
end function
end module