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

second initial import ;-)

This commit is contained in:
Philip Eisenlohr 2007-03-20 13:50:57 +00:00
parent 823c64c28a
commit 987c3b9d89
6 changed files with 3241 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

883
CPFEM.f90 Normal file
View File

@ -0,0 +1,883 @@
! ---------------------------
MODULE CPFEM
! ---------------------------
! *** CPFEM engine ***
use prec, only: pRe,pIn
implicit none
! ****************************************************************
! *** General variables for the material behaviour calculation ***
! ****************************************************************
real(pRe), allocatable :: CPFEM_stress_all (:,:,:)
real(pRe), allocatable :: CPFEM_jacobi_all (:,:,:,:)
real(pRe), allocatable :: CPFEM_results (:,:,:,:)
real(pRe), allocatable :: CPFEM_thickness (:,:)
real(pRe), allocatable :: CPFEM_ini_ori (:,:,:,:)
real(pRe), allocatable :: CPFEM_sigma_old (:,:,:,:)
real(pRe), allocatable :: CPFEM_sigma_new (:,:,:,:)
real(pRe), allocatable :: CPFEM_Fp_old (:,:,:,:,:)
real(pRe), allocatable :: CPFEM_Fp_new (:,:,:,:,:)
real(pRe), allocatable :: CPFEM_tauc_slip_old(:,:,:,:)
real(pRe), allocatable :: CPFEM_tauc_slip_new(:,:,:,:)
real(pRe), allocatable :: CPFEM_g_old (:,:,:,:)
real(pRe), allocatable :: CPFEM_g_new (:,:,:,:)
real(pRe), allocatable :: CPFEM_jaco_old (:,:,:,:)
real(pRe), allocatable :: CPFEM_mat (:,:)
CONTAINS
!***********************************************************************
!*** This routine allocates the arrays defined in module mpie ***
!*** and initializes them ***
!***********************************************************************
subroutine ALLOCATION(mpie_numel,mpie_nip)
use prec, only: pRe,pIn
use IO, only: _error
use math
use mesh
use constitutive
implicit none
integer(pIn) i
! *** mpie.marc parameters ***
allocate(CPFEM_stress_all(6,mesh_Nelems,mesh_Nips))
allocate(CPFEM_jacobi_all(6,6,mesh_Nelems,mesh_Nips))
CPFEM_stress_all=0.0_pRe
CPFEM_jacobi_all=0.0_pRe
! *** User defined results ***
allocate(CPFEM_results(constitutive_Nresults,
& constitutive_maxNgrains,
& mesh_Nelems,mesh_Nips))
CPFEM_results=0.0_pRe
! *** Relative sheet thickness ***
allocate(CPFEM_thickness(mesh_Nelems,mesh_Nips))
CPFEM_thickness=0.0_pRe
! *** Initial orientations ***
allocate(CPFEM_ini_ori(3,constitutive_maxNgrains,mesh_Nelems,mesh_Nips))
CPFEM_ini_ori=0.0_pRe
! *** Second Piola-Kirchoff stress tensor at (t=t0) and (t=t1) ***
allocate(CPFEM_sigma_old(6,constitutive_maxNgrains,mesh_Nelems,mesh_Nips))
allocate(CPFEM_sigma_new(6,constitutive_maxNgrains,mesh_Nelems,mesh_Nips))
CPFEM_sigma_old=0.0_pRe
CPFEM_sigma_new=0.0_pRe
! *** Plastic deformation gradient at (t=t0) and (t=t1) ***
allocate(CPFEM_Fp_old(3,3,constitutive_maxNgrains,mesh_Nelems,mesh_Nips))
allocate(CPFEM_Fp_new(3,3,constitutive_maxNgrains,mesh_Nelems,mesh_Nips))
CPFEM_Fp_old=0.0_pRe
CPFEM_Fp_new=0.0_pRe
do i=1,3
CPFEM_Fp_old(i,i,:,:,:)=1.0_pRe
CPFEM_Fp_new(i,i,:,:,:)=1.0_pRe
enddo
! QUESTION: would it be wise to outsource these to _constitutive_ ??
! *** Slip resistances at (t=t0) and (t=t1) ***
allocate(CPFEM_tauc_slip_old(nslip,constitutive_maxNgrains,mesh_Nelems,
& mesh_Nips))
allocate(CPFEM_tauc_slip_new(nslip,constitutive_maxNgrains,mesh_Nelems,
& mesh_Nips))
CPFEM_tauc_slip_old=0.0_pRe
CPFEM_tauc_slip_new=0.0_pRe
! *** Cumulative shear at (t=t0) and (t=t1) ***
! QUESTION which nslip to use here ?!?
allocate(CPFEM_g_old(nslip,constitutive_maxNgrains,mesh_Nelems,mesh_Nips))
allocate(CPFEM_g_new(nslip,constitutive_maxNgrains,mesh_Nelems,mesh_Nips))
CPFEM_g_old=0.0_pRe
CPFEM_g_new=0.0_pRe
! *** Old jacobian (consistent tangent) ***
allocate(CPFEM_jaco_old(6,6,mesh_Nelems,mesh_Nips))
! *** Output to MARC output file ***
write(6,*)
write(6,*) 'Arrays allocated:'
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_thickness: ', shape(CPFEM_thickness)
write(6,*) 'CPFEM_ini_ori: ', shape(CPFEM_ini_ori)
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_tauc_slip_old: ', shape(CPFEM_tauc_slip_old)
write(6,*) 'CPFEM_tauc_slip_new: ', shape(CPFEM_tauc_slip_new)
write(6,*) 'CPFEM_g_old: ', shape(CPFEM_g_old)
write(6,*) 'CPFEM_g_new: ', shape(CPFEM_g_new)
write(6,*) 'CPFEM_jaco_old: ', shape(CPFEM_jaco_old)
write(6,*)
call flush(6)
return
end
subroutine CPFEM_general_material(
& CPFEM_s, ! Stress vector
& CPFEM_d, ! Jacobi matrix (consistent tangent)
& CPFEM_ndi, ! Dimension
& CPFEM_ffn, ! Deformation gradient at begin of increment
& CPFEM_ffn1, ! Deformation gradient at end of increment
& CPFEM_inc, ! Increment number
& CPFEM_subinc, ! Subincrement number
& CPFEM_cn, ! Cycle number
& CPFEM_tinc, ! Time increment (dt)
& CPFEM_timefactor, ! Factor for timestep correction
! & mesh_Nelems, ! Number of elements in mesh
! & CPFEM_nip, ! Maximum number of integration points per element
& CPFEM_en, ! Element number
& CPFEM_in, ! Integration point number
& CPFEM_mn, ! Material number
& CPFEM_dimStress ! Dimension of stress/strain vector
&)
!***********************************************************************
!*** This routine calculates the material behaviour ***
!***********************************************************************
use prec, only: pRe,pIn
use IO, only _error
use math
use mesh
use constitutive
implicit none
! *** Definition of variables ***
integer(pIn) CPFEM_ndi,CPFEM_inc,CPFEM_subinc,CPFEM_cn,
& CPFEM_en,CPFEM_in,CPFEM_mn,CPFEM_dimStress
real(pRe) CPFEM_timefactor,CPFEM_tinc,CPFEM_s(CPFEM_dimStress),
& CPFEM_d(CPFEM_dimStress,CPFEM_dimStress),
& CPFEM_ffn(3,3),CPFEM_ffn1(3,3)
! QUESTION which nslip to use?
real(pRe) Fp_old(3,3),tauc_slip_old(nslip),
& tauc_slip_new(nslip),g_old(nslip),
& g_new(nslip),Tstar_v(6),
& Fp_new(3,3),cs(6),phi1mis(2),PHImis(2),phi2mis(2),
& cd(6,6),ori_mat(3,3),hh6(6,6)
integer(pIn) jpara,nori
real(pRe) phi1,PHI,phi2,scatter,vf,alpha1,alpha2,beta1,
& beta2,phi1_s,PHI_s,phi2_s,p10,P0,p20,p11,P1,p21,
& dgmax,dgmaxc,orimis
integer(pIn) i,iori,iconv,ising,icut
! *** Numerical parameters ***
! *** How often the jacobian is recalculated ***
integer (pIn), parameter :: ijaco=1_pIn
! *** Reference shear rate for the calculation of CPFEM_timefactor ***
real (pRe), parameter :: dgs=0.01_pRe
! *** Initialization step ***
if (CPFEM_first_call==1_pIn) then
call INITIALIZATION(mesh_Nelems,CPFEM_nip)
CPFEM_first_call=0_pIn
endif
! *** Case of a new increment ***
if (CPFEM_inc.NE.CPFEM_inc_old) then
CPFEM_sigma_old=CPFEM_sigma_new
CPFEM_Fp_old=CPFEM_Fp_new
CPFEM_tauc_slip_old=CPFEM_tauc_slip_new
CPFEM_g_old=CPFEM_g_new
CPFEM_inc_old=CPFEM_inc
CPFEM_subinc_old=1_pIn
CPFEM_timefactor_max=0.0_pRe
endif
! *** case of a new subincrement:update starting with subinc 2 ***
if (CPFEM_subinc.GT.CPFEM_subinc_old) then
CPFEM_sigma_old=CPFEM_sigma_new
CPFEM_Fp_old=CPFEM_Fp_new
CPFEM_tauc_slip_old=CPFEM_tauc_slip_new
CPFEM_g_old=CPFEM_g_new
CPFEM_subinc_old=CPFEM_subinc
endif
! *** Flag for recalculation of jacobian ***
jpara=1_pIn
! ************************************
! *** Orientation initialization ***
! ************************************
! *** Number of components per state ***
nori=CPFEM_mat(CPFEM_mn,1)
if (CPFEM_inc==0_pIn) then
! *** Three dimensional stress state ***
if (CPFEM_ndi.NE.3_pIn) then
call CPFEM_error(300)
endif
if ((CPFEM_en==1_pIn).AND.(CPFEM_in==1_pIn)) then
write(6,*) 'MPIE Material Routine Ver. 0.1 by L. Hantcherli'
write(6,*)
write(6,*) 'Orientation initialization'
call flush(6)
endif
i=1
do while (i.LE.nori)
! *** Direct ODF sampling ***
if (CPFEM_mat(CPFEM_mn,2)==2) then
call CPFEM_odf_ori(CPFEM_cko(CPFEM_mn,:,:,:,:),
& CPFEM_odfmax(CPFEM_mn),phi1,PHI,phi2)
else
! *** Gauss/Spherical component ***
if (CPFEM_mat(CPFEM_mn,7*i-4)==1) then
phi1=CPFEM_mat(CPFEM_mn,7*i-3)
PHI=CPFEM_mat(CPFEM_mn,7*i-2)
phi2=CPFEM_mat(CPFEM_mn,7*i-1)
scatter=CPFEM_mat(CPFEM_mn,7*i+1)
! *** Random orientation to this component to represent ***
! *** random fraction of texture using halton series ***
if (phi1==400.0) then
call CPFEM_halton_ori(phi1,PHI,phi2,scatter)
! *** ELSE modify orientation to represent gauss distribution ***
else if (scatter.GT.0.1) then
call CPFEM_gauss(phi1,PHI,phi2,scatter)
endif
! *** Fiber component ***
else if (CPFEM_mat(CPFEM_mn,7*i-4)==2) then
alpha1=CPFEM_mat(CPFEM_mn,7*i-3)
alpha2=CPFEM_mat(CPFEM_mn,7*i-2)
beta1=CPFEM_mat(CPFEM_mn,7*i-1)
beta2=CPFEM_mat(CPFEM_mn,7*i)
scatter=CPFEM_mat(CPFEM_mn,7*i+1)
! *** Random orientation to this component to represent ***
! *** random fraction of texture using random numbers ***
if (alpha1==400.0) then
call CPFEM_random_ori(phi1,PHI,phi2,scatter)
! *** ELSE calculate orientation to represent fiber component ***
else if (scatter.GT.0.1) then
call CPFEM_fiber(alpha1,alpha2,beta1,beta2,
& scatter,phi1,PHI,phi2)
endif
else
call CPFEM_error(510)
endif
endif
CPFEM_ini_ori(1,i,CPFEM_en,CPFEM_in)=phi1
CPFEM_ini_ori(2,i,CPFEM_en,CPFEM_in)=PHI
CPFEM_ini_ori(3,i,CPFEM_en,CPFEM_in)=phi2
! *** Orientation matrix ***
call CPFEM_euldreh(phi1,PHI,phi2,ori_mat)
CPFEM_Fp_old(:,:,i,CPFEM_en,CPFEM_in)=ori_mat
i=i+1
! *** If symmetric component, creation of additional three orientations ***
if (CPFEM_mat(CPFEM_mn,2)==1) then
! *** First one ***
phi1_s=180.0_pRe-phi1
if (phi1_s.LT.0.0_pRe) phi1_s=phi1_s+360.0_pRe
PHI_s=180.0_pRe-PHI
if (PHI_s.LT.0.0_pRe) PHI_s=PHI_s+360.0_pRe
phi2_s=phi2+180.0_pRe
if (phi2_s.GT.360.0_pRe) phi2_s=phi2_s-360.0_pRe
CPFEM_ini_ori(1,i,CPFEM_en,CPFEM_in)=phi1_s
CPFEM_ini_ori(2,i,CPFEM_en,CPFEM_in)=PHI_s
CPFEM_ini_ori(3,i,CPFEM_en,CPFEM_in)=phi2_s
! *** Orientation matrix for initial orientation ***
call CPFEM_euldreh(phi1_s,PHI_s,phi2_s,ori_mat)
CPFEM_Fp_old(:,:,i,CPFEM_en,CPFEM_in)=ori_mat
i=i+1
! *** Second one ***
phi1_s=360.0_pRe-phi1
PHI_s=180.0_pRe-PHI
if (PHI_s.LT.0.0_pRe) PHI_s=PHI_s+360.0_pRe
phi2_s=phi2+180.0_pRe
if (phi2_s.GT.360.0_pRe) phi2_s=phi2_s-360.0_pRe
CPFEM_ini_ori(1,i,CPFEM_en,CPFEM_in)=phi1_s
CPFEM_ini_ori(2,i,CPFEM_en,CPFEM_in)=PHI_s
CPFEM_ini_ori(3,i,CPFEM_en,CPFEM_in)=phi2_s
! *** Orientation matrix for initial orientation ***
call CPFEM_euldreh(phi1_s,PHI_s,phi2_s,ori_mat)
CPFEM_Fp_old(:,:,i,CPFEM_en,CPFEM_in)=ori_mat
i=i+1
! *** Third one ***
phi1_s=phi1+180.0_pRe
if (phi1_s.GT.360.0_pRe) phi1_s=phi1_s-360.0_pRe
PHI_s=PHI
phi2_s=phi2
CPFEM_ini_ori(1,i,CPFEM_en,CPFEM_in)=phi1_s
CPFEM_ini_ori(2,i,CPFEM_en,CPFEM_in)=PHI_s
CPFEM_ini_ori(3,i,CPFEM_en,CPFEM_in)=phi2_s
! *** Orientation matrix for initial orientation ***
call CPFEM_euldreh(phi1_s,PHI_s,phi2_s,ori_mat)
CPFEM_Fp_old(:,:,i,CPFEM_en,CPFEM_in)=ori_mat
i=i+1
else if ((CPFEM_mat(CPFEM_mn,2).NE.0).AND.
& (CPFEM_mat(CPFEM_mn,2).NE.2)) then
call CPFEM_error(520)
endif
enddo
CPFEM_tauc_slip_old(:,:,CPFEM_en,CPFEM_in)=s0_slip
endif
! ************************************
! *** CP-FEM Calculation ***
! ************************************
! *** Reinitialization of stress and consistent tangent ***
CPFEM_s=0
CPFEM_d=0
! *** Loop over all the components ***
do iori=1,nori
! *** Initialization of the matrices for t=t0 ***
Fp_old=CPFEM_Fp_old(:,:,iori,CPFEM_en,CPFEM_in)
tauc_slip_old=CPFEM_tauc_slip_old(:,iori,CPFEM_en,CPFEM_in)
tauc_slip_new=tauc_slip_old
g_old=CPFEM_g_old(:,iori,CPFEM_en,CPFEM_in)
Tstar_v=CPFEM_sigma_old(:,iori,CPFEM_en,CPFEM_in)
p10=CPFEM_ini_ori(1,iori,CPFEM_en,CPFEM_in)
P0=CPFEM_ini_ori(2,iori,CPFEM_en,CPFEM_in)
p20=CPFEM_ini_ori(3,iori,CPFEM_en,CPFEM_in)
vf=CPFEM_mat(CPFEM_mn,7*iori+2)
! *** Calculation of the solution at t=t1 ***
if (modulo(CPFEM_cn,ijaco).EQ.0) then
call CPFEM_stress(CPFEM_tinc,CPFEM_ffn,CPFEM_ffn1,Fp_old,Fp_new,
& g_old,g_new,tauc_slip_old,
& tauc_slip_new,
& Tstar_v,cs,cd,p11,P1,p21,dgmaxc,1,iconv,ising,
& icut,CPFEM_en,CPFEM_in,CPFEM_inc)
! *** Evaluation of ising ***
! *** ising=2 => singular matrix in jacobi calculation ***
! *** => use old jacobi ***
if (ising==2) then
jpara=0
endif
! *** Calculation of the consistent tangent ***
CPFEM_d=CPFEM_d+vf*cd
else
call CPFEM_stress(CPFEM_tinc,CPFEM_ffn,CPFEM_ffn1,Fp_old,Fp_new,
& g_old,g_new,tauc_slip_old,
& tauc_slip_new,
& Tstar_v,cs,hh6,p11,P1,p21,dgmaxc,0,iconv,
& ising,icut,CPFEM_en,CPFEM_in,CPFEM_inc)
jpara=0
endif
! *** Cases of unsuccessful calculations ***
! *** Evaluation od 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 CPFEM_error(700)
CPFEM_timefactor=1.e5_pRe
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 CPFEM_error(600)
CPFEM_timefactor=1.e5_pRe
return
endif
! *** Evaluation of iconv ***
! *** iconv!=0 => no convergence ***
if (iconv==1) then
write(6,*) 'Inner loop did not converged!'
write(6,*) 'Integration point: ',CPFEM_in
write(6,*) 'Element:',CPFEM_en
call CPFEM_error(600)
CPFEM_timefactor=1.e5_pRe
return
else
if (iconv==2) then
write(6,*) 'Outer loop did not converged!'
write(6,*) 'Integration point: ',CPFEM_in
write(6,*) 'Element: ',CPFEM_en
call CPFEM_error(600)
CPFEM_timefactor=1.e5_pRe
return
endif
endif
! *** Update the differents matrices for t=t1 ***
CPFEM_Fp_new(:,:,iori,CPFEM_en,CPFEM_in)=Fp_new
CPFEM_tauc_slip_new(:,iori,CPFEM_en,CPFEM_in)=tauc_slip_new
CPFEM_g_new(:,iori,CPFEM_en,CPFEM_in)=g_new
CPFEM_sigma_new(:,iori,CPFEM_en,CPFEM_in)=Tstar_v
! *** Calculation of the misorientation ***
phi1mis(1)=p10
PHImis(1)=P0
phi2mis(1)=p20
phi1mis(2)=p11
PHImis(2)=P1
phi2mis(2)=p21
call CPFEM_misori(phi1mis,PHImis,phi2mis,orimis)
! *** Update the results plotted in MENTAT ***
CPFEM_results(1,iori,CPFEM_en,CPFEM_in)=p11
CPFEM_results(2,iori,CPFEM_en,CPFEM_in)=P1
CPFEM_results(3,iori,CPFEM_en,CPFEM_in)=p21
CPFEM_results(4,iori,CPFEM_en,CPFEM_in)=orimis
CPFEM_results(5,iori,CPFEM_en,CPFEM_in)=sum(g_new)
CPFEM_results(7,iori,CPFEM_en,CPFEM_in)=sum(tauc_slip_new)/nslip
CPFEM_results(21,iori,CPFEM_en,CPFEM_in)=vf
! *** Evaluation of the maximum shear ***
dgmax=max(dgmax,dgmaxc)
! *** Evaluation of the average Cauchy stress ***
CPFEM_s=CPFEM_s+vf*cs
enddo
! *** End of the loop over the components ***
! *************************************
! *** End of the CP-FEM Calculation ***
! *************************************
! *** Approximate relative element thickness ***
call CPFEM_thick(CPFEM_ffn1,CPFEM_en,CPFEM_in)
! *** Restoration of the old jacobian if necessary ***
if (jpara==0) then
CPFEM_d=CPFEM_jaco_old(:,:,CPFEM_en,CPFEM_in)
else
! *** Store the new jacobian ***
CPFEM_jaco_old(:,:,CPFEM_en,CPFEM_in)=CPFEM_d
endif
! *** Calculate timefactor ***
CPFEM_timefactor=dgmax/dgs
return
end
subroutine CPFEM_stress(
&CPFEM_tinc,
&CPFEM_ffn,
&CPFEM_ffn1,
&Fp_old,
&Fp_new,
&g_old,
&g_new,
&tauc_slip_old,
&tauc_slip_new,
&Tstar_v,
&cs,
&dcs_de,
&phi1,
&PHI,
&phi2,
&dgmaxc,
&isjaco,
&iconv,
&ising,
&icut,
&CPFEM_en,
&CPFEM_in,
&CPFEM_inc
&)
c********************************************************************
c This routine calculates the stress for a single component
c and manages the independent time incrmentation
c********************************************************************
use mpie
use prec, only: pRe,pIn
implicit none
! *** Definition of variables ***
integer(pIn) isjaco,iconv,ising,icut,CPFEM_en,CPFEM_in,CPFEM_inc
real(pRe) CPFEM_tinc,CPFEM_ffn(3,3),CPFEM_ffn1(3,3),Fp_old(3,3),
& Fp_new(3,3),g_old(nslip),g_new(nslip),
& tauc_slip_old(nslip),tauc_slip_new(nslip),
& Tstar_v(6),
& cs(6),dcs_de(6,6),phi1,PHI,phi2,dgmaxc
integer(pIn) jcut
real(pRe) Tstar_v_h(6),tauc_slip_new_h(nslip),
& dt_i,delta_Fg(3,3),Fg_i(3,3),
& tauc_slip_new_i(nslip),time,mm(6,6)
! *** Numerical parameters ***
integer(pIn), parameter :: ncut=7_pIn
icut=0
! *** First attempt to calculate Tstar and tauc with initial timestep ***
Tstar_v_h=Tstar_v
tauc_slip_new_h=tauc_slip_new
call CPFEM_stress_int(CPFEM_tinc,CPFEM_ffn,CPFEM_ffn1,Fp_old,Fp_new,
& g_old,g_new,tauc_slip_old,
& tauc_slip_new,
& Tstar_v,cs,dcs_de,phi1,PHI,phi2,dgmaxc,
& isjaco,iconv,ising,CPFEM_en,CPFEM_in,CPFEM_inc)
if ((iconv==0).AND.(ising==0)) then
return
endif
! *** Calculation of stress and resistences with a cut timestep ***
! *** when first try did not converge ***
jcut=1_pIn
dt_i=0.5*CPFEM_tinc
delta_Fg=0.5*(CPFEM_ffn1-CPFEM_ffn)
Fg_i=CPFEM_ffn+delta_Fg
Tstar_v=Tstar_v_h
tauc_slip_new_i=tauc_slip_new_h
! *** Start time ***
time=dt_i
do while (time.LE.CPFEM_tinc)
call CPFEM_stress_int(time,CPFEM_ffn,Fg_i,Fp_old,Fp_new,g_old,
& g_new,tauc_slip_old,
& tauc_slip_new_i,
& Tstar_v,cs,mm,phi1,PHI,
& phi2,dgmaxc,0_pIn,iconv,ising,CPFEM_en,
& CPFEM_in,CPFEM_inc)
if ((iconv==0).AND.(ising==0)) then
time=time+dt_i
Fg_i=Fg_i+delta_Fg
Tstar_v_h=Tstar_v
tauc_slip_new_h=tauc_slip_new_i
else
jcut=jcut+1
if (jcut.GT.ncut) then
icut=1
return
endif
dt_i=0.5*dt_i
time=time-dt_i
delta_Fg=0.5*delta_Fg
Fg_i=Fg_i-delta_Fg
Tstar_v=Tstar_v_h
tauc_slip_new_i=tauc_slip_new_h
endif
enddo
! *** Final calculation of stress and resistences withb full timestep ***
tauc_slip_new=tauc_slip_new_i
call CPFEM_stress_int(CPFEM_tinc,CPFEM_ffn,CPFEM_ffn1,Fp_old,Fp_new,
& g_old,g_new,tauc_slip_old,
& tauc_slip_new,
& Tstar_v,cs,dcs_de,phi1,PHI,phi2,dgmaxc,
& isjaco,iconv,ising,CPFEM_en,CPFEM_in,CPFEM_inc)
return
end
subroutine CPFEM_stress_int(
&dt, ! Time increment
&Fg_old, ! Old global deformation gradient
&Fg_new, ! New global deformation gradient
&Fp_old, ! Old plastic deformation gradient
&Fp_new, ! New plastic deformation gradient
&g_old, ! Old cumulative plastic strain of a slip system
&g_new, ! New cumulative plastic strain of a slip system
&tauc_slip_old, ! Old resistence of a slip system
&tauc_slip_new, ! New resistence of a slip system
&Tstar_v, ! Second Piola-Kirschoff stress tensor
&cs, ! Cauchy stress vector
&dcs_de, ! Consistent tangent
&phi1, ! Euler angle phi1
&PHI, ! Euler angle PHI
&phi2, ! Euler angle phi2
&dgmaxc,
&isjaco,
&iconv,
&ising,
&CPFEM_en,
&CPFEM_in,
&CPFEM_inc
&)
c********************************************************************
c This routine calculates the stress for a single component
c it is based on the paper by Kalidindi et al.:
c J. Mech. Phys, Solids Vol. 40, No. 3, pp. 537-569, 1992
c it is modified to use anisotropic elasticity matrix
c********************************************************************
use mpie
use prec
implicit none
! *** Definition of variables ***
integer(pIn) isjaco,iconv,ising,CPFEM_en,CPFEM_in,CPFEM_inc
real(pRe) dt,Fg_old(3,3),Fg_new(3,3),Fp_old(3,3),Fp_new(3,3),
& g_old(nslip),g_new(nslip),
& tauc_slip_old(nslip),tauc_slip_new(nslip),
& Tstar_v(6),
& cs(6),dcs_de(6,6),phi1,PHI,phi2,dgmaxc
integer(pIn) ic
real(pRe) gdot_slip(nslip),Fe(3,3),R(3,3),
& U(3,3),de(3,3),tauc2(nslip),Fp2(3,3),
& sgm2(6),cs1(6),dF(3,3),Fg2(3,3),dev(6)
! *** Numerical parameters ***
real(pRe), parameter :: pert_ct=1.0e-5_pRe
! *** Error treatment ***
dgmaxc=0
iconv=0
ising=0
! *********************************************
! *** Calculation of the new Cauchy stress ***
! *********************************************
! *** Call Newton-Raphson method ***
call NEWTON_RAPHSON(dt,Fg_old,Fg_new,Fp_old,Fp_new,Fe,gdot_slip,
& tauc_slip_old,tauc_slip_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)
! *** Evaluation of the maximum slip shear ***
dgmaxc=maxval(abs(gdot_slip*dt))
g_new=g_old+abs(gdot_slip)*dt
! *** Choice of the calculation of the consistent tangent ***
if (isjaco==0) then
return
endif
! *********************************************
! *** 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
dev=0
if(ic.le.3) dev(ic)=pert_ct
if(ic.gt.3) dev(ic)=pert_ct/2
call CPFEM_conv6to33(dev,de)
dF=matmul(de,Fg_old)
Fg2=Fg_new+dF
sgm2=Tstar_v
tauc2=tauc_slip_new
! *** Calculation of the perturbated Cauchy stress ***
call NEWTON_RAPHSON(dt,Fg_old,Fg2,Fp_old,Fp2,Fe,gdot_slip,
& tauc_slip_old,tauc2,
& sgm2,cs1,iconv,ising)
! *** Consistent tangent ***
dcs_de(:,ic)=(cs1-cs)/pert_ct
enddo
return
end
subroutine NEWTON_RAPHSON(
&dt,
&Fg_old,
&Fg_new,
&Fp_old,
&Fp_new,
&Fe,
&gdot_slip,
&tauc_slip_old,
&tauc_slip_new,
&Tstar_v,
&cs,
&iconv,
&ising
&)
!***********************************************************************
!*** NEWTON-RAPHSON Calculation ***
!***********************************************************************
use mpie
use prec
implicit none
! *** Definition of variables ***
integer(pIn) isjaco,iconv,ising,CPFEM_en,CPFEM_in,CPFEM_inc
real(pRe) dt,Fg_old(3,3),Fg_new(3,3),Fp_old(3,3),Fp_new(3,3),
& g_old(nslip),g_new(nslip),
& tauc_slip_old(nslip),tauc_slip_new(nslip),
& Tstar_v(6),cs(6),dcs_de(6,6),phi1,PHI,phi2,dgmaxc
integer(pIn) i,j,k,iouter,iinner,ijac,ic
real(pRe) invFp_old(3,3),det,A(3,3),Estar0_v(6),Tstar0_v(6),
& mm(3,3),mm1(3,3),vv(6),Dslip(6,nslip),
& tau_slip(nslip),gdot_slip(nslip),
& R1(6),norm1,Tstar_v_per(6),R1_per(6),
& Jacobi(6,6),invJacobi(6,6),dTstar_v(6),R2(nslip),
& dtauc_slip(nslip),norm2,dLp(3,3),
& Estar(3,3),Estar_v(6),invFp_new(3,3),
& invFp2(3,3),Lp(3,3),Fe(3,3),
& R(3,3),U(3,3),dgdot_dtaucslip(nslip)
real(pRe) de(3,3),dev(6),tauc2(nslip),fp2(3,3),
& sgm2(6),cs1(6),df(3,3),
& fg2(3,3),tauc_old(nslip),crite,tol_in,tol_out
! *** Numerical parameters ***
integer(pIn), parameter :: nouter=50
real(pRe), parameter :: tol_outer=1.0e-4_pRe
integer(pIn), parameter :: ninner=2000
real(pRe), parameter :: tol_inner=1.0e-3_pRe
real(pRe), parameter :: eta=13.7_pRe
integer(pIn), parameter :: numerical=0
real(pRe), parameter :: pert_nr=1.0e-8_pRe
crite=eta*s0_slip/n_slip
! *** Tolerences ***
tol_in=tol_inner*s0_slip
tol_out=tol_outer*s0_slip
! *** Error treatment ***
dgmaxc=0
iconv=0
ising=0
! *** Calculation of Fp_old(-1) ***
invFp_old=Fp_old
call invert(invFp_old,3,0,0,det,3)
if (det==0.0_pRe) then
ising=1
return
endif
! *** Calculation of A and T*0 (see Kalidindi) ***
A=matmul(transpose(matmul(Fg_new,invFp_old)),
& matmul(Fg_new,invFp_old))
call CPFEM_conv33to6((A-I3)/2,Estar0_v)
Tstar0_v=matmul(Cslip_66,Estar0_v)
! *** Calculation of Dslip (see Kalidindi) ***
do i=1,nslip
mm=matmul(A,Sslip(i,:,:))
mm1=(mm+transpose(mm))/2
vv = math_33to6(mm1)
Dslip(:,i)=matmul(Cslip_66,vv)
enddo
! *** Second level of iterative procedure: Resistences ***
do iouter=1,nouter
! *** First level of iterative procedure: Stresses ***
do iinner=1,ninner
! *** Calculation of gdot_slip ***
do i=1,nslip
tau_slip(i)=dot_product(Tstar_v,Sslip_v(i,:))
enddo
call slip_rate(tau_slip,tauc_slip_new,gdot_slip,
& dgdot_dtaucslip)
! *** Evaluation of Tstar and Gn (see Kalidindi) ***
vv=0
do i=1,nslip
vv=vv-gdot_slip(i)*Dslip(:,i)
enddo
R1=Tstar_v-Tstar0_v-vv*dt
norm1=maxval(abs(R1))
if (norm1.LT.tol_in) then
goto 100
endif
! *** Jacobi Calculation ***
if (numerical==1) then
! *** Perturbation method ***
else
! *** Analytical Calculation ***
Jacobi=0
do i=1,nslip
do j=1,6
do k=1,6
Jacobi(j,k)=Jacobi(j,k)
& +Dslip(j,i)*Sslip_v(i,k)*dgdot_dtaucslip(i)
enddo
enddo
enddo
Jacobi=Jacobi*dt
do i=1,6
Jacobi(i,i)=1.0_pRe+Jacobi(i,i)
enddo
endif
! *** End of the Jacobi calculation ***
! *** Inversion of the Jacobi matrix ***
invJacobi=Jacobi
call invert(invJacobi,6,0,0,det,6)
if (det==0.0_pRe) then
do i=1,6
Jacobi(i,i)=1.05d0*maxval(Jacobi(i,:))
enddo
invJacobi=Jacobi
call invert(invJacobi,6,0,0,det,6)
if (det==0.0_pRe) then
ising=1
return
endif
endif
dTstar_v=matmul(invJacobi,R1)
! *** Correction (see Kalidindi) ***
do i=1,6
if (abs(dTstar_v(i)).GT.crite) then
dTstar_v(i)=sign(crite,dTstar_v(i))
endif
enddo
Tstar_v=Tstar_v-dTstar_v
enddo
iconv=1
return
! *** End of the first level of iterative procedure ***
100 continue
call hardening(tauc_slip_new,gdot_slip,dtauc_slip)
! *** Arrays of residuals ***
R2=tauc_slip_new-tauc_slip_old-dtauc_slip*dt
norm2=maxval(abs(R2))
if (norm2.LT.tol_out) then
goto 200
endif
tauc_slip_new=tauc_slip_old+dtauc_slip*dt
enddo
iconv=2
return
! *** End of the second level of iterative procedure ***
200 continue
call plastic_vel_grad(dt,tau_slip,tauc_slip_new,Lp)
! *** Calculation of Fp(t+dt) (see Kalidindi) ***
dLp=I3+Lp*dt
Fp_new=matmul(dLp,Fp_old)
call CPFEM_determ(Fp_new,det)
Fp_new=Fp_new/det**(1.0_pRe/3.0_pRe)
! *** Calculation of F*(t+dt) (see Kalidindi) ***
invFp_new=Fp_new
call invert(invFp_new,3,0,0,det,3)
if (det==0.0_pRe) then
ising=1
return
endif
Fe=matmul(Fg_new,invFp_new)
! *** Calculation of Estar ***
Estar=0.5_pRe*(matmul(transpose(Fe),Fe)-I3)
call CPFEM_conv33to6(Estar,Estar_v)
! *** Calculation of the Cauchy stress ***
call cauchy_stress(Estar_v,Fe,cs)
return
end
end module

