223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

implemented yield stop criteria

This commit is contained in:
Fengbo Han 2017-07-27 16:21:02 +02:00
parent a4ceebc00b
commit 5cedba0721
4 changed files with 290 additions and 6 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

85
src/DAMASK_spectral.f90
View File

@ -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)') ' ###########################################################################'

4
src/crystallite.f90 Normal file → Executable file
View File

@ -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

65
src/math.f90 Normal file → Executable file
View File

@ -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
!--------------------------------------------------------------------------------------------------

142
src/spectral_utilities.f90 Normal file → Executable file
View File

@ -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