Herezh_dev/comportement/thermique/Taux_crista/Hoffman2.cc

// This file is part of the Herezh++ application.
//
// The finite element software Herezh++ is dedicated to the field
// of mechanics for large transformations of solid structures.
// It is developed by Gérard Rio (APP: IDDN.FR.010.0106078.000.R.P.2006.035.20600)
// INSTITUT DE RECHERCHE DUPUY DE LÔME (IRDL) <https://www.irdl.fr/>.
//
// Herezh++ is distributed under GPL 3 license ou ultérieure.
//
// Copyright (C) 1997-2022 Université Bretagne Sud (France)
// AUTHOR : Gérard Rio
// E-MAIL  : gerardrio56@free.fr
//
// This program is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License,
// or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty
// of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// See the GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program. If not, see <https://www.gnu.org/licenses/>.
//
// For more information, please consult: <https://herezh.irdl.fr/>.

#include "Hoffman2.h"
#include <cmath>

#include "ConstMath.h"
#include "ConstPhysico.h"

//==================== cas de la class de sauvegarde SaveResul ===================

// constructeur  par défaut à ne pas utiliser
Hoffman2::SaveCrista_Hoffman2::SaveCrista_Hoffman2() :
	I_Kcinetique(0.),I_Kcinetique_t(0.)
  { }; 

// le constructeur courant
Hoffman2::SaveCrista_Hoffman2::SaveCrista_Hoffman2
                        (const double & cristal_t,const double & cristal):
	I_Kcinetique(cristal) ,I_Kcinetique_t(cristal_t) 
	 {}; 

// constructeur de copie 
Hoffman2::SaveCrista_Hoffman2::SaveCrista_Hoffman2
                                         (const Hoffman2::SaveCrista_Hoffman2& sav ):
	I_Kcinetique(sav.I_Kcinetique) ,I_Kcinetique_t(sav.I_Kcinetique_t)
	 { }; 
	   
     //============= lecture écriture dans base info ==========
     
// cas donne le niveau de la récupération
// = 1 : on récupère tout
// = 2 : on récupère uniquement les données variables (supposées comme telles)
void Hoffman2::SaveCrista_Hoffman2::Lecture_base_info 
                                                 (ifstream& ent,const int )
  {  // ici toutes les données sont toujours a priori variables
     // ou en tout cas pour les méthodes appelées, elles sont gérées par le paramètre: cas
     string toto; ent >> toto;
     #ifdef MISE_AU_POINT	 	 
      if (toto != "S_R_Hoffman")
          { cout << "\n erreur en lecture du conteneur pour la loi thermo physique de Tait"
                 << " on attendait S_R_Hoffman et on a lue " << toto
                 << " \n Hoffman2::SaveCrista_Hoffman2::Lecture_base_info(..";
            Sortie(1);     
          }
     #endif
//     ent >> I_Kcinetique_t >> I_Kcinetique;
     ent >> I_Kcinetique_t ; I_Kcinetique = I_Kcinetique_t;
    }; 

// cas donne le niveau de sauvegarde
// = 1 : on sauvegarde tout
// = 2 : on sauvegarde uniquement les données variables (supposées comme telles)
void Hoffman2::SaveCrista_Hoffman2::Ecriture_base_info
                                                 (ofstream& sort,const int )
   { // ici toutes les données sont toujours a priori variables
     // ou en tout cas pour les méthodes appelées, elles sont gérées par le paramètre: cas
//     sort << "\n S_R_Hoffman  " << I_Kcinetique_t << " " << I_Kcinetique << " ";
     sort << "\n S_R_Hoffman  " << I_Kcinetique_t << " " ;
	 }; 
	        
// affichage des infos
void Hoffman2::SaveCrista_Hoffman2::Affiche()	
   { cout << "\n S_R_Hoffman  I_Kcinetique_t= " << I_Kcinetique_t << " I_Kcinetique= " << I_Kcinetique << " ";
	 }; 
// idem sur un ofstream
void Hoffman2::SaveCrista_Hoffman2::Affiche(ofstream& sort)			        
   { sort << "\n S_R_Hoffman  I_Kcinetique_t= " << I_Kcinetique_t << " I_Kcinetique= " << I_Kcinetique << " ";
	 }; 	        

