From a3e1b712404ff658d3964a74b878203bf1354aec Mon Sep 17 00:00:00 2001 From: Jean-Michel Scherer Date: Wed, 22 Jun 2022 11:22:19 -0700 Subject: [PATCH 01/16] Implement the rotation of the stiffness tensor from the local frame (yaml input) to the global frame. The rotation matrix follows the Bunge (passive) convention (ZXZ) and can be defined from 3 angles (phi1, Phi, phi2) or three vectors (V1,V2,V3). --- MEF90/m_MEF90_Materials.F90 | 141 +++++++++++++++++++++++++++++++++--- 1 file changed, 130 insertions(+), 11 deletions(-) diff --git a/MEF90/m_MEF90_Materials.F90 b/MEF90/m_MEF90_Materials.F90 index a5789a3085..c1b22b5cb2 100644 --- a/MEF90/m_MEF90_Materials.F90 +++ b/MEF90/m_MEF90_Materials.F90 @@ -7,13 +7,14 @@ Module m_MEF90_Materials_Types Type MEF90HookesLaw2D Type(Tens4OS2D) :: fullTensor + Type(Tens4OS3D) :: fullTensorLocal,fullTensor3D PetscReal :: lambda,mu,YoungsModulus,PoissonRatio,BulkModulus PetscEnum :: type PetscBool :: isPlaneStress End Type MEF90HookesLaw2D Type MEF90HookesLaw3D - Type(Tens4OS3D) :: fullTensor + Type(Tens4OS3D) :: fullTensor,fullTensorLocal PetscReal :: lambda,mu,YoungsModulus,PoissonRatio,BulkModulus PetscEnum :: type #if (PETSC_SIZEOF_INT == 4) @@ -25,6 +26,17 @@ Module m_MEF90_Materials_Types #endif End Type MEF90HookesLaw3D + Type MEF90RotationMatrix3D + ! The Bunge (passive) convention is used. The rotation matrix transforms a vector from the global frame to the local frame. + ! - rotation of a vector: [V_local]_i = [R]_ij . [V_global]_j + ! - rotation of a 2nd order tensor: [M_local]_ij = [R]_ik . [M_global]_kl . [R^T]_lj + ! - rotation of a fourth order tensor: [C_local]_ijkl = [R]_ip . [R]_jq . [R]_kr . [R]_ls . [C_global]_pqrs + Type(MAT3D) :: fullTensor + PetscReal :: phi1,Phi,phi2 + Type(Vect3D) :: V1,V2,V3 + PetscBool :: fromEuler + End Type MEF90RotationMatrix3D + Type MEF90MatProp2D_Type PetscReal :: Density ! rho PetscReal :: FractureToughness ! Gc @@ -62,6 +74,7 @@ Module m_MEF90_Materials_Types PetscReal :: drivingForceGamma ! gamma parameter in Drucker-Prager driving force PetscBool :: isLinearIsotropicHardening PetscBool :: isNoPlCoupling + Type(MEF90RotationMatrix3D) :: RotationMatrix ! rotation matrix from the global frame to the material frame: X_local = R . X_global Character(len=MEF90_MXSTRLEN) :: Name End Type MEF90MatProp2D_Type @@ -102,6 +115,7 @@ Module m_MEF90_Materials_Types PetscReal :: drivingForceGamma ! gamma parameter in Drucker-Prager driving force PetscBool :: isLinearIsotropicHardening PetscBool :: isNoPlCoupling + Type(MEF90RotationMatrix3D) :: RotationMatrix ! rotation matrix from the global frame to the material frame: X_local = R . X_global Character(len=MEF90_MXSTRLEN) :: Name End Type MEF90MatProp3D_Type @@ -123,6 +137,18 @@ Module m_MEF90_Materials_Types Tens4OS2D(1.09890_Kr,0.32967_Kr,0.00000_Kr, & ! HookesLaw XXXX,XXYY,XXXY 1.09890_Kr,0.00000_Kr, & ! YYYY,YYXY 0.38462_Kr), & ! XYXY + Tens4OS3D(1.34615_Kr,0.57692_Kr,0.57692_Kr,0.00000_Kr,0.00000_Kr,0.00000_Kr, & ! HookesLaw Local XXXX,XXYY,XXZZ,XXYZ,XXXZ,XXXY + 1.34615_Kr,0.57692_Kr,0.00000_Kr,0.00000_Kr,0.00000_Kr, & ! YYYY,YYZZ,YYYZ,YYXZ,YYXY + 1.34615_Kr,0.00000_Kr,0.00000_Kr,0.00000_Kr, & ! ZZZZ ZZYZ,ZZXZ,ZZXY + 0.38462_Kr,0.00000_Kr,0.00000_Kr, & ! YXYX,YZXZ,YZXY + 0.38462_Kr,0.00000_Kr, & ! XZXZ,XZXY + 0.38462_Kr), & ! + Tens4OS3D(1.34615_Kr,0.57692_Kr,0.57692_Kr,0.00000_Kr,0.00000_Kr,0.00000_Kr, & ! HookesLaw 3D XXXX,XXYY,XXZZ,XXYZ,XXXZ,XXXY + 1.34615_Kr,0.57692_Kr,0.00000_Kr,0.00000_Kr,0.00000_Kr, & ! YYYY,YYZZ,YYYZ,YYXZ,YYXY + 1.34615_Kr,0.00000_Kr,0.00000_Kr,0.00000_Kr, & ! ZZZZ ZZYZ,ZZXZ,ZZXY + 0.38462_Kr,0.00000_Kr,0.00000_Kr, & ! YXYX,YZXZ,YZXY + 0.38462_Kr,0.00000_Kr, & ! XZXZ,XZXY + 0.38462_Kr), & ! 0.0_Kr,0.0_Kr,1.0_Kr,.3_Kr,0.0_Kr, & ! lambda, mu, E, nu, kappa (lambda, mu, kappa will be recomputed) MEF90HookesLawTypeIsotropic, & ! type .FALSE.), & ! isPlaneStress @@ -155,6 +181,9 @@ Module m_MEF90_Materials_Types 0.0_Kr, & ! drivingForceGamma .FALSE., & ! isLinearIsotropicHardening .FALSE., & ! isNoPlCoupling + MEF90RotationMatrix3D(MEF90Mat3DIdentity,0.0_Kr,0.0_Kr,0.0_Kr, & ! RotationMatrix, phi1, Phi, phi2, + Vect3D(1.0_Kr,0.0_Kr,0.0_Kr),Vect3D(0.0_Kr,1.0_Kr,0.0_Kr), & ! V1, V2, + Vect3D(0.0_Kr,0.0_Kr,1.0_Kr),.False.), & ! V3, fromEuler "MEF90Mathium2D") Type(MEF90MatProp3D_Type),Parameter :: MEF90Mathium3D = MEF90MatProp3D_Type( & @@ -165,12 +194,18 @@ Module m_MEF90_Materials_Types MEF90MatS3DIdentity, & ! ThermalConductivity MEF90MatS3DIdentity, & ! LinearThermalExpansion MEF90HookesLaw3D( & - Tens4OS3D(1.34615_Kr,0.57692_Kr,0.57692_Kr,0.00000_Kr,0.00000_Kr,0.00000_Kr, & ! XXXX,XXYY,XXZZ,XXYZ,XXXZ,XXXY - 1.34615_Kr,0.57692_Kr,0.00000_Kr,0.00000_Kr,0.00000_Kr, & ! YYYY,YYZZ,YYYZ,YYXZ,YYXY - 1.34615_Kr,0.00000_Kr,0.00000_Kr,0.00000_Kr, & ! ZZZZ ZZYZ,ZZXZ,ZZXY - 0.38462_Kr,0.00000_Kr,0.00000_Kr, & ! YXYX,YZXZ,YZXY - 0.38462_Kr,0.00000_Kr, & ! XZXZ,XZXY - 0.38462_Kr), & ! XYXY + Tens4OS3D(1.34615_Kr,0.57692_Kr,0.57692_Kr,0.00000_Kr,0.00000_Kr,0.00000_Kr, & ! HookesLaw XXXX,XXYY,XXZZ,XXYZ,XXXZ,XXXY + 1.34615_Kr,0.57692_Kr,0.00000_Kr,0.00000_Kr,0.00000_Kr, & ! YYYY,YYZZ,YYYZ,YYXZ,YYXY + 1.34615_Kr,0.00000_Kr,0.00000_Kr,0.00000_Kr, & ! ZZZZ ZZYZ,ZZXZ,ZZXY + 0.38462_Kr,0.00000_Kr,0.00000_Kr, & ! YXYX,YZXZ,YZXY + 0.38462_Kr,0.00000_Kr, & ! XZXZ,XZXY + 0.38462_Kr), & ! XYXY + Tens4OS3D(1.34615_Kr,0.57692_Kr,0.57692_Kr,0.00000_Kr,0.00000_Kr,0.00000_Kr, & ! HookesLaw Local XXXX,XXYY,XXZZ,XXYZ,XXXZ,XXXY + 1.34615_Kr,0.57692_Kr,0.00000_Kr,0.00000_Kr,0.00000_Kr, & ! YYYY,YYZZ,YYYZ,YYXZ,YYXY + 1.34615_Kr,0.00000_Kr,0.00000_Kr,0.00000_Kr, & ! ZZZZ ZZYZ,ZZXZ,ZZXY + 0.38462_Kr,0.00000_Kr,0.00000_Kr, & ! YXYX,YZXZ,YZXY + 0.38462_Kr,0.00000_Kr, & ! XZXZ,XZXY + 0.38462_Kr), & ! 0.0_Kr,0.0_Kr,1.0_Kr,.3_Kr,0.0_Kr, & ! lambda, mu, E, nu, kappa (lambda, mu, kappa will be recomputed) MEF90HookesLawTypeIsotropic, & ! type #if (PETSC_SIZEOF_INT == 4) @@ -206,6 +241,9 @@ Module m_MEF90_Materials_Types 2, & ! drivingForcep .FALSE., & ! isLinearIsotropicHardening .FALSE., & ! isNoPlCoupling + MEF90RotationMatrix3D(MEF90Mat3DIdentity,0.0_Kr,0.0_Kr,0.0_Kr, & ! RotationMatrix, phi1, Phi, phi2, + Vect3D(1.0_Kr,0.0_Kr,0.0_Kr),Vect3D(0.0_Kr,1.0_Kr,0.0_Kr), & ! V1, V2, + Vect3D(0.0_Kr,0.0_Kr,1.0_Kr),.False.), & ! V3, fromEuler "MEF90Mathium3D") End Module m_MEF90_Materials_Types @@ -233,6 +271,8 @@ Subroutine PetscBagGetDataMEF90MatProp2D(bag,data,ierr) PetscBag :: bag type(MEF90MatProp2D_Type),pointer :: data PetscErrorCode :: ierr + PetscReal :: normV1,normV2,normV3 + Type(Tens4OS3D) :: HookesLaw3D Call PetscBagGetData(bag,data,ierr) Select case(data%HookesLaw%type) @@ -249,6 +289,38 @@ Subroutine PetscBagGetDataMEF90MatProp2D(bag,data,ierr) data%HookesLaw%BulkModulus = data%HookesLaw%lambda + data%HookesLaw%mu End If End Select + If (data%RotationMatrix%fromEuler) Then + data%RotationMatrix%fullTensor%XX = cos(data%RotationMatrix%phi1)*cos(data%RotationMatrix%phi2)-sin(data%RotationMatrix%phi1)*sin(data%RotationMatrix%phi2)*cos(data%RotationMatrix%Phi) + data%RotationMatrix%fullTensor%XY = sin(data%RotationMatrix%phi1)*cos(data%RotationMatrix%phi2)+cos(data%RotationMatrix%phi1)*sin(data%RotationMatrix%phi2)*cos(data%RotationMatrix%Phi) + data%RotationMatrix%fullTensor%YX = -cos(data%RotationMatrix%phi1)*sin(data%RotationMatrix%phi2)-sin(data%RotationMatrix%phi1)*cos(data%RotationMatrix%phi2)*cos(data%RotationMatrix%Phi) + data%RotationMatrix%fullTensor%YY = -sin(data%RotationMatrix%phi1)*sin(data%RotationMatrix%phi2)+cos(data%RotationMatrix%phi1)*cos(data%RotationMatrix%phi2)*cos(data%RotationMatrix%Phi) + data%RotationMatrix%fullTensor%XZ = sin(data%RotationMatrix%phi2)*sin(data%RotationMatrix%Phi) + data%RotationMatrix%fullTensor%YZ = cos(data%RotationMatrix%phi2)*sin(data%RotationMatrix%Phi) + data%RotationMatrix%fullTensor%ZX = sin(data%RotationMatrix%phi1)*sin(data%RotationMatrix%Phi) + data%RotationMatrix%fullTensor%ZY = -cos(data%RotationMatrix%phi1)*sin(data%RotationMatrix%Phi) + data%RotationMatrix%fullTensor%ZZ = cos(data%RotationMatrix%Phi) + Else + normV1 = SQRT(data%RotationMatrix%V1%X*data%RotationMatrix%V1%X + data%RotationMatrix%V1%Y*data%RotationMatrix%V1%Y + data%RotationMatrix%V1%Z*data%RotationMatrix%V1%Z) + normV2 = SQRT(data%RotationMatrix%V2%X*data%RotationMatrix%V2%X + data%RotationMatrix%V2%Y*data%RotationMatrix%V2%Y + data%RotationMatrix%V2%Z*data%RotationMatrix%V2%Z) + normV3 = SQRT(data%RotationMatrix%V3%X*data%RotationMatrix%V3%X + data%RotationMatrix%V3%Y*data%RotationMatrix%V3%Y + data%RotationMatrix%V3%Z*data%RotationMatrix%V3%Z) + data%RotationMatrix%fullTensor%XX = data%RotationMatrix%V1%X / normV1 + data%RotationMatrix%fullTensor%YX = data%RotationMatrix%V1%Y / normV1 + data%RotationMatrix%fullTensor%ZX = data%RotationMatrix%V1%Z / normV1 + data%RotationMatrix%fullTensor%XY = data%RotationMatrix%V2%X / normV2 + data%RotationMatrix%fullTensor%YY = data%RotationMatrix%V2%Y / normV2 + data%RotationMatrix%fullTensor%ZY = data%RotationMatrix%V2%Z / normV2 + data%RotationMatrix%fullTensor%XZ = data%RotationMatrix%V3%X / normV3 + data%RotationMatrix%fullTensor%YZ = data%RotationMatrix%V3%Y / normV3 + data%RotationMatrix%fullTensor%ZZ = data%RotationMatrix%V3%Z / normV3 + End If + HookesLaw3D = Tens4OSTransform(data%HookesLaw%fullTensorLocal,transpose(data%RotationMatrix%fullTensor)) + data%HookesLaw%fullTensor3D = HookesLaw3D + data%HookesLaw%fullTensor%XXXX = HookesLaw3D%XXXX + data%HookesLaw%fullTensor%XXYY = HookesLaw3D%XXYY + data%HookesLaw%fullTensor%XXXY = HookesLaw3D%XXXY + data%HookesLaw%fullTensor%YYYY = HookesLaw3D%YYYY + data%HookesLaw%fullTensor%YYXY = HookesLaw3D%YYXY + data%HookesLaw%fullTensor%XYXY = HookesLaw3D%XYXY End Subroutine PetscBagGetDataMEF90MatProp2D End Module m_MEF90_Materials_Interface2D @@ -276,6 +348,7 @@ Subroutine PetscBagGetDataMEF90MatProp3D(bag,data,ierr) PetscBag :: bag type(MEF90MatProp3D_Type),pointer :: data PetscErrorCode :: ierr + PetscReal :: normV1,normV2,normV3 Call PetscBagGetData(bag,data,ierr) Select case(data%HookesLaw%type) @@ -286,6 +359,31 @@ Subroutine PetscBagGetDataMEF90MatProp3D(bag,data,ierr) data%HookesLaw%mu = data%HookesLaw%YoungsModulus / (1.0_Kr + data%HookesLaw%PoissonRatio) * .5_Kr data%HookesLaw%BulkModulus = data%HookesLaw%lambda + 2.0_Kr * data%HookesLaw%mu / 3.0_Kr End Select + If (data%RotationMatrix%fromEuler) Then + data%RotationMatrix%fullTensor%XX = cos(data%RotationMatrix%phi1)*cos(data%RotationMatrix%phi2)-sin(data%RotationMatrix%phi1)*sin(data%RotationMatrix%phi2)*cos(data%RotationMatrix%Phi) + data%RotationMatrix%fullTensor%XY = sin(data%RotationMatrix%phi1)*cos(data%RotationMatrix%phi2)+cos(data%RotationMatrix%phi1)*sin(data%RotationMatrix%phi2)*cos(data%RotationMatrix%Phi) + data%RotationMatrix%fullTensor%YX = -cos(data%RotationMatrix%phi1)*sin(data%RotationMatrix%phi2)-sin(data%RotationMatrix%phi1)*cos(data%RotationMatrix%phi2)*cos(data%RotationMatrix%Phi) + data%RotationMatrix%fullTensor%YY = -sin(data%RotationMatrix%phi1)*sin(data%RotationMatrix%phi2)+cos(data%RotationMatrix%phi1)*cos(data%RotationMatrix%phi2)*cos(data%RotationMatrix%Phi) + data%RotationMatrix%fullTensor%XZ = sin(data%RotationMatrix%phi2)*sin(data%RotationMatrix%Phi) + data%RotationMatrix%fullTensor%YZ = cos(data%RotationMatrix%phi2)*sin(data%RotationMatrix%Phi) + data%RotationMatrix%fullTensor%ZX = sin(data%RotationMatrix%phi1)*sin(data%RotationMatrix%Phi) + data%RotationMatrix%fullTensor%ZY = -cos(data%RotationMatrix%phi1)*sin(data%RotationMatrix%Phi) + data%RotationMatrix%fullTensor%ZZ = cos(data%RotationMatrix%Phi) + Else + normV1 = SQRT(data%RotationMatrix%V1%X*data%RotationMatrix%V1%X + data%RotationMatrix%V1%Y*data%RotationMatrix%V1%Y + data%RotationMatrix%V1%Z*data%RotationMatrix%V1%Z) + normV2 = SQRT(data%RotationMatrix%V2%X*data%RotationMatrix%V2%X + data%RotationMatrix%V2%Y*data%RotationMatrix%V2%Y + data%RotationMatrix%V2%Z*data%RotationMatrix%V2%Z) + normV3 = SQRT(data%RotationMatrix%V3%X*data%RotationMatrix%V3%X + data%RotationMatrix%V3%Y*data%RotationMatrix%V3%Y + data%RotationMatrix%V3%Z*data%RotationMatrix%V3%Z) + data%RotationMatrix%fullTensor%XX = data%RotationMatrix%V1%X / normV1 + data%RotationMatrix%fullTensor%YX = data%RotationMatrix%V1%Y / normV1 + data%RotationMatrix%fullTensor%ZX = data%RotationMatrix%V1%Z / normV1 + data%RotationMatrix%fullTensor%XY = data%RotationMatrix%V2%X / normV2 + data%RotationMatrix%fullTensor%YY = data%RotationMatrix%V2%Y / normV2 + data%RotationMatrix%fullTensor%ZY = data%RotationMatrix%V2%Z / normV2 + data%RotationMatrix%fullTensor%XZ = data%RotationMatrix%V3%X / normV3 + data%RotationMatrix%fullTensor%YZ = data%RotationMatrix%V3%Y / normV3 + data%RotationMatrix%fullTensor%ZZ = data%RotationMatrix%V3%Z / normV3 + End If + data%HookesLaw%fullTensor = Tens4OSTransform(data%HookesLaw%fullTensorLocal,transpose(data%RotationMatrix%fullTensor)) End Subroutine PetscBagGetDataMEF90MatProp3D End Module m_MEF90_Materials_Interface3D @@ -396,13 +494,14 @@ Subroutine PetscBagRegisterMEF90MatProp2D(bag,name,prefix,default,ierr) Call PetscBagRegisterEnum(bag,matprop%HookesLaw%type,MEF90HookesLawTypeList,default%HookesLaw%type,'hookeslaw_type','Type of Hooke''s law',ierr);CHKERRQ(ierr) Select case(matprop%HookesLaw%type) Case (MEF90HookesLawTypeFull) - matprop%HookesLaw%fullTensor = default%HookesLaw%fullTensor - Call PetscBagRegisterRealArray(bag,matprop%HookesLaw%fullTensor,6,'HookesLaw_tensor','[N.m^(-2)] (A) Hooke''s law',ierr) + matprop%HookesLaw%fullTensorLocal = default%HookesLaw%fullTensorLocal + Call PetscBagRegisterRealArray(bag,matprop%HookesLaw%fullTensorLocal,21,'HookesLaw_tensor','[N.m^(-2)] (A) Hooke''s law in the local frame',ierr) Case(MEF90HookesLawTypeIsotropic) Call PetscBagRegisterReal(bag,matprop%HookesLaw%YoungsModulus,default%HookesLaw%YoungsModulus,'hookeslaw_YoungsModulus','[N.m^(-2)] (E) Young''s Modulus',ierr) Call PetscBagRegisterReal(bag,matprop%HookesLaw%PoissonRatio,default%HookesLaw%PoissonRatio,'hookeslaw_PoissonRatio','[] (nu) Poisson Modulus',ierr) Call PetscBagRegisterBool(bag,matprop%HookesLaw%isPlaneStress,default%HookesLaw%isPlaneStress,'hookeslaw_planeStress','Use plane stress elasticity',ierr);CHKERRQ(ierr) matprop%HookesLaw%fulltensor = -1.D+30 + matprop%HookesLaw%fulltensorLocal = -1.D+30 End Select Call PetscBagRegisterReal(bag,matprop%internalLength,default%internalLength,'internalLength','[m] (l) Internal Length',ierr) @@ -438,6 +537,16 @@ Subroutine PetscBagRegisterMEF90MatProp2D(bag,name,prefix,default,ierr) Call PetscBagRegisterBool(bag,matprop%isLinearIsotropicHardening,default%isLinearIsotropicHardening,'isLinearIsotropicHardening','[bool] Plasticity with Linear Isotropic Hardening',ierr);CHKERRQ(ierr) Call PetscBagRegisterBool(bag,matprop%isNoPlCoupling,default%isNoPlCoupling,'isNoPlCoupling','[bool] Coupling between damage and plastic dissipation',ierr);CHKERRQ(ierr) + Call PetscBagRegisterReal(bag,matprop%RotationMatrix%phi1,default%RotationMatrix%phi1,'RotationMatrix_phi1','[radians] (phi1) First Bunge-Euler angle',ierr) + Call PetscBagRegisterReal(bag,matprop%RotationMatrix%Phi,default%RotationMatrix%Phi,'RotationMatrix_Phi','[radians] (Phi) Second Bunge-Euler angle',ierr) + Call PetscBagRegisterReal(bag,matprop%RotationMatrix%phi2,default%RotationMatrix%phi2,'RotationMatrix_phi2','[radians] (phi2) Third Bunge-Euler angle',ierr) + matprop%RotationMatrix%V1 = default%RotationMatrix%V1 + Call PetscBagRegisterRealArray(bag,matprop%RotationMatrix%V1,3,'RotationMatrix_V1','[] (V1) First column of the rotation matrix',ierr) + matprop%RotationMatrix%V2 = default%RotationMatrix%V2 + Call PetscBagRegisterRealArray(bag,matprop%RotationMatrix%V2,3,'RotationMatrix_V2','[] (V2) Second column of the rotation matrix',ierr) + matprop%RotationMatrix%V3 = default%RotationMatrix%V3 + Call PetscBagRegisterRealArray(bag,matprop%RotationMatrix%V3,3,'RotationMatrix_V3','[] (V3) Third column of the rotation matrix',ierr) + Call PetscBagRegisterBool(bag,matprop%RotationMatrix%fromEuler,default%RotationMatrix%fromEuler,'RotationMatrix_fromEuler','Define rotation matrix from Bunge-Euler angles',ierr);CHKERRQ(ierr) !Call PetscBagSetFromOptions(bag,ierr) End Subroutine PetscBagRegisterMEF90MatProp2D @@ -474,8 +583,8 @@ Subroutine PetscBagRegisterMEF90MatProp3D(bag,name,prefix,default,ierr) Call PetscBagRegisterEnum(bag,matprop%HookesLaw%type,MEF90HookesLawTypeList,default%HookesLaw%type,'hookeslaw_type','Type of Hooke''s law',ierr);CHKERRQ(ierr) Select case(matprop%HookesLaw%type) Case (MEF90HookesLawTypeFull) - matprop%HookesLaw%fullTensor = default%HookesLaw%fullTensor - Call PetscBagRegisterRealArray(bag,matprop%HookesLaw%fullTensor,21,'HookesLaw','[N.m^(-2)] (A) Hooke''s law',ierr) + matprop%HookesLaw%fullTensorLocal = default%HookesLaw%fullTensorLocal + Call PetscBagRegisterRealArray(bag,matprop%HookesLaw%fullTensorLocal,21,'HookesLaw_tensor','[N.m^(-2)] (A) Hooke''s law in the local frame',ierr) Case(MEF90HookesLawTypeIsotropic) Call PetscBagRegisterReal(bag,matprop%HookesLaw%YoungsModulus,default%HookesLaw%YoungsModulus,'hookeslaw_YoungsModulus','[N.m^(-2)] (E) Young''s Modulus',ierr) Call PetscBagRegisterReal(bag,matprop%HookesLaw%PoissonRatio,default%HookesLaw%PoissonRatio,'hookeslaw_PoissonRatio','[] (nu) Poisson Modulus',ierr) @@ -514,6 +623,16 @@ Subroutine PetscBagRegisterMEF90MatProp3D(bag,name,prefix,default,ierr) Call PetscBagRegisterBool(bag,matprop%isLinearIsotropicHardening,default%isLinearIsotropicHardening,'isLinearIsotropicHardening','[bool] Plasticity with Linear Isotropic Hardening',ierr);CHKERRQ(ierr) Call PetscBagRegisterBool(bag,matprop%isNoPlCoupling,default%isNoPlCoupling,'isNoPlCoupling','[bool] Coupling between damage and plastic dissipation',ierr);CHKERRQ(ierr) + Call PetscBagRegisterReal(bag,matprop%RotationMatrix%phi1,default%RotationMatrix%phi1,'RotationMatrix_phi1','[radians] (phi1) First Bunge-Euler angle',ierr) + Call PetscBagRegisterReal(bag,matprop%RotationMatrix%Phi,default%RotationMatrix%Phi,'RotationMatrix_Phi','[radians] (Phi) Second Bunge-Euler angle',ierr) + Call PetscBagRegisterReal(bag,matprop%RotationMatrix%phi2,default%RotationMatrix%phi2,'RotationMatrix_phi2','[radians] (phi2) Third Bunge-Euler angle',ierr) + matprop%RotationMatrix%V1 = default%RotationMatrix%V1 + Call PetscBagRegisterRealArray(bag,matprop%RotationMatrix%V1,3,'RotationMatrix_V1','[] (V1) First column of the rotation matrix',ierr) + matprop%RotationMatrix%V2 = default%RotationMatrix%V2 + Call PetscBagRegisterRealArray(bag,matprop%RotationMatrix%V2,3,'RotationMatrix_V2','[] (V2) Second column of the rotation matrix',ierr) + matprop%RotationMatrix%V3 = default%RotationMatrix%V3 + Call PetscBagRegisterRealArray(bag,matprop%RotationMatrix%V3,3,'RotationMatrix_V3','[] (V3) Third column of the rotation matrix',ierr) + Call PetscBagRegisterBool(bag,matprop%RotationMatrix%fromEuler,default%RotationMatrix%fromEuler,'RotationMatrix_fromEuler','Define rotation matrix from Bunge-Euler angles',ierr);CHKERRQ(ierr) !Call PetscBagSetFromOptions(bag,ierr) End Subroutine PetscBagRegisterMEF90MatProp3D From a65890758fa4246480bb33fabcfbde062510bd6b Mon Sep 17 00:00:00 2001 From: Jean-Michel Scherer Date: Wed, 22 Jun 2022 11:23:08 -0700 Subject: [PATCH 02/16] fix MPI types --- MEF90/m_MEF90_EXO.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/MEF90/m_MEF90_EXO.F90 b/MEF90/m_MEF90_EXO.F90 index 4bf91693b1..e2776ea49f 100644 --- a/MEF90/m_MEF90_EXO.F90 +++ b/MEF90/m_MEF90_EXO.F90 @@ -235,7 +235,7 @@ Subroutine EXOGetCellSetElementType_Scal(MEF90Ctx,elemType,ierr) If (MEF90Ctx%rank == 0) Then Call EXGELB(exoID,setID(set),EXOelemType,junk1,junk2,junk3,exoerr) EndIf - Call MPI_Bcast(EXOelemType,MXSTLN,MPI_CHAR,0,MEF90Ctx%comm,ierr) + Call MPI_Bcast(EXOelemType,MXSTLN,MPI_CHARACTER,0,MEF90Ctx%comm,ierr) Call EXO2MEF90ElementType_Scal(EXOelemType,numDim,elemType(set),ierr) End Do DeAllocate(setID) @@ -291,7 +291,7 @@ Subroutine EXOGetCellSetElementType_Vect(MEF90Ctx,elemType,ierr) If (MEF90Ctx%rank == 0) Then Call EXGELB(exoID,setID(set),EXOelemType,junk1,junk2,junk3,exoerr) EndIf - Call MPI_Bcast(EXOelemType,MXSTLN,MPI_CHAR,0,MEF90Ctx%comm,ierr) + Call MPI_Bcast(EXOelemType,MXSTLN,MPI_CHARACTER,0,MEF90Ctx%comm,ierr) Call EXO2MEF90ElementType_Vect(EXOelemType,numDim,elemType(set),ierr) End Do DeAllocate(setID) @@ -347,7 +347,7 @@ Subroutine EXOGetCellSetElementType_Elast(MEF90Ctx,elemType,ierr) If (MEF90Ctx%rank == 0) Then Call EXGELB(exoID,setID(set),EXOelemType,junk1,junk2,junk3,exoerr) EndIf - Call MPI_Bcast(EXOelemType,MXSTLN,MPI_CHAR,0,MEF90Ctx%comm,ierr) + Call MPI_Bcast(EXOelemType,MXSTLN,MPI_CHARACTER,0,MEF90Ctx%comm,ierr) Call EXO2MEF90ElementType_Elast(EXOelemType,numDim,elemType(set),ierr) End Do DeAllocate(setID) From 613beefde92b5a7cfdba6c4091ba542506d68835 Mon Sep 17 00:00:00 2001 From: Jean-Michel Scherer Date: Wed, 22 Jun 2022 11:23:38 -0700 Subject: [PATCH 03/16] make prints compatible with python3 --- bin/exo2exo.py | 39 ++++++++++++++++++++------------------- 1 file changed, 20 insertions(+), 19 deletions(-) diff --git a/bin/exo2exo.py b/bin/exo2exo.py index 6bd8d582e8..c562743d7e 100755 --- a/bin/exo2exo.py +++ b/bin/exo2exo.py @@ -7,6 +7,7 @@ import os.path import numpy as np import pymef90 +import ctypes def exo2exo(fin,fout): @@ -18,7 +19,7 @@ def exo2exo(fin,fout): cell2D = ("TRI","TRI3","TRIANGLE","TRISHELL","TRISHELL3","TRI6","TRISHELL6","QUAD","QUAD4","SHELL","SHELL4","QUAD9","SHELL9") cell3D = ("TETRA","TETRA4","TETRA10","HEX","HEX8","HEX27") - print "Opening {0}".format(fin) + print("Opening {0}".format(fin)) e = exo.exodus(fin,"r") dim = e.num_dimensions() nVertices = e.num_nodes() @@ -26,18 +27,18 @@ def exo2exo(fin,fout): nCSets = e.num_blks() nVSets = e.num_node_sets() - print "\tnumber of dimensions: {0}".format(dim) - print "\tnumber of vertices: {0}".format(nVertices) - print "\tnumber of cells: {0}".format(nCells) - print "\tnumber of cell sets: {0}".format(nCSets) - print "\tnumber of vertex sets: {0}".format(nVSets) + print( "\tnumber of dimensions: {0}".format(dim)) + print( "\tnumber of vertices: {0}".format(nVertices)) + print( "\tnumber of cells: {0}".format(nCells)) + print( "\tnumber of cell sets: {0}".format(nCSets)) + print( "\tnumber of vertex sets: {0}".format(nVSets)) ### Reorder cell sets ### List cells of coDIm0 first then cells of coDim 1 ### - print "Reordering cell sets by increasing co-dimension" + print( "Reordering cell sets by increasing co-dimension") cellsType = [] for set in range(nCSets): setID = e.get_elem_blk_ids()[set] @@ -63,7 +64,7 @@ def exo2exo(fin,fout): ### Find hanging nodes ### ### 1. find largest vertex index in all connectivity tables - print "Removing hanging nodes" + print( "Removing hanging nodes") minV = 100000000 maxV = 0 for set in range(nCSets): @@ -72,7 +73,7 @@ def exo2exo(fin,fout): for v in connect[0]: maxV = max(maxV,v) minV = min(minV,v) - print "\tVertex range: {0}/{1}".format(minV,maxV) + print( "\tVertex range: {0}/{1}".format(minV,maxV)) listedVertices = [False,]*maxV for set in range(nCSets): @@ -85,7 +86,7 @@ def exo2exo(fin,fout): vertexReordering = np.array(range(maxV),dtype=int) for v in range(len(listedVertices)): if not listedVertices[v]: - print "\tvertex {0} is missing.".format(v) + print( "\tvertex {0} is missing.".format(v)) vertexReordering[v:] = vertexReordering[v:]-1 vertexReordering[v] = -1 numMissing += 1 @@ -100,9 +101,9 @@ def exo2exo(fin,fout): ### X,Y,Z=e.get_coords() eout.put_coord_names(["x","y","z"][0:dim]) - fixedX = np.empty(nVertices,dtype=exo.c_double) - fixedY = np.empty(nVertices,dtype=exo.c_double) - fixedZ = np.empty(nVertices,dtype=exo.c_double) + fixedX = np.empty(nVertices,dtype=ctypes.c_double) #exo.c_double) + fixedY = np.empty(nVertices,dtype=ctypes.c_double) #exo.c_double) + fixedZ = np.empty(nVertices,dtype=ctypes.c_double) #exo.c_double) for i in range(nVertices): if listedVertices[i]: @@ -128,9 +129,9 @@ def exo2exo(fin,fout): usedID.append(setFixedID) cellType,numCells,numVertexPerCell,numAttr = e.elem_blk_info(setID) print('Assigning ID {0:4d} to cell set "{1}". \tmef90/vDef name will be cs{0:04}'.format(setFixedID,setName)) - print "\tNumber of cells: {0}".format(numCells) - print "\tCell type: {0}".format(cellType) - print "\tNumber of vertex per cell: {0}".format(numVertexPerCell) + print( "\tNumber of cells: {0}".format(numCells)) + print( "\tCell type: {0}".format(cellType)) + print( "\tNumber of vertex per cell: {0}".format(numVertexPerCell)) eout.put_elem_blk_info(setFixedID,cellType,numCells,numVertexPerCell,1) eout.put_elem_blk_name(setFixedID,setName) @@ -177,18 +178,18 @@ def exo2exo(fin,fout): if __name__ == '__main__': if not (len(sys.argv) == 3): - print "usage: {0} ".format(sys.argv[0]) + print("usage: {0} ".format(sys.argv[0])) sys.exit(-1) fin = sys.argv[-2] fout = sys.argv[-1] if not os.path.isfile(fin): - print "usage: {0} ".format(sys.argv[0]) + print( "usage: {0} ".format(sys.argv[0])) sys.exit(-1) if os.path.exists(fout): if pymef90.confirm("ExodusII file {0} already exists. Overwrite?".format(fout)): os.remove(fout) else: - print "bye!" + print( "bye!") sys.exit(-1) exo2exo(fin,fout) From d64abaff945ec80d5b4d3fbbb3ffc2bc6a98272b Mon Sep 17 00:00:00 2001 From: Jean-Michel Scherer Date: Wed, 22 Jun 2022 11:31:18 -0700 Subject: [PATCH 04/16] add operator for double contraction of 4th order tensors --- MEF90/m_MEF90_LinAlg.F90 | 72 ++++++++++++++++++++++++++++++++++++++++ 1 file changed, 72 insertions(+) diff --git a/MEF90/m_MEF90_LinAlg.F90 b/MEF90/m_MEF90_LinAlg.F90 index e63aa4c69b..3aec00e21e 100644 --- a/MEF90/m_MEF90_LinAlg.F90 +++ b/MEF90/m_MEF90_LinAlg.F90 @@ -2843,6 +2843,78 @@ Function Tens4OS3DTransform(T,M) Call PetscLogFlops(flops,ierr);CHKERRQ(ierr) End Function Tens4OS3DTransform + Function Tens4OS2DXTens4OS2D(T1,T2) + !!! Compute the double contraction of 4th order tensors + !!! i.e. C_{ijkl} = T1_{ijpq}:T2_{pqkl} + !!! (c) 2022 Jean-Michel Scherer scherer@caltech.edu + Type(Tens4OS2D),Intent(IN) :: T1,T2 + + Type(Tens4OS2D) :: Tens4OS2DXTens4OS2D + + PetscReal,Dimension(2,2,2,2) :: TT1,TT2,C + Integer :: i,j,k,l + Integer :: p,q + PetscLogDouble :: flops + PetscErrorCode :: ierr + + TT1 = T1 + TT2 = T2 + C = 0.0_Kr + Do i = 1,2 + Do j = 1,2 + Do k = 1,2 + Do l = 1,2 + Do p = 1,2 + Do q = 1,2 + C(i,j,k,l) = C(i,j,k,l) + TT1(i,j,p,q) * TT2(p,q,k,l) + End Do + End Do + End Do + End Do + End Do + End Do + + Tens4OS2DXTens4OS2D = C + flops = 192 ! 2**6 * 3 + Call PetscLogFlops(flops,ierr);CHKERRQ(ierr) + End Function Tens4OS2DXTens4OS2D + + Function Tens4OS3DXTens4OS3D(T1,T2) + !!! Compute the double contraction of 4th order tensors + !!! i.e. C_{ijkl} = T1_{ijpq}:T2_{pqkl} + !!! (c) 2022 Jean-Michel Scherer scherer@caltech.edu + Type(Tens4OS3D),Intent(IN) :: T1,T2 + + Type(Tens4OS3D) :: Tens4OS3DXTens4OS3D + + PetscReal,Dimension(3,3,3,3) :: TT1,TT2,C + Integer :: i,j,k,l + Integer :: p,q + PetscLogDouble :: flops + PetscErrorCode :: ierr + + TT1 = T1 + TT2 = T2 + C = 0.0_Kr + Do i = 1,3 + Do j = 1,3 + Do k = 1,3 + Do l = 1,3 + Do p = 1,3 + Do q = 1,3 + C(i,j,k,l) = C(i,j,k,l) + TT1(i,j,p,q) * TT2(p,q,k,l) + End Do + End Do + End Do + End Do + End Do + End Do + + Tens4OS3DXTens4OS3D = C + flops = 2187 ! 3**6 * 3 + Call PetscLogFlops(flops,ierr);CHKERRQ(ierr) + End Function Tens4OS3DXTens4OS3D + Function Tens4OS2DSquareRoot(T) Type(Tens4OS2D),Intent(IN) :: T Type(Tens4OS2D) :: Tens4OS2DSquareRoot From eb4dfb74531d40bcc2a1cc55dec066b55ebe7a1d Mon Sep 17 00:00:00 2001 From: Jean-Michel Scherer Date: Wed, 22 Jun 2022 11:59:07 -0700 Subject: [PATCH 05/16] fix incompatibility issue in Hookeslaw2DSum --- m_DefMech/m_MEF90_DefMechSplitNone.F90 | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/m_DefMech/m_MEF90_DefMechSplitNone.F90 b/m_DefMech/m_MEF90_DefMechSplitNone.F90 index 3cb1aad75f..d62d036184 100644 --- a/m_DefMech/m_MEF90_DefMechSplitNone.F90 +++ b/m_DefMech/m_MEF90_DefMechSplitNone.F90 @@ -96,7 +96,7 @@ Subroutine D2EEDNone(self,Strain,HookesLaw,D2EEDPlus,D2EEDMinus) D2EEDPlus = HookesLaw - D2EEDMinus%type = MEF90HookesLawTypeIsotropic + D2EEDMinus%type = HookesLaw%type #if MEF90_DIM==2 D2EEDMinus%isPlaneStress = HookesLaw%isPlaneStress #endif @@ -104,6 +104,18 @@ Subroutine D2EEDNone(self,Strain,HookesLaw,D2EEDPlus,D2EEDMinus) D2EEDMinus%PoissonRatio = 0.0_Kr D2EEDMinus%lambda = 0.0_Kr D2EEDMinus%mu = 0.0_Kr +#if MEF90_DIM==2 + D2EEDMinus%fullTensor = 0.0_Kr * HookesLaw%fullTensor !Tens4OS2D(0.00000_Kr,0.00000_Kr,0.00000_Kr, & ! HookesLaw XXXX,XXYY,XXXY + ! 0.00000_Kr,0.00000_Kr, & ! YYYY,YYXY + ! 0.00000_Kr) & ! XYXY +#elif MEF90_DIM==3 + D2EEDMinus%fullTensor = 0.0_Kr * HookesLaw%fullTensor !Tens4OS3D(0.00000_Kr,0.00000_Kr,0.00000_Kr,0.00000_Kr,0.00000_Kr,0.00000_Kr, & ! XXXX,XXYY,XXZZ,XXYZ,XXXZ,XXXY + !0.00000_Kr,0.00000_Kr,0.00000_Kr,0.00000_Kr,0.00000_Kr, & ! YYYY,YYZZ,YYYZ,YYXZ,YYXY + ! 0.00000_Kr,0.00000_Kr,0.00000_Kr,0.00000_Kr, & ! ZZZZ ZZYZ,ZZXZ,ZZXY + ! 0.00000_Kr,0.00000_Kr,0.00000_Kr, & ! YXYX,YZXZ,YZXY + ! 0.00000_Kr,0.00000_Kr, & ! XZXZ,XZXY + ! 0.00000_Kr) &! XYXY +#endif End Subroutine D2EEDNone End Module MEF90_APPEND(m_MEF90_DefMechSplitNone,MEF90_DIM)D From 9b4440256634860713691b8443798c4cba603856 Mon Sep 17 00:00:00 2001 From: Jean-Michel Scherer Date: Wed, 22 Jun 2022 12:08:44 -0700 Subject: [PATCH 06/16] fix link with snlp --- Makefile.include | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Makefile.include b/Makefile.include index 0e94d8b3b1..4e8b920c52 100644 --- a/Makefile.include +++ b/Makefile.include @@ -8,8 +8,8 @@ endif MEF90_INCLUDE=${MEF90_DIR}/objs/${PETSC_ARCH} ifneq (${SNLP_DIR},) - FC_FLAGS+= -DMEF90_HAVE_SNLP -I${SNLP_DIR}/include - C_FLAGS+= -I${SNLP_DIR}/include + FC_FLAGS+= -DMEF90_HAVE_SNLP -I${SNLP_DIR}/include -I${SNLP_DIR}/obj + C_FLAGS+= -I${SNLP_DIR}/include -I${SNLP_DIR}/obj SNLP_LIB:= -L${SNLP_DIR}/lib -lsnlp -lsnlpF90 -Wl,-rpath,${SNLP_DIR}/lib endif From 641f06aa9989d5acf8e55eec7167dd3cd238235a Mon Sep 17 00:00:00 2001 From: Jean-Michel Scherer Date: Fri, 24 Jun 2022 11:54:30 -0700 Subject: [PATCH 07/16] Implement crystal plasticity in a rate-dependent setting (viscous Norton-like evolution equation). Single slip and Body-Centered Cubic (BCC) behaviours are considered. --- MEF90/m_MEF90_Materials.F90 | 26 ++ m_DefMech/m_MEF90_DefMechCtx.F90 | 14 +- m_DefMech/m_MEF90_DefMechPlasticity.F90 | 356 +++++++++++++++++++++++- 3 files changed, 386 insertions(+), 10 deletions(-) diff --git a/MEF90/m_MEF90_Materials.F90 b/MEF90/m_MEF90_Materials.F90 index c1b22b5cb2..72ce95b1f2 100644 --- a/MEF90/m_MEF90_Materials.F90 +++ b/MEF90/m_MEF90_Materials.F90 @@ -75,6 +75,10 @@ Module m_MEF90_Materials_Types PetscBool :: isLinearIsotropicHardening PetscBool :: isNoPlCoupling Type(MEF90RotationMatrix3D) :: RotationMatrix ! rotation matrix from the global frame to the material frame: X_local = R . X_global + PetscBool :: isViscousPlasticity ! boolean telling if crystal plasticity is viscous or rate-independent + PetscReal :: ViscosityGamma0 ! viscosity reference slip rate + PetscReal :: ViscosityN ! viscosity exponent + PetscReal :: Viscositydt ! time step size Character(len=MEF90_MXSTRLEN) :: Name End Type MEF90MatProp2D_Type @@ -116,6 +120,10 @@ Module m_MEF90_Materials_Types PetscBool :: isLinearIsotropicHardening PetscBool :: isNoPlCoupling Type(MEF90RotationMatrix3D) :: RotationMatrix ! rotation matrix from the global frame to the material frame: X_local = R . X_global + PetscBool :: isViscousPlasticity ! boolean telling if crystal plasticity is viscous or rate-independent + PetscReal :: ViscosityGamma0 ! viscosity reference slip rate + PetscReal :: ViscosityN ! viscosity exponent + PetscReal :: Viscositydt ! time step size Character(len=MEF90_MXSTRLEN) :: Name End Type MEF90MatProp3D_Type @@ -184,6 +192,10 @@ Module m_MEF90_Materials_Types MEF90RotationMatrix3D(MEF90Mat3DIdentity,0.0_Kr,0.0_Kr,0.0_Kr, & ! RotationMatrix, phi1, Phi, phi2, Vect3D(1.0_Kr,0.0_Kr,0.0_Kr),Vect3D(0.0_Kr,1.0_Kr,0.0_Kr), & ! V1, V2, Vect3D(0.0_Kr,0.0_Kr,1.0_Kr),.False.), & ! V3, fromEuler + .FALSE., & ! isViscousPlasticity + 1.0_Kr, & ! ViscosityGamma0 + 1.0_Kr, & ! ViscosityN + 1.0_Kr, & ! Viscositydt "MEF90Mathium2D") Type(MEF90MatProp3D_Type),Parameter :: MEF90Mathium3D = MEF90MatProp3D_Type( & @@ -244,6 +256,10 @@ Module m_MEF90_Materials_Types MEF90RotationMatrix3D(MEF90Mat3DIdentity,0.0_Kr,0.0_Kr,0.0_Kr, & ! RotationMatrix, phi1, Phi, phi2, Vect3D(1.0_Kr,0.0_Kr,0.0_Kr),Vect3D(0.0_Kr,1.0_Kr,0.0_Kr), & ! V1, V2, Vect3D(0.0_Kr,0.0_Kr,1.0_Kr),.False.), & ! V3, fromEuler + .FALSE., & ! isViscousPlasticity + 1.0_Kr, & ! ViscosityGamma0 + 1.0_Kr, & ! ViscosityN + 1.0_Kr, & ! Viscositydt "MEF90Mathium3D") End Module m_MEF90_Materials_Types @@ -547,6 +563,11 @@ Subroutine PetscBagRegisterMEF90MatProp2D(bag,name,prefix,default,ierr) matprop%RotationMatrix%V3 = default%RotationMatrix%V3 Call PetscBagRegisterRealArray(bag,matprop%RotationMatrix%V3,3,'RotationMatrix_V3','[] (V3) Third column of the rotation matrix',ierr) Call PetscBagRegisterBool(bag,matprop%RotationMatrix%fromEuler,default%RotationMatrix%fromEuler,'RotationMatrix_fromEuler','Define rotation matrix from Bunge-Euler angles',ierr);CHKERRQ(ierr) + + Call PetscBagRegisterBool(bag,matprop%isViscousPlasticity,default%isViscousPlasticity,'isViscousPlasticity','[bool] Viscous plastic potential',ierr);CHKERRQ(ierr) + Call PetscBagRegisterReal(bag,matprop%ViscosityGamma0,default%ViscosityGamma0,'ViscosityGamma0','[s^(-1)] Reference plastic deformation rate',ierr);CHKERRQ(ierr) + Call PetscBagRegisterReal(bag,matprop%ViscosityN,default%ViscosityN,'ViscosityN','[unit-less] Viscosity exponent',ierr);CHKERRQ(ierr) + Call PetscBagRegisterReal(bag,matprop%Viscositydt,default%Viscositydt,'Viscositydt','[s] Viscosity time step size',ierr);CHKERRQ(ierr) !Call PetscBagSetFromOptions(bag,ierr) End Subroutine PetscBagRegisterMEF90MatProp2D @@ -633,6 +654,11 @@ Subroutine PetscBagRegisterMEF90MatProp3D(bag,name,prefix,default,ierr) matprop%RotationMatrix%V3 = default%RotationMatrix%V3 Call PetscBagRegisterRealArray(bag,matprop%RotationMatrix%V3,3,'RotationMatrix_V3','[] (V3) Third column of the rotation matrix',ierr) Call PetscBagRegisterBool(bag,matprop%RotationMatrix%fromEuler,default%RotationMatrix%fromEuler,'RotationMatrix_fromEuler','Define rotation matrix from Bunge-Euler angles',ierr);CHKERRQ(ierr) + + Call PetscBagRegisterBool(bag,matprop%isViscousPlasticity,default%isViscousPlasticity,'isViscousPlasticity','[bool] Viscous plastic potential',ierr);CHKERRQ(ierr) + Call PetscBagRegisterReal(bag,matprop%ViscosityGamma0,default%ViscosityGamma0,'ViscosityGamma0','[s^(-1)] Reference plastic deformation rate',ierr);CHKERRQ(ierr) + Call PetscBagRegisterReal(bag,matprop%ViscosityN,default%ViscosityN,'ViscosityN','[unit-less] Viscosity exponent',ierr);CHKERRQ(ierr) + Call PetscBagRegisterReal(bag,matprop%Viscositydt,default%Viscositydt,'Viscositydt','[s] Viscosity time step size',ierr);CHKERRQ(ierr) !Call PetscBagSetFromOptions(bag,ierr) End Subroutine PetscBagRegisterMEF90MatProp3D diff --git a/m_DefMech/m_MEF90_DefMechCtx.F90 b/m_DefMech/m_MEF90_DefMechCtx.F90 index 4047d7f30e..650a4d8d4c 100644 --- a/m_DefMech/m_MEF90_DefMechCtx.F90 +++ b/m_DefMech/m_MEF90_DefMechCtx.F90 @@ -289,9 +289,11 @@ Module m_MEF90_DefMechCtx MEF90DefMech_plasticityTypeVonMises1D, & MEF90DefMech_plasticityTypeHillPlaneTheory, & MEF90DefMech_PlasticityTypeGreen, & - MEF90DefMech_PlasticityTypeGurson + MEF90DefMech_PlasticityTypeGurson, & + MEF90DefMech_PlasticityTypeCrystalSingleSlip, & + MEF90DefMech_PlasticityTypeCrystalBCC End Enum - Character(len = MEF90_MXSTRLEN),Dimension(13),protected :: MEF90DefMech_plasticityTypeList + Character(len = MEF90_MXSTRLEN),Dimension(15),protected :: MEF90DefMech_plasticityTypeList Enum,bind(c) enumerator :: MEF90DefMech_unilateralContactTypeNone = 0, & @@ -377,9 +379,11 @@ Subroutine MEF90DefMechCtxInitialize_Private(ierr) MEF90DefMech_plasticityTypeList(8) = 'HillPlaneTheory' MEF90DefMech_plasticityTypeList(9) = 'Green' MEF90DefMech_plasticityTypeList(10) = 'Gurson' - MEF90DefMech_plasticityTypeList(11) = 'MEF90DefMech_plasticityType' - MEF90DefMech_plasticityTypeList(12) = '_MEF90DefMech_plasticityType' - MEF90DefMech_plasticityTypeList(13) = '' + MEF90DefMech_plasticityTypeList(11) = 'CrystalSingleSlip' + MEF90DefMech_plasticityTypeList(12) = 'CrystalBCC' + MEF90DefMech_plasticityTypeList(13) = 'MEF90DefMech_plasticityType' + MEF90DefMech_plasticityTypeList(14) = '_MEF90DefMech_plasticityType' + MEF90DefMech_plasticityTypeList(15) = '' MEF90DefMech_unilateralContactTypeList(1) = 'None' MEF90DefMech_unilateralContactTypeList(2) = 'HydrostaticDeviatoric' diff --git a/m_DefMech/m_MEF90_DefMechPlasticity.F90 b/m_DefMech/m_MEF90_DefMechPlasticity.F90 index ba0e7c4214..0e581d89a6 100644 --- a/m_DefMech/m_MEF90_DefMechPlasticity.F90 +++ b/m_DefMech/m_MEF90_DefMechPlasticity.F90 @@ -42,6 +42,12 @@ Module MEF90_APPEND(m_MEF90_DefMechPlasticity,MEF90_DIM)D real(Kind = Kr) :: Phi real(Kind = Kr) :: delta logical(Kind = Kr) :: isNoPlCoupling + Type(MEF90RotationMatrix3D) :: RotationMatrix3D + logical(Kind = Kr) :: isViscousPlasticity + real(Kind = Kr) :: ViscosityGamma0 + real(Kind = Kr) :: ViscosityN + real(Kind = Kr) :: Viscositydt + real(Kind = Kr) :: viscouscumulatedDissipatedPlasticEnergyVariation end type MEF90DefMechPlasticityCtx contains @@ -68,6 +74,9 @@ Module MEF90_APPEND(m_MEF90_DefMechPlasticity,MEF90_DIM)D #define FHG_GURSON MEF90_APPEND(fhg_Gurson,MEF90_DIM)D +#define FHG_CRYSTALSINGLESLIP MEF90_APPEND(fhg_CrystalSingleSlip,MEF90_DIM)D + +#define FHG_CRYSTALBCC MEF90_APPEND(fhg_CrystalBCC,MEF90_DIM)D #undef __FUNCT__ #define __FUNCT__ "FHG_NONE" @@ -706,6 +715,300 @@ subroutine FHG_GURSON(x,f,h,g,myctx) bind(c) end subroutine FHG_GURSON +#undef __FUNCT__ +#define __FUNCT__ "FHG_CRYSTALSINGLESLIP" +!!! +!!! +!!! fhg: Single slip Crystal Plasticity: m=[100] and n=[010] +!!! +!!! (c) 2022 Jean-Michel Scherer, Caltech , +!!! Blaise Bourdin, LSU, +!!! +!!! + + subroutine FHG_CRYSTALSINGLESLIP(x,f,h,g,myctx) bind(c) + use,intrinsic :: iso_c_binding + use m_MEF90 + + real(kind=c_double) :: x(*) + real(kind=c_double) :: f(*) + real(kind=c_double) :: h(*) + real(kind=c_double) :: g(*) + real(kind = Kr) :: StiffnessA,StiffnessB + + type(c_ptr),intent(in),value :: myctx + type(MEF90DefMechPlasticityCtx),pointer :: myctx_ptr + type(MEF90_MATS) :: xMatS + type(MEF90_MATS) :: InelasticStrain, PlasticStrainFlow, Stress + type(MatS3D) :: Strain3D, PlasticStrainFlow3D, Stress3D, Stress3DCrystal, MatrixMu, TotalPlasticIncrement, TotalPlasticIncrementCrystal + type(MAT3D) :: MatrixR + real(Kind = Kr),dimension(1) :: ResolvedShearStress,PlasticSlipIncrement + integer(Kind = Kr) :: s + type(Vect3D) :: n + type(Vect3D) :: m + real(Kind = Kr) :: normS, dt, fixed_point_norm + real(Kind = Kr),dimension(3,1) :: n_s = reshape((/ 0., 1., 0./), (/3,1/)) !!! slip planes normals + real(Kind = Kr),dimension(3,1) :: m_s = reshape((/ 1., 0., 0./), (/3,1/)) !!! slip directions + real(Kind = Kr) :: lambda,mu,E,nu + + !!! Casting x into a MEF90_MATS + xMatS = x(1:SIZEOFMEF90_MATS) + !!! This is the fortran equivalent of casting ctx into a c_ptr + call c_f_pointer(myctx,myctx_ptr) + + !!! Select which softening young model + if (myctx_ptr%CoefficientLinSoft==0) then + StiffnessA = (1.0_Kr - myctx_ptr%Damage)**2 + myctx_ptr%residualStiffness + StiffnessB = (1.0_Kr - myctx_ptr%Damage)**myctx_ptr%DuctileCouplingPower + myctx_ptr%residualStiffness + else + StiffnessA = ( (1.0_Kr - myctx_ptr%Damage)**2 /( 1.0_Kr + ( myctx_ptr%CoefficientLinSoft - 1.0_Kr )*(1.0_Kr - (1.0_Kr - myctx_ptr%Damage)**2 ) ) ) + myctx_ptr%residualStiffness + StiffnessB = (1.0_Kr - myctx_ptr%Damage)**myctx_ptr%DuctileCouplingPower + myctx_ptr%residualStiffness + endif + + if (myctx_ptr%isNoPlCoupling .eqv. .true.) then + StiffnessB = 1.0_Kr + else + StiffnessB = ( (1.0_Kr-myctx_ptr%residualYieldStress)*StiffnessB + myctx_ptr%residualYieldStress ) + endif + + PlasticStrainFlow = xMatS-myctx_ptr%PlasticStrainOld + Stress = (myctx_ptr%HookesLaw*(myctx_ptr%InelasticStrain-xMatS)) !*StiffnessA + +#if MEF90_DIM==2 + !!! If plane strain + E = myctx_ptr%HookesLaw%YoungsModulus + nu = myctx_ptr%HookesLaw%PoissonRatio + mu = E / (1.0_Kr + nu) * .5_Kr + lambda = E * nu / (1.0_Kr + nu) / (1 - 2.0_Kr * nu) + + Strain3D = 0.0_Kr + Strain3D%XX = myctx_ptr%InelasticStrain%XX + Strain3D%YY = myctx_ptr%InelasticStrain%YY + Strain3D%XY = myctx_ptr%InelasticStrain%XY + + PlasticStrainFlow3D = 0.0_Kr + PlasticStrainFlow3D%XX = xMatS%XX - myctx_ptr%PlasticStrainOld%XX + PlasticStrainFlow3D%YY = xMatS%YY - myctx_ptr%PlasticStrainOld%YY + PlasticStrainFlow3D%XY = xMatS%XY - myctx_ptr%PlasticStrainOld%XY + PlasticStrainFlow3D%ZZ = -( PlasticStrainFlow3D%XX + PlasticStrainFlow3D%YY ) + + Stress3D = 0.0_Kr + Stress3D%XX = Stress%XX + Stress3D%XY = Stress%XY + Stress3D%YY = Stress%YY + Stress3D%ZZ = lambda*(Trace(Strain3D)) + 2*mu*(Strain3D%ZZ+Trace(xMatS)) +#elif MEF90_DIM==3 + Stress3D = Stress + Strain3D = myctx_ptr%InelasticStrain + PlasticStrainFlow3D = xMatS - myctx_ptr%PlasticStrainOld +#endif + + Stress3DCrystal = MatRaRt(Stress3D,myctx_ptr%RotationMatrix3D%fullTensor) + + TotalPlasticIncrementCrystal = 0.0_Kr + TotalPlasticIncrement = 0.0_Kr + myctx_ptr%viscouscumulatedDissipatedPlasticEnergyVariation = 0.0_Kr + + normS = (2.0_Kr*SQRT(1.0_Kr)) + + dt = myctx_ptr%Viscositydt + + f(1) = 0.5_Kr * StiffnessA * Stress .DotP. (myctx_ptr%InelasticStrain-xMatS) + Do s = 1,1 + m = m_s(:,s) + n = n_s(:,s) + MatrixMu = ((m .TensP. n) + (n .TensP. m)) / normS + ResolvedShearStress(s) = Stress3DCrystal .DotP. MatrixMu + if (myctx_ptr%isViscousPlasticity) then + PlasticSlipIncrement(s) = dt * myctx_ptr%ViscosityGamma0 * SIGN(1.0_Kr, ResolvedShearStress(s)) *& + & MAX( (ABS(StiffnessA*ResolvedShearStress(s)) - StiffnessB*myctx_ptr%YieldTau0) / myctx_ptr%YieldTau0 , 0. )**myctx_ptr%ViscosityN + TotalPlasticIncrementCrystal = TotalPlasticIncrementCrystal + (PlasticSlipIncrement(s) * MatrixMu) + !myctx_ptr%plasticSlipsVariation(s) = PlasticSlipIncrement(s) + myctx_ptr%viscouscumulatedDissipatedPlasticEnergyVariation = myctx_ptr%viscouscumulatedDissipatedPlasticEnergyVariation + ResolvedShearStress(s)*PlasticSlipIncrement(s) + else + print *, "Rate-independent crystal plasticity is not implemented." + !PlasticSlipIncrement(s) = dt * myctx_ptr%eta * SIGN(1.0_Kr, ResolvedShearStress(s)) *& + ! & MAX( (ABS(StiffnessA*ResolvedShearStress(s)) - StiffnessB*myctx_ptr%YieldTau0) / myctx_ptr%YieldTau0 , 0. )**myctx_ptr%ViscosityN + !TotalPlasticIncrementCrystal = TotalPlasticIncrementCrystal + (PlasticSlipIncrement(s) * MatrixMu) + !myctx_ptr%viscouscumulatedDissipatedPlasticEnergyVariation = myctx_ptr%viscouscumulatedDissipatedPlasticEnergyVariation + ResolvedShearStress(s)*PlasticSlipIncrement(s) + endif + End Do + f(1) = f(1) + StiffnessB*myctx_ptr%viscouscumulatedDissipatedPlasticEnergyVariation + TotalPlasticIncrement = MatRtaR(TotalPlasticIncrementCrystal,myctx_ptr%RotationMatrix3D%fullTensor) +#if MEF90_DIM==2 + h(1) = PlasticStrainFlow3D%XX - TotalPlasticIncrement%XX + h(2) = PlasticStrainFlow3D%XY - TotalPlasticIncrement%XY + h(3) = PlasticStrainFlow3D%YY - TotalPlasticIncrement%YY +#elif MEF90_DIM==3 + h(1) = NORM(PlasticStrainFlow3D - TotalPlasticIncrement) +#endif + end subroutine FHG_CRYSTALSINGLESLIP + + +#undef __FUNCT__ +#define __FUNCT__ "FHG_CRYSTALBCC" +!!! +!!! +!!! fhg: BCC Crystal Plasticity (octahedral): equivalent to FCC in small strains +!!! +!!! (c) 2022 Jean-Michel Scherer, Caltech , +!!! Blaise Bourdin, LSU, +!!! +!!! + + subroutine FHG_CRYSTALBCC(x,f,h,g,myctx) bind(c) + use,intrinsic :: iso_c_binding + use m_MEF90 + + real(kind=c_double) :: x(*) + real(kind=c_double) :: f(*) + real(kind=c_double) :: h(*) + real(kind=c_double) :: g(*) + real(kind = Kr) :: StiffnessA,StiffnessB + + type(c_ptr),intent(in),value :: myctx + type(MEF90DefMechPlasticityCtx),pointer :: myctx_ptr + type(MEF90_MATS) :: xMatS + type(MEF90_MATS) :: InelasticStrain,PlasticStrainFlow,Stress + type(MatS3D) :: Strain3D,PlasticStrainFlow3D,Stress3D,Stress3DCrystal,TotalPlasticIncrement,TotalPlasticIncrementCrystal !, MatrixMu + type(MatS3D) :: PlasticStrain3D,PlasticStrain3DCrystal + type(MatS3D),dimension(12) :: MatrixMu + type(MAT3D) :: MatrixR + real(Kind = Kr),dimension(12) :: ResolvedShearStress,PlasticSlipIncrement + real(Kind = Kr),dimension(12) :: PlasticSlips + integer(Kind = Kr) :: s,k,active + type(Vect3D) :: n + type(Vect3D) :: m + real(Kind = Kr) :: normS,dt,CRSS + !!! slip systems nomenclature: + !!! Bd,B4 Ba,B2 Bc,B5 Db,D4 Dc,D1 Da,D6 + !!! Ab,A2 Ad,A6 Ac,A3 Cb,C5 Ca,C3 Cd,C1 + real(Kind = Kr),dimension(3,12) :: n_s = reshape((/-1., 0., 1., 0.,-1., 1., -1., 1., 0., -1., 0., 1., 0., 1., 1., 1., 1., 0., & + & 0.,-1., 1., 1., 1., 0., 1., 0., 1., -1., 1., 0., 1., 0., 1., 0., 1., 1. /), (/3,12/)) !!! slip planes normals + real(Kind = Kr),dimension(3,12) :: m_s = reshape((/ 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,-1., 1., 1.,-1., 1., 1.,-1., 1., & + &-1., 1., 1., -1., 1., 1., -1., 1., 1., -1.,-1., 1., -1.,-1., 1., -1.,-1., 1. /), (/3,12/)) !!! slip directions + real(Kind = Kr) :: lambda,mu,E,nu + + !!! Casting x into a MEF90_MATS + xMatS = x(1:SIZEOFMEF90_MATS) + !!! This is the fortran equivalent of casting ctx into a c_ptr + call c_f_pointer(myctx,myctx_ptr) + + !!! Select which softening young model + if (myctx_ptr%CoefficientLinSoft==0) then + StiffnessA = (1.0_Kr - myctx_ptr%Damage)**2 + myctx_ptr%residualStiffness + StiffnessB = (1.0_Kr - myctx_ptr%Damage)**myctx_ptr%DuctileCouplingPower + myctx_ptr%residualStiffness + else + StiffnessA = ( (1.0_Kr - myctx_ptr%Damage)**2 /( 1.0_Kr + ( myctx_ptr%CoefficientLinSoft - 1.0_Kr )*(1.0_Kr - (1.0_Kr - myctx_ptr%Damage)**2 ) ) ) + myctx_ptr%residualStiffness + StiffnessB = (1.0_Kr - myctx_ptr%Damage)**myctx_ptr%DuctileCouplingPower + myctx_ptr%residualStiffness + endif + + if (myctx_ptr%isNoPlCoupling) then + StiffnessB = 1.0_Kr + else + StiffnessB = ( (1.0_Kr-myctx_ptr%residualYieldStress)*StiffnessB + myctx_ptr%residualYieldStress ) + endif + + PlasticStrainFlow = xMatS-myctx_ptr%PlasticStrainOld + Stress = (myctx_ptr%HookesLaw*(myctx_ptr%InelasticStrain-xMatS)) + +#if MEF90_DIM==2 + !!! If plane strain + E = myctx_ptr%HookesLaw%YoungsModulus + nu = myctx_ptr%HookesLaw%PoissonRatio + mu = E / (1.0_Kr + nu) * .5_Kr + lambda = E * nu / (1.0_Kr + nu) / (1 - 2.0_Kr * nu) + + Strain3D = 0.0_Kr + Strain3D%XX = myctx_ptr%InelasticStrain%XX + Strain3D%YY = myctx_ptr%InelasticStrain%YY + Strain3D%XY = myctx_ptr%InelasticStrain%XY + + PlasticStrainFlow3D = 0.0_Kr + PlasticStrainFlow3D%XX = xMatS%XX - myctx_ptr%PlasticStrainOld%XX + PlasticStrainFlow3D%YY = xMatS%YY - myctx_ptr%PlasticStrainOld%YY + PlasticStrainFlow3D%XY = xMatS%XY - myctx_ptr%PlasticStrainOld%XY + PlasticStrainFlow3D%ZZ = -( PlasticStrainFlow3D%XX + PlasticStrainFlow3D%YY ) + + Stress3D = 0.0_Kr + Stress3D%XX = Stress%XX + Stress3D%XY = Stress%XY + Stress3D%YY = Stress%YY + Stress3D%ZZ = lambda*(Trace(Strain3D)) + 2*mu*(Strain3D%ZZ+Trace(xMatS)) + + PlasticStrain3D = 0.0_Kr + PlasticStrain3D%XX = xMatS%XX + PlasticStrain3D%YY = xMatS%YY + PlasticStrain3D%XY = xMatS%XY + PlasticStrain3D%ZZ = -(PlasticStrain3D%XX + PlasticStrain3D%YY) +#elif MEF90_DIM==3 + Stress3D = Stress + Strain3D = myctx_ptr%InelasticStrain + PlasticStrainFlow3D = xMatS - myctx_ptr%PlasticStrainOld + PlasticStrain3D = xMatS +#endif + + Stress3DCrystal = MatRaRt(Stress3D,myctx_ptr%RotationMatrix3D%fullTensor) + + normS = (2.0_Kr*SQRT(6.0_Kr)) + Do s=1,12 + m = m_s(:,s) + n = n_s(:,s) + MatrixMu(s) = ((m .TensP. n) + (n .TensP. m)) / normS + end do + + !if ( .NOT. (myctx_ptr%YieldQ==0.0_Kr) ) then + ! PlasticStrain3DCrystal = MatRaRt(PlasticStrain3D,myctx_ptr%RotationMatrix3D%fullTensor) + ! PlasticSlips(s) = PlasticStrain3DCrystal .DotP. MatrixMu(s) + !end if + + TotalPlasticIncrementCrystal = 0.0_Kr + TotalPlasticIncrement = 0.0_Kr + myctx_ptr%viscouscumulatedDissipatedPlasticEnergyVariation = 0.0_Kr + + dt = myctx_ptr%Viscositydt + active = 0 + + f(1) = 0.5_Kr * StiffnessA * Stress .DotP. (myctx_ptr%InelasticStrain-xMatS) + Do s = 1,12 + ResolvedShearStress(s) = Stress3DCrystal .DotP. MatrixMu(s) + CRSS = myctx_ptr%YieldTau0 + !if ( .NOT. (myctx_ptr%YieldQ==0.0_Kr) ) then + ! Do k=1,12 + ! CRSS = CRSS + myctx_ptr%YieldQ * myctx_ptr%InteractionMatrix%him(s,k) * (1.0_Kr - EXP(-myctx_ptr%Yieldb * ABS(PlasticSlips(k)) )) + ! end do + !end if + if (myctx_ptr%isViscousPlasticity) then + PlasticSlipIncrement(s) = dt * myctx_ptr%ViscosityGamma0 * SIGN(1.0_Kr, ResolvedShearStress(s)) *& + & MAX( (ABS(StiffnessA*ResolvedShearStress(s)) - StiffnessB*CRSS) / myctx_ptr%YieldTau0 , 0. )**myctx_ptr%ViscosityN + TotalPlasticIncrementCrystal = TotalPlasticIncrementCrystal + (PlasticSlipIncrement(s) * MatrixMu(s)) + !myctx_ptr%plasticSlipsVariation(s) = PlasticSlipIncrement(s) + myctx_ptr%viscouscumulatedDissipatedPlasticEnergyVariation = myctx_ptr%viscouscumulatedDissipatedPlasticEnergyVariation + ResolvedShearStress(s)*PlasticSlipIncrement(s) + print *,'CRSS = ',CRSS,' ', 'taus = ',ResolvedShearStress(s) + print *,'gammas = ',PlasticSlipIncrement(s) + else + print *, "Rate-independent crystal plasticity is not implemented." + !PlasticSlipIncrement(s) = dt * myctx_ptr%eta * SIGN(1.0_Kr, ResolvedShearStress(s)) *& + ! & MAX( (ABS(StiffnessA*ResolvedShearStress(s)) - StiffnessB*CRSS) / myctx_ptr%YieldTau0 , 0. )**myctx_ptr%ViscosityN + !TotalPlasticIncrementCrystal = TotalPlasticIncrementCrystal + (PlasticSlipIncrement(s) * MatrixMu(s)) + !myctx_ptr%plasticSlipsVariation(s) = PlasticSlipIncrement(s) + !myctx_ptr%viscouscumulatedDissipatedPlasticEnergyVariation = myctx_ptr%viscouscumulatedDissipatedPlasticEnergyVariation + ResolvedShearStress(s)*PlasticSlipIncrement(s) + endif + if ((ABS(StiffnessA*ResolvedShearStress(s)) - StiffnessB*CRSS)>0) then + active = active + 1 + end if + End Do + f(1) = f(1) + StiffnessB * myctx_ptr%viscouscumulatedDissipatedPlasticEnergyVariation + TotalPlasticIncrement = MatRtaR(TotalPlasticIncrementCrystal,myctx_ptr%RotationMatrix3D%fullTensor) +#if MEF90_DIM==2 + h(1) = PlasticStrainFlow3D%XX - TotalPlasticIncrement%XX + h(2) = PlasticStrainFlow3D%XY - TotalPlasticIncrement%XY + h(3) = PlasticStrainFlow3D%YY - TotalPlasticIncrement%YY +#elif MEF90_DIM==3 + h(1) = NORM(PlasticStrainFlow3D - TotalPlasticIncrement) +#endif + end subroutine FHG_CRYSTALBCC + #undef __FUNCT__ #define __FUNCT__ "MEF90DefMechPlasticStrainUpdate" @@ -891,6 +1194,44 @@ Subroutine MEF90DefMechPlasticStrainUpdate(MEF90DefMechCtx,plasticStrain,x,Plast snlp_m = 0 snlp_p = 0 snlp_ctx = c_loc(PlasticityCtx) + + case(MEF90DefMech_plasticityTypeCrystalSingleSlip) + snlp_Dfhg = c_null_funptr + snlp_fhg = c_funloc(FHG_CRYSTALSINGLESLIP) + snlp_n = SIZEOFMEF90_MATS +#if MEF90_DIM==2 + if (matPropSet%HookesLaw%isPlaneStress) then + write(*,*) "Plane stress CrystalSingleSlip is not implemented" + end if + snlp_m = 3 +#elif MEF90_DIM==3 + snlp_m = 1 +#endif + if (matpropSet%isViscousPlasticity) then + snlp_p = 0 + else + snlp_p = 0 + end if + snlp_ctx = c_loc(PlasticityCtx) + + case(MEF90DefMech_plasticityTypeCrystalBCC) + snlp_Dfhg = c_null_funptr + snlp_fhg = c_funloc(FHG_CRYSTALBCC) + snlp_n = SIZEOFMEF90_MATS +#if MEF90_DIM==2 + if (matPropSet%HookesLaw%isPlaneStress) then + write(*,*) "Plane stress CrystalBCC is not implemented" + end if + snlp_m = 3 +#elif MEF90_DIM==3 + snlp_m = 1 +#endif + if (matpropSet%isViscousPlasticity) then + snlp_p = 0 + else + snlp_p = 0 + end if + snlp_ctx = c_loc(PlasticityCtx) case default Print*,__FUNCT__,': Unimplemented plasticity Type',cellSetOptions%PlasticityType STOP @@ -923,9 +1264,11 @@ Subroutine MEF90DefMechPlasticStrainUpdate(MEF90DefMechCtx,plasticStrain,x,Plast PlasticityCtx%phi1 = matpropSet%phi1 PlasticityCtx%phi2 = matpropSet%phi2 PlasticityCtx%Phi = matpropSet%Phi - - - + PlasticityCtx%RotationMatrix3D = matpropSet%RotationMatrix + PlasticityCtx%isViscousPlasticity = matpropSet%isViscousPlasticity + PlasticityCtx%ViscosityGamma0 = matpropSet%ViscosityGamma0 + PlasticityCtx%ViscosityN = matpropSet%ViscosityN + PlasticityCtx%Viscositydt = matpropSet%Viscositydt #if MEF90_DIM == 2 @@ -999,8 +1342,11 @@ Subroutine MEF90DefMechPlasticStrainUpdate(MEF90DefMechCtx,plasticStrain,x,Plast Print*,__FUNCT__,': Unimplemented damage Type, only AT1Elastic and AT2Elastic implement',cellSetOptions%damageType STOP End Select - cumulatedDissipatedPlasticEnergyVariationLoc(1) = Stiffness * ( PlasticityCtx%HookesLaw * ( PlasticityCtx%InelasticStrain - PlasticStrainMatS ) ) .dotP. ( PlasticStrainMatS - PlasticityCtx%plasticStrainOld ) - + if (.NOT. PlasticityCtx%isViscousPlasticity) then + cumulatedDissipatedPlasticEnergyVariationLoc(1) = Stiffness * ( PlasticityCtx%HookesLaw * ( PlasticityCtx%InelasticStrain - PlasticStrainMatS ) ) .dotP. ( PlasticStrainMatS - PlasticityCtx%plasticStrainOld ) + else + cumulatedDissipatedPlasticEnergyVariationLoc(1) = Stiffness * PlasticityCtx%viscousCumulatedDissipatedPlasticEnergyVariation + endif #if MEF90_DIM == 2 if (PlasticityCtx%isPlaneStress .eqv. .FALSE.) then From 8e22bfe6edd171cba701267a74ed47173a33e315 Mon Sep 17 00:00:00 2001 From: Jean-Michel Scherer Date: Fri, 24 Jun 2022 19:31:43 -0700 Subject: [PATCH 08/16] make surfing BCs implementation match the notation used in the papers --- bin/vDefSurfingBC.py | 22 +++++++++++++++------- 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/bin/vDefSurfingBC.py b/bin/vDefSurfingBC.py index 7f32df7ada..62edbf7b42 100755 --- a/bin/vDefSurfingBC.py +++ b/bin/vDefSurfingBC.py @@ -16,12 +16,14 @@ def parse(args=None): parser.add_argument("--time_numstep",type=int,help="number of time step",default=11) parser.add_argument("--E",type=float,help="Youngs modulus",default=1.) parser.add_argument("--nu",type=float,help="Poisson Ratio",default=.3) + parser.add_argument("--Gc",type=float,help="Critical Energy Release Rate",default=1.) parser.add_argument("--ampl",type=float,help="Amplification",default=1.) parser.add_argument("--initialpos",type=float,nargs=3,help="initial crack tip postion",default=[0.,0,0]) parser.add_argument("--cs",type=int,nargs='*',help="list of cell sets where surfing boundary displacements are applied",default=[]) parser.add_argument("--vs",type=int,nargs='*',help="list of vertex sets where surfing boundary displacements are applied",default=[]) parser.add_argument("--plasticity",default=False,action="store_true",help="Add extended variables for plasticity related fields") parser.add_argument("--force",action="store_true",default=False,help="Overwrite existing files without prompting") + parser.add_argument("--planestrain",action="store_true",default=False,help="Apply plane strain surfing BCs: kappa=3-4*nu") return parser.parse_args() def exoformat(e,plasticity=False): @@ -67,11 +69,17 @@ def cart2polar(x, y): theta = np.arctan2(y, x) return r, theta -def surfingBC(e,t,Xc,cslist,vslist,E,nu,ampl): +def surfingBC(e,t,Xc,cslist,vslist,E,nu,Gc,ampl,planestrain): import numpy as np - kappa = (3.0-nu)/(1.0+nu) + if planestrain: + kappa = (3.0-nu)/(1.0+nu) + Ep = E / (1. - nu**2) + else: + kappa = 3.0-4.0*nu + Ep = E mu = E / (1. + nu) * .5 + KI = np.sqrt(Gc * Ep) dim = e.num_dimensions() X,Y,Z=e.get_coords() @@ -85,16 +93,16 @@ def surfingBC(e,t,Xc,cslist,vslist,E,nu,ampl): for v in vertices: r,theta = cart2polar(X[v-1]-Xc[0]-t,Y[v-1]-Xc[1]) z = Z[v-1]-Xc[2] - U[0,v-1] = ampl * np.sqrt(r / np.pi * .5) / mu * .5 * np.cos(theta * .5) * (kappa - np.cos(theta)) - U[1,v-1] = ampl * np.sqrt(r / np.pi * .5) / mu * .5 * np.sin(theta * .5) * (kappa - np.cos(theta)) + U[0,v-1] = ampl * KI * np.sqrt(r / np.pi * .5) / mu * .5 * np.cos(theta * .5) * (kappa - np.cos(theta)) + U[1,v-1] = ampl * KI * np.sqrt(r / np.pi * .5) / mu * .5 * np.sin(theta * .5) * (kappa - np.cos(theta)) U[2,v-1] = 0.0 for set in vslist: for v in e.get_node_set_nodes(set): r,theta = cart2polar(X[v-1]-Xc[0]-t,Y[v-1]-Xc[1]) z = Z[v-1] - U[0,v-1] = ampl * np.sqrt(r / np.pi * .5) / mu * .5 * np.cos(theta * .5) * (kappa - np.cos(theta)) - U[1,v-1] = ampl * np.sqrt(r / np.pi * .5) / mu * .5 * np.sin(theta * .5) * (kappa - np.cos(theta)) + U[0,v-1] = ampl * KI * np.sqrt(r / np.pi * .5) / mu * .5 * np.cos(theta * .5) * (kappa - np.cos(theta)) + U[1,v-1] = ampl * KI * np.sqrt(r / np.pi * .5) / mu * .5 * np.sin(theta * .5) * (kappa - np.cos(theta)) U[2,v-1] = 0.0 return U @@ -132,7 +140,7 @@ def main(): for t in np.linspace(options.time_min,options.time_max,options.time_numstep): print ("writing step {0}, t = {1:0.4f}".format(step+1,t)) exoout.put_time(step+1,t) - U = surfingBC(exoout,t,options.initialpos,options.cs,options.vs,options.E,options.nu,options.ampl) + U = surfingBC(exoout,t,options.initialpos,options.cs,options.vs,options.E,options.nu,options.Gc,options.ampl,options.planestrain) X,Y,Z=exoout.get_coords() exoout.put_node_variable_values("Displacement_X",step+1,U[0,:]) exoout.put_node_variable_values("Displacement_Y",step+1,U[1,:]) From 76bb4872844768d03a8413e6435a43eb3967cecf Mon Sep 17 00:00:00 2001 From: Jean-Michel Scherer Date: Thu, 30 Jun 2022 18:03:45 -0700 Subject: [PATCH 09/16] remove debugging comments --- m_DefMech/m_MEF90_DefMechPlasticity.F90 | 2 -- 1 file changed, 2 deletions(-) diff --git a/m_DefMech/m_MEF90_DefMechPlasticity.F90 b/m_DefMech/m_MEF90_DefMechPlasticity.F90 index 0e581d89a6..96917e321c 100644 --- a/m_DefMech/m_MEF90_DefMechPlasticity.F90 +++ b/m_DefMech/m_MEF90_DefMechPlasticity.F90 @@ -984,8 +984,6 @@ subroutine FHG_CRYSTALBCC(x,f,h,g,myctx) bind(c) TotalPlasticIncrementCrystal = TotalPlasticIncrementCrystal + (PlasticSlipIncrement(s) * MatrixMu(s)) !myctx_ptr%plasticSlipsVariation(s) = PlasticSlipIncrement(s) myctx_ptr%viscouscumulatedDissipatedPlasticEnergyVariation = myctx_ptr%viscouscumulatedDissipatedPlasticEnergyVariation + ResolvedShearStress(s)*PlasticSlipIncrement(s) - print *,'CRSS = ',CRSS,' ', 'taus = ',ResolvedShearStress(s) - print *,'gammas = ',PlasticSlipIncrement(s) else print *, "Rate-independent crystal plasticity is not implemented." !PlasticSlipIncrement(s) = dt * myctx_ptr%eta * SIGN(1.0_Kr, ResolvedShearStress(s)) *& From 99043c147b430a4b2dbb3b34290118e61153b496 Mon Sep 17 00:00:00 2001 From: Jean-Michel Scherer Date: Sat, 2 Jul 2022 18:25:23 -0700 Subject: [PATCH 10/16] launch visit without window when computing J-integral --- bin/JIntegral2Dx2017.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/bin/JIntegral2Dx2017.py b/bin/JIntegral2Dx2017.py index e99ada2c95..0a5cbb7df6 100755 --- a/bin/JIntegral2Dx2017.py +++ b/bin/JIntegral2Dx2017.py @@ -1,3 +1,5 @@ +#!/usr/bin/env python2 + from visit import * def parse(args=None): import argparse @@ -28,6 +30,8 @@ def plot(opts): import shutil import math + AddArgument("-nowin") + Launch() ## ## Open the database ## From 0a57dbe2358ad5cdefd39824d0b74fe171a30f2b Mon Sep 17 00:00:00 2001 From: Jean-Michel Scherer Date: Mon, 4 Jul 2022 17:17:39 -0700 Subject: [PATCH 11/16] fix bug on boolean isNoPlCoupling which was not initialized (randomly True or False) --- m_DefMech/m_MEF90_DefMechPlasticity.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/m_DefMech/m_MEF90_DefMechPlasticity.F90 b/m_DefMech/m_MEF90_DefMechPlasticity.F90 index 96917e321c..176371ff3e 100644 --- a/m_DefMech/m_MEF90_DefMechPlasticity.F90 +++ b/m_DefMech/m_MEF90_DefMechPlasticity.F90 @@ -1244,6 +1244,7 @@ Subroutine MEF90DefMechPlasticStrainUpdate(MEF90DefMechCtx,plasticStrain,x,Plast PlasticityCtx%residualStiffness = matpropSet%residualStiffness PlasticityCtx%YieldStress = matpropSet%YieldStress PlasticityCtx%DuctileCouplingPower = matpropSet%DuctileCouplingPower + PlasticityCtx%isNoPlCoupling = matpropSet%isNoPlCoupling PlasticityCtx%CoefficientDruckerPrager = matpropSet%CoefficientDruckerPrager PlasticityCtx%CoefficientCapModel0 = matpropSet%CoefficientCapModel0 PlasticityCtx%CoefficientCapModel1 = matpropSet%CoefficientCapModel1 From d5da9ed61d0d9dfb357e6020cbe494fdc6cd2ca0 Mon Sep 17 00:00:00 2001 From: Jean-Michel Scherer Date: Mon, 4 Jul 2022 17:18:22 -0700 Subject: [PATCH 12/16] fix implementation of coupled damage-plasticity after commits from may/june 2020 (ATmodel and split operator) --- m_DefMech/m_MEF90_DefMechAssembly.F90 | 177 +++++++++++++++++--------- 1 file changed, 118 insertions(+), 59 deletions(-) diff --git a/m_DefMech/m_MEF90_DefMechAssembly.F90 b/m_DefMech/m_MEF90_DefMechAssembly.F90 index 4ab5a43cd0..2b553a786c 100644 --- a/m_DefMech/m_MEF90_DefMechAssembly.F90 +++ b/m_DefMech/m_MEF90_DefMechAssembly.F90 @@ -1212,6 +1212,8 @@ Subroutine MEF90DefMechStress(displacement,MEF90DefMechCtx,stress,ierr) PetscReal,Dimension(:),Pointer :: stressCellPtr PetscReal :: damageGauss,temperatureGauss PetscInt :: iDoF1,iDoF2,iGauss,numDofDisplacement,numDofDamage,numGauss + PetscReal,Dimension(:),Pointer :: plasticStrainLoc + Type(MEF90_MATS) :: plasticStrainCell Call PetscBagGetDataMEF90DefMechCtxGlobalOptions(MEF90DefMechCtx%GlobalOptionsBag,GlobalOptions,ierr);CHKERRQ(ierr) @@ -1280,6 +1282,13 @@ Subroutine MEF90DefMechStress(displacement,MEF90DefMechCtx,stress,ierr) Call SectionRealRestrictClosure(displacementSec,MEF90DefMechCtx%DMVect,cellID(cell),elemDisplacementType%numDof,displacementDof,ierr);CHKERRQ(ierr) Call SectionRealRestrictClosure(damageSec,MEF90DefMechCtx%DMScal,cellID(cell),elemDamageType%numDof,damageDof,ierr);CHKERRQ(ierr) + If (Associated(MEF90DefMechCtx%plasticStrain)) Then + Call SectionRealRestrict(plasticStrainSec,cellID(cell),plasticStrainLoc,ierr);CHKERRQ(ierr) + plasticStrainCell = plasticStrainLoc + Else + plasticStrainCell = 0.0_Kr + End If + !!! Loop over Gauss point, evaluate fields and assemble the cell contribution to the residual numDofDisplacement = size(elemDisplacement(cell)%BF,1) numDofDamage = size(elemDamage(cell)%BF,1) @@ -1306,7 +1315,8 @@ Subroutine MEF90DefMechStress(displacement,MEF90DefMechCtx,stress,ierr) inelasticStrainGauss = inelasticStrainGauss - (temperatureGauss * matpropSet%LinearThermalExpansion) End If ! Associated(temperatureDof) - Call Split%DEED(inelasticStrainGauss,matpropSet%HookesLaw,stressGaussPlus,stressGaussMinus) + ! Call Split%DEED(inelasticStrainGauss,matpropSet%HookesLaw,stressGaussPlus,stressGaussMinus) + Call Split%DEED(inelasticStrainGauss - plasticStrainCell,matpropSet%HookesLaw,stressGaussPlus,stressGaussMinus) If (cellIsElastic) Then stressCell = stressCell + (stressGaussPlus + stressGaussMinus) * elemDisplacement(cell)%Gauss_C(iGauss) Else @@ -1314,7 +1324,7 @@ Subroutine MEF90DefMechStress(displacement,MEF90DefMechCtx,stress,ierr) End If cellSize = cellSize + elemDisplacement(cell)%Gauss_C(iGauss) End Do ! iGauss - + stressCell = stressCell / cellSize stressCellPtr = stressCell Call SectionRealUpdate(stressSec,cellID(cell),stressCellPtr,INSERT_VALUES,ierr);CHKERRQ(ierr) @@ -1519,15 +1529,16 @@ Subroutine MEF90DefMechOperatorDamage(snesDamage,damage,residual,MEF90DefMechCt PetscErrorCode,Intent(OUT) :: ierr Type(SectionReal) :: displacementSec,residualSec - Type(SectionReal) :: plasticStrainSec,temperatureSec,damageSec,damagePreviousStepSec + Type(SectionReal) :: plasticStrainSec,temperatureSec,damageSec,damagePreviousStepSec,cumulatedDissipatedPlasticEnergySec Type(SectionReal) :: boundaryDamageSec,forceSec,pressureForceSec Type(SectionReal) :: CrackPressureSec PetscReal,Dimension(:),Pointer :: damagePtr,residualPtr PetscReal,Dimension(:),Pointer :: boundaryDamagePtr PetscReal,Dimension(:),Pointer :: displacementDof,damageDof,damagePreviousStepDof,temperatureDof,residualLoc - PetscReal,Dimension(:),Pointer :: plasticStrainLoc + PetscReal,Dimension(:),Pointer :: plasticStrainLoc,cumulatedDissipatedPlasticEnergyLoc PetscReal,Dimension(:),Pointer :: CrackPressureLoc Type(MEF90_MATS) :: plasticStrainCell + PetscReal :: cumulatedDissipatedPlasticEnergyCell Type(IS) :: VertexSetGlobalIS,CellSetGlobalIS,setIS,setISdof,bcIS PetscInt,Dimension(:),Pointer :: setID,cellID PetscInt,Dimension(:),Pointer :: setIdx,setdofIdx @@ -1548,10 +1559,11 @@ Subroutine MEF90DefMechOperatorDamage(snesDamage,damage,residual,MEF90DefMechCt Class(MEF90_DEFMECHSPLIT),Allocatable :: Split !Class(MEF90_DefMechDrivingForce_Type), Allocatable :: drivingForce PetscBool :: cellIsElastic - PetscReal :: damageGauss,damageDampingGauss,temperatureGauss,EEDGaussPlus,EEDGaussMinus,C1,C3,cellDampingCoefficient + PetscReal :: damageGauss,damageDampingGauss,temperatureGauss,EEDGaussPlus,EEDGaussMinus,C1,C3,cellDampingCoefficient,elasticEnergyDensityGauss Type(MEF90_VECT) :: gradDamageGauss,displacementGauss,C4 Type(MEF90_MATS) :: inelasticStrainGauss,C2 PetscInt :: iDoF1,iDoF2,iGauss,numDofDisplacement,numDofDamage,numGauss + PetscReal :: N Call SNESGetDM(snesDamage,MEF90DefMechCtx%DM,ierr);CHKERRQ(ierr) @@ -1599,6 +1611,13 @@ Subroutine MEF90DefMechOperatorDamage(snesDamage,damage,residual,MEF90DefMechCt Else PlasticStrainSec%v = 0 End If + + If (Associated(MEF90DefMechCtx%cumulatedDissipatedPlasticEnergy)) Then + Call SectionRealDuplicate(MEF90DefMechCtx%cellDMScalSec,cumulatedDissipatedPlasticEnergySec,ierr);CHKERRQ(ierr) + Call SectionRealToVec(cumulatedDissipatedPlasticEnergySec,MEF90DefMechCtx%CellDMScalScatter,SCATTER_REVERSE,MEF90DefMechCtx%cumulatedDissipatedPlasticEnergy,ierr);CHKERRQ(ierr) + Else + cumulatedDissipatedPlasticEnergySec%v = 0 + End If !!! get IS for cell sets Call DMMeshGetDimension(MEF90DefMechCtx%DM,dim,ierr);CHKERRQ(ierr) @@ -1625,6 +1644,7 @@ Subroutine MEF90DefMechOperatorDamage(snesDamage,damage,residual,MEF90DefMechCt C1 = matpropSet%fractureToughness / ATModel%cw * 0.25_Kr / matpropSet%internalLength C2 = matpropSet%fractureToughness / ATModel%cw * 0.5_Kr * matpropSet%internalLength * matpropSet%toughnessAnisotropyMatrix cellDampingCoefficient = globalOptions%dampingCoefficientDamage * MEF90DefMechCtx%timeStep + N = matpropSet%DuctileCouplingPower !!! Allocate storage for fields at dof and Gauss points !!! Leaving plasticstrain aside until I can remember how it is supposed to be dealt with @@ -1669,7 +1689,14 @@ Subroutine MEF90DefMechOperatorDamage(snesDamage,damage,residual,MEF90DefMechCt Else plasticStrainCell = 0.0_Kr End If - + + If (Associated(MEF90DefMechCtx%cumulatedDissipatedPlasticEnergy)) Then + Call SectionRealRestrict(cumulatedDissipatedPlasticEnergySec,cellID(cell),cumulatedDissipatedPlasticEnergyLoc,ierr);CHKERRQ(ierr) + cumulatedDissipatedPlasticEnergyCell = cumulatedDissipatedPlasticEnergyLoc(1) + Else + cumulatedDissipatedPlasticEnergyCell = 0.0_Kr + End If + !!! Loop over Gauss point, evaluate fields and assemble the cell contribution to the residual numDofDisplacement = size(elemDisplacement(cell)%BF,1) numDofDamage = size(elemDamage(cell)%BF,1) @@ -1681,7 +1708,7 @@ Subroutine MEF90DefMechOperatorDamage(snesDamage,damage,residual,MEF90DefMechCt damageGauss = 0.0_Kr gradDamageGauss = 0.0_Kr Do iDof1 = 1,numDofDamage - damageGauss = damageGauss + damageDof(iDoF1) * elemDamage(cell)%BF(iDoF1,iGauss) + damageGauss = damageGauss + damageDof(iDoF1) * elemDamage(cell)%BF(iDoF1,iGauss) gradDamageGauss = gradDamageGauss + damageDof(iDoF1) * elemDamage(cell)%Grad_BF(iDoF1,iGauss) End Do ! iDof1 @@ -1699,13 +1726,13 @@ Subroutine MEF90DefMechOperatorDamage(snesDamage,damage,residual,MEF90DefMechCt End Do ! iDof1 inelasticStrainGauss = inelasticStrainGauss - (temperatureGauss * matpropSet%LinearThermalExpansion) End If ! Associated(temperatureDof) -#if MEF90_DIM == 2 +! #if MEF90_DIM == 2 !!! We need something along these lines - !!! Adding terms in planestrain for plasticity with tr(p) = 0 - ! If (.NOT. matProp%HookesLaw%isPlaneStress) Then - ! stressGauss = stressGauss + stiffness * ( matProp%HookesLaw%lambda*trace(plasticStrainCell)*MEF90MatS2DIdentity ) - ! End If -#endif +! !!! Adding terms in planestrain for plasticity with tr(p) = 0 +! ! If (.NOT. matProp%HookesLaw%isPlaneStress) Then +! ! stressGauss = stressGauss + stiffness * ( matProp%HookesLaw%lambda*trace(plasticStrainCell)*MEF90MatS2DIdentity ) +! ! End If +! #endif If (.NOT. cellIsElastic) Then Call Split%EED(inelasticStrainGauss - plasticStrainCell,matpropSet%HookesLaw,EEDGaussPlus,EEDGaussMinus) @@ -1713,7 +1740,7 @@ Subroutine MEF90DefMechOperatorDamage(snesDamage,damage,residual,MEF90DefMechCt EEDGaussPlus = 0.0_Kr End If - C3 = ATModel%Da(damageGauss) * EEDGaussPlus + C1 * ATModel%Dw(damageGauss) + C3 = ATModel%Da(damageGauss) * EEDGaussPlus + C1 * ATModel%Dw(damageGauss) ! ? + N * cumulatedDissipatedPlasticEnergyCell * (1.0_Kr - damageGauss)**(N - 1.0_Kr) ) ? C4 = C2 * gradDamageGauss + CrackPressureCell * displacementGauss Do iDoF2 = 1,numDofDamage residualLoc(iDoF2) = residualLoc(iDoF2) + elemDamage(cell)%Gauss_C(iGauss) * ( & @@ -1738,28 +1765,37 @@ Subroutine MEF90DefMechOperatorDamage(snesDamage,damage,residual,MEF90DefMechCt ! Call localOperatorFunction2(residualLoc,damageDof,displacementDof,nullPtr,nullPtr,temperatureDof,plasticStrainCell,cumulatedDissipatedPlasticEnergyCell,matpropSet,elemDisplacement(cell),elemDamage(cell),CrackPressureCell) ! End If !!! Add plastic Dissipation here -! #if MEF90_DIM == 2 -! !!! Adding terms in planestrain for plasticity with tr(p) = 0 -! If (.NOT. matProp%HookesLaw%isPlaneStress) Then -! Stress_ZZ_planeStrain = ( matprop%HookesLaw%YoungsModulus - 2.0_Kr*matprop%HookesLaw%PoissonRatio*matprop%HookesLaw%mu )*trace(plasticStrainCell) + matprop%HookesLaw%lambda*trace(inelasticStrainGauss) -! elasticEnergyDensityGauss = elasticEnergyDensityGauss + ( Stress_ZZ_planeStrain + matprop%HookesLaw%lambda*trace(inelasticStrainGauss - plasticStrainCell) ) * trace(plasticStrainCell) -! End If -! #endif -! If (N==1.0_Kr) Then -! Do iDoF2 = 1,numDofDamage -! residualLoc(iDoF2) = residualLoc(iDoF2) + elemDamage%Gauss_C(iGauss) * ( & -! ( elasticEnergyDensityGauss * (damageGauss - 1.0_Kr) & -! - cumulatedDissipatedPlasticEnergyCell + C1 ) * elemDamage%BF(iDoF2,iGauss) + & -! (( C2 * matprop%toughnessAnisotropyMatrix * gradientDamageGauss + Displacement*CrackPressureCell ) .dotP. elemDamage%Grad_BF(iDoF2,iGauss)) ) -! End Do -! Else -! Do iDoF2 = 1,numDofDamage -! residualLoc(iDoF2) = residualLoc(iDoF2) + elemDamage%Gauss_C(iGauss) * ( & -! ( elasticEnergyDensityGauss * (damageGauss - 1.0_Kr) & -! - ( N * cumulatedDissipatedPlasticEnergyCell * (1.0_Kr - damageGauss)**(N - 1.0_Kr) ) + C1 ) * elemDamage%BF(iDoF2,iGauss) + & -! (( C2 * matprop%toughnessAnisotropyMatrix * gradientDamageGauss + Displacement*CrackPressureCell ) .dotP. elemDamage%Grad_BF(iDoF2,iGauss)) ) -! End Do -! EndIf +!#if MEF90_DIM == 2 +! !!! Adding terms in planestrain for plasticity with tr(p) = 0 +! If (.NOT. matProp%HookesLaw%isPlaneStress) Then +! Stress_ZZ_planeStrain = ( matprop%HookesLaw%YoungsModulus - 2.0_Kr*matprop%HookesLaw%PoissonRatio*matprop%HookesLaw%mu )*trace(plasticStrainCell) + matprop%HookesLaw%lambda*trace(inelasticStrainGauss) +! elasticEnergyDensityGauss = elasticEnergyDensityGauss + ( Stress_ZZ_planeStrain + matprop%HookesLaw%lambda*trace(inelasticStrainGauss - plasticStrainCell) ) * trace(plasticStrainCell) +! End If +!#endif + + If (N==1.0_Kr) Then + Do iDoF2 = 1,numDofDamage + residualLoc(iDoF2) = residualLoc(iDoF2) + elemDamage(cell)%Gauss_C(iGauss) * ( & + ( & + ! elasticEnergyDensityGauss * (damageGauss - 1.0_Kr) & + - cumulatedDissipatedPlasticEnergyCell & + !+ C1 + ) ) * elemDamage(cell)%BF(iDoF2,iGauss) + ! + & + ! (( C2 * matpropset%toughnessAnisotropyMatrix * gradDamageGauss + Displacement*CrackPressureCell ) .dotP. elemDamage%Grad_BF(iDoF2,iGauss)) ) + End Do + Else + Do iDoF2 = 1,numDofDamage + residualLoc(iDoF2) = residualLoc(iDoF2) + elemDamage(cell)%Gauss_C(iGauss) * ( & + ( & + ! elasticEnergyDensityGauss * (damageGauss - 1.0_Kr) & + - ( N * cumulatedDissipatedPlasticEnergyCell * (1.0_Kr - damageGauss)**(N - 1.0_Kr) ) & + !+ C1 + ) ) * elemDamage(cell)%BF(iDoF2,iGauss) + ! + & + ! (( C2 * matpropset%toughnessAnisotropyMatrix * gradDamageGauss + Displacement*CrackPressureCell ) .dotP. elemDamage%Grad_BF(iDoF2,iGauss)) ) + End Do + EndIf End Do ! iGauss Call SectionRealUpdateClosure(residualSec,MEF90DefMechCtx%DMScal,cellID(cell),residualLoc,ADD_VALUES,ierr);CHKERRQ(iErr) @@ -1889,10 +1925,11 @@ Subroutine MEF90DefMechBilinearFormDamage(snesDamage,damage,A,M,flg,MEF90DefMech Type(MEF90DefMechCtx_Type),Intent(IN) :: MEF90DefMechCtx PetscErrorCode,Intent(OUT) :: ierr - Type(SectionReal) :: displacementSec,damageSec,temperatureSec,plasticStrainSec + Type(SectionReal) :: displacementSec,damageSec,temperatureSec,plasticStrainSec,cumulatedDissipatedPlasticEnergySec PetscReal,Dimension(:),Pointer :: displacementDOF,damageDof,temperatureDof,nullPtr - PetscReal,Dimension(:),Pointer :: plasticStrainLoc + PetscReal,Dimension(:),Pointer :: plasticStrainLoc,cumulatedDissipatedPlasticEnergyLoc Type(MEF90_MATS) :: plasticStrainCell + PetscReal :: cumulatedDissipatedPlasticEnergyCell Type(IS) :: VertexSetGlobalIS,CellSetGlobalIS,setIS,setISdof PetscInt,dimension(:),Pointer :: setID PetscInt :: set,QuadratureOrder @@ -1913,10 +1950,11 @@ Subroutine MEF90DefMechBilinearFormDamage(snesDamage,damage,A,M,flg,MEF90DefMech Class(MEF90_DEFMECHSPLIT),Allocatable :: Split !Class(MEF90_DefMechDrivingForce_Type), Allocatable :: drivingForce PetscBool :: cellIsElastic - PetscReal :: damageGauss,damageDampingGauss,temperatureGauss,EEDGaussPlus,EEDGaussMinus,C1,C3,cellDampingCoefficient,BF1 + PetscReal :: damageGauss,damageDampingGauss,temperatureGauss,EEDGaussPlus,EEDGaussMinus,C1,C3,cellDampingCoefficient,BF1,elasticEnergyDensityGauss Type(MEF90_VECT) :: gradDamageGauss,displacementGauss,grad_BF1 Type(MEF90_MATS) :: inelasticStrainGauss,C2 PetscInt :: iDoF1,iDoF2,iGauss,numDofDisplacement,numDofDamage,numGauss + PetscReal :: N Call PetscBagGetDataMEF90DefMechCtxGlobalOptions(MEF90DefMechCtx%GlobalOptionsBag,GlobalOptions,ierr);CHKERRQ(ierr) Call MatZeroEntries(A,ierr);CHKERRQ(ierr) @@ -1943,7 +1981,14 @@ Subroutine MEF90DefMechBilinearFormDamage(snesDamage,damage,A,M,flg,MEF90DefMech Else PlasticStrainSec%v = 0 End If - + + If (Associated(MEF90DefMechCtx%cumulatedDissipatedPlasticEnergy)) Then + Call SectionRealDuplicate(MEF90DefMechCtx%CellDMScalSec,cumulatedDissipatedPlasticEnergySec,ierr);CHKERRQ(ierr) + Call SectionRealToVec(cumulatedDissipatedPlasticEnergySec,MEF90DefMechCtx%CellDMScalScatter,SCATTER_REVERSE,MEF90DefMechCtx%cumulatedDissipatedPlasticEnergy,ierr);CHKERRQ(ierr) + Else + cumulatedDissipatedPlasticEnergySec%v = 0 + End If + !!! get IS for cell sets Call DMMeshGetDimension(MEF90DefMechCtx%DM,dim,ierr);CHKERRQ(ierr) Call DMmeshGetLabelIdIS(MEF90DefMechCtx%DM,'Cell Sets',CellSetGlobalIS,ierr);CHKERRQ(ierr) @@ -1965,6 +2010,7 @@ Subroutine MEF90DefMechBilinearFormDamage(snesDamage,damage,A,M,flg,MEF90DefMech C1 = matpropSet%fractureToughness / ATModel%cw * 0.25_Kr / matpropSet%internalLength C2 = matpropSet%fractureToughness / ATModel%cw * 0.5_Kr * matpropSet%internalLength * matpropSet%toughnessAnisotropyMatrix cellDampingCoefficient = globalOptions%dampingCoefficientDamage * MEF90DefMechCtx%timeStep + N = matpropSet%DuctileCouplingPower !!! Allocate storage for fields at dof and Gauss points Allocate(Aloc(ElemDamageType%numDof,ElemDamageType%numDof)) @@ -1991,7 +2037,13 @@ Subroutine MEF90DefMechBilinearFormDamage(snesDamage,damage,A,M,flg,MEF90DefMech Else plasticStrainCell = 0.0_Kr End If - + If (Associated(MEF90DefMechCtx%cumulatedDissipatedPlasticEnergy)) Then + Call SectionRealRestrict(cumulatedDissipatedPlasticEnergySec,cellID(cell),cumulatedDissipatedPlasticEnergyLoc,ierr);CHKERRQ(ierr) + cumulatedDissipatedPlasticEnergyCell = cumulatedDissipatedPlasticEnergyLoc(1) + Else + cumulatedDissipatedPlasticEnergyCell = 0.0_Kr + End If + !!! Loop over Gauss point, evaluate fields and assemble the cell contribution to the bilinear form numDofDisplacement = size(elemDisplacement(cell)%BF,1) numDofDamage = size(elemDamage(cell)%BF,1) @@ -2032,6 +2084,8 @@ Subroutine MEF90DefMechBilinearFormDamage(snesDamage,damage,A,M,flg,MEF90DefMech Else EEDGaussPlus = 0.0_Kr End If + + ! elasticEnergyDensityGauss = (matpropset%HookesLaw * (inelasticStrainGauss - plasticStrainCell)) .DotP. (inelasticStrainGauss - plasticStrainCell) C3 = ATModel%D2a(damageGauss) * EEDGaussPlus + C1 * ATModel%D2w(damageGauss) Do iDoF1 = 1,numDofDamage @@ -2061,28 +2115,33 @@ Subroutine MEF90DefMechBilinearFormDamage(snesDamage,damage,A,M,flg,MEF90DefMech ! #endif ! !!! This is really twice the elastic energy density -! If ((N == 1.0_Kr) .OR. (cumulatedDissipatedPlasticEnergyCell == 0.0_Kr)) Then -! Do iDoF1 = 1,numDofDamage -! Do iDoF2 = 1,numDofDamage -! Aloc(iDoF2,iDoF1) = Aloc(iDoF2,iDoF1) + elemDamage%Gauss_C(iGauss) * ( & -! elasticEnergyDensityGauss * elemDamage%BF(iDoF1,iGauss) * elemDamage%BF(iDoF2,iGauss) + & -! C2 * ( (matprop%toughnessAnisotropyMatrix * elemDamage%Grad_BF(iDoF1,iGauss)) .dotP. elemDamage%Grad_BF(iDoF2,iGauss))) -! End Do -! End Do -! Else + If ((N == 1.0_Kr) .OR. (cumulatedDissipatedPlasticEnergyCell == 0.0_Kr)) Then + Do iDoF1 = 1,numDofDamage + BF1 = elemDamage(cell)%BF(iDoF1,iGauss) + Do iDoF2 = 1,numDofDamage + Aloc(iDoF2,iDoF1) = Aloc(iDoF2,iDoF1) !+ elemDamage(cell)%Gauss_C(iGauss) * ( & + !elasticEnergyDensityGauss * BF1 * elemDamage(cell)%BF(iDoF2,iGauss) ) + ! + & + ! C2 * ( (matpropset%toughnessAnisotropyMatrix * elemDamage%Grad_BF(iDoF1,iGauss)) .dotP. elemDamage%Grad_BF(iDoF2,iGauss))) + End Do + End Do + Else ! damageGauss = 0.0_Kr ! Do iDoF1 = 1,numDofDamage ! damageGauss = damageGauss + elemDamage%BF(iDoF1,iGauss) * xDof(iDoF1) ! End Do -! Do iDoF1 = 1,numDofDamage -! Do iDoF2 = 1,numDofDamage -! Aloc(iDoF2,iDoF1) = Aloc(iDoF2,iDoF1) + elemDamage%Gauss_C(iGauss) * ( & -! (elasticEnergyDensityGauss + & -! ( N * (N - 1.0_Kr) *( (1.0_Kr-damageGauss)**(N - 2.0_Kr) ) ) * cumulatedDissipatedPlasticEnergyCell ) * elemDamage%BF(iDoF1,iGauss) * elemDamage%BF(iDoF2,iGauss) + & -! C2 * ((matprop%toughnessAnisotropyMatrix * elemDamage%Grad_BF(iDoF1,iGauss)) .dotP. elemDamage%Grad_BF(iDoF2,iGauss))) -! End Do -! End Do -! EndIf + Do iDoF1 = 1,numDofDamage + BF1 = elemDamage(cell)%BF(iDoF1,iGauss) + Do iDoF2 = 1,numDofDamage + Aloc(iDoF2,iDoF1) = Aloc(iDoF2,iDoF1) + elemDamage(cell)%Gauss_C(iGauss) * ( & + ( & + ! elasticEnergyDensityGauss + & + ( N * (N - 1.0_Kr) *( (1.0_Kr-damageGauss)**(N - 2.0_Kr) ) ) * cumulatedDissipatedPlasticEnergyCell ) * BF1 * elemDamage(cell)%BF(iDoF2,iGauss) ) + ! + & + ! C2 * ((matpropset%toughnessAnisotropyMatrix * elemDamage%Grad_BF(iDoF1,iGauss)) .dotP. elemDamage%Grad_BF(iDoF2,iGauss))) + End Do + End Do + EndIf End Do ! iGauss Call DMmeshAssembleMatrix(A,MEF90DefMechCtx%DM,damageSec,cellID(cell),ALoc,ADD_VALUES,ierr);CHKERRQ(ierr) From ac463313ca8d2343d91c5a1bb1b59575dce97816 Mon Sep 17 00:00:00 2001 From: Jean-Michel Scherer Date: Tue, 5 Jul 2022 15:50:30 -0700 Subject: [PATCH 13/16] implement rate-independent crystal plasticity based on the Gambin-Arminjon regularization (L-inf norm regularized by L-p norm): see for instance Manik, Asadkandi and Holmedal (CMAME, 2022) --- MEF90/m_MEF90_Materials.F90 | 8 +++++ m_DefMech/m_MEF90_DefMechPlasticity.F90 | 46 +++++++++++++++---------- 2 files changed, 35 insertions(+), 19 deletions(-) diff --git a/MEF90/m_MEF90_Materials.F90 b/MEF90/m_MEF90_Materials.F90 index 72ce95b1f2..737d4fb0cf 100644 --- a/MEF90/m_MEF90_Materials.F90 +++ b/MEF90/m_MEF90_Materials.F90 @@ -79,6 +79,7 @@ Module m_MEF90_Materials_Types PetscReal :: ViscosityGamma0 ! viscosity reference slip rate PetscReal :: ViscosityN ! viscosity exponent PetscReal :: Viscositydt ! time step size + PetscReal :: m ! equivalent stress exponent for rate-independent crystal plasticity Character(len=MEF90_MXSTRLEN) :: Name End Type MEF90MatProp2D_Type @@ -124,6 +125,7 @@ Module m_MEF90_Materials_Types PetscReal :: ViscosityGamma0 ! viscosity reference slip rate PetscReal :: ViscosityN ! viscosity exponent PetscReal :: Viscositydt ! time step size + PetscReal :: m ! equivalent stress exponent for rate-independent crystal plasticity Character(len=MEF90_MXSTRLEN) :: Name End Type MEF90MatProp3D_Type @@ -196,6 +198,7 @@ Module m_MEF90_Materials_Types 1.0_Kr, & ! ViscosityGamma0 1.0_Kr, & ! ViscosityN 1.0_Kr, & ! Viscositydt + 1.0_Kr, & ! m "MEF90Mathium2D") Type(MEF90MatProp3D_Type),Parameter :: MEF90Mathium3D = MEF90MatProp3D_Type( & @@ -260,6 +263,7 @@ Module m_MEF90_Materials_Types 1.0_Kr, & ! ViscosityGamma0 1.0_Kr, & ! ViscosityN 1.0_Kr, & ! Viscositydt + 1.0_Kr, & ! m "MEF90Mathium3D") End Module m_MEF90_Materials_Types @@ -568,6 +572,8 @@ Subroutine PetscBagRegisterMEF90MatProp2D(bag,name,prefix,default,ierr) Call PetscBagRegisterReal(bag,matprop%ViscosityGamma0,default%ViscosityGamma0,'ViscosityGamma0','[s^(-1)] Reference plastic deformation rate',ierr);CHKERRQ(ierr) Call PetscBagRegisterReal(bag,matprop%ViscosityN,default%ViscosityN,'ViscosityN','[unit-less] Viscosity exponent',ierr);CHKERRQ(ierr) Call PetscBagRegisterReal(bag,matprop%Viscositydt,default%Viscositydt,'Viscositydt','[s] Viscosity time step size',ierr);CHKERRQ(ierr) + + Call PetscBagRegisterReal(bag,matprop%m,default%m,'m','[unit-less] Equivalent stress exponent for rate-independent crystal plasticity',ierr);CHKERRQ(ierr) !Call PetscBagSetFromOptions(bag,ierr) End Subroutine PetscBagRegisterMEF90MatProp2D @@ -659,6 +665,8 @@ Subroutine PetscBagRegisterMEF90MatProp3D(bag,name,prefix,default,ierr) Call PetscBagRegisterReal(bag,matprop%ViscosityGamma0,default%ViscosityGamma0,'ViscosityGamma0','[s^(-1)] Reference plastic deformation rate',ierr);CHKERRQ(ierr) Call PetscBagRegisterReal(bag,matprop%ViscosityN,default%ViscosityN,'ViscosityN','[unit-less] Viscosity exponent',ierr);CHKERRQ(ierr) Call PetscBagRegisterReal(bag,matprop%Viscositydt,default%Viscositydt,'Viscositydt','[s] Viscosity time step size',ierr);CHKERRQ(ierr) + + Call PetscBagRegisterReal(bag,matprop%m,default%m,'m','[unit-less] Equivalent stress exponent for rate-independent crystal plasticity',ierr);CHKERRQ(ierr) !Call PetscBagSetFromOptions(bag,ierr) End Subroutine PetscBagRegisterMEF90MatProp3D diff --git a/m_DefMech/m_MEF90_DefMechPlasticity.F90 b/m_DefMech/m_MEF90_DefMechPlasticity.F90 index 176371ff3e..a6dfb8df49 100644 --- a/m_DefMech/m_MEF90_DefMechPlasticity.F90 +++ b/m_DefMech/m_MEF90_DefMechPlasticity.F90 @@ -48,6 +48,7 @@ Module MEF90_APPEND(m_MEF90_DefMechPlasticity,MEF90_DIM)D real(Kind = Kr) :: ViscosityN real(Kind = Kr) :: Viscositydt real(Kind = Kr) :: viscouscumulatedDissipatedPlasticEnergyVariation + real(Kind = Kr) :: m end type MEF90DefMechPlasticityCtx contains @@ -870,8 +871,8 @@ subroutine FHG_CRYSTALBCC(x,f,h,g,myctx) bind(c) type(MEF90DefMechPlasticityCtx),pointer :: myctx_ptr type(MEF90_MATS) :: xMatS type(MEF90_MATS) :: InelasticStrain,PlasticStrainFlow,Stress - type(MatS3D) :: Strain3D,PlasticStrainFlow3D,Stress3D,Stress3DCrystal,TotalPlasticIncrement,TotalPlasticIncrementCrystal !, MatrixMu - type(MatS3D) :: PlasticStrain3D,PlasticStrain3DCrystal + type(MatS3D) :: Strain3D,PlasticStrainFlow3D,Stress3D,Stress3DCrystal,TotalPlasticIncrement,TotalPlasticIncrementCrystal + type(MatS3D) :: PlasticStrain3D,PlasticStrain3DCrystal,PlasticStrainFlow3DCrystal type(MatS3D),dimension(12) :: MatrixMu type(MAT3D) :: MatrixR real(Kind = Kr),dimension(12) :: ResolvedShearStress,PlasticSlipIncrement @@ -887,7 +888,7 @@ subroutine FHG_CRYSTALBCC(x,f,h,g,myctx) bind(c) & 0.,-1., 1., 1., 1., 0., 1., 0., 1., -1., 1., 0., 1., 0., 1., 0., 1., 1. /), (/3,12/)) !!! slip planes normals real(Kind = Kr),dimension(3,12) :: m_s = reshape((/ 1., 1., 1., 1., 1., 1., 1., 1., 1., 1.,-1., 1., 1.,-1., 1., 1.,-1., 1., & &-1., 1., 1., -1., 1., 1., -1., 1., 1., -1.,-1., 1., -1.,-1., 1., -1.,-1., 1. /), (/3,12/)) !!! slip directions - real(Kind = Kr) :: lambda,mu,E,nu + real(Kind = Kr) :: lambda,mu,E,nu,taueq !!! Casting x into a MEF90_MATS xMatS = x(1:SIZEOFMEF90_MATS) @@ -949,6 +950,7 @@ subroutine FHG_CRYSTALBCC(x,f,h,g,myctx) bind(c) #endif Stress3DCrystal = MatRaRt(Stress3D,myctx_ptr%RotationMatrix3D%fullTensor) + PlasticStrainFlow3DCrystal = MatRaRt(PlasticStrainFlow3D,myctx_ptr%RotationMatrix3D%fullTensor) normS = (2.0_Kr*SQRT(6.0_Kr)) Do s=1,12 @@ -968,8 +970,7 @@ subroutine FHG_CRYSTALBCC(x,f,h,g,myctx) bind(c) dt = myctx_ptr%Viscositydt active = 0 - - f(1) = 0.5_Kr * StiffnessA * Stress .DotP. (myctx_ptr%InelasticStrain-xMatS) + taueq = 0.0_Kr Do s = 1,12 ResolvedShearStress(s) = Stress3DCrystal .DotP. MatrixMu(s) CRSS = myctx_ptr%YieldTau0 @@ -985,26 +986,28 @@ subroutine FHG_CRYSTALBCC(x,f,h,g,myctx) bind(c) !myctx_ptr%plasticSlipsVariation(s) = PlasticSlipIncrement(s) myctx_ptr%viscouscumulatedDissipatedPlasticEnergyVariation = myctx_ptr%viscouscumulatedDissipatedPlasticEnergyVariation + ResolvedShearStress(s)*PlasticSlipIncrement(s) else - print *, "Rate-independent crystal plasticity is not implemented." - !PlasticSlipIncrement(s) = dt * myctx_ptr%eta * SIGN(1.0_Kr, ResolvedShearStress(s)) *& - ! & MAX( (ABS(StiffnessA*ResolvedShearStress(s)) - StiffnessB*CRSS) / myctx_ptr%YieldTau0 , 0. )**myctx_ptr%ViscosityN - !TotalPlasticIncrementCrystal = TotalPlasticIncrementCrystal + (PlasticSlipIncrement(s) * MatrixMu(s)) - !myctx_ptr%plasticSlipsVariation(s) = PlasticSlipIncrement(s) - !myctx_ptr%viscouscumulatedDissipatedPlasticEnergyVariation = myctx_ptr%viscouscumulatedDissipatedPlasticEnergyVariation + ResolvedShearStress(s)*PlasticSlipIncrement(s) + taueq = taueq + ABS( (StiffnessA * ResolvedShearStress(s)) / CRSS ) ** myctx_ptr%m + !myctx_ptr%PlasticSlipsVariation(s) = SIGN(1.0_KR, ResolvedShearStress(s)) * (Stress3DCrystal .DotP. PlasticStrainFlow3DCrystal) * (ABS( (StiffnessA * ResolvedShearStress(s)) / CRSS ) ** (myctx_ptr%m - 1.0_Kr)) / CRSS endif if ((ABS(StiffnessA*ResolvedShearStress(s)) - StiffnessB*CRSS)>0) then active = active + 1 end if End Do - f(1) = f(1) + StiffnessB * myctx_ptr%viscouscumulatedDissipatedPlasticEnergyVariation - TotalPlasticIncrement = MatRtaR(TotalPlasticIncrementCrystal,myctx_ptr%RotationMatrix3D%fullTensor) + if (myctx_ptr%isViscousPlasticity) then + f(1) = (0.5_Kr * StiffnessA * Stress .DotP. (myctx_ptr%InelasticStrain-xMatS)) + (StiffnessB * myctx_ptr%viscouscumulatedDissipatedPlasticEnergyVariation) + TotalPlasticIncrement = MatRtaR(TotalPlasticIncrementCrystal,myctx_ptr%RotationMatrix3D%fullTensor) #if MEF90_DIM==2 - h(1) = PlasticStrainFlow3D%XX - TotalPlasticIncrement%XX - h(2) = PlasticStrainFlow3D%XY - TotalPlasticIncrement%XY - h(3) = PlasticStrainFlow3D%YY - TotalPlasticIncrement%YY + h(1) = PlasticStrainFlow3D%XX - TotalPlasticIncrement%XX + h(2) = PlasticStrainFlow3D%XY - TotalPlasticIncrement%XY + h(3) = PlasticStrainFlow3D%YY - TotalPlasticIncrement%YY #elif MEF90_DIM==3 - h(1) = NORM(PlasticStrainFlow3D - TotalPlasticIncrement) + h(1) = NORM(PlasticStrainFlow3D - TotalPlasticIncrement) #endif + else + f(1) = ( (myctx_ptr%HookesLaw *(xMatS-myctx_ptr%PlasticStrainOld)) .DotP. (xMatS-myctx_ptr%PlasticStrainOld) ) + g(1) = ((taueq) ** (1.0_Kr/myctx_ptr%m)) - StiffnessB + h(1) = Trace(xMatS) + end if end subroutine FHG_CRYSTALBCC @@ -1220,14 +1223,18 @@ Subroutine MEF90DefMechPlasticStrainUpdate(MEF90DefMechCtx,plasticStrain,x,Plast if (matPropSet%HookesLaw%isPlaneStress) then write(*,*) "Plane stress CrystalBCC is not implemented" end if - snlp_m = 3 + if (matpropSet%isViscousPlasticity) then + snlp_m = 3 + else + snlp_m = 1 + end if #elif MEF90_DIM==3 snlp_m = 1 #endif if (matpropSet%isViscousPlasticity) then snlp_p = 0 else - snlp_p = 0 + snlp_p = 1 end if snlp_ctx = c_loc(PlasticityCtx) case default @@ -1268,6 +1275,7 @@ Subroutine MEF90DefMechPlasticStrainUpdate(MEF90DefMechCtx,plasticStrain,x,Plast PlasticityCtx%ViscosityGamma0 = matpropSet%ViscosityGamma0 PlasticityCtx%ViscosityN = matpropSet%ViscosityN PlasticityCtx%Viscositydt = matpropSet%Viscositydt + PlasticityCtx%m = matpropSet%m #if MEF90_DIM == 2 From 4e4754d8c98a5b21b04dce369e1f7092524c7184 Mon Sep 17 00:00:00 2001 From: Jean-Michel Scherer Date: Wed, 6 Jul 2022 10:00:32 -0700 Subject: [PATCH 14/16] fix memory leak --- m_DefMech/m_MEF90_DefMechAssembly.F90 | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) diff --git a/m_DefMech/m_MEF90_DefMechAssembly.F90 b/m_DefMech/m_MEF90_DefMechAssembly.F90 index 2b553a786c..2395287363 100644 --- a/m_DefMech/m_MEF90_DefMechAssembly.F90 +++ b/m_DefMech/m_MEF90_DefMechAssembly.F90 @@ -1807,6 +1807,11 @@ Subroutine MEF90DefMechOperatorDamage(snesDamage,damage,residual,MEF90DefMechCt If (Associated(MEF90DefMechCtx%CrackPressure)) Then Call SectionRealRestore(CrackPressureSec,cellID(cell),CrackPressureLoc,ierr);CHKERRQ(ierr) EndIf + + If (Associated(MEF90DefMechCtx%cumulatedDissipatedPlasticEnergy)) Then + Call SectionRealRestore(cumulatedDissipatedPlasticEnergySec,cellID(cell),cumulatedDissipatedPlasticEnergyLoc,ierr);CHKERRQ(ierr) + EndIf + End Do ! cell DeAllocate(displacementDof) DeAllocate(damageDof) @@ -1899,6 +1904,10 @@ Subroutine MEF90DefMechOperatorDamage(snesDamage,damage,residual,MEF90DefMechCt Call SectionRealDestroy(temperatureSec,ierr);CHKERRQ(ierr) End If + If (Associated(MEF90DefMechCtx%cumulatedDissipatedPlasticEnergy)) Then + Call SectionRealDestroy(cumulatedDissipatedPlasticEnergySec,ierr);CHKERRQ(ierr) + End If + Call SectionRealDestroy(damageSec,ierr);CHKERRQ(ierr) If (globalOptions%dampingCoefficientDamage * MEF90DefMechCtx%timeStep /= 0.0_Kr) Then Call SectionRealDestroy(damagePreviousStepSec,ierr);CHKERRQ(ierr) @@ -2149,6 +2158,10 @@ Subroutine MEF90DefMechBilinearFormDamage(snesDamage,damage,A,M,flg,MEF90DefMech If (Associated(MEF90DefMechCtx%plasticStrain)) Then Call SectionRealRestore(plasticStrainSec,cellID(cell),plasticStrainLoc,ierr);CHKERRQ(ierr) End If + + If (Associated(MEF90DefMechCtx%cumulatedDissipatedPlasticEnergy)) Then + Call SectionRealRestore(cumulatedDissipatedPlasticEnergySec,cellID(cell),cumulatedDissipatedPlasticEnergyLoc,ierr);CHKERRQ(ierr) + End If End Do ! cell DeAllocate(displacementDof) DeAllocate(damageDof) @@ -2216,6 +2229,10 @@ Subroutine MEF90DefMechBilinearFormDamage(snesDamage,damage,A,M,flg,MEF90DefMech If (Associated(MEF90DefMechCtx%plasticStrain)) Then Call SectionRealDestroy(plasticStrainSec,ierr);CHKERRQ(ierr) End If + + If (Associated(MEF90DefMechCtx%cumulatedDissipatedPlasticEnergy)) Then + Call SectionRealDestroy(cumulatedDissipatedPlasticEnergySec,ierr);CHKERRQ(ierr) + End If flg = SAME_NONZERO_PATTERN End Subroutine MEF90DefMechBilinearFormDamage From 7472672d8f55393417ea3b53d20a4b3e80fb9158 Mon Sep 17 00:00:00 2001 From: Jean-Michel Scherer Date: Wed, 6 Jul 2022 17:55:55 -0700 Subject: [PATCH 15/16] the trace of the 3D total plastic strain tensor must be zero --- m_DefMech/m_MEF90_DefMechPlasticity.F90 | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/m_DefMech/m_MEF90_DefMechPlasticity.F90 b/m_DefMech/m_MEF90_DefMechPlasticity.F90 index a6dfb8df49..121d87e46a 100644 --- a/m_DefMech/m_MEF90_DefMechPlasticity.F90 +++ b/m_DefMech/m_MEF90_DefMechPlasticity.F90 @@ -987,15 +987,18 @@ subroutine FHG_CRYSTALBCC(x,f,h,g,myctx) bind(c) myctx_ptr%viscouscumulatedDissipatedPlasticEnergyVariation = myctx_ptr%viscouscumulatedDissipatedPlasticEnergyVariation + ResolvedShearStress(s)*PlasticSlipIncrement(s) else taueq = taueq + ABS( (StiffnessA * ResolvedShearStress(s)) / CRSS ) ** myctx_ptr%m + PlasticSlipIncrement(s) = SIGN(1.0_KR, ResolvedShearStress(s)) * (Stress3DCrystal .DotP. PlasticStrainFlow3DCrystal) * (ABS( (StiffnessA * ResolvedShearStress(s)) / CRSS ) ** (myctx_ptr%m - 1.0_Kr)) / CRSS + TotalPlasticIncrementCrystal = TotalPlasticIncrementCrystal + (PlasticSlipIncrement(s) * MatrixMu(s)) + !myctx_ptr%plasticSlipsVariation(s) = PlasticSlipIncrement(s) !myctx_ptr%PlasticSlipsVariation(s) = SIGN(1.0_KR, ResolvedShearStress(s)) * (Stress3DCrystal .DotP. PlasticStrainFlow3DCrystal) * (ABS( (StiffnessA * ResolvedShearStress(s)) / CRSS ) ** (myctx_ptr%m - 1.0_Kr)) / CRSS endif if ((ABS(StiffnessA*ResolvedShearStress(s)) - StiffnessB*CRSS)>0) then active = active + 1 end if End Do + TotalPlasticIncrement = MatRtaR(TotalPlasticIncrementCrystal,myctx_ptr%RotationMatrix3D%fullTensor) if (myctx_ptr%isViscousPlasticity) then f(1) = (0.5_Kr * StiffnessA * Stress .DotP. (myctx_ptr%InelasticStrain-xMatS)) + (StiffnessB * myctx_ptr%viscouscumulatedDissipatedPlasticEnergyVariation) - TotalPlasticIncrement = MatRtaR(TotalPlasticIncrementCrystal,myctx_ptr%RotationMatrix3D%fullTensor) #if MEF90_DIM==2 h(1) = PlasticStrainFlow3D%XX - TotalPlasticIncrement%XX h(2) = PlasticStrainFlow3D%XY - TotalPlasticIncrement%XY @@ -1006,7 +1009,7 @@ subroutine FHG_CRYSTALBCC(x,f,h,g,myctx) bind(c) else f(1) = ( (myctx_ptr%HookesLaw *(xMatS-myctx_ptr%PlasticStrainOld)) .DotP. (xMatS-myctx_ptr%PlasticStrainOld) ) g(1) = ((taueq) ** (1.0_Kr/myctx_ptr%m)) - StiffnessB - h(1) = Trace(xMatS) + h(1) = Trace(TotalPlasticIncrementCrystal) end if end subroutine FHG_CRYSTALBCC From 988d3a191deeadd84ce62ab6f67d8b1038bd36ca Mon Sep 17 00:00:00 2001 From: Jean-Michel Scherer Date: Thu, 28 Jul 2022 15:23:50 -0700 Subject: [PATCH 16/16] compute dissipation with plastic strain only --- m_DefMech/m_MEF90_DefMechPlasticity.F90 | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/m_DefMech/m_MEF90_DefMechPlasticity.F90 b/m_DefMech/m_MEF90_DefMechPlasticity.F90 index 121d87e46a..c1be340b2f 100644 --- a/m_DefMech/m_MEF90_DefMechPlasticity.F90 +++ b/m_DefMech/m_MEF90_DefMechPlasticity.F90 @@ -1007,10 +1007,13 @@ subroutine FHG_CRYSTALBCC(x,f,h,g,myctx) bind(c) h(1) = NORM(PlasticStrainFlow3D - TotalPlasticIncrement) #endif else - f(1) = ( (myctx_ptr%HookesLaw *(xMatS-myctx_ptr%PlasticStrainOld)) .DotP. (xMatS-myctx_ptr%PlasticStrainOld) ) + f(1) = PlasticStrainFlow .DotP. PlasticStrainFlow !( (myctx_ptr%HookesLaw *(xMatS-myctx_ptr%PlasticStrainOld)) .DotP. (xMatS-myctx_ptr%PlasticStrainOld) ) g(1) = ((taueq) ** (1.0_Kr/myctx_ptr%m)) - StiffnessB h(1) = Trace(TotalPlasticIncrementCrystal) end if + !print *,"f(1) = ", f(1) + !print *,"h(1) = ", h(1) + !print *,"g(1) = ", g(1) end subroutine FHG_CRYSTALBCC @@ -1296,6 +1299,7 @@ Subroutine MEF90DefMechPlasticStrainUpdate(MEF90DefMechCtx,plasticStrain,x,Plast elemDisplacement,elemDisplacementType,elemScal,elemScalType,ierr) Allocate(damageloc(elemScalType%numDof)) Do cell = 1,size(cellID) + !print *,"cell = ", cell !! actualiser le ctx ( HookesLaw ,InelasticStrainSec, plasticStrainStrainSec, plasticStrainOldSec ) Call SectionRealRestrict(plasticStrainSec,cellID(cell),plasticStrainLoc,ierr);CHKERRQ(ierr) Call SectionRealRestrict(plasticStrainOldSec,cellID(cell),plasticStrainOldLoc,ierr);CHKERRQ(ierr)