diff --git a/NonLinearSolver/materialLaw/mlawSMP.cpp b/NonLinearSolver/materialLaw/mlawSMP.cpp index a58b2f6fd916378232faff6452ab6cda3dc588e7..b3d5dbfa5acf1a945e77b97e7f27b162bdb3533d 100644 --- a/NonLinearSolver/materialLaw/mlawSMP.cpp +++ b/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.; - } */ - - double _htg =1.e-9;//pow(0.001,_r1); + } +#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.); + } + 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 + + }