// Surcharge de l'operateur =  
CristaliniteAbstraite::SaveCrista& Hoffman2::SaveCrista_Hoffman2::operator= 
          (const CristaliniteAbstraite::SaveCrista& a)
 {  // récup de la bonne classe
    Hoffman2::SaveCrista_Hoffman2& b =  *((Hoffman2::SaveCrista_Hoffman2*) &a);
    I_Kcinetique = b.I_Kcinetique;
    I_Kcinetique_t = b.I_Kcinetique_t;
    return *this; 
  };

//==================== fin du cas de la class de sauvegarde SaveResul ============

Hoffman2::Hoffman2():
  n(ConstMath::Pi),invt012(ConstMath::Pi)
  ,Ustar(ConstMath::Pi),Tinf(ConstMath::Pi)
  ,Kg(ConstMath::Pi),Tm(ConstMath::Pi),alphaP(ConstMath::Pi)
  ,betaP(ConstMath::Pi),X_inf(ConstMath::Pi)
  ,algo_integ(new Algo_Integ1D(3)) // par défaut on utilise 3 points de Gauss
  ,Te_t(),Te_tdt(),Pr_t(),Pr_tdt() // variables de travail
    { };
    
// constructeur de copie
Hoffman2::Hoffman2(const Hoffman2& co ):
  n(co.n),invt012(co.invt012)
  ,Ustar(co.Ustar),Tinf(co.Tinf)
  ,Kg(co.Kg),Tm(co.Tm),alphaP(co.alphaP)
  ,betaP(co.betaP),X_inf(co.X_inf)
  ,algo_integ(new Algo_Integ1D(*(co.algo_integ)))
  ,Te_t(),Te_tdt(),Pr_t(),Pr_tdt() // variables de travail
    {};
    
// Lecture des donnees de la classe sur fichier
 void Hoffman2::LectureDonneesLoiCrista(UtilLecture * entreePrinc,LesCourbes1D& ,LesFonctions_nD& lesFonctionsnD)
{   // ---- lecture  du tag d'entête  -----
    string nom;
    *(entreePrinc->entree) >> nom; // lecture du tag 
    if(nom != "coefficientsLoiHoffman_")
         { cout << "\n erreur en lecture du tag du debut des parametres "
                << " on attendait la chaine : coefficientsLoiHoffman_ et on a lue: " << nom;
           cout << "\n Hoffman2::LectureDonneesLoiCrista(  "
                << "(UtilLecture * entreePrinc) " << endl ;
           entreePrinc->MessageBuffer("** erreur1 de lecture des parametres de la loi Hoffman2 **");
           throw (UtilLecture::ErrNouvelleDonnee(-1));
           Sortie(1);          
          };
    entreePrinc->NouvelleDonnee(); // prepa du flot de lecture
    
    // --- lecture de l'ensemble des paramètres  ----
    MotCle motCle; // ref aux mots cle
    while (strstr(entreePrinc->tablcar,"fin_coefficientsLoiHoffman_")==0) 
       {
        // si on a  un mot clé global dans la ligne courante c-a-d dans tablcar --> erreur
        if ( motCle.SimotCle(entreePrinc->tablcar))
         { cout << "\n erreur de lecture des parametre on n'a pas trouve le mot cle "
                << " fin_coefficientsLoiHoffman_ et par contre la ligne courante contient un mot cle global  ";
           entreePrinc->MessageBuffer("** erreur2  de lecture des parametres de la loi Hoffman2 **");
           throw (UtilLecture::ErrNouvelleDonnee(-1));
           Sortie(1);     
         };
         
        // lecture d'un mot clé
        *(entreePrinc->entree) >> nom;
        
        if ((entreePrinc->entree)->rdstate() == 0) 
          {} // lecture normale
        #ifdef ENLINUX 
        else  if ((entreePrinc->entree)->fail())
           // on a atteind la fin de la ligne et on appelle un nouvel enregistrement 
          {   entreePrinc->NouvelleDonnee();  // lecture d'un nouvelle enregistrement
              *(entreePrinc->entree) >>nom;
           }
        #else
        else  if ((entreePrinc->entree)->eof())
          // la lecture est bonne mais on a atteind la fin de la ligne
          { if(nom != "fin_coefficientsLoiHoffman_")
              {entreePrinc->NouvelleDonnee(); *(entreePrinc->entree) >> nom;};
          } 
        #endif  
        else // cas d'une erreur de lecture
         { cout << "\n erreur de lecture inconnue  ";
           entreePrinc->MessageBuffer("** erreur3  de lecture des parametres de la loi Hoffman2 **");
           throw (UtilLecture::ErrNouvelleDonnee(-1));
           Sortie(1);     
         };
        // dans le cas où le fait d'avoir lue une nouvelle ligne fait que l'on a fini la lecture on stop 
        if (nom == "fin_coefficientsLoiHoffman_") 
          {break;}
        // sinon c'est ok et  on passe en revue tous les paramètres
        else if(nom == "n=")            {*(entreePrinc->entree) >> n;}
        else if (nom == "invt012=")     {*(entreePrinc->entree) >> invt012;}
        else if (nom == "Ustar=")  {*(entreePrinc->entree) >> Ustar;}
        else if (nom == "Tinf=")   {*(entreePrinc->entree) >> Tinf;}
        else if (nom == "Kg=")     {*(entreePrinc->entree) >> Kg;}
        else if (nom == "Tm=")     {*(entreePrinc->entree) >> Tm;}
        else if (nom == "alphaP=") {*(entreePrinc->entree) >> alphaP;}
        else if (nom == "betaP=")  {*(entreePrinc->entree) >> betaP;}
        else if (nom == "X_inf=")  {*(entreePrinc->entree) >> X_inf;}
        else
         { cout << "\n erreur: lecture d'un identificateur inconnue, on a lue:  " << nom ;
           entreePrinc->MessageBuffer("** erreur4  de lecture des parametres de la loi Hoffman2 **");
           throw (UtilLecture::ErrNouvelleDonnee(-1));
           Sortie(1);     
         };
        
       }; //-- fin du while
     
    entreePrinc->NouvelleDonnee(); // prepa du flot de lecture
   };

