!!!!! IMPORTANT !!!!!

All subroutines were committed at once:
- constitutive_pheno works
- constitutive_dislo without twinning also works

This release should serve as reference
This commit is contained in:
Luc Hantcherli 2007-12-14 13:36:04 +00:00
parent 3ef451824c
commit fe70a82d6d
11 changed files with 72 additions and 91 deletions

View File

@ -116,7 +116,6 @@
real(pReal) ffn(3,3), ffn1(3,3), Temperature, CPFEM_dt, CPFEM_stress(CPFEM_ngens), CPFEM_jaco(CPFEM_ngens,CPFEM_ngens)
logical CPFEM_stress_recovery
!
Temperature = 293.0_pReal
! calculate only every second cycle
if(mod(CPFEM_cn,2)==0) then
! really calculate only in first call of new cycle and when in stress recovery
@ -151,7 +150,6 @@ if(mod(CPFEM_cn,2)==0) then
! this shall be done in a parallel loop in the future
do e=1,mesh_NcpElems
do i=1,FE_Nips(FE_mapElemtype(mesh_element(2,e)))
write(6,*) 'Debut', e,i
call CPFEM_stressIP(CPFEM_cn, CPFEM_dt, i, e)
enddo
enddo
@ -330,8 +328,6 @@ if(mod(CPFEM_cn,2)==0) then
real(pReal), dimension(6,6) :: dcs_de
real(pReal), dimension(constitutive_Nstatevars(grain,CPFEM_in,cp_en)) :: state_old,state_new,state_pert
write(6,*)
write(6,*) 'SOLUTION'
call CPFEM_timeIntegration(msg,Fp_new,Fe_new,Tstar_v,state_new, & ! def gradients and PK2 at end of time step
dt,cp_en,CPFEM_in,grain,Fg_new,Fp_old,state_old)
@ -346,9 +342,7 @@ if(mod(CPFEM_cn,2)==0) then
Fg_pert = Fg_new+matmul(E_pert,Fg_old) ! perturbated Fg
Tstar_v_pert = Tstar_v ! initial guess from end of time step
state_pert = state_new ! initial guess from end of time step
write(6,*)
write(6,*) 'CONSISTENT TANGENT', i
state_pert = state_new ! initial guess from end of time step
call CPFEM_timeIntegration(msg,Fp_pert,Fe_pert,Tstar_v_pert,state_pert, &
dt,cp_en,CPFEM_in,grain,Fg_pert,Fp_old,state_old)
if (msg /= 'ok') then
@ -360,14 +354,6 @@ if(mod(CPFEM_cn,2)==0) then
enddo
endif
write(6,*) 'OPERATEUR TANGENTE'
write(6,*) dcs_de(1,:)
write(6,*) dcs_de(2,:)
write(6,*) dcs_de(3,:)
write(6,*) dcs_de(4,:)
write(6,*) dcs_de(5,:)
write(6,*) dcs_de(6,:)
return
END SUBROUTINE
@ -394,7 +380,7 @@ if(mod(CPFEM_cn,2)==0) then
use prec
use constitutive, only: constitutive_Nstatevars,&
constitutive_homogenizedC,constitutive_dotState,constitutive_LpAndItsTangent,&
constitutive_Microstructure
constitutive_Microstructure
use math
implicit none
@ -437,10 +423,10 @@ state: do ! outer iteration: state
endif
call constitutive_Microstructure(state_new,CPFEM_Temperature(CPFEM_in,cp_en),grain,CPFEM_in,cp_en)
C_66 = constitutive_HomogenizedC(state_new, grain, CPFEM_in, cp_en)
iStress = 0_pInt
stress: do ! inner iteration: stress
iStress = iStress+1
write(6,*) 'istate,istress', istate, istress
if (iStress > nStress) then ! too many loops required
msg = 'limit stress iteration'
return
@ -462,9 +448,6 @@ stress: do ! inner iteration: stress
dTstar_v=0.5*dTstar_v
cycle
endif
write(6,*) 'Tstar_v', Tstar_v
write(6,*) 'dTstar_v', dTstar_v
!write(6,*) 'norm stress', maxval(abs(Rstress/maxval(abs(Tstar_v))))
if (iStress > 1 .and. &
(maxval(abs(Tstar_v)) < abstol_Stress .or. maxval(abs(Rstress/maxval(abs(Tstar_v)))) < reltol_Stress)) exit stress
@ -480,9 +463,9 @@ stress: do ! inner iteration: stress
enddo
enddo
enddo
enddo
Jacobi = math_identity2nd(6) + 0.5_pReal*dt*matmul(C_66,math_Mandel3333to66(LTL))
j = 0_pInt
enddo
Jacobi = math_identity2nd(6) + 0.5_pReal*dt*matmul(C_66,math_Mandel3333to66(LTL))
j = 0_pInt
call math_invert6x6(Jacobi,invJacobi,dummy,failed)
do while (failed .and. j <= nReg)
forall (i=1:6) Jacobi(i,i) = 1.05_pReal*maxval(Jacobi(i,:)) ! regularization
@ -496,10 +479,14 @@ stress: do ! inner iteration: stress
dTstar_v = matmul(invJacobi,Rstress) ! correction to Tstar
Rstress_old=Rstress
Tstar_v = Tstar_v-dTstar_v
enddo stress
Tstar_v = 0.5_pReal*matmul(C_66,math_Mandel33to6(matmul(transpose(B),AB)-math_I3))
!if ((printer==1_pInt).AND.(CPFEM_in==1_pInt).AND.(cp_en==1_pInt)) then
!write(6,'(A10, 12ES12.3)') 'state_new', state_new
!write(6,'(A10, 6ES12.3)') 'Tstar_v', Tstar_v
!endif
dstate = dt*constitutive_dotState(Tstar_v,state_new,CPFEM_Temperature(CPFEM_in,cp_en),grain,CPFEM_in,cp_en) ! evolution of microstructure
Rstate = state_new - (state_old+dstate)
RstateS = 0.0_pReal
@ -507,12 +494,9 @@ stress: do ! inner iteration: stress
RstateS(i) = Rstate(i)/state_new(i)
state_new = state_old+dstate
write(6,*) 'norm state', maxval(abs(RstateS))
if (maxval(abs(RstateS)) < reltol_State) exit state
enddo state
! write(999,*) 'Tstar_v raus', Tstar_v
! write(999,*)
invFp_new = matmul(invFp_old,B)
call math_invert3x3(invFp_new,Fp_new,det,failed)

