// 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) . // // 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 . // // For more information, please consult: . //#include "Debug.h" # include using namespace std; //introduces namespace std #include #include #include "Sortie.h" #include "ConstMath.h" #include "Prandtl_Reuss1D.h" // ========== fonctions pour la classe de sauvegarde des résultats ========= // affectation Loi_comp_abstraite::SaveResul & Prandtl_Reuss1D::SaveResulPrandtl_Reuss1D::operator = ( const Loi_comp_abstraite::SaveResul & a) { Prandtl_Reuss1D::SaveResulPrandtl_Reuss1D& sav = *((Prandtl_Reuss1D::SaveResulPrandtl_Reuss1D*) &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_Reuss1D::SaveResulPrandtl_Reuss1D::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_Reuss1D::SaveResulPrandtl_Reuss1D::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_Reuss1D::SaveResulPrandtl_Reuss1D::Affiche() const { cout << "\n SaveResulPrandtl_Reuss1D: " ; 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_Reuss1D::Prandtl_Reuss1D () : // Constructeur par defaut Loi_comp_abstraite(PRANDTL_REUSS1D,CAT_MECANIQUE,1),E(0.),nu(-2.) ,f_ecrouissage(NULL) ,tolerance_plas(1.e-6),nb_boucle_maxi(100),nb_sous_increment(4) {}; // Constructeur de copie Prandtl_Reuss1D::Prandtl_Reuss1D (const Prandtl_Reuss1D& 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) ,nb_sous_increment(loi.nb_sous_increment) { }; Prandtl_Reuss1D::~Prandtl_Reuss1D () // Destructeur { if (f_ecrouissage != NULL) if (f_ecrouissage->NomCourbe() == "_") delete f_ecrouissage; }; // Lecture des donnees de la classe sur fichier void Prandtl_Reuss1D::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_Reuss1D::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_Reuss1D::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_Reuss1D::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_Reuss1D::Affiche() const { cout << " \n loi_de_comportement PRANDTL_REUSS1D " << " \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_Reuss1D::Info_commande_LoisDeComp(UtilLecture& entreePrinc) { ofstream & sort = *(entreePrinc.Commande_pointInfo()); // pour simplifier sort << "\n# ....... loi_de_comportement PRANDTL_REUSS1D ........" << "\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_Reuss1D::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_Reuss1D::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 & 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() != 3) { cout << "\nErreur : la dimension devrait etre 3 !\n"; cout << " Prandtl_Reuss1D::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_Reuss1D::Calcul_SigmaHH\n"; Sortie(1); }; #endif /* const Tenseur1BB & epsBB = *((Tenseur1BB*) &epsBB_); // passage en dim 3 const Tenseur1BB & delta_epsBB = *((Tenseur1BB*) &delta_epsBB_); // passage en dim 3 const Tenseur1HH & gijHH = *((Tenseur1HH*) &gijHH_); // " " " " const Tenseur1BB & gijBB = *((Tenseur1BB*) &gijBB_); // " " " " Tenseur1HH & sigHH = *((Tenseur1HH*) &sigHH_); // " " " " Tenseur1HH & sigHH_i = *((Tenseur1HH*) &sigHH_t); // " " " " SaveResulPrandtl_Reuss1D & save_resul = *((SaveResulPrandtl_Reuss1D*) saveResul); // le tenseur des contraintes initiale en mixte Tenseur1BH 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 Tenseur1BB eps_elas_nBB = epsBB - save_resul.def_plasBB_t; //def car utile pour la plasticité Tenseur1BH 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(); Tenseur1BH sigBH = (Ieps * coef1) * IdBH3 + coef2 * eps_elasBH ; // contrainte en mixte // b) calcul de la nouvelle contrainte équivalente double Isig = sigBH.Trace(); Tenseur1BH 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 Tenseur1BB deps_plasBB; // def de l'incrément de la déformation plastique en BB Tenseur1BH 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 Tenseur1BB S_BB; // déviateur en BB Tenseur1BB 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 Tenseur1BB& def_plasBB_t = save_resul.def_plasBB_t; // """ Tenseur1BB& 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 Tenseur1BH eps_elas_barre_BH = eps_elasBH - (un_tiers * Ieps) *IdBH3; // puis en deux fois covariant Tenseur1BB 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 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 = 4.*lambda*omega *sqrt(c_c); 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 = 3.*c_c*omega2 - 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_Reuss1D::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 = -6.*c_c*omega2; // 2 - calcul de la variation de la déformation plastique cumulée double der_eps_plas_lambda = 2.*omega*(1.-2.*lambda*omega)* sqrt(c_c); // 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; } // 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_Reuss1D::Calcul_DsigmaHH_tdt (TenseurHH & sigHH_t,TenseurBB& ,DdlElement & tab_ddl ,BaseB& ,TenseurBB & gijBB_t,TenseurHH & gijHH_t, BaseB& giB_tdt,Tableau & d_giB_tdt,BaseH& giH_tdt,Tableau & d_giH_tdt, TenseurBB & epsBB_tdt,Tableau & d_epsBB,TenseurBB & delta_epsBB_, TenseurBB & gijBB_tdt ,TenseurHH & gijHH_tdt, Tableau & d_gijBB_tdt, Tableau & d_gijHH_tdt,double& ,double& , Vecteur& ,TenseurHH& sigHH_tdt,Tableau & 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() != 1) { cout << "\nErreur : la dimension devrait etre 1 !\n"; cout << " Prandtl_Reuss1D::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_Reuss1D::Calcul_DsigmaHH_tdt\n"; Sortie(1); }; #endif const Tenseur1BB & epsBB = *((Tenseur1BB*) &epsBB_tdt); // passage en dim 1 const Tenseur1BB & delta_epsBB = *((Tenseur1BB*) &delta_epsBB_); // passage en dim 1 const Tenseur1HH & gijHH = *((Tenseur1HH*) &gijHH_tdt); // " " " " const Tenseur1BB & gijBB = *((Tenseur1BB*) &gijBB_tdt); // " " " " Tenseur1HH & sigHH = *((Tenseur1HH*) &sigHH_tdt); // " " " " Tenseur1HH & sigHH_i = *((Tenseur1HH*) &sigHH_t); // " " " " SaveResulPrandtl_Reuss1D & save_resul = *((SaveResulPrandtl_Reuss1D*) saveResul); // pour la variation du tenseur des contraintes int nbddl = d_gijBB_tdt.Taille(); // le tenseur des contraintes initiale en mixte Tenseur1BH 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 Tenseur1BB eps_elasBB = epsBB - save_resul.def_plasBB_t; // deformation en mixte Tenseur1BH eps_elasBH = eps_elasBB * gijHH; // deformation en mixte Tenseur1BH sigBH = E * eps_elasBH; sigHH = gijHH * sigBH; // b) test et orientation ad hoc if (f_ecrouissage->Valeur(save_resul.epsilon_barre_t) >= Dabs(sigBH(1,1))) // cas ou l'élasticité est confirmée { // la variation de la contrainte c'est en fait la variation du delta contrainte 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 () Tenseur1HH & dsigHH = *((Tenseur1HH*) (d_sigHH(i))); // passage en dim 1 const Tenseur1BB & d_gijBB = *((Tenseur1BB*)(d_gijBB_tdt(i))); // passage en dim 1 const Tenseur1HH & dgijHH = *((Tenseur1HH*)(d_gijHH_tdt(i))) ; // pour simplifier l'ecriture // pour chacun des ddl on calcul les tenseurs derivees dsigHH = dgijHH * sigBH + gijHH * E * (0.5 * d_gijBB * gijHH + epsBB * dgijHH ); } } else // cas ou l'on est en élastoplasticité { // la procédure de calcul est de type newton // calcul de coefficients double K = E/(1.-2.*nu); double deux_G = E/(1.+ nu); const int nbddl_def = 1; // le nombre de ddl de déformation en mixte Tenseur1BH 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 double & epsilon_barre = save_resul.epsilon_barre; // pour simplifier l'écriture double & epsilon_barre_t = save_resul.epsilon_barre_t; // pour simplifier l'écriture Tenseur1BB& def_plasBB_t = save_resul.def_plasBB_t; // """ Tenseur1BB& 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 deux_G_carre = deux_G * deux_G; const double deux_G_cube = deux_G * deux_G_carre; const double racine_deux_tiers = sqrt(deux_tiers); 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 // acroissement de la déformation équivalente plastique double delta_eps_equi = 0; // constante c double c_c = eps_elasBH && eps_elasBH; // 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; double alphaa,alphaa2; double omega,omega2,omega3; double sig_equi; Courbe1D::ValDer valder; // ==== algorithme de recherche du multiplicateur plastique, // dans le cas où un passage n'est pas convergent on décompose l'incrément // en sous incréments, ceci jusqu'à convergence double coeff_increment = 1.; // init, a priori un seul incrément bool convergence_lambda = true; // init int nb_sous_inc = 1; bool fin_plastique = false; // while ((!fin_plastique) && (nb_sous_inc <= nb_sous_increment)) {int nb_iter = 1; bool fin_inc_plastique = false; // pour la fin de la boucle suivante while ((!fin_inc_plastique) && (nb_iter <= nb_boucle_maxi)) { un_sur_1_plus_2_G_lambda = 1. / (1. + deux_G*lambda); un_sur_1_plus_2_G_lambda2 = un_sur_1_plus_2_G_lambda * un_sur_1_plus_2_G_lambda; alphaa = un_sur_1_plus_2_G_lambda*deux_G; alphaa2 = un_sur_1_plus_2_G_lambda2 * deux_G_carre; omega = alphaa*K/(alphaa+2.*K); omega2 = omega*omega; omega3 = omega*omega2; // a moins que lambda soit très grand on considère qu'omega est positif // on ne teste pas le signe de lambda dans la boucle seulement en sortie // ce qui permet de garder le caratère le plus quadratique et continu possible delta_eps_equi = 2.*lambda*omega * Dabs(eps_elasBH(1,1)); 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 = 3.*c_c*omega2 - 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_inc_plastique = true; if (lambda < 0) cout << "\n *** erreur : lambda plastique négatif " << "\n Prandtl_Reuss1D::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 = -12.*c_c*omega3; // 2 - calcul de la variation de la déformation plastique cumulée double der_eps_plas_lambda = 2.*omega*(1.-2.*lambda*omega)* sqrt(c_c); //eps_elasBH(1,1); // 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_inc_plastique) { // pas de convergence, on diminue la sous_incrémentation si c'est possible cout << "\n attention non convergence sur la plasticité " << endl; } // { cout << "\n attention non convergence sur la plasticité " << endl; // } } // fin de la recherche de lambda // nouvelle valeur de la contrainte en mixte sigBH = 3. * omega * eps_elasBH; // calcul de l'incrément de déformation plastique deps_plasBH = deux_tiers *lambda * sigBH;// delta def plastique en BH // passage en 2 fois contravariants sigHH = gijHH * sigBH; // mise à jour de la déformation plastique def_plasBB = def_plasBB_t + deps_plasBH * gijBB ; // maintenant on s'occupe de la raideur tangente // 0- calcul de la variation de c par rapport à epsilonBarre_ij_BH // à lambda constant Tenseur1BH DC_Depsij_BH = 2.*eps_elasBH ; // ========= petite vérif ========= Tenseur3BH S_BHverif;S_BHverif.Coor(1,1)=2.; S_BHverif.Coor(2,2) = -1.;S_BHverif.Coor(3,3) = -1.; S_BHverif *= sigBH(1,1)/3.; double f_verif = un_demi*(S_BHverif && S_BHverif)- sig_equi * sig_equi * un_tiers; Tenseur3BH delta_eps_plasBHverif = lambda * S_BHverif; double delta_epsB_plasverif = sqrt(deux_tiers* (delta_eps_plasBHverif && delta_eps_plasBHverif)); //=========== pour le calcul de la dérivée de C de manière numérique ============= /* double peti= 1.E-10; valder = f_ecrouissage->Valeur_Et_derivee(epsilon_barre); sig_equi = valder.valeur; Tenseur1BH DC_Depsij_BH_ver; Tenseur1BH eps_elasBH_ver = eps_elasBH; eps_elasBH_ver(1,1) += peti; double c_c_ver = (eps_elasBH_ver)&&(eps_elasBH_ver); DC_Depsij_BH_ver(1,1) = (c_c_ver - c_c)/peti; */ //=========== fin du calcul de la dérivée ========================================= // 1- calcul de la variation de f par rapport à epsilon_ij_BH // à lambda constant Tenseur1BH Df_Depsij_BH (DC_Depsij_BH); Df_Depsij_BH *= (3.*omega2 - deux_tiers * sig_equi * lambda * omega / sqrt(c_c) * valder.derivee); //=========== pour le calcul de la dérivée de DF de manière numérique ============= /* eps_elasBH_ver = eps_elasBH; Tenseur1BH Df_Depsij_BH_ver; eps_elasBH_ver(1,1) += peti; c_c_ver = (eps_elasBH_ver)&&(eps_elasBH_ver); double delta_eps_equi_ver = 2.*lambda*omega * sqrt(c_c_ver); double epsilon_barre_ver = epsilon_barre_t + // def plastique cumulée delta_eps_equi_ver; // nouvelle valeur de la variation de sigma_barre // et valeur de la tangente de la courbe d'écrouissage Courbe1D::ValDer valder_ver = f_ecrouissage->Valeur_Et_derivee(epsilon_barre_ver); double sig_equi_ver = valder_ver.valeur; // calcul du résidu double res_plas_ver = 3.*c_c_ver*omega2 - sig_equi_ver * sig_equi_ver * un_tiers; Df_Depsij_BH_ver(1,1) = res_plas_ver/peti; */ //=========== fin du calcul de la dérivée ========================================= // 2- calcul de la variation de lambda par rapport à epsilon_ij_BH Tenseur1BH Dlambda_Depsij_BH = (-1./Df_Dlambda) * Df_Depsij_BH; //=========== pour le calcul de la dérivée de DF de manière numérique ============= /* Tenseur1BH Dlambda_Depsij_BH_ver; double Df_Dlambda_ver,delta_lambda_ver; double un_sur_1_plus_2_G_lambda_ver ,un_sur_1_plus_2_G_lambda2_ver; double lambda_ver = lambda; double alphaa_ver,alphaa2_ver; double omega_ver,omega2_ver,omega3_ver; eps_elasBH_ver = eps_elasBH; eps_elasBH_ver(1,1) += peti; c_c_ver = (eps_elasBH_ver)&&(eps_elasBH_ver); nb_iter = 1; 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_ver = 1. / (1. + deux_G*lambda_ver); un_sur_1_plus_2_G_lambda2_ver = un_sur_1_plus_2_G_lambda_ver * un_sur_1_plus_2_G_lambda_ver; alphaa_ver = un_sur_1_plus_2_G_lambda_ver*deux_G; alphaa2_ver = un_sur_1_plus_2_G_lambda2_ver * deux_G_carre; omega_ver = alphaa_ver*K/(alphaa_ver+2.*K); omega2_ver = omega_ver*omega_ver;omega3_ver = omega_ver*omega2_ver; delta_eps_equi_ver = 2.*lambda_ver*omega_ver * sqrt(c_c_ver); epsilon_barre_ver = epsilon_barre_t + // def plastique cumulée delta_eps_equi_ver; // nouvelle valeur de la variation de sigma_barre // et valeur de la tangente de la courbe d'écrouissage valder_ver = f_ecrouissage->Valeur_Et_derivee(epsilon_barre_ver); sig_equi_ver = valder_ver.valeur; // calcul du résidu res_plas_ver = 3.*c_c_ver*omega2_ver - sig_equi_ver * sig_equi_ver * 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_ver < 0) cout << "\n *** erreur : lambda plastique négatif " << "\n Prandtl_Reuss1D::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_ver = -12.*c_c*omega3_ver; // 2 - calcul de la variation de la déformation plastique cumulée double der_eps_plas_lambda_ver = 2.*omega_ver*(1.-2.*lambda_ver*omega_ver)* sqrt(c_c_ver); // 3 - calcul de la variation du second terme de f double D2_Dlambda_ver = - deux_tiers * sig_equi_ver * valder_ver.derivee * der_eps_plas_lambda_ver; // 4 - calcul de la variation de f Df_Dlambda_ver = D1_Dlambda_ver + D2_Dlambda_ver; // calcul de l'incrément de lambda delta_lambda_ver = - res_plas_ver / Df_Dlambda_ver; // nouvelle valeur de lambda lambda_ver += delta_lambda_ver; // 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 de la dérivée numérique de lambda Dlambda_Depsij_BH_ver(1,1) = (lambda_ver - lambda)/ peti; */ //=========== fin du calcul de la dérivée ========================================= // 3- calcul de la variation de sigma par rapport à epsilon_ij_BH Tenseur1BH D_sig_BH_DepsijBH_BH; D_sig_BH_DepsijBH_BH.Coor(1,1) = 3.*omega*(1.-2.*omega*eps_elasBH(1,1)*Dlambda_Depsij_BH(1,1)); //=========== pour le calcul de la dérivée de S_BH par rapport à epsBB de manière numérique ============= /* Tenseur1BH D_sig_BH_DepsijBH_BH_ver; Tenseur1BH sigBH_ver; eps_elasBH_ver = eps_elasBH; eps_elasBH_ver(1,1) += peti; c_c_ver = (eps_elasBH_ver)&&(eps_elasBH_ver); nb_iter = 1; 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_ver = 1. / (1. + deux_G*lambda_ver); un_sur_1_plus_2_G_lambda2_ver = un_sur_1_plus_2_G_lambda_ver * un_sur_1_plus_2_G_lambda_ver; alphaa_ver = un_sur_1_plus_2_G_lambda_ver*deux_G; alphaa2_ver = un_sur_1_plus_2_G_lambda2_ver * deux_G_carre; omega_ver = alphaa_ver*K/(alphaa_ver+2.*K); omega2_ver = omega_ver*omega_ver;omega3_ver = omega_ver*omega2_ver; // nouvelle valeur de la contrainte en mixte sigBH_ver = 3. * omega_ver * eps_elasBH_ver; delta_eps_equi_ver = 2.*lambda_ver*omega_ver * sqrt(c_c_ver); epsilon_barre_ver = epsilon_barre_t + // def plastique cumulée delta_eps_equi_ver; // nouvelle valeur de la variation de sigma_barre // et valeur de la tangente de la courbe d'écrouissage valder_ver = f_ecrouissage->Valeur_Et_derivee(epsilon_barre_ver); sig_equi_ver = valder_ver.valeur; // calcul du résidu res_plas_ver = 3.*c_c_ver*omega2_ver - sig_equi_ver * sig_equi_ver * 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_ver < 0) cout << "\n *** erreur : lambda plastique négatif " << "\n Prandtl_Reuss1D::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_ver = -12.*c_c*omega3_ver; // 2 - calcul de la variation de la déformation plastique cumulée double der_eps_plas_lambda_ver = 2.*omega_ver*(1.-2.*lambda_ver*omega_ver)* sqrt(c_c_ver); // 3 - calcul de la variation du second terme de f double D2_Dlambda_ver = - deux_tiers * sig_equi_ver * valder_ver.derivee * der_eps_plas_lambda_ver; // 4 - calcul de la variation de f Df_Dlambda_ver = D1_Dlambda_ver + D2_Dlambda_ver; // calcul de l'incrément de lambda delta_lambda_ver = - res_plas_ver / Df_Dlambda_ver; // nouvelle valeur de lambda lambda_ver += delta_lambda_ver; // 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 de la dérivée numérique de lambda D_sig_BH_DepsijBH_BH_ver(1,1) = (sigBH_ver(1,1) - sigBH(1,1))/ peti; double dersig1_ver =(3.*omega* eps_elasBH_ver(1,1) - sigBH(1,1))/ peti; double dersig1 = 3.*omega; double deromega = -2.*omega2*Dlambda_Depsij_BH_ver(1,1); double deromega_ver = (omega2-omega2)/peti; double dersig2 = -6.* omega2*eps_elasBH(1,1) * Dlambda_Depsij_BH(1,1); double dersig2_ver = (3.*omega_ver*eps_elasBH(1,1)-sigBH(1,1))/peti; double dersig = dersig1 + dersig2; */ // D_sig_BH_DepsijBH_BH(1,1) = 3.*omega*(1.-2.*omega*eps_elasBH(1,1)*Dlambda_Depsij_BH(1,1)); //=========== fin du calcul de la dérivée ========================================= // 4- calcul de la variation de sigmaHH par rapport au ddl for (int iddl = 1; iddl<= nbddl; iddl++) { // on fait de faire uniquement une égalité d'adresse et de ne pas utiliser // le constructeur d'ou la profusion d'* et de () Tenseur1HH & dsigHH = *((Tenseur1HH*) (d_sigHH(iddl))); // passage en dim 1 const Tenseur1BB & d_gijBB = *((Tenseur1BB*)(d_gijBB_tdt(iddl))); // passage en dim 1 const Tenseur1HH & dgijHH = *((Tenseur1HH*)(d_gijHH_tdt(iddl))) ; // pour simplifier l'ecriture // pour chacun des ddl on calcul les tenseurs derivees // 4.1- calcul de la variation de la déformation élastique / aux ddl Tenseur1BB d_eps_BB = 0.5 * d_gijBB; Tenseur1BH d_eps_elasBH = d_eps_BB * gijHH + epsBB * dgijHH; // 4.2- calcul de la variation de sigma en mixte Tenseur1BH dsigBH = D_sig_BH_DepsijBH_BH * d_eps_elasBH; // 4.3- en deux fois contravariant dsigHH = dgijHH * sigBH + gijHH * dsigBH; } } 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_Reuss1D::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_Reuss1D_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_Reuss1D::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); }; /* // **********calcul d'une dérivée numérique------------- double peti= 1.E-10; double lambda_ver = lambda+peti; double un_sur_1_plus_2_G_lambda_ver = 1. / (1. + deux_G*lambda_ver); double un_sur_1_plus_2_G_lambda2_ver = un_sur_1_plus_2_G_lambda_ver * un_sur_1_plus_2_G_lambda_ver; double alphaa_ver = un_sur_1_plus_2_G_lambda_ver*deux_G; double alphaa2_ver = un_sur_1_plus_2_G_lambda2_ver * deux_G_carre; double omega_ver = alphaa_ver*K/(alphaa_ver+2.*K); double omega2_ver = omega_ver*omega_ver; // a moins que lambda soit très grand on considère qu'omega est positif // par contre lambda peut-être négatif // if (lambda >= 0) double delta_eps_equi_ver = 2.*lambda_ver*omega_ver * eps_elasBH(1,1); // delta_eps_equi = 2.*lambda*omega * sqrt(c_c); double epsilon_barre_ver = epsilon_barre_t + // def plastique cumulée delta_eps_equi_ver; // nouvelle valeur de la variation de sigma_barre // et valeur de la tangente de la courbe d'écrouissage Courbe1D::ValDer valder_ver = f_ecrouissage->Valeur_Et_derivee(epsilon_barre_ver); double sig_equi_ver = valder_ver.valeur; // calcul du résidu double res_plas_ver = 3.*c_c*omega2_ver - sig_equi_ver * sig_equi_ver * un_tiers; double der1 = 3.*c_c*(omega2_ver-omega2)/peti; double dereps_barre = (epsilon_barre_ver -epsilon_barre)/peti; double der2 = un_tiers*(- sig_equi_ver * sig_equi_ver - - sig_equi * sig_equi)/peti; double der = (res_plas_ver - res_plas)/peti; double delta_lambda_ver = - res_plas/der ; // **********fin du calcul de la dérivée numérique------- */