From 796fe0e8e92c16ffad7549ca836e17c2a1a2779b Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Tue, 20 Mar 2007 13:38:35 +0000 Subject: [PATCH] initial import --- CPFEM.f90 | 883 ++++++++++++++++++++++++++++ IO.f90 | 254 ++++++++ constitutive.f90 | 423 ++++++++++++++ math.f90 | 1460 ++++++++++++++++++++++++++++++++++++++++++++++ mesh.f90 | 210 +++++++ prec.f90 | 11 + 6 files changed, 3241 insertions(+) create mode 100644 CPFEM.f90 create mode 100644 IO.f90 create mode 100644 constitutive.f90 create mode 100644 math.f90 create mode 100644 mesh.f90 create mode 100644 prec.f90 diff --git a/CPFEM.f90 b/CPFEM.f90 new file mode 100644 index 000000000..8d2d3df9f --- /dev/null +++ b/CPFEM.f90 @@ -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 \ No newline at end of file diff --git a/IO.f90 b/IO.f90 new file mode 100644 index 000000000..72acd1462 --- /dev/null +++ b/IO.f90 @@ -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