View File

@ -499,10 +499,10 @@
END FUNCTION
!********************************************************************
!*********************************************************************
! read consecutive lines of ints concatenated by "c" as last char
! or range of values a "to" b
!********************************************************************
!*********************************************************************
FUNCTION IO_continousIntValues (unit,maxN,lookupName,lookupMap,lookupMaxN)
use prec, only: pReal,pInt

View File

@ -349,7 +349,7 @@ do while(.true.)
case ('burgers')
material_bg(section)=IO_floatValue(line,positions,2)
write(6,*) 'burgers', material_bg(section)
case ('Qedge')
case ('qedge')
material_Qedge(section)=IO_floatValue(line,positions,2)
write(6,*) 'Qedge', material_Qedge(section)
case ('tau0')
@ -925,11 +925,15 @@ constitutive_rho_p=matmul(constitutive_Pparallel(1:material_Nslip(matID),1:mater
do i=1,material_Nslip(matID)
constitutive_passing_stress(i)=material_tau0(matID)+material_c1(matID)*material_Gmod(matID)*material_bg(matID)*&
sqrt(constitutive_rho_p(i))
constitutive_jump_width(i)=material_c2(matID)/sqrt(constitutive_rho_f(i))
constitutive_activation_volume(i)=material_c3(matID)*constitutive_jump_width(i)*material_bg(matID)**2.0_pReal
constitutive_rho_m(i)=(2.0_pReal*Kb*Tp*sqrt(constitutive_rho_p(i)))/&
(material_c1(matID)*material_c3(matID)*material_Gmod(matID)*constitutive_jump_width(i)*&
material_bg(matID)**3.0_pReal)
constitutive_g0_slip(i)=constitutive_rho_m(i)*material_bg(matID)*attack_frequency*constitutive_jump_width(i)*&
exp(-(material_Qedge(matID)+constitutive_passing_stress(i)*constitutive_activation_volume(i))/&
(Kb*Tp))
@ -949,15 +953,6 @@ do i=1,material_Ntwin(matID)
constitutive_twin_volume(i)=(pi/6.0_pReal)*material_StackSize(matID)*constitutive_twin_mfp(i)**2.0_pReal
enddo
!write(6,*) 'rho_f', constitutive_rho_f
!write(6,*) 'rho_p', constitutive_rho_p
!write(6,*) 'passing', constitutive_passing_stress
!write(6,*) 'jump_width', constitutive_jump_width
!write(6,*) 'activation_volume', constitutive_activation_volume
!write(6,*) 'Temperature', Tp
!write(6,*) 'rho_m', constitutive_rho_m
!write(6,*) 'g0_slip', constitutive_g0_slip
return
end subroutine
@ -994,18 +989,19 @@ real(pReal), dimension(constitutive_Nstatevars(ipc,ip,el)) :: state
matID = constitutive_matID(ipc,ip,el)
startIdxTwin = material_Nslip(matID)
!* Recompute arrays from the microstructure (may be not needed)
call constitutive_Microstructure(state,Tp,ipc,ip,el)
!* Calculation of Lp - slip
Ftwin = sum(state((startIdxTwin+1):(startIdxTwin+material_Ntwin(matID))))
Lp = 0.0_pReal
do i=1,material_Nslip(matID)
constitutive_tau_slip(i)=dot_product(Tstar_v,crystal_Sslip_v(:,i,material_CrystalStructure(matID)))
constitutive_gdot_slip(i)=constitutive_g0_slip(i)*sinh((abs(constitutive_tau_slip(i))*&
constitutive_activation_volume(i))/(Kb*Tp))*sign(1.0_pReal,constitutive_tau_slip(i))
constitutive_dgdot_dtauslip(i)=((constitutive_g0_slip(i)*constitutive_activation_volume(i))/(Kb*Tp))*&
cosh((abs(constitutive_tau_slip(i))*constitutive_activation_volume(i))/(Kb*Tp))
if (abs(constitutive_tau_slip(i))<constitutive_passing_stress(i)) then
constitutive_gdot_slip(i) = 0.0_pReal
constitutive_dgdot_dtauslip(i) = 0.0_pReal
else
constitutive_gdot_slip(i)=constitutive_g0_slip(i)*sinh(constitutive_tau_slip(i)*constitutive_activation_volume(i)/Kb/Tp)
constitutive_dgdot_dtauslip(i)=constitutive_g0_slip(i)*constitutive_activation_volume(i)/Kb/Tp*&
cosh(constitutive_tau_slip(i)*constitutive_activation_volume(i)/Kb/Tp)
endif
Lp=Lp+(1.0_pReal-Ftwin)*constitutive_gdot_slip(i)*crystal_Sslip(:,:,i,material_CrystalStructure(matID))
enddo
@ -1038,8 +1034,6 @@ do i=1,material_Ntwin(matID)
crystal_Stwin(:,:,i,material_CrystalStructure(matID))
enddo
!write(6,*) 'constitutive_gdot_slip'
!write(6,*) constitutive_gdot_slip
!* Calculation of the tangent of Lp
dLp_dTstar=0.0_pReal
@ -1100,9 +1094,6 @@ real(pReal), dimension(constitutive_Nstatevars(ipc,ip,el)) :: constitutive_dotSt
!* Get the material-ID from the triplet(ipc,ip,el)
matID = constitutive_matID(ipc,ip,el)
!* Recompute arrays from the microstructure (may be not needed)
call constitutive_Microstructure(state,Tp,ipc,ip,el)
!* Hardening of each system
do i=1,material_Nslip(matID)
tau_slip = dot_product(Tstar_v,crystal_Sslip_v(:,i,material_CrystalStructure(matID)))

View File

@ -84,9 +84,9 @@ real(pReal), dimension(:,:,:) , allocatable :: constitutive_texVolFrac
!* State variables *
!************************************
integer(pInt) constitutive_maxNstatevars
integer(pInt), dimension(:,:,:), allocatable :: constitutive_Nstatevars
real(pReal), dimension(:,:,:,:), allocatable :: constitutive_state_old
real(pReal), dimension(:,:,:,:), allocatable :: constitutive_state_new
integer(pInt), dimension(:,:,:) , allocatable :: constitutive_Nstatevars
real(pReal), dimension(:,:,:,:) , allocatable :: constitutive_state_old
real(pReal), dimension(:,:,:,:) , allocatable :: constitutive_state_new
!************************************
!* Hardening matrices *

View File

@ -310,7 +310,7 @@ implicit none
!* Definition of variables
integer(pInt) i,j,k,l
real(pReal) invNorm
real(pReal) norm_d,norm_t,norm_n
!* Iteration over the crystal structures
do l=1,crystal_MaxCrystalStructure
@ -320,14 +320,16 @@ do l=1,crystal_MaxCrystalStructure
crystal_st(1,k,l)=crystal_sn(2,k,l)*crystal_sd(3,k,l)-crystal_sn(3,k,l)*crystal_sd(2,k,l)
crystal_st(2,k,l)=crystal_sn(3,k,l)*crystal_sd(1,k,l)-crystal_sn(1,k,l)*crystal_sd(3,k,l)
crystal_st(3,k,l)=crystal_sn(1,k,l)*crystal_sd(2,k,l)-crystal_sn(2,k,l)*crystal_sd(1,k,l)
norm_d=dsqrt(crystal_sd(1,k,l)**2+crystal_sd(2,k,l)**2+crystal_sd(3,k,l)**2)
norm_t=dsqrt(crystal_st(1,k,l)**2+crystal_st(2,k,l)**2+crystal_st(3,k,l)**2)
norm_n=dsqrt(crystal_sn(1,k,l)**2+crystal_sn(2,k,l)**2+crystal_sn(3,k,l)**2)
crystal_sd(:,k,l)=crystal_sd(:,k,l)/norm_d
crystal_st(:,k,l)=crystal_st(:,k,l)/norm_t
crystal_sn(:,k,l)=crystal_sn(:,k,l)/norm_n
!* Defintion of Schmid matrix
forall (i=1:3,j=1:3)
crystal_Sslip(i,j,k,l)=crystal_sd(i,k,l)*crystal_sn(j,k,l)
endforall
!* Normalization of Schmid matrix
invNorm=dsqrt(1.0_pReal/((crystal_sn(1,k,l)**2+crystal_sn(2,k,l)**2+crystal_sn(3,k,l)**2)*&
(crystal_sd(1,k,l)**2+crystal_sd(2,k,l)**2+crystal_sd(3,k,l)**2)))
crystal_Sslip(:,:,k,l)=crystal_Sslip(:,:,k,l)*invNorm
!* Vectorization of normalized Schmid matrix
crystal_Sslip_v(1,k,l)=crystal_Sslip(1,1,k,l)
crystal_Sslip_v(2,k,l)=crystal_Sslip(2,2,k,l)
@ -344,16 +346,18 @@ do l=1,crystal_MaxCrystalStructure
crystal_tt(1,k,l)=crystal_tn(2,k,l)*crystal_td(3,k,l)-crystal_tn(3,k,l)*crystal_td(2,k,l)
crystal_tt(2,k,l)=crystal_tn(3,k,l)*crystal_td(1,k,l)-crystal_tn(1,k,l)*crystal_td(3,k,l)
crystal_tt(3,k,l)=crystal_tn(1,k,l)*crystal_td(2,k,l)-crystal_tn(2,k,l)*crystal_td(1,k,l)
norm_d=dsqrt(crystal_td(1,k,l)**2+crystal_td(2,k,l)**2+crystal_td(3,k,l)**2)
norm_t=dsqrt(crystal_tt(1,k,l)**2+crystal_tt(2,k,l)**2+crystal_tt(3,k,l)**2)
norm_n=dsqrt(crystal_tn(1,k,l)**2+crystal_tn(2,k,l)**2+crystal_tn(3,k,l)**2)
crystal_td(:,k,l)=crystal_td(:,k,l)/norm_d
crystal_tt(:,k,l)=crystal_tt(:,k,l)/norm_t
crystal_tn(:,k,l)=crystal_tn(:,k,l)/norm_n
!* Defintion of Schmid matrix and transformation matrices
crystal_Qtwin(:,:,k,l)=-math_I3
forall (i=1:3,j=1:3)
crystal_Stwin(i,j,k,l)=crystal_td(i,k,l)*crystal_tn(j,k,l)
crystal_Qtwin(i,j,k,l)=crystal_Qtwin(i,j,k,l)+2*crystal_tn(i,k,l)*crystal_tn(j,k,l)
endforall
!* Normalization of Schmid matrix
invNorm=dsqrt(1.0_pReal/((crystal_tn(1,k,l)**2+crystal_tn(2,k,l)**2+crystal_tn(3,k,l)**2)*&
(crystal_td(1,k,l)**2+crystal_td(2,k,l)**2+crystal_td(3,k,l)**2)))
crystal_Stwin(:,:,k,l)=crystal_Stwin(:,:,k,l)*invNorm
!* Vectorization of normalized Schmid matrix
crystal_Stwin_v(1,k,l)=crystal_Stwin(1,1,k,l)
crystal_Stwin_v(2,k,l)=crystal_Stwin(2,2,k,l)

View File

@ -114,7 +114,7 @@
a(k,i) = a(k,j)
a(k,j) = tmp
enddo
else ! if they do cross, exchange left value with pivot and return with the partition index
else ! if they do cross, exchange left value with pivot and return with the partition index
do k = 1,d
tmp = a(k,istart)
a(k,istart) = a(k,j)

View File

@ -1,28 +1,28 @@
<materials>
[Aluminium]
[TWIP steel FeMnC]
crystal_structure 1
Nslip 12
Ntwin 0
## Elastic constants
# Unit in [Pa]
C11 106.75e9
C12 60.41e9
C44 28.34e9
C11 245.0e9
C12 105.0e9
C44 65.0e9
## Parameters for phenomenological modeling
# Unit in [Pa]
s0_slip 31.0e6
s0_slip 85.0e6
gdot0_slip 0.001
n_slip 20
h0 75.0e6
s_sat 63.0e6
w0 2.25
n_slip 100.0
h0 355.0e6
s_sat 265.0e6
w0 1.0
# Self and latent hardening coefficients
hardening_coefficients 1.0 1.4
## Parameters for dislocation-based modeling
# Initial dislocation density [m]²
rho0 1.5e11
rho0 2.8e13
# Burgers vector [m]
burgers 2.86e-10
# Activation energy for dislocation glide [J/K]
@ -34,31 +34,33 @@ c1 0.1
# Jump width adjustment
c2 2.0
# Activation volume adjustment
c3 0.4
c3 1.2
# Dislocation storage adjustment
c4 20.0
# = c4(Anxin)*c2(Anxin) !!!!!!
c4 14.29
# Athermal annihilation adjustment
c7 20.0
# Dislocation interaction coefficients
interaction_coefficients 1.0 2.2 3.0 1.6 3.8 4.5
# Twin parameters
# Grain boundaries storage adjustment
c5 1.0e100
# Twin boundaries storage adjustment
c6 1.0e100
# Athermal annihilation adjustment
c7 10.0
# Dislocation interaction coefficients
interaction_coefficients 1.0 2.2 3.0 1.6 3.8 4.5
# Twin parameters
# grain size, average size of stacks of twins [m]
grain_size 1.5e-5
stack_size 5.0e-8
# stacking fault energy
stacking_fault_energy 2.0e-2
# Twin reference [?], twin resistance [Pa], twin sensitivity
twin_ref 1.0e-15
twin_res 150.0e6
twin_sens 10.0
twin_reference 1.0e-15
twin_resistance 150.0e6
twin_sensitivity 10.0
<textures>
[cube SX]
symmetry monoclinic /orthorhombic
Ngrains 2 /4
(gauss) phi1 0.0 phi 0.0 phi2 0.0 scatter 0.0 fraction 1.0
symmetry no /monoclinic /orthorhombic
Ngrains 1 /2 /4
(gauss) phi1 0.0 phi 45.0 phi2 0.0 scatter 0.0 fraction 1.0

View File

@ -830,7 +830,7 @@ matchFace: do j = 1,FE_NfaceNodes(-neighbor,t) ! count over nodes on matc
write (6,*)
write (fmt,"(a,i3,a,a)") "(",1+mesh_maxValStateVar(1),"(i8)",")"
do i=1,mesh_maxValStateVar(1) ! loop over all (possibly assigned) materials
do i=1,mesh_maxValStateVar(1) ! loop over all (possibly assigned) materials
write (6,fmt) i,mesh_MatTex(i,:) ! loop over all (possibly assigned) textures
enddo
write (6,*)

View File

@ -34,7 +34,7 @@
include "crystal.f90"
include "constitutive.f90"
include "CPFEM.f90"
!
SUBROUTINE hypela2(d,g,e,de,s,t,dt,ngens,n,nn,kcus,matus,ndi,&
nshear,disp,dispt,coord,ffn,frotn,strechn,eigvn,ffn1,&
frotn1,strechn1,eigvn1,ncrd,itel,ndeg,ndm,&

View File

@ -119,7 +119,7 @@
!2 continue
!3 continue
!
!
use prec, only: pReal,pInt
use CPFEM, only: CPFEM_general
implicit real(pReal) (a-h,o-z)

View File

@ -9,7 +9,7 @@
integer, parameter :: pInt = 4
! *** Numerical parameters ***
! *** How frequently the jacobian is recalculated ***
integer (pInt), parameter :: ijaco = 5_pInt
integer (pInt), parameter :: ijaco = 1_pInt
! *** Maximum number of internal cutbacks in time step ***
integer(pInt), parameter :: nCutback = 7_pInt
! *** Maximum number of regularization attempts for Jacobi inversion ***