mg_calcul_fatigue.cpp

Aller à la documentation de ce fichier.

//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//  MAGiC
//####//  Jean Christophe Cuilliere et Vincent FRANCOIS
//####//  Departement de Genie Mecanique - UQTR
//####//------------------------------------------------------------
//####//  MAGIC est un projet de recherche de l equipe ERICCA
//####//  du departement de genie mecanique de l Universite du Quebec a Trois Rivieres
//####//  http://www.uqtr.ca/ericca
//####//  http://www.uqtr.ca/
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//
//####//       mg_calcul_fatigue.cpp
//####//
//####//------------------------------------------------------------
//####//------------------------------------------------------------
//####//    COPYRIGHT 2000-2024
//####//  jeu 13 jun 2024 11:58:57 EDT
//####//------------------------------------------------------------
//####//------------------------------------------------------------
#include "gestionversion.h"
#include "mg_calcul_fatigue.h"
#include "mg_file.h"
#include "mg_volume.h"
#include "fem_solution.h"
#include "ot_parametres.h"
#include <fct_generateur_calibrage.h>
#include "fem_solution.h"
#include <cmath>
#include "ot_parametres.h"
#include <limits>
 
DAMAGE_FEM3::DAMAGE_FEM3(FEM_ELEMENT3* ele3)
{
  m_ele3=ele3;
}
 
DAMAGE_FEM3::~DAMAGE_FEM3(void)
{
 
}
 
FEM_ELEMENT3* DAMAGE_FEM3::get_fem_element3(void)
{
  return m_ele3;
}
 
void DAMAGE_FEM3::change_fem_element3(FEM_ELEMENT3* ele3)
{
  m_ele3=ele3;
}
 
double DAMAGE_FEM3::get_damage(void)
{
  return m_damage;
}
 
void DAMAGE_FEM3::change_damage(double damage)
{
  m_damage=damage;
}
 
int DAMAGE_FEM3::get_niveau_maille(void)
{
  return m_niveau_maille;
}
 
void DAMAGE_FEM3::change_niveau_maille(int niveau)
{
  m_niveau_maille=niveau;
}
 
 
MG_CALCUL_FATIGUE::MG_CALCUL_FATIGUE():affichageactif(0)
{
}
 
MG_CALCUL_FATIGUE::MG_CALCUL_FATIGUE(MG_CALCUL_FATIGUE& mdd):affichageactif(mdd.affichageactif)
{
}
 
MG_CALCUL_FATIGUE::~MG_CALCUL_FATIGUE()
{
}
 
void MG_CALCUL_FATIGUE::active_affichage(void (*fonc)(char*))
{
affiche=fonc;
affichageactif=1;
}
 
char MG_CALCUL_FATIGUE::get_seprateur_decimale(void)
{
OT_PARAMETRES param;
param.ajouter("Separateur_decimale",",",OT_PARAMETRES::STRING,"Seprateur décimale pour l'exportation des fichiers excel");
char fichier[500];
sprintf(fichier,"%s/.magic",getenv("HOME"));
param.lire(fichier);
param.enregistrer(fichier);
char c=(param.get_nom("Separateur_decimale").c_str())[0];
return c;
}
 
