add energy integral

NonLinearSolver/materialLaw/mlawHyperelastic.cpp

@@ -98,7 +98,7 @@ void mlawHyperViscoElastic::setViscoElasticData(const std::string filename){
   fclose(file);
 };
 
-void mlawHyperViscoElastic::evaluatePhiPCorrection(double tr, const STensor3 &dev, double &A_v, double &dA_vdE, double &B_d, STensor3 &dB_vddev) const
+void mlawHyperViscoElastic::evaluatePhiPCorrection(double tr, const STensor3 &dev, double &A_v, double &dA_vdE, double &intA, double &B_d, STensor3 &dB_vddev, double &intB) const
 {
   int method =1;
   if(method ==0)
@@ -114,18 +114,37 @@ void mlawHyperViscoElastic::evaluatePhiPCorrection(double tr, const STensor3 &de
  
     dB_vddev=dev;
     dB_vddev*=2.*getPiDevCorrection()*getDevCorrection()*pow(exp(getThetaDevCorrection()*dev.dotprod())-1.,getPiDevCorrection()-1.)* exp(getThetaDevCorrection()*dev.dotprod())*getThetaDevCorrection();
+    Msg::Error("mlawHyperViscoElastic::evaluatePhiPCorrection: need to evaluate the volume energy");
   }
   else
   {
     A_v = getVolumeCorrection()*(tanh(getXiVolumeCorrection()/3.*tr*tr-getZetaVolumeCorrection())+tanh(getZetaVolumeCorrection()));  
     dA_vdE = getVolumeCorrection()*getXiVolumeCorrection()*2./3.*(1.-tanh(getXiVolumeCorrection()/3.*tr*tr-getZetaVolumeCorrection())*tanh(getXiVolumeCorrection()/3.*tr*tr-getZetaVolumeCorrection()));
   
+    intA=0.;
+    if (getXiVolumeCorrection()!=0.)
+    {
+      double u=tr*tr;
+      intA=0.5*getVolumeCorrection()*(3./getXiVolumeCorrection()*(log(cosh(getXiVolumeCorrection()/3.*u-getZetaVolumeCorrection())) -log(cosh(-getZetaVolumeCorrection())) ) +tanh(getZetaVolumeCorrection())*u ) ;
+    
+    }
   
     B_d = getDevCorrection()*(tanh(getThetaDevCorrection()*dev.dotprod()-getPiDevCorrection())+tanh(getPiDevCorrection()));
     STensorOperation::zero(dB_vddev);
  
     dB_vddev=dev;
     dB_vddev*=2.*getPiDevCorrection()*getDevCorrection()*(1.-tanh(getThetaDevCorrection()*dev.dotprod()-getPiDevCorrection())*tanh(getThetaDevCorrection()*dev.dotprod()-getPiDevCorrection()) );
+    intB=0.;
+    if(getThetaDevCorrection()!=0)
+    {
+      double u=dev.dotprod();
+      intB=0.5*getDevCorrection()*(1./getThetaDevCorrection()*(log(cosh(getThetaDevCorrection()*u-getPiDevCorrection())) -log(cosh(-getPiDevCorrection())) ) +tanh(getPiDevCorrection())*u ) ;
+    
+    }
+    
+    
+    
+    
     
   }
 }
@@ -140,10 +159,10 @@ double mlawHyperViscoElastic::deformationEnergy(const IPHyperViscoElastic& q) co
     static STensor3 devEe;
     STensorOperation::decomposeDevTr(q._Ee,devEe,trEe);
     
-    double Av, dAv,B_d;
+    double Av, dAv,B_d, intA,intB;
     static STensor3 dB_d;
