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

changed interfacing to marc due to Mandel notation 

set relative convergence limits
This commit is contained in:
Franz Roters 2007-03-28 12:49:12 +00:00
parent 2be98a39ac
commit 05db614589
2 changed files with 47 additions and 41 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

67
trunk/CPFEM.f90
View File

@ -1,4 +1,4 @@
! last modified 26.03.07
! last modified 28.03.07
! ---------------------------
MODULE CPFEM
! ---------------------------
@ -399,7 +399,7 @@
!********************************************************************
use prec, only: pReal,pInt
use constitutive, only: constitutive_Nstatevars
use math, only: math_6to33
use math, only: math_Mandel6to33
implicit none
!
! *** Definition of variables ***
@ -447,7 +447,7 @@
dev=0
if(ic<=3) dev(ic) = pert_ct
if(ic>3) dev(ic) = pert_ct/2
dF=matmul(math_6to33(dev),Fg_old)
dF=matmul(math_Mandel6to33(dev),Fg_old)
Fg2=Fg_new+dF
sgm2=Tstar_v
state2=state_new
@ -493,17 +493,18 @@
real(pReal) state_old(constitutive_Nstatevars(iori, CPFEM_in, cp_en)), state_new(constitutive_Nstatevars(iori, CPFEM_in, cp_en))
real(pReal) Tstar_v(6), cs(6)
! *** Local variables ***
real(pReal) crite, tol_in, tol_out, invFp_old(3,3), det, A(3,3), C_66(6,6), Lp(3,3), dLp(3,3,6)
real(pReal) I3tLp(3,3), help(3,3), help1(6), Tstar0_v(6), R1(6), norm1, tdLp(3,3)
real(pReal) dstate(constitutive_Nstatevars(iori, CPFEM_in, cp_en)), R2(6), norm2, invFp_new(3,3), Estar(3,3)
real(pReal) Estar_v(6), Jacobi(6,6), invJacobi(6,6), dTstar_v(6)
integer(pInt) iouter, iinner , dummy, err, i, j, k
real(pReal) tol_in, tol_out, 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), R1s(6), norm1, tdLp(3,3)
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)), norm2, invFp_new(3,3), Estar(3,3)
real(pReal) Estar_v(6), Jacobi(6,6), invJacobi(6,6), dTstar_v(6), help2(6,6)
integer(pInt) iouter, iinner , dummy, err, i, j, k, l, m
! *** Numerical parameters ***
integer(pInt), parameter :: nouter = 50_pInt
real(pReal), parameter :: tol_outer = 1.0e-4_pReal
integer(pInt), parameter :: ninner = 2000_pInt
real(pReal), parameter :: tol_inner = 1.0e-3_pReal
real(pReal), parameter :: eta = 13.7_pReal
real(pReal), parameter :: crite = 0.3_pReal
! crite=eta*constitutive_s0_slip/constitutive_n_slip !ÄÄÄ
!
@ -528,6 +529,7 @@
A = matmul(Fg_new,invFp_old)
A = matmul(transpose(A), A)
C_66=constitutive_HomogenizedC(iori, CPFEM_in, cp_en) !ÄÄÄ
Tstar_v=matmul(C_66, math_Mandel33to6(A-math_I3))
!
! *** Second level of iterative procedure: Resistences ***
do iouter=1,nouter
@ -535,26 +537,24 @@
do iinner=1,ninner
!
! *** Calculation of gdot_slip ***
call constitutive_LpAndItsTangent(Lp, dLp, iori, CPFEM_in, cp_en)
call constitutive_LpAndItsTangent(Tstar_v, iori, CPFEM_in, cp_en, Lp, dLp)
I3tLp = math_I3-dt*Lp
help=matmul(transpose(I3tLp),matmul(A, I3tLp))-math_I3
Tstar0_v = 0.5_pReal * matmul(C_66, math_33to6(help))
R1=Tstar_v-Tstar0_v
norm1=maxval(abs(R1))
if (norm1<tol_in) goto 100
Tstar0_v = 0.5_pReal * matmul(C_66, math_Mandel33to6(help))
R1=Tstar_v-Tstar0_v
R1s=0
forall(i=1:6, Tstar_v(i)/=0) R1s(i)=R1(i)/Tstar_v(i)
norm1=maxval(abs(R1s))
if (norm1<tol_inner) goto 100
!
! *** Jacobi Calculation ***
help=matmul(A, I3tLp)
! MISSING 1..6 outer loop, inner sum required
help1=0
do k=1,6
do j=1,3
do i=1,6
help1(k)=help1(k)+dLp(j,i,k)*help(i,j)+dLp(i,j,k)*help(j,i)
enddo
enddo
enddo
! Jacobi= matmul(C_66, help1) + math_identity(6)
forall(i=1:3, j=1:3, k=1:3, l=1:3,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)
help2=math_Mandel3333to66(help1)
Jacobi= 0.5_pReal*matmul(C_66, help2) + math_identity2nd(6)
call math_invert6x6(Jacobi, invJacobi, dummy, err) !ÄÄÄ
if (err==1_pInt) then
forall (i=1:6) Jacobi(i,i)=1.05d0*maxval(Jacobi(i,:)) ! regularization
@ -566,12 +566,13 @@
endif
dTstar_v=matmul(invJacobi,R1) ! correction to Tstar
! *** Correction (see Kalidindi) ***
do i=1,6
if (abs(dTstar_v(i))>crite) then
dTstar_v(i)=sign(crite,dTstar_v(i))
endif
enddo
! *** Correction (see Kalidindi) ***
forall(i=1:6, abs(dTstar_v(i)) > crite*Tstar_v(i).AND. Tstar_v(i)/=0)&
! do i=1,6
! if (abs(dTstar_v(i))> crite*Tstar_v(i)) then
dTstar_v(i)=sign(crite*Tstar_v(i),dTstar_v(i))
! endif
! enddo
Tstar_v=Tstar_v-dTstar_v
!
enddo
@ -582,8 +583,10 @@
100 dstate=dt*constitutive_dotState(Tstar_v, iori, CPFEM_in, cp_en)
! *** Arrays of residuals ***
R2=state_new-state_old-dstate
norm2=maxval(abs(R2))
if (norm2<tol_out) goto 200
R2s=0
forall(i=1:constitutive_Nstatevars(iori, CPFEM_in, cp_en), state_new(i)/=0) R2s(i)=R2(i)/state_new(i)
norm2=maxval(abs(R2s))
if (norm2<tol_outer) goto 200
state_new=state_old+dstate
enddo
iconv=2
@ -624,8 +627,8 @@
real(pReal), dimension(3,3) :: cs_33
! *** Calculation of Estar ***
det = math_det3x3(Fe)
PK = math_6to33(PK_v)
PK = math_Mandel6to33(PK_v)
cs_33 = matmul(matmul(Fe,PK),transpose(Fe))/det
CPFEM_cauchy_stress = math_33to6(cs_33)
CPFEM_cauchy_stress = math_Mandel33to6(cs_33)
end function
end module

