change haredening law

NonLinearSolver/materialLaw/mlawSMP.cpp

@@ -1791,32 +1791,45 @@ void mlawSMP:: computeSaphi(const double depsilon1 ,const double Sa0,const doubl
 			    const double dh,double &phiastar0,double &Sastar0,double &Sa1,double &Sastar1,double &phia1,double &phiastar1, 
 			    double &dSa1depsilon, double &dSa1dTg, double &dSa1dT)const
 {
-  /* if(T<=Tg  && depsilon1>0.)
+	
+  double _htg =1.e-9;//pow(0.001,_r1);
+#if 0
+  if(T<=Tg  && depsilon1>0.)
   {
-    phiastar1=_Z1*pow(1.-T/Tg,_r1)*pow(depsilon1/(dh*_epsilonr),_s1);
+    phiastar1=_Z1*(pow(1.-T/Tg,_r1)+_htg)*pow(depsilon1/(dh*_epsilonr),_s1);
   }
-  else if (T>Tg  && depsilon1>0.)//check in the serevern ....old and not making difference
+  else if (T>Tg  && depsilon1>0.)//check in the serever
   { 
-    phiastar1=_Z1*pow(T/Tg-1,_r1)*pow(depsilon1/(dh*_epsilonr),_s1);
+    phiastar1=_Z1*_htg*pow(depsilon1/(dh*_epsilonr),_s1);
   }
   else
   {
     phiastar1=0.;
-  }	*/
+  }	
+#endif
+#if 1
+  double tmp=(pow(1.+depsilon1/(dh*_epsilonr),_s1)-1.+tanh(depsilon1/(dh*_epsilonr)));
+  double Dtmp=1./(dh*_epsilonr)*(_s1*pow(1.+depsilon1/(dh*_epsilonr),_s1-1.)+1./cosh(depsilon1/(dh*_epsilonr))/cosh(depsilon1/(dh*_epsilonr)));
+  bool oldlaw=true;
+  if(oldlaw)
+  {
+    tmp=pow(depsilon1/(dh*_epsilonr),_s1);
+    Dtmp=_s1/(dh*_epsilonr)*pow(depsilon1/(dh*_epsilonr),_s1-1.);
+  }
 
-  double _htg =1.e-9;//pow(0.001,_r1);
   if(T<=Tg  && depsilon1>0.)
   {
-    phiastar1=_Z1*(pow(1.-T/Tg,_r1)+_htg)*pow(depsilon1/(dh*_epsilonr),_s1);
+    phiastar1=_Z1*(pow(1.-T/Tg,_r1)+_htg)*tmp;
   }
   else if (T>Tg  && depsilon1>0.)//check in the serever
   {
-    phiastar1=_Z1*_htg*pow(depsilon1/(dh*_epsilonr),_s1);
+    phiastar1=_Z1*_htg*tmp;
   }
   else
   {
     phiastar1=0.;
   }
+#endif
   
   phia1=(phia0 +_g1*_be1* phiastar1*depsilon1+_g1*(1.-_be1)*(phiastar0-phia0)*depsilon1)/(1.+_be1*_g1*depsilon1);	
   Sastar1=_b1*(phiastar1-phia1);
@@ -1831,6 +1844,7 @@ void mlawSMP:: computeSaphi(const double depsilon1 ,const double Sa0,const doubl
   dSa1depsilon+=(pow(_be1,2.)*_ha1*(1.+_be1*_g1*depsilon1)+pow(_be1,2.)*_g1*(1.+_be1*_ha1*depsilon1))*_ha1*_b1*phia0*depsilon1/(pow((1.+_be1*_ha1*depsilon1),2.)*pow((1.+_be1*_g1*depsilon1),2.));
   dSa1depsilon+=(pow(_be1,2.)*_ha1*(1.+_be1*_g1*depsilon1)+pow(_be1,2.)*_g1*(1.+_be1*_ha1*depsilon1))*_ha1*_b1*(_g1*(1.-_be1)*(phiastar0-phia0)*depsilon1)*depsilon1/(pow((1.+_be1*_ha1*depsilon1),2.)*pow((1.+_be1*_g1*depsilon1),2.));
      
+#if 0
   if(T<=Tg  && depsilon1>0.)
   {
     dSa1depsilon+=(_be1*_ha1*_b1*_Z1*(pow(1.-T/Tg,_r1)+_htg)*pow(depsilon1/(dh*_epsilonr),_s1))/((1.+_be1*_ha1*depsilon1)*(1.+_be1*_g1*depsilon1));
@@ -1854,6 +1868,34 @@ void mlawSMP:: computeSaphi(const double depsilon1 ,const double Sa0,const doubl
   {
 	  
   }
+#endif
+#if 1
+  if(T<=Tg  && depsilon1>0.)
+  {
+    dSa1depsilon+=(_be1*_ha1*_b1*_Z1*(pow(1.-T/Tg,_r1)+_htg)*tmp)/((1.+_be1*_ha1*depsilon1)*(1.+_be1*_g1*depsilon1));
+    dSa1depsilon-=(pow(_be1,2.)*_ha1*(1.+_be1*_g1*depsilon1)+pow(_be1,2.)*_g1*(1.+_be1*_ha1*depsilon1))*_ha1*_b1*_Z1*(pow(1.-T/Tg,_r1)+_htg)*tmp*depsilon1/(pow((1.+_be1*_ha1*depsilon1),2.)*pow((1.+_be1*_g1*depsilon1),2.));
+	    
+    dSa1depsilon+=_be1*_ha1*_b1*_Z1*(pow(1.-T/Tg,_r1)+_htg)*Dtmp*depsilon1/((1.+_be1*_ha1*depsilon1)*(1.+_be1*_g1*depsilon1));
+     
+    dSa1dTg=(_r1*_be1*_ha1*_b1*_Z1*T/(Tg*Tg)*(pow(1.-T/Tg,_r1-1.)+_htg)*tmp*depsilon1)/((1.+_be1*_ha1*depsilon1)*(1.+_be1*_g1*depsilon1));	
+    dSa1dT=-_r1*_be1*_ha1*_b1*_Z1/Tg*(pow(1.-T/Tg,_r1-1.)+_htg)*tmp*depsilon1/((1.+_be1*_ha1*depsilon1)*(1.+_be1*_g1*depsilon1));
+  }
+  else if (T>Tg  && depsilon1>0.)
+  { 
+    dSa1depsilon+=(_be1*_ha1*_b1*_Z1*(_htg)*tmp)/((1.+_be1*_ha1*depsilon1)*(1.+_be1*_g1*depsilon1));
+    dSa1depsilon-=(pow(_be1,2.)*_ha1*(1.+_be1*_g1*depsilon1)+pow(_be1,2.)*_g1*(1.+_be1*_ha1*depsilon1))*_ha1*_b1*_Z1*(_htg)*tmp*depsilon1/(pow((1.+_be1*_ha1*depsilon1),2.)*pow((1.+_be1*_g1*depsilon1),2.));	    
+    dSa1depsilon+=_be1*_ha1*_b1*_Z1*(_htg)*Dtmp*depsilon1/((1.+_be1*_ha1*depsilon1)*(1.+_be1*_g1*depsilon1)); 
+    dSa1dTg=0.;	
+    dSa1dT =0.;
+	  
+  }
+  else
+  {
+	  
+  }
+#endif
+
+
 }