From 7d05845d5d3586928526cf7b11ba73b3a450eb77 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Mon, 12 Oct 2009 16:26:16 +0000 Subject: [PATCH] # now $ID aware # LF as line separation and some tidy up --- code/constitutive_dislotwin.f90 | 886 +++++++++++++++++--------------- 1 file changed, 467 insertions(+), 419 deletions(-) diff --git a/code/constitutive_dislotwin.f90 b/code/constitutive_dislotwin.f90 index 12d6e5747..9a22d7655 100644 --- a/code/constitutive_dislotwin.f90 +++ b/code/constitutive_dislotwin.f90 @@ -1,4 +1,4 @@ -!* $Id: constitutive_dislotwin.f90 412 2009-09-18 15:37:14Z MPIE\c.kords $ +!* $Id$ !************************************ !* Module: CONSTITUTIVE * !************************************ @@ -12,16 +12,16 @@ implicit none !* Lists of states and physical parameters character(len=*), parameter :: constitutive_dislotwin_label = 'dislotwin' character(len=18), dimension(2), parameter:: constitutive_dislotwin_listBasicSlipStates = (/'rhoEdge ', & - 'rhoEdgeDip'/) -character(len=18), dimension(1), parameter:: constitutive_dislotwin_listBasicTwinStates = (/'twinFraction'/) + 'rhoEdgeDip'/) +character(len=18), dimension(1), parameter:: constitutive_dislotwin_listBasicTwinStates = (/'twinFraction'/) character(len=18), dimension(4), parameter:: constitutive_dislotwin_listDependentSlipStates =(/'invLambdaSlip ', & - 'invLambdaSlipTwin', & - 'meanFreePathSlip ', & - 'tauSlipThreshold '/) + 'invLambdaSlipTwin', & + 'meanFreePathSlip ', & + 'tauSlipThreshold '/) character(len=18), dimension(4), parameter:: constitutive_dislotwin_listDependentTwinStates =(/'invLambdaTwin ', & - 'meanFreePathTwin', & - 'tauTwinThreshold', & - 'twinVolume '/) + 'meanFreePathTwin', & + 'tauTwinThreshold', & + 'twinVolume '/) real(pReal), parameter :: kB = 1.38e-23_pReal ! Boltzmann constant in J/Kelvin !* Definition of global variables @@ -78,16 +78,16 @@ real(pReal), dimension(:,:), allocatable :: constitutive_dislotwin constitutive_dislotwin_twinsizePerTwinFamily, & ! twin thickness [m] for each twin family and instance constitutive_dislotwin_twinsizePerTwinSystem, & ! twin thickness [m] for each twin system and instance constitutive_dislotwin_CLambdaSlipPerSlipFamily, & ! Adj. parameter for distance between 2 forest dislocations for each slip family and instance - constitutive_dislotwin_CLambdaSlipPerSlipSystem, & ! Adj. parameter for distance between 2 forest dislocations for each slip system and instance - constitutive_dislotwin_interactionSlipSlip, & ! coefficients for slip-slip interaction for each interaction type and instance + constitutive_dislotwin_CLambdaSlipPerSlipSystem, & ! Adj. parameter for distance between 2 forest dislocations for each slip system and instance + constitutive_dislotwin_interactionSlipSlip, & ! coefficients for slip-slip interaction for each interaction type and instance constitutive_dislotwin_interactionSlipTwin, & ! coefficients for slip-twin interaction for each interaction type and instance constitutive_dislotwin_interactionTwinSlip, & ! coefficients for twin-slip interaction for each interaction type and instance constitutive_dislotwin_interactionTwinTwin ! coefficients for twin-twin interaction for each interaction type and instance real(pReal), dimension(:,:,:), allocatable :: constitutive_dislotwin_interactionMatrixSlipSlip, & ! interaction matrix of the different slip systems for each instance constitutive_dislotwin_interactionMatrixSlipTwin, & ! interaction matrix of slip systems with twin systems for each instance constitutive_dislotwin_interactionMatrixTwinSlip, & ! interaction matrix of twin systems with slip systems for each instance - constitutive_dislotwin_interactionMatrixTwinTwin, & ! interaction matrix of the different twin systems for each instance - constitutive_dislotwin_forestProjectionEdge ! matrix of forest projections of edge dislocations for each instance + constitutive_dislotwin_interactionMatrixTwinTwin, & ! interaction matrix of the different twin systems for each instance + constitutive_dislotwin_forestProjectionEdge ! matrix of forest projections of edge dislocations for each instance CONTAINS !**************************************** !* - constitutive_dislotwin_init @@ -122,53 +122,52 @@ character(len=1024) line !write(6,*) !write(6,'(a20,a20,a12)') '<<<+- constitutive_',constitutive_dislotwin_label,' init -+>>>' -!write(6,*) '$Id: constitutive_dislotwin.f90 412 2009-09-18 15:37:14Z MPIE\c.kords $' +!write(6,*) '$Id$' !write(6,*) maxNinstance = count(phase_constitution == constitutive_dislotwin_label) if (maxNinstance == 0) return - + !* Space allocation for global variables allocate(constitutive_dislotwin_sizeDotState(maxNinstance)) allocate(constitutive_dislotwin_sizeState(maxNinstance)) allocate(constitutive_dislotwin_sizePostResults(maxNinstance)) allocate(constitutive_dislotwin_sizePostResult(maxval(phase_Noutput),maxNinstance)) -allocate(constitutive_dislotwin_output(maxval(phase_Noutput),maxNinstance)) -constitutive_dislotwin_sizeDotState = 0_pInt -constitutive_dislotwin_sizeState = 0_pInt -constitutive_dislotwin_sizePostResults = 0_pInt -constitutive_dislotwin_sizePostResult = 0_pInt -constitutive_dislotwin_output = '' +allocate(constitutive_dislotwin_output(maxval(phase_Noutput),maxNinstance)) +constitutive_dislotwin_sizeDotState = 0_pInt +constitutive_dislotwin_sizeState = 0_pInt +constitutive_dislotwin_sizePostResults = 0_pInt +constitutive_dislotwin_sizePostResult = 0_pInt +constitutive_dislotwin_output = '' allocate(constitutive_dislotwin_structureName(maxNinstance)) allocate(constitutive_dislotwin_structure(maxNinstance)) -allocate(constitutive_dislotwin_Nslip(lattice_maxNslipFamily,maxNinstance)) +allocate(constitutive_dislotwin_Nslip(lattice_maxNslipFamily,maxNinstance)) allocate(constitutive_dislotwin_Ntwin(lattice_maxNtwinFamily,maxNinstance)) -allocate(constitutive_dislotwin_slipFamily(lattice_maxNslip,maxNinstance)) +allocate(constitutive_dislotwin_slipFamily(lattice_maxNslip,maxNinstance)) allocate(constitutive_dislotwin_twinFamily(lattice_maxNtwin,maxNinstance)) allocate(constitutive_dislotwin_slipSystemLattice(lattice_maxNslip,maxNinstance)) allocate(constitutive_dislotwin_twinSystemLattice(lattice_maxNtwin,maxNinstance)) allocate(constitutive_dislotwin_totalNslip(maxNinstance)) allocate(constitutive_dislotwin_totalNtwin(maxNinstance)) -constitutive_dislotwin_structureName = '' -constitutive_dislotwin_structure = 0_pInt -constitutive_dislotwin_Nslip = 0_pInt -constitutive_dislotwin_Ntwin = 0_pInt -constitutive_dislotwin_slipFamily = 0_pInt +constitutive_dislotwin_structureName = '' +constitutive_dislotwin_structure = 0_pInt +constitutive_dislotwin_Nslip = 0_pInt +constitutive_dislotwin_Ntwin = 0_pInt +constitutive_dislotwin_slipFamily = 0_pInt constitutive_dislotwin_twinFamily = 0_pInt constitutive_dislotwin_slipSystemLattice = 0.0_pReal -constitutive_dislotwin_twinSystemLattice = 0.0_pReal -constitutive_dislotwin_totalNslip = 0_pInt -constitutive_dislotwin_totalNtwin = 0_pInt - -allocate(constitutive_dislotwin_CoverA(maxNinstance)) -allocate(constitutive_dislotwin_C11(maxNinstance)) -allocate(constitutive_dislotwin_C12(maxNinstance)) -allocate(constitutive_dislotwin_C13(maxNinstance)) -allocate(constitutive_dislotwin_C33(maxNinstance)) -allocate(constitutive_dislotwin_C44(maxNinstance)) -allocate(constitutive_dislotwin_Gmod(maxNinstance)) -allocate(constitutive_dislotwin_CAtomicVolume(maxNinstance)) +constitutive_dislotwin_twinSystemLattice = 0.0_pReal +constitutive_dislotwin_totalNslip = 0_pInt +constitutive_dislotwin_totalNtwin = 0_pInt +allocate(constitutive_dislotwin_CoverA(maxNinstance)) +allocate(constitutive_dislotwin_C11(maxNinstance)) +allocate(constitutive_dislotwin_C12(maxNinstance)) +allocate(constitutive_dislotwin_C13(maxNinstance)) +allocate(constitutive_dislotwin_C33(maxNinstance)) +allocate(constitutive_dislotwin_C44(maxNinstance)) +allocate(constitutive_dislotwin_Gmod(maxNinstance)) +allocate(constitutive_dislotwin_CAtomicVolume(maxNinstance)) allocate(constitutive_dislotwin_D0(maxNinstance)) allocate(constitutive_dislotwin_Qsd(maxNinstance)) allocate(constitutive_dislotwin_GrainSize(maxNinstance)) @@ -176,19 +175,19 @@ allocate(constitutive_dislotwin_p(maxNinstance)) allocate(constitutive_dislotwin_q(maxNinstance)) allocate(constitutive_dislotwin_MaxTwinFraction(maxNinstance)) allocate(constitutive_dislotwin_r(maxNinstance)) -allocate(constitutive_dislotwin_CEdgeDipMinDistance(maxNinstance)) -allocate(constitutive_dislotwin_Cmfptwin(maxNinstance)) +allocate(constitutive_dislotwin_CEdgeDipMinDistance(maxNinstance)) +allocate(constitutive_dislotwin_Cmfptwin(maxNinstance)) allocate(constitutive_dislotwin_Cthresholdtwin(maxNinstance)) -allocate(constitutive_dislotwin_relevantRho(maxNinstance)) -allocate(constitutive_dislotwin_Cslip_66(6,6,maxNinstance)) -allocate(constitutive_dislotwin_Cslip_3333(3,3,3,3,maxNinstance)) -constitutive_dislotwin_CoverA = 0.0_pReal -constitutive_dislotwin_C11 = 0.0_pReal -constitutive_dislotwin_C12 = 0.0_pReal -constitutive_dislotwin_C13 = 0.0_pReal -constitutive_dislotwin_C33 = 0.0_pReal -constitutive_dislotwin_C44 = 0.0_pReal -constitutive_dislotwin_Gmod = 0.0_pReal +allocate(constitutive_dislotwin_relevantRho(maxNinstance)) +allocate(constitutive_dislotwin_Cslip_66(6,6,maxNinstance)) +allocate(constitutive_dislotwin_Cslip_3333(3,3,3,3,maxNinstance)) +constitutive_dislotwin_CoverA = 0.0_pReal +constitutive_dislotwin_C11 = 0.0_pReal +constitutive_dislotwin_C12 = 0.0_pReal +constitutive_dislotwin_C13 = 0.0_pReal +constitutive_dislotwin_C33 = 0.0_pReal +constitutive_dislotwin_C44 = 0.0_pReal +constitutive_dislotwin_Gmod = 0.0_pReal constitutive_dislotwin_CAtomicVolume = 0.0_pReal constitutive_dislotwin_D0 = 0.0_pReal constitutive_dislotwin_Qsd = 0.0_pReal @@ -197,15 +196,14 @@ constitutive_dislotwin_p = 0.0_pReal constitutive_dislotwin_q = 0.0_pReal constitutive_dislotwin_MaxTwinFraction = 0.0_pReal constitutive_dislotwin_r = 0.0_pReal -constitutive_dislotwin_CEdgeDipMinDistance = 0.0_pReal -constitutive_dislotwin_Cmfptwin = 0.0_pReal +constitutive_dislotwin_CEdgeDipMinDistance = 0.0_pReal +constitutive_dislotwin_Cmfptwin = 0.0_pReal constitutive_dislotwin_Cthresholdtwin = 0.0_pReal -constitutive_dislotwin_relevantRho = 0.0_pReal -constitutive_dislotwin_Cslip_66 = 0.0_pReal -constitutive_dislotwin_Cslip_3333 = 0.0_pReal - +constitutive_dislotwin_relevantRho = 0.0_pReal +constitutive_dislotwin_Cslip_66 = 0.0_pReal +constitutive_dislotwin_Cslip_3333 = 0.0_pReal allocate(constitutive_dislotwin_rhoEdge0(lattice_maxNslipFamily,maxNinstance)) -allocate(constitutive_dislotwin_rhoEdgeDip0(lattice_maxNslipFamily,maxNinstance)) +allocate(constitutive_dislotwin_rhoEdgeDip0(lattice_maxNslipFamily,maxNinstance)) allocate(constitutive_dislotwin_burgersPerSlipFamily(lattice_maxNslipFamily,maxNinstance)) allocate(constitutive_dislotwin_burgersPerTwinFamily(lattice_maxNtwinFamily,maxNinstance)) allocate(constitutive_dislotwin_QedgePerSlipFamily(lattice_maxNslipFamily,maxNinstance)) @@ -213,25 +211,26 @@ allocate(constitutive_dislotwin_v0PerSlipFamily(lattice_maxNslipFamily,maxNinsta allocate(constitutive_dislotwin_Ndot0PerTwinFamily(lattice_maxNtwinFamily,maxNinstance)) allocate(constitutive_dislotwin_twinsizePerTwinFamily(lattice_maxNtwinFamily,maxNinstance)) allocate(constitutive_dislotwin_CLambdaSlipPerSlipFamily(lattice_maxNslipFamily,maxNinstance)) -constitutive_dislotwin_rhoEdge0 = 0.0_pReal -constitutive_dislotwin_rhoEdgeDip0 = 0.0_pReal -constitutive_dislotwin_burgersPerSlipFamily = 0.0_pReal -constitutive_dislotwin_burgersPerTwinFamily = 0.0_pReal -constitutive_dislotwin_QedgePerSlipFamily = 0.0_pReal +constitutive_dislotwin_rhoEdge0 = 0.0_pReal +constitutive_dislotwin_rhoEdgeDip0 = 0.0_pReal +constitutive_dislotwin_burgersPerSlipFamily = 0.0_pReal +constitutive_dislotwin_burgersPerTwinFamily = 0.0_pReal +constitutive_dislotwin_QedgePerSlipFamily = 0.0_pReal constitutive_dislotwin_v0PerSlipFamily = 0.0_pReal constitutive_dislotwin_Ndot0PerTwinFamily = 0.0_pReal -constitutive_dislotwin_twinsizePerTwinFamily = 0.0_pReal +constitutive_dislotwin_twinsizePerTwinFamily = 0.0_pReal constitutive_dislotwin_CLambdaSlipPerSlipFamily = 0.0_pReal - allocate(constitutive_dislotwin_interactionSlipSlip(lattice_maxNinteraction,maxNinstance)) allocate(constitutive_dislotwin_interactionSlipTwin(lattice_maxNinteraction,maxNinstance)) allocate(constitutive_dislotwin_interactionTwinSlip(lattice_maxNinteraction,maxNinstance)) -allocate(constitutive_dislotwin_interactionTwinTwin(lattice_maxNinteraction,maxNinstance)) -constitutive_dislotwin_interactionSlipSlip = 0.0_pReal -constitutive_dislotwin_interactionSlipTwin = 0.0_pReal +allocate(constitutive_dislotwin_interactionTwinTwin(lattice_maxNinteraction,maxNinstance)) +constitutive_dislotwin_interactionSlipSlip = 0.0_pReal +constitutive_dislotwin_interactionSlipTwin = 0.0_pReal constitutive_dislotwin_interactionTwinSlip = 0.0_pReal -constitutive_dislotwin_interactionTwinTwin = 0.0_pReal - +constitutive_dislotwin_interactionTwinTwin = 0.0_pReal + + + !* Readout data from material.config file rewind(file) line = '' @@ -272,41 +271,52 @@ do ! read thru sections of case ('c44') constitutive_dislotwin_C44(i) = IO_floatValue(line,positions,2) case ('nslip') - forall (j = 1:lattice_maxNslipFamily) constitutive_dislotwin_Nslip(j,i) = IO_intValue(line,positions,1+j) + forall (j = 1:lattice_maxNslipFamily) & + constitutive_dislotwin_Nslip(j,i) = IO_intValue(line,positions,1+j) case ('ntwin') - forall (j = 1:lattice_maxNtwinFamily) constitutive_dislotwin_Ntwin(j,i) = IO_intValue(line,positions,1+j) - case ('rhoedge0') - forall (j = 1:lattice_maxNslipFamily) constitutive_dislotwin_rhoEdge0(j,i) = IO_floatValue(line,positions,1+j) - case ('rhoedgedip0') - forall (j = 1:lattice_maxNslipFamily) constitutive_dislotwin_rhoEdgeDip0(j,i) = IO_floatValue(line,positions,1+j) + forall (j = 1:lattice_maxNtwinFamily) & + constitutive_dislotwin_Ntwin(j,i) = IO_intValue(line,positions,1+j) + case ('rhoedge0') + forall (j = 1:lattice_maxNslipFamily) & + constitutive_dislotwin_rhoEdge0(j,i) = IO_floatValue(line,positions,1+j) + case ('rhoedgedip0') + forall (j = 1:lattice_maxNslipFamily) & + constitutive_dislotwin_rhoEdgeDip0(j,i) = IO_floatValue(line,positions,1+j) case ('slipburgers') - forall (j = 1:lattice_maxNslipFamily) constitutive_dislotwin_burgersPerSlipFamily(j,i) = IO_floatValue(line,positions,1+j) - case ('twinburgers') - forall (j = 1:lattice_maxNtwinFamily) constitutive_dislotwin_burgersPerTwinFamily(j,i) = IO_floatValue(line,positions,1+j) + forall (j = 1:lattice_maxNslipFamily) & + constitutive_dislotwin_burgersPerSlipFamily(j,i) = IO_floatValue(line,positions,1+j) + case ('twinburgers') + forall (j = 1:lattice_maxNtwinFamily) & + constitutive_dislotwin_burgersPerTwinFamily(j,i) = IO_floatValue(line,positions,1+j) case ('qedge') - forall (j = 1:lattice_maxNslipFamily) constitutive_dislotwin_QedgePerSlipFamily(j,i) = IO_floatValue(line,positions,1+j) - case ('v0') - forall (j = 1:lattice_maxNslipFamily) constitutive_dislotwin_v0PerSlipFamily(j,i) = IO_floatValue(line,positions,1+j) - case ('ndot0') - forall (j = 1:lattice_maxNtwinFamily) constitutive_dislotwin_Ndot0PerTwinFamily(j,i) = IO_floatValue(line,positions,1+j) - case ('twinsize') - forall (j = 1:lattice_maxNtwinFamily) constitutive_dislotwin_twinsizePerTwinFamily(j,i) = IO_floatValue(line,positions,1+j) + forall (j = 1:lattice_maxNslipFamily) & + constitutive_dislotwin_QedgePerSlipFamily(j,i) = IO_floatValue(line,positions,1+j) + case ('v0') + forall (j = 1:lattice_maxNslipFamily) & + constitutive_dislotwin_v0PerSlipFamily(j,i) = IO_floatValue(line,positions,1+j) + case ('ndot0') + forall (j = 1:lattice_maxNtwinFamily) & + constitutive_dislotwin_Ndot0PerTwinFamily(j,i) = IO_floatValue(line,positions,1+j) + case ('twinsize') + forall (j = 1:lattice_maxNtwinFamily) & + constitutive_dislotwin_twinsizePerTwinFamily(j,i) = IO_floatValue(line,positions,1+j) case ('clambdaslip') - forall (j = 1:lattice_maxNslipFamily) constitutive_dislotwin_CLambdaSlipPerSlipFamily(j,i) = IO_floatValue(line,positions,1+j) - case ('grainsize') + forall (j = 1:lattice_maxNslipFamily) & + constitutive_dislotwin_CLambdaSlipPerSlipFamily(j,i) = IO_floatValue(line,positions,1+j) + case ('grainsize') constitutive_dislotwin_GrainSize(i) = IO_floatValue(line,positions,2) case ('maxtwinfraction') constitutive_dislotwin_MaxTwinFraction(i) = IO_floatValue(line,positions,2) - case ('pexponent') - constitutive_dislotwin_p(i) = IO_floatValue(line,positions,2) - case ('qexponent') - constitutive_dislotwin_q(i) = IO_floatValue(line,positions,2) - case ('rexponent') - constitutive_dislotwin_r(i) = IO_floatValue(line,positions,2) - case ('d0') - constitutive_dislotwin_D0(i) = IO_floatValue(line,positions,2) - case ('qsd') - constitutive_dislotwin_Qsd(i) = IO_floatValue(line,positions,2) + case ('pexponent') + constitutive_dislotwin_p(i) = IO_floatValue(line,positions,2) + case ('qexponent') + constitutive_dislotwin_q(i) = IO_floatValue(line,positions,2) + case ('rexponent') + constitutive_dislotwin_r(i) = IO_floatValue(line,positions,2) + case ('d0') + constitutive_dislotwin_D0(i) = IO_floatValue(line,positions,2) + case ('qsd') + constitutive_dislotwin_Qsd(i) = IO_floatValue(line,positions,2) case ('relevantrho') constitutive_dislotwin_relevantRho(i) = IO_floatValue(line,positions,2) case ('cmfptwin') @@ -319,7 +329,7 @@ do ! read thru sections of constitutive_dislotwin_CAtomicVolume(i) = IO_floatValue(line,positions,2) case ('interactionslipslip') forall (j = 1:lattice_maxNinteraction) & - constitutive_dislotwin_interactionSlipSlip(j,i) = IO_floatValue(line,positions,1+j) + constitutive_dislotwin_interactionSlipSlip(j,i) = IO_floatValue(line,positions,1+j) case ('interactionsliptwin') forall (j = 1:lattice_maxNinteraction) & constitutive_dislotwin_interactionSlipTwin(j,i) = IO_floatValue(line,positions,1+j) @@ -345,16 +355,16 @@ enddo if (sum(constitutive_dislotwin_Ntwin(:,i)) < 0_pInt) call IO_error(225) !*** do f = 1,lattice_maxNslipFamily if (constitutive_dislotwin_Nslip(f,i) > 0_pInt) then - if (constitutive_dislotwin_rhoEdge0(f,i) < 0.0_pReal) call IO_error(220) - if (constitutive_dislotwin_rhoEdgeDip0(f,i) < 0.0_pReal) call IO_error(220) - if (constitutive_dislotwin_burgersPerSlipFamily(f,i) <= 0.0_pReal) call IO_error(221) - if (constitutive_dislotwin_v0PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(226) + if (constitutive_dislotwin_rhoEdge0(f,i) < 0.0_pReal) call IO_error(220) + if (constitutive_dislotwin_rhoEdgeDip0(f,i) < 0.0_pReal) call IO_error(220) + if (constitutive_dislotwin_burgersPerSlipFamily(f,i) <= 0.0_pReal) call IO_error(221) + if (constitutive_dislotwin_v0PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(226) endif enddo do f = 1,lattice_maxNtwinFamily if (constitutive_dislotwin_Nslip(f,i) > 0_pInt) then - if (constitutive_dislotwin_burgersPerTwinFamily(f,i) <= 0.0_pReal) call IO_error(221) !*** - if (constitutive_dislotwin_Ndot0PerTwinFamily(f,i) < 0.0_pReal) call IO_error(226) !*** + if (constitutive_dislotwin_burgersPerTwinFamily(f,i) <= 0.0_pReal) call IO_error(221) !*** + if (constitutive_dislotwin_Ndot0PerTwinFamily(f,i) < 0.0_pReal) call IO_error(226) !*** endif enddo ! if (any(constitutive_dislotwin_interactionSlipSlip(1:maxval(lattice_interactionSlipSlip(:,:,myStructure)),i) < 1.0_pReal)) call IO_error(229) @@ -373,7 +383,8 @@ enddo !* Allocation of variables whose size depends on the total number of active slip systems maxTotalNslip = maxval(constitutive_dislotwin_totalNslip) -maxTotalNtwin = maxval(constitutive_dislotwin_totalNtwin) +maxTotalNtwin = maxval(constitutive_dislotwin_totalNtwin) + allocate(constitutive_dislotwin_burgersPerSlipSystem(maxTotalNslip, maxNinstance)) allocate(constitutive_dislotwin_burgersPerTwinSystem(maxTotalNtwin, maxNinstance)) @@ -439,29 +450,30 @@ do i = 1,maxNinstance !* Determine size of postResults array do o = 1,maxval(phase_Noutput) select case(constitutive_dislotwin_output(o,i)) - case('edge_density', & - 'dipole_density', & - 'shear_rate_slip', & - 'mfp_slip', & - 'resolved_stress_slip', & - 'threshold_stress_slip' & - ) - mySize = constitutive_dislotwin_totalNslip(i) - case('twin_fraction', & - 'shear_rate_twin', & - 'mfp_twin', & - 'resolved_stress_twin', & - 'threshold_stress_twin' & - ) - mySize = constitutive_dislotwin_totalNtwin(i) - case default - mySize = 0_pInt + case('edge_density', & + 'dipole_density', & + + 'shear_rate_slip', & + 'mfp_slip', & + 'resolved_stress_slip', & + 'threshold_stress_slip' & + ) + mySize = constitutive_dislotwin_totalNslip(i) + case('twin_fraction', & + 'shear_rate_twin', & + 'mfp_twin', & + 'resolved_stress_twin', & + 'threshold_stress_twin' & + ) + mySize = constitutive_dislotwin_totalNtwin(i) + case default + mySize = 0_pInt end select - if (mySize > 0_pInt) then ! any meaningful output found - constitutive_dislotwin_sizePostResult(o,i) = mySize - constitutive_dislotwin_sizePostResults(i) = constitutive_dislotwin_sizePostResults(i) + mySize - endif + if (mySize > 0_pInt) then ! any meaningful output found + constitutive_dislotwin_sizePostResult(o,i) = mySize + constitutive_dislotwin_sizePostResults(i) = constitutive_dislotwin_sizePostResults(i) + mySize + endif enddo !* Elasticity matrix and shear modulus according to material.config @@ -496,16 +508,16 @@ do i = 1,maxNinstance do j=1,lattice_maxNtwinFamily do k=1,constitutive_dislotwin_Ntwin(j,i) do l=1,3 ; do m=1,3 ; do n=1,3 ; do o=1,3 ; do p=1,3 ; do q=1,3 ; do r=1,3 ; do s=1,3 - constitutive_dislotwin_Ctwin_3333(l,m,n,o,sum(constitutive_dislotwin_Nslip(1:j-1,i))+k,i) = & - constitutive_dislotwin_Ctwin_3333(l,m,n,o,sum(constitutive_dislotwin_Nslip(1:j-1,i))+k,i) + & - constitutive_dislotwin_Cslip_3333(p,q,r,s,i)*& - lattice_Qtwin(l,p,sum(lattice_NslipSystem(1:j-1,myStructure))+k,myStructure)* & - lattice_Qtwin(m,q,sum(lattice_NslipSystem(1:j-1,myStructure))+k,myStructure)* & - lattice_Qtwin(n,r,sum(lattice_NslipSystem(1:j-1,myStructure))+k,myStructure)* & - lattice_Qtwin(o,s,sum(lattice_NslipSystem(1:j-1,myStructure))+k,myStructure) - enddo ; enddo ; enddo ; enddo ; enddo ; enddo ; enddo ; enddo + constitutive_dislotwin_Ctwin_3333(l,m,n,o,sum(constitutive_dislotwin_Nslip(1:j-1,i))+k,i) = & + constitutive_dislotwin_Ctwin_3333(l,m,n,o,sum(constitutive_dislotwin_Nslip(1:j-1,i))+k,i) + & + constitutive_dislotwin_Cslip_3333(p,q,r,s,i)*& + lattice_Qtwin(l,p,sum(lattice_NslipSystem(1:j-1,myStructure))+k,myStructure)* & + lattice_Qtwin(m,q,sum(lattice_NslipSystem(1:j-1,myStructure))+k,myStructure)* & + lattice_Qtwin(n,r,sum(lattice_NslipSystem(1:j-1,myStructure))+k,myStructure)* & + lattice_Qtwin(o,s,sum(lattice_NslipSystem(1:j-1,myStructure))+k,myStructure) + enddo ; enddo ; enddo ; enddo ; enddo ; enddo ; enddo ; enddo constitutive_dislotwin_Ctwin_66(:,:,k,i) = math_Mandel3333to66(constitutive_dislotwin_Ctwin_3333(:,:,:,:,k,i)) - enddo + enddo enddo !* Burgers vector, dislocation velocity prefactor, mean free path prefactor and minimum dipole distance for each slip system @@ -514,7 +526,7 @@ do i = 1,maxNinstance constitutive_dislotwin_burgersPerSlipSystem(s,i) = constitutive_dislotwin_burgersPerSlipFamily(f,i) constitutive_dislotwin_QedgePerSlipSystem(s,i) = constitutive_dislotwin_QedgePerSlipFamily(f,i) constitutive_dislotwin_v0PerSlipSystem(s,i) = constitutive_dislotwin_v0PerSlipFamily(f,i) - constitutive_dislotwin_CLambdaSlipPerSlipSystem(s,i) = constitutive_dislotwin_CLambdaSlipPerSlipFamily(f,i) + constitutive_dislotwin_CLambdaSlipPerSlipSystem(s,i) = constitutive_dislotwin_CLambdaSlipPerSlipFamily(f,i) enddo !* Burgers vector, nucleation rate prefactor and twin size for each twin system @@ -531,7 +543,8 @@ do i = 1,maxNinstance constitutive_dislotwin_interactionMatrixSlipSlip(s1,s2,i) = & constitutive_dislotwin_interactionSlipSlip(lattice_interactionSlipSlip(constitutive_dislotwin_slipSystemLattice(s1,i), & constitutive_dislotwin_slipSystemLattice(s2,i), & - myStructure),i) + myStructure),i) + enddo; enddo do s1 = 1,constitutive_dislotwin_totalNslip(i) @@ -576,7 +589,8 @@ function constitutive_dislotwin_stateInit(myInstance) !********************************************************************* !* initial microstructural state * !********************************************************************* -use prec, only: pReal,pInt +use prec, only: pReal,pInt + use math, only: pi use lattice, only: lattice_maxNslipFamily,lattice_maxNtwinFamily implicit none @@ -586,24 +600,31 @@ integer(pInt) :: myInstance real(pReal), dimension(constitutive_dislotwin_sizeState(myInstance)) :: constitutive_dislotwin_stateInit !* Local variables integer(pInt) s0,s1,s,t,f,ns,nt -real(pReal), dimension(constitutive_dislotwin_totalNslip(myInstance)) :: rhoEdge0,rhoEdgeDip0,invLambdaSlip0,MeanFreePathSlip0,tauSlipThreshold0 -real(pReal), dimension(constitutive_dislotwin_totalNtwin(myInstance)) :: MeanFreePathTwin0,TwinVolume0 +real(pReal), dimension(constitutive_dislotwin_totalNslip(myInstance)) :: rhoEdge0, & + rhoEdgeDip0, & + invLambdaSlip0, & + MeanFreePathSlip0, & + tauSlipThreshold0 -ns = constitutive_dislotwin_totalNslip(myInstance) +real(pReal), dimension(constitutive_dislotwin_totalNtwin(myInstance)) :: MeanFreePathTwin0,TwinVolume0 + +ns = constitutive_dislotwin_totalNslip(myInstance) nt = constitutive_dislotwin_totalNtwin(myInstance) constitutive_dislotwin_stateInit = 0.0_pReal - -!* Initialize basic slip state variables -s1 = 0_pInt -do f = 1,lattice_maxNslipFamily - s0 = s1 + 1_pInt - s1 = s0 + constitutive_dislotwin_Nslip(f,myInstance) - 1_pInt - do s = s0,s1 - rhoEdge0(s) = constitutive_dislotwin_rhoEdge0(f,myInstance) - rhoEdgeDip0(s) = constitutive_dislotwin_rhoEdgeDip0(f,myInstance) - enddo -enddo -constitutive_dislotwin_stateInit(1:ns) = rhoEdge0 + +!* Initialize basic slip state variables + +s1 = 0_pInt +do f = 1,lattice_maxNslipFamily + s0 = s1 + 1_pInt + s1 = s0 + constitutive_dislotwin_Nslip(f,myInstance) - 1_pInt + do s = s0,s1 + rhoEdge0(s) = constitutive_dislotwin_rhoEdge0(f,myInstance) + rhoEdgeDip0(s) = constitutive_dislotwin_rhoEdgeDip0(f,myInstance) + enddo +enddo + +constitutive_dislotwin_stateInit(1:ns) = rhoEdge0 constitutive_dislotwin_stateInit(ns+1:2*ns) = rhoEdgeDip0 !* Initialize dependent slip microstructural variables @@ -614,14 +635,16 @@ constitutive_dislotwin_stateInit(2*ns+nt+1:3*ns+nt) = invLambdaSlip0 forall (s = 1:ns) & MeanFreePathSlip0(s) = & -constitutive_dislotwin_GrainSize(myInstance)/(1.0_pReal+invLambdaSlip0(s)*constitutive_dislotwin_GrainSize(myInstance)) +constitutive_dislotwin_GrainSize(myInstance)/(1.0_pReal+invLambdaSlip0(s)*constitutive_dislotwin_GrainSize(myInstance)) + constitutive_dislotwin_stateInit(4*ns+2*nt+1:5*ns+2*nt) = MeanFreePathSlip0 forall (s = 1:ns) & tauSlipThreshold0(s) = & constitutive_dislotwin_Gmod(myInstance)*constitutive_dislotwin_burgersPerSlipSystem(s,myInstance)* & sqrt(dot_product(rhoEdge0,constitutive_dislotwin_interactionMatrixSlipSlip(1:ns,s,myInstance))) -constitutive_dislotwin_stateInit(5*ns+3*nt+1:6*ns+3*nt) = tauSlipThreshold0 +constitutive_dislotwin_stateInit(5*ns+3*nt+1:6*ns+3*nt) = tauSlipThreshold0 + !* Initialize dependent twin microstructural variables forall (t = 1:nt) & @@ -632,16 +655,17 @@ forall (t = 1:nt) & TwinVolume0(t) = & (pi/6.0_pReal)*constitutive_dislotwin_twinsizePerTwinSystem(t,myInstance)*MeanFreePathTwin0(t)**(2.0_pReal) constitutive_dislotwin_stateInit(6*ns+4*nt+1:6*ns+5*nt) = TwinVolume0 - -!write(6,*) '#STATEINIT#' -!write(6,*) -!write(6,'(a,/,4(3(f30.20,x)/))') 'RhoEdge',rhoEdge0 -!write(6,'(a,/,4(3(f30.20,x)/))') 'RhoEdgedip',rhoEdgeDip0 -!write(6,'(a,/,4(3(f30.20,x)/))') 'invLambdaSlip',invLambdaSlip0 -!write(6,'(a,/,4(3(f30.20,x)/))') 'MeanFreePathSlip',MeanFreePathSlip0 -!write(6,'(a,/,4(3(f30.20,x)/))') 'tauSlipThreshold', tauSlipThreshold0 -!write(6,'(a,/,4(3(f30.20,x)/))') 'MeanFreePathTwin', MeanFreePathTwin0 -!write(6,'(a,/,4(3(f30.20,x)/))') 'TwinVolume', TwinVolume0 + + +!write(6,*) '#STATEINIT#' +!write(6,*) +!write(6,'(a,/,4(3(f30.20,x)/))') 'RhoEdge',rhoEdge0 +!write(6,'(a,/,4(3(f30.20,x)/))') 'RhoEdgedip',rhoEdgeDip0 +!write(6,'(a,/,4(3(f30.20,x)/))') 'invLambdaSlip',invLambdaSlip0 +!write(6,'(a,/,4(3(f30.20,x)/))') 'MeanFreePathSlip',MeanFreePathSlip0 +!write(6,'(a,/,4(3(f30.20,x)/))') 'tauSlipThreshold', tauSlipThreshold0 +!write(6,'(a,/,4(3(f30.20,x)/))') 'MeanFreePathTwin', MeanFreePathTwin0 +!write(6,'(a,/,4(3(f30.20,x)/))') 'TwinVolume', TwinVolume0 return end function @@ -690,13 +714,15 @@ myInstance = phase_constitutionInstance(material_phase(g,ip,el)) ns = constitutive_dislotwin_totalNslip(myInstance) nt = constitutive_dislotwin_totalNtwin(myInstance) -!* Total twin volume fraction -sumf = sum(state(g,ip,el)%p((2*ns+1):(2*ns+nt))) ! safe for nt == 0 +!* Total twin volume fraction + +sumf = sum(state(g,ip,el)%p((2*ns+1):(2*ns+nt))) ! safe for nt == 0 !* Homogenized elasticity matrix constitutive_dislotwin_homogenizedC = (1.0_pReal-sumf)*constitutive_dislotwin_Cslip_66(:,:,myInstance) do i=1,nt - constitutive_dislotwin_homogenizedC = & + constitutive_dislotwin_homogenizedC = & + constitutive_dislotwin_homogenizedC + state(g,ip,el)%p(2*ns+i)*constitutive_dislotwin_Ctwin_66(:,:,i,myInstance) enddo @@ -735,7 +761,8 @@ myInstance = phase_constitutionInstance(material_phase(g,ip,el)) myStructure = constitutive_dislotwin_structure(myInstance) ns = constitutive_dislotwin_totalNslip(myInstance) nt = constitutive_dislotwin_totalNtwin(myInstance) -!* State: 1 : ns rho_edge +!* State: 1 : ns rho_edge + !* State: ns+1 : 2*ns rho_dipole !* State: 2*ns+1 : 2*ns+nt f !* State: 2*ns+nt+1 : 3*ns+nt 1/lambda_slip @@ -745,83 +772,93 @@ nt = constitutive_dislotwin_totalNtwin(myInstance) !* State: 5*ns+2*nt+1 : 5*ns+3*nt mfp_twin !* State: 5*ns+3*nt+1 : 6*ns+3*nt threshold_stress_slip !* State: 6*ns+3*nt+1 : 6*ns+4*nt threshold_stress_twin -!* State: 6*ns+4*nt+1 : 6*ns+5*nt twin volume +!* State: 6*ns+4*nt+1 : 6*ns+5*nt twin volume -!* Total twin volume fraction -sumf = sum(state(g,ip,el)%p((2*ns+1):(2*ns+nt))) ! safe for nt == 0 - -!* Stacking fault energy -sfe = 0.0002_pReal*Temperature-0.0396_pReal + +!* Total twin volume fraction + +sumf = sum(state(g,ip,el)%p((2*ns+1):(2*ns+nt))) ! safe for nt == 0 + + +!* Stacking fault energy + +sfe = 0.0002_pReal*Temperature-0.0396_pReal !* rescaled twin volume fraction for topology forall (t = 1:nt) & -fOverStacksize(t) = & -state(g,ip,el)%p(2*ns+t)/constitutive_dislotwin_twinsizePerTwinSystem(t,myInstance) + fOverStacksize(t) = & + state(g,ip,el)%p(2*ns+t)/constitutive_dislotwin_twinsizePerTwinSystem(t,myInstance) !* 1/mean free distance between 2 forest dislocations seen by a moving dislocation forall (s = 1:ns) & -state(g,ip,el)%p(2*ns+nt+s) = & -sqrt(dot_product((state(g,ip,el)%p(1:ns)+state(g,ip,el)%p(ns+1:2*ns)),constitutive_dislotwin_forestProjectionEdge(1:ns,s,myInstance)))/ & -constitutive_dislotwin_CLambdaSlipPerSlipSystem(s,myInstance) + state(g,ip,el)%p(2*ns+nt+s) = & + sqrt(dot_product((state(g,ip,el)%p(1:ns)+state(g,ip,el)%p(ns+1:2*ns)),& + constitutive_dislotwin_forestProjectionEdge(1:ns,s,myInstance)))/ & + constitutive_dislotwin_CLambdaSlipPerSlipSystem(s,myInstance) !* 1/mean free distance between 2 twin stacks from different systems seen by a moving dislocation !$OMP CRITICAL (evilmatmul) state(g,ip,el)%p((3*ns+nt+1):(4*ns+nt)) = 0.0_pReal -if (nt > 0_pInt) state(g,ip,el)%p((3*ns+nt+1):(4*ns+nt)) = & -matmul(constitutive_dislotwin_interactionMatrixSlipTwin(1:ns,1:nt,myInstance),fOverStacksize(1:nt))/(1.0_pReal-sumf) +if (nt > 0_pInt) & + state(g,ip,el)%p((3*ns+nt+1):(4*ns+nt)) = & + matmul(constitutive_dislotwin_interactionMatrixSlipTwin(1:ns,1:nt,myInstance),fOverStacksize(1:nt))/(1.0_pReal-sumf) !$OMP END CRITICAL (evilmatmul) !* 1/mean free distance between 2 twin stacks from different systems seen by a growing twin !$OMP CRITICAL (evilmatmul) -if (nt > 0_pInt) state(g,ip,el)%p((4*ns+nt+1):(4*ns+2*nt)) = & -matmul(constitutive_dislotwin_interactionMatrixTwinTwin(1:nt,1:nt,myInstance),fOverStacksize(1:nt))/(1.0_pReal-sumf) +if (nt > 0_pInt) & + state(g,ip,el)%p((4*ns+nt+1):(4*ns+2*nt)) = & + matmul(constitutive_dislotwin_interactionMatrixTwinTwin(1:nt,1:nt,myInstance),fOverStacksize(1:nt))/(1.0_pReal-sumf) !$OMP END CRITICAL (evilmatmul) !* mean free path between 2 obstacles seen by a moving dislocation do s = 1,ns if (nt > 0_pInt) then - state(g,ip,el)%p(4*ns+2*nt+s) = & - constitutive_dislotwin_GrainSize(myInstance)/(1.0_pReal+constitutive_dislotwin_GrainSize(myInstance)*& - (state(g,ip,el)%p(2*ns+nt+s)+state(g,ip,el)%p(3*ns+nt+s))) + state(g,ip,el)%p(4*ns+2*nt+s) = & + constitutive_dislotwin_GrainSize(myInstance)/(1.0_pReal+constitutive_dislotwin_GrainSize(myInstance)*& + (state(g,ip,el)%p(2*ns+nt+s)+state(g,ip,el)%p(3*ns+nt+s))) else - state(g,ip,el)%p(4*ns+s) = & - constitutive_dislotwin_GrainSize(myInstance)/& - (1.0_pReal+constitutive_dislotwin_GrainSize(myInstance)*(state(g,ip,el)%p(2*ns+s))) + state(g,ip,el)%p(4*ns+s) = & + constitutive_dislotwin_GrainSize(myInstance)/& + (1.0_pReal+constitutive_dislotwin_GrainSize(myInstance)*(state(g,ip,el)%p(2*ns+s))) endif enddo !* mean free path between 2 obstacles seen by a growing twin forall (t = 1:nt) & -state(g,ip,el)%p(5*ns+2*nt+t) = & -(constitutive_dislotwin_Cmfptwin(myInstance)*constitutive_dislotwin_GrainSize(myInstance))/& -(1.0_pReal+constitutive_dislotwin_GrainSize(myInstance)*state(g,ip,el)%p(4*ns+nt+t)) + state(g,ip,el)%p(5*ns+2*nt+t) = & + (constitutive_dislotwin_Cmfptwin(myInstance)*constitutive_dislotwin_GrainSize(myInstance))/& + (1.0_pReal+constitutive_dislotwin_GrainSize(myInstance)*state(g,ip,el)%p(4*ns+nt+t)) !* threshold stress for dislocation motion forall (s = 1:ns) & -state(g,ip,el)%p(5*ns+3*nt+s) = & -constitutive_dislotwin_Gmod(myInstance)*constitutive_dislotwin_burgersPerSlipSystem(s,myInstance)*& -sqrt(dot_product(state(g,ip,el)%p(1:ns),constitutive_dislotwin_interactionMatrixSlipSlip(1:ns,s,myInstance))) + state(g,ip,el)%p(5*ns+3*nt+s) = & + constitutive_dislotwin_Gmod(myInstance)*constitutive_dislotwin_burgersPerSlipSystem(s,myInstance)*& + sqrt(dot_product(state(g,ip,el)%p(1:ns),constitutive_dislotwin_interactionMatrixSlipSlip(1:ns,s,myInstance))) + !* threshold stress for growing twin forall (t = 1:nt) & -state(g,ip,el)%p(6*ns+3*nt+t) = & -constitutive_dislotwin_Cthresholdtwin(myInstance)*& -(sfe/(3.0_pReal*constitutive_dislotwin_burgersPerTwinSystem(t,myInstance))+& -(constitutive_dislotwin_burgersPerTwinSystem(t,myInstance)*constitutive_dislotwin_Gmod(myInstance))/& -state(g,ip,el)%p(5*ns+2*nt+t)) + state(g,ip,el)%p(6*ns+3*nt+t) = & + constitutive_dislotwin_Cthresholdtwin(myInstance)*& + (sfe/(3.0_pReal*constitutive_dislotwin_burgersPerTwinSystem(t,myInstance))+& + (constitutive_dislotwin_burgersPerTwinSystem(t,myInstance)*constitutive_dislotwin_Gmod(myInstance))/& + state(g,ip,el)%p(5*ns+2*nt+t)) !* final twin volume after growth forall (t = 1:nt) & -state(g,ip,el)%p(6*ns+4*nt+t) = & -(pi/6.0_pReal)*constitutive_dislotwin_twinsizePerTwinSystem(t,myInstance)*state(g,ip,el)%p(5*ns+2*nt+t)**(2.0_pReal) - -!if ((ip==1).and.(el==1)) then -! write(6,*) '#MICROSTRUCTURE#' -! write(6,*) -! write(6,'(a,/,4(3(f10.4,x)/))') 'rhoEdge',state(g,ip,el)%p(1:ns)/1e9 -! write(6,'(a,/,4(3(f10.4,x)/))') 'rhoEdgeDip',state(g,ip,el)%p(ns+1:2*ns)/1e9 -! write(6,'(a,/,4(3(f10.4,x)/))') 'Fraction',state(g,ip,el)%p(2*ns+1:2*ns+nt) -!endif + state(g,ip,el)%p(6*ns+4*nt+t) = & + (pi/6.0_pReal)*constitutive_dislotwin_twinsizePerTwinSystem(t,myInstance)*state(g,ip,el)%p(5*ns+2*nt+t)**(2.0_pReal) + + +!if ((ip==1).and.(el==1)) then +! write(6,*) '#MICROSTRUCTURE#' +! write(6,*) +! write(6,'(a,/,4(3(f10.4,x)/))') 'rhoEdge',state(g,ip,el)%p(1:ns)/1e9 +! write(6,'(a,/,4(3(f10.4,x)/))') 'rhoEdgeDip',state(g,ip,el)%p(ns+1:2*ns)/1e9 +! write(6,'(a,/,4(3(f10.4,x)/))') 'Fraction',state(g,ip,el)%p(2*ns+1:2*ns+nt) +!endif + return end subroutine @@ -846,7 +883,7 @@ use math, only: math_Plain3333to99 use mesh, only: mesh_NcpElems,mesh_maxNips use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance use lattice, only: lattice_Sslip,lattice_Sslip_v,lattice_Stwin,lattice_Stwin_v,lattice_maxNslipFamily,lattice_maxNtwinFamily, & - lattice_NslipSystem,lattice_NtwinSystem,lattice_shearTwin + lattice_NslipSystem,lattice_NtwinSystem,lattice_shearTwin implicit none !* Input-Output variables @@ -871,8 +908,9 @@ myStructure = constitutive_dislotwin_structure(myInstance) ns = constitutive_dislotwin_totalNslip(myInstance) nt = constitutive_dislotwin_totalNtwin(myInstance) -!* Total twin volume fraction -sumf = sum(state(g,ip,el)%p((2*ns+1):(2*ns+nt))) ! safe for nt == 0 +!* Total twin volume fraction + +sumf = sum(state(g,ip,el)%p((2*ns+1):(2*ns+nt))) ! safe for nt == 0 Lp = 0.0_pReal dLp_dTstar3333 = 0.0_pReal @@ -887,37 +925,38 @@ do f = 1,lattice_maxNslipFamily ! loop over all do i = 1,constitutive_dislotwin_Nslip(f,myInstance) ! process each (active) slip system in family j = j+1_pInt - !* Calculation of Lp - - !* Resolved shear stress on slip system - tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure)) - !* Stress ratios - StressRatio_p = (abs(tau_slip(j))/state(g,ip,el)%p(5*ns+3*nt+j))**constitutive_dislotwin_p(myInstance) - StressRatio_pminus1 = (abs(tau_slip(j))/state(g,ip,el)%p(5*ns+3*nt+j))**(constitutive_dislotwin_p(myInstance)-1.0_pReal) - !* Boltzmann ratio - BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(f,myInstance)/(kB*Temperature) - !* Initial shear rates - DotGamma0 = & - state(g,ip,el)%p(j)*constitutive_dislotwin_burgersPerSlipSystem(f,myInstance)*constitutive_dislotwin_v0PerSlipSystem(f,myInstance) - - !* Shear rates due to slip - gdot_slip(j) = DotGamma0*exp(-BoltzmannRatio*(1-StressRatio_p)**constitutive_dislotwin_q(myInstance))*sign(1.0_pReal,tau_slip(j)) - - !* Derivatives of shear rates - dgdot_dtauslip(j) = & - ((gdot_slip(j)*BoltzmannRatio*& - constitutive_dislotwin_p(myInstance)*constitutive_dislotwin_q(myInstance))/state(g,ip,el)%p(5*ns+3*nt+j))*& - StressRatio_pminus1*(1-StressRatio_p)**(constitutive_dislotwin_q(myInstance)-1.0_pReal) - + !* Calculation of Lp + !* Resolved shear stress on slip system + tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure)) + + !* Stress ratios + StressRatio_p = (abs(tau_slip(j))/state(g,ip,el)%p(5*ns+3*nt+j))**constitutive_dislotwin_p(myInstance) + StressRatio_pminus1 = (abs(tau_slip(j))/state(g,ip,el)%p(5*ns+3*nt+j))**(constitutive_dislotwin_p(myInstance)-1.0_pReal) + !* Boltzmann ratio + BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(f,myInstance)/(kB*Temperature) + !* Initial shear rates + DotGamma0 = & + state(g,ip,el)%p(j)*constitutive_dislotwin_burgersPerSlipSystem(f,myInstance)*& + constitutive_dislotwin_v0PerSlipSystem(f,myInstance) + + !* Shear rates due to slip + gdot_slip(j) = DotGamma0*exp(-BoltzmannRatio*(1-StressRatio_p)**constitutive_dislotwin_q(myInstance))*sign(1.0_pReal,tau_slip(j)) + + !* Derivatives of shear rates + dgdot_dtauslip(j) = & + ((gdot_slip(j)*BoltzmannRatio*& + constitutive_dislotwin_p(myInstance)*constitutive_dislotwin_q(myInstance))/state(g,ip,el)%p(5*ns+3*nt+j))*& + StressRatio_pminus1*(1-StressRatio_p)**(constitutive_dislotwin_q(myInstance)-1.0_pReal) + !* Plastic velocity gradient for dislocation glide Lp = Lp + (1.0_pReal - sumf)*gdot_slip(j)*lattice_Sslip(:,:,index_myFamily+i,myStructure) - - !* Calculation of the tangent of Lp + + !* Calculation of the tangent of Lp forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLp_dTstar3333(k,l,m,n) = & - dLp_dTstar3333(k,l,m,n) + dgdot_dtauslip(j)*& - lattice_Sslip(k,l,index_myFamily+i,myStructure)*& - lattice_Sslip(m,n,index_myFamily+i,myStructure) + dLp_dTstar3333(k,l,m,n) = & + dLp_dTstar3333(k,l,m,n) + dgdot_dtauslip(j)*& + lattice_Sslip(k,l,index_myFamily+i,myStructure)*& + lattice_Sslip(m,n,index_myFamily+i,myStructure) enddo enddo @@ -930,42 +969,44 @@ do f = 1,lattice_maxNtwinFamily ! loop over all do i = 1,constitutive_dislotwin_Ntwin(f,myInstance) ! process each (active) slip system in family j = j+1_pInt - !* Calculation of Lp - - !* Resolved shear stress on twin system - tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,myStructure)) - !* Stress ratios - StressRatio_r = (state(g,ip,el)%p(6*ns+3*nt+j)/tau_twin(j))**constitutive_dislotwin_r(myInstance) - - !* Shear rates and their derivatives due to twin - if ( tau_twin(j) > 0.0_pReal ) then - gdot_twin(j) = & - (constitutive_dislotwin_MaxTwinFraction(myInstance)-sumf)*lattice_shearTwin(index_myFamily+i,myStructure)*& - state(g,ip,el)%p(6*ns+4*nt+j)*constitutive_dislotwin_Ndot0PerTwinSystem(f,myInstance)*exp(-StressRatio_r) - - dgdot_dtautwin(j) = ((gdot_twin(j)*constitutive_dislotwin_r(myInstance))/tau_twin(j))*StressRatio_r - endif - + !* Calculation of Lp + !* Resolved shear stress on twin system + + tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,myStructure)) + !* Stress ratios + StressRatio_r = (state(g,ip,el)%p(6*ns+3*nt+j)/tau_twin(j))**constitutive_dislotwin_r(myInstance) + + !* Shear rates and their derivatives due to twin + if ( tau_twin(j) > 0.0_pReal ) then + gdot_twin(j) = & + (constitutive_dislotwin_MaxTwinFraction(myInstance)-sumf)*lattice_shearTwin(index_myFamily+i,myStructure)*& + state(g,ip,el)%p(6*ns+4*nt+j)*constitutive_dislotwin_Ndot0PerTwinSystem(f,myInstance)*exp(-StressRatio_r) + dgdot_dtautwin(j) = ((gdot_twin(j)*constitutive_dislotwin_r(myInstance))/tau_twin(j))*StressRatio_r + endif + + !* Plastic velocity gradient for mechanical twinning Lp = Lp + gdot_twin(j)*lattice_Stwin(:,:,index_myFamily+i,myStructure) !* Calculation of the tangent of Lp forall (k=1:3,l=1:3,m=1:3,n=1:3) & - dLp_dTstar3333(k,l,m,n) = & - dLp_dTstar3333(k,l,m,n) + dgdot_dtautwin(j)*& - lattice_Stwin(k,l,index_myFamily+i,myStructure)*& - lattice_Stwin(m,n,index_myFamily+i,myStructure) + dLp_dTstar3333(k,l,m,n) = & + dLp_dTstar3333(k,l,m,n) + dgdot_dtautwin(j)*& + lattice_Stwin(k,l,index_myFamily+i,myStructure)*& + lattice_Stwin(m,n,index_myFamily+i,myStructure) enddo enddo -dLp_dTstar = math_Plain3333to99(dLp_dTstar3333) - -!if ((ip==1).and.(el==1)) then -! write(6,*) '#MICROSTRUCTURE#' -! write(6,*) -! write(6,'(a,/,4(3(f10.4,x)/))') 'rhoEdge',state(g,ip,el)%p(1:12) -! write(6,'(a,/,4(3(f10.4,x)/))') 'rhoEdgeDip',state(g,ip,el)%p(13:24) -!endif +dLp_dTstar = math_Plain3333to99(dLp_dTstar3333) + + +!if ((ip==1).and.(el==1)) then +! write(6,*) '#MICROSTRUCTURE#' +! write(6,*) +! write(6,'(a,/,4(3(f10.4,x)/))') 'rhoEdge',state(g,ip,el)%p(1:12) +! write(6,'(a,/,4(3(f10.4,x)/))') 'rhoEdgeDip',state(g,ip,el)%p(13:24) +!endif + return end subroutine @@ -984,7 +1025,8 @@ function constitutive_dislotwin_dotState(Tstar_v,Temperature,state,g,ip,el) !* OUTPUT: * !* - constitutive_dotState : evolution of state variable * !********************************************************************* -use prec, only: pReal,pInt,p_vec +use prec, only: pReal,pInt,p_vec + use math, only: pi use mesh, only: mesh_NcpElems,mesh_maxNips use material, only: homogenization_maxNgrains,material_phase, phase_constitutionInstance @@ -995,15 +1037,18 @@ implicit none !* Input-Output variables integer(pInt), intent(in) :: g,ip,el real(pReal), intent(in) :: Temperature -real(pReal), dimension(6), intent(in) :: Tstar_v +real(pReal), dimension(6), intent(in) :: Tstar_v + type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: state real(pReal), dimension(constitutive_dislotwin_sizeDotState(phase_constitutionInstance(material_phase(g,ip,el)))) :: & constitutive_dislotwin_dotState !* Local variables integer(pInt) MyInstance,MyStructure,ns,nt,f,i,j,k,index_myFamily -real(pReal) sumf,StressRatio_p,StressRatio_pminus1,BoltzmannRatio,DotGamma0,EdgeDipMinDistance,AtomicVolume,VacancyDiffusion,StressRatio_r +real(pReal) sumf,StressRatio_p,StressRatio_pminus1,BoltzmannRatio,DotGamma0,& + EdgeDipMinDistance,AtomicVolume,VacancyDiffusion,StressRatio_r real(pReal), dimension(constitutive_dislotwin_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: & -gdot_slip,tau_slip,DotRhoMultiplication,EdgeDipDistance,DotRhoEdgeEdgeAnnihilation,DotRhoEdgeDipAnnihilation,& +gdot_slip,tau_slip,DotRhoMultiplication,EdgeDipDistance,DotRhoEdgeEdgeAnnihilation,DotRhoEdgeDipAnnihilation,& + ClimbVelocity,DotRhoEdgeDipClimb,DotRhoDipFormation real(pReal), dimension(constitutive_dislotwin_totalNtwin(phase_constitutionInstance(material_phase(g,ip,el)))) :: gdot_twin,tau_twin @@ -1013,8 +1058,9 @@ MyStructure = constitutive_dislotwin_structure(myInstance) ns = constitutive_dislotwin_totalNslip(myInstance) nt = constitutive_dislotwin_totalNtwin(myInstance) -!* Total twin volume fraction -sumf = sum(state(g,ip,el)%p((2*ns+1):(2*ns+nt))) ! safe for nt == 0 +!* Total twin volume fraction + +sumf = sum(state(g,ip,el)%p((2*ns+1):(2*ns+nt))) ! safe for nt == 0 constitutive_dislotwin_dotState = 0.0_pReal @@ -1025,66 +1071,73 @@ do f = 1,lattice_maxNslipFamily ! loop over all index_myFamily = sum(lattice_NslipSystem(1:f-1,MyStructure)) ! at which index starts my family do i = 1,constitutive_dislotwin_Nslip(f,myInstance) ! process each (active) slip system in family j = j+1_pInt - - !* Resolved shear stress on slip system + + + !* Resolved shear stress on slip system tau_slip(j) = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure)) - !* Stress ratios - StressRatio_p = (abs(tau_slip(j))/state(g,ip,el)%p(5*ns+3*nt+j))**constitutive_dislotwin_p(myInstance) - StressRatio_pminus1 = (abs(tau_slip(j))/state(g,ip,el)%p(5*ns+3*nt+j))**(constitutive_dislotwin_p(myInstance)-1.0_pReal) - !* Boltzmann ratio - BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(f,myInstance)/(kB*Temperature) - !* Initial shear rates - DotGamma0 = & - state(g,ip,el)%p(j)*constitutive_dislotwin_burgersPerSlipSystem(f,myInstance)*constitutive_dislotwin_v0PerSlipSystem(f,myInstance) + !* Stress ratios + StressRatio_p = (abs(tau_slip(j))/state(g,ip,el)%p(5*ns+3*nt+j))**constitutive_dislotwin_p(myInstance) + StressRatio_pminus1 = (abs(tau_slip(j))/state(g,ip,el)%p(5*ns+3*nt+j))**(constitutive_dislotwin_p(myInstance)-1.0_pReal) + !* Boltzmann ratio + BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(f,myInstance)/(kB*Temperature) + !* Initial shear rates + DotGamma0 = & + state(g,ip,el)%p(j)*constitutive_dislotwin_burgersPerSlipSystem(f,myInstance)*& + constitutive_dislotwin_v0PerSlipSystem(f,myInstance) - !* Shear rates due to slip + !* Shear rates due to slip gdot_slip(j) = DotGamma0*exp(-BoltzmannRatio*(1-StressRatio_p)**constitutive_dislotwin_q(myInstance))*sign(1.0_pReal,tau_slip(j)) - + + !* Multiplication - DotRhoMultiplication(j) = abs(gdot_slip(j))/(constitutive_dislotwin_burgersPerSlipSystem(f,myInstance)*state(g,ip,el)%p(4*ns+2*nt+j)) - - !* Dipole formation - if (tau_slip(j) == 0.0_pReal) then - DotRhoDipFormation(j) = 0.0_pReal - else - EdgeDipDistance(j) = & - (3.0_pReal*constitutive_dislotwin_Gmod(myInstance)*constitutive_dislotwin_burgersPerSlipSystem(f,myInstance))/& - (16.0_pReal*pi*abs(tau_slip(j))) - DotRhoDipFormation(j) = & - ((2.0_pReal*EdgeDipDistance(j))/constitutive_dislotwin_burgersPerSlipSystem(f,myInstance))*& - state(g,ip,el)%p(j)*abs(gdot_slip(j)) - endif - - !* Spontaneous annihilation of 2 single edge dislocations - EdgeDipMinDistance = & - constitutive_dislotwin_CEdgeDipMinDistance(myInstance)*constitutive_dislotwin_burgersPerSlipSystem(f,myInstance) - DotRhoEdgeEdgeAnnihilation(j) = & - ((2.0_pReal*EdgeDipMinDistance)/constitutive_dislotwin_burgersPerSlipSystem(f,myInstance))*& - state(g,ip,el)%p(j)*abs(gdot_slip(j)) - - !* Spontaneous annihilation of a single edge dislocation with a dipole constituent - DotRhoEdgeDipAnnihilation(j) = & - ((2.0_pReal*EdgeDipMinDistance)/constitutive_dislotwin_burgersPerSlipSystem(f,myInstance))*& - state(g,ip,el)%p(ns+j)*abs(gdot_slip(j)) - - !* Dislocation dipole climb - AtomicVolume = & - constitutive_dislotwin_CAtomicVolume(myInstance)*constitutive_dislotwin_burgersPerSlipSystem(f,myInstance)**(3.0_pReal) - VacancyDiffusion = & - constitutive_dislotwin_D0(myInstance)*exp(-constitutive_dislotwin_Qsd(myInstance)/(kB*Temperature)) - ClimbVelocity(j) = & - ((3.0_pReal*constitutive_dislotwin_Gmod(myInstance)*VacancyDiffusion*AtomicVolume)/(2.0_pReal*pi*kB*Temperature))*& - (1/(EdgeDipDistance(j)+EdgeDipMinDistance)) - DotRhoEdgeDipClimb(j) = & - (4.0_pReal*ClimbVelocity(j)*state(g,ip,el)%p(ns+j))/(EdgeDipDistance(j)+EdgeDipMinDistance) - - !* Edge dislocation density rate of change - constitutive_dislotwin_dotState(j) = & - DotRhoMultiplication(j)-DotRhoDipFormation(j)-DotRhoEdgeEdgeAnnihilation(j) - - !* Edge dislocation dipole density rate of change - constitutive_dislotwin_dotState(ns+j) = & - DotRhoDipFormation(j)-DotRhoEdgeDipAnnihilation(j)-DotRhoEdgeDipClimb(j) + DotRhoMultiplication(j) = abs(gdot_slip(j))/& + (constitutive_dislotwin_burgersPerSlipSystem(f,myInstance)*state(g,ip,el)%p(4*ns+2*nt+j)) + + !* Dipole formation + + if (tau_slip(j) == 0.0_pReal) then + DotRhoDipFormation(j) = 0.0_pReal + else + EdgeDipDistance(j) = & + (3.0_pReal*constitutive_dislotwin_Gmod(myInstance)*constitutive_dislotwin_burgersPerSlipSystem(f,myInstance))/& + (16.0_pReal*pi*abs(tau_slip(j))) + DotRhoDipFormation(j) = & + ((2.0_pReal*EdgeDipDistance(j))/constitutive_dislotwin_burgersPerSlipSystem(f,myInstance))*& + state(g,ip,el)%p(j)*abs(gdot_slip(j)) + endif + + !* Spontaneous annihilation of 2 single edge dislocations + EdgeDipMinDistance = & + constitutive_dislotwin_CEdgeDipMinDistance(myInstance)*constitutive_dislotwin_burgersPerSlipSystem(f,myInstance) + DotRhoEdgeEdgeAnnihilation(j) = & + ((2.0_pReal*EdgeDipMinDistance)/constitutive_dislotwin_burgersPerSlipSystem(f,myInstance))*& + state(g,ip,el)%p(j)*abs(gdot_slip(j)) + + !* Spontaneous annihilation of a single edge dislocation with a dipole constituent + DotRhoEdgeDipAnnihilation(j) = & + ((2.0_pReal*EdgeDipMinDistance)/constitutive_dislotwin_burgersPerSlipSystem(f,myInstance))*& + state(g,ip,el)%p(ns+j)*abs(gdot_slip(j)) + + !* Dislocation dipole climb + AtomicVolume = & + constitutive_dislotwin_CAtomicVolume(myInstance)*constitutive_dislotwin_burgersPerSlipSystem(f,myInstance)**(3.0_pReal) + + VacancyDiffusion = & + constitutive_dislotwin_D0(myInstance)*exp(-constitutive_dislotwin_Qsd(myInstance)/(kB*Temperature)) + ClimbVelocity(j) = & + ((3.0_pReal*constitutive_dislotwin_Gmod(myInstance)*VacancyDiffusion*AtomicVolume)/(2.0_pReal*pi*kB*Temperature))*& + (1/(EdgeDipDistance(j)+EdgeDipMinDistance)) + DotRhoEdgeDipClimb(j) = & + (4.0_pReal*ClimbVelocity(j)*state(g,ip,el)%p(ns+j))/(EdgeDipDistance(j)+EdgeDipMinDistance) + + !* Edge dislocation density rate of change + constitutive_dislotwin_dotState(j) = & + DotRhoMultiplication(j)-DotRhoDipFormation(j)-DotRhoEdgeEdgeAnnihilation(j) + + + !* Edge dislocation dipole density rate of change + constitutive_dislotwin_dotState(ns+j) = & + DotRhoDipFormation(j)-DotRhoEdgeDipAnnihilation(j)-DotRhoEdgeDipClimb(j) enddo enddo @@ -1096,35 +1149,39 @@ do f = 1,lattice_maxNtwinFamily ! loop over all index_myFamily = sum(lattice_NtwinSystem(1:f-1,MyStructure)) ! at which index starts my family do i = 1,constitutive_dislotwin_Ntwin(f,myInstance) ! process each (active) twin system in family j = j+1_pInt - - !* Resolved shear stress on twin system - tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,myStructure)) - !* Stress ratios - StressRatio_r = (state(g,ip,el)%p(6*ns+3*nt+j)/tau_twin(j))**constitutive_dislotwin_r(myInstance) - - !* Shear rates and their derivatives due to twin - if ( tau_twin(j) > 0.0_pReal ) then - gdot_twin(j) = & - (constitutive_dislotwin_MaxTwinFraction(myInstance)-sumf)*& - state(g,ip,el)%p(6*ns+4*nt+j)*constitutive_dislotwin_Ndot0PerTwinSystem(f,myInstance)*exp(-StressRatio_r) - endif + + !* Resolved shear stress on twin system + tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,myStructure)) + !* Stress ratios + StressRatio_r = (state(g,ip,el)%p(6*ns+3*nt+j)/tau_twin(j))**constitutive_dislotwin_r(myInstance) + + !* Shear rates and their derivatives due to twin + if ( tau_twin(j) > 0.0_pReal ) then + gdot_twin(j) = & + (constitutive_dislotwin_MaxTwinFraction(myInstance)-sumf)*& + state(g,ip,el)%p(6*ns+4*nt+j)*constitutive_dislotwin_Ndot0PerTwinSystem(f,myInstance)*exp(-StressRatio_r) + endif + enddo enddo - -!if ((ip==1).and.(el==1)) then -! 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)/))') 'RhoEdge',state(g,ip,el)%p(1:ns) -! write(6,'(a,/,4(3(f30.20,x)/))') 'Threshold Slip', state(g,ip,el)%p(5*ns+3*nt+1:6*ns+3*nt) -! write(6,'(a,/,4(3(f30.20,x)/))') 'Multiplication',DotRhoMultiplication -! write(6,'(a,/,4(3(f30.20,x)/))') 'DipFormation',DotRhoDipFormation -! write(6,'(a,/,4(3(f30.20,x)/))') 'SingleSingle',DotRhoEdgeEdgeAnnihilation -! write(6,'(a,/,4(3(f30.20,x)/))') 'SingleDipole',DotRhoEdgeDipAnnihilation -! write(6,'(a,/,4(3(f30.20,x)/))') 'DipClimb',DotRhoEdgeDipClimb -!endif - + + +!if ((ip==1).and.(el==1)) then +! 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)/))') 'RhoEdge',state(g,ip,el)%p(1:ns) +! write(6,'(a,/,4(3(f30.20,x)/))') 'Threshold Slip', state(g,ip,el)%p(5*ns+3*nt+1:6*ns+3*nt) +! write(6,'(a,/,4(3(f30.20,x)/))') 'Multiplication',DotRhoMultiplication +! write(6,'(a,/,4(3(f30.20,x)/))') 'DipFormation',DotRhoDipFormation +! write(6,'(a,/,4(3(f30.20,x)/))') 'SingleSingle',DotRhoEdgeEdgeAnnihilation +! write(6,'(a,/,4(3(f30.20,x)/))') 'SingleDipole',DotRhoEdgeDipAnnihilation +! write(6,'(a,/,4(3(f30.20,x)/))') 'DipClimb',DotRhoEdgeDipClimb +!endif + + + return end function @@ -1174,7 +1231,7 @@ use prec, only: pReal,pInt,p_vec use mesh, only: mesh_NcpElems,mesh_maxNips use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance,phase_Noutput use lattice, only: lattice_Sslip_v,lattice_Stwin_v,lattice_maxNslipFamily,lattice_maxNtwinFamily, & - lattice_NslipSystem,lattice_NtwinSystem,lattice_shearTwin + lattice_NslipSystem,lattice_NtwinSystem,lattice_shearTwin implicit none !* Definition of variables @@ -1206,39 +1263,36 @@ do o = 1,phase_Noutput(material_phase(g,ip,el)) case ('edge_density') constitutive_dislotwin_postResults(c+1:c+ns) = state(g,ip,el)%p(1:ns) c = c + ns - - case ('dipole_density') - constitutive_dislotwin_postResults(c+1:c+ns) = state(g,ip,el)%p(ns+1:2*ns) - c = c + ns - + case ('dipole_density') + constitutive_dislotwin_postResults(c+1:c+ns) = state(g,ip,el)%p(ns+1:2*ns) + c = c + ns case ('shear_rate_slip') j = 0_pInt do f = 1,lattice_maxNslipFamily ! loop over all slip families index_myFamily = sum(lattice_NslipSystem(1:f-1,myStructure)) ! at which index starts my family do i = 1,constitutive_dislotwin_Nslip(f,myInstance) ! process each (active) slip system in family j = j + 1_pInt - - !* Resolved shear stress on slip system - tau = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure)) - !* Stress ratios - StressRatio_p = (abs(tau)/state(g,ip,el)%p(5*ns+3*nt+j))**constitutive_dislotwin_p(myInstance) - StressRatio_pminus1 = (abs(tau)/state(g,ip,el)%p(5*ns+3*nt+j))**(constitutive_dislotwin_p(myInstance)-1.0_pReal) - !* Boltzmann ratio - BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(f,myInstance)/(kB*Temperature) - !* Initial shear rates - DotGamma0 = & - state(g,ip,el)%p(j)*constitutive_dislotwin_burgersPerSlipSystem(f,myInstance)*constitutive_dislotwin_v0PerSlipSystem(f,myInstance) - - !* Shear rates due to slip - constitutive_dislotwin_postResults(c+j) = & - DotGamma0*exp(-BoltzmannRatio*(1-StressRatio_p)**constitutive_dislotwin_q(myInstance))*sign(1.0_pReal,tau) - enddo ; enddo - c = c + ns + !* Resolved shear stress on slip system + tau = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure)) + !* Stress ratios + StressRatio_p = (abs(tau)/state(g,ip,el)%p(5*ns+3*nt+j))**constitutive_dislotwin_p(myInstance) + StressRatio_pminus1 = (abs(tau)/state(g,ip,el)%p(5*ns+3*nt+j))**(constitutive_dislotwin_p(myInstance)-1.0_pReal) + !* Boltzmann ratio + BoltzmannRatio = constitutive_dislotwin_QedgePerSlipSystem(f,myInstance)/(kB*Temperature) + !* Initial shear rates + DotGamma0 = & + state(g,ip,el)%p(j)*constitutive_dislotwin_burgersPerSlipSystem(f,myInstance)* & + constitutive_dislotwin_v0PerSlipSystem(f,myInstance) + + !* Shear rates due to slip + constitutive_dislotwin_postResults(c+j) = & + DotGamma0*exp(-BoltzmannRatio*(1-StressRatio_p)**constitutive_dislotwin_q(myInstance))*sign(1.0_pReal,tau) + enddo ; enddo + c = c + ns case ('mfp_slip') constitutive_dislotwin_postResults(c+1:c+ns) = state(g,ip,el)%p((4*ns+2*nt+1):(5*ns+2*nt)) c = c + ns - case ('resolved_stress_slip') j = 0_pInt do f = 1,lattice_maxNslipFamily ! loop over all slip families @@ -1246,64 +1300,58 @@ do o = 1,phase_Noutput(material_phase(g,ip,el)) do i = 1,constitutive_dislotwin_Nslip(f,myInstance) ! process each (active) slip system in family j = j + 1_pInt constitutive_dislotwin_postResults(c+j) = dot_product(Tstar_v,lattice_Sslip_v(:,index_myFamily+i,myStructure)) - enddo; enddo + enddo; enddo c = c + ns - case ('threshold_stress_slip') constitutive_dislotwin_postResults(c+1:c+ns) = state(g,ip,el)%p((5*ns+3*nt+1):(6*ns+3*nt)) c = c + ns - case ('twin_fraction') constitutive_dislotwin_postResults(c+1:c+nt) = state(g,ip,el)%p((2*ns+1):(2*ns+nt)) c = c + nt - case ('shear_rate_twin') - if (nt > 0_pInt) then - j = 0_pInt - do f = 1,lattice_maxNtwinFamily ! loop over all twin families - index_myFamily = sum(lattice_NtwinSystem(1:f-1,myStructure)) ! at which index starts my family - do i = 1,constitutive_dislotwin_Ntwin(f,myInstance) ! process each (active) twin system in family + if (nt > 0_pInt) then + j = 0_pInt + do f = 1,lattice_maxNtwinFamily ! loop over all twin families + index_myFamily = sum(lattice_NtwinSystem(1:f-1,myStructure)) ! at which index starts my family + do i = 1,constitutive_dislotwin_Ntwin(f,myInstance) ! process each (active) twin system in family j = j + 1_pInt - !* Resolved shear stress on twin system - tau = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,myStructure)) - !* Stress ratios - StressRatio_r = (state(g,ip,el)%p(6*ns+3*nt+j)/tau)**constitutive_dislotwin_r(myInstance) - - !* Shear rates and their derivatives due to twin - if ( tau > 0.0_pReal ) then - constitutive_dislotwin_postResults(c+j) = & - (constitutive_dislotwin_MaxTwinFraction(myInstance)-sumf)*& - state(g,ip,el)%p(6*ns+4*nt+j)*constitutive_dislotwin_Ndot0PerTwinSystem(f,myInstance)*exp(-StressRatio_r) - endif - enddo ; enddo - endif - c = c + nt + !* Resolved shear stress on twin system + tau = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,myStructure)) + !* Stress ratios + StressRatio_r = (state(g,ip,el)%p(6*ns+3*nt+j)/tau)**constitutive_dislotwin_r(myInstance) + !* Shear rates and their derivatives due to twin + if ( tau > 0.0_pReal ) then + constitutive_dislotwin_postResults(c+j) = & + (constitutive_dislotwin_MaxTwinFraction(myInstance)-sumf)*& + state(g,ip,el)%p(6*ns+4*nt+j)*constitutive_dislotwin_Ndot0PerTwinSystem(f,myInstance)*exp(-StressRatio_r) + endif + + enddo ; enddo + endif + c = c + nt case ('mfp_twin') constitutive_dislotwin_postResults(c+1:c+nt) = state(g,ip,el)%p((5*ns+2*nt+1):(5*ns+3*nt)) c = c + nt - case ('resolved_stress_twin') - if (nt > 0_pInt) then - j = 0_pInt - do f = 1,lattice_maxNtwinFamily ! loop over all slip families - index_myFamily = sum(lattice_NtwinSystem(1:f-1,myStructure)) ! at which index starts my family - do i = 1,constitutive_dislotwin_Ntwin(f,myInstance) ! process each (active) slip system in family + if (nt > 0_pInt) then + j = 0_pInt + do f = 1,lattice_maxNtwinFamily ! loop over all slip families + index_myFamily = sum(lattice_NtwinSystem(1:f-1,myStructure)) ! at which index starts my family + do i = 1,constitutive_dislotwin_Ntwin(f,myInstance) ! process each (active) slip system in family j = j + 1_pInt constitutive_dislotwin_postResults(c+j) = dot_product(Tstar_v,lattice_Stwin_v(:,index_myFamily+i,myStructure)) - enddo; enddo - endif + enddo; enddo + endif c = c + nt - case ('threshold_stress_twin') constitutive_dislotwin_postResults(c+1:c+nt) = state(g,ip,el)%p((6*ns+3*nt+1):(6*ns+4*nt)) c = c + nt - end select enddo return end function -END MODULE \ No newline at end of file +END MODULE