Merge pull request #12393 from KratosMultiphysics/poro/feature-postpr…

…ocess-tools [Poro] Add new post-processing variables
KratosMultiphysics · May 21, 2024 · 8956b64 · 8956b64
2 parents 8419608 + c23f58f
commit 8956b64
Show file tree

Hide file tree

Showing 12 changed files with 221 additions and 29 deletions.
diff --git a/...custom_constitutive/interface_element_laws/elastoplastic_mohr_coulomb_cohesive_2D_law.cpp b/...custom_constitutive/interface_element_laws/elastoplastic_mohr_coulomb_cohesive_2D_law.cpp
@@ -45,6 +45,7 @@ void ElastoPlasticMohrCoulombCohesive2DLaw::InitializeMaterial( const Properties
     mPlasticStrainVector[1]    = 0.0;
     mOldPlasticStrainVector[0] = 0.0;
     mOldPlasticStrainVector[1] = 0.0;
+    mSlipTendency = 0.0;
 }
 
 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

diff --git a/...custom_constitutive/interface_element_laws/elastoplastic_mohr_coulomb_cohesive_3D_law.cpp b/...custom_constitutive/interface_element_laws/elastoplastic_mohr_coulomb_cohesive_3D_law.cpp
@@ -97,6 +97,7 @@ void ElastoPlasticMohrCoulombCohesive3DLaw::InitializeMaterial( const Properties
     mOldPlasticStrainVector[0] = 0.0;
     mOldPlasticStrainVector[1] = 0.0;
     mOldPlasticStrainVector[2] = 0.0;
+    mSlipTendency = 0.0;
 }
 
 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
@@ -181,6 +182,18 @@ void ElastoPlasticMohrCoulombCohesive3DLaw::FinalizeMaterialResponseCauchy (Para
 
 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
 
+double& ElastoPlasticMohrCoulombCohesive3DLaw::GetValue( const Variable<double>& rThisVariable, double& rValue )
+{
+    if (rThisVariable == SLIP_TENDENCY) 
+    {
+        rValue = mSlipTendency;
+    }
+
+    return rValue;
+}
+
+//----------------------------------------------------------------------------------------
+
 Vector& ElastoPlasticMohrCoulombCohesive3DLaw::GetValue( const Variable<Vector>& rThisVariable, Vector& rValue )
 {
     rValue = mPlasticStrainVector;
@@ -269,7 +282,10 @@ double ElastoPlasticMohrCoulombCohesive3DLaw::GetShearResultantStressVector(Vect
 // This method computes the traction vector and the plastic multiplier based on a Backward-Euler integration scheme.
 // The implementation does not consider any softening or hardening rule.
 
-void ElastoPlasticMohrCoulombCohesive3DLaw::ReturnMapping(Vector& rStressVector, Matrix& rConstitutiveMatrix, Vector& TrialStressVector, Matrix& ElasticConstitutiveMatrix, ConstitutiveLawVariables& rVariables, ElastoPlasticConstitutiveLawVariables& rEPlasticVariables, Parameters& rValues)
+void ElastoPlasticMohrCoulombCohesive3DLaw::ReturnMapping(Vector& rStressVector, Matrix& rConstitutiveMatrix, 
+                                                          Vector& TrialStressVector, Matrix& ElasticConstitutiveMatrix, 
+                                                          ConstitutiveLawVariables& rVariables, ElastoPlasticConstitutiveLawVariables& rEPlasticVariables, 
+                                                          Parameters& rValues)
 {
     // Get the size of the strain vector
     const Vector& StrainVector = rValues.GetStrainVector();
@@ -298,7 +314,10 @@ void ElastoPlasticMohrCoulombCohesive3DLaw::ReturnMapping(Vector& rStressVector,
     rStressVector = TrialStressVector;             
 
     // Get the shear resultant
-    double ts = this->GetShearResultantStressVector(rStressVector);                
+    double ts = this->GetShearResultantStressVector(rStressVector);
+
+    // Get the normal component of the stress vector
+    double tn = rStressVector[VoigtSize-1];              
 
     // Compute the normal to the plastic potential surface (np) and its derivative wrt to the stress vector
     this->DerivativesPlasticPotentialSurface(rStressVector, rVariables, rEPlasticVariables, rValues);
@@ -327,6 +346,14 @@ void ElastoPlasticMohrCoulombCohesive3DLaw::ReturnMapping(Vector& rStressVector,
             rConstitutiveMatrix(VoigtSize-1, VoigtSize-1) = 0.0;
         }  
 
+        // Compute the slip tendency
+        if (rEPlasticVariables.YieldFunction_MC == 0.0) {
+            mSlipTendency = 1.0;        
+        }
+        else {
+            mSlipTendency = abs(ts) / abs(c - tn*tanPhi);
+        }
+
     } else if ((rEPlasticVariables.YieldFunction_MC > 0.0) && (rEPlasticVariables.YieldFunction_CutOff < 0.0)){ // ------ Return to the Mohr-Coulomb surface
 
         // Result from the product between n^T * Tel * np
@@ -352,6 +379,9 @@ void ElastoPlasticMohrCoulombCohesive3DLaw::ReturnMapping(Vector& rStressVector,
             noalias(rConstitutiveMatrix) = ElasticConstitutiveMatrix - outer_prod(Tel_np,Tel_n) / n_Tel_np;
         }
 
+        // Compute the slip tendency
+        mSlipTendency = 1.0;
+
     }else{ // ------------------------------------------------------------------------------------------- Return to the point which both surfaces intersect
 
         // Compute the plastic multipliers
@@ -371,6 +401,9 @@ void ElastoPlasticMohrCoulombCohesive3DLaw::ReturnMapping(Vector& rStressVector,
         if(Options.Is(ConstitutiveLaw::COMPUTE_CONSTITUTIVE_TENSOR)){
             this->ConstitutiveMatrixInstersectionYieldSurfaces(rStressVector,rConstitutiveMatrix, rVariables);
         }
+
+        // Compute the slip tendency
+        mSlipTendency = 1.0;
     }
 
 }

diff --git a/...custom_constitutive/interface_element_laws/elastoplastic_mohr_coulomb_cohesive_3D_law.hpp b/...custom_constitutive/interface_element_laws/elastoplastic_mohr_coulomb_cohesive_3D_law.hpp
@@ -73,6 +73,8 @@ class KRATOS_API(POROMECHANICS_APPLICATION) ElastoPlasticMohrCoulombCohesive3DLa
     void FinalizeMaterialResponseCauchy(Parameters & rValues) override;
 
 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
+    double& GetValue( const Variable<double>& rThisVariable, double& rValue ) override;
+
     Vector& GetValue( const Variable<Vector>& rThisVariable, Vector& rValue ) override;
 
     void SetValue( const Variable<Vector>& rThisVariable,
@@ -108,11 +110,13 @@ class KRATOS_API(POROMECHANICS_APPLICATION) ElastoPlasticMohrCoulombCohesive3DLa
         // Vector normal to the plastic potential surface
         Vector np_MC;
         Vector np_TC;
+
     };
 
     // Member Variables
-    Vector mPlasticStrainVector;    
+    Vector mPlasticStrainVector;
     Vector mOldPlasticStrainVector;
+    double mSlipTendency;
 
 //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------