diff --git a/src/DAMASK_spectral.f90 b/src/DAMASK_spectral.f90 index dfa1746b2..bffc89790 100755 --- a/src/DAMASK_spectral.f90 +++ b/src/DAMASK_spectral.f90 @@ -75,7 +75,8 @@ program DAMASK_spectral FIELD_UNDEFINED_ID, & FIELD_MECH_ID, & FIELD_THERMAL_ID, & - FIELD_DAMAGE_ID + FIELD_DAMAGE_ID, & + utilities_calcPlasticity use spectral_mech_Basic use spectral_mech_AL use spectral_mech_Polarisation @@ -153,7 +154,21 @@ program DAMASK_spectral MPI_finalize, & MPI_allreduce, & PETScFinalize - +!-------------------------------------------------------------------------------------------------- +! variables related to stop criterion for yielding + real(pReal) :: plasticWorkOld, plasticWorkNew, & ! plastic work + eqTotalStrainOld, eqTotalStrainNew, & ! total equivalent strain + eqPlasticStrainOld, eqPlasticStrainNew, & ! total equivalent plastic strain + eqStressOld, eqStressNew , & ! equivalent stress + yieldStopValue + real(pReal), dimension(3,3) :: yieldStress,yieldStressOld,yieldStressNew, plasticStrainOld, plasticStrainNew, plasticStrainRate + integer(pInt) :: yieldResUnit = 0_pInt + character(len=13) :: stopFlag + logical :: yieldStop, yieldStopSatisfied + ! logical :: & + ! stop_totalStrain, & !< stop criterion + ! stop_plasticStrain, & + ! stop_plasticWork !-------------------------------------------------------------------------------------------------- ! init DAMASK (all modules) call CPFEM_initAll(el = 1_pInt, ip = 1_pInt) @@ -210,6 +225,8 @@ program DAMASK_spectral !-------------------------------------------------------------------------------------------------- ! reading the load case and assign values to the allocated data structure + yieldStop = .False. + yieldStopSatisfied = .False. rewind(FILEUNIT) do line = IO_read(FILEUNIT) @@ -284,6 +301,18 @@ program DAMASK_spectral temp_valueVector(j) = IO_floatValue(line,chunkPos,i+j) enddo loadCases(currentLoadCase)%rotation = math_plain9to33(temp_valueVector) + case('totalStrain', 'totalstrain') + yieldStop = .True. + stopFlag = 'totalStrain' + yieldStopValue = IO_floatValue(line,chunkPos,i+1_pInt) + case('plasticStrain', 'plasticstrain') + yieldStop = .True. + stopFlag = 'plasticStrain' + yieldStopValue = IO_floatValue(line,chunkPos,i+1_pInt) + case('plasticWork', 'plasticwork') + yieldStop = .True. + stopFlag = 'plasticWork' + yieldStopValue = IO_floatValue(line,chunkPos,i+1_pInt) end select enddo; enddo close(FILEUNIT) @@ -658,10 +687,60 @@ program DAMASK_spectral time = time + timeinc guess = .true. endif forwarding + + yieldCheck: if(yieldStop) then ! check if it yields or satisfies the certain stop condition + yieldStressOld = yieldStressNew + plasticStrainOld = plasticStrainNew + eqStressOld = eqStressNew + eqTotalStrainOld = eqTotalStrainNew + eqPlasticStrainOld = eqPlasticStrainNew + plasticWorkOld = plasticWorkNew + + call utilities_calcPlasticity(yieldStressNew, plasticStrainNew, eqStressNew, eqTotalStrainNew, eqPlasticStrainNew, plasticWorkNew) + + if(stopFlag .eq. 'totalStrain') then + if(eqTotalStrainNew > yieldStopValue) then + yieldStress = yieldStressOld * (eqTotalStrainNew - yieldStopValue)/(eqTotalStrainNew - eqTotalStrainOld) & ! linear interpolation of stress values + + yieldStressNew * (yieldStopValue - eqTotalStrainOld)/(eqTotalStrainNew - eqTotalStrainOld) + plasticStrainRate = (plasticStrainNew - plasticStrainOld)/(time - time0) ! calculate plastic strain rate + yieldStopSatisfied = .True. + endif + elseif(stopFlag .eq. 'plasticStrain') then + if(eqPlasticStrainNew > yieldStopValue) then + yieldStress = yieldStressOld * (eqPlasticStrainNew - yieldStopValue)/(eqPlasticStrainNew - eqPlasticStrainOld) & + + yieldStressNew * (yieldStopValue - eqPlasticStrainOld)/(eqPlasticStrainNew - eqPlasticStrainOld) + plasticStrainRate = (plasticStrainNew - plasticStrainOld)/(time - time0) + yieldStopSatisfied = .True. + endif + elseif(stopFlag .eq. 'plasticWork') then + if(plasticWorkNew > yieldStopValue) then + yieldStress = yieldStressOld * (plasticWorkNew - yieldStopValue)/(plasticWorkNew - plasticWorkOld) & + + yieldStressNew * (yieldStopValue - plasticWorkOld)/(plasticWorkNew - plasticWorkOld) + plasticStrainRate = (plasticStrainNew - plasticStrainOld)/(time - time0) + yieldStopSatisfied = .True. + endif + endif + endif yieldCheck + + if (yieldStopSatisfied) then ! when yield, write the yield stress and strain rate to file and quit the job + if (worldrank == 0) then + open(newunit=yieldResUnit,file=trim(getSolverWorkingDirectoryName())//trim(getSolverJobName())//& + '.yield',form='FORMATTED',status='REPLACE') + do i = 1_pInt,3_pInt + write(yieldResUnit,*) (yieldStress(i,j), j=1,3) + enddo + do i = 1_pInt,3_pInt + write(yieldResUnit,*) (plasticStrainRate(i,j), j=1,3) + enddo + close(yieldResUnit) + call quit(0_pInt) + endif + endif enddo incLooping enddo loadCaseLooping - + + !-------------------------------------------------------------------------------------------------- ! report summary of whole calculation write(6,'(/,a)') ' ###########################################################################' diff --git a/src/crystallite.f90 b/src/crystallite.f90 old mode 100644 new mode 100755 index c5bd4d979..d9f01e3d1 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -55,14 +55,14 @@ module crystallite crystallite_Li0, & !< intermediate velocitiy grad at start of FE inc crystallite_partionedLi0,& !< intermediate velocity grad at start of homog inc crystallite_Fe, & !< current "elastic" def grad (end of converged time step) - crystallite_P !< 1st Piola-Kirchhoff stress per grain + crystallite_P, & !< 1st Piola-Kirchhoff stress per grain + crystallite_subF !< def grad to be reached at end of crystallite inc real(pReal), dimension(:,:,:,:,:), allocatable, private :: & crystallite_subFe0,& !< "elastic" def grad at start of crystallite inc crystallite_invFp, & !< inverse of current plastic def grad (end of converged time step) crystallite_subFp0,& !< plastic def grad at start of crystallite inc crystallite_invFi, & !< inverse of current intermediate def grad (end of converged time step) crystallite_subFi0,& !< intermediate def grad at start of crystallite inc - crystallite_subF, & !< def grad to be reached at end of crystallite inc crystallite_subF0, & !< def grad at start of crystallite inc crystallite_subLp0,& !< plastic velocity grad at start of crystallite inc crystallite_subLi0,& !< intermediate velocity grad at start of crystallite inc diff --git a/src/math.f90 b/src/math.f90 old mode 100644 new mode 100755 index 3a0533ddf..937f00d7a --- a/src/math.f90 +++ b/src/math.f90 @@ -145,6 +145,7 @@ module math math_sampleGaussVar, & math_symmetricEulers, & math_eigenvectorBasisSym33, & + math_eigenvectorBasisSym33_log, & math_eigenvectorBasisSym, & math_eigenValuesVectorsSym33, & math_eigenValuesVectorsSym, & @@ -2064,6 +2065,70 @@ function math_eigenvectorBasisSym33(m) end function math_eigenvectorBasisSym33 +!-------------------------------------------------------------------------------------------------- +!> @brief logarithm eigenvector basis of symmetric 33 matrix m +!-------------------------------------------------------------------------------------------------- +function math_eigenvectorBasisSym33_log(m) + + implicit none + real(pReal), dimension(3,3) :: math_eigenvectorBasisSym33_log + real(pReal), dimension(3) :: invariants, values + real(pReal), dimension(3,3), intent(in) :: m + real(pReal) :: P, Q, rho, phi + real(pReal), parameter :: TOL=1.e-14_pReal + real(pReal), dimension(3,3,3) :: N, EB + + invariants = math_invariantsSym33(m) + EB = 0.0_pReal + + P = invariants(2)-invariants(1)**2.0_pReal/3.0_pReal + Q = -2.0_pReal/27.0_pReal*invariants(1)**3.0_pReal+product(invariants(1:2))/3.0_pReal-invariants(3) + + threeSimilarEigenvalues: if(all(abs([P,Q]) < TOL)) then + values = invariants(1)/3.0_pReal +! this is not really correct, but at least the basis is correct + EB(1,1,1)=1.0_pReal + EB(2,2,2)=1.0_pReal + EB(3,3,3)=1.0_pReal + else threeSimilarEigenvalues + rho=sqrt(-3.0_pReal*P**3.0_pReal)/9.0_pReal + phi=acos(math_limit(-Q/rho*0.5_pReal,-1.0_pReal,1.0_pReal)) + values = 2.0_pReal*rho**(1.0_pReal/3.0_pReal)* & + [cos(phi/3.0_pReal), & + cos((phi+2.0_pReal*PI)/3.0_pReal), & + cos((phi+4.0_pReal*PI)/3.0_pReal) & + ] + invariants(1)/3.0_pReal + N(1:3,1:3,1) = m-values(1)*math_I3 + N(1:3,1:3,2) = m-values(2)*math_I3 + N(1:3,1:3,3) = m-values(3)*math_I3 + twoSimilarEigenvalues: if(abs(values(1)-values(2)) < TOL) then + EB(1:3,1:3,3)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,2))/ & + ((values(3)-values(1))*(values(3)-values(2))) + EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,3) + elseif(abs(values(2)-values(3)) < TOL) then twoSimilarEigenvalues + EB(1:3,1:3,1)=math_mul33x33(N(1:3,1:3,2),N(1:3,1:3,3))/ & + ((values(1)-values(2))*(values(1)-values(3))) + EB(1:3,1:3,2)=math_I3-EB(1:3,1:3,1) + elseif(abs(values(3)-values(1)) < TOL) then twoSimilarEigenvalues + EB(1:3,1:3,2)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,3))/ & + ((values(2)-values(1))*(values(2)-values(3))) + EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,2) + else twoSimilarEigenvalues + EB(1:3,1:3,1)=math_mul33x33(N(1:3,1:3,2),N(1:3,1:3,3))/ & + ((values(1)-values(2))*(values(1)-values(3))) + EB(1:3,1:3,2)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,3))/ & + ((values(2)-values(1))*(values(2)-values(3))) + EB(1:3,1:3,3)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,2))/ & + ((values(3)-values(1))*(values(3)-values(2))) + endif twoSimilarEigenvalues + endif threeSimilarEigenvalues + + math_eigenvectorBasisSym33_log = log(sqrt(values(1))) * EB(1:3,1:3,1) & + + log(sqrt(values(2))) * EB(1:3,1:3,2) & + + log(sqrt(values(3))) * EB(1:3,1:3,3) + +end function math_eigenvectorBasisSym33_log + !-------------------------------------------------------------------------------------------------- !> @brief rotational part from polar decomposition of 33 tensor m !-------------------------------------------------------------------------------------------------- diff --git a/src/spectral_utilities.f90 b/src/spectral_utilities.f90 old mode 100644 new mode 100755 index 0773a9065..d0d720155 --- a/src/spectral_utilities.f90 +++ b/src/spectral_utilities.f90 @@ -145,7 +145,8 @@ module spectral_utilities FIELD_UNDEFINED_ID, & FIELD_MECH_ID, & FIELD_THERMAL_ID, & - FIELD_DAMAGE_ID + FIELD_DAMAGE_ID, & + utilities_calcPlasticity private :: & utilities_getFreqDerivative @@ -1040,6 +1041,145 @@ subroutine utilities_constitutiveResponse(F_lastInc,F,timeinc, & end subroutine utilities_constitutiveResponse +!-------------------------------------------------------------------------------------------------- +!> @brief calculates yield stress, plastic strain, total strain and their equivalent values +!-------------------------------------------------------------------------------------------------- +subroutine utilities_calcPlasticity(yieldStress, plasticStrain, eqStress, eqTotalStrain, eqPlasticStrain, plasticWork) + use crystallite, only: & + crystallite_Fp, & + crystallite_P, & + crystallite_subF + use material, only: & + homogenization_maxNgrains + use mesh, only: & + mesh_maxNips,& + mesh_NcpElems + use math, only: & + math_det33, & + math_inv33, & + math_mul33x33, & + math_trace33, & + math_transpose33, & + math_identity2nd, & + math_crossproduct, & + math_eigenvectorBasisSym, & + math_eigenvectorBasisSym33, & + math_eigenvectorBasisSym33_log, & + math_eigenValuesVectorsSym33 + + implicit none + + real(pReal), intent(out) :: eqStress, eqTotalStrain, eqPlasticStrain, plasticWork + real(pReal), dimension(3,3),intent(out) :: yieldStress, plasticStrain + real(pReal), dimension(3,3) :: cauchy, cauchy_av, P, Vp_av, V_total_av !< average + real(pReal), dimension(3) :: Values, S + real(pReal), dimension(3,3) :: Vectors, diag + real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & + Bp, Vp, F, F_temp, Fp, U, VT, R, V, V_total + real(pReal), dimension(15) :: WORK !< previous deformation gradient + integer(pInt) :: INFO, i, j, k, l, ierr + real(pReal) :: stress, stress_av, strain_total, strain_total_av, strain_plastic, strain_plastic_av, wgtm + + external :: dgesvd + + wgtm = 1.0/real(mesh_NcpElems*mesh_maxNips*homogenization_maxNgrains,pReal) + Vp_av = 0.0_pReal + V_total_av = 0.0_pReal + strain_plastic_av = 0.0_pReal + strain_total_av = 0.0_pReal + cauchy_av = 0.0_pReal + diag = 0.0_pReal + do k = 1_pInt, mesh_NcpElems; do j = 1_pInt, mesh_maxNips; do i = 1_pInt,homogenization_maxNgrains + F(1:3,1:3,i,j,k) = crystallite_subF(1:3,1:3,i,j,k) + Fp(1:3,1:3,i,j,k) = crystallite_Fp(1:3,1:3,i,j,k) + + P = crystallite_P(1:3,1:3,i,j,k) + cauchy = 1.0/math_det33(F(1:3,1:3,i,j,k))*math_mul33x33(P,transpose(F(1:3,1:3,i,j,k))) + cauchy_av = cauchy_av + cauchy + stress = Mises(cauchy, 'stress') + stress_av = stress_av + stress + + Bp(1:3,1:3,i,j,k) = math_mul33x33(Fp(1:3,1:3,i,j,k),math_transpose33(Fp(1:3,1:3,i,j,k))) ! plastic part of left Cauchy–Green deformation tensor + Vp(1:3,1:3,i,j,k) = math_eigenvectorBasisSym33_log(Bp(1:3,1:3,i,j,k)) + strain_plastic = Mises(Vp(1:3,1:3,i,j,k), 'strain') + strain_plastic_av = strain_plastic_av + strain_plastic + Vp_av = Vp_av + Vp(1:3,1:3,i,j,k) + + F_temp(1:3,1:3,i,j,k) = F(1:3,1:3,i,j,k) + call dgesvd ('A', 'A', 3, 3, F_temp(1:3,1:3,i,j,k), 3, & + S, U(1:3,1:3,i,j,k), 3, VT(1:3,1:3,i,j,k), 3, WORK, 15, INFO) ! singular value decomposition + R(1:3,1:3,i,j,k) = math_mul33x33(U(1:3,1:3,i,j,k), VT(1:3,1:3,i,j,k)) ! rotation of polar decomposition + V(1:3,1:3,i,j,k) = math_mul33x33(F(1:3,1:3,i,j,k),math_inv33(R(1:3,1:3,i,j,k))) + + call math_eigenValuesVectorsSym33(V(1:3,1:3,i,j,k),Values,Vectors) + do l = 1_pInt, 3_pInt + if (Values(l) < 0.0_pReal) then + Values(l) = -Values(l) + Vectors(1:3, l) = -Vectors(1:3, l) + endif + Values(l) = log(Values(l)) + diag(l,l) = Values(l) + enddo + if (dot_product(Vectors(1:3,1),Vectors(1:3,2)) /= 0) then + Vectors(1:3,2) = math_crossproduct(Vectors(1:3,3), Vectors(1:3,1)) + Vectors(1:3,2) = Vectors(1:3,2)/sqrt(dot_product(Vectors(1:3,2),Vectors(1:3,2))) + endif + if (dot_product(Vectors(1:3,2),Vectors(1:3,3)) /= 0) then + Vectors(1:3,3) = math_crossproduct(Vectors(1:3,1), Vectors(1:3,2)) + Vectors(1:3,3) = Vectors(1:3,3)/sqrt(dot_product(Vectors(1:3,3),Vectors(1:3,3))) + endif + if (dot_product(Vectors(1:3,3),Vectors(1:3,1)) /= 0) then + Vectors(1:3,1) = math_crossproduct(Vectors(1:3,2), Vectors(1:3,3)) + Vectors(1:3,1) = Vectors(1:3,1)/sqrt(dot_product(Vectors(1:3,1),Vectors(1:3,1))) + endif + + V_total(1:3,1:3,i,j,k) = REAL(math_mul33x33(Vectors, math_mul33x33(diag, transpose(Vectors)))) + strain_total = Mises(V_total(1:3,1:3,i,j,k), 'strain') + strain_total_av = strain_total_av + strain_total + + enddo; enddo; enddo + + + yieldStress = cauchy_av * wgtm + call MPI_Allreduce(MPI_IN_PLACE,yieldStress,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + plasticStrain = Vp_av * wgtm + call MPI_Allreduce(MPI_IN_PLACE,plasticStrain,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + + plasticWork = plasticWork + 0.5*(eqStress + stress_av * wgtm) * (strain_total_av * wgtm - eqTotalStrain) + call MPI_Allreduce(MPI_IN_PLACE,plasticWork,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + eqStress = stress_av * wgtm + call MPI_Allreduce(MPI_IN_PLACE,eqStress,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + eqTotalStrain = strain_total_av * wgtm + call MPI_Allreduce(MPI_IN_PLACE,eqTotalStrain,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + eqPlasticStrain = strain_plastic_av * wgtm + call MPI_Allreduce(MPI_IN_PLACE,eqPlasticStrain,1,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) + +end subroutine utilities_calcPlasticity + +!-------------------------------------------------------------------------------------------------- +!> @brief vonMises equivalent values for symmetric part of requested strains and/or stresses. +!-------------------------------------------------------------------------------------------------- +function Mises(m, mType) + use math, only: & + math_trace33, & + math_identity2nd + + implicit none + real(pReal), intent(in), dimension(3,3) :: m + character(len=*), intent(in) :: mType + real(pReal) :: trace, Mises + real(pReal), dimension(3,3) :: eye, dev, symdev + + eye = math_identity2nd(3) + trace = math_trace33(m) + dev = m - trace/3.0*eye + symdev = 0.5*(dev+transpose(dev)) + if (mType .eq. 'stress') then + Mises = sqrt(3.0/2.0*sum(symdev*transpose(symdev))) + else + Mises = sqrt(2.0/3.0*sum(symdev*transpose(symdev))) + endif +end function Mises !-------------------------------------------------------------------------------------------------- !> @brief calculates forward rate, either guessing or just add delta/timeinc