int MG_CALCUL_FATIGUE::cycle_jump(MG_GESTIONNAIRE &gest,MG_GEOMETRIE *geo,int &Njump,char *fichierparam,int numchamps,int N,int niveaumax,char *file,std::vector<FEM_ELEMENT3*> *lst,std::vector<double>& tableau_D,int nombre_classe,double &max_damage_element,FILE *strength,double x,double y,double z,char *fichiercalcul,int pos,double val_pos,int num_noeud)
{
double alpha,beta,epsilonint,sigma,epsilon,D,E,val_CM,varD,coef1,coef2,coef3,coef4,coef5,rapport,sigmamax,Njump_Gauss,deltaD,coef_model;
int nomb,l,coord,col1,s,ll,l1,lo;
l=0;
ll=0;
l1=0;
s=0;
coord=0;
lo=0;
OT_PARAMETRES* parametre=new OT_PARAMETRES;
parametre->lire(fichiercalcul);
sigmamax= parametre->get_valeur((char*)"sigmamax");
alpha= parametre->get_valeur((char*)"alpha");
beta= parametre->get_valeur((char*)"beta");
epsilonint= parametre->get_valeur((char*)"epsilonint");
E= parametre->get_valeur((char*)"E");
nomb= parametre->get_valeur((char*)"nomb");
char ligne[300];
char ligne1[300];
double coef[10];
FILE *coefficient_optim = NULL;
coefficient_optim= fopen((char*)"resistance.coe", "r");
while (fgets(ligne,100,coefficient_optim) !=NULL)
{
sscanf (ligne,"%s %lf",ligne1,&coef_model);
coef[lo]=coef_model;
lo++;
}
fclose(coefficient_optim); 
FEM_MAILLAGE* fem;
FEM_SOLUTION* sol;
for(long i=0;i<geo->get_nb_mg_volume();i++)
{
recupere_maillage(gest,&fem);  
OT_PARAMETRES param_aster; 
param_aster.lire(fichierparam);
int k=param_aster.get_valeur("It_resultat");
recupere_resultat(gest,fem,&sol,k);
LISTE_FEM_ELEMENT3::iterator it;
int element=fem->get_nb_fem_element3();
int nombre_Gauss_max=15;
std::vector<double> tableau_Njump_Gauss;
std::vector<double> tableau_varD;
std::vector<double> tableau_freq;
int taille=element*nombre_Gauss_max;
tableau_Njump_Gauss.reserve(taille);
tableau_varD.reserve(taille);
tableau_freq.reserve(taille);
long int id=get_id_element_face(fem,x,y,z,pos,val_pos,num_noeud);
 for (FEM_ELEMENT3* ele=fem->get_premier_element3(it);ele!=NULL;ele=fem->get_suivant_element3(it))
 {
      unsigned long id1=ele->get_id();
      int nbgauss=sol->get_nb_gauss(l);
      for (int r=0;r<nbgauss;r++)
        {   
         if (N==1) 
         {
           D=0;
           tableau_D[s]=0;  
         }  
         else {D=tableau_D[s];}
        double valeur=sol->lire(l,numchamps,coord,r);
      //************************Modele d'endommagement******************************//
        rapport=fabs(valeur/sigmamax);
        //varD=coef[0]*rapport*exp(-coef[1]*(D/sqrt(rapport)))+coef[2]*D*pow(rapport,2)*(1+exp(coef[4]*(rapport-coef[3])));
        // varD=coef[0]*rapport*exp(-coef[1]*(D/sqrt(rapport)))+coef[2]*D*pow(rapport,2)*(1+exp(coef[3]*rapport));
        varD=coef[0]*rapport*exp(-coef[1]*(D/sqrt(rapport)))+coef[2]*D*pow(rapport,coef[3]);
        //varD=rapport*(coef[0]*exp(-coef[1]*(D/sqrt(rapport)))+coef[2]*D*exp(-1/sqrt(rapport)));
        tableau_varD[s]=varD;
     //************************Modele d'endommagement******************************//
        if (D==0) 
        {
         // deltaD=1E-09;
          deltaD=1E-10; 
        }
        if (D>0 && D<=0.2) 
        {
          //deltaD=0.5*D;
          deltaD=0.2*D;
            //deltaD=0.9*D;
        }
        if (D>0.2)
        {
         // deltaD=0.1;
          deltaD=0.1*D;
           // deltaD=0.05*D;
        }
        Njump_Gauss=round(deltaD/varD);
        tableau_Njump_Gauss[s]=Njump_Gauss;
        s++;
        }
        if (id1==id)
        {
        double valeur_element=sol->get_moyenne_volumique_elementaire_champs(numchamps,ele,l,coord);
        //BOITE_3D boite3d = ele->get_boite_3D();
        //fprintf(strength,"%d %lf centre_boite_3d (%lf %lf %lf)\n",N,valeur_element,boite3d.get_xcentre(),boite3d.get_ycentre(),boite3d.get_zcentre());
        fprintf(strength,"%d %lf\n",N,valeur_element);
        }
        l++;
 } 
 
 
double max_Njump_Gauss,freq_relative,freq_cumulative,low_freq,classe;
int co1,nomb_jump;
co1=0;
freq_cumulative=0;
tableau_freq[0]=0;
for(int j=0;j<s;j++)
{ 
  if(tableau_Njump_Gauss[j]>co1)
  {
    max_Njump_Gauss=tableau_Njump_Gauss[j];
    co1=tableau_Njump_Gauss[j];
  }
  else 
  {
    max_Njump_Gauss=co1;  
  }
}
 
 
classe=max_Njump_Gauss/nombre_classe;  
 
 
for(int m=0;m<nombre_classe;m++)
{
  nomb_jump=0;
  for(int jj=0;jj<s;jj++)
  {
    if (((classe*m)<=tableau_Njump_Gauss[jj]) && (tableau_Njump_Gauss[jj]<(classe*(m+1))))
    {
      nomb_jump++;
    }
    if ((m==nombre_classe-1) && (tableau_Njump_Gauss[jj]==max_Njump_Gauss))    //inclure les max_Njump_Gauss dans la derniere classe
    {
      nomb_jump++;   
    }
   
  }
  freq_relative=(double)nomb_jump/s;
  freq_cumulative=freq_cumulative+freq_relative; 
  tableau_freq[m+1]=freq_cumulative;
}
low_freq=0.1*freq_cumulative;
for(int m=0;m<nombre_classe;m++)
{
   if ((low_freq>tableau_freq[m]) && (low_freq<tableau_freq[m+1]))
   {
    Njump=(int)round(((m+1)*classe)-((tableau_freq[m+1]-low_freq)*((((m+1)*classe)-(m*classe))/(tableau_freq[m+1]-tableau_freq[m])))+1);      //interpolation lineaire
   }
}
int n;
double damage_element,co2,co3,min_damage_element;
n=0;
co2=0;
co3=1;
char message1[500];
char message2[500];
sprintf(message1,"endommagement.sol");
sprintf(message2,"Nom");
FEM_SOLUTION* sol1=new FEM_SOLUTION(fem,1,message1,fem->get_nb_fem_element3(),message2,MAGIC::ENTITE_SOLUTION::ENTITE_ELEMENT3_GAUSS,MAGIC::TYPE_SOLUTION::SCALAIRE);
sol1->change_legende(0,"Fatigue_damage");
std::vector<DAMAGE_FEM3*> lst_damage_fem3;
for (FEM_ELEMENT3* ele=fem->get_premier_element3(it);ele!=NULL;ele=fem->get_suivant_element3(it))
 {
      int nbgauss=sol->get_nb_gauss(ll);
      for (int r=0;r<nbgauss;r++)
        {  
        tableau_D[n]=tableau_D[n]+(tableau_varD[n]*Njump);
        sol1->ecrire(tableau_D[n],ll,0);
        n++;
        }
        damage_element=sol1->get_moyenne_volumique_elementaire_champs(numchamps,ele,ll,coord);
         if(damage_element>co2)
         {
          max_damage_element=damage_element;
          co2=damage_element;
         }
         else 
         {
          max_damage_element=co2;
         }
         if(damage_element<co3)
         {
          min_damage_element=damage_element;
          co3=damage_element;
         }
         else 
         {
          min_damage_element=co3;
         }
        DAMAGE_FEM3 *damagefem3=new DAMAGE_FEM3(ele);
        damagefem3->change_damage(damage_element);
        lst_damage_fem3.push_back(damagefem3);
        ll++;       
 }
delete sol1;
std::vector<FEM_ELEMENT3*> *lstniveau=new std::vector<FEM_ELEMENT3*>[niveaumax+1];
std::vector<DAMAGE_FEM3*>::iterator it_damage_ele;
int l2=0;
for(it_damage_ele=lst_damage_fem3.begin();it_damage_ele!=lst_damage_fem3.end();it_damage_ele++)
{
  DAMAGE_FEM3 *damagefem3=lst_damage_fem3[l2];
  damagefem3->change_niveau_maille((int)((damagefem3->get_damage()-min_damage_element)*niveaumax/(max_damage_element-min_damage_element))+1);
  if ((damagefem3->get_niveau_maille())>niveaumax) damagefem3->change_niveau_maille(niveaumax);
  lstniveau[damagefem3->get_niveau_maille()].push_back(damagefem3->get_fem_element3()); 
  lst[damagefem3->get_niveau_maille()]=lstniveau[damagefem3->get_niveau_maille()];
  l2++;
}
delete [] lstniveau;
double Dzone; 
MG_GEOM_FONCTION *fonction;
double coo[2];
if (max_damage_element<=0.99999)
{
for(int l3=1;l3<=niveaumax;l3++)
{
  Dzone=((((l3*(max_damage_element-min_damage_element))/niveaumax)+(((l3-1)*(max_damage_element-min_damage_element))/niveaumax))/2)+min_damage_element;
  val_CM=geo->get_mg_geom_fonction((int)l3-1)->get_id();
  fonction = geo->get_mg_geom_fonctionid((int)val_CM);
  for(int k=0;k<nomb;k++)
    {
      epsilon=epsilonint+(k*0.0010);
      sigma=exp(alpha*log(epsilon)+log(fabs(E*(1-Dzone)))+beta);     //enlever la valeur absolu quand Dzone sera correct et inferieur a 1
      coo[0]=epsilon;
      coo[1]=sigma;
      fonction->change_point(k,coo);
      gest.enregistrer(file);
    }  
}
 
}
}
return 0;
}
int MG_CALCUL_FATIGUE::recupere_maillage(MG_GESTIONNAIRE &gest,FEM_MAILLAGE **femmesh)
{
 
FEM_MAILLAGE *fem=gest.get_fem_maillage(0);
*femmesh=fem;
 
return 0;
}
 
