// FICHIER : Loi_newton3D.cc
// CLASSE : Loi_newton3D
// 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 "Debug.h"
# include <iostream>
using namespace std; //introduces namespace std
#include <math.h>
#include <stdlib.h>
#include "Sortie.h"
#include "TypeConsTens.h"
#include "ParaGlob.h"
#include "ConstMath.h"
#include "CharUtil.h"
#include "Loi_newton3D.h"
Loi_newton3D::Loi_newton3D () : // Constructeur par defaut
Loi_comp_abstraite(NEWTON3D,CAT_THERMO_MECANIQUE,3),mu(-ConstMath::trespetit),xn(1.),simple(true)
,mu_temperature(NULL),xn_temperature(NULL),Deps0(0.01)
{ };
// Constructeur de copie
Loi_newton3D::Loi_newton3D (const Loi_newton3D& loi) :
Loi_comp_abstraite(loi),mu(loi.mu),xn(loi.xn),simple(loi.simple)
,mu_temperature(loi.mu_temperature),xn_temperature(loi.xn_temperature)
,Deps0(loi.Deps0)
{// on regarde s'il s'agit d'une courbe locale ou d'une courbe globale
if (mu_temperature != NULL)
if (mu_temperature->NomCourbe() == "_")
mu_temperature = Courbe1D::New_Courbe1D(*(loi.mu_temperature));
if (xn_temperature != NULL)
if (xn_temperature->NomCourbe() == "_")
xn_temperature = Courbe1D::New_Courbe1D(*(loi.xn_temperature));;
};
Loi_newton3D::~Loi_newton3D ()
// Destructeur
{ if (mu_temperature != NULL)
if (mu_temperature->NomCourbe() == "_") delete mu_temperature;
if (xn_temperature != NULL)
if (xn_temperature->NomCourbe() == "_") delete xn_temperature;
};
// Lecture des donnees de la classe sur fichier
void Loi_newton3D::LectureDonneesParticulieres (UtilLecture * entreePrinc,LesCourbes1D& lesCourbes1D
,LesFonctions_nD& lesFonctionsnD)
{ string nom;
// lecture de la viscosité
*(entreePrinc->entree) >> nom;
if (nom != "mu=")
{ cout << "\n erreur en lecture de la viscosite, on aurait du lire le mot mu=";
entreePrinc->MessageBuffer("**erreur1 Loi_newton3D::LectureDonneesParticulieres (...**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
};
// on regarde si la viscosité est thermo dépendante
if(strstr(entreePrinc->tablcar,"mu_thermo_dependant_")!=0)
{ thermo_dependant=true;
*(entreePrinc->entree) >> nom;
if (nom != "mu_thermo_dependant_")
{ cout << "\n erreur en lecture de la thermodependance de mu, on aurait du lire le mot cle mu_thermo_dependant_"
<< " suivi du nom d'une courbe de charge ou de la courbe elle meme ";
entreePrinc->MessageBuffer("**erreur2 Loi_newton3D::LectureDonneesParticulieres (...**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
}
// lecture de la loi d'évolution en fonction de la température
*(entreePrinc->entree) >> nom;
// on regarde si la courbe existe, si oui on récupère la référence
if (lesCourbes1D.Existe(nom))
{ mu_temperature = lesCourbes1D.Trouve(nom);
}
else
{ // sinon il faut la lire maintenant
string non_courbe("_");
mu_temperature = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom.c_str()));
// lecture de la courbe
mu_temperature->LectDonnParticulieres_courbes (non_courbe,entreePrinc);
};
// prepa du flot de lecture
if(strstr(entreePrinc->tablcar,"fin_coeff_NEWTON3D")==0) entreePrinc->NouvelleDonnee();
}
else
{ // lecture de mu
*(entreePrinc->entree) >> mu ;
};
// on regarde s'il y a un coefficient non linéaire
simple = true; // par défaut
if (strstr(entreePrinc->tablcar,"xn=")!=NULL)
{ *(entreePrinc->entree) >> nom ;
simple = false;
#ifdef MISE_AU_POINT
if (nom != "xn=")
{ cout << "\n erreur en lecture de la loi de Newton, on attendait xn= suivi d'un nombre "
<< " ou du mot cle xn_thermo_dependant_ et ensuite une courbe et on a lu: "
<< nom <<" ";
entreePrinc->MessageBuffer("**erreur3 Loi_newton3D::LectureDonneesParticulieres (...**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
}
#endif
// on regarde si le coefficient non linéaire est thermo dépendant
if(strstr(entreePrinc->tablcar,"xn_thermo_dependant_")!=0)
{ thermo_dependant=true; *(entreePrinc->entree) >> nom;
if (nom != "xn_thermo_dependant_")
{ cout << "\n erreur en lecture de la thermodependance de xn, on aurait du lire le mot cle xn_thermo_dependant_"
<< " suivi du nom d'une courbe de charge ou de la courbe elle meme ";
entreePrinc->MessageBuffer("**erreur4 Loi_newton3D::LectureDonneesParticulieres (...**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
}
// lecture de la loi d'évolution du coefficient non linéaire en fonction de la température
*(entreePrinc->entree) >> nom;
// on regarde si la courbe existe, si oui on récupère la référence
if (lesCourbes1D.Existe(nom)) { xn_temperature = lesCourbes1D.Trouve(nom);}
else { // sinon il faut la lire maintenant
string non_courbe("_");
xn_temperature = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom.c_str()));
// lecture de la courbe
xn_temperature->LectDonnParticulieres_courbes (non_courbe,entreePrinc);
}
// prepa du flot de lecture
if(strstr(entreePrinc->tablcar,"fin_coeff_NEWTON3D")==0) entreePrinc->NouvelleDonnee();
}
else
{ // lecture du coeff
*(entreePrinc->entree) >> xn ;
};
};
// on regarde s'il y a un décalage à l'origine
simple = true; // par défaut
if (strstr(entreePrinc->tablcar,"Deps0=")!=NULL)
{ *(entreePrinc->entree) >> nom ;
#ifdef MISE_AU_POINT
if (nom != "Deps0=")
{ cout << "\n erreur en lecture de la loi de Newton, on attendait Deps0= suivi d'un nombre "
<< " et on a lu: " << nom <<" ";
entreePrinc->MessageBuffer("**erreur5 Loi_newton3D::LectureDonneesParticulieres (...**");
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
}
#endif
// lecture du coeff
*(entreePrinc->entree) >> Deps0 ;
};
// prepa du flot de lecture
if(strstr(entreePrinc->tablcar,"fin_coeff_NEWTON3D")==0) entreePrinc->NouvelleDonnee();
// appel au niveau de la classe mère
Loi_comp_abstraite::Lecture_type_deformation_et_niveau_commentaire
(*entreePrinc,lesFonctionsnD);
};
// affichage de la loi
void Loi_newton3D::Affiche() const
{ cout << " \n loi de comportement Newton 3D ";
if ( mu_temperature != NULL) { cout << " viscosite thermo dependant "
<< " courbe mu=f(T): " << mu_temperature->NomCourbe() <<" ";}
else { cout << " viscosite mu= " << mu ;}
if (!simple)
{ if (xn_temperature != NULL) { cout << " coef non lineaire thermo dependant "
<< " courbe xn=f(T): " << xn_temperature->NomCourbe() <<" ";}
else { cout << " coef non lin xn= " << xn ;}
};
cout << " Deps0= "<<Deps0<<endl;
// appel de la classe mère
Loi_comp_abstraite::Affiche_don_classe_abstraite();
};
// affichage et definition interactive des commandes particulières à chaques lois
void Loi_newton3D::Info_commande_LoisDeComp(UtilLecture& entreePrinc)
{ ofstream & sort = *(entreePrinc.Commande_pointInfo()); // pour simplifier
cout << "\n definition standart (rep o) ou exemples exhaustifs (rep n'importe quoi) ? ";
string rep = "_";
// procédure de lecture avec prise en charge d'un retour chariot
rep = lect_return_defaut(true,"o");
if (mu == -ConstMath::trespetit)
{ // on initialise à une valeur arbitraire
mu = 0.15;}
sort << "\n# ......................... loi de comportement Newton 3D ............................"
<< "\n# | viscosite | et eventuellement une puissance | puis eventuellement |"
<< "\n# | | pour une evolution non lineaire | une vitesse limite |"
<< "\n# | mu (obligatoire) | xn (facultatif) | inferieur (par defaut: 0.01)|"
<< "\n#......................................................................................"
<< "\n mu= " << setprecision(8) << mu << " xn= " << xn << " Deps0= "<<Deps0
<< "\n fin_coeff_NEWTON3D " << endl;
if ((rep != "o") && (rep != "O" ) && (rep != "0") )
{ sort
<< "\n# \n# l'equation constitutive est : S = mu * ((D_barre:D_barre+Deps0)^(xn/2) * D_barre "
<< "\n# \n# chaque parametre peut etre remplace par une fonction dependante de la temperature "
<< "\n# pour ce faire on utilise un mot cle puis une nom de courbe ou la courbe directement comme avec "
<< "\n# les autre loi de comportement "
<< "\n# exemple pour la viscosite: mu= mu_thermo_dependant_ courbe2 "
<< "\n# exemple pour la coef non lineaire: xn= xn_thermo_dependant_ courbe3 "
<< "\n# IMPORTANT: a chaque fois qu'il y a une thermodependence, il faut passer une ligne apres la description"
<< "\n# de la grandeur thermodependante, mais pas de passage à la ligne si se n'est pas thermo dependant "
<< "\n# la derniere ligne doit contenir uniquement le mot cle: fin_coeff_NEWTON3D "
<< endl;
};
// appel de la classe mère
Loi_comp_abstraite::Info_commande_don_LoisDeComp(entreePrinc);
};
// test si la loi est complete
int Loi_newton3D::TestComplet()
{ int ret = LoiAbstraiteGeneral::TestComplet();
if ((mu_temperature == NULL) && (mu == -ConstMath::trespetit))
{ cout << " \n la viscosite n'est pas defini pour la loi " << Nom_comp(id_comp)
<< '\n';
ret = 0;
}
return ret;
};
//----- lecture écriture de restart -----
// 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 Loi_newton3D::Lecture_base_info_loi(istream& ent,const int cas,LesReferences& lesRef,LesCourbes1D& lesCourbes1D
,LesFonctions_nD& lesFonctionsnD)
{ string nom;
if (cas == 1)
{ ent >> nom;
if (nom != "NEWTON3D")
{ cout << "\n erreur en lecture de la loi : Loi_newton3D, on attendait le mot cle : NEWTON3D "
<< "\n Loi_newton3D::Lecture_base_info_loi(...";
Sortie(1);
}
// ensuite normalement il n'y a pas de pb de lecture puisque c'est écrit automatiquement (sauf en debug)
// viscosité
bool test; ent >> nom >> test;
if (!test)
{ ent >> mu;
if (mu_temperature != NULL) {if (mu_temperature->NomCourbe() == "_") delete mu_temperature; mu_temperature = NULL;};
}
else
{ ent >> nom; mu_temperature = lesCourbes1D.Lecture_pour_base_info(ent,cas,mu_temperature); };
// le coef non linéaire
ent >> nom >> nom >> simple >> test;
if (!test)
{ ent >> xn;
if (xn_temperature != NULL) {if (xn_temperature->NomCourbe() == "_") delete xn_temperature; xn_temperature = NULL;};
}
else
{ ent >> nom; xn_temperature = lesCourbes1D.Lecture_pour_base_info(ent,cas,xn_temperature); };
ent >> nom >> Deps0;
};
Loi_comp_abstraite::Lecture_don_base_info(ent,cas,lesRef,lesCourbes1D,lesFonctionsnD);
};
// cas donne le niveau de sauvegarde
// = 1 : on sauvegarde tout
// = 2 : on sauvegarde uniquement les données variables (supposées comme telles)
void Loi_newton3D::Ecriture_base_info_loi(ostream& sort,const int cas)
{ if (cas == 1)
{ sort << " NEWTON3D " ;
sort << "\n viscosite ";
if (mu_temperature == NULL)
{ sort << false << " " << mu << " ";}
else
{ sort << true << " fonction_mu_temperature ";
LesCourbes1D::Ecriture_pour_base_info(sort,cas,mu_temperature);
};
sort << "\n non_lineaire " << " simple " << simple << " ";
if (xn_temperature == NULL)
{ sort << false << " " << xn << " ";}
else
{ sort << true << " fonction_xn_temperature ";
LesCourbes1D::Ecriture_pour_base_info(sort,cas,xn_temperature);
};
sort << " Deps0= "<<Deps0<<" ";
};
// appel de la classe mère
Loi_comp_abstraite::Ecriture_don_base_info(sort,cas);
};
// calcul d'un module d'young équivalent à la loi, ceci pour un
// chargement nul
double Loi_newton3D::Module_young_equivalent(Enum_dure temps,const Deformation & ,SaveResul * )
{ // ici on ne tiens pas compte de la non linéarité éventuelle, car on ne connait
// pas la vitesse de déformation !
/* const VariablesTemps& vartemps = ParaGlob::Variables_de_temps();
double temps = vartemps.TempsCourant();
if (temps <= ConstMath::trespetit)
// cas du temps nul
return ConstMath::tresgrand;
else
// cas du temps non nul
return mu/temps;*/
// non en fait il faudrait plus d'info, donc a minima on ramène mu, à charge aux prog appelant
// de gérer les conséquences
return 0.; //mu;
// en fait après réflexion, il n'y a pas de module d'young équivalent
// après re-re réflexion c'est le cas précédent qui est bon !! enfin je pense !
// mais le pb est qu'ici la partie newton ne s'intéresse qu'au cisaillement et que le temps critique
// intéresse aussi la partie volumique (peut-être même que la partie volumique) pour les ondes de compression
// or pour les ondes de compression, newton doit donner 0 !! car il intervient pas !!
// return ConstMath::trespetit;
};
// ========== codage des METHODES VIRTUELLES protegees:================
// calcul des contraintes a t+dt
void Loi_newton3D::Calcul_SigmaHH (TenseurHH& ,TenseurBB& DepsBB_,DdlElement & ,
TenseurBB & ,TenseurHH & ,BaseB& ,BaseH& ,TenseurBB& ,
TenseurBB& , TenseurBB& ,
TenseurHH & gijHH_,Tableau <TenseurBB *>& ,double& ,double& ,TenseurHH & sigHH_
,EnergieMeca & energ,const EnergieMeca & energ_t,double& module_compressibilite,double& module_cisaillement
,const Met_abstraite::Expli_t_tdt& )
{
#ifdef MISE_AU_POINT
if (DepsBB_.Dimension() != 3)
{ cout << "\nErreur : la dimension devrait etre 3 !\n";
cout << " Loi_newton3D::Calcul_SigmaHH\n";
Sortie(1);
};
#endif
const Tenseur3BB & DepsBB = *((Tenseur3BB*) &DepsBB_); // passage en dim 1
const Tenseur3HH & gijHH = *((Tenseur3HH*) &gijHH_); // " " " "
Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_); // " " " "
// cas de la thermo dépendance, on calcul les grandeurs
if (mu_temperature != NULL) mu = mu_temperature->Valeur(*temperature);
if (xn_temperature != NULL) xn = xn_temperature->Valeur(*temperature);
Tenseur3BH DepsBH = DepsBB * gijHH;
// calcul de la trace
double IDeps = DepsBH.Trace();
// le déviateur
Tenseur3BH Deps_barre_BH = DepsBH - (IDeps/3.) * IdBH3;
double d2= Deps_barre_BH && Deps_barre_BH; // inter
if (!simple && (d2>=ConstMath::trespetit) ) // cas d'une viscosité non linéaire et une vitesse non nulle
{ // cas d'une viscosité linéaire
Tenseur3BH sigBH = mu * Deps_barre_BH;
sigHH = gijHH * sigBH;
}
else
{ // cas d'une viscosité non linéaire
Tenseur3BH sigBH = (mu * pow((d2+Deps0),(0.5*xn))) * Deps_barre_BH;
sigHH = gijHH * sigBH;
};
// traitement des énergies
energ.Inita(0.);
// on considère une vitesse de déformation constante sur le pas de temps
double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant(); // recup de l'incrément de temps
energ.ChangeDissipationVisqueuse(energ_t.DissipationVisqueuse()+(sigHH && DepsBB)*deltat);
// --- débug -----
// if (energ.DissipationVisqueuse() > (1.e5))
// { cout << "\n attention l'energie visqueuse est bizarement grande !!! " << energ.DissipationVisqueuse();
// cout << "\n Loi_newton3D::Calcul_SigmaHH(... ";
// Sortie(1);
// };
// --- fin débug ----
////debug
//double toto = (sigHH && DepsBB)*deltat;
//cout << toto << " " << energ_t.DissipationVisqueuse() << " " << energ.DissipationVisqueuse() << " " << endl;
//if (toto < 0.)
// { cout << "\n debug: Loi_newton3D::Calcul_SigmaHH "
// << " d visqueux = "<< toto << endl;
//
// };
////fin debug
// on libère les tenseurs intermédiaires
LibereTenseur();
};
// calcul des contraintes a t+dt et de ses variations
void Loi_newton3D::Calcul_DsigmaHH_tdt (TenseurHH& ,TenseurBB& DepsBB_,DdlElement & tab_ddl
,BaseB& ,TenseurBB & ,TenseurHH & ,
BaseB& ,Tableau <BaseB> & ,BaseH& ,Tableau <BaseH> & ,
TenseurBB & ,Tableau <TenseurBB *>& d_epsBB,TenseurBB & ,
TenseurBB & ,TenseurHH & gijHH_tdt,
Tableau <TenseurBB *>& d_gijBB_tdt,
Tableau <TenseurHH *>& d_gijHH_tdt,double& , double& ,
Vecteur& ,TenseurHH& sigHH_tdt,Tableau <TenseurHH *>& d_sigHH
,EnergieMeca & energ,const EnergieMeca & energ_t,double& module_compressibilite,double& module_cisaillement
,const Met_abstraite::Impli& )
{
#ifdef MISE_AU_POINT
if (DepsBB_.Dimension() != 3)
{ cout << "\nErreur : la dimension devrait etre 3 !\n";
cout << " Loi_newton3D::Calcul_DsigmaHH_tdt\n";
Sortie(1);
};
if (tab_ddl.NbDdl() != d_gijBB_tdt.Taille())
{ cout << "\nErreur : le nb de ddl est != de la taille de d_gijBB_tdt !\n";
cout << " Loi_newton3D::Calcul_SDsigmaHH_tdt\n";
Sortie(1);
};
#endif
const Tenseur3BB & DepsBB = *((Tenseur3BB*) &DepsBB_); // passage en dim 1
const Tenseur3HH & gijHH = *((Tenseur3HH*) &gijHH_tdt); // " " " "
Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_tdt); // " " " "
// cas de la thermo dépendance, on calcul les grandeurs
if (mu_temperature != NULL) mu = mu_temperature->Valeur(*temperature);
if (xn_temperature != NULL) xn = xn_temperature->Valeur(*temperature);
Tenseur3BH DepsBH = DepsBB * gijHH;
// calcul de la trace
double IDeps = DepsBH.Trace();
// le déviateur
Tenseur3BH Deps_barre_BH = DepsBH - (IDeps/3.) * IdBH3;
double d2= Deps_barre_BH && Deps_barre_BH; // inter, nombre positif
Tenseur3BH sigBH = mu * Deps_barre_BH; // cas d'une viscosité linéaire
if (!simple && (d2>=ConstMath::trespetit) ) // cas d'une viscosité non linéaire et une vitesse non nulle
{sigBH *= pow((d2+Deps0),(0.5*xn)) ;}
sigHH = gijHH * sigBH;
// cout << "\n d2= " << d2 << " sigB= " << sqrt(3./2.*sigBH.II()) << " ";
int nbddl = d_gijBB_tdt.Taille();
// recup de l'incrément de temps
double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();
double unSurDeltat;
if (Abs(deltat) >= ConstMath::trespetit)
{unSurDeltat = 1./deltat;}
else
// si l'incrément de temps est tres petit on fait comme pour le module d'young équivalent
// on considère une raideur petite (pas complèment nul quand même pour éviter une singularité)
// { unSurDeltat = Signe(deltat)*ConstMath::trespetit;}; // ce qui donnera un terme très petit
// si l'incrément de temps est tres petit on remplace 1/deltat par un nombre tres grand
{ // un pas de temps doit être positif !! or certaine fois il peut y avoir des pb
if (unSurDeltat < 0)
{ cout << "\n le pas de temps est négatif !! "; };
unSurDeltat = ConstMath::tresgrand;
};
// boucle sur les ddl
Tenseur3BH dDepsBH,d_Deps_barre_BH; // tenseurs de travail
for (int i = 1; i<= nbddl; i++)
{// on fait uniquement une égalité d'adresse de manière à ne pas utiliser
// le constructeur d'ou la profusion d'* et de ()
Tenseur3HH & dsigHH = *((Tenseur3HH*) (d_sigHH(i))); // passage en dim 3
// const Tenseur3BB & dgijBB = *((Tenseur3BB*)(d_gijBB_tdt(i))); // passage en dim 3
const Tenseur3HH & dgijHH = *((Tenseur3HH*)(d_gijHH_tdt(i))) ; // pour simplifier l'ecriture
const Tenseur3BB & depsBB = *((Tenseur3BB *) (d_epsBB(i))); // "
// variation de la vitesse de déformation
dDepsBH = ( unSurDeltat) * depsBB * gijHH_tdt + DepsBB * dgijHH;
// variation de la trace
double d_IDeps = dDepsBH.Trace();
// variation du déviateur
d_Deps_barre_BH = dDepsBH - (d_IDeps/3.) * IdBH3;
if ((simple)||(d2<=ConstMath::trespetit)) // cas d'une viscosité linéaire ou vitesse nulle
{ dsigHH = dgijHH * sigBH + gijHH * mu * d_Deps_barre_BH;}
else // cas d'une viscosité non linéaire
{ double dDHB_DHB = ((d_Deps_barre_BH && Deps_barre_BH) + (Deps_barre_BH && d_Deps_barre_BH)) ;
//d_Deps_barre_BH && Deps_barre_BH;
dsigHH = dgijHH * sigBH
+ gijHH * (
(mu * xn * 0.5 * pow((d2+Deps0),(xn*0.5-1.))*dDHB_DHB) * Deps_barre_BH
+ (mu * pow((d2+Deps0),(0.5*xn))) * d_Deps_barre_BH
);
}
}; //-- fin de la boucle sur les ddl
// traitement des énergies
energ.Inita(0.);
// on considère une vitesse de déformation constante sur le pas de temps
energ.ChangeDissipationVisqueuse(energ_t.DissipationVisqueuse()+(sigHH && DepsBB)*deltat);
// on libère les tenseurs intermédiaires
LibereTenseur();
};
// calcul des contraintes et ses variations par rapport aux déformations a t+dt
// en_base_orthonormee: le tenseur de contrainte en entrée est en orthonormee
// le tenseur de déformation et son incrémentsont également en orthonormees
// si = false: les bases transmises sont utilisées
// ex: contient les éléments de métrique relativement au paramétrage matériel = X_(0)^a
void Loi_newton3D::Calcul_dsigma_deps (bool en_base_orthonormee, TenseurHH & ,TenseurBB& DepsBB_
,TenseurBB & epsBB_tdt,TenseurBB & ,double& ,double&
,TenseurHH& sigHH_tdt,TenseurHHHH& d_sigma_deps_
,EnergieMeca & energ,const EnergieMeca & energ_t,double& module_compressibilite,double& module_cisaillement
,const Met_abstraite::Umat_cont& )
{
#ifdef MISE_AU_POINT
if (epsBB_tdt.Dimension() != 3)
{ cout << "\nErreur : la dimension devrait etre 3 !\n";
cout << " Loi_newton3D::Calcul_dsigma_deps\n";
Sortie(1);
};
#endif
const Tenseur3BB & epsBB = *((Tenseur3BB*) &epsBB_tdt); // passage en dim 3
Tenseur3HH & sigHH = *((Tenseur3HH*) &sigHH_tdt); // " " " "
Tenseur3BB& DepsBB = *((Tenseur3BB*) &DepsBB_);
// cas de la thermo dépendance, on calcul les grandeurs
if (mu_temperature != NULL) mu = mu_temperature->Valeur(*temperature);
if (xn_temperature != NULL) xn = xn_temperature->Valeur(*temperature);
// cas du tenseur des contraintes
Tenseur3BH DepsBH = DepsBB.MonteDernierIndice();
// calcul de la trace
double IDeps = DepsBH.Trace();
// le déviateur
Tenseur3BH Deps_barre_BH = DepsBH - (IDeps/3.) * IdBH3;
Tenseur3BB Deps_barre_BB = Deps_barre_BH * IdBB3;
Tenseur3HH Deps_barre_HH = IdHH3 * Deps_barre_BH ;
double d2= Deps_barre_BH && Deps_barre_BH; // inter, nombre positif
Tenseur3BH sigBH = mu * Deps_barre_BH; // cas d'une viscosité linéaire
if (!simple && (d2>=ConstMath::trespetit) ) // cas d'une viscosité non linéaire et une vitesse non nulle
{sigBH *= pow((d2+Deps0),(0.5*xn)) ;}
sigHH = IdHH3 * sigBH;
// recup de l'incrément de temps
double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();
double unSurDeltat;
if (Abs(deltat) >= ConstMath::trespetit)
{unSurDeltat = 1./deltat;}
else
// si l'incrément de temps est tres petit on remplace 1/deltat par un nombre tres grand
{ // un pas de temps doit être positif !! or certaine fois il peut y avoir des pb
if (unSurDeltat < 0)
{ cout << "\n le pas de temps est négatif !! "; };
unSurDeltat = ConstMath::tresgrand;
};
// cas le la variation du tenseur des contraintes par rapport aux déformations
Tenseur3HHHH & d_sigma_depsHHHH = *((Tenseur3HHHH*) &d_sigma_deps_);
// on commence par calculer la variation de la vitesse de déformation barre / à eps
// (d_D_ij)barre / d_eps_kl = 1/deltat (delta_i^k * delta_j^l
// - 1/3 * delta_ij * g^nm * delta_m^k * delta_n^l)
// ici g^nm = delta^nm, et g^nm * delta_m^k * delta_n^l = id * id = id = delta^kl
static const double untiers = 1./3.;
Tenseur3HHHH d_D_barre_deps_HHHH =
((Tenseur3BBBB*) &(unSurDeltat * (PIdBBBB3 - untiers * IdBBBB3)))->Monte4Indices();
// puis calcul en fonction du test sur "simple"
if ((simple)||(d2<=ConstMath::trespetit)) // cas d'une viscosité linéaire ou vitesse nulle
{d_sigma_depsHHHH = mu * d_D_barre_deps_HHHH;}
else // cas d'une viscosité non linéaire
{ Tenseur3HH d_D2_HH= (d_D_barre_deps_HHHH && Deps_barre_BB)
+ (Deps_barre_BB && d_D_barre_deps_HHHH);
d_sigma_depsHHHH = mu * (d_D_barre_deps_HHHH
+ (xn * 0.5 * pow((d2+Deps0),(xn*0.5-1.)))
* Tenseur3HHHH::Prod_tensoriel(Deps_barre_HH,d_D2_HH));
};
// traitement des énergies
energ.Inita(0.);
// on considère une vitesse de déformation constante sur le pas de temps
energ.ChangeDissipationVisqueuse(energ_t.DissipationVisqueuse()+(sigHH && DepsBB)*deltat);
// on libère les tenseurs intermédiaires
LibereTenseur(); LibereTenseurQ();
};