From e562df35a984a8b040b91fcea602c47fa090bba8 Mon Sep 17 00:00:00 2001 From: Alankar Alankar Date: Wed, 29 Sep 2010 06:35:08 +0000 Subject: [PATCH] Added more state dependent variables in post results. --- code/constitutive_titanmod.f90 | 176 ++++++++++++++++++++++----------- 1 file changed, 117 insertions(+), 59 deletions(-) diff --git a/code/constitutive_titanmod.f90 b/code/constitutive_titanmod.f90 index d7d945720..193ddc427 100644 --- a/code/constitutive_titanmod.f90 +++ b/code/constitutive_titanmod.f90 @@ -12,14 +12,18 @@ character(len=18), dimension(2), parameter:: constitutive_titanmod_listBasicSlip 'rho_screw'/) character(len=18), dimension(1), parameter:: constitutive_titanmod_listBasicTwinStates = (/'gdot_twin'/) -character(len=18), dimension(8), parameter:: constitutive_titanmod_listDependentSlipStates =(/'segment_edge', & +character(len=18), dimension(12), parameter:: constitutive_titanmod_listDependentSlipStates =(/'segment_edge', & 'segment_screw', & 'resistance_edge', & 'resistance_screw', & 'tau_slip', & 'gdot_slip', & 'velocity_edge', & - 'velocity_screw' & + 'velocity_screw', & + 'gdot_slip_edge', & + 'gdot_slip_screw', & + 'stressratio_edge_p', & + 'stressratio_screw_p' & /) character(len=18), dimension(2), parameter:: constitutive_titanmod_listDependentTwinStates =(/'twin_fraction', & @@ -49,7 +53,8 @@ real(pReal), dimension(:), allocatable :: constitutive_titanmod_ constitutive_titanmod_C13, & ! C13 element in elasticity matrix constitutive_titanmod_C33, & ! C33 element in elasticity matrix constitutive_titanmod_C44, & ! C44 element in elasticity matrix - constitutive_titanmod_debyfrequency, & !Deby frequency + constitutive_titanmod_debyefrequency, & !Debye frequency + constitutive_titanmod_kinkf0, & !Debye frequency constitutive_titanmod_Gmod, & ! shear modulus constitutive_titanmod_CAtomicVolume, & ! atomic volume in Bugers vector unit constitutive_titanmod_dc, & ! prefactor for self-diffusion coefficient @@ -213,7 +218,8 @@ allocate(constitutive_titanmod_C12(maxNinstance)) allocate(constitutive_titanmod_C13(maxNinstance)) allocate(constitutive_titanmod_C33(maxNinstance)) allocate(constitutive_titanmod_C44(maxNinstance)) -allocate(constitutive_titanmod_debyfrequency(maxNinstance)) +allocate(constitutive_titanmod_debyefrequency(maxNinstance)) +allocate(constitutive_titanmod_kinkf0(maxNinstance)) allocate(constitutive_titanmod_Gmod(maxNinstance)) allocate(constitutive_titanmod_CAtomicVolume(maxNinstance)) allocate(constitutive_titanmod_dc(maxNinstance)) @@ -233,7 +239,8 @@ constitutive_titanmod_C12 = 0.0_pReal constitutive_titanmod_C13 = 0.0_pReal constitutive_titanmod_C33 = 0.0_pReal constitutive_titanmod_C44 = 0.0_pReal -constitutive_titanmod_debyfrequency = 0.0_pReal +constitutive_titanmod_debyefrequency = 0.0_pReal +constitutive_titanmod_kinkf0 = 0.0_pReal constitutive_titanmod_Gmod = 0.0_pReal constitutive_titanmod_CAtomicVolume = 0.0_pReal constitutive_titanmod_dc = 0.0_pReal @@ -367,9 +374,12 @@ do ! read thru sections of case ('c44') constitutive_titanmod_C44(i) = IO_floatValue(line,positions,2) write(6,*) tag,constitutive_titanmod_C44(i) - case ('debyfrequency') - constitutive_titanmod_debyfrequency(i) = IO_floatValue(line,positions,2) - write(6,*) tag,constitutive_titanmod_debyfrequency(i) + case ('debyefrequency') + constitutive_titanmod_debyefrequency(i) = IO_floatValue(line,positions,2) + write(6,*) tag,constitutive_titanmod_debyefrequency(i) + case ('kinkf0') + constitutive_titanmod_kinkf0(i) = IO_floatValue(line,positions,2) + write(6,*) tag,constitutive_titanmod_kinkf0(i) case ('nslip') forall (j = 1:lattice_maxNslipFamily) & constitutive_titanmod_Nslip(j,i) = IO_intValue(line,positions,1+j) @@ -697,7 +707,21 @@ do i = 1,maxNinstance constitutive_titanmod_sizeDotState(i)+ & size(constitutive_titanmod_listDependentSlipStates)*ns+size(constitutive_titanmod_listDependentTwinStates)*nt write(6,*) 'Determined size of state and dot state' - !* Determine size of postResults array + !* Determine size of postResults array + +! 'segment_edge' +! 'segment_screw' +! 'resistance_edge' +! 'resistance_screw' +! 'tau_slip' +! 'gdot_slip' +! 'velocity_edge' +! 'velocity_screw' +! 'gdot_slip_edge' +! 'gdot_slip_screw' +! 'StressRatio_edge_p' +! 'StressRatio_screw_p' + do o = 1,maxval(phase_Noutput) select case(constitutive_titanmod_output(o,i)) case('rhoedge', & @@ -710,7 +734,10 @@ do i = 1,maxNinstance 'gdot_slip', & 'velocity_edge', & 'velocity_screw', & - 'total_density' & + 'gdot_slip_edge', & + 'gdot_slip_screw', & + 'stressratio_edge_p', & + 'stressratio_screw_p' & ) mySize = constitutive_titanmod_totalNslip(i) case('twin_fraction', & @@ -1097,7 +1124,7 @@ nt = constitutive_titanmod_totalNtwin(myInstance) ! Need to update this list !* State: 1 : ns rho_edge !* State: ns+1 : 2*ns rho_screw -!* State: 2*ns+1 : 2*ns+nt f +!* State: 2*ns+1 : 2*ns+nt gamma_twin !* State: 2*ns+nt+1 : 3*ns+nt 1/lambda_slip !* State: 3*ns+nt+1 : 4*ns+nt 1/lambda_sliptwin !* State: 4*ns+nt+1 : 4*ns+2*nt 1/lambda_twin @@ -1207,10 +1234,10 @@ integer(pInt) myInstance,myStructure,ns,nt,f,i,j,k,l,m,n,index_myFamily real(pReal) sumf,StressRatio_edge_p,minusStressRatio_edge_p,StressRatio_edge_pminus1,StressRatio_screw_p, & StressRatio_screw_pminus1, StressRatio_r,BoltzmannRatio,DotGamma0, minusStressRatio_screw_p,gdotTotal, & screwvelocity_prefactor,twinStressRatio_p,twinminusStressRatio_p,twinStressRatio_pminus1, & - twinStressRatio_r, twinDotGamma0 + twinStressRatio_r, twinDotGamma0,BoltzmannRatioscrew real(pReal), dimension(3,3,3,3) :: dLp_dTstar3333 real(pReal), dimension(constitutive_titanmod_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: & - gdot_slip,dgdot_dtauslip,tau_slip, edge_velocity, screw_velocity + gdot_slip,dgdot_dtauslip,tau_slip, edge_velocity, screw_velocity,gdot_slip_edge,gdot_slip_screw real(pReal), dimension(constitutive_titanmod_totalNtwin(phase_constitutionInstance(material_phase(g,ip,el)))) :: & gdot_twin,dgdot_dtautwin,tau_twin, twinedge_velocity, twinscrew_velocity,volumefraction_pertwinsystem @@ -1235,6 +1262,8 @@ dLp_dTstar = 0.0_pReal !* Dislocation glide part gdot_slip = 0.0_pReal +gdot_slip_edge = 0.0_pReal +gdot_slip_screw = 0.0_pReal dgdot_dtauslip = 0.0_pReal j = 0_pInt do f = 1,lattice_maxNslipFamily ! loop over all slip families @@ -1245,25 +1274,25 @@ do f = 1,lattice_maxNslipFamily ! loop over all !* Calculation of Lp !* Resolved shear stress on slip system tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure)) - state(g,ip,el)%p(9*ns+3*nt+j)=tau_slip(j) !************************************************* - if(myStructure==3.and.j>3.and.j<25) then ! only for prismatic and pyr systems in hex - screwvelocity_prefactor=constitutive_titanmod_debyfrequency(myInstance)* & +! if(myStructure==3.and.j>3.and.j<13) then ! only for prismatic and pyr systems in hex +! if(myStructure==3) then ! only for prismatic and pyr systems in hex + screwvelocity_prefactor=constitutive_titanmod_debyefrequency(myInstance)* & state(g,ip,el)%p(3*ns+nt+j)*(constitutive_titanmod_burgersPerSlipSystem(j,myInstance)/ & constitutive_titanmod_kinkcriticallength_PerSlipSystem(j,myInstance))**2 - else - screwvelocity_prefactor=constitutive_titanmod_v0s_PerSlipSystem(j,myInstance) - endif +! else +! screwvelocity_prefactor=constitutive_titanmod_v0s_PerSlipSystem(j,myInstance) +! endif !* Stress ratio for edge StressRatio_edge_p = ((abs(tau_slip(j)))/ & ( constitutive_titanmod_tau0e_PerSlipSystem(j,myInstance)+state(g,ip,el)%p(4*ns+nt+j)) & )**constitutive_titanmod_pe_PerSlipSystem(j,myInstance) - !* Stress ratio for screw + !* Stress ratio for screw ! No slip resistance for screw dislocations, only Peierls stress StressRatio_screw_p = ((abs(tau_slip(j)))/ & - ( constitutive_titanmod_tau0s_PerSlipSystem(j,myInstance)+state(g,ip,el)%p(5*ns+nt+j)) & + ( constitutive_titanmod_tau0s_PerSlipSystem(j,myInstance)) & )**constitutive_titanmod_ps_PerSlipSystem(j,myInstance) if((1.0_pReal-StressRatio_edge_p)>0.001_pReal) then @@ -1289,21 +1318,35 @@ do f = 1,lattice_maxNslipFamily ! loop over all !* Boltzmann ratio BoltzmannRatio = constitutive_titanmod_f0_PerSlipSystem(j,myInstance)/(kB*Temperature) + !* Boltzmann ratio for screw + BoltzmannRatioscrew = constitutive_titanmod_kinkf0(myInstance)/(kB*Temperature) + edge_velocity(j) =constitutive_titanmod_v0e_PerSlipSystem(j,myInstance)*exp(-BoltzmannRatio* & (minusStressRatio_edge_p)** & constitutive_titanmod_qe_PerSlipSystem(j,myInstance)) screw_velocity(j) =screwvelocity_prefactor * & ! there is no v0 for screw now because it is included in the prefactor - exp(-BoltzmannRatio*(minusStressRatio_screw_p)** & + exp(-BoltzmannRatioscrew*(minusStressRatio_screw_p)** & constitutive_titanmod_qs_PerSlipSystem(j,myInstance)) - !* Shear rates due to slip - gdot_slip(j) = constitutive_titanmod_burgersPerSlipSystem(j,myInstance)*(state(g,ip,el)%p(j)* & - edge_velocity(j)+state(g,ip,el)%p(ns+j) * screw_velocity(j))* sign(1.0_pReal,tau_slip(j)) + !* Shear rates due to edge slip + gdot_slip_edge(j) = constitutive_titanmod_burgersPerSlipSystem(j,myInstance)*(state(g,ip,el)%p(j)* & + edge_velocity(j))* sign(1.0_pReal,tau_slip(j)) + !* Shear rates due to screw slip + gdot_slip_screw(j) = constitutive_titanmod_burgersPerSlipSystem(j,myInstance)*(state(g,ip,el)%p(ns+j) * & + screw_velocity(j))* sign(1.0_pReal,tau_slip(j)) + !Total shear rate + + gdot_slip(j) = gdot_slip_edge(j) + gdot_slip_screw(j) state(g,ip,el)%p(6*ns+3*nt+j)=gdot_slip(j) state(g,ip,el)%p(7*ns+3*nt+j)=edge_velocity(j) state(g,ip,el)%p(8*ns+3*nt+j)=screw_velocity(j) + state(g,ip,el)%p(9*ns+3*nt+j)=tau_slip(j) + state(g,ip,el)%p(10*ns+3*nt+j)=gdot_slip_edge(j) + state(g,ip,el)%p(11*ns+3*nt+j)=gdot_slip_screw(j) + state(g,ip,el)%p(12*ns+3*nt+j)=StressRatio_edge_p + state(g,ip,el)%p(13*ns+3*nt+j)=StressRatio_screw_p !* Derivatives of shear rates dgdot_dtauslip(j) = constitutive_titanmod_burgersPerSlipSystem(j,myInstance)*(( & @@ -1405,7 +1448,7 @@ do f = 1,lattice_maxNtwinFamily ! loop over all gdot_twin(j) =sign(1.0_pReal,tau_twin(j))*constitutive_titanmod_twingamma0_PerTwinSystem(j,myInstance)* & exp(-BoltzmannRatio*(twinminusStressRatio_p)**constitutive_titanmod_twinq_PerTwinSystem(j,myInstance)) - state(g,ip,el)%p(6*ns+2*nt+j)=nt + state(g,ip,el)%p(6*ns+2*nt+j)=gdot_twin(j) !* Derivatives of shear rates in twin @@ -1485,7 +1528,7 @@ constitutive_titanmod_dotState integer(pInt) MyInstance,MyStructure,ns,nt,f,i,j,k,index_myFamily,s,t real(pReal) sumf,StressRatio_edge_p,minusStressRatio_edge_p,StressRatio_pminus1,BoltzmannRatio,DotGamma0,& EdgeDipMinDistance,AtomicVolume,VacancyDiffusion,StressRatio_r,StressRatio_screw_p,minusStressRatio_screw_p, & - screwvelocity_prefactor,twinStressRatio_p,twinminusStressRatio_p,twinStressRatio_pminus1, & + twinStressRatio_p,twinminusStressRatio_p,twinStressRatio_pminus1, & twinDotGamma0 real(pReal), dimension(constitutive_titanmod_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: & gdot_slip,tau_slip,DotRhoEdgeGeneration,EdgeDipDistance,DotRhoEdgeAnnihilation,DotRhoScrewAnnihilation,& @@ -1517,17 +1560,17 @@ constitutive_titanmod_dotState = 0.0_pReal j = j+1_pInt !* Multiplication of edge dislocations - DotRhoEdgeGeneration(j) = 2.0_pReal*(state(g,ip,el)%p(ns+j)*state(g,ip,el)%p(8*ns+3*nt+j)/state(g,ip,el)%p(3*ns+nt+j)) + DotRhoEdgeGeneration(j) = (state(g,ip,el)%p(ns+j)*state(g,ip,el)%p(8*ns+3*nt+j)/state(g,ip,el)%p(3*ns+nt+j)) !* Multiplication of screw dislocations - DotRhoScrewGeneration(j) = 2.0_pReal*(state(g,ip,el)%p(j)*state(g,ip,el)%p(7*ns+3*nt+j)/state(g,ip,el)%p(2*ns+nt+j)) + DotRhoScrewGeneration(j) = (state(g,ip,el)%p(j)*state(g,ip,el)%p(7*ns+3*nt+j)/state(g,ip,el)%p(2*ns+nt+j)) !* Annihilation of edge dislocations DotRhoEdgeAnnihilation(j) = -((state(g,ip,el)%p(j))**2)* & - constitutive_titanmod_capre_PerSlipSystem(j,myInstance)*state(g,ip,el)%p(7*ns+3*nt+j) + constitutive_titanmod_capre_PerSlipSystem(j,myInstance)*state(g,ip,el)%p(7*ns+3*nt+j)/2.0_pReal !* Annihilation of screw dislocations DotRhoScrewAnnihilation(j) = -((state(g,ip,el)%p(ns+j))**2)* & - constitutive_titanmod_caprs_PerSlipSystem(j,myInstance)*state(g,ip,el)%p(8*ns+3*nt+j) + constitutive_titanmod_caprs_PerSlipSystem(j,myInstance)*state(g,ip,el)%p(8*ns+3*nt+j)/2.0_pReal !* Edge dislocation density rate of change constitutive_titanmod_dotState(j) = & @@ -1581,15 +1624,10 @@ do f = 1,lattice_maxNtwinFamily ! loop over all !* Boltzmann ratio BoltzmannRatio = constitutive_titanmod_twinf0_PerTwinSystem(j,myInstance)/(kB*Temperature) - -! if (tau_slip(j) == 0.0_pReal) then -! edge_velocity(j) = 0.0_pReal -! screw_velocity(j) = 0.0_pReal -! else + gdot_twin(j) =constitutive_titanmod_twingamma0_PerTwinSystem(j,myInstance)*exp(-BoltzmannRatio* & (twinminusStressRatio_p)** & constitutive_titanmod_twinq_PerTwinSystem(j,myInstance))*sign(1.0_pReal,tau_twin(j)) -! endif constitutive_titanmod_dotState(2*ns+j)=gdot_twin(j) @@ -1598,16 +1636,11 @@ enddo !write(6,*) '#DOTSTATE#' !write(6,*) -!write(6,'(a,/,4(3(f30.20,x)/))') 'tau slip',tau_slip -!write(6,'(a,/,4(3(f30.20,x)/))') 'gamma slip',gdot_slip -!write(6,'(a,/,4(3(f30.20,x)/))') 'rho_edge',state(g,ip,el)%p(1:ns) -!write(6,'(a,/,4(3(f30.20,x)/))') 'Threshold Slip Edge', state(g,ip,el)%p(5*ns+3*nt+1:6*ns+3*nt) -!write(6,'(a,/,4(3(f30.20,x)/))') 'Threshold Slip Screw', state(g,ip,el)%p(6*ns+3*nt+1:7*ns+3*nt) !write(6,'(a,/,4(3(f30.20,x)/))') 'EdgeGeneration',DotRhoEdgeGeneration !write(6,'(a,/,4(3(f30.20,x)/))') 'ScrewGeneration',DotRhoScrewGeneration !write(6,'(a,/,4(3(f30.20,x)/))') 'EdgeAnnihilation',DotRhoEdgeAnnihilation !write(6,'(a,/,4(3(f30.20,x)/))') 'ScrewAnnihilation',DotRhoScrewAnnihilation -!write(6,'(a,/,4(3(f30.20,x)/))') 'DipClimb',DotRhoEdgeDipClimb + return end function @@ -1668,8 +1701,7 @@ real(pReal), intent(in) :: dt,Temperature real(pReal), dimension(6), intent(in) :: Tstar_v type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state integer(pInt) myInstance,myStructure,ns,nt,f,o,i,c,j,index_myFamily -real(pReal) sumf,tau,StressRatio_edge_p,StressRatio_screw_p,StressRatio_pminus1,BoltzmannRatio,DotGamma0,StressRatio_r, & - gdot_slip,dgdot_dtauslip +real(pReal) sumf real(pReal), dimension(constitutive_titanmod_sizePostResults(phase_constitutionInstance(material_phase(g,ip,el)))) :: & constitutive_titanmod_postResults real(pReal), dimension(constitutive_titanmod_totalNtwin(phase_constitutionInstance(material_phase(g,ip,el)))) :: & @@ -1705,18 +1737,6 @@ do o = 1,phase_Noutput(material_phase(g,ip,el)) case ('rhoscrew') constitutive_titanmod_postResults(c+1:c+ns) = state(g,ip,el)%p(ns+1:2*ns) c = c + ns - case ('gdot_slip') - constitutive_titanmod_postResults(c+1:c+ns) = state(g,ip,el)%p((6*ns+3*nt+1):(7*ns+3*nt)) - c = c + ns - case ('gdot_twin') - constitutive_titanmod_postResults(c+1:c+nt) = state(g,ip,el)%p((6*ns+2*nt+1):(6*ns+3*nt)) - c = c + nt - case ('velocity_edge') - constitutive_titanmod_postResults(c+1:c+ns) = state(g,ip,el)%p((7*ns+3*nt+1):(8*ns+3*nt)) - c = c + ns - case ('velocity_screw') - constitutive_titanmod_postResults(c+1:c+ns) = state(g,ip,el)%p((8*ns+3*nt+1):(9*ns+3*nt)) - c = c + ns case ('segment_edge') constitutive_titanmod_postResults(c+1:c+ns) = state(g,ip,el)%p((2*ns+nt+1):(3*ns+nt)) c = c + ns @@ -1730,11 +1750,49 @@ do o = 1,phase_Noutput(material_phase(g,ip,el)) constitutive_titanmod_postResults(c+1:c+ns) = state(g,ip,el)%p((5*ns+nt+1):(6*ns+nt)) c = c + ns case ('tau_slip') - constitutive_titanmod_postResults(c+1:c+ns) = state(g,ip,el)%p((9*ns+3*nt+1):(10*ns+3*nt)) - c=c + ns + constitutive_titanmod_postResults(c+1:c+ns) = abs(state(g,ip,el)%p((9*ns+3*nt+1):(10*ns+3*nt))) + c = c + ns + case ('gdot_slip') + constitutive_titanmod_postResults(c+1:c+ns) = abs(state(g,ip,el)%p((6*ns+3*nt+1):(7*ns+3*nt))) + c = c + ns + case ('velocity_edge') + constitutive_titanmod_postResults(c+1:c+ns) = state(g,ip,el)%p((7*ns+3*nt+1):(8*ns+3*nt)) + c = c + ns + case ('velocity_screw') + constitutive_titanmod_postResults(c+1:c+ns) = state(g,ip,el)%p((8*ns+3*nt+1):(9*ns+3*nt)) + c = c + ns + case ('gdot_slip_edge') + constitutive_titanmod_postResults(c+1:c+ns) = abs(state(g,ip,el)%p((10*ns+3*nt+1):(11*ns+3*nt))) + c = c + ns + case ('gdot_slip_screw') + constitutive_titanmod_postResults(c+1:c+ns) = abs(state(g,ip,el)%p((11*ns+3*nt+1):(12*ns+3*nt))) + c = c + ns + case ('stressratio_edge_p') + constitutive_titanmod_postResults(c+1:c+ns) = abs(state(g,ip,el)%p((12*ns+3*nt+1):(13*ns+3*nt))) + c = c + ns + case ('stressratio_screw_p') + constitutive_titanmod_postResults(c+1:c+ns) = abs(state(g,ip,el)%p((13*ns+3*nt+1):(14*ns+3*nt))) + c = c + ns + case ('gdot_twin') + constitutive_titanmod_postResults(c+1:c+nt) = abs(state(g,ip,el)%p((6*ns+2*nt+1):(6*ns+3*nt))) + c = c + nt case ('twin_fraction') constitutive_titanmod_postResults(c+1:c+nt) = abs(volumefraction_pertwinsystem(1:nt)) c = c + nt + ! 'rhoedge', & + ! 'rhoscrew', & + ! 'segment_edge', & + ! 'segment_screw', & + ! 'resistance_edge', & + ! 'resistance_screw', & + ! 'tau_slip', & + ! 'gdot_slip', & + ! 'velocity_edge', & + ! 'velocity_screw', & + ! 'gdot_slip_edge' & + ! 'gdot_slip_screw' & + ! 'StressRatio_edge_p' & + ! 'StressRatio_screw_p' end select enddo