int MG_CALCUL_FATIGUE::recupere_resultat(MG_GESTIONNAIRE &gest,FEM_MAILLAGE *fem,FEM_SOLUTION **solmesh,int k)
{  
FEM_SOLUTION *sol=gest.get_fem_solution(k);
*solmesh=sol;
return 0;
}
 
long int MG_CALCUL_FATIGUE::get_id_element(FEM_MAILLAGE *fem, double x, double y, double z)
{
  long id=-1;
  double distance_min=std::numeric_limits<double>::max();
  LISTE_FEM_ELEMENT3::iterator it_ele;
  for(FEM_ELEMENT3* ele=fem->get_premier_element3(it_ele);ele!=NULL;ele=fem->get_suivant_element3(it_ele))
  {
    BOITE_3D boite = ele->get_boite_3D();
    double centre[3];
    boite.get_centre(centre);
    double distance = sqrt((centre[0]-x)*(centre[0]-x)+(centre[1]-y)*(centre[1]-y)+(centre[2]-z)*(centre[2]-z));
    if(distance<distance_min)
    {
      distance_min=distance;
      id = ele->get_id();
    }
  }
  return id;
}
 
long int MG_CALCUL_FATIGUE::get_id_element_face(FEM_MAILLAGE *fem, double x, double y, double z, int position, double valeur_position, int nomb_noeud)
{
  long id=-1;
  double distance_min=std::numeric_limits<double>::max();
  LISTE_FEM_ELEMENT3::iterator it_ele;
  for(FEM_ELEMENT3* ele=fem->get_premier_element3(it_ele);ele!=NULL;ele=fem->get_suivant_element3(it_ele))
  {
    BOITE_3D boite = ele->get_boite_3D();
    double centre[3];
    boite.get_centre(centre);
    double distance = sqrt((centre[0]-x)*(centre[0]-x)+(centre[1]-y)*(centre[1]-y)+(centre[2]-z)*(centre[2]-z));
    FEM_NOEUD *noeud;
    double *corr;
    int cmax=ele->get_nb_fem_noeud();
    int compteur=0;
    for(int c=0;c<cmax;c++)
        {
         noeud=ele->get_fem_noeud(c); 
         corr=noeud->get_coord();
         if (corr[position]==valeur_position){compteur++;};
        }
    if((distance<distance_min) && (compteur==nomb_noeud))
        {
         distance_min=distance;
         id = ele->get_id();
        }
  }
  return id;
}
 
