correct plastic dissipation

NonLinearSolver/materialLaw/mlawAnisotropicPlasticity.cpp

@@ -752,8 +752,8 @@ bool mlawAnisotropicPlasticity::predictorCorector(const STensor3& F0,         //
   q->getRefToPlasticPower() = 0.;
     // plastic power (Wp1- Wp0)/dt
   if (eps-eps0 > 0.){
-    q->getRefToPlasticEnergy() += DeltaHatP*R;
-    q->getRefToPlasticPower() = DeltaHatP/this->getTimeStep();
+    q->getRefToPlasticEnergy() += detF*DeltaHatP*R;
+    q->getRefToPlasticPower() = detF*DeltaHatP*R/this->getTimeStep();
   }
 
   // irreversible energy
@@ -800,8 +800,8 @@ bool mlawAnisotropicPlasticity::predictorCorector(const STensor3& F0,         //
       }
 
       STensor3& dplasticpowerDF = q->getRefToDPlasticPowerDF();
-      dplasticpowerDF = dpdF;
-      dplasticpowerDF *= (R+H*(eps-eps0))/this->getTimeStep();
+      dplasticpowerDF = q->getConstRefToDPlasticEnergyDF();
+      dplasticpowerDF *= (1./this->getTimeStep());
     }
     else
     {
@@ -878,7 +878,7 @@ bool mlawAnisotropicPlasticity::predictorCorector(const STensor3& F0,         //
         
         // plastic energy DplasticEnegy = R*DeltaHatP
         STensor3& DplasticEnergyDF = q->getRefToDPlasticEnergyDF();
-        STensorOperation::scale(DDeltaHatPDEepr,H*DeltaHatP+R,DplasticEnergyDF);
+        STensorOperation::scale(DDeltaHatPDEepr,detF*(H*DeltaHatP+R),DplasticEnergyDF);
       }
     
       static STensor43 EprToF;
@@ -907,7 +907,18 @@ bool mlawAnisotropicPlasticity::predictorCorector(const STensor3& F0,         //
       DplEnergyDEpr = q->getConstRefToDPlasticEnergyDF();
       STensor3& DplEnergyDF = q->getRefToDPlasticEnergyDF();
       STensorOperation::multSTensor3STensor43(DplEnergyDEpr,EprToF,DplEnergyDF);
-
+      if (_stressFormulation == CORO_CAUCHY)
+      {
+        static STensor3 invF;
+        STensorOperation::inverseSTensor3(F, invF);
+        for (int i=0; i<3; i++)
+        {
+          for (int j=0; j< 3; j++)
+          {
+            DplEnergyDF(i,j) += detF*invF(j,i)*DeltaHatP*R;
+          }
+        }
+      }
       if (this->getMacroSolver()->withPathFollowing()){
         // irreversible energy
         if ((this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::PLASTIC_ENERGY) or

NonLinearSolver/materialLaw/mlawJ2linear.cpp

@@ -582,8 +582,8 @@ void mlawJ2linear::predictorCorector(const STensor3& F0,         // initial defo
   q->getRefToPlasticEnergy() = q0->plasticEnergy();
     // plastic power (Wp1- Wp0)/dt
   if (Deps > 0.){
-    q->getRefToPlasticEnergy() += Deps*Sy;
-    q->getRefToPlasticPower() = Deps*Sy/this->getTimeStep();
+    q->getRefToPlasticEnergy() += Jformule*Deps*Sy;
+    q->getRefToPlasticPower() = Jformule*Deps*Sy/this->getTimeStep();
   }
   else
   {
@@ -631,8 +631,8 @@ void mlawJ2linear::predictorCorector(const STensor3& F0,         // initial defo
       }
 
       STensor3& dplasticpowerDF = q->getRefToDPlasticPowerDF();
-      dplasticpowerDF = dpdF;
-      dplasticpowerDF *= (Sy+H*Deps)/this->getTimeStep();
+      dplasticpowerDF = q->getConstRefToDIrreversibleEnergyDF();
+      dplasticpowerDF *= (1./this->getTimeStep());
     }
     else{
       static STensor43 DcorKirDEepr;
@@ -882,24 +882,31 @@ void mlawJ2linear::predictorCorector(const STensor3& F0,         // initial defo
       }
 
 
-      // plastic power
-      STensor3& dplasticpowerDF = q->getRefToDPlasticPowerDF();
-      if (Deps > 0.){
-        dplasticpowerDF = dpdF;
-        dplasticpowerDF *= (Sy+H*Deps)/this->getTimeStep();
-      }
-      else{
-        STensorOperation::zero(dplasticpowerDF);
-      }
       
       STensor3& DplasticEnergyDF = q->getRefToDPlasticEnergyDF();
       if (Deps > 0.){
         DplasticEnergyDF = dpdF;
-        DplasticEnergyDF *= (Sy+H*Deps);
+        DplasticEnergyDF *= Jformule*(Sy+H*Deps);
+        if (_stressFormulation == CORO_CAUCHY)
+        {
+          static STensor3 Finv;
+          STensorOperation::inverseSTensor3(F, Finv);
+          for (int i=0; i<3; i++)
+          {
+            for (int j=0; j<3; j++)
+            {
+              DplasticEnergyDF(i,j) += Jformule*Finv(j,i)*Deps*Sy;
+            }
+          }
+        }
       }
       else{
         STensorOperation::zero(DplasticEnergyDF);
       }
+      // plastic power
+      STensor3& dplasticpowerDF = q->getRefToDPlasticPowerDF();
+      dplasticpowerDF = DplasticEnergyDF;
+      dplasticpowerDF *= (1./this->getTimeStep());
 
       // irreversible energy
       STensor3& DirrEnegDF = q->getRefToDIrreversibleEnergyDF();
@@ -919,8 +926,7 @@ void mlawJ2linear::predictorCorector(const STensor3& F0,         // initial defo
       }
       else if ((this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::PLASTIC_ENERGY) or
                (this->getMacroSolver()->getPathFollowingLocalIncrementType() == pathFollowingManager::DISSIPATION_ENERGY)) {
-        DirrEnegDF = dpdF;
-        DirrEnegDF*= (Sy+H*Deps);
+        DirrEnegDF = q->getConstRefToDPlasticEnergyDF();
       }
       else{
         STensorOperation::zero(DirrEnegDF);