// affichage de la loi
void Hoffman2::Affiche() const
 { cout << "\n# ----- loi de Hoffman2 permettant le calcul de K(T,P)  -------------------------------"
        << "\n     n= "      << n      << " invt012= "    << invt012    << " "
        << "\n     Ustar= "  << Ustar  << " Tinf= "  << Tinf  << " Kg= "  << Kg  << " Tm= "  << Tm << " "
        << "\n     alphaP= " << alphaP << " betaP= " << betaP << " X_inf= " << X_inf << " " ;
	}; 
    
// affichage et definition interactive des paramètres 
void Hoffman2::Info_commande_LoisCrista(UtilLecture& entreePrinc)
 {  ofstream & sort = *(entreePrinc.Commande_pointInfo()); // pour simplifier
 
    n= 3;             invt012= 1.528e8;       
    Ustar= 6284.;     Tinf= 215.;       Kg= 3.81e5;   Tm= 445.; 
    alphaP= 0; betaP= 0; X_inf= 0.443;     
    
    sort << "\n# ----- loi de Hoffman2 permettant le calcul de K(T,P)  ------------------------------"
         << "\n# |la loi necessite 11 coefficients  n,invt012,Ustar,Tinf,Kg,Tm,alphaP,betaP,X_inf   |"
         << "\n# |   le calcul de K(T,P) s'effectue alors suivant la formule:                       |"
         << "\n# |   K(T,P) = (log(2))**(1/n) * invt012 * exp(- Ustar/(R*(Tcorr - Tinf)))           |"
         << "\n# |         * exp(- Kg/(Tcorr * (Tm - Tcorr) * f))                                   |"   
         << "\n# |   avec les relations:                                                            |"
         << "\n# |       f=2*Tcorr/(Tcorr + Tm) ;                                                   |"
         << "\n# |       DTcorr_Pres = P*(alphaP + betaP*P);                                        |"
         << "\n# |       Tcorr       = T - DTcorr_Pres + 273.;                                      |"
         << "\n# |   et:  P la pression relative, T la temperature celsius                          |"
         << "\n# -----------------------------------------------------------------------------------";
         
    sort << "\n#  definition des coefficients, ils peuvent être donnés dans un ordre quelconque et "
         << "\n#  sur une ou plusieurs lignes "
         << "\n\n  coefficientsLoiHoffman_ "
         << "\n     n= "      << n      << " invt012= "    << invt012    << " "
         << "\n     Ustar= "  << Ustar  << " Tinf= "  << Tinf  << " Kg= "  << Kg  << " Tm= "  << Tm << " "
         << "\n     alphaP= " << alphaP << " betaP= " << betaP << " X_inf= " << X_inf << " " ;

    sort << "\n     fin_coefficientsLoiHoffman_       \n "; 
    sort << "\n#   la derniere ligne doit contenir uniquement le mot cle:     fin_coefficientsLoiHoffman_ "
         << endl;     
  };  		  	  

