rio/Herezh_dev
1
0
Fork 0
Code Issues Pull requests Projects Releases Wiki Activity
Herezh_dev/comportement/plasticite/Prandtl_Reuss2D_D.cc

691 lines
33 KiB
C++
Raw Normal View History

ajout des lois de comportement
2021-09-23 11:21:15 +02:00
// FICHIER : Loi_iso_elas.cp
// CLASSE : Loi_iso_elas
// 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-2021 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 "ConstMath.h"
#include "Prandtl_Reuss2D_D.h"
// ========== fonctions pour la classe de sauvegarde des résultats =========
// affectation
Loi_comp_abstraite::SaveResul &
Prandtl_Reuss2D_D::SaveResulPrandtl_Reuss2D_D::operator = ( const Loi_comp_abstraite::SaveResul & a)
{ Prandtl_Reuss2D_D::SaveResulPrandtl_Reuss2D_D& sav
= *((Prandtl_Reuss2D_D::SaveResulPrandtl_Reuss2D_D*) &a);
// données protégées
epsilon_barre = sav.epsilon_barre;
def_plasBB = sav.def_plasBB;
epsilon_barre_t = sav.epsilon_barre_t;
def_plasBB_t = sav.def_plasBB_t;
return *this;
};
//------- 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 Prandtl_Reuss2D_D::SaveResulPrandtl_Reuss2D_D::Lecture_base_info
(ifstream& ent,const int )
{ // ici toutes les données sont toujours a priori variables
string toto;
ent >> toto >> epsilon_barre_t;
ent >> toto >> def_plasBB_t ;
};
// cas donne le niveau de sauvegarde
// = 1 : on sauvegarde tout
// = 2 : on sauvegarde uniquement les données variables
//(supposées comme telles)
void Prandtl_Reuss2D_D::SaveResulPrandtl_Reuss2D_D::Ecriture_base_info
(ofstream& sort,const int )
{ // ici toutes les données sont toujours a priori variables
sort << " epsb_t " << epsilon_barre_t << " " ;
sort << " def_plasBB_t " << def_plasBB_t << " ";
};
// affichage à l'écran des infos
void Prandtl_Reuss2D_D::SaveResulPrandtl_Reuss2D_D::Affiche() const
{ cout << "\n SaveResulPrandtl_Reuss2D_D: " ;
cout << "\n epsilon_barre= " << epsilon_barre << " def_plasBB= " << def_plasBB
<< " epsilon_barre_t= " << epsilon_barre_t << " def_plasBB_t= " << def_plasBB_t;
cout << " ";
};
// ========== fin des fonctions pour la classe de sauvegarde des résultats =========
Prandtl_Reuss2D_D::Prandtl_Reuss2D_D () : // Constructeur par defaut
Loi_comp_abstraite(PRANDTL_REUSS2D_D,CAT_MECANIQUE,2),E(0.),nu(-2.)
,f_ecrouissage(NULL)
,tolerance_plas(1.e-6),nb_boucle_maxi(100)
{};
// Constructeur de copie
Prandtl_Reuss2D_D::Prandtl_Reuss2D_D (const Prandtl_Reuss2D_D& loi) :
Loi_comp_abstraite(loi),E(loi.E),nu(loi.nu)
,f_ecrouissage(Courbe1D::New_Courbe1D(*(loi.f_ecrouissage)))
,tolerance_plas(loi.tolerance_plas),nb_boucle_maxi(loi.nb_boucle_maxi)
{ };
Prandtl_Reuss2D_D::~Prandtl_Reuss2D_D ()
// Destructeur
{ if (f_ecrouissage != NULL)
if (f_ecrouissage->NomCourbe() == "_") delete f_ecrouissage;
};
// Lecture des donnees de la classe sur fichier
void Prandtl_Reuss2D_D::LectureDonneesParticulieres (UtilLecture * entreePrinc,LesCourbes1D& lesCourbes1D
,LesFonctions_nD& lesFonctionsnD)
{ // lecture du module d'young et du coefficient de poisson
if(strstr(entreePrinc->tablcar,"E=")== NULL)
{ cout << "\n erreur en lecture du module d'young "
<< " on attendait la chaine : E= ";
cout << "\n Prandtl_Reuss2D_D::LectureDonneesParticulieres "
<< "(UtilLecture * entreePrinc) " << endl ;
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
}
if(strstr(entreePrinc->tablcar,"nu=")== NULL)
{ cout << "\n erreur en lecture du coefficient de poisson "
<< " on attendait la chaine : nu= ";
cout << "\n Prandtl_Reuss2D_D::LectureDonneesParticulieres "
<< "(UtilLecture * entreePrinc) " << endl ;
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
}
string nom,toto;
*(entreePrinc->entree) >> nom >> E >> nom >> nu;
// lecture de la loi d'écrouissage
entreePrinc->NouvelleDonnee(); // lecture d'une nouvelle ligne
if(strstr(entreePrinc->tablcar,"loi_ecrouissage")== NULL)
{ cout << "\n erreur en lecture de la loi d'écrouissage "
<< " on attendait la chaine : loi_ecrouissage";
cout << "\n Prandtl_Reuss2D_D::LectureDonneesParticulieres "
<< "(UtilLecture * entreePrinc) " << endl ;
throw (UtilLecture::ErrNouvelleDonnee(-1));
Sortie(1);
}
*(entreePrinc->entree) >> toto >> nom;
// on regarde si la courbe existe, si oui on récupère la référence
if (lesCourbes1D.Existe(nom))
{ f_ecrouissage = lesCourbes1D.Trouve(nom);
}
else
{ // sinon il faut la lire maintenant
string non_courbe("_");
f_ecrouissage = Courbe1D::New_Courbe1D(non_courbe,Id_Nom_Courbe1D (nom.c_str()));
// lecture de la courbe
f_ecrouissage->LectDonnParticulieres_courbes (non_courbe,entreePrinc);
}
// appel au niveau de la classe mère
Loi_comp_abstraite::Lecture_type_deformation_et_niveau_commentaire
(*entreePrinc,lesFonctionsnD);
};
// affichage de la loi
void Prandtl_Reuss2D_D::Affiche() const
{ cout << " \n loi_de_comportement PRANDTL_REUSS "
<< " \n E= " << E << " nu= " << nu ;
cout << " \n loi_ecrouissage " ;
f_ecrouissage->Affiche();
// appel de la classe mère
Loi_comp_abstraite::Affiche_don_classe_abstraite();
};
// affichage et definition interactive des commandes particulières à chaques lois
void Prandtl_Reuss2D_D::Info_commande_LoisDeComp(UtilLecture& entreePrinc)
{ ofstream & sort = *(entreePrinc.Commande_pointInfo()); // pour simplifier
sort << "\n# ....... loi_de_comportement PRANDTL_REUSS 2D en deformation plane ........"
<< "\n# module d'young : coefficient de poisson "
<< "\n E= " << setprecision(6) << E << " nu= " << setprecision(6) << nu
<< "\n# on doit maintenant definir le nom d'une courbe 1D deja defini au debut du fichier .info,"
<< "\n# qui donnera la courbe d'ecrouissabe sigmabarre = f(epsilonbarre): par exemple "
<< "\n# fonction1 ou alors a la suite definir directement la courbe (cf. def de courbe) "
<< "\n# sans un nom de reference " << endl;
// appel de la classe mère
Loi_comp_abstraite::Info_commande_don_LoisDeComp(entreePrinc);
};
// test si la loi est complete
int Prandtl_Reuss2D_D::TestComplet()
{ int ret = LoiAbstraiteGeneral::TestComplet();
if (E == 0.)
{ cout << " \n le module d'young n'est pas défini pour la loi " << Nom_comp(id_comp)
<< '\n';
ret = 0;
}
if (nu == -2.)
{ cout << " \n le coefficient de poisson n'est pas défini pour la loi " << Nom_comp(id_comp)
<< '\n';
ret = 0;
}
if ( f_ecrouissage == NULL)
{ cout << " \n la fonction d'écrouissage n'est pas défini pour la loi "
<< Nom_comp(id_comp)
<< '\n';
ret = 0;
}
return ret;
};
// ========== codage des METHODES VIRTUELLES protegees:================
// calcul des contraintes
void Prandtl_Reuss2D_D::Calcul_SigmaHH (TenseurHH & sigHH_t,TenseurBB& ,DdlElement & tab_ddl,
TenseurBB & gijBB_t,TenseurHH & gijHH_t,BaseB& giB,BaseH& gi_H,TenseurBB& epsBB_,
TenseurBB& delta_epsBB_,TenseurBB & gijBB_,
TenseurHH & gijHH_,Tableau <TenseurBB *>& d_gijBB_,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 (epsBB_.Dimension() != 2)
{ cout << "\nErreur : la dimension devrait etre 2 !\n";
cout << " Prandtl_Reuss2D_D::Calcul_SigmaHH\n";
Sortie(1);
};
if (tab_ddl.NbDdl() != d_gijBB_.Taille())
{ cout << "\nErreur : le nb de ddl est != de la taille de d_gijBB_ !\n";
cout << " Prandtl_Reuss2D_D::Calcul_SigmaHH\n";
Sortie(1);
};
#endif
const Tenseur2BB & epsBB = *((Tenseur2BB*) &epsBB_); // passage en dim 2
const Tenseur2BB & delta_epsBB = *((Tenseur2BB*) &delta_epsBB_); // passage en dim 2
const Tenseur2HH & gijHH = *((Tenseur2HH*) &gijHH_); // " " " "
const Tenseur2BB & gijBB = *((Tenseur2BB*) &gijBB_); // " " " "
Tenseur2HH & sigHH = *((Tenseur2HH*) &sigHH_); // " " " "
Tenseur2HH & sigHH_i = *((Tenseur2HH*) &sigHH_t); // " " " "
SaveResulPrandtl_Reuss2D_D & save_resul = *((SaveResulPrandtl_Reuss2D_D*) saveResul);
//================== arrêt ici du travail ========================
// le tenseur des contraintes initiale en mixte
Tenseur2BH sigBH_i = gijBB * sigHH_i ;
// tout d'abord on considère que l'incrément est purement élastique et on
// regarde si la surface de plasticité n'a pas augmenté
// a) calcul du tenseur élastique résiduel
Tenseur2BB eps_elas_nBB = epsBB - save_resul.def_plasBB_t; //def car utile pour la plasticité
Tenseur2BH eps_elasBH = eps_elas_nBB * gijHH; // deformation en mixte
// calcul des coefficients
double coef1 = (E*nu)/((1.-2.*nu)*(1.+nu));
double coef2 = E/(1.+ nu);
// calcul du deviateur des deformations élastiques
double Ieps = eps_elasBH.Trace();
Tenseur2BH sigBH = (Ieps * coef1) * IdBH3 + coef2 * eps_elasBH ; // contrainte en mixte
// b) calcul de la nouvelle contrainte équivalente
double Isig = sigBH.Trace();
Tenseur2BH S_BH = sigBH - Isig * IdBH3/3.; // le déviateur
double sig_equi=sqrt(3./2. * S_BH && S_BH);
// c) test et orientation ad hoc
if (f_ecrouissage->Valeur(save_resul.epsilon_barre_t) >= sig_equi)
// cas ou l'élasticité est confirmée
{ // passage en 2fois contravariants
sigHH = gijHH * sigBH;
}
else
// cas ou l'on est en élastoplasticité
{ // la procédure de calcul est de type newton
const int nbddl_def = 6; // le nombre de ddl de déformation
Tenseur2BB deps_plasBB; // def de l'incrément de la déformation plastique en BB
Tenseur2BH deps_plasBH; // def de l'incrément de la déformation plastique en BH
// delta de lambda d'une itération à l'autre en BB
double delta_lambda; // initialisation à zéro
double lambda = 0.; // le lambda résultant
double res_plas; // résidu de l'équation sur la surface plastique
Tenseur2BB S_BB; // déviateur en BB
Tenseur2BB sigBB; // contrainte en BB
double & epsilon_barre = save_resul.epsilon_barre; // pour simplifier l'écriture
double & epsilon_barre_t = save_resul.epsilon_barre_t; // pour simplifier l'écriture
Tenseur2BB& def_plasBB_t = save_resul.def_plasBB_t; // """
Tenseur2BB& def_plasBB = save_resul.def_plasBB; // """
const double deux_tiers = 2./3.;
const double un_tiers = 1./3.;
const double un_demi = 1./2.;
const double coef2_carre = coef2 * coef2;
const double coef2_cube = coef2 * coef2_carre;
const double racine_deux_tiers = sqrt(deux_tiers);
// le déviateur de la déformation élastique
Tenseur2BH eps_elas_barre_BH = eps_elasBH - (un_tiers * Ieps) *IdBH3;
// puis en deux fois covariant
Tenseur2BB eps_elas_barre_BB = eps_elas_barre_BH * gijBB;
epsilon_barre = epsilon_barre_t; // init
// def du sigma équivalent initial
double sig_equi_i = f_ecrouissage->Valeur(epsilon_barre_t);
sigBH = sigBH_i; // init
Tableau2 <int> IJK = OrdreContrainte(nbddl_def); // pour la transformation 6 -> (i,j)
// ijk(n,1) correspond au premier indice i, et ijk(n,2) correspond au deuxième indice j
// acroissement de la déformation équivalente plastique
double delta_eps_equi = 0;
// constante c
double c_c = eps_elas_barre_BH && eps_elas_barre_BH;
// pour le calcul de la variation de f c-a-d le résidu
double Df_Dlambda ;
// variables qui sont utilisées après la boucle
double un_sur_1_plus_2_G_lambda
,un_sur_1_plus_2_G_lambda2,un_sur_1_plus_2_G_lambda3;
Courbe1D::ValDer valder;
int nb_iter = 1;
bool fin_plastique = false; // pour la fin de la boucle suivante
while ((!fin_plastique) && (nb_iter <= nb_boucle_maxi))
{
un_sur_1_plus_2_G_lambda = 1. / (1. + coef2*lambda);
un_sur_1_plus_2_G_lambda2 =
un_sur_1_plus_2_G_lambda * un_sur_1_plus_2_G_lambda;
un_sur_1_plus_2_G_lambda3 =
un_sur_1_plus_2_G_lambda * un_sur_1_plus_2_G_lambda2;
// nouveau déviateur de contrainte en mixte
S_BH = (coef2 * un_sur_1_plus_2_G_lambda) * eps_elas_barre_BH;
// en deux fois covariant
S_BB = S_BH * gijBB;
// calcul de l'incrément de déformation plastique
deps_plasBB = lambda * S_BB;
deps_plasBH = deps_plasBB * gijHH;// delta def plastique en BH
def_plasBB = def_plasBB_t + deps_plasBB; // deformation plastique
delta_eps_equi = sqrt(deux_tiers * (deps_plasBH && deps_plasBH));
epsilon_barre = epsilon_barre_t + // def plastique cumulée
delta_eps_equi;
// nouvelle valeur de la variation de sigma_barre
// et valeur de la tangente de la courbe d'écrouissage
valder = f_ecrouissage->Valeur_Et_derivee(epsilon_barre);
sig_equi = valder.valeur;
// calcul du résidu
res_plas = coef2_carre * c_c * un_sur_1_plus_2_G_lambda2 * un_demi
- sig_equi * sig_equi * un_tiers;
// test d'arrêt, pas pour la première itération car il faut
// le calcul de pas mal de grandeur pour le calcul de variation qui
// suit la boucle plastique
if (( Dabs(res_plas) < tolerance_plas) && (nb_iter!= 1))
{fin_plastique = true;
if (lambda < 0)
cout << "\n *** erreur : lambda plastique négatif "
<< "\n Prandtl_Reuss2D_D::Calcul_DsigmaHH_tdt (.. ";
break;
}
// cas où l'on doit faire une itération supplémentaire
// calcul de la variation du résidu par rapport à lamda
// 1 - calcul de la variation du premier terme
double D1_Dlambda = - un_sur_1_plus_2_G_lambda3 * c_c * coef2_cube;
// 2 - calcul de la variation de la déformation plastique cumulée
double der_eps_plas_lambda = racine_deux_tiers * coef2
* sqrt(c_c) * un_sur_1_plus_2_G_lambda2;
// 3 - calcul de la variation du second terme de f
double D2_Dlambda = - deux_tiers * sig_equi * valder.derivee * der_eps_plas_lambda;
// 4 - calcul de la variation de f
Df_Dlambda = D1_Dlambda + D2_Dlambda;
// calcul de l'incrément de lambda
delta_lambda = - res_plas / Df_Dlambda;
// nouvelle valeur de lambda
lambda += delta_lambda;
// incrémentation de la boucle
nb_iter++;
}
// sortie de la boucle on vérifie que la convergence est ok sinon message
if (!fin_plastique)
{ cout << "\n attention non convergence sur la plasticité " << endl;
}
// calcul du tenseur des contraintes
sigBH = S_BH + (un_tiers*Ieps * E/(1.-2.*nu)) * IdBH3;
// passage en 2 fois contravariants
sigHH = gijHH * sigBH;
// passage aussi en 2 fois covariants (utilisé pour les variations)
sigBB = sigBH * gijBB;
}
LibereTenseur();
};
// calcul des contraintes a t+dt et de ses variations
void Prandtl_Reuss2D_D::Calcul_DsigmaHH_tdt (TenseurHH & sigHH_t,TenseurBB& ,DdlElement & tab_ddl
,BaseB& ,TenseurBB & gijBB_t,TenseurHH & gijHH_t,
BaseB& giB_tdt,Tableau <BaseB> & d_giB_tdt,BaseH& giH_tdt,Tableau <BaseH> & d_giH_tdt,
TenseurBB & epsBB_tdt,Tableau <TenseurBB *>& d_epsBB,TenseurBB & delta_epsBB_,
TenseurBB & gijBB_tdt ,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 (epsBB_tdt.Dimension() != 3)
{ cout << "\nErreur : la dimension devrait etre 3 !\n";
cout << " Prandtl_Reuss2D_D::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 << " Prandtl_Reuss2D_D::Calcul_DsigmaHH_tdt\n";
Sortie(1);
};
#endif
const Tenseur2BB & epsBB = *((Tenseur2BB*) &epsBB_tdt); // passage en dim 3
const Tenseur2BB & delta_epsBB = *((Tenseur2BB*) &delta_epsBB_); // passage en dim 3
const Tenseur2HH & gijHH = *((Tenseur2HH*) &gijHH_tdt); // " " " "
const Tenseur2BB & gijBB = *((Tenseur2BB*) &gijBB_tdt); // " " " "
Tenseur2HH & sigHH = *((Tenseur2HH*) &sigHH_tdt); // " " " "
Tenseur2HH & sigHH_i = *((Tenseur2HH*) &sigHH_t); // " " " "
SaveResulPrandtl_Reuss2D_D & save_resul = *((SaveResulPrandtl_Reuss2D_D*) saveResul);
// pour la variation du tenseur des contraintes
int nbddl = d_gijBB_tdt.Taille();
// le tenseur des contraintes initiale en mixte
Tenseur2BH sigBH_i = gijBB * sigHH_i ;
// tout d'abord on considère que l'incrément est purement élastique et on
// regarde si la surface de plasticité n'a pas augmenté
// a) calcul du tenseur élastique résiduel
Tenseur2BB eps_elas_nBB = epsBB - save_resul.def_plasBB_t; //def car utile pour la plasticité
Tenseur2BH eps_elasBH = eps_elas_nBB * gijHH; // deformation en mixte
// calcul des coefficients
double coef1 = (E*nu)/((1.-2.*nu)*(1.+nu));
double coef2 = E/(1.+ nu);
// calcul du deviateur des deformations élastiques
double Ieps = eps_elasBH.Trace();
Tenseur2BH sigBH = (Ieps * coef1) * IdBH3 + coef2 * eps_elasBH ; // contrainte en mixte
// b) calcul de la nouvelle contrainte équivalente
double Isig = sigBH.Trace();
Tenseur2BH S_BH = sigBH - Isig * IdBH3/3.; // le déviateur
double sig_equi=sqrt(3./2. * S_BH && S_BH);
// c) test et orientation ad hoc
if (f_ecrouissage->Valeur(save_resul.epsilon_barre_t) >= sig_equi)
// cas ou l'élasticité est confirmée
{ // passage en 2foi contravariants
sigHH = gijHH * sigBH;
// la variation de la contrainte c'est en fait la variation du delta contrainte
for (int i = 1; i<= nbddl; i++)
{ // on fait de faire uniquement une égalité d'adresse et de ne pas utiliser
// le constructeur d'ou la profusion d'* et de ()
Tenseur2HH & dsigHH = *((Tenseur2HH*) (d_sigHH(i))); // passage en dim 3
const Tenseur2BB & d_gijBB = *((Tenseur2BB*)(d_gijBB_tdt(i))); // passage en dim 3
const Tenseur2HH & dgijHH = *((Tenseur2HH*)(d_gijHH_tdt(i))) ; // pour simplifier l'ecriture
// pour chacun des ddl on calcul les tenseurs derivees
Tenseur2BB d_eps_elasBB = 0.5 * d_gijBB;
Tenseur2BH d_eps_elasBH = eps_elas_nBB * dgijHH + d_eps_elasBB * gijHH ;
double dIeps = d_eps_elasBH.Trace();
dsigHH = dgijHH * sigBH + gijHH *
((dIeps * coef1) * IdBH3 + coef2 * d_eps_elasBH);
}
}
else
// cas ou l'on est en élastoplasticité
{ // la procédure de calcul est de type newton
const int nbddl_def = 6; // le nombre de ddl de déformation
Tenseur2BB deps_plasBB; // def de l'incrément de la déformation plastique en BB
Tenseur2BH deps_plasBH; // def de l'incrément de la déformation plastique en BH
// delta de lambda d'une itération à l'autre en BB
double delta_lambda; // initialisation à zéro
double lambda = 0.; // le lambda résultant
double res_plas; // résidu de l'équation sur la surface plastique
Tenseur2BB S_BB; // déviateur en BB
Tenseur2BB sigBB; // contrainte en BB
double & epsilon_barre = save_resul.epsilon_barre; // pour simplifier l'écriture
double & epsilon_barre_t = save_resul.epsilon_barre_t; // pour simplifier l'écriture
Tenseur2BB& def_plasBB_t = save_resul.def_plasBB_t; // """
Tenseur2BB& def_plasBB = save_resul.def_plasBB; // """
const double deux_tiers = 2./3.;
const double un_tiers = 1./3.;
const double un_demi = 1./2.;
const double coef2_carre = coef2 * coef2;
const double coef2_cube = coef2 * coef2_carre;
const double racine_deux_tiers = sqrt(deux_tiers);
// le déviateur de la déformation élastique
Tenseur2BH eps_elas_barre_BH = eps_elasBH - (un_tiers * Ieps) *IdBH3;
// puis en deux fois covariant
Tenseur2BB eps_elas_barre_BB = eps_elas_barre_BH * gijBB;
epsilon_barre = epsilon_barre_t; // init
// def du sigma équivalent initial
double sig_equi_i = f_ecrouissage->Valeur(epsilon_barre_t);
sigBH = sigBH_i; // init
Tableau2 <int> IJK = OrdreContrainte(nbddl_def); // pour la transformation 6 -> (i,j)
// ijk(n,1) correspond au premier indice i, et ijk(n,2) correspond au deuxième indice j
// acroissement de la déformation équivalente plastique
double delta_eps_equi = 0;
// constante c
double c_c = eps_elas_barre_BH && eps_elas_barre_BH;
// pour le calcul de la variation de f c-a-d le résidu
double Df_Dlambda ;
// variables qui sont utilisées après la boucle
double un_sur_1_plus_2_G_lambda
,un_sur_1_plus_2_G_lambda2,un_sur_1_plus_2_G_lambda3;
Courbe1D::ValDer valder;
int nb_iter = 1;
bool fin_plastique = false; // pour la fin de la boucle suivante
while ((!fin_plastique) && (nb_iter <= nb_boucle_maxi))
{
un_sur_1_plus_2_G_lambda = 1. / (1. + coef2*lambda);
un_sur_1_plus_2_G_lambda2 =
un_sur_1_plus_2_G_lambda * un_sur_1_plus_2_G_lambda;
un_sur_1_plus_2_G_lambda3 =
un_sur_1_plus_2_G_lambda * un_sur_1_plus_2_G_lambda2;
// nouveau déviateur de contrainte en mixte
S_BH = (coef2 * un_sur_1_plus_2_G_lambda) * eps_elas_barre_BH;
// en deux fois covariant
S_BB = S_BH * gijBB;
// calcul de l'incrément de déformation plastique
deps_plasBB = lambda * S_BB;
deps_plasBH = deps_plasBB * gijHH;// delta def plastique en BH
def_plasBB = def_plasBB_t + deps_plasBB; // deformation plastique
delta_eps_equi = sqrt(deux_tiers * (deps_plasBH && deps_plasBH));
epsilon_barre = epsilon_barre_t + // def plastique cumulée
delta_eps_equi;
// nouvelle valeur de la variation de sigma_barre
// et valeur de la tangente de la courbe d'écrouissage
valder = f_ecrouissage->Valeur_Et_derivee(epsilon_barre);
sig_equi = valder.valeur;
// calcul du résidu
res_plas = coef2_carre * c_c * un_sur_1_plus_2_G_lambda2 * un_demi
- sig_equi * sig_equi * un_tiers;
// if (res_plas < 0)
// cout << "\n *** erreur : res_plas négatif "
// << "\n Prandtl_Reuss2D_D::Calcul_DsigmaHH_tdt (.. " << endl;
// test d'arrêt, pas pour la première itération car il faut
// le calcul de pas mal de grandeur pour le calcul de variation qui
// suit la boucle plastique
if (( Dabs(res_plas) < tolerance_plas) && (nb_iter!= 1))
{fin_plastique = true;
if (lambda < 0)
cout << "\n *** erreur : lambda plastique négatif "
<< "\n Prandtl_Reuss2D_D::Calcul_DsigmaHH_tdt (.. ";
break;
}
// cas où l'on doit faire une itération supplémentaire
// calcul de la variation du résidu par rapport à lamda
// 1 - calcul de la variation du premier terme
double D1_Dlambda = - un_sur_1_plus_2_G_lambda3 * c_c * coef2_cube;
// 2 - calcul de la variation de la déformation plastique cumulée
double der_eps_plas_lambda = racine_deux_tiers * coef2
* sqrt(c_c) * un_sur_1_plus_2_G_lambda2;
// 3 - calcul de la variation du second terme de f
double D2_Dlambda = - deux_tiers * sig_equi * valder.derivee * der_eps_plas_lambda;
// 4 - calcul de la variation de f
Df_Dlambda = D1_Dlambda + D2_Dlambda;
// calcul de l'incrément de lambda
delta_lambda = - res_plas / Df_Dlambda;
// nouvelle valeur de lambda
lambda += delta_lambda;
// incrémentation de la boucle
nb_iter++;
}
// sortie de la boucle on vérifie que la convergence est ok sinon message
if (!fin_plastique)
{ cout << "\n attention non convergence sur la plasticité " << endl;
}
// calcul du tenseur des contraintes
sigBH = S_BH + (un_tiers*Ieps * E/(1.-2.*nu)) * IdBH3;
// passage en 2 fois contravariants
sigHH = gijHH * sigBH;
// passage aussi en 2 fois covariants (utilisé pour les variations)
sigBB = sigBH * gijBB;
// maintenant on s'occupe de la raideur tangente
// 1- calcul de la variation de f par rapport à epsilonBarre_ij_BB
// à lambda constant
Tenseur2HH Df_DepsBij_HH(eps_elas_barre_BB(1,1),eps_elas_barre_BB(2,2),
eps_elas_barre_BB(2,1));
Df_DepsBij_HH *= (coef2_carre * un_sur_1_plus_2_G_lambda2
- coef2 * deux_tiers * racine_deux_tiers
* un_sur_1_plus_2_G_lambda * sig_equi * lambda / sqrt(c_c)
* valder.derivee);
// 2- calcul de la variation de lambda par rapport à epsilonBarre_ij_BB
Tenseur2HH Dlambda_DepsBij_HH = (-1./Df_Dlambda) * Df_DepsBij_HH;
// 3- calcul de la variation de S_BB par rapport à epsilonBarre_ij_BB
// normalement c'est un tenseur du 4ième ordre mais que l'on remplace par un
// tableau de tenseur du second ordre de la variance de S, ici BB
Tableau <Tenseur2BB> D_S_BB_DepsBijBB(nbddl_def);
double a1 = - coef2_carre * un_sur_1_plus_2_G_lambda2;
double a2 = coef2 * un_sur_1_plus_2_G_lambda;
// le premier terme
for (int e4=1; e4<= nbddl_def; e4++)
for (int k=1; k<= 3;k++)
for (int m=1;m<=k;m++)
D_S_BB_DepsBijBB(e4).Coor(k,m) =
a1 * eps_elas_barre_BB(k,m)* Dlambda_DepsBij_HH(IJK(e4,1),IJK(e4,2));
// le second terme qui aurait du être : + a2 * IdBH3(k,IJK(e4,1)) * IdBH3(m,IJK(e4,2));
// mais cela ne fonctionne pas car e4 = 4 5 6 correspond à 1,2 2,3 et 1,3 mais il manque
// les termes symétriques
D_S_BB_DepsBijBB(1).Coor(1,1) += a2; D_S_BB_DepsBijBB(4).Coor(1,2) += a2;
D_S_BB_DepsBijBB(2).Coor(2,2) += a2; D_S_BB_DepsBijBB(5).Coor(2,3) += a2;
D_S_BB_DepsBijBB(3).Coor(3,3) += a2; D_S_BB_DepsBijBB(6).Coor(1,3) += a2;
// 4- calcul de la variation par rapport à epsilon_ij_BB
// Tableau <Tenseur2BB> D_S_BB_DepsijBB(nbddl_def) = D_S_BB_DepsBijBB(ff);
Tableau <Tenseur2BB> D_S_BB_DepsijBB(D_S_BB_DepsBijBB);
for (int e=1; e<= nbddl_def; e++)
{for (int ff=1; ff<= 3; ff++) // tout d'abord les termes ii
D_S_BB_DepsijBB(e) -= D_S_BB_DepsBijBB(ff) *
( un_tiers * gijHH(IJK(e,1),IJK(e,2)) * gijBB(IJK(ff,1),IJK(ff,2)));
for (int ff=4; ff<= nbddl_def; ff++) // puis les termes ij (i différent de j)
// il y en a deux à chaque fois
D_S_BB_DepsijBB(e) -= D_S_BB_DepsBijBB(ff) *
( deux_tiers * gijHH(IJK(e,1),IJK(e,2)) * gijBB(IJK(ff,1),IJK(ff,2)));
}
// 5- calcul de la variation de sigmaBB par rapport à epsilon_ij_BB
Tableau <Tenseur2BB> D_sig_BB_DepsijBB(nbddl_def);
double a3 = E/(3.*(1.-2.*nu));
for (int e=1; e<= nbddl_def; e++)
D_sig_BB_DepsijBB(e) = D_S_BB_DepsijBB(e)
+ (a3 * gijHH(IJK(e,1),IJK(e,2))) * gijBB;
// 6- calcul de la variation de sigmaHH par rapport au ddl
for (int iddl = 1; iddl<= nbddl; iddl++)
{ // on fait uniquement une égalité d'adresse pour ne pas utiliser
// le constructeur d'ou la profusion d'* et de ()
Tenseur2HH & dsigHH = *((Tenseur2HH*) (d_sigHH(iddl))); // passage en dim 3
const Tenseur2BB & d_gijBB = *((Tenseur2BB*)(d_gijBB_tdt(iddl))); // passage en dim 3
const Tenseur2HH & dgijHH = *((Tenseur2HH*)(d_gijHH_tdt(iddl))) ; // pour simplifier l'ecriture
// pour chacun des ddl on calcul les tenseurs derivees
// def de la variation de la déformation
Tenseur2BB d_eps_elasBB = 0.5 * d_gijBB;
// 7- calcul de la variation de sigmaBB par rapport au ddl
Tenseur2BB dsigBB;
for (int e1=1;e1<= 3; e1++) // tout d'abord les termes ii
dsigBB += D_sig_BB_DepsijBB(e1) * d_eps_elasBB(IJK(e1,1),IJK(e1,2));
for (int e1=4;e1<= nbddl_def; e1++) // puis les termes ij (i différent de j)
//qui doivent être doublée car il y en a deux, exemple 12 puis 21
dsigBB += 2. * D_sig_BB_DepsijBB(e1) * d_eps_elasBB(IJK(e1,1),IJK(e1,2));
// 8- calcul maintenant en deux fois contravariants
dsigHH = gijHH * (dsigBB * gijHH) + gijHH * (sigBB * dgijHH)
+ dgijHH * sigBH;
}
}
LibereTenseur();
};
//----- 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 Prandtl_Reuss2D_D::Lecture_base_info_loi(ifstream& ent,const int cas,LesReferences& lesRef,LesCourbes1D& lesCourbes1D
,LesFonctions_nD& lesFonctionsnD)
{ string toto;
if (cas == 1)
{ ent >> toto >> E >> toto >> nu;
// la courbe d'écrouissage
ent >> toto;
if (toto != "f_ecrouissage")
{ cout << "\n erreur en lecture de la fonction d'ecrouissage, on attendait f_ecrouissage et on a lue " << toto
<< "\n Prandtl_Reuss2D_D::Lecture_base_info_loi(...";
Sortie(1);
};
f_ecrouissage = lesCourbes1D.Lecture_pour_base_info(ent,cas,f_ecrouissage);
// lecture des tol
ent >> toto >> tolerance_plas
>> toto >> nb_boucle_maxi ;
}
// appel de la classe mère
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 Prandtl_Reuss2D_D::Ecriture_base_info_loi(ofstream& sort,const int cas)
{ if (cas == 1)
{ sort << " module_d'young " << E << " nu " << nu ;
// la courbe d'écrouissage
sort << " \n f_ecrouissage ";
LesCourbes1D::Ecriture_pour_base_info(sort,cas,f_ecrouissage);
sort << "\n tolerance_algorithme " << tolerance_plas
<< " nb_boucle_maxi " << nb_boucle_maxi << " ";
}
// appel de la classe mère
Loi_comp_abstraite::Ecriture_don_base_info(sort,cas);
};
Reference in a new issue Copy permalink