254
IO.f90 Normal file
View File

@ -0,0 +1,254 @@
!##############################################################
MODULE IO
!##############################################################
CONTAINS
!---------------------------
! function IO_open_file(unit,relPath)
! function IO_open_inputFile(unit)
! function IO_stringPos(line,N)
! function IO_stringValue(line,positions,pos)
! function IO_floatValue(line,positions,pos)
! function IO_intValue(line,positions,pos)
! function IO_lowercase(line)
! subroutine IO_error(ID)
!---------------------------
!********************************************************************
! open existing file to given unit
! path to file is relative to working directory
!********************************************************************
logical FUNCTION IO_open_file(unit,relPath)
use prec, only: pInt
implicit none
character(len=*), parameter :: pathSep = achar(47)//achar(92) ! /, \
character(len=*) relPath
integer(pInt) unit
character(256) path
inquire(6, name=path) ! determine outputfile
open(unit,status='old',err=100,file=path(1:scan(path,pathSep,back=.true.))//relPath)
IO_open_file = .true.
return
100 IO_open_file = .false.
return
END FUNTION
!********************************************************************
! open FEM inputfile to given unit
!********************************************************************
logical FUNCTION IO_open_inputFile(unit)
use prec, only: pReal, pInt
implicit none
character(256) outName
integer(pInt) unit, extPos
character(3) ext
inquire(6, name=outName) ! determine outputfileName
extPos = len_trim(outName)-2
if(outName(extPos:extPos+2)=='out') then
ext='dat' ! MARC
else
ext='inp' ! ABAQUS
end if
open(unit,status='old',err=100,file=outName(1:extPos-1)//ext)
IO_open_inputFile = .true.
return
100 IO_open_inputFile = .false.
return
END FUNCTION
!********************************************************************
! locate at most N space-separated parts in line
! return array containing number of parts found and
! their left/right positions to be used by IO_xxxVal
!********************************************************************
FUNCTION IO_stringPos (line,N)
use prec, only: pReal,pInt
implicit none
character(len=*) line
character(len=*), parameter :: sep=achar(32)//achar(9) ! whitespaces
integer(pInt) N, part
integer(pInt) IO_stringPos(1+N*2)
IO_stringPos = -1
IO_stringPos(1) = 0
part = 1
do while ((N<1 .or. part<=N) .and. verify(line(IO_stringPos(part*2-1)+1:),sep)>0)
IO_stringPos(part*2) = IO_stringPos(part*2-1)+verify(line(IO_stringPos(part*2-1)+1:),sep)
IO_stringPos(part*2+1) = IO_stringPos(part*2)+scan(line(IO_stringPos(part*2):),sep)-2
part = part+1
end do
IO_stringPos(1) = part-1
return
END FUNCTION
!********************************************************************
! read string value at pos from line
!********************************************************************
FUNCTION IO_stringValue (line,positions,pos)
use prec, only: pReal,pInt
implicit none
character(len=*) line
integer(pInt) positions(*),pos
character(len=1+positions(pos*2+1)-positions(pos*2)) IO_stringValue
IO_stringValue = line(positions(pos*2):positions(pos*2+1))
return
END FUNCTION
!********************************************************************
! read float value at pos from line
!********************************************************************
FUNCTION IO_floatValue (line,positions,pos)
use prec, only: pReal,pInt
implicit none
character(len=*) line
real(pReal) IO_floatValue
integer(pInt) positions(*),pos
READ(UNIT=line(positions(pos*2):positions(pos*2+1)),ERR=100,FMT='(F)') IO_floatValue
return
100 IO_floatValue = -1.0_pReal
return
END FUNCTION
!********************************************************************
! read int value at pos from line
!********************************************************************
FUNCTION IO_intValue (line,positions,pos)
use prec, only: pReal,pInt
implicit none
character(len=*) line
integer(pInt) IO_intValue
integer(pInt) positions(*),pos
READ(UNIT=line(positions(pos*2):positions(pos*2+1)),ERR=100,FMT='(I)') IO_intValue
return
100 IO_intValue = -1_pInt
return
END FUNCTION
!********************************************************************
! change character in line to lower case
!********************************************************************
FUNCTION IO_lowercase (line)
use prec, only: pInt
implicit none
character (len=*) line
character (len=len(line)) IO_lowercase
integer(pInt) i
IO_lowercase = line
forall (i=1:len(line),64<ichar(line(i:i)).and.ichar(line(i:i))<91) IO_lowercase(i:i)=achar(ichar(line(i:i))+32)
return
END FUNCTION
!********************************************************************
! in place change character in line to lower case
!********************************************************************
SUBROUTINE IO_lowercaseInplace (line)
use prec, only: pInt
implicit none
character (len=*) line
integer(pInt) i
forall (i=1:len(line),64<ichar(line(i:i)).and.ichar(line(i:i))<91) line(i:i)=achar(ichar(line(i:i))+32)
return
END SUBROUTINE
!********************************************************************
! write error statements to standard out
! and terminate the Marc run with exit #9xxx
! in ABAQUS either time step is reduced or execution terminated
!********************************************************************
SUBROUTINE IO_error(ID)
use prec, only: pInt
implicit none
integer(pInt) ID
character(len=80) msg
select case (ID)
case (100)
msg='File material.mpie can not be opened'
case (110)
msg='File material99.mpie can not be opened'
case (120)
msg='File with c-coefficience can not be opened'
case (130)
msg='File with single orientations can not be opened'
case (200)
msg='Error reading from file material.mpie'
case (210)
msg='Error reading from file material99.mpie'
case (220)
msg='Error reading from file containing c-coefficiences'
case (230)
msg='Error reading from file containing single orientations'
case (300)
msg='This material can only be used with &
&elements with three direct stress components'
case (400)
msg='Unknown alloy number specified'
case (500)
msg='Unknown lattice number specified'
case (510)
msg='Unknown component type specified'
case (520)
msg='Unknown component symmetry specified'
case (600)
msg='Stress iteration did not converge'
case (700)
msg='Singular matrix in stress iteration'
case default
msg='Unknown error number'
end select
write(6,*) 'MPIE Material Routine Ver. 0.7 by Dr. F. Roters'
write(6,*)
write(6,*) msg
call flush(6)
! call quit(9000+ID)
! ABAQUS returns in some cases
return
END SUBROUTINE
END MODULE IO

423
constitutive.f90 Normal file
View File

@ -0,0 +1,423 @@
! ---------------------------
MODULE constitutive
! ---------------------------
! *** constitutive equations ***
use prec, only: pRe,pIn
implicit none
! ***************************
! *** Material parameters ***
! ***************************
real(pRe), allocatable :: Cslip_66(:,:,:),Cslip_3333(:,:,:,:,:)
real(pRe), allocatable :: s0_slip(:),gdot0_slip(:)
real(pRe), allocatable :: h0(:),w0(:),s_sat(:),q0(:),n_slip(:)
real(pRe), allocatable :: hardening_matrix(:,:,:)
character*80, allocatable :: TCfile(:), ODFfile(:)
real(pRe), parameter :: latent=1.4_pRe
integer(pIn), parameter :: Nslip(3)
integer(pIn) Nmats
real(pRe) sn(3,48,3),sd(3,48,3)
real(pRe) Sslip(3,48,3,3),Sslip_v(3,48,6)
! *** Vectors n and d for each fcc slip systems ***
! MISSING needs to be generalized to fcc and bcc (and hcp?)
! 1: fcc, 2: bcc, 3: hcp
! the respective crystal structure has to be defined
! via material parameter 'crystal_structure' in [material]
data sd( 1,:)/ 0, 1,-1/ ; data sn( 1,:)/ 1, 1, 1/
data sd( 2,:)/-1, 0, 1/ ; data sn( 2,:)/ 1, 1, 1/
data sd( 3,:)/ 1,-1, 0/ ; data sn( 3,:)/ 1, 1, 1/
data sd( 4,:)/ 0,-1,-1/ ; data sn( 4,:)/-1,-1, 1/
data sd( 5,:)/ 1, 0, 1/ ; data sn( 5,:)/-1,-1, 1/
data sd( 6,:)/-1, 1, 0/ ; data sn( 6,:)/-1,-1, 1/
data sd( 7,:)/ 0,-1, 1/ ; data sn( 7,:)/ 1,-1,-1/
data sd( 8,:)/-1, 0,-1/ ; data sn( 8,:)/ 1,-1,-1/
data sd( 9,:)/ 1, 1, 0/ ; data sn( 9,:)/ 1,-1,-1/
data sd(10,:)/ 0, 1, 1/ ; data sn(10,:)/-1, 1,-1/
data sd(11,:)/ 1, 0,-1/ ; data sn(11,:)/-1, 1,-1/
data sd(12,:)/-1,-1, 0/ ; data sn(12,:)/-1, 1,-1/
contains
! **************************************
! *** module Init ***
! **************************************
subroutine constitutive_init()
call constitutive_calc_SchmidM()
call constitutive_calc_hardeningM()
call constitutive_parse_materialDat()
end subroutine
! **************************************
! *** Calculation of Schmid matrices ***
! **************************************
subroutine constitutive_calc_SchmidM()
use prec, only: pRe,pIn
implicit none
integer(pIn) i,j,k,l
real(pRe) invNorm
do j=1,3 ! iterate over crystal system
do i=1,Nslip(j) ! iterate over slip systems
do k=1,3
do l=1,3
Sslip(j,i,k,l)=sd(j,i,k)*sn(j,i,l)
enddo
enddo
invNorm = dsqrt(1.0_pRe/
& (sn(j,i,1)**2+sn(j,1,2)**2+sn(j,i,3)**2)/
& (sd(j,i,1)**2+sd(j,1,2)**2+sd(j,i,3)**2))
Sslip(j,i,:,:) = Sslip(j,i,:,:)*invNorm
Sslip_v(j,i,1)=Sslip(j,i,1,1)
Sslip_v(j,i,2)=Sslip(j,i,2,2)
Sslip_v(j,i,3)=Sslip(j,i,3,3)
Sslip_v(j,i,4)=Sslip(j,i,1,2)+Sslip(j,i,2,1)
Sslip_v(j,i,5)=Sslip(j,i,2,3)+Sslip(j,i,3,2)
Sslip_v(j,i,6)=Sslip(j,i,1,3)+Sslip(j,i,3,1)
enddo
enddo
end subroutine
! ****************************************
! *** Hardening matrix (see Kalidindi) ***
! ****************************************
subroutine constitutive_calc_hardeningM()
use prec, only: pRe,pIn
implicit none
integer(pIn) i,j,k,l
! MISSING iteration over crystal systems
! PE does not understand the j,k looping
hardening_matrix=latent
do i=1,10,3
do j=1,3
do k=1,3
hardening_matrix(i-1+j,i-1+k)=1.0_ZdRe
enddo
enddo
enddo
! ****************************************
! *** Reading 'material.mpie' ***
! ****************************************
subroutine constitutive_parse_materialDat()
use prec, only: pRe,pIn
implicit none
character*80 line
integer(pIn) i,j,k,l,positions(4)
! MISSING: needs to be 2 pass
! first pass to count Nmats and allocate
! 2nd pass to read actual parameters
write(6,*) '## constitutive_parse_materialDat ##'
write(6,*)
constitutive_Nmats = 1
open(200,FILE='material.mpie',ACTION='READ',STATUS='OLD',ERR=100)
read(200,610,ERR=200,END=200) line
IF( line(1:1).ne.'[' )THEN
WRITE(6,*) 'Problem with mat file: no mat. in 1st line'
ELSE
WRITE(6,*) 'Reading mat. data'
DO WHILE( .true. )
READ(200,610,END=220) line
IF( line(1:1).eq.'[' )THEN
constitutive_Nmats = constitutive_Nmats+1
ELSE
positions = IO_stringPos(line,2) ! parse 2 parts
SELECT CASE (IO_stringValue(line,positions,1))
CASE ('s0_slip')
s0_slip(mat) = IO_floatValue(line,positions,2)
CASE ('g0_slip')
g0_slip(mat) = IO_floatValue(line,positions,2)
CASE ('n_slip')
n_slip(mat) = IO_intValue(line,positions,2)
CASE ('h0')
h0(mat) = IO_floatValue(line,positions,2)
CASE ('w0')
w0(mat) = IO_floatValue(line,positions,2)
CASE ('tauc_sat')
tauc_sat(mat) = IO_floatValue(line,positions,2)
CASE ('C11')
C11(mat) = IO_floatValue(line,positions,2)
CASE ('C12')
C12(mat) = IO_floatValue(line,positions,2)
CASE ('C44')
C44(mat) = IO_floatValue(line,positions,2)
CASE ('TCfile')
TCfile(mat) = IO_stringValue(line,positions,2)
CASE ('ODFfile')
ODFfile(mat) = IO_stringValue(line,positions,2)
CASE ('Ngrains')
Ngrains(mat) = IO_intValue(line,positions,2)
CASE DEFAULT
WRITE(6,*) 'Unknown mat. parameter ',line
END IF
END DO
END IF
220 continue
close(200)
! ** Defintion of stiffness matrices **
! MISSING: this needs to be iterated over the materials
Cslip_66 = 0.0_pRe
do i=1,3
do j=1,3
Cslip_66(i,j) = C12
enddo
Cslip_66(i,i) = C11
Cslip_66(i+3,i+3) = C44
enddo
Cslip_3333(:,:,:,:) = math_66to3333(Cslip_66(:,:))
! *** Transformation to get the MARC order ***
! *** 11,22,33,12,23,13 ***
! MISSING this should be outsourced to FEM-spec
temp=Cslip_66(4,:)
Cslip_66(4,:)=Cslip_66(6,:)
Cslip_66(6,:)=Cslip_66(5,:)
Cslip_66(5,:)=temp
temp=Cslip_66(:,4)
Cslip_66(:,4)=2.0d0*Cslip_66(:,6)
Cslip_66(:,6)=2.0d0*Cslip_66(:,5)
Cslip_66(:,5)=2.0d0*temp
! *** Output to MARC output file ***
write(6,*) 'Material data:'
write(6,*) 'Slip parameter:(s0_slip,g0_slip,n_slip)'
write(6,*) s0_slip,g0_slip,n_slip
write(6,*) 'Slip hardening parameter:(h0,tauc_sat,w0)'
write(6,*) h0,tauc_sat,w0
write(6,*) 'Elasticity matrix:'
write(6,*) Cslip_66(1,:)
write(6,*) Cslip_66(2,:)
write(6,*) Cslip_66(3,:)
write(6,*) Cslip_66(4,:)/2.0d0
write(6,*) Cslip_66(5,:)/2.0d0
write(6,*) Cslip_66(6,:)/2.0d0
write(6,*)
call flush(6)
! END OF MISSING mat iterations
return
100 call _error(110)
200 call _error(210)
end
subroutine READ_ORIENTATIONS
!***********************************************************************
!*** This routine reads orientations from 'orientations.mpie' ***
!***********************************************************************
use mpie
use Zahlendarstellung, only: ZdRe,ZdIn
implicit none
! *** Definition of variables ***
integer(ZdIn) i,j
! *** Read 'orientations.mpie' file ***
open(100,FILE='orientations.mpie',ACTION='READ',STATUS='OLD',
& ERR=100)
read(100,*,ERR=200,END=200)
! *** Read number of states, maximum of components over the states ***
read(100,*,ERR=200,END=200) mpie_nmat,mpie_norimx
! *** Allocate memory for the arrays ***
allocate(mpie_mat(mpie_nmat,2+7*mpie_norimx))
allocate(mpie_cko(mpie_nmat,4:35,3,0:35,2))
allocate(mpie_ckofile(mpie_nmat,80))
allocate(mpie_odfmax(mpie_nmat))
mpie_mat=0.0_ZdRe
mpie_cko=0.0_ZdRe
mpie_ckofile=''
mpie_odfmax=0.0_ZdRe
! *** Read the different states ***
do i=1,mpie_nmat
read(100,*,ERR=200,END=200)
! *** Number of component and symmetry ***
read(100,*,ERR=200,END=200) mpie_mat(i,1),mpie_mat(i,2)
! *** If symmetry = 2, use direct ODF sampling,i.e. read coefficience ***
if (mpie_mat(i,2)==2_ZdIn) then
read(100,'(80A)',ERR=200,END=201) mpie_ckofile(i,:)
201 call mpie_read_ckofile(mpie_cko(i,:,:,:,:),
& mpie_ckofile(i,:))
call mpie_odf_max(mpie_cko(i,:,:,:,:),mpie_odfmax(i))
! *** Set volume fraction to inverse of orientation number for each orientation ***
do j=1,int(mpie_mat(i,1),ZdIn)
mpie_mat(i,2+7*j)=1/mpie_mat(i,1)
enddo
else
! *** Read for every component: ***
! *** gauss: euler angles (phi1, PHI, phi2), dummy, scatter, volume fraction ***
! *** fiber: alpha1, alpha2, beta1, beta2, scatter, volume fraction ***
do j=1,int(mpie_mat(i,1),ZdIn)
read(100,*,ERR=200,END=200) mpie_mat(i,7*j-4),
& mpie_mat(i,7*j-3),mpie_mat(i,7*j-2),
& mpie_mat(i,7*j-1),mpie_mat(i,7*j),
& mpie_mat(i,7*j+1),mpie_mat(i,7*j+2)
enddo
endif
enddo
close(100)
! *** Output to MARC output file ***
write(6,*) 'MPIE Material Routine Ver. 0.1 by L. Hantcherli'
write(6,*)
write(6,*) 'Orientations data:'
write(6,*) 'Number of materials: ', mpie_nmat
write(6,*) 'Maximum number of components: ', mpie_norimx
write(6,*)
do i=1,mpie_nmat
write(6,*) 'State', i
if (mpie_mat(i,2)==2_ZdIn) then
write(6,*) mpie_ckofile(i,:),mpie_mat(i,9),mpie_odfmax(i)
else
write(6,*) mpie_mat(i,:)
endif
write(6,*)
enddo
call flush(6)
return
100 call _error(100)
200 call _error(200)
end
subroutine slip_rate (tau_slip,tauc_slip_new,gdot_slip,
& dgdot_dtaucslip)
C ********************************************************************
C Subroutine contains the constitutive equation for the slip
C rate on each slip system
C Input: tau_slip : shear stress on each slip system
C tauc_slip_new : critical shear stress on each slip system
C Output: gdot_slip : slip rate on each slip system
C dgdot_dtaucslip: derivative of slip rate on each slip system
C ********************************************************************
use mpie
use Zahlendarstellung
implicit none
real(ZdRe) tau_slip(nslip),tauc_slip_new(nslip),
& gdot_slip(nslip),dgdot_dtaucslip(nslip)
integer(ZdIn) i
do i=1,nslip
gdot_slip(i)=g0_slip*(abs(tau_slip(i))/tauc_slip_new(i))
& **n_slip*sign(1.0_ZdRe,tau_slip(i))
dgdot_dtaucslip(i)=g0_slip*(abs(tau_slip(i))/tauc_slip_new(i))
& **(n_slip-1) *n_slip/tauc_slip_new(i)
enddo
return
end
subroutine hardening (tauc_slip_new,gdot_slip,dtauc_slip)
C *********************************************************************
C Subroutine calculates the increment in critical shear stress due
C to plastic deformation on each slip system
C Input: tauc_slip_new :critical shear stress needed for slip on each
C slip system
C gdot_slip :slip rate on each slip system
C Output: dtauc_slip :increment of hardening due to slip on each
C slip system
C Local : selfhr
C *********************************************************************
use mpie
use Zahlendarstellung
implicit none
real(ZdRe) tauc_slip_new(nslip),gdot_slip(nslip),
& dtauc_slip(nslip)
real(ZdRe) selfhr(nslip)
integer(ZdIn) i
do i=1,nslip
selfhr(i)=h0*(1.0_ZdRe-tauc_slip_new(i)/
& tauc_sat)**w0
& *abs(gdot_slip(i))
enddo
dtauc_slip=matmul(hardening_matrix,selfhr)
return
end
subroutine plastic_vel_grad(dt,tau_slip,tauc_slip_new,Lp)
C *************************************************************
C Subroutine calculates the plastic velocity gradient given the
C slip rates
C Input: dt : time step
C tau_slip : shear stress on each slip system on each
C slip system
C tauc_slip_new : critical shear stress needed for slip on each
C slip system
C Output: Lp : plastic velocity gradient
C gdot_slip : slip rate on each slip system
C *************************************************************
use mpie
use Zahlendarstellung
implicit none
real(ZdRe) dt,tau_slip(nslip),tauc_slip_new(nslip),
& Lp(3,3),gdot_slip(nslip)
integer(ZdIn) i
Lp=0
do i=1,nslip
gdot_slip(i)=g0_slip*(abs(tau_slip(i))/tauc_slip_new(i))
& **n_slip*sign(1.0_ZdRe,tau_slip(i))
Lp=Lp+gdot_slip(i)*Sslip(i,:,:)
enddo
return
end
function CPFEM_Cauchy(Estar_v,Fe,C66)
C ***************************************************************
C Subroutine calculates the cauchy from the elastic strain tensor
C Input: Estar_v : elastic strain tensor (in vector form)
C Fe : elastic deformation gradient
C C66 : Stiffness Tensor
C Output: cs : cauchy stress
C Local: Tstar_v,Tstar,mm,det
C ***************************************************************
use math
use prec
implicit none
real(pRe) Estar_v(6),Fe(3,3),C66(6,6),CPFEM_Cauchy(6)
real(pRe) det,mm(3,3),Tstar(3,3)
integer(pIn) i
det = math_det(Fe)
Tstar = math_6to33(matmul(C66,Estar_v))
mm=matmul(matmul(Fe,Tstar),transpose(Fe))/det
CPFEM_Cauchy = math_33to6(mm)
return
end function
end module

1460
math.f90 Normal file
View File

File diff suppressed because it is too large Load Diff

210
mesh.f90 Normal file
View File

@ -0,0 +1,210 @@
!##############################################################
MODULE mesh
!##############################################################
use prec, only: pRe,pIn
implicit none
! ---------------------------
! _Nelems : total number of elements in mesh
! _Nnodes : total number of nodes in mesh
! _maxNnodes : max number of nodes in any element
! _maxNips : max number of IPs in any element
! _element : type, material, node indices
! _node : x,y,z coordinates (initially!)
! _nodeIndex : count of elements containing node,
! [element_num, node_index], ...
! _envIP : 6 neighboring IPs as [element_num, IP_index]
! order is +x, +y,+z, -x, -y, -z in local coord
! ---------------------------
integer(pIn) mesh_Nelems, mesh_Nnodes, mesh_maxNnodes,mesh_maxNips
integer(pIn), allocatable :: mesh_element (:,:)
integer(pIn), allocatable :: mesh_nodeIndex (:,:)
integer(pIn), allocatable :: mesh_envIP (:,:)
real(pRe), allocatable :: mesh_node (:,:)
CONTAINS
! ---------------------------
! subroutine mesh_init()
! subroutine mesh_parse_inputFile()
! ---------------------------
! ***********************************************************
! initialization
! ***********************************************************
SUBROUTINE mesh_init ()
mesh_Nelems = 0_pIn
mesh_Nnodes = 0_pIn
mesh_maxNips = 0_pIn
mesh_maxNnodes = 0_pIn
call mesh_parse_inputFile ()
END SUBROUTINE
! ***********************************************************
! parsing of input file
! ***********************************************************
FUNCTION mesh_parse_inputFile ()
use prec, only: pRe,pIn
use IO
implicit none
logical mesh_parse_inputFile
integer(pIn) i,j,positions(10*2+1)
integer(pIn) elem_num,elem_type,Nnodes,node_num,num_ip,mat,tp(70,2)
! Set a format to read the entire line (max. len is 80 characters)
610 FORMAT(A80)
if (.not. IO_open_inputFile(600)) then
mesh_parse_inputFile = .false.
return
endif
do while(.true.)
read(600,610,end=620) line
positions = IO_stringPos(line,3)
select case (IO_stringValue(line,positions,1)
!-----------------------------------
case ('sizing')
!-----------------------------------
mesh_Nelems = IO_intValue(line,positions,2)
mesh_Nnodes = IO_intValue(line,positions,3)
!-----------------------------------
case ('elements')
!-----------------------------------
select case (IO_intValue(line,positions,2)) ! elem type
case (3) ! 2D Triangle
mesh_maxNips = max(3,mesh_maxNips)
mesh_maxNnodes = max(3,mesh_maxNnodes)
case (6) ! 2D Quad.
mesh_maxNips = max(4,mesh_maxNips)
mesh_maxNnodes = max(4,mesh_maxNnodes)
case (7) ! 3D hexahedral
mesh_maxNips = max(8,mesh_maxNips)
mesh_maxNnodes = max(8,mesh_maxNnodes)
case default
mesh_maxNips = max(8,mesh_maxNips)
mesh_maxNnodes = max(8,mesh_maxNnodes)
end select
!-----------------------------------
case ('connectivity')
!-----------------------------------
allocate (mesh_element(mesh_Nelems,2+mesh_maxNips))
allocate (mesh_nodeIndex (mesh_Nnodes,1+mesh_maxNnodes*2)
allocate (mesh_envIP(mesh_Nelems,mesh_maxNips,6,2))
mesh_element = 0_pIn
mesh_nodeIndex = 0_pIn
mesh_envIP = 0_pIn
! MISSING: setting up of envIP
read(600,610,end=620) line ! skip line ??
do i=1,mesh_Nelems
read(600,610,end=620) line
positions = IO_stringPos(line,0) ! find all chunks
elem_num = IO_intValue(line,positions,1)
elem_type = IO_intValue(line,positions,2)
select case (elem_type)
case (3) ! 2D Triangle
Nnodes = 3
case (6) ! 2D Quad.
Nnodes = 4
case (7) ! 3D hexahedral
Nnodes = 8
case default
Nnodes = 8
end select
mesh_element(elem_num,1) = elem_type
do j=1,Nnodes ! store all node indices
node_num = IO_intValue(line,positions,2+j)
mesh_element(elem_num,1+j) = node_num
mesh_nodeIndex(node_num,1) = mesh_nodeIndex(node_num,1)+1 ! inc count
mesh_nodeIndex(node_num,mesh_nodeIndex(node_num,1)*2 ) = elem_num
mesh_nodeIndex(node_num,mesh_nodeIndex(node_num,1)*2+1) = j
end do
end do
!-----------------------------------
case ('coordinates')
!-----------------------------------
allocate (mesh_node(mesh_Nnodes,3)) ! x,y,z, per node
read(600,610,end=620) line ! skip line ??
do i=1,mesh_Nnodes
read(600,610,end=620) line
positions = IO_stringPos(line,0) ! find all (4) chunks
node_num = IO_intValue(line,positions,1)
do j=1,3 ! store x,y,z coordinates
mesh_node(node_num,j) = IO_floatValue(line,positions,1+j)
end do
end do
!-----------------------------------
case ('hypoelastic')
!-----------------------------------
! ***************************************************
! Search for key word "hypoelastic".
! This section contains the # of materials and
! the element range of each material
! ***************************************************
ELSE IF( line(1:11).eq.'hypoelastic' )THEN
mat=0
flag=0
DO WHILE( line(1:8).ne.'geometry' )
READ(600,610,END=620) line
i=1
DO WHILE( i.le.len(line)-8 )
IF( line(i:i+2).eq.'mat' )THEN
mat=mat+1
flag=1
END IF
i=i+1
END DO
IF( flag.eq.1 )THEN
flag=0
READ(600,610,END=620) line
READ(600,610,END=620) line
i=1
DO WHILE( line(i:i).eq.' ')
i=i+1
END DO
left=i
DO WHILE( line(i:i).ne.' ')
i=i+1
END DO
right=i-1
READ(UNIT=line(left:right), FMT=' (I5) ') tp(mat,1)
DO WHILE( (line(i:i).eq.' ').or.
& (line(i:i).eq.'t').or.
& (line(i:i).eq.'o') )
i=i+1
END DO
left=i
DO WHILE( line(i:i).ne.' ')
i=i+1
END DO
right=i-1
READ(UNIT=line(left:right), FMT=' (I5) ') tp(mat,2)
END IF
END DO
WRITE(6,*) 'mat: ',mat,' ',tp(1,1),' ',tp(1,2)
end select
END DO
! Code jumps to 620 when it reaches the end of the file
620 continue
WRITE(6,*) 'Finished with .dat file'
CALL FLUSH(6)
END FUNCTION
END MODULE mesh

11
prec.f90 Normal file
View File

@ -0,0 +1,11 @@
!##############################################################
MODULE prec
!##############################################################
implicit none
integer, parameter :: pReal = 8
integer, parameter :: pInz = 4
END MODULE prec