// calcul de la fonction K(T,P)
double  Hoffman2::fct_KT(const double& P, const double& T ) const
	{ // calcul des variables intermédiaires
	  double DTcorr_Pres = P*(alphaP + betaP*P);
	  double Tcorr       = T - DTcorr_Pres + 273.;
	  double f           = 2. * Tcorr/(Tcorr + Tm);
	  // calcul de KT
	  double KT = 0.;
	  if ((Tcorr > Tinf) && (Tcorr < Tm))
	     KT = pow((log(2.)),(1./n)) * invt012 * exp(- Ustar/( ConstPhysico::R *(Tcorr - Tinf))) 
                  * exp(- Kg/(Tcorr * (Tm - Tcorr) * f)); 
      // retour
      return KT;              
	};
	
// calcul du taux de cristalinité
// P_t, P : pression à t et, t+deltat (= actuelle)
// T_t, T : température à t et, t+deltat (=actuelle)
double  Hoffman2::Cristalinite(const double& P_t, const double& T_t,
                       CristaliniteAbstraite::SaveCrista * save_donTP
                       ,const double& P, const double& T,Enum_dure temps)
  { saveCrista =  save_donTP; // mise à jour au cas où des méthodes attachées utiliseraient le pointeur de la 
                              // classe générique
    Hoffman2::SaveCrista_Hoffman2* saveCristaHoffman2 = ((Hoffman2::SaveCrista_Hoffman2*) save_donTP);
    
    // on calcul par intégration avec la formule de Nakamura, la valeur de la cristalinité courante
    // pour cela on considère une intégration explicite en 1 pas de temps à partir de la forme différentielle
    // de l'équation cf. formule 9.12 de la thèse de Mikaël, page 139
    // on pourrait également utiliser une méthode de Runge, ce qui est possible si l'on veut plus de précision
//    double cris = saveCristaHoffman2->I_Kcinetique_t; // pour simplifier
    // recup de l'incrément de temps
	double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();
	double temps_courant = ParaGlob::Variables_de_temps().TempsCourant();

    // calcul de l'accroissement de la cristalinité sur le pas de temps par linéarisation de la forme différentielle
// non ne fonctionne pas car initialement la cristalinité = 0, et avec une formule purement explicite
// ce qui est le cas ici, la valeur du delta_cris = 0 donc la fonction n'arrive jamais à décoller !!
//    double delta_cris = deltat * (n * fct_KT(P,T) * (1.-cris)  
//                                    * pow( log( 1./(1.-cris) )
//                                           ,(n-1.)/n 
//                                          )
//                                   );
//    // mise à jour de la cristalinité
//    saveCristaHoffman2->cristalinite = cris + delta_cris;

     // utilisation d'une intégration  avec la méthode de Gauss
     // deux cas suivant que l'on calcul a t ou tdt
     double crista = 0.;
     if (temps == TEMPS_tdt)
     	{// dans ce cas on fait un calcul avec cumul à partir du temps précédent et du résultat
     	 // sauvegardé, qui est supposé à t
     	 // on renseigne les variables qui servent pour le calcul de la fonction a intégrer
     	 Te_t=T_t; Te_tdt=T;Pr_t=P_t; Pr_tdt=P;
     	 // on intègre la constante cinétique
         saveCristaHoffman2->I_Kcinetique = saveCristaHoffman2->I_Kcinetique_t +
                algo_integ ->IntegGauss(*this,&Hoffman2::fct_KT_,deltat);
//  cout << "\n I_Kcinetique= " << saveCristaHoffman2->I_Kcinetique;  //----- débug       
         // on calcul la cristalinité  
         crista = X_inf *(1.-exp(-pow(saveCristaHoffman2->I_Kcinetique,n)));    
     	}
     else // sinon pour l'instant ce n'est pas implanté, car il faut réfléchir
          // à la méthodologie d'intégration, pas difficile mais il vaut mieux le faire au moment
          // ou cela se présentera
        { cout << " \n   attention, cas non implante !!! pour l'instant (le demander si pb) "
               << " cas d'un calcul a TEMPS_t "
               << "\n Hoffman2::Cristalinite(..." ;
          Sortie(1);
        };
      
//--- debug
//  cout << fct_KT(P,T) << " ";
//  cout << delta_cris << " ";
//--- fin debug
                                                                    
    // retour
    return crista;
	};

	