21
trunk/mpie_cpfem_marc.f90
View File

@ -119,7 +119,8 @@
!
!
use prec, only: pReal,pInt
use CPFEM, only : CPFEM_stress_all, CPFEM_jaco_old
use CPFEM, only: CPFEM_stress_all, CPFEM_jaco_old
use math, only: invnrmMandel, nrmMandel
implicit real(pReal) (a-h,o-z)
!
! Marc common blocks are in fixed format so they have to be pasted in here beware of changes in newer Marc versions
@ -156,10 +157,7 @@
!
! call general material routine only in increment 0 and for lovl==6 (stress recovery)
! subroutine cpfem_general(mpie_s, mpie_d, mpie_ndi,
! 1 mpie_ffn, mpie_ffn1, mpie_cn, mpie_tinc,
! 2 mpie_timefactor, mpie_numel, mpie_nip, mpie_en,
! 3 mpie_in, mpie_mn, mpie_dimension, state_var)
! subroutine cpfem_general(mpie_ffn, mpie_ffn1, mpie_cn, mpie_tinc, mpie_enp, mpie_in)
!********************************************************************
! This routine calculates the material behaviour
!********************************************************************
@ -169,15 +167,20 @@
! mpie_tinc time increment
! mpie_en element number
! mpie_in intergration point number
! mpie_dimension dimension of stress/strain vector
!********************************************************************
cp_en=mesh_FEasCP('elem', n(1))
if ((lovl==6).or.(inc==0)) then
call cpfem_general(ffn, ffn1, inc, incsub, ncycle, timinc, cp_en, nn)
endif
! return stress and jacobi
s(1:ngens)=CPFEM_stress_all(1:ngens, nn, cp_en)
d(1:ngens,1:ngens)=CPFEM_jaco_old(1:ngens,1:ngens, nn, cp_en)
! return stress and jacobi
! Mandel: 11, 22, 33, 12, 23, 13
! Marc: 11, 22, 33, 12, 23, 13
s(1:ngens)=invnrmMandel(1:ngens)*CPFEM_stress_all(1:ngens, nn, cp_en)
d(1:ngens,1:ngens)=CPFEM_jaco_old(1:ngens,1:ngens, nn, cp_en)
forall(i=1:ngens) d(1:ngens,i)=d(1:ngens,i)*invnrmMandel(1:ngens)
d(1:ngens,1:ngens)=transpose(d(1:ngens,1:ngens))
forall(i=1:ngens) d(1:ngens,i)=d(1:ngens,i)*nrmMandel(1:ngens)
d(1:ngens,1:ngens)=transpose(d(1:ngens,1:ngens))
return
end