void MG_CALCUL_FATIGUE::modification_courbe_exp_sidolo(FILE *strength)
{
 char ligne2[300];
 char ligne3[300];
 char mot1[300];
 char mot2[300];
 int nbligne=0;
 int c,val3;
 double val1,val2,val4,val5;
 FILE *strength_exp = NULL;
 strength_exp = fopen((char*)"resistanceavantmodif.exp", "r");
 while((c=fgetc(strength_exp)) != EOF)
 {
  if (c == '\n') {nbligne++;}
 }
 double vec[nbligne];
 double vec1[nbligne];
 int f=0;
 fclose(strength_exp);
 strength_exp = fopen((char*)"resistanceavantmodif.exp", "r");
 fgets(ligne2,100,strength_exp);
 sscanf (ligne2,"%s %s",mot1,mot2);
 while (fgets(ligne2,100,strength_exp) !=NULL)
 {
   sscanf (ligne2,"%lf %lf",&val1,&val2);
   vec[f]=val1;
   vec1[f]=val2;
   f++;
 }
 fclose(strength_exp);
 FILE *strength_exp_modif = NULL;
 strength_exp_modif = fopen((char*)"resistance.exp", "w");
 strength = fopen((char*)"resistancenum.txt", "r"); 
 fprintf(strength_exp_modif,"%s %s\n",mot1,mot2);
 while (fgets(ligne3,100,strength) !=NULL)
  {
   sscanf (ligne3,"%d %lf",&val3,&val4); 
   double val=(double)val3;
   for(f=0;f<nbligne;f++)
   {
     if ((val3>vec[f]) && (val3<vec[f+1]))
     {
      val5=((vec1[f]-vec1[f+1])/(vec[f]-vec[f+1]))*(val3-vec[f+1])+vec1[f+1];      //interpolation lineaire
      fprintf(strength_exp_modif,"%lf %lf\n",val,val5);
     }
     if (val3==vec[f]) 
     {
      val5=vec1[f];
      fprintf(strength_exp_modif,"%lf %lf\n",val,val5);
     }
  
   }
  }
 fclose(strength);
 fclose(strength_exp_modif);
}
 