// calcul du taux de cristalinité à tdt à partir des grandeurs à tdt
// P : pression à  t+deltat (= actuelle)
// T : température à  t+deltat (=actuelle)
double  Hoffman2::Cristalinite(CristaliniteAbstraite::SaveCrista * save_donTP
                         ,const double& , const double& )
  { saveCrista =  save_donTP; // mise à jour au cas où des méthodes attachées utiliseraient le pointeur de la 
                              // classe générique
    Hoffman2::SaveCrista_Hoffman2* saveCristaHoffman2 = ((Hoffman2::SaveCrista_Hoffman2*) save_donTP);
    
    // on calcul la cristalinité pour la valeur de la constante cinétique actuelle
    double crista = X_inf *(1.-exp(-pow(saveCristaHoffman2->I_Kcinetique,n)));    
                                                                    
    // retour
    return crista;
	};

//----- lecture écriture de restart spécifique aux données de la classe -----
// cas donne le niveau de la récupération
// = 1 : on récupère tout
// = 2 : on récupère uniquement les données variables (supposées comme telles)
void Hoffman2::Lecture_don_base_info(ifstream& ent,const int cas ,LesCourbes1D& ,LesFonctions_nD& )
 { // ici on ne lit que dans le cas = 1
   if (cas == 1)
   	{ string nom;
   	  ent >> nom;
   	  if (nom != "Cristallinite_Hoffman2")
   	   { cout << "\n erreur en lecture des donnees de la loi Hoffman2, on attendait le mot cle: "
   	          << " Cristallinite_Hoffman2 et on a lu: " << nom ;
   	     Sortie(1); 
   	   };
   	  // lecture des données
      ent >> nom >> n >> nom >> invt012 ;
      ent >> nom >> Ustar >> nom >> Tinf >> nom >> Kg >> nom >> Tm ;
      ent >> nom >> alphaP >> nom >> betaP >> nom >> X_inf;
   	};
 };	
 
// cas donne le niveau de sauvegarde
// = 1 : on sauvegarde tout
// = 2 : on sauvegarde uniquement les données variables (supposées comme telles)
void Hoffman2::Ecriture_don_base_info(ofstream& sort ,const int cas) const 
 { // ici on ne sauvegarde que dans le cas = 1
   if (cas == 1)
   	{ sort << "\n Cristallinite_Hoffman2 "
           << "\n     n= "      << n      << " invt012= "    << invt012    << " "
           << "\n     Ustar= "  << Ustar  << " Tinf= "  << Tinf  << " Kg= "  << Kg  << " Tm= "  << Tm << " "
           << "\n     alphaP= " << alphaP << " betaP= " << betaP << " X_inf= " << X_inf << " " ;
   	};
 };	
    
// --------------------------- METHODES PROTEGEES : ---------------------

// calcul de la fonction K(t) en fonction d'un paramètre qui varie de -1 à 1 et des variables
// définit par ailleurs: Te_t,Te_tdt,Pr_t,Pr_tdt
double  Hoffman2::fct_KT_(const double& thetai )
  { // la pression et la température sont interpolées linéairement sur le pas de temps
    double P = 0.5 * ((1.-thetai)*Pr_t + (1.+thetai)*Pr_tdt);
    double T = 0.5 * ((1.-thetai)*Te_t + (1.+thetai)*Te_tdt);
	return fct_KT(P,T);	
//    return 0.1;  //--- débug
	};