clean and check implementation for offset method 0 and 1

NonLinearSolver/materialLaw/mlawNonLocalPorousCoupled.cpp

@@ -368,10 +368,15 @@ void mlawNonLocalPorousCoupledLaw::checkCoalescence(IPNonLocalPorosity* q1, cons
   // Coupling management
   if (_withCrackTransition){
     // If crack transition is used, 
-      // At the interface: coalescence is used as crack insertion criterion
+      // At the interface: coalescence is used, namely as crack insertion criterion
       // In the bulk, nothing is done (== always Gurson)
-    if (q1->getLocation() == IPVariable::INTERFACE_MINUS or q1->getLocation() == IPVariable::INTERFACE_PLUS)
-    {
+    if (q1->getLocation() == IPVariable::INTERFACE_MINUS or q1->getLocation() == IPVariable::INTERFACE_PLUS){
+      // At the interface:
+      // Check if transition has already occurs or not
+      if (!q0Thom->getCoalescenceOnsetFlag()){
+        // If Coalescence has not yet occured at least once
+        // Check if the normal traction force is sufficient to obtain coalescence in normal direction of the interface
+      
         // Get normal traction on the interface from tau stress
         double sigmax = 0.;
         static SVector3 Norm;
@@ -384,47 +389,40 @@ void mlawNonLocalPorousCoupledLaw::checkCoalescence(IPNonLocalPorosity* q1, cons
             }
           }
           // update maximal stress -- will be removed soon
-        q1Thom->getRefToMaximalTractionStress()= q0Thom->getMaximalTractionStress();
-        if (sigmax > q1Thom->getMaximalTractionStress()){
-          q1Thom->getRefToMaximalTractionStress() = sigmax;
-        }
+          //q1Thom->getRefToMaximalTractionStress()= q0Thom->getMaximalTractionStress();
+          //if (sigmax > q1Thom->getMaximalTractionStress()){
+          //  q1Thom->getRefToMaximalTractionStress() = sigmax;
+          //}
         }
         else{
           Msg::Fatal("mlawNonLocalPorousCoupledLaw::updatePlasticState: zero normal");
         }
         
         // Check if the normal traction force is sufficient to obtain coalescence in normal direction of the interface
-      if (!q0Thom->getCoalescenceOnsetFlag())
-      {
-        // Check Thomason criterion:
+        // and manage offset
+          // Why is the offset needed?
+          // We want to ensure a continuity with the Thomason yield surface at coalescence onset
+          // The offset correspond to a factor to obtain fT = 0 with the actual stress state
+            // -actual stress state:  yieldThomason = sthingOfStress - Cft*R/R0
+            // -correction to satidfy fT=0 with actual stress:  yieldThomason + Cft*R/R0 - CftOffset*Cft*R/R0  = 0
+            // -solving for offset gives: CftOffset = (yieldThomason + Cft*R/R0)/(Cft*R/R0) = 1+yieldThomason*R0/(Cft*R)
+        
+        // Compute Thomason yield surface
         double yieldThomason = _mlawCoales->yieldFunction(kcorEq,pcor,R,yieldfV,q0,q1,T); 
-        /* check this value */
-        if(_CfTOffsetMethod == 2)
-        {
+        
+        // Offset management
+        if(_CfTOffsetMethod == 2){
+          // Method 2: Offset is always computed to ensure contuinity if coalescence is forced (offset can be <1)
           q1Thom->getRefToCrackOffsetOnCft() = 1. + yieldThomason*R0/(Cft*R);
         }
+        
         if ((sigmax - Cft*R > _tol) and (yieldThomason > _tol)){
           // If coalescence occurs, update Cft offset value and others variables
-          //
-          // why offset? because:
-          // we want to follow the continuity of the Thomason yeild surface after satisfying coalescence
-          // the value of offset can be larger than 1
-          // this value is found by resolving the equation
-          // before coalescence yieldThomason = sthingOfStress - Cft*R/R0
-          // at coalescence  yieldThomason + Cft*R/R0 - CftOffset*Cft*R/R0  = 0
-          // leading to offset = (yieldThomason + Cft*R/R0)/(Cft*R/R0) = 1+yieldThomason*R0/(Cft*R)
-
-          // Compute Cft Offset
-          if (_CfTOffsetMethod == 1){
-            q1Thom->getRefToCrackOffsetOnCft() = 1. + yieldThomason*R0/(Cft*R);
-            Msg::Info("coalescence occurs for the first time: Offset =%e",q1Thom->getRefToCrackOffsetOnCft());
-          }
-          else if (_CfTOffsetMethod == 0){
-            q1Thom->getRefToCrackOffsetOnCft() = 1.;
-          }
           
+          // Update onset coalescence variables
+            // Activate coalescence onset
           q1Thom->getRefToCoalescenceOnsetFlag() = true;
-          // Set onset variables to their actual value
+            // Set onset variables to the actual values (yield porosity and geometrical parameters)
           q1Thom->getRefToPorosityAtCoalescenceOnset() = q1->getYieldPorosity();
           q1Thom->getRefToLigamentRatioAtCoalescenceOnset() = q1Thom->getLigamentRatio();
           q1Thom->getRefToAspectRatioAtCoalescenceOnset() = q1Thom->getAspectRatio();
@@ -432,42 +430,54 @@ void mlawNonLocalPorousCoupledLaw::checkCoalescence(IPNonLocalPorosity* q1, cons
           
           q1Thom->getRefToCoalescenceActiveFlag() = true;
           q1Thom->getRefToAccelerateRate() = 1.;
+          
+          // Offset management
+          if (_CfTOffsetMethod == 1){
+            // Method 1: offset computed at coalescence onset only
+            q1Thom->getRefToCrackOffsetOnCft() = 1. + yieldThomason*R0/(Cft*R);
+            Msg::Info("Coalescence occurs for the first time: Offset =%e",q1Thom->getRefToCrackOffsetOnCft());
+          }
+          else if (_CfTOffsetMethod == 0){
+            // Method 0: Offset is not used: always equal to 1.
+            q1Thom->getRefToCrackOffsetOnCft() = 1.;
+          }
+
+
         }
         else{
           // Coalescence never occurs and is not occuring at the next step
           q1Thom->getRefToCoalescenceOnsetFlag() = false;
-          // Reset onset variables to their default value
-            // Porosity
+          // Reset onset variables to their default value (yield porosity and geometrical parameters)
           q1Thom->getRefToPorosityAtCoalescenceOnset() = 0.;
-            // Geometrical parameters
           q1Thom->getRefToLigamentRatioAtCoalescenceOnset() = 0.;
           q1Thom->getRefToAspectRatioAtCoalescenceOnset() = 0.;
           q1Thom->getRefToShapeFactorAtCoalescenceOnset() = 0.;
 
           q1Thom->getRefToCoalescenceActiveFlag() = false;
-
           q1Thom->getRefToAccelerateRate() = 1.;
-          q1Thom->getRefToCrackOffsetOnCft() = 1.; // offset equal to 1 by default
+          
+          // Offset management
+          if (_CfTOffsetMethod == 1 or _CfTOffsetMethod == 0){
+            q1Thom->getRefToCrackOffsetOnCft() = 1.;
+          }
         }
       }
       else{
         // Coalescence has occured at least once
         q1Thom->getRefToCoalescenceOnsetFlag() = true;
-        // Update onset value
-          // Porosity
+        // Update onset value (yield porosity and geometrical parameters)
         q1Thom->getRefToPorosityAtCoalescenceOnset() = q0Thom->getPorosityAtCoalescenceOnset();
-          // Geometrical parameters
         q1Thom->getRefToLigamentRatioAtCoalescenceOnset() = q0Thom->getLigamentRatioAtCoalescenceOnset();
         q1Thom->getRefToAspectRatioAtCoalescenceOnset() = q0Thom->getAspectRatioAtCoalescenceOnset();
         q1Thom->getRefToShapeFactorAtCoalescenceOnset() = q0Thom->getShapeFactorAtCoalescenceOnset();
-          // Cft onset
+        
+        // Cft onset - Method 0-1-2: keep it always constant once the coalescence occurs
         q1Thom->getRefToCrackOffsetOnCft() = q0Thom->getCrackOffsetOnCft();
 
 
         // Determine which mode should be used for the next time step.
-        if (q0Thom->getCoalescenceActiveFlag())
-        {
-          // If Thomason is/was used for this time step : check if Gurson is better
+        if (q0Thom->getCoalescenceActiveFlag()){
+          // If Thomason is/was used for this time step : check if Gurson is more restrictive
           double yieldGurson = _mlawGrowth->yieldFunction(kcorEq,pcor,R,yieldfV,q0,q1,T);
           if(yieldGurson > _tol){
             Msg::Info("Change yield surface to Gurson fG =%e",yieldGurson);
@@ -480,9 +490,8 @@ void mlawNonLocalPorousCoupledLaw::checkCoalescence(IPNonLocalPorosity* q1, cons
           }
 
         }
-        else
-        {
-          // If Gurson was/is used for this time step : check if Thomason is better
+        else{
+          // If Gurson was/is used for this time step : check if Thomason is more restrictive
           double yieldThomason = _mlawCoales->yieldFunction(kcorEq,pcor,R,yieldfV,q0,q1,T);
           if(yieldThomason > _tol ){
             Msg::Info("Change yield surface to Thomason fT = %e",yieldThomason);