void MG_CALCUL_FATIGUE::parametres_calcul(OT_PARAMETRES* params,char* fichier_params)
{
  params->ajouter("Nmax","5000000.0",OT_PARAMETRES::DOUBLE,"Nombre de cycles maximal");
  params->ajouter("nombreclasse","800.0",OT_PARAMETRES::DOUBLE,"Nombre de classe dans le calcul de la frequence cumulative");
  params->ajouter("valx","63.275",OT_PARAMETRES::DOUBLE,"coordonnee 1 du point localisant un element 3D");  
  params->ajouter("valy","-6.22",OT_PARAMETRES::DOUBLE,"coordonnee 2 du point localisant un element 3D");
  params->ajouter("valz","0.0",OT_PARAMETRES::DOUBLE,"coordonnee 3 du point localisant un element 3D");
  params->ajouter("sigmamax","34.35E+06",OT_PARAMETRES::DOUBLE,"Contrainte maximale dans un essai de traction");
  params->ajouter("alpha","0.6473738171",OT_PARAMETRES::DOUBLE,"parametre 1 du modele qui definie la courbe de traction");
  params->ajouter("beta","-2.096345371",OT_PARAMETRES::DOUBLE,"parametre 2 du modele qui definie la courbe de comportement");
  params->ajouter("nomb","13.0",OT_PARAMETRES::DOUBLE,"Nombre de points dans la courbe de comportement");
  params->ajouter("epsilonint","0.002",OT_PARAMETRES::DOUBLE,"Premiere deformation dans la courbe de comportement");
  params->ajouter("E","3.38E+09",OT_PARAMETRES::DOUBLE,"Module Young");
  params->ajouter("nbrniveau","4.0",OT_PARAMETRES::DOUBLE,"Nombre de niveau");
  params->ajouter("position","2.0",OT_PARAMETRES::DOUBLE,"Coordonnee du noeud selon laxe perpendiculaire au plan qui contient une face de l element se trouvant dans la boite");
  params->ajouter("valeur_position","0.0",OT_PARAMETRES::DOUBLE,"Position du plan contenant une face de l element se trouvant dans la boite");
  params->ajouter("nomb_noeud","6.0",OT_PARAMETRES::DOUBLE,"Nombre de noeuds appartenant au plan qui contient une face de l element se trouvant dans la boite");
}