-    evaluatePhiPCorrection(trEe, devEe, Av, dAv, B_d, dB_d);
-    Psy = getUpdatedBulkModulus(&q)*(1.+Av)*(0.5*trEe*trEe)+getUpdatedShearModulus(&q)*(1.+B_d)*STensorOperation::doubledot(devEe,devEe); //this is not correct we should do an integral
+    evaluatePhiPCorrection(trEe, devEe, Av, dAv, intA, B_d, dB_d, intB);
+    Psy = getUpdatedBulkModulus(&q)*(0.5*trEe*trEe+intA)+2.*getUpdatedShearModulus(&q)*(0.5*STensorOperation::doubledot(devEe,devEe)+intB); //this is not correct we should do an integral
 																// 
 																//    
     for (int i=0; i<_N; i++){
@@ -417,7 +436,7 @@ void mlawHyperViscoElastic::updateViscoElasticFlow(const IPHyperViscoElastic *q0
 };
 void mlawHyperViscoElastic::viscoElasticPredictor(const STensor3& Ee, const STensor3& Ee0,
           const IPHyperViscoElastic *q0, IPHyperViscoElastic *q1,
-          double& Ke, double& Ge,     double& A_v,  double& B_d, double& dApr, STensor3& dB_vddev) const{
+          double& Ke, double& Ge,     double& A_v,  double& B_d, double& dApr, STensor3& dB_vddev, double &intA, double &intB) const{
   if ((_Ki.size() > 0) or (_Gi.size() > 0)){
     static STensor3 DE, devDE;
     static double trDE;
@@ -455,7 +474,7 @@ void mlawHyperViscoElastic::viscoElasticPredictor(const STensor3& Ee, const STen
       double p, trEe;
       STensorOperation::decomposeDevTr(Ee,devE,trEe);
 
-      evaluatePhiPCorrection(trEe, devE, A_v, dApr, B_d, dB_vddev);
+      evaluatePhiPCorrection(trEe, devE, A_v, dApr, intA, B_d, dB_vddev, intB);
 
       devK = devE;
       devK *= (2.*getUpdatedShearModulus(q1)*(1.+B_d));  // deviatoric part
@@ -547,7 +566,7 @@ void mlawHyperViscoElastic::viscoElasticPredictor(const STensor3& Ee, const STen
     STensorOperation::decomposeDevTr(Ee,devE,trEe);
 
     STensorOperation::zero(dB_vddev);
-    evaluatePhiPCorrection(trEe, devE, A_v, dApr, B_d, dB_vddev);
+    evaluatePhiPCorrection(trEe, devE, A_v, dApr, intA, B_d, dB_vddev, intB);
             
     devK = devE;
     devK *= (2.*getUpdatedShearModulus(q1)*(1.+B_d));  // deviatoric part
@@ -602,10 +621,10 @@ void mlawHyperViscoElastic::predictorCorrector_ViscoElastic(const STensor3& F0,
   STensorOperation::decomposeDevTr(E,devE,trE);
   
   double Ke, Ge;
-  double A_v_pr=0., B_d_pr = 0., dApr=0.;
+  double A_v_pr=0., B_d_pr = 0., dApr=0., intApr=0., intBpr=0.;
   static STensor3 dB_devpr;
   STensorOperation::zero(dB_devpr);
-  viscoElasticPredictor(E,q0->_Ee,q0,q1,Ke,Ge,A_v_pr, B_d_pr, dApr, dB_devpr);
+  viscoElasticPredictor(E,q0->_Ee,q0,q1,Ke,Ge,A_v_pr, B_d_pr, dApr, dB_devpr,intApr,intBpr);
 
   const STensor3& corKir = q1->getConstRefToCorotationalKirchhoffStress();
   static STensor3 secondPK;
@@ -894,10 +913,10 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
   
   
   
-  double A_v_pr, B_d_pr, dApr;
+  double A_v_pr=0., B_d_pr=0., dApr=0.,intApr=0.,intBpr=0.;
   static STensor3 dB_devpr;
   STensorOperation::zero(dB_devpr);
-  viscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,A_v_pr, B_d_pr, dApr, dB_devpr);
+  viscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,A_v_pr, B_d_pr, dApr, dB_devpr, intApr, intBpr);
   
   STensorOperation::decomposeDevTr(Ee,devEepr,trEepr); 
   //std::cout << "trEepr = "<<trEepr <<std::endl;
@@ -980,7 +999,7 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
   
   double A_v=0.;
   double B_d = 0.;
-  double dAdtrEe=0.;
+  double dAdtrEe=0., intA=0., intB=0.;
   static STensor3 dBddevEe;
   STensorOperation::zero(dBddevEe);
   
@@ -1078,7 +1097,7 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
   
         while (1){
 
-        	evaluatePhiPCorrection(trEe, devEe, A_v, dAdtrEe, B_d, dBddevEe);
+        	evaluatePhiPCorrection(trEe, devEe, A_v, dAdtrEe, intA, B_d, dBddevEe, intB);
         	J=getUpdatedBulkModulus(q1)*A_v*trEe-ptilde_Vol;
         	dJdptilde_Vol = getUpdatedBulkModulus(q1)*( dAdtrEe*dtrEedptilde*trEe+A_v*dtrEedptilde)-1.;
        
@@ -1121,7 +1140,7 @@ void mlawPowerYieldHyper::predictorCorrector_nonAssociatedFlow(const STensor3& F
           
         while (1){
   
-          	evaluatePhiPCorrection(trEe, devEe, A_v, dAdtrEe, B_d, dBddevEe);
+          	evaluatePhiPCorrection(trEe, devEe, A_v, dAdtrEe, intA, B_d, dBddevEe, intB);
 
                 for (int i=0; i<3; i++){
         	        for (int j=0; j<3; j++){
@@ -1735,10 +1754,10 @@ void mlawPowerYieldHyper::predictorCorrector_associatedFlow(const STensor3& F, c
   
   
   
-  double A_v_pr=0., B_d_pr = 0., dApr=0.;
+  double A_v_pr=0., B_d_pr = 0., dApr=0., intApr=0., intBpr=0.;
   static STensor3 dB_devpr;
   STensorOperation::zero(dB_devpr);
-  viscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,A_v_pr, B_d_pr, dApr, dB_devpr);
+  viscoElasticPredictor(Ee,q0->_Ee,q0,q1,Ke,Ge,A_v_pr, B_d_pr, dApr, dB_devpr, intApr, intBpr);
 
   static STensor3 devKpr; // dev corotational kirchoff stress predictor
   static double ppr; // pressure predictor

NonLinearSolver/materialLaw/mlawHyperelastic.h

@@ -59,14 +59,14 @@ class mlawHyperViscoElastic : public materialLaw{
     void updateViscoElasticFlow(const IPHyperViscoElastic *q0, IPHyperViscoElastic *q1, double& Ke, double& Ge) const;
     void viscoElasticPredictor(const STensor3& Ee, const STensor3& Ee0,
               const IPHyperViscoElastic *q0, IPHyperViscoElastic *q1,
-              double& Ke, double& Ge,     double& A_v,  double& B_d, double& dAprd, STensor3& dB_vddev) const;
+              double& Ke, double& Ge,     double& A_v,  double& B_d, double& dAprd, STensor3& dB_vddev, double &intA, double &intB) const;
 
     void isotropicHookTensor(const double K, const double G, STensor43& L) const;
 
     void predictorCorrector_ViscoElastic(const STensor3& F0, const STensor3& F, STensor3&P, const IPHyperViscoElastic *q0_, IPHyperViscoElastic *q1,
                                   const bool stiff, STensor43& Tangent) const;
 
-    void evaluatePhiPCorrection(double trEe, const STensor3 &devEpr,  double &A_v, double &dAprdEepr, double &B_d, STensor3 &dB_vddevEe) const;
+    void evaluatePhiPCorrection(double trEe, const STensor3 &devEpr,  double &A_v, double &dAprdEepr, double &intA, double &B_d, STensor3 &dB_vddevEe, double &intB) const;
     //void evaluateA_dA(double trEe, double &A_v, double &dAprdEepr) const;
     
   #endif // SWIG