From fe70a82d6dc860b4fa9812d31effbefa78342ec3 Mon Sep 17 00:00:00 2001 From: Luc Hantcherli Date: Fri, 14 Dec 2007 13:36:04 +0000 Subject: [PATCH] !!!!! IMPORTANT !!!!! All subroutines were committed at once: - constitutive_pheno works - constitutive_dislo without twinning also works This release should serve as reference --- trunk/CPFEM.f90 | 38 ++++++++------------------ trunk/IO.f90 | 4 +-- trunk/constitutive_dislo.f90 | 35 +++++++++--------------- trunk/constitutive_pheno.f90 | 6 ++--- trunk/crystal.f90 | 22 ++++++++------- trunk/math.f90 | 2 +- trunk/mattex.mpie | 48 +++++++++++++++++---------------- trunk/mesh.f90 | 2 +- trunk/mpie_cpfem_marc2005r3.f90 | 2 +- trunk/mpie_cpfem_marc2007r1.f90 | 2 +- trunk/prec.f90 | 2 +- 11 files changed, 72 insertions(+), 91 deletions(-) diff --git a/trunk/CPFEM.f90 b/trunk/CPFEM.f90 index 63b04fd14..33ef9deee 100644 --- a/trunk/CPFEM.f90 +++ b/trunk/CPFEM.f90 @@ -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) diff --git a/trunk/IO.f90 b/trunk/IO.f90 index baba4a9a6..95a7452bb 100644 --- a/trunk/IO.f90 +++ b/trunk/IO.f90 @@ -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 diff --git a/trunk/constitutive_dislo.f90 b/trunk/constitutive_dislo.f90 index 82cd87fc8..823355e6f 100644 --- a/trunk/constitutive_dislo.f90 +++ b/trunk/constitutive_dislo.f90 @@ -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)) -[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 [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 diff --git a/trunk/mesh.f90 b/trunk/mesh.f90 index 84ed3035e..ccbce5503 100644 --- a/trunk/mesh.f90 +++ b/trunk/mesh.f90 @@ -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,*) diff --git a/trunk/mpie_cpfem_marc2005r3.f90 b/trunk/mpie_cpfem_marc2005r3.f90 index c712ba784..cfc1761db 100644 --- a/trunk/mpie_cpfem_marc2005r3.f90 +++ b/trunk/mpie_cpfem_marc2005r3.f90 @@ -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,& diff --git a/trunk/mpie_cpfem_marc2007r1.f90 b/trunk/mpie_cpfem_marc2007r1.f90 index 489693f94..987ef360e 100644 --- a/trunk/mpie_cpfem_marc2007r1.f90 +++ b/trunk/mpie_cpfem_marc2007r1.f90 @@ -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) diff --git a/trunk/prec.f90 b/trunk/prec.f90 index 7e3ca1211..4e781356e 100644 --- a/trunk/prec.f90 +++ b/trunk/prec.f90 @@ -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 ***