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

3077 lines
158 KiB
C++
Raw Permalink Normal View History

ajout des lois de comportement
2021-09-23 11:21:15 +02:00
// FICHIER : LoiContraintesPlanesDouble.cp
// CLASSE : LoiContraintesPlanesDouble
// 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.
//
mise à jour 7.010
2023-05-03 17:23:49 +02:00
// Copyright (C) 1997-2022 Université Bretagne Sud (France)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
// 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 "LesLoisDeComp.h"
# include <iostream>
using namespace std; //introduces namespace std
#include <math.h>
#include <stdlib.h>
#include "Sortie.h"
#include "TypeQuelconqueParticulier.h"
#include "TypeConsTens.h"
#include "NevezTenseurQ.h"
#include "Util.h"
#include "ExceptionsLoiComp.h"
#include "MotCle.h"
#include "MathUtil.h"
#include "MathUtil2.h"
#include "LoiContraintesPlanesDouble.h"
// calcul des contraintes et ses variations par rapport aux déformations a t+dt
// ceci dans le cas où l'axe de traction simple est quelconque
// ici on suppose que:
// - tous les calculs sont en dim 3
// - l'axe 3 est normal aux deux autres (ceci pour Vi et pour gi)
//
// ViB et ViH : correspond à la base de traction: les vecteurs sont identiques
// par contre leurs composantes sont différentes car exprimés dans g^i et g_i
// telle que \vec ViB(1) = \vec ViH(1) = l'axe de traction
// La base Vi est supposée orthonormée
// 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 !!** non pas tout le temps a revoir
void LoiContraintesPlanesDouble::Calcul_dsigma_deps
(bool en_base_orthonormee,const BaseB * ViB,const BaseH * ViH
,const TenseurHH & sigHH_t,const TenseurBB& DepsBB
,const TenseurBB & epsBB_tdt,const TenseurBB & delta_epsBB,double& jacobien_0,double& jacobien
,TenseurHH& sigHH_tdt,TenseurHHHH& d_sigma_deps
,EnergieMeca & energ,const EnergieMeca &
,double& module_compressibilite,double& module_cisaillement
,const Met_abstraite::Umat_cont& ex)
{ // récup du conteneur spécifique
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
module_compressibilite=module_cisaillement=0.; // init
energ.Inita(0.); // initialisation des énergies mises en jeux
//récup de eps_mecaBB_t
// 1) dans le repère gi
Tenseur3BB& eps_mecaBB_t = *((Tenseur3BB*) (save_resul.eps_mecaBB_t));
Tenseur3BB& eps_mecaBB = *((Tenseur3BB*) (save_resul.eps_mecaBB));
// 2) dans le dernier repère V_a: a utiliser en sortie post, pas pour le calcul !!
Tenseur3BB& eps_P_mecaBB_t = *((Tenseur3BB*) (save_resul.eps_P_mecaBB_t));
Tenseur3BB& eps_P_mecaBB = *((Tenseur3BB*) (save_resul.eps_P_mecaBB));
// on récupère la dimension de l'espace de travail qui doit-être de dimension 3
// ainsi que le nombre de vecteur: pour l'instant 3
#ifdef MISE_AU_POINT
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{if (Permet_affichage() > 5)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
cout << "\n LoiContraintesPlanesDouble::Calcul_dsigma_deps(...";
if ((ViB->Dimension() != 3) || (ViH->Dimension() != 3))
{ cout << "\n***Erreur : la dimension de la base ViB = "<< ViB->Dimension()
<< " ou la dimension de la base ViH = " << ViH->Dimension()
<< " est different de 3, on ne peut pas continuer ! "
<< "\n LoiContraintesPlanesDouble::Calcul_dsigma_deps(...";
Sortie(1);
};
if ((ViB->NbVecteur() != 3) || (ViH->NbVecteur() != 3))
{ cout << "\n***Erreur : le nombre de vecteur de la base ViB = "
<< ViB->NbVecteur()
<< " ou le nombre de vecteur de la base ViH = " << ViH->NbVecteur()
<< " est different de 3, on ne peut pas continuer ! "
<< "\n LoiContraintesPlanesDouble::Calcul_dsigma_deps(...";
Sortie(1);
};
// idem pour l'Umat: là on ne test que les vecteurs à 0, pour alléger le test
if (ex.giB_0->Dimension() != 3)
{ cout << "\n***Erreur : la dimension de la base ex.giB_0 = "
<< ex.giB_0->Dimension()
<< " du conteneur Umat, est different de 3, on ne peut pas continuer ! "
<< "\n LoiContraintesPlanesDouble::Calcul_dsigma_deps(...";
Sortie(1);
};
if (ex.giB_0->NbVecteur() != 3)
{ cout << "\n***Erreur : le nombre de vecteur de la base ex.giB_0 = "
<< ex.giB_0->NbVecteur()
<< " du conteneur Umat est different de 3, on ne peut pas continuer ! "
<< "\n LoiContraintesPlanesDouble::Calcul_dsigma_deps(...";
Sortie(1);
};
}
#endif
// les coordonnées des vecteurs \vec V_e sont exprimées
// dans la base de travail g_i finale on a donc
// Vi_B(i) = {\beta'}_i^{.j} \hat{\vec g}_j = betaP(i,j) * giB_tdt(j)
// i indice de ligne, j indice de colonne
// \vec V_i = \vec V^i : les vecteurs sont égaux mais pas les composantes !
// 1 -- passage de la base \vec V_e à la base \hat{\vec g}_i
// et \vec V^e à la base \hat{\vec g}^i
// calcul des grandeurs gamma et beta:
// \hat{\vec g}^i = \gamma^i_{~.e}~\vec V^e
// \hat{\vec g}_i = \beta_i^{~.e}~\vec V_e
// 2 -- passage de la base \hat{\vec g}_i à la base \vec V_e
// et \hat{\vec g}^i à la base \vec V^e
//
// Vi_B(i) = {\beta'}_i^{.j} \hat{\vec g}_j
// Vi_H(i) = {\gamma'}^i_{.j} \hat{\vec g}^j
// gamma3D(i,e) = \gamma^i_{~.e} = \hat{\vec g}^i . \vec V_e
// ce qui donne que [gamma3D(i,e)] = [{\beta'}_i^{.j}]^T
for (int i=1;i<4;i++) // calcul de la matrice beta'
for (int e=1;e<4;e++) // au lieu de passer par un produit scalaire on récupère
{betaP_3D(i,e) = (*ViH)(i)(e); // directement les coordonnées
gammaP_3D(i,e) = (*ViB)(i)(e); // " "
};
// \beta_i^{~.e} = \hat{\vec g}_i . \vec V^e
// c'est aussi l'inverse de la transposé de gamma3D
// on choisit de calculer avec la formule
gamma3D = betaP_3D.Transpose(); // ok
beta3D = gammaP_3D.Transpose(); // betaP_3D.Inverse();
// gammaP_3D = beta3D.Transpose();
#ifdef MISE_AU_POINT
mise à jour 7.010
2023-05-03 17:23:49 +02:00
if (Permet_affichage() > 6)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
{cout << "\n Vi_B(1)= "<<(*ViB)(1) << " Vi_B(2)= "<<(*ViB)(2) << " Vi_B(3)= "<<(*ViB)(3);
cout << "\n Vi_H(1)= "<<(*ViH)(1) << " Vi_H(2)= "<<(*ViH)(2) << " Vi_H(3)= "<<(*ViH)(3);
cout << "\n Vi_B(i) = {beta'}_i^{.j} hat{vec g}_j --> beta'= ";betaP_3D.Affiche();
cout << "\n Vi_H(i) = {gamma'}^i_{.j} hat{vec g}^j--> gamma'= ";gammaP_3D.Affiche();
cout << "\n hat{vec g}^i = gamma^i_{~.e}~vec V^e --> gamma= ";gamma3D.Affiche();
cout << "\n hat{vec g}_i = beta_i^{~.e}~vec V_e--> beta= ";beta3D.Affiche();
};
#endif
// on alimente les pointeurs de tenseurs d'entrée de la cinématique
epsBB_tdt_cin = ((Tenseur3BB*) &epsBB_tdt); // passage en dim 3
DepsBB_cin = ((Tenseur3BB*) &DepsBB); // passage en dim 3
delta_epsBB_cin = ((Tenseur3BB*) &delta_epsBB); // passage en dim 3
#ifdef MISE_AU_POINT
mise à jour 7.010
2023-05-03 17:23:49 +02:00
if (Permet_affichage() > 6)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
{cout << "\n epsBB_tdt_cin "<< (*epsBB_tdt_cin) << "\n DepsBB_cin "<< (*DepsBB_cin)
<< "\n delta_epsBB_cin "<< (*delta_epsBB_cin) << flush;
// en base orthonormee
Tenseur3BB epsBB_tdt_cin_inter;
epsBB_tdt_cin->BaseAbsolue(epsBB_tdt_cin_inter,(* ex.giH_tdt));
cout << "\n epsBB_tdt_cin= en absolue 3D : "<<epsBB_tdt_cin_inter;
};
#endif
// --- initialisation des déformations mécaniques dans le repère V_a
eps_P_BB = eps_mecaBB_t; // init: là on est dans le repère g^i
// on passe ensuite dans le repère V_a qui est donc le nouveau repère, pour ce passage
eps_P_BB.ChBase(betaP_3D);
#ifdef MISE_AU_POINT
mise à jour 7.010
2023-05-03 17:23:49 +02:00
if (Permet_affichage() > 6)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
{cout << "\n eps_mecaBB_t= "; eps_mecaBB_t.Ecriture(cout);
cout << "\n eps_P_BB apres change base donc en V_a "<< eps_P_BB;
};
#endif
// on calcul la déformation cinématique dans la direction V_1
// \epsilon'_{11}(t+\Delta t) = \epsilon_{kl}(t+\Delta t) ~\gamma^k_{~.1} ~\gamma^l_{~.1}
// idem pour l'accroissement de déformation et la vitesse de déformation
eps_P_BB.Coor(1,1) = 0.; Deps_P_BB.Coor(1,1) = 0.; delta_eps_P_BB.Coor(1,1) = 0.;
for (int k=1;k<4;k++) // calcul de la matrice gamma
for (int l=1;l<4;l++)
{eps_P_BB.Coor(1,1) += (*epsBB_tdt_cin)(k,l) * gamma3D(k,1) * gamma3D(l,1);
Deps_P_BB.Coor(1,1) += (*DepsBB_cin)(k,l) * gamma3D(k,1) * gamma3D(l,1);
delta_eps_P_BB.Coor(1,1) += (*delta_epsBB_cin)(k,l) * gamma3D(k,1) * gamma3D(l,1);
};
#ifdef MISE_AU_POINT
mise à jour 7.010
2023-05-03 17:23:49 +02:00
if (Permet_affichage() > 6)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
{cout << "\n eps_P_BB avec la cinematique "<< eps_P_BB;
cout << "\n Deps_P_BB avec la cinematique "<< Deps_P_BB;
cout << "\n delta_eps_P_BB avec la cinematique "<< delta_eps_P_BB << flush;
};
#endif
// on affecte les contraintes a t: en fait pour les lois incrémentales, les infos à t
// sont déjà stocké, donc cela ne sert à rien, mais pour une loi élastique par exemple ce n'est pas le cas
sig_HH_t_3D = sigHH_t;
// concernant la métrique relative au repère de travail, on utilise celle passée en paramètre
(*umat_cont_3D)= ex;
// concernant le jacobien, le jacobien du conteneur est supposé être celui 1D
sauve_jacobien1D_tdt = *(ex.jacobien_tdt);
// --- pour les contraintes passage en 3D
// bool plusZero = true;
// --- pour la cinématique
// passage des informations spécifique à la loi le_SaveResul
lois_interne->IndiqueSaveResult(save_resul.le_SaveResul);
lois_interne->IndiquePtIntegMecaInterne(ptintmeca_en_cours);// idem pour ptintmeca
lois_interne->IndiqueDef_en_cours(def_en_cours); // idem pour def en coursVhVg
// la base de traction
ViB_3D = ViB; ViH_3D = ViH;
// calcul d'élements intermédiaires
V1V1_BB = Tenseur3BB::Prod_tensoriel((*ViB_3D)(1),(*ViB_3D)(1));
#ifdef MISE_AU_POINT
if (Permet_affichage() > 7)
{ cout << "\n V1V1_BB: "<< V1V1_BB;
}
#endif
// pour (m,n) = (2,2); (3,3); (1,2), et pour (g,h) = (2,2); (3,3); et (1,2)
// on va utiliser des matrices intermédiaires tels que:
// in = 1 pour (2,2) , 2 pour (3,3) et 3 pour (1,2)
for (int gh=1;gh<4;gh++)
{//if (gh <3)
{VhVg_BB(gh) = Tenseur3BB::Prod_tensoriel((*ViB_3D)(jnd(gh)),(*ViB_3D)(ind(gh)));}
// else // donne des tenseurs non symétriques pour (1,2) -> on symétrise
// {VhVg_BB(gh) = 0.5* Tenseur3BB::Prod_tensoriel((*ViB_3D)(ind(gh)),(*ViB_3D)(jnd(gh)))
// + 0.5* Tenseur3BB::Prod_tensoriel((*ViB_3D)(jnd(gh)),(*ViB_3D)(ind(gh)));
// }
#ifdef MISE_AU_POINT
if (Permet_affichage() > 7)
{ cout << "\n VhVg_BB("<<gh<<")= "<< VhVg_BB(gh);
}
#endif
};
// on sauvegarde les dimensions transverse, car elles vont être modifiées et si on doit
// repartir au départ, ce ne sera plus possible sans cette sauvegarde
double sauve_hsurh0=save_resul.hsurh0;
double sauve_bsurb0=save_resul.bsurb0;
// on initialise les jacobiens à 0 et t
jacobien_0_3D = jacobien_0; // à t=0 -> save_resul.hsurh0 * save_resul.bsurb0 == 1;
jacobien_t_3D = (*ex.jacobien_t) * save_resul.h_tsurh0 * save_resul.b_tsurb0;
// choix entre calcul avec une boucle locale de Newton on non.
// Dans le premier cas il nous faut la version d_sig/d_eps
// au lieu de d_sig/d_ddl
bool mauvaise_convergence = false; // init
bool newtonlocal = false;
switch (type_de_contrainte)
{case NEWTON_LOCAL: // cas de la version d_sig/d_eps
{ newtonlocal=true; // pour la gestion du cas par défaut
// initialisation du tenseur final de contrainte
sig_HH_3D = sigHH_tdt; // c-a-d de la solution: peut-être pas nécessaire ***
// passage des informations liées à la déformation de 1 vers 3, et variation de volume
Tableau <TenseurBB *>* toto = NULL; // pour dire que ce n'est pas attribué
// on appel la procédure de résolution de sig33(..)=0
if ((sortie_post)&&(save_resul.indicateurs_resolution.Taille()!= 2)) // dimensionnement éventuelle de la sortie d'indicateurs
save_resul.indicateurs_resolution.Change_taille(2);
// ----- pour ce faire on appelle une méthode de recherche de zero
// on initialise les inconnues à 0. On pourrait partir des def_P sauvegardés
// mais si on change d'orientation dans la pratique la convergence est moins bonne ??
val_initiale_3D.Zero(); // = save_resul.def_P_t;
// val_initiale_3D = save_resul.def_P_t;
// on impose que les grandeurs soient dans les limites admises
double& hsurh0 = save_resul.h_tsurh0;double& bsurb0 = save_resul.b_tsurb0; // init
// la méthode va mettre à jour les variations bsurb0 et hsurh0 et le jacobien
if (Calcul_et_Limitation2_h_b(hsurh0,eps_P_BB,bsurb0,jacobien_tdt_3D,jacobien_0_3D))
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{if (Permet_affichage() > 3)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
cout << "\n pb limitation val initiales: LoiContraintesPlanesDouble::Calcul_dsigma_deps(... "
<< flush;
};
//// ---- debut
// // on vérifie que les anciennes transversales ne sont pas négatives
// if ((val_initiale(1) < 0.) || (val_initiale(2) < 0.) )
// { cout << "\n debug LoiContraintesPlanesDouble::Calcul_DsigmaHH_tdt : "
// << " initialement "
// << "\n save_resul.b_tsurb0 = " << save_resul.b_tsurb0
// << ", save_resul.h_tsurh0 " << save_resul.h_tsurh0 << endl;
// Sortie(1);
// };
////------ fin debug
residu_3D.Zero(); // init du résultat à 0., mais c'est une grandeur de sortie donc cela ne sert à rien
derResidu_3D.Initialise(0.); // init de la matrice dérivée à 0.
try // on met le bloc sous surveillance
{
// résolution de l'équation constitutive d'avancement discrétisée en euler implicite
int nb_incr_total,nb_iter_total; // variables intermédiaires d'indicateurs
// dans le cas où on utilise une précision qui est pilotée
{ // opération de transmission de la métrique: encapsulé ici
const Met_abstraite::Impli* ex_impli = NULL;
const Met_abstraite::Expli_t_tdt* ex_expli_tdt = NULL;
const Met_abstraite::Umat_cont* ex_expli = &ex;
if (fct_tolerance_residu != NULL)
{// ici on utilise les variables connues aux pti, ou calculées à partir de
// on commence par récupérer les conteneurs des grandeurs à fournir
List_io <Ddl_enum_etendu>& li_enu_scal = fct_tolerance_residu->Li_enu_etendu_scalaire();
List_io <TypeQuelconque >& li_quelc = fct_tolerance_residu->Li_equi_Quel_evolue();
bool absolue = true; // on se place systématiquement en absolu
// on va utiliser la méhode Valeur_multi_interpoler_ou_calculer
// pour les grandeurs strictement scalaire
Tableau <double> val_ddl_enum(Valeur_multi_interpoler_ou_calculer
(absolue,TEMPS_tdt,li_enu_scal,ex_impli,ex_expli_tdt,ex_expli,NULL)
);
// on utilise la méthode Valeurs_Tensorielles_interpoler_ou_calculer
// pour les Coordonnees et Tenseur
Valeurs_Tensorielles_interpoler_ou_calculer
(absolue,TEMPS_tdt,li_quelc,ex_impli,ex_expli_tdt,ex_expli,NULL);
// calcul de la valeur et retour dans tab_ret
Tableau <double> & tab_val = fct_tolerance_residu->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
#ifdef MISE_AU_POINT
if (tab_val.Taille() != 1)
{ cout << "\nErreur : la fonction nD de pilotage de la precision "
<< " doit calculer un scalaire or le tableau de retour est de taille "
<< tab_val.Taille() << " ce n'est pas normal !\n";
cout << " LoiContraintesPlanesDouble::Calcul_dsigma_deps\n";
Sortie(1);
};
#endif
// on récupère le premier élément du tableau uniquement
double tol = tab_val(1);
#ifdef MISE_AU_POINT
mise à jour 7.010
2023-05-03 17:23:49 +02:00
if (Permet_affichage() > 5)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
cout << "\n Newton_prec: tol_= "<< tol;
#endif
alg_zero.Modif_prec_res_abs(tol);
};
if (fct_tolerance_residu_rel != NULL)
{// ici on utilise les variables connues aux pti, ou calculées à partir de
// on commence par récupérer les conteneurs des grandeurs à fournir
List_io <Ddl_enum_etendu>& li_enu_scal = fct_tolerance_residu_rel->Li_enu_etendu_scalaire();
List_io <TypeQuelconque >& li_quelc = fct_tolerance_residu_rel->Li_equi_Quel_evolue();
bool absolue = true; // on se place systématiquement en absolu
// on va utiliser la méhode Valeur_multi_interpoler_ou_calculer
// pour les grandeurs strictement scalaire
Tableau <double> val_ddl_enum(Valeur_multi_interpoler_ou_calculer
(absolue,TEMPS_tdt,li_enu_scal,ex_impli,ex_expli_tdt,ex_expli,NULL)
);
// on utilise la méthode Valeurs_Tensorielles_interpoler_ou_calculer
// pour les Coordonnees et Tenseur
Valeurs_Tensorielles_interpoler_ou_calculer
(absolue,TEMPS_tdt,li_quelc,ex_impli,ex_expli_tdt,ex_expli,NULL);
// calcul de la valeur et retour dans tab_ret
Tableau <double> & tab_val = fct_tolerance_residu_rel->Valeur_FnD_Evoluee(&val_ddl_enum,&li_enu_scal,&li_quelc,NULL,NULL);
#ifdef MISE_AU_POINT
if (tab_val.Taille() != 1)
{ cout << "\nErreur : la fonction nD de pilotage de la precision relative "
<< " doit calculer un scalaire or le tableau de retour est de taille "
<< tab_val.Taille() << " ce n'est pas normal !\n";
cout << " LoiContraintesPlanesDouble::Calcul_dsigma_deps\n";
Sortie(1);
};
#endif
// on récupère le premier élément du tableau uniquement
double tol_rel = tab_val(1);
#ifdef MISE_AU_POINT
mise à jour 7.010
2023-05-03 17:23:49 +02:00
if (Permet_affichage() > 5)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
cout << ", tol_relative= "<< tol_rel;
#endif
alg_zero.Modif_prec_res_rel(tol_rel);
};
};
bool conver=alg_zero.Newton_raphson
(*this,&LoiContraintesPlanesDouble::Residu_constitutif_3D
,&LoiContraintesPlanesDouble::Mat_tangente_constitutif_3D
,val_initiale_3D,racine_3D,derResidu_3D,nb_incr_total,nb_iter_total
,maxi_delta_var_eps_sur_iter_pour_Newton);
if(sortie_post) // sauvegarde éventuelle des indicateurs
{save_resul.indicateurs_resolution(1)+=nb_incr_total;
save_resul.indicateurs_resolution(2)+=nb_iter_total;
};
// on vérifie qu'il n'y a pas de pb de convergence
double absracinemax=racine_3D.Max_val_abs();
if ((!conver) || (!isfinite(absracinemax)) || (isnan(absracinemax)) )
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{ if (Permet_affichage() > 0)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
{ cout << "\n**>> non convergence sur l'algo de la resolution de (sigma.V_2=0, sigma.V_3=0),\n";
if (ptintmeca_en_cours != NULL) ptintmeca_en_cours->Signature();
cout << ", nb_incr_total=" << nb_incr_total
<< ", nb_iter_total=" << nb_iter_total<<", val obtenues:"
<< "\n eps_P(2,2)= " << racine_3D(1) << ", "
<< " eps_P(3,3)= " << racine_3D(2) << ", "<< " eps_P(1,2)= " << racine_3D(3)
<< "\n precedentes: "<<eps_mecaBB_t(2,2)<<", "<<eps_mecaBB_t(3,3)<<", "<<eps_mecaBB_t(1,2);
cout << "\n dernier_residu: "<< residu_3D(1) <<", "<<residu_3D(2) <<", "<<residu_3D(3)
<< "\n LoiContraintesPlanesDouble::Calcul_dsigma_deps (...";
};
// on garde en mémoire
mauvaise_convergence=true;
};
}
catch (ErrNonConvergence_Newton erreur)
{
mise à jour 7.010
2023-05-03 17:23:49 +02:00
if (Permet_affichage() > 0)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
{ cout << "\n erreur de non convergence avec l'algo de la resolution de"
<< " (sig'22, sig'33, sig'12)=(0,0,0) "
<< " on obtient une valeur infinie ou NaN "
<< "\n dernier_residu: "<< residu_3D(1) <<", "<<residu_3D(2) <<", "<<residu_3D(3)
<< "\n LoiContraintesPlanesDouble::Calcul_dsigma_deps (...";
};
// on garde en mémoire
mauvaise_convergence=true;
}
catch (ErrSortieFinale)
// cas d'une direction voulue vers la sortie
// on relance l'interuption pour le niveau supérieur
{ ErrSortieFinale toto;
throw (toto);
}
catch ( ... )
{ // dans le cas d'une erreur inconnue, on génère également une exception
mise à jour 7.010
2023-05-03 17:23:49 +02:00
if (Permet_affichage() > 0)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
{ cout << "\n erreur non identifiee sur l'algo de la resolution de "
<< " (sig'22, sig'33, sig'12)=(0,0,0) "
<< "\n dernier_residu: "<< residu_3D(1) <<", "<<residu_3D(2) <<", "<<residu_3D(3)
<< "\n LoiContraintesPlanesDouble::Calcul_dsigma_deps (...";
};
// on garde en mémoire
mauvaise_convergence=true;
};
if (!mauvaise_convergence)
{ // arrivée ici, cela veut dire que tout à bien fonctionné
save_resul.def_P = racine_3D; // récup de la solution dans le repère V_a
// on doit sauvegarder les def méca dans le repère g^i, ce sont directement
// les dernières composantes utilisées pour l'appel de la loi 3D
eps_mecaBB = eps_BB_3D; // dans g^i
eps_P_mecaBB = eps_P_BB; // dans V_a pour le post traitement
// on met à jour les modules
module_compressibilite= module_compressibilite_3D;
module_cisaillement= module_cisaillement_3D;
// retour du tenseurs résultat
sigHH_tdt = sig_HH_3D;
// calcul de l'opérateur tangent final: ici en 3D
// d sigma^{ij} / d epsilon_{kl}
Calcul2_dsigma_deps(d_sigma_deps);
#ifdef MISE_AU_POINT
mise à jour 7.010
2023-05-03 17:23:49 +02:00
if (Permet_affichage() > 6)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
{ cout << "\nLoiContraintesPlanesDouble::Calcul_dsigma_deps: bonne convergence de Newton ";
TenseurHH* ptHH = NevezTenseurHH(sig_HH_3D);
sig_HH_3D.BaseAbsolue(*ptHH,giB_tdt_3D);
cout << "\n sigma apres CP2: ";
ptHH->Ecriture(cout);
delete ptHH;
};
if (Permet_affichage() > 5)
{ cout << "\n def_P(2,2) " << save_resul.def_P(1)
<< " def_P(3,3) " << save_resul.def_P(2)
<< " def_P(1,2) " << save_resul.def_P(3)
<< "\n en g^i: eps_mecaBB.Coor(2,2) " << eps_mecaBB.Coor(2,2)
<< " eps_mecaBB.Coor(3,3) " << eps_mecaBB.Coor(3,3)
<< " eps_mecaBB.Coor(1,2) " << eps_mecaBB.Coor(1,2)
<< "\n en V_a: eps_P_mecaBB.Coor(2,2) " << eps_P_mecaBB.Coor(2,2)
<< " eps_P_mecaBB.Coor(3,3) " << eps_P_mecaBB.Coor(3,3)
<< " eps_P_mecaBB.Coor(1,2) " << eps_P_mecaBB.Coor(1,2)
<< "\n module_compressibilite= " << module_compressibilite
<< " module_cisaillement= " << module_cisaillement;
// on affiche les épaisseur et largeur
double& hsurh0 = save_resul.hsurh0;double& bsurb0 = save_resul.bsurb0;
cout << "\n hsurh0= "<< hsurh0 << " bsurb0= " << bsurb0 << flush;
};
if (Permet_affichage() > 8)
{ cout << "\n dans le repere locale d_sigma_deps_3D= \n";
int e=1;
for (int i=1;i<4;i++) for (int j=1;j<4;j++) for (int k=1;k<4;k++)for (int l=1;l<4;l++,e++)
{ cout << "("<<i<<","<<j<<","<<k<<","<<l<<")= "<<d_sigma_deps_3D(i,j,k,l) << " ; ";
if (e>6) {cout << "\n"; e=1;}
};
Tenseur3HHHH inter_HHHH;
d_sigma_deps_3D.ChangeBase(inter_HHHH,giB_tdt_3D);
cout << "\n dans le repere orthonormee d_sigma_deps_3D= \n";
e=1;
for (int i=1;i<4;i++) for (int j=1;j<4;j++) for (int k=1;k<4;k++)for (int l=1;l<4;l++,e++)
{ cout << "("<<i<<","<<j<<","<<k<<","<<l<<")= "<<inter_HHHH(i,j,k,l) << " ; ";
if (e>6) {cout << "\n"; e=1;}
};
};
#endif
break;
};
// sinon on ne fait pas de break, donc on continue avec la méthode de perturbation
};
default: // cas de la version avec perturbation, on encore du cas où la méthode de Newton n'a pas convergé
{
// initialisation du tenseurs contrainte
sig_HH_3D.Inita(0.);
// dans le cas ou on est passé par la méthode de newton et que cela n'a pas
// convergé, on doit réinitialiser les variables de départ
mise à jour 7.010
2023-05-03 17:23:49 +02:00
if ((newtonlocal) && (Permet_affichage() > 3))
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
cout << "\n calcul direct sans Newton ";
save_resul.hsurh0 = sauve_hsurh0;
save_resul.bsurb0 = sauve_bsurb0;
// appel du calcul en perturbation
Cal_avec_perturbation2();
// on met à jour les modules
module_compressibilite= module_compressibilite_3D;
module_cisaillement= module_cisaillement_3D;
if (mauvaise_convergence)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{if (Permet_affichage() > 0)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
{ cout << ": valeur final de h/h_0(t+dt)= " << save_resul.hsurh0
<< " et de b/b_0(t+dt)= " << save_resul.bsurb0;
};
// on annulle les dérivées des épaisseurs
save_resul.d_hsurh0.Zero();save_resul.d_bsurb0.Zero();
};
}
break;
};
// sig_HH_3D.Ecriture(cout);
(*save_resul.l_sigoHH) = sig_HH_3D; // sauvegarde en locale
// ---calcul des invariants de déformation et de vitesse de déformation
Calcul_invariants_et_def_cumul();
energ = save_resul.l_energ; // récup des énergies
// --- on remet à jour éventuellement l'épaisseur
// dans le cas d'une contrainte par perturbation
// if ((type_de_contrainte == PERTURBATION )||mauvaise_convergence)
if (type_de_contrainte == PERTURBATION )
// calcul de la déformation d'épaisseur correspondant à la condition de contraintes planes
{ //Tenseur1BH sigBH = gijBB_tdt * sigHH_tdt;
Tenseur3BH sigBH = gijBB_tdt_3D * sig_HH_3D;
Tenseur3BH sigBH_t = gijBB_t_3D * sig_HH_t_3D;
Calcul_eps_trans_parVarVolume2(sigBH_t,jacobien_0,module_compressibilite,jacobien,sigBH,*(ex.jacobien_t));
mise à jour 7.010
2023-05-03 17:23:49 +02:00
if ((Permet_affichage() > 0) && mauvaise_convergence)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
{ cout << "\n valeur final de h/h_0(t+dt)= " << save_resul.hsurh0
<< " et de b/b_0(t+dt)= " << save_resul.bsurb0;
};
};
LibereTenseur();
LibereTenseurQ();
};
// fonction interne utilisée par les classes dérivées de Loi_comp_abstraite
// pour répercuter les modifications de la température
// ici utiliser pour modifier la température des lois élémentaires
// l'Enum_dure: indique quel est la température courante : 0 t ou tdt
void LoiContraintesPlanesDouble::RepercuteChangeTemperature(Enum_dure temps)
{ lois_interne->temperature_0 = this->temperature_0;
lois_interne->temperature_t = this->temperature_t;
lois_interne->temperature_tdt = this->temperature_tdt;
lois_interne->dilatation=dilatation;
// on répercute également les déformations thermiques, qui ne sont utilisées
// telles quelles que pour certaines lois: ex: loi hyper-élastique
if (dilatation)
{// a- dimensionnement des tenseurs intermédiaires
int dim_tens = epsBB_therm->Dimension();
// -- cas de la déformation
if (lois_interne->epsBB_therm == NULL) { lois_interne->epsBB_therm = NevezTenseurBB(dim_tens);}
else if (lois_interne->epsBB_therm->Dimension() != dim_tens)
{ delete lois_interne->epsBB_therm;lois_interne->epsBB_therm = NevezTenseurBB(dim_tens);};
// -- cas de la vitesse de déformation
if (lois_interne->DepsBB_therm == NULL) { lois_interne->DepsBB_therm = NevezTenseurBB(dim_tens);}
else if (lois_interne->DepsBB_therm->Dimension() != dim_tens)
{ delete lois_interne->DepsBB_therm;lois_interne->DepsBB_totale = NevezTenseurBB(dim_tens);};
// b- affectation des tenseurs
(*lois_interne->epsBB_therm)=(*epsBB_therm);
(*lois_interne->DepsBB_therm)=(*DepsBB_therm);
};
// puis en interne
lois_interne->RepercuteChangeTemperature(temps);
switch (temps)
{ case TEMPS_0:
{lois_interne->temperature = &lois_interne->temperature_0;
break;
}
case TEMPS_t:
{lois_interne->temperature = &lois_interne->temperature_t;
break;
}
case TEMPS_tdt:
{lois_interne->temperature = &lois_interne->temperature_tdt;
break;
}
default:
{ cout << "\n erreur, cas de temps non prevu !! "
<< "\n LoiContraintesPlanesDouble::RepercuteChangeTemperature(...";
Sortie(1);
};
};
};
void LoiContraintesPlanesDouble::CalculGrandeurTravail
(const PtIntegMecaInterne& ptintmec
,const Deformation & def,Enum_dure temps,const ThermoDonnee& dTP
,const Met_abstraite::Impli* ex_impli
,const Met_abstraite::Expli_t_tdt* ex_expli_tdt
,const Met_abstraite::Umat_cont* ex_umat
,const List_io<Ddl_etendu>* exclure_dd_etend
,const List_io<const TypeQuelconque *>* exclure_Q
)
{ // récup des infos spécifiques
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
// recup de l'incrément de temps
double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();
double unSurDeltat=0;
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
{ // non 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;
};
// par défaut on active les invariants et on intègre les infos def en 3D
if (!calcul_en_3D_via_direction_quelconque)
{ //dans le cas classique, on dispose de la variation d'épaisseur et de largeur pour déterminer les def
// ptintmec concerne des tenseurs d'ordres 1 donc il faut transformer en 3D
bool plusZero=true; // on ajoute des 0 pour les grandeurs manquantes dans une première étape
// on commence par transférer le tout
ptintmeca.Affectation_1D_a_3D(ptintmec,plusZero);
// ensuite on traite les points particuliers
// pour tout ce qui est contrainte: l'ajout de 0 est ok
// pour les déformations on met à jour en fonction des grandeurs sauvegardées
double eps33=0.; double eps33_t=0.; //init par défaut
double eps22=0.; double eps22_t=0.; //init par défaut
// -- maintenant on tient compte de l'élongation d'épaisseur
switch (type_de_deformation)
{case DEFORMATION_STANDART : case DEFORMATION_POUTRE_PLAQUE_STANDART :
// cas d'une déformation d'Almansi
{ // dans le repère local: epsBB33 = 1/2 * (h^2 - 1.), or h0=1. donc : epsBB33 = 1/2 * ((h/h0)^2 - 1.)
eps33 = 0.5 * (save_resul.hsurh0 * save_resul.hsurh0 - 1.);
eps33_t = 0.5 * (save_resul.h_tsurh0 * save_resul.h_tsurh0 - 1.);
eps22 = 0.5 * (save_resul.bsurb0 * save_resul.bsurb0 - 1.);
eps22_t = 0.5 * (save_resul.b_tsurb0 * save_resul.b_tsurb0 - 1.);
};
break;
case DEFORMATION_LOGARITHMIQUE : case DEF_CUMUL_CORROTATIONNEL : case DEFORMATION_CUMU_LOGARITHMIQUE :
// cas d'une def logarithmique ou une approximation
{ eps33 = log(save_resul.hsurh0);
eps33_t = log(save_resul.h_tsurh0);
eps22 = log(save_resul.bsurb0);
eps22_t = log(save_resul.b_tsurb0);
};
break;
default :
cout << "\nErreur : type de deformation qui n'est pas actuellement pris en compte, type= "
<< Nom_type_deformation(type_de_deformation);
cout << "\n LoiContraintesPlanesDouble::CalculGrandeurTravail \n";
Sortie(1);
};
(*ptintmeca.EpsBB()).Coor(3,3) = eps33;
(*ptintmeca.DeltaEpsBB()).Coor(3,3) = eps33 - eps33_t;
(*ptintmeca.EpsBB()).Coor(2,2) = eps22;
(*ptintmeca.DeltaEpsBB()).Coor(2,2) = eps22 - eps22_t;
(*ptintmeca.DepsBB()).Coor(3,3) = (eps33 - eps33_t) * unSurDeltat;
(*ptintmeca.DepsBB()).Coor(2,2) = (eps22 - eps22_t) * unSurDeltat;
}
else
// dans le cas où on travaille dans le repère gi en 3D on dispose directement des infos de déformation
{//récup de eps_mecaBB
Tenseur3BB& eps_mecaBB_t = *((Tenseur3BB*) (save_resul.eps_mecaBB_t));
Tenseur3BB& eps_mecaBB = *((Tenseur3BB*) (save_resul.eps_mecaBB));
// on sauve ce qu'il y a dans la direction 11 pour éviter les erreurs d'arrondis
double eps11= (*ptintmeca.EpsBB())(1,1);
double Deps11= (*ptintmeca.DepsBB())(1,1);
double delta_eps11= (*ptintmeca.DeltaEpsBB())(1,1);
// les manips entre tenseurs
(*ptintmeca.EpsBB()) = eps_mecaBB;
(*ptintmeca.DeltaEpsBB()) = eps_mecaBB - eps_mecaBB_t;
(*ptintmeca.DepsBB()) = (*ptintmeca.DeltaEpsBB()) * unSurDeltat;
// récup des valeurs sauvegardées (normalement ne devrait pas changer...)
(*ptintmeca.EpsBB()).Coor(1,1) = eps11;
(*ptintmeca.DeltaEpsBB()).Coor(1,1) = delta_eps11;
(*ptintmeca.DepsBB()).Coor(1,1) = Deps11;
};
// on met à jour les invariants 3D, on récupère ce qui a été sauvegardé
if (ptintmeca.Statut_Invariants_deformation()) // le conteneur des invariants a pu être effacé par Affectation_1D_a_3D
ptintmeca.EpsInvar() = save_resul.epsInvar;
if (ptintmeca.Statut_Invariants_vitesseDeformation()) // le conteneur des invariants a pu être effacé par Affectation_1D_a_3D
ptintmeca.DepsInvar() = save_resul.depsInvar;
// idem au niveau des déformations cumulées
ptintmeca.Deformation_equi() = save_resul.def_equi;
ptintmeca.Deformation_equi_t() = save_resul.def_equi_t;
// passage des informations spécifique à la loi le_SaveResul
lois_interne->IndiqueSaveResult(save_resul.le_SaveResul);
lois_interne->IndiquePtIntegMecaInterne(ptintmeca_en_cours);// idem pour ptintmeca
lois_interne->IndiqueDef_en_cours(def_en_cours); // idem pour def en cours
lois_interne->CalculGrandeurTravail(ptintmeca,def,temps,dTP
,ex_impli,ex_expli_tdt,ex_umat,exclure_dd_etend,exclure_Q); // passage à la loi 3D
};
// passage des grandeurs métriques de l'ordre 1 à 3, cas implicite
void LoiContraintesPlanesDouble::Passage_metrique_ordre1_vers_3(const Met_abstraite::Impli& ex)
{// on s'occupe du redimensionnement éventuel
// on affecte et/ou on redimensionne éventuellement les tableaux fonction du nombre de ddl pour le passage 3D
// la partie dépendant des vitesses: entre accolades pour pouvoir fermer
{if (ex.gradVmoyBB_t != NULL) {impli_3D->gradVmoyBB_t= gradVmoyBB_t_3D_P = &gradVmoyBB_t_3D;};
if (ex.gradVmoyBB_tdt != NULL) {impli_3D->gradVmoyBB_tdt=gradVmoyBB_tdt_3D_P = &gradVmoyBB_tdt_3D;};
if (ex.gradVBB_tdt != NULL) {impli_3D->gradVBB_tdt=gradVBB_tdt_3D_P = &gradVBB_tdt_3D;};
// on s'occupe tout d'abord des tableaux de vitesses qui peuvent ne pas exister
if (ex.d_gradVmoyBB_t != NULL) // cas où il existe
{ int tail=ex.d_gradVmoyBB_t->Taille();
if (d_gradVmoyBB_t_3D_P == NULL) {d_gradVmoyBB_t_3D_P = new Tableau<TenseurBB *> (tail);
impli_3D->d_gradVmoyBB_t = d_gradVmoyBB_t_3D_P;}
else {d_gradVmoyBB_t_3D_P->Change_taille(tail);};
d_gradVmoyBB_t_3D.Change_taille(tail);
for (int i=1;i<= tail;i++)
(*d_gradVmoyBB_t_3D_P)(i) = &(d_gradVmoyBB_t_3D(i));
};
if (ex.d_gradVmoyBB_tdt != NULL) // cas où il existe
{ int tail=ex.d_gradVmoyBB_tdt->Taille();
if (d_gradVmoyBB_tdt_3D_P == NULL) {d_gradVmoyBB_tdt_3D_P = new Tableau<TenseurBB *> (tail);
impli_3D->d_gradVmoyBB_tdt = d_gradVmoyBB_tdt_3D_P;}
else {d_gradVmoyBB_tdt_3D_P->Change_taille(tail);};
d_gradVmoyBB_tdt_3D.Change_taille(tail);
for (int i=1;i<= tail;i++)
(*d_gradVmoyBB_tdt_3D_P)(i) = &(d_gradVmoyBB_tdt_3D(i));
};
if (ex.d_gradVBB_t != NULL) // cas où il existe
{ int tail=ex.d_gradVBB_t->Taille();
if (d_gradVBB_t_3D_P == NULL) {d_gradVBB_t_3D_P = new Tableau<TenseurBB *> (tail);
impli_3D->d_gradVBB_t = d_gradVBB_t_3D_P;}
else {d_gradVBB_t_3D_P->Change_taille(tail);};
d_gradVBB_t_3D.Change_taille(tail);
for (int i=1;i<= tail;i++)
(*d_gradVBB_t_3D_P)(i) = &(d_gradVBB_t_3D(i));
};
if (ex.d_gradVBB_tdt != NULL) // cas où il existe
{ int tail=ex.d_gradVBB_tdt->Taille();
if (d_gradVBB_tdt_3D_P == NULL) {d_gradVBB_tdt_3D_P = new Tableau<TenseurBB *> (tail);
impli_3D->d_gradVBB_tdt = d_gradVBB_tdt_3D_P;}
else {d_gradVBB_tdt_3D_P->Change_taille(tail);};
d_gradVBB_tdt_3D.Change_taille(tail);
for (int i=1;i<= tail;i++)
(*d_gradVBB_tdt_3D_P)(i) = &(d_gradVBB_tdt_3D(i));
};
}; // fin de la partie dédiée à la vitesse
// maintenant on s'occupe uniquement du redimensionnement des tableaux restants
// -- on s'occupe des tableaux nécessaire à la métrique
int ta_ex_d_giB_tdt = ex.d_giB_tdt->Taille(); // la dimension
if (d_giB_tdt_3D.Taille() != ta_ex_d_giB_tdt)
{ d_giB_tdt_3D.Change_taille(ta_ex_d_giB_tdt);};
int ta_ex_d_giH_tdt = ex.d_giH_tdt->Taille(); // la dimension
if (d_giH_tdt_3D.Taille() != ta_ex_d_giH_tdt)
{ d_giH_tdt_3D.Change_taille(ta_ex_d_giH_tdt);};
int ta_ex_d_gijBB_tdt = ex.d_gijBB_tdt->Taille(); // la dimension
if (d_gijBB_tdt_3D.Taille() != ta_ex_d_gijBB_tdt)
{ d_gijBB_tdt_3D.Change_taille(ta_ex_d_gijBB_tdt); d_gijBB_tdt_3D_P.Change_taille(ta_ex_d_gijBB_tdt);
for (int i=1;i<=ta_ex_d_gijBB_tdt;i++) {d_gijBB_tdt_3D_P(i) = &(d_gijBB_tdt_3D(i));};
};
int ta_ex_d_gijHH_tdt = ex.d_gijHH_tdt->Taille();
if (d_gijHH_tdt_3D.Taille() != ta_ex_d_gijHH_tdt)
{ d_gijHH_tdt_3D.Change_taille(ta_ex_d_gijHH_tdt); d_gijHH_tdt_3D_P.Change_taille(ta_ex_d_gijHH_tdt);
for (int i=1;i<=ta_ex_d_gijHH_tdt;i++) {d_gijHH_tdt_3D_P(i) = &(d_gijHH_tdt_3D(i));};
};
if( d_jacobien_tdt_3D.Taille() != ex.d_jacobien_tdt->Taille())
{d_jacobien_tdt_3D.Change_taille(ex.d_jacobien_tdt->Taille());}
// le tableau à double dimension
if (ex.d2_gijBB_tdt != NULL)
{ int taille1=ex.d2_gijBB_tdt->Taille1(); int taille2=ex.d2_gijBB_tdt->Taille2();
if (d2_gijBB_tdt_3D_P == NULL)
{ d2_gijBB_tdt_3D_P = impli_3D->d2_gijBB_tdt = new Tableau2<TenseurBB *> (taille1,taille2);
d2_gijBB_tdt_3D.Change_taille(taille1,taille2);
for (int i=1;i<=taille1;i++) for (int j=1;j<= taille2;j++)
{ (*d2_gijBB_tdt_3D_P)(i,j) = &d2_gijBB_tdt_3D(i,j);};
}
else if ((taille1 != d2_gijBB_tdt_3D.Taille1()) || (taille2 != d2_gijBB_tdt_3D.Taille2()))
{ d2_gijBB_tdt_3D_P->Change_taille(taille1,taille2);
d2_gijBB_tdt_3D.Change_taille(taille1,taille2);
for (int i=1;i<=taille1;i++) for (int j=1;j<= taille2;j++)
{ (*d2_gijBB_tdt_3D_P)(i,j) = &d2_gijBB_tdt_3D(i,j);};
};
};
// récup du conteneur spécifique
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
// si les grandeurs à 0 ont déjà été calculées, on les sauvegarde d'abord car elles vont
// être écrasées
// on commence par recopier les grandeurs de l'ordre 1 à 3
bool plusZero = true; // on complète avec des 0 dans un premier temps
int type_recopie = 0; // init: on recopie les grandeurs de 0,t et tdt
if ( save_resul.meti_ref_00.giB_ != NULL)
type_recopie = 1; // = 1 -> on transfert les grandeurs à t et tdt
impli_3D->Passage_de_Ordre1_vers_Ordre3(ex,plusZero,type_recopie);
#ifdef MISE_AU_POINT
if (Permet_affichage() > 8)
{ cout << "\n LoiContraintesPlanesDouble::Passage_metrique_ordre1_vers_3("
<< "\n 1) base giB_tdt_3D(1) "<< giB_tdt_3D(1)
<< "\n base giB_tdt_3D(2) "<< giB_tdt_3D(2)
<< "\n base giB_tdt_3D(2) "<< giB_tdt_3D(3) << flush;
};
#endif
// maintenant on s'occupe de mettre à jour les grandeurs manquantes
// - les bases naturelles: les vecteurs normaux dans les directions 2 et 3 sont normées
// et sont identiques pour les bases naturelles et duales
// Dans le cas d'un comportement isotrope, on choisit un repère par défaut
// si les repères associés ne sont pas définis on les définit
if ( save_resul.meti_ref_00.giB_ == NULL)
{ Deformation::Stmet inter(3,3); // création d'un conteneur bien dimensionné
save_resul.meti_ref_00 = save_resul.meti_ref_t = save_resul.meti_ref_tdt = inter;
// on s'occupe de la création du premier repère
CoordonneeB & giB0_1 = (*(save_resul.meti_ref_00.giB_)).CoordoB(1); // pour simplifier
giB0_1 = giB_0_3D(1); // on récupère la direction du premier vecteur
giB0_1.Normer(); // on le norme
// on crée 2 vecteurs normés, perpendiculaires à ce premier vecteur
CoordonneeB & giB0_2 = (*(save_resul.meti_ref_00.giB_)).CoordoB(2); // pour simplifier
CoordonneeB & giB0_3 = (*(save_resul.meti_ref_00.giB_)).CoordoB(3); // pour simplifier
// deux cas suivant que l'on dispose d'une orientation supplémentaire ou non
if (ex.giB_0->NbVecteur() > 1)
{ // a priori le second vecteur de la base est un vecteur d'un plan de ref
giB0_3 = (Util::ProdVec_coorB(giB_0_3D(1), (*(ex.giB_0))(2))).Normer();
giB0_2 = Util::ProdVec_coorB(giB0_3,giB0_1);
}
else
{ // cas où on n'a pas d'indication sup: -> on definit des vecteurs a priori
MathUtil2::Def_vecteurs_plan(giB0_1,giB0_2,giB0_3);
};
////--- debug
// cout << "\n debug LoiContraintesPlanesDouble::Passage_metrique_ordre1_vers_3(";
// cout << "\n giB0_1: "; giB0_1.Affiche();
// cout << "\n giB0_2: "; giB0_2.Affiche();
// cout << "\n giB0_3: "; giB0_3.Affiche();
//// --- fin debug
giB_0_3D.CoordoB(2)=giB0_2;giB_0_3D.CoordoB(3)=giB0_3;
giH_0_3D.CoordoH(2) = giB_0_3D(2).Bas_haut();
giH_0_3D.CoordoH(3) = giB_0_3D(3).Bas_haut();
};
////--- debug
// cout << "\n debug LoiContraintesPlanesDouble::Passage_metrique_ordre1_vers_3(";
// cout << "\n giH_0(2): "; giH_0_3D(2).Affiche();
// cout << "\n giH_0(3): "; giH_0_3D(3).Affiche();
//// --- fin debug
// on met à jour les deux autres directions
// *** pour l'instant on fait la même chose que pour t= 0, pour les deux autres temps
// *** , mais il faudra surement améliorer
CoordonneeB & giBt_1 = (*(save_resul.meti_ref_t.giB_)).CoordoB(1); // pour simplifier
giBt_1 = giB_t_3D(1); // on récupère la direction du premier vecteur
giBt_1.Normer(); // on le norme
// on crée 2 vecteurs normés, perpendiculaires à ce premier vecteur
CoordonneeB & giBt_2 = (*(save_resul.meti_ref_t.giB_)).CoordoB(2); // pour simplifier
CoordonneeB & giBt_3 = (*(save_resul.meti_ref_t.giB_)).CoordoB(3); // pour simplifier
// deux cas suivant que l'on dispose d'une orientation supplémentaire ou non
if (ex.giB_t->NbVecteur() > 1)
{ // a priori le second vecteur de la base est un vecteur d'un plan de ref
giBt_3 = (Util::ProdVec_coorB(giB_t_3D(1), (*(ex.giB_t))(2))).Normer();
giBt_2 = Util::ProdVec_coorB(giBt_3,giBt_1);
}
else
{ // cas où on n'a pas d'indication sup: -> on definit des vecteurs a priori
MathUtil2::Def_vecteurs_plan(giBt_1,giBt_2,giBt_3);
};
giB_t_3D.CoordoB(2)=giBt_2;giB_t_3D.CoordoB(3)=giBt_3;
giH_t_3D.CoordoH(2) = giB_t_3D(2).Bas_haut();
giH_t_3D.CoordoH(3) = giB_t_3D(3).Bas_haut();
#ifdef MISE_AU_POINT
if (Permet_affichage() > 8)
{ cout << "\n 2) base giB_tdt_3D(1) "<< giB_tdt_3D(1)
<< "\n base giB_tdt_3D(2) "<< giB_tdt_3D(2)
<< "\n base giB_tdt_3D(2) "<< giB_tdt_3D(3) << flush;
};
#endif
CoordonneeB & giBtdt_1 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(1); // pour simplifier
giBtdt_1 = giB_tdt_3D(1); // on récupère la direction du premier vecteur
giBtdt_1.Normer(); // on le norme
// on crée 2 vecteurs normés, perpendiculaires à ce premier vecteur
CoordonneeB & giBtdt_2 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(2); // pour simplifier
CoordonneeB & giBtdt_3 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(3); // pour simplifier
// deux cas suivant que l'on dispose d'une orientation supplémentaire ou non
if (ex.giB_tdt->NbVecteur() > 1)
{ // a priori le second vecteur de la base est un vecteur d'un plan de ref
giBtdt_3 = (Util::ProdVec_coorB(giB_tdt_3D(1), (*(ex.giB_tdt))(2))).Normer();
giBtdt_2 = Util::ProdVec_coorB(giBtdt_3,giBtdt_1);
}
else
{ // cas où on n'a pas d'indication sup: -> on definit des vecteurs a priori
MathUtil2::Def_vecteurs_plan(giBtdt_1,giBtdt_2,giBtdt_3);
};
giB_tdt_3D.CoordoB(2)=giBtdt_2;giB_tdt_3D.CoordoB(3)=giBtdt_3;
giH_tdt_3D.CoordoH(2) = giB_tdt_3D(2).Bas_haut();
giH_tdt_3D.CoordoH(3) = giB_tdt_3D(3).Bas_haut();
#ifdef MISE_AU_POINT
if (Permet_affichage() > 8)
{ cout << "\n 3) base giB_tdt_3D(1) "<< giB_tdt_3D(1)
<< "\n base giB_tdt_3D(2) "<< giB_tdt_3D(2)
<< "\n base giB_tdt_3D(2) "<< giB_tdt_3D(3) << flush;
};
#endif
// ici on calcul les vecteurs
// >>>> et la première partie de leurs variations, fonction du premier vecteur uniquement
// - les tenseurs métriques: a priori 1 pour les directions 2 et 3, initialement
gijBB_0_3D.Coor(2,2)=gijHH_0_3D.Coor(2,2)=gijBB_t_3D.Coor(2,2)=gijHH_t_3D.Coor(2,2)
=gijBB_tdt_3D.Coor(2,2)=gijHH_tdt_3D.Coor(2,2)=1.;
gijBB_0_3D.Coor(3,3)=gijHH_0_3D.Coor(3,3)=gijBB_t_3D.Coor(3,3)=gijHH_t_3D.Coor(3,3)
=gijBB_tdt_3D.Coor(3,3)=gijHH_tdt_3D.Coor(3,3)=1.;
// lors d'un premier appel, les variables save_resul.h_tsurh0 et save_resul.hsurh0 et idem pour l'épaisseur
// ne sont pas affectés car elles sont issues d'un premier calcul, dans ce cas on les mets arbitrairement à 1
// ensuite elles évolueront, idem pour la largeur
bool premier_passage = false;
if (save_resul.h_tsurh0 == 0.)
{ save_resul.h_tsurh0 = 1.; save_resul.hsurh0 = 1.;
save_resul.b_tsurb0 = 1.; save_resul.bsurb0 = 1.;
premier_passage = true;
};
if (save_resul.hsurh0 == 0.) // peut arriver dans le cas où on a une erreur dès le début
{ save_resul.hsurh0 = save_resul.h_tsurh0;}; // on réinitialise
if (save_resul.bsurb0 == 0.) // peut arriver dans le cas où on a une erreur dès le début
{ save_resul.bsurb0 = save_resul.b_tsurb0;}; // on réinitialise
// on calcul les vecteur 2 et 3 naturels et duals
giB_t_3D.CoordoB(2) *= save_resul.b_tsurb0; giH_t_3D.CoordoH(2) /= save_resul.b_tsurb0;
giB_tdt_3D.CoordoB(2) *= save_resul.bsurb0; giH_tdt_3D.CoordoH(2) /= save_resul.bsurb0;
giB_t_3D.CoordoB(3) *= save_resul.h_tsurh0; giH_t_3D.CoordoH(3) /= save_resul.h_tsurh0;
giB_tdt_3D.CoordoB(3) *= save_resul.hsurh0; giH_tdt_3D.CoordoH(3) /= save_resul.hsurh0;
#ifdef MISE_AU_POINT
if (Permet_affichage() > 8)
{ cout << "\n 4) base giB_tdt_3D(1) "<< giB_tdt_3D(1)
<< "\n base giB_tdt_3D(2) "<< giB_tdt_3D(2)
<< "\n base giB_tdt_3D(2) "<< giB_tdt_3D(3) << flush;
};
#endif
// - et les tenseurs métriques à t et tdt, comme les vecteurs en épaisseurs et largeurs
// sont normées initialement, seule l'intensité
// des normales est à mettre à jour compte tenue de la variation d'épaisseurs
double h_tsurh0_carre = save_resul.h_tsurh0 * save_resul.h_tsurh0;
gijBB_t_3D.Coor(3,3)=h_tsurh0_carre; gijHH_t_3D.Coor(3,3)= 1./h_tsurh0_carre;
double hsurh0_carre = save_resul.hsurh0 * save_resul.hsurh0;
gijBB_tdt_3D.Coor(3,3) = hsurh0_carre;gijHH_tdt_3D.Coor(3,3) = 1. / hsurh0_carre;
double b_tsurb0_carre = save_resul.b_tsurb0 * save_resul.b_tsurb0;
gijBB_t_3D.Coor(2,2)=b_tsurb0_carre; gijHH_t_3D.Coor(2,2)= 1./b_tsurb0_carre;
double bsurb0_carre = save_resul.bsurb0 * save_resul.bsurb0;
gijBB_tdt_3D.Coor(2,2) = bsurb0_carre;gijHH_tdt_3D.Coor(2,2) = 1. / bsurb0_carre;
// -- cas des variations
if (save_resul.d_hsurh0.Taille() != ta_ex_d_giB_tdt)
{save_resul.d_hsurh0.Change_taille(ta_ex_d_giB_tdt); // donc à 0 la première fois
save_resul.d_bsurb0.Change_taille(ta_ex_d_giB_tdt);
};
// - les variations du vecteur normal
#ifdef MISE_AU_POINT
if (ex.d_giB_tdt->Taille() != ex.d_giH_tdt->Taille())
{ cout << "\n **** sans doute une erreur: (ex.d_giB_tdt->Taille() "<< ex.d_giB_tdt->Taille()
<< " != ex.d_giH_tdt->Taille()) " << ex.d_giH_tdt->Taille()
<< "\n LoiDeformationsPlanes::Passage_metrique_ordre1_vers_3(const Met_abstraite::Impli& ex)";
Sortie(1);
};
#endif
// maintenant on s'occupe de la seconde partie des variations
// >>>> fonction du changement d'épaisseur et de largeur
if (!premier_passage)
{Vecteur& d_hsurh0 = save_resul.d_hsurh0;
Vecteur& d_bsurb0 = save_resul.d_bsurb0;
for (int iddl=1;iddl<= ta_ex_d_giB_tdt;iddl++)
{ double d_hsurh0_iddl = d_hsurh0(iddl);
// gib(3) = h/h0 * N -> d_gib(3) = d (h/h0 * N) = d_(h/h0) * N + h/h0 * d_N
// actuellement, c'est la partie d_N qui est stockée dans d_gib(3), on rajoute donc la partie restante
d_giB_tdt_3D(iddl).CoordoB(3) *= save_resul.hsurh0;
d_giB_tdt_3D(iddl).CoordoB(3) += d_hsurh0_iddl * giB_normer_3_tdt_3D_sauve; // base naturelle
// gih(3) = 1/giB(3) -> d_gih(3)(3) = - d_giB(3)(3) / (giB(3)*giB(3))
d_giH_tdt_3D(iddl).CoordoH(3) = (d_giB_tdt_3D(iddl)(3) / (-hsurh0_carre)).Bas_haut() ; // base duale
// gijBB(3)= giB(3) * giB(3) -> d_gijBB(3)= 2. * d_giB(3) * giB(3)
d_gijBB_tdt_3D(iddl).Coor(3,3) = 2. * (d_giB_tdt_3D(iddl)(3).ScalBB(giB_tdt_3D(3)));
// gijHH(3)= giH(3) * giH(3) -> d_gijHH(3)= 2. * d_giH(3) * giH(3)
d_gijHH_tdt_3D (iddl).Coor(3,3) = 2. * (d_giH_tdt_3D(iddl)(3).ScalHH(giH_tdt_3D(3)));
double d_bsurb0_iddl = d_bsurb0(iddl);
// gib(2) = b/b0 * N -> d_gib(2) = d (b/b0 * N) = d_(b/b0) * N + b/b0 * d_N
// actuellement, c'est la partie d_N qui est stockée dans d_gib(3), on rajoute donc la partie restante
d_giB_tdt_3D(iddl).CoordoB(2) *= save_resul.bsurb0;
d_giB_tdt_3D(iddl).CoordoB(2) += d_bsurb0_iddl * giB_normer_3_tdt_3D_sauve; // base naturelle
// gih(2) = 1/giB(2) -> d_gih(2)(2) = - d_giB(3)(3) / (giB(3)*giB(3))
d_giH_tdt_3D(iddl).CoordoH(2) = (d_giB_tdt_3D(iddl)(2) / (-bsurb0_carre)).Bas_haut() ; // base duale
// gijBB(2)= giB(2) * giB(2) -> d_gijBB(2)= 2. * d_giB(2) * giB(2)
d_gijBB_tdt_3D(iddl).Coor(2,2) = 2. * (d_giB_tdt_3D(iddl)(2).ScalBB(giB_tdt_3D(2)));
// gijHH(2)= giH(2) * giH(2) -> d_gijHH(2)= 2. * d_giH(2) * giH(2)
d_gijHH_tdt_3D (iddl).Coor(2,2) = 2. * (d_giH_tdt_3D(iddl)(2).ScalHH(giH_tdt_3D(2)));
};
};
};
// passage des grandeurs métriques de l'ordre 1 à 3: cas explicite
void LoiContraintesPlanesDouble::Passage_metrique_ordre1_vers_3(const Met_abstraite::Expli_t_tdt& ex)
{ // on s'occupe du redimensionnement éventuel
// on affecte et/ou on redimensionne éventuellement les tableaux fonction du nombre de ddl pour le passage 3D
// la partie dépendant des vitesses: entre accolades pour pouvoir fermer
{if (ex.gradVmoyBB_t != NULL) {expli_3D->gradVmoyBB_t= gradVmoyBB_t_3D_P = &gradVmoyBB_t_3D;};
if (ex.gradVmoyBB_tdt != NULL) {expli_3D->gradVmoyBB_tdt=gradVmoyBB_tdt_3D_P = &gradVmoyBB_tdt_3D;};
if (ex.gradVBB_tdt != NULL) {expli_3D->gradVBB_tdt=gradVBB_tdt_3D_P = &gradVBB_tdt_3D;};
}; // fin de la partie dédiée à la vitesse
// la dimension des tableaux : on considère que tous les tableaux doivent avoir la même dimension
int ta_ex_d_gijBB_tdt = ex.d_gijBB_tdt->Taille(); // la dimension
// maintenant on s'occupe uniquement du redimensionnement des tableaux restants
int ta_d_gijBB_tdt_3D = d_gijBB_tdt_3D.Taille();
if (ta_d_gijBB_tdt_3D != ta_ex_d_gijBB_tdt)
{// cela veut dire que tous les tableaux sont mal dimensionnés
// -- on s'occupe des tableaux nécessaire à la métrique
d_gijBB_tdt_3D.Change_taille(ta_ex_d_gijBB_tdt); d_gijBB_tdt_3D_P.Change_taille(ta_ex_d_gijBB_tdt);
for (int i=1;i<=ta_ex_d_gijBB_tdt;i++) {d_gijBB_tdt_3D_P(i) = &(d_gijBB_tdt_3D(i));};
};
// récup du conteneur spécifique
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
// on commence par recopier les grandeurs de l'ordre 1 à 3
bool plusZero = true; // on complète avec des 0 dans un premier temps
int type_recopie = 0; // init: on recopie les grandeurs de 0,t et tdt
if ( save_resul.meti_ref_00.giB_ != NULL)
type_recopie = 1; // = 1 -> on transfert les grandeurs à t et tdt
expli_3D->Passage_de_Ordre1_vers_Ordre3(ex,plusZero,type_recopie);
// maintenant on s'occupe de mettre à jour les grandeurs manquantes
// - les bases naturelles: les vecteurs normaux dans les directions 2 et 3 sont normées
// et sont identiques pour les bases naturelles et duales
// Dans le cas d'un comportement isotrope, on choisit un repère par défaut
// si les repères associés ne sont définis on les définit
if ( save_resul.meti_ref_00.giB_ == NULL)
{ Deformation::Stmet inter(3,3); // création d'un conteneur bien dimensionné
save_resul.meti_ref_00 = save_resul.meti_ref_t = save_resul.meti_ref_tdt = inter;
// on s'occupe de la création du premier repère
CoordonneeB & giB0_1 = (*(save_resul.meti_ref_00.giB_)).CoordoB(1); // pour simplifier
giB0_1 = giB_0_3D(1); // on récupère la direction du premier vecteur
giB0_1.Normer(); // on le norme
// on crée 2 vecteurs normés, perpendiculaires à ce premier vecteur
CoordonneeB & giB0_2 = (*(save_resul.meti_ref_00.giB_)).CoordoB(2); // pour simplifier
CoordonneeB & giB0_3 = (*(save_resul.meti_ref_00.giB_)).CoordoB(3); // pour simplifier
// Def_vecteurs_plan(giB0_1,giB0_2,giB0_3);
giB_0_3D.CoordoB(2)=giB0_2;giB_0_3D.CoordoB(3)=giB0_3;
};
// on met à jour les deux autres directions
// *** pour l'instant on fait la même chose que pour t= 0, pour les deux autres temps
// *** , mais il faudra surement améliorer
CoordonneeB & giBt_1 = (*(save_resul.meti_ref_t.giB_)).CoordoB(1); // pour simplifier
giBt_1 = giB_t_3D(1); // on récupère la direction du premier vecteur
giBt_1.Normer(); // on le norme
// on crée 2 vecteurs normés, perpendiculaires à ce premier vecteur
CoordonneeB & giBt_2 = (*(save_resul.meti_ref_t.giB_)).CoordoB(2); // pour simplifier
CoordonneeB & giBt_3 = (*(save_resul.meti_ref_t.giB_)).CoordoB(3); // pour simplifier
// Def_vecteurs_plan(giBt_1,giBt_2,giBt_3);
giB_t_3D.CoordoB(2)=giBt_2;giB_t_3D.CoordoB(3)=giBt_3;
CoordonneeB & giBtdt_1 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(1); // pour simplifier
giBtdt_1 = giB_tdt_3D(1); // on récupère la direction du premier vecteur
giBtdt_1.Normer(); // on le norme
// on crée 2 vecteurs normés, perpendiculaires à ce premier vecteur
CoordonneeB & giBtdt_2 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(2); // pour simplifier
CoordonneeB & giBtdt_3 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(3); // pour simplifier
// Def_vecteurs_plan(giBtdt_1,giBtdt_2,giBtdt_3);
giB_tdt_3D.CoordoB(2)=giBtdt_2;giB_tdt_3D.CoordoB(3)=giBtdt_3;
// - les tenseurs métriques: a priori 1 pour les directions 2 et 3, initialement
gijBB_0_3D.Coor(2,2)=gijHH_0_3D.Coor(2,2)=gijBB_t_3D.Coor(2,2)=gijHH_t_3D.Coor(2,2)
=gijBB_tdt_3D.Coor(2,2)=gijHH_tdt_3D.Coor(2,2)=1.;
gijBB_0_3D.Coor(3,3)=gijHH_0_3D.Coor(3,3)=gijBB_t_3D.Coor(3,3)=gijHH_t_3D.Coor(3,3)
=gijBB_tdt_3D.Coor(3,3)=gijHH_tdt_3D.Coor(3,3)=1.;
// lors d'un premier appel, les variables save_resul.h_tsurh0 et save_resul.hsurh0
// ne sont pas affecté car elles sont issues d'un premier calcul, dans ce cas on les mets arbitrairement à 1
// ensuite elles évolueront, idem pour la largeur
bool premier_passage = false;
if (save_resul.h_tsurh0 == 0.)
{ save_resul.h_tsurh0 = 1.; save_resul.hsurh0 = 1.;
save_resul.b_tsurb0 = 1.; save_resul.bsurb0 = 1.;
premier_passage = true;
};
if (save_resul.hsurh0 == 0.) // peut arriver dans le cas où on a une erreur dès le début
{ save_resul.hsurh0 = save_resul.h_tsurh0;}; // on réinitialise
if (save_resul.bsurb0 == 0.) // peut arriver dans le cas où on a une erreur dès le début
{ save_resul.bsurb0 = save_resul.b_tsurb0;}; // on réinitialise
// on calcul les vecteur 2 et 3 naturels et duals
giB_t_3D.CoordoB(2) *= save_resul.b_tsurb0; giH_t_3D.CoordoH(2) /= save_resul.b_tsurb0;
giB_tdt_3D.CoordoB(2) *= save_resul.bsurb0; giH_tdt_3D.CoordoH(2) /= save_resul.bsurb0;
giB_t_3D.CoordoB(3) *= save_resul.h_tsurh0; giH_t_3D.CoordoH(3) /= save_resul.h_tsurh0;
giB_tdt_3D.CoordoB(3) *= save_resul.hsurh0; giH_tdt_3D.CoordoH(3) /= save_resul.hsurh0;
// - et les tenseurs métriques à t et tdt, comme les vecteurs en épaisseurs et largeurs
// sont normées initialement, seule l'intensité
// des normales est à mettre à jour compte tenue de la variation d'épaisseurs
double b_tsurb0_carre = save_resul.b_tsurb0 * save_resul.b_tsurb0;
gijBB_t_3D.Coor(2,2)=b_tsurb0_carre; gijHH_t_3D.Coor(2,2)= 1./b_tsurb0_carre;
double bsurb0_carre = save_resul.bsurb0 * save_resul.bsurb0;
gijBB_tdt_3D.Coor(2,2) = bsurb0_carre;gijHH_tdt_3D.Coor(2,2) = 1. / bsurb0_carre;
double h_tsurh0_carre = save_resul.h_tsurh0 * save_resul.h_tsurh0;
gijBB_t_3D.Coor(3,3)=h_tsurh0_carre; gijHH_t_3D.Coor(3,3)= 1./h_tsurh0_carre;
double hsurh0_carre = save_resul.hsurh0 * save_resul.hsurh0;
gijBB_tdt_3D.Coor(3,3) = hsurh0_carre;gijHH_tdt_3D.Coor(3,3) = 1. / hsurh0_carre;
};
// passage des grandeurs métriques de l'ordre 1 à 3: cas Umat
void LoiContraintesPlanesDouble::Passage_metrique_ordre1_vers_3(const Met_abstraite::Umat_cont& ex)
{ // on s'occupe du redimensionnement éventuel
// on affecte et/ou on redimensionne éventuellement les tableaux fonction du nombre de ddl pour le passage 3D
// la partie dépendant des vitesses: entre accolades pour pouvoir fermer
{if (ex.gradVmoyBB_t != NULL) {umat_cont_3D->gradVmoyBB_t= gradVmoyBB_t_3D_P = &gradVmoyBB_t_3D;};
if (ex.gradVmoyBB_tdt != NULL) {umat_cont_3D->gradVmoyBB_tdt=gradVmoyBB_tdt_3D_P = &gradVmoyBB_tdt_3D;};
if (ex.gradVBB_tdt != NULL) {umat_cont_3D->gradVBB_tdt=gradVBB_tdt_3D_P = &gradVBB_tdt_3D;};
}; // fin de la partie dédiée à la vitesse
// récup du conteneur spécifique
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
// on commence par recopier les grandeurs de l'ordre 1 à 3
bool plusZero = true; // on complète avec des 0 dans un premier temps
int type_recopie = 0; // init: on recopie les grandeurs de 0,t et tdt
if ( save_resul.meti_ref_00.giB_ != NULL)
type_recopie = 1; // = 1 -> on transfert les grandeurs à t et tdt
umat_cont_3D->Passage_de_Ordre1_vers_Ordre3(ex,plusZero,type_recopie);
// maintenant on s'occupe de mettre à jour les grandeurs manquantes
// - les bases naturelles: les vecteurs normaux dans les directions 2 et 3 sont normées
// et sont identiques pour les bases naturelles et duales
// Dans le cas d'un comportement isotrope, on choisit un repère par défaut
// si les repères associés ne sont définis on les définit
if ( save_resul.meti_ref_00.giB_ == NULL)
{ Deformation::Stmet inter(3,3); // création d'un conteneur bien dimensionné
save_resul.meti_ref_00 = save_resul.meti_ref_t = save_resul.meti_ref_tdt = inter;
// on s'occupe de la création du premier repère
CoordonneeB & giB0_1 = (*(save_resul.meti_ref_00.giB_)).CoordoB(1); // pour simplifier
giB0_1 = giB_0_3D(1); // on récupère la direction du premier vecteur
giB0_1.Normer(); // on le norme
// on crée 2 vecteurs normés, perpendiculaires à ce premier vecteur
CoordonneeB & giB0_2 = (*(save_resul.meti_ref_00.giB_)).CoordoB(2); // pour simplifier
CoordonneeB & giB0_3 = (*(save_resul.meti_ref_00.giB_)).CoordoB(3); // pour simplifier
MathUtil2::Def_vecteurs_plan(giB0_1,giB0_2,giB0_3);
// comme la base est orthonormée, les giH sont identiques aux giB
(*(save_resul.meti_ref_00.giH_)).Affectation_trans_variance((*(save_resul.meti_ref_00.giB_)));
// on affecte les vecteurs qui seront transmis via l'umat
giB_0_3D.CoordoB(2)=giB0_2;giB_0_3D.CoordoB(3)=giB0_3;
BaseH& a0H = (*(save_resul.meti_ref_00.giH_)); // pour simplifier
giH_0_3D.CoordoH(1)=a0H(1);giH_0_3D.CoordoH(2)=a0H(2);giH_0_3D.CoordoH(3)=a0H(3);
};
////--- debug
// cout << "\n debug LoiContraintesPlanesDouble::Passage_metrique_ordre1_vers_3(";
// cout << "\n giH_0(2): "; giH_0_3D(2).Affiche();
// cout << "\n giH_0(3): "; giH_0_3D(3).Affiche();
//// --- fin debug
// on met à jour les deux autres directions
// *** pour l'instant on fait la même chose que pour t= 0, pour les deux autres temps
// *** , mais il faudra surement améliorer
CoordonneeB & giBt_1 = (*(save_resul.meti_ref_t.giB_)).CoordoB(1); // pour simplifier
giBt_1 = giB_t_3D(1); // on récupère la direction du premier vecteur
giBt_1.Normer(); // on le norme
// on crée 2 vecteurs normés, perpendiculaires à ce premier vecteur
CoordonneeB & giBt_2 = (*(save_resul.meti_ref_t.giB_)).CoordoB(2); // pour simplifier
CoordonneeB & giBt_3 = (*(save_resul.meti_ref_t.giB_)).CoordoB(3); // pour simplifier
MathUtil2::Def_vecteurs_plan(giBt_1,giBt_2,giBt_3);
// comme la base est orthonormée, les giH sont identiques aux giB
// ensuite plus bas, les normes seront changées pour tenir compte des variations d'épaisseur et de largeur
(*(save_resul.meti_ref_t.giH_)).Affectation_trans_variance((*(save_resul.meti_ref_t.giB_)));
// on affecte les vecteurs qui seront transmis via l'umat
giB_t_3D.CoordoB(2)=giBt_2;giB_t_3D.CoordoB(3)=giBt_3;
BaseH& a_t_H = (*(save_resul.meti_ref_t.giH_)); // pour simplifier
giH_t_3D.CoordoH(1)=a_t_H(1);giH_t_3D.CoordoH(2)=a_t_H(2);giH_t_3D.CoordoH(3)=a_t_H(3);
//--- cas à tdt
CoordonneeB & giBtdt_1 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(1); // pour simplifier
giBtdt_1 = giB_tdt_3D(1); // on récupère la direction du premier vecteur
giBtdt_1.Normer(); // on le norme
// on crée 2 vecteurs normés, perpendiculaires à ce premier vecteur
CoordonneeB & giBtdt_2 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(2); // pour simplifier
CoordonneeB & giBtdt_3 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(3); // pour simplifier
MathUtil2::Def_vecteurs_plan(giBtdt_1,giBtdt_2,giBtdt_3);
giB_tdt_3D.CoordoB(2)=giBtdt_2;giB_tdt_3D.CoordoB(3)=giBtdt_3;
// comme la base est orthonormée ici, les giH sont identiques aux giB
// ensuite plus bas, les normes seront changées pour tenir compte des variations d'épaisseur et de largeur
(*(save_resul.meti_ref_tdt.giH_)).Affectation_trans_variance((*(save_resul.meti_ref_tdt.giB_)));
// on affecte les vecteurs qui seront transmis via l'umat
giB_tdt_3D.CoordoB(2)=giBtdt_2;giB_tdt_3D.CoordoB(3)=giBtdt_3;
BaseH& a_tdt_H = (*(save_resul.meti_ref_tdt.giH_)); // pour simplifier
giH_tdt_3D.CoordoH(1)=a_tdt_H(1);giH_tdt_3D.CoordoH(2)=a_tdt_H(2);giH_tdt_3D.CoordoH(3)=a_tdt_H(3);
// - les tenseurs métriques: a priori 1 pour les directions 2 et 3, initialement
gijBB_0_3D.Coor(2,2)=gijHH_0_3D.Coor(2,2)=gijBB_t_3D.Coor(2,2)=gijHH_t_3D.Coor(2,2)
=gijBB_tdt_3D.Coor(2,2)=gijHH_tdt_3D.Coor(2,2)=1.;
gijBB_0_3D.Coor(3,3)=gijHH_0_3D.Coor(3,3)=gijBB_t_3D.Coor(3,3)=gijHH_t_3D.Coor(3,3)
=gijBB_tdt_3D.Coor(3,3)=gijHH_tdt_3D.Coor(3,3)=1.;
// lors d'un premier appel, les variables save_resul.h_tsurh0 et save_resul.hsurh0
// ne sont pas affecté car elles sont issues d'un premier calcul, dans ce cas on les mets arbitrairement à 1
// ensuite elles évolueront, idem pour la largeur
bool premier_passage = false;
if (save_resul.h_tsurh0 == 0.)
{ save_resul.h_tsurh0 = 1.; save_resul.hsurh0 = 1.;
premier_passage = true;
};
if (save_resul.hsurh0 == 0.) // peut arriver dans le cas où on a une erreur dès le début
{ save_resul.hsurh0 = save_resul.h_tsurh0;}; // on réinitialise
// on peut avoir la même chose en b, indépendamment, lorsque l'on passe de contrainte plane
// à doublement plane
if (save_resul.b_tsurb0 == 0.)
{ save_resul.b_tsurb0 = 1.; save_resul.bsurb0 = 1.;
premier_passage = true;
};
if (save_resul.bsurb0 == 0.) // peut arriver dans le cas où on a une erreur dès le début
{ save_resul.bsurb0 = save_resul.b_tsurb0;}; // on réinitialise
// on calcul les vecteur 2 et 3 naturels et duals
giB_t_3D.CoordoB(2) *= save_resul.b_tsurb0; giH_t_3D.CoordoH(2) /= save_resul.b_tsurb0;
giB_tdt_3D.CoordoB(2) *= save_resul.bsurb0; giH_tdt_3D.CoordoH(2) /= save_resul.bsurb0;
giB_t_3D.CoordoB(3) *= save_resul.h_tsurh0; giH_t_3D.CoordoH(3) /= save_resul.h_tsurh0;
giB_tdt_3D.CoordoB(3) *= save_resul.hsurh0; giH_tdt_3D.CoordoH(3) /= save_resul.hsurh0;
// - et les tenseurs métriques à t et tdt, comme les vecteurs en épaisseurs et largeurs
// sont normées initialement, seule l'intensité
// des normales est à mettre à jour compte tenue de la variation d'épaisseurs
double b_tsurb0_carre = save_resul.b_tsurb0 * save_resul.b_tsurb0;
gijBB_t_3D.Coor(2,2)=b_tsurb0_carre; gijHH_t_3D.Coor(2,2)= 1./b_tsurb0_carre;
double bsurb0_carre = save_resul.bsurb0 * save_resul.bsurb0;
gijBB_tdt_3D.Coor(2,2) = bsurb0_carre;gijHH_tdt_3D.Coor(2,2) = 1. / bsurb0_carre;
double h_tsurh0_carre = save_resul.h_tsurh0 * save_resul.h_tsurh0;
gijBB_t_3D.Coor(3,3)=h_tsurh0_carre; gijHH_t_3D.Coor(3,3)= 1./h_tsurh0_carre;
double hsurh0_carre = save_resul.hsurh0 * save_resul.hsurh0;
gijBB_tdt_3D.Coor(3,3) = hsurh0_carre;gijHH_tdt_3D.Coor(3,3) = 1. / hsurh0_carre;
};
// passage des informations liées à la déformation de 1 vers 3, et variation de volume
// utilisée pour les méthodes 1 et 2
// entre autre: calcul du jacobien 3D en fonction de l'élongation 1D et des variations h/h0 et b/b0
void LoiContraintesPlanesDouble::Passage_deformation_volume_ordre1_vers_3
// (int cas, TenseurBB& DepsBB,TenseurBB & epsBB_tdt,Tableau <TenseurBB *>* d_epsBB
(TenseurBB& DepsBB,TenseurBB & epsBB_tdt,Tableau <TenseurBB *>* d_epsBB
,TenseurBB & delta_epsBB,const double& jacobien_0,double& jacobien
,Vecteur* d_jacobien_tdt,const double& jacobien_t)
{ // au début on complète avec des 0, puis dans un second temps on tient compte de la déformation d'épaisseur eps33
bool plusZero = true;
Deps_BB_3D.Affectation_trans_dimension(*((Tenseur1BB*) &DepsBB),plusZero);
eps_BB_3D.Affectation_trans_dimension(*((Tenseur1BB*) &epsBB_tdt),plusZero);
delta_eps_BB_3D.Affectation_trans_dimension(*((Tenseur1BB*) &delta_epsBB),plusZero);
// récup du conteneur spécifique
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
if(d_epsBB != NULL)
{ int taille = d_epsBB->Taille();
// redimensionnement éventuel
int ta_d_eps_BB_3D = d_eps_BB_3D.Taille();
if (ta_d_eps_BB_3D != taille)
{ // cela veut dire que tous les tableaux sont mal dimensionnés
d_eps_BB_3D.Change_taille(taille); d_eps_BB_3D_P.Change_taille(taille);
for (int i=1;i<=taille;i++) {d_eps_BB_3D_P(i) = &(d_eps_BB_3D(i));}
};
// passage 1D 3D
for (int i=1;i<=taille;i++)
{ d_eps_BB_3D(i).Affectation_trans_dimension(*((Tenseur1BB*) (*d_epsBB)(i)),plusZero);}
if (save_resul.d_hsurh0.Taille() != taille)
save_resul.d_hsurh0.Change_taille(taille); // donc à 0 la première fois
if (save_resul.d_bsurb0.Taille() != taille)
save_resul.d_bsurb0.Change_taille(taille); // donc à 0 la première fois
};
// -- calcul du jacobien 3D
/* la partie cas = 2: ne fonctionne pas bien le choix est fait
d'utiliser les variations d'épaisseur et de largeur
// dans le cas d'un calcul dans une direction quelconque
// on utilise la déformation mécanique sauvegardée pour calculer les jacobiens
if (cas == 2)
{ //récup de eps_mecaBB_t
Tenseur3BB& eps_mecaBB_t = *((Tenseur3BB*) (save_resul.eps_mecaBB_t));
// --- on s'occupe des jacobiens
// à t le jacobien est calculée à partir de la déformation supposée d'Almansi
// ici il s'agit du jacobien relatif au comportement mécanique
// |J| = sqrt(det(2 * eps_ij_t + g_ij(0)))
// calcul de la matrice correspondante à (2 * eps_ij_t + g_ij(0))
gij_meca_BB_t = (2. * eps_mecaBB_t + (*(umat_cont_3D->gijBB_0)));
gij_meca_BB_t.Matrice_composante(mat_inter);
// calcul du jacobien
jacobien_t_3D = sqrt(Dabs(mat_inter.Determinant()));
// on fait de même pour le jacobien initial ce qui permet d'être cohérent
// avec la métrique initiale et la déformation au sens d'Almansi
(umat_cont_3D->gijBB_0)->Matrice_composante(mat_inter);
jacobien_0_3D = sqrt(Dabs(mat_inter.Determinant()));
// le jacobien à tdt sera calculé pendant la résolution
}
else if (cas == 1)
*/
{// première phase on se contente de recopier, ce qui suppose une épaisseur = 1.
jacobien_0_3D = jacobien_0; // -> partie 1D
jacobien_tdt_3D = sauve_jacobien1D_tdt = jacobien;
// ici, l'épaisseur initiale h0 est toujours 1., donc h/h0 = également h et idem pour b/b0
jacobien_tdt_3D *= save_resul.hsurh0 * save_resul.bsurb0;
// même travail pour t
jacobien_t_3D = jacobien_t;
jacobien_t_3D *= save_resul.h_tsurh0 * save_resul.b_tsurb0;
};
// recup de l'incrément de temps
double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();
double unSurDeltat=0;
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
{ // non 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 des def, delta def et D
if (!calcul_en_3D_via_direction_quelconque)
{ //dans le cas classique, on dispose de la variation d'épaisseur et de largeur pour déterminer les def
switch (type_de_deformation)
{case DEFORMATION_STANDART : case DEFORMATION_POUTRE_PLAQUE_STANDART :
// cas d'une déformation d'Almansi
{ // epsBB33 = 1/2 * (1. - (h0/h)^2), en orthonormee
// dans le repère local: epsBB33 = 1/2 * (h^2 - 1.), or h0=1. donc : epsBB33 = 1/2 * ((h/h0)^2 - 1.)
eps_BB_3D.Coor(3,3) = 0.5 * (save_resul.hsurh0 * save_resul.hsurh0 - 1.);
delta_eps_BB_3D.Coor(3,3) = 0.5 * (save_resul.hsurh0 * save_resul.hsurh0 - save_resul.h_tsurh0 * save_resul.h_tsurh0);
// idem pour la largeur
eps_BB_3D.Coor(2,2) = 0.5 * (save_resul.bsurb0 * save_resul.bsurb0 - 1.);
delta_eps_BB_3D.Coor(2,2) = 0.5 * (save_resul.bsurb0 * save_resul.bsurb0 - save_resul.b_tsurb0 * save_resul.b_tsurb0);
Deps_BB_3D.Coor(3,3) = delta_eps_BB_3D(3,3)* unSurDeltat;
Deps_BB_3D.Coor(2,2) = delta_eps_BB_3D(2,2)* unSurDeltat;
// dans le cas où on calcule les variations, on intègre également les variations de eps33 et eps22
if (d_epsBB != NULL)
{ Vecteur& d_hsurh0 = save_resul.d_hsurh0;
Vecteur& d_bsurb0 = save_resul.d_bsurb0;
int taille = d_epsBB->Taille();
for (int i=1;i<=taille;i++)
{d_eps_BB_3D(i).Coor(3,3)= save_resul.hsurh0 * d_hsurh0(i);
d_eps_BB_3D(i).Coor(2,2)= save_resul.bsurb0 * d_bsurb0(i);
}
};
};
break;
case DEFORMATION_LOGARITHMIQUE : case DEF_CUMUL_CORROTATIONNEL : case DEFORMATION_CUMU_LOGARITHMIQUE :
// cas d'une def logarithmique ou une approximation
{ eps_BB_3D.Coor(3,3) = log(save_resul.hsurh0);
delta_eps_BB_3D.Coor(3,3) = log(save_resul.hsurh0) - log(save_resul.h_tsurh0);
// idem pour la largeur
eps_BB_3D.Coor(2,2) = log(save_resul.bsurb0);
delta_eps_BB_3D.Coor(2,2) = log(save_resul.bsurb0) - log(save_resul.b_tsurb0);
Deps_BB_3D.Coor(3,3) = delta_eps_BB_3D(3,3)* unSurDeltat;
Deps_BB_3D.Coor(2,2) = delta_eps_BB_3D(2,2)* unSurDeltat;
// dans le cas où on calcule les variations, on intègre également les variations de eps33 et eps22
if (d_epsBB != NULL)
{ Vecteur& d_hsurh0 = save_resul.d_hsurh0;
Vecteur& d_bsurb0 = save_resul.d_bsurb0;
int taille = d_epsBB->Taille();
// d_eps_BB_3D(3,3) = d_hsurh0 * exp(save_resul.hsurh0)
for (int i=1;i<=taille;i++)
{d_eps_BB_3D(i).Coor(3,3)= d_hsurh0(i) * eps_BB_3D(3,3);
d_eps_BB_3D(i).Coor(2,2)= d_bsurb0(i) * eps_BB_3D(2,2);
}
};
};
break;
default :
cout << "\nErreur : type de deformation qui n'est pas actuellement pris en compte, type= "
<< Nom_type_deformation(type_de_deformation);
cout << "\n LoiContraintesPlanesDouble::Passage_deformation_volume_ordre1_vers_3 \n";
Sortie(1);
};
}
else
// dans le cas où on travaille dans le repère gi en 3D on dispose directement des infos de déformation
{//récup de eps_mecaBB
Tenseur3BB& eps_mecaBB_t = *((Tenseur3BB*) (save_resul.eps_mecaBB_t));
Tenseur3BB& eps_mecaBB = *((Tenseur3BB*) (save_resul.eps_mecaBB));
// on sauve ce qu'il y a dans la direction 11 pour éviter les erreurs d'arrondis
double eps11= (*ptintmeca.EpsBB())(1,1);
double Deps11= (*ptintmeca.DepsBB())(1,1);
double delta_eps11= (*ptintmeca.DeltaEpsBB())(1,1);
// les manips entre tenseurs
(*ptintmeca.EpsBB()) = eps_mecaBB;
(*ptintmeca.DeltaEpsBB()) = eps_mecaBB - eps_mecaBB_t;
(*ptintmeca.DepsBB()) = (*ptintmeca.DeltaEpsBB()) * unSurDeltat;
// récup des valeurs sauvegardées (normalement ne devrait pas changer...)
(*ptintmeca.EpsBB()).Coor(1,1) = eps11;
(*ptintmeca.DeltaEpsBB()).Coor(1,1) = delta_eps11;
(*ptintmeca.DepsBB()).Coor(1,1) = Deps11;
};
if(d_jacobien_tdt != NULL)
{ Vecteur& d_hsurh0 = save_resul.d_hsurh0;
Vecteur& d_bsurb0 = save_resul.d_bsurb0;
// jacobien_tdt_3D = hsurh0 * jacobien_1D
// --> d_jacobien_tdt_3D = d_hsurh0 * jacobien_1D + hsurh0 * d_jacobien_1D
d_jacobien_tdt_3D = d_hsurh0 * save_resul.bsurb0 * jacobien
+ save_resul.hsurh0 * d_bsurb0 * jacobien
+ save_resul.hsurh0 * save_resul.bsurb0 * (*d_jacobien_tdt);
};
//---debug
//cout << "\n LoiContraintesPlanesDouble::Passage_deformation_volume_ordre1_vers_3, eps_BB_3D(3,3) = "
// << eps_BB_3D(3,3) << " jacobien_tdt_3D= " << jacobien_tdt_3D;
// --- fin debug
};
// mise à jour des informations liées à la déformation de 1 vers 3
// uniquement méthode 1
void LoiContraintesPlanesDouble::Mise_a_jour_deformations_et_Jacobien_en_3D
()
{ // récup du conteneur spécifique
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
// -- calcul du jacobien 3D
// première phase on se contente de recopier, ce qui suppose une épaisseur = 1. et une largeur = 1.
jacobien_tdt_3D = sauve_jacobien1D_tdt * save_resul.hsurh0 * save_resul.bsurb0;
// -- maintenant on tient compte de la déformation d'épaisseur
switch (type_de_deformation)
{case DEFORMATION_STANDART : case DEFORMATION_POUTRE_PLAQUE_STANDART :
// cas d'une déformation d'Almansi
{ // epsBB33 = 1/2 * (1. - (h0/h)^2), en orthonormee
// dans le repère local: epsBB33 = 1/2 * (h^2 - 1.), or h0=1. donc : epsBB33 = 1/2 * ((h/h0)^2 - 1.)
eps_BB_3D.Coor(3,3) = 0.5 * (save_resul.hsurh0 * save_resul.hsurh0 - 1.);
delta_eps_BB_3D.Coor(3,3) = 0.5 * (save_resul.hsurh0 * save_resul.hsurh0 - save_resul.h_tsurh0 * save_resul.h_tsurh0);
// idem pour la direction 2
eps_BB_3D.Coor(2,2) = 0.5 * (save_resul.bsurb0 * save_resul.bsurb0 - 1.);
delta_eps_BB_3D.Coor(2,2) = 0.5 * (save_resul.bsurb0 * save_resul.bsurb0 - save_resul.b_tsurb0 * save_resul.b_tsurb0);
};
break;
case DEFORMATION_LOGARITHMIQUE : case DEF_CUMUL_CORROTATIONNEL : case DEFORMATION_CUMU_LOGARITHMIQUE :
// cas d'une def logarithmique ou une approximation
{ eps_BB_3D.Coor(3,3) = log(save_resul.hsurh0);
delta_eps_BB_3D.Coor(3,3) = log(save_resul.hsurh0) - log(save_resul.h_tsurh0);
// idem pour la direction 2
eps_BB_3D.Coor(2,2) = log(save_resul.bsurb0);
delta_eps_BB_3D.Coor(2,2) = log(save_resul.bsurb0) - log(save_resul.b_tsurb0);
};
break;
default :
cout << "\nErreur : type de deformation qui n'est pas actuellement pris en compte, type= "
<< Nom_type_deformation(type_de_deformation);
cout << "\n LoiContraintesPlanesDouble::Passage_deformation_volume_ordre1_vers_3 \n";
Sortie(1);
};
// -- cas de la vitesse de déformation
// recup de l'incrément de temps
double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();
double unSurDeltat=0;
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
{ // non 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;
};
Deps_BB_3D.Coor(2,2) = delta_eps_BB_3D(2,2) * unSurDeltat;
Deps_BB_3D.Coor(3,3) = delta_eps_BB_3D(3,3) * unSurDeltat;
};
// calcul de la déformation d'épaisseur et de largeur
// cas de la méthode 1
void LoiContraintesPlanesDouble::Calcul_eps_trans_parVarVolume1(Tenseur3BH& sigBH_t,double& jaco_1D_0,const double& module_compressibilite,double& jaco_1D_tdt
,Tenseur3BH& sigBH,double& jaco_1D_t)
{ // récup des infos spécifiques
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
// on utilise la relation: (V-V_0)/V = trace(sigma)/3 /K_moy
// avec trace(sigma) en 1D = idem en 3D car on est en contraintes doublement planes
// jacobien 1D = la longueur donc le volume = jacobien 1D * h * b -> (jaco1D_tdt * h * b -jaco1D_0*h_0*b_0) = delta_V
// on considère de plus que l'on est isotrope donc h/h_0 = b/b_0
// en appelant B= trace(sigma)/3 /K_moy on a:
// (h/h_0)^2 = jaco1D_0 / jaco1D_tdt * 1./(1.-B) d'où la déformation logarithmique
if (Dabs(module_compressibilite)>0.) // si le module est nul, cela signifie qu'il n'est pas disponible
{ //double B = (sigBH(1,1)+sigBH(2,2)) / (3. * module_compressibilite );
double log_VsurV0 = (sigBH.Trace()) / (3. * module_compressibilite );
double exp_log_VsurV0 = exp(log_VsurV0);
double hSurh_0 = sqrt((jaco_1D_0 / jaco_1D_tdt) * exp_log_VsurV0) ;
// on test les mini-maxi
if (hSurh_0 < mini_hsurh0)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{if (Permet_affichage() > 3)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
{cout << "\n *** attention hsurh0 " << hSurh_0 << "est inferieur a la limite fixee "<<mini_hsurh0
<< " on le ramene a la limite "
<< "\n LoiContraintesPlanesDouble::Calcul_eps_trans_parVarVolume(... " << endl;
};
hSurh_0 = mini_hsurh0;
};
if (hSurh_0 > maxi_hsurh0)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{if (Permet_affichage() > 3)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
{cout << "\n *** attention bsurb0 " << maxi_hsurh0 << "est superieur a la limite fixee "<<maxi_hsurh0
<< " on le ramene a la limite "
<< "\n LoiContraintesPlanesDouble::Calcul_eps_trans_parVarVolume(... " << endl;
};
hSurh_0 = maxi_hsurh0;
};
// on met à jour épaisseur et largeur d'une même valeur
save_resul.hsurh0 = hSurh_0; // élongation d'épaisseur
save_resul.bsurb0 = hSurh_0; // élongation de largeur
////--- debug
// save_resul.hsurh0 = save_resul.h_tsurh0; // élongation d'épaisseur
// save_resul.bsurb0 = save_resul.b_tsurb0; // élongation d'épaisseur
////-----fin debug
}
else
{ save_resul.hsurh0 = save_resul.bsurb0 = 1.;};
// ---- la suite est un essai, mais ne semble pas non plus fonctionner avec la loi critère
// car le jacobien à t n'est pas calculé ?? ce qui n'est pas normal
// ticket #120
// // a priori cela ne marche pas bien: on va tenter avec le delta jacobien
// // on utilise la relation: (V(tdt)-V(t))/V(tdt) = trace(delta_sigma)/3 /K_moy
// // avec trace(delta_sigma) en 1D = idem en 3D car on est en contraintes doublement planes
// // jacobien 1D = la longueur donc le volume = jacobien 1D * h * b -> (jaco1D_tdt * h * b -jaco1D_t*h_t*b_t) = delta_V
// // on considère de plus que l'on est isotrope donc delta_h/h = delta_b/b d'où h/h_t = b/b_t
// // en appelant B= trace(delta_sigma)/3 /K / (h_t*b_t) on a:
// // (h/h_t)^2 = jaco1D_t / (jaco1D_tdt * (1.-B))
// if (Dabs(module_compressibilite)>0.) // si le module est nul, cela signifie qu'il n'est pas disponible
// { //double B = (sigBH(1,1)+sigBH(2,2)) / (3. * module_compressibilite );
// double B = (sigBH(1,1)-sigBH_t(1,1)) / (3. * module_compressibilite );
// double hSurh_t = sqrt((jaco_1D_t / jaco_1D_tdt) / (1.-B)) ;
//
// save_resul.hsurh0 = hSurh_t * save_resul.h_tsurh0;
// save_resul.bsurb0 = hSurh_t * save_resul.b_tsurb0; // élongation d'épaisseur
// }
// else
// { save_resul.hsurh0 = save_resul.bsurb0 = 1.;};
};
// calcul de la déformation d'épaisseur et de largeur
// cas de la méthode 2
void LoiContraintesPlanesDouble::Calcul_eps_trans_parVarVolume2(Tenseur3BH& sigBH_t,double& jaco_1D_0,const double& module_compressibilite,double& jaco_1D_tdt
,Tenseur3BH& sigBH,double& jaco_1D_t)
{ // récup des infos spécifiques
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
// on utilise la relation: (V-V_0)/V = trace(sigma)/3 /K_moy
// avec trace(sigma) en 1D = idem en 3D car on est en contraintes doublement planes
// jacobien 1D = la longueur donc le volume = jacobien 1D * h * b -> (jaco1D_tdt * h * b -jaco1D_0*h_0*b_0) = delta_V
// on considère de plus que l'on est isotrope donc h/h_0 = b/b_0
// en appelant B= trace(sigma)/3 /K_moy on a:
// (h/h_0)^2 = jaco1D_0 / jaco1D_tdt * 1./(1.-B) d'où la déformation logarithmique
if (Dabs(module_compressibilite)>0.) // si le module est nul, cela signifie qu'il n'est pas disponible
{ //double B = (sigBH(1,1)+sigBH(2,2)) / (3. * module_compressibilite );
double log_VsurV0 = (sigBH.Trace()) / (3. * module_compressibilite );
double exp_log_VsurV0 = exp(log_VsurV0);
double hSurh_0 = sqrt((jaco_1D_0 / jaco_1D_tdt) * exp_log_VsurV0) ;
// on test les mini-maxi
if (hSurh_0 < mini_hsurh0)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{if (Permet_affichage() > 3)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
{cout << "\n *** attention hsurh0 " << hSurh_0 << "est inferieur a la limite fixee "<<mini_hsurh0
<< " on le ramene a la limite "
<< "\n LoiContraintesPlanesDouble::Calcul_eps_trans_parVarVolume2(... " << endl;
};
hSurh_0 = mini_hsurh0;
};
if (hSurh_0 > maxi_hsurh0)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{if (Permet_affichage() > 3)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
{cout << "\n *** attention bsurb0 " << maxi_hsurh0 << "est superieur a la limite fixee "<<maxi_hsurh0
<< " on le ramene a la limite "
<< "\n LoiContraintesPlanesDouble::Calcul_eps_trans_parVarVolume1(... " << endl;
};
hSurh_0 = maxi_hsurh0;
};
// on met à jour épaisseur et largeur d'une même valeur
save_resul.hsurh0 = hSurh_0; // élongation d'épaisseur
save_resul.bsurb0 = hSurh_0; // élongation de largeur
////--- debug
// save_resul.hsurh0 = save_resul.h_tsurh0; // élongation d'épaisseur
// save_resul.bsurb0 = save_resul.b_tsurb0; // élongation d'épaisseur
////-----fin debug
}
else
{ save_resul.hsurh0 = save_resul.bsurb0 = 1.;};
// ---- la suite est un essai, mais ne semble pas non plus fonctionner avec la loi critère
// car le jacobien à t n'est pas calculé ?? ce qui n'est pas normal
// ticket #120
// // a priori cela ne marche pas bien: on va tenter avec le delta jacobien
// // on utilise la relation: (V(tdt)-V(t))/V(tdt) = trace(delta_sigma)/3 /K_moy
// // avec trace(delta_sigma) en 1D = idem en 3D car on est en contraintes doublement planes
// // jacobien 1D = la longueur donc le volume = jacobien 1D * h * b -> (jaco1D_tdt * h * b -jaco1D_t*h_t*b_t) = delta_V
// // on considère de plus que l'on est isotrope donc delta_h/h = delta_b/b d'où h/h_t = b/b_t
// // en appelant B= trace(delta_sigma)/3 /K / (h_t*b_t) on a:
// // (h/h_t)^2 = jaco1D_t / (jaco1D_tdt * (1.-B))
// if (Dabs(module_compressibilite)>0.) // si le module est nul, cela signifie qu'il n'est pas disponible
// { //double B = (sigBH(1,1)+sigBH(2,2)) / (3. * module_compressibilite );
// double B = (sigBH(1,1)-sigBH_t(1,1)) / (3. * module_compressibilite );
// double hSurh_t = sqrt((jaco_1D_t / jaco_1D_tdt) / (1.-B)) ;
//
// save_resul.hsurh0 = hSurh_t * save_resul.h_tsurh0;
// save_resul.bsurb0 = hSurh_t * save_resul.b_tsurb0; // élongation d'épaisseur
// }
// else
// { save_resul.hsurh0 = save_resul.bsurb0 = 1.;};
};
// calcul des déformations d'épaisseur et de largeur et des variations
// méthode 1
void LoiContraintesPlanesDouble::Calcul_d_eps_trans_parVarVolume(double& jaco_1D_0,const double& module_compressibilite,double& jaco_1D_tdt
,TenseurHH& sigHH_,Vecteur& d_jaco_1D
,Tableau <TenseurHH *>& d_sigHH,Tableau <TenseurBB *>& d_gijBB_tdt,TenseurBB & gijBB_)
{ // récup des infos spécifiques
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
// on utilise la relation: log(V/V0) = trace(sigma)/3 /K_moy
// avec trace(sigma) en 1D = idem en 3D car on est en contraintes planes double
// jacobien 1D = la longueur, donc le volume: V = jacobien 1D * h * b -> (jaco1D_tdt * h * b)/(jaco1D_0*h_0*b_0) = V/V0
// en appelant log_VsurV0= trace(sigma)/3 /K_moy on a:
// h/h_0 * b/b_0= jaco1D_0 / jaco1D_tdt * exp(log_VsurV0)
// ici sans info sup, on considère une loi isotrope d'où h/h_0 = b/b_0
// -> h/h_0 = (jaco1D_0 / jaco1D_tdt * exp(log_VsurV0))^(0.5)
if (Dabs(module_compressibilite)>0.) // si le module est nul, cela signifie qu'il n'est pas disponible
{ const Tenseur1BB & gijBB = *((Tenseur1BB*) &gijBB_); // pour simplifier
const Tenseur1HH & sigHH = *((Tenseur1HH*) &sigHH_); // pour simplifier
Tenseur1BH sigBH = gijBB * sigHH;
double log_VsurV0 = (sigBH.Trace()) / (3. * module_compressibilite );
double exp_log_VsurV0 = exp(log_VsurV0);
double hSurh_0 = sqrt((jaco_1D_0 / jaco_1D_tdt) * exp_log_VsurV0) ;
// on test les mini-maxi
if (hSurh_0 < mini_hsurh0)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{if (Permet_affichage() > 3)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
{cout << "\n *** attention hsurh0 " << hSurh_0 << "est inferieur a la limite fixee "<<mini_hsurh0
<< " on le ramene a la limite "
<< "\n LoiContraintesPlanesDouble::Calcul_d_eps_trans_parVarVolume(... " << endl;
};
hSurh_0 = mini_hsurh0;
};
if (hSurh_0 > maxi_hsurh0)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{if (Permet_affichage() > 3)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
{cout << "\n *** attention bsurb0 " << maxi_hsurh0 << "est superieur a la limite fixee "<<maxi_hsurh0
<< " on le ramene a la limite "
<< "\n LoiContraintesPlanesDouble::Calcul_d_eps_trans_parVarVolume(... " << endl;
};
hSurh_0 = maxi_hsurh0;
};
// on met à jour épaisseur et largeur d'une même valeur
save_resul.hsurh0 = hSurh_0; // élongation d'épaisseur
save_resul.bsurb0 = hSurh_0; // élongation de largeur
// maintenant le cas de la variation de h/h0
// on tient compte de la variation relative à la dépendance du jacobien et à la trace de sigma
// a) vérification et modif éventuelle de la taille
Vecteur& d_hsurh0 = save_resul.d_hsurh0;
// b) calcul de la variation
//2.* h/h_0 * d h/h0 = - djaco_1D * (jaco_1D_0 / ((jaco_1D_tdt)^2)) * exp(log_VsurV0)
// + d_(sigBH.Trace()) * (jaco_1D_0 / jaco_1D_tdt) * exp(log_VsurV0) / (3. * module_compressibilite )
int taille = d_hsurh0.Taille();
double fact1 = - jaco_1D_0 / (jaco_1D_tdt*jaco_1D_tdt) * exp_log_VsurV0 / (2. * hSurh_0);
double fact2 = (jaco_1D_0 / jaco_1D_tdt) * exp_log_VsurV0 / (3. * module_compressibilite / (2. * hSurh_0));
////---debug
//cout << "\n Calcul_d_eps33_parVarVolume() module_compressibilite= " << module_compressibilite
// << " log_VsurV0= "<< log_VsurV0
// << " sigBH.Trace()= " << sigBH.Trace()
// << " hSurh_0 =" << hSurh_0
// << "\n fact1 =" << fact1 << " fact2 =" << fact2;
//// --- fin debug
for (int i=1;i<=taille;i++)
{
Tenseur1HH & dsigHH = *((Tenseur1HH*) (d_sigHH(i))); // passage en dim 1
const Tenseur1BB & dgijBB = *((Tenseur1BB*)(d_gijBB_tdt(i))) ; // pour simplifier l'ecriture
Tenseur1BH d_sigBH = gijBB * dsigHH + dgijBB * sigHH;
d_hsurh0(i) = fact1 * d_jaco_1D(i) + fact2 * (d_sigBH.Trace());
//---debug
//cout << "\n Calcul_d_eps33_parVarVolume() "
// << " d_hsurh0("<< i<<")= "<< d_hsurh0(i);
// --- fin debug
};
// en isotrope les axes h et b sont équivalents
save_resul.d_bsurb0 = save_resul.d_hsurh0;
}
else
{ save_resul.hsurh0 = save_resul.bsurb0 = 1.;
save_resul.d_hsurh0.Zero();save_resul.d_bsurb0.Zero();
};
};
// calcul de la variation des déformations d'hauteur et de largeur, par rapport à la def_11, à cause
// des contraintes planes doubles: première méthode
void LoiContraintesPlanesDouble::Calcul_d_eps_eg_11()
{ // on prépare la résolution
// ,d_sig_ef_gh(2,2),d_eps_ef_11(2)
// la matrice
double a = d_sig_ef_gh(1,1) = d_sigma_deps_3D(2,2,2,2);
double b = d_sig_ef_gh(1,2) = d_sigma_deps_3D(2,2,3,3);
double c = d_sig_ef_gh(2,1) = d_sigma_deps_3D(3,3,2,2);
double d = d_sig_ef_gh(2,2) = d_sigma_deps_3D(3,3,3,3);
double determinant = a*d - b*c;
// //----- debug
// cout << "\n debug LoiContraintesPlanesDouble::Calcul_d_eps_eg_11() "<< "\n d_sig_ef_gh= ";
// d_sig_ef_gh.Affiche();
// //---fin debug
if (Dabs(determinant) > ConstMath::petit)
{ // le second membre
d_eps_ef_11(1)=-d_sigma_deps_3D(2,2,1,1);
d_eps_ef_11(2)=-d_sigma_deps_3D(3,3,1,1);
// //----- debug
// cout << "\n second membre= "<< d_eps_ef_11(1) << ", " << d_eps_ef_11(2);
// //---fin debug
// résolution
d_sig_ef_gh.Resol_systID(d_eps_ef_11);
// //----- debug
// cout << "\n solution= "<< d_eps_ef_11(1) << ", " << d_eps_ef_11(2);
// // vérif
// a = d_sig_ef_gh(1,1) = d_sigma_deps_3D(2,2,2,2);
// b = d_sig_ef_gh(1,2) = d_sigma_deps_3D(2,2,3,3);
// c = d_sig_ef_gh(2,1) = d_sigma_deps_3D(3,3,2,2);
// d = d_sig_ef_gh(2,2) = d_sigma_deps_3D(3,3,3,3);
// Vecteur truc(2);
// truc = d_sig_ef_gh * d_eps_ef_11;
// cout << "\n second membre correspondant "<< truc(1) << ", " << truc(2)<< flush;
// //---fin debug
}
else // sinon on met à 0
{d_eps_ef_11.Inita(0.);};
};
// calcul des invariants de déformation et de vitesse de déformation, et les def cumulées
// correspondant aux cas 3D: valable pour les 2 méthodes
void LoiContraintesPlanesDouble::Calcul_invariants_et_def_cumul()
{ // récup des infos spécifiques
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
// recup de l'incrément de temps
double deltat=ParaGlob::Variables_de_temps().IncreTempsCourant();
double unSurDeltat=0;
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
{ // non 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;
};
// ---calcul des invariants de déformation et de vitesse de déformation
Tenseur3BH espBH = eps_BB_3D * gijHH_tdt_3D;
Tenseur3BH delta_espBH = delta_eps_BB_3D * gijHH_tdt_3D;
save_resul.epsInvar(1) = espBH.Trace();
save_resul.epsInvar(2) = espBH.II();
save_resul.epsInvar(3) = espBH.Det();
save_resul.depsInvar(1) = unSurDeltat * delta_espBH.Trace();
save_resul.depsInvar(2) = unSurDeltat * delta_espBH.II();
save_resul.depsInvar(3) = unSurDeltat * delta_espBH.Det();
// cas des grandeurs cumulées (cf. la classe Deformation)
// tableau relatif aux différentes grandeurs de type def scalaires équivalentes
// def_equi(1) = deformation cumulée = somme des sqrt(2./3. * (delta_eps_barre_BH && delta_eps_barre_BH)) ;
// def_equi(2) = deformation duale de la contrainte de mises = sqrt(2./3. * (eps_barre_BH && eps_barre_BH)) ;
// def_equi(3) = niveau maxi atteind par def_equi(2)
// def_equi(4) = delta def cumulée = sqrt(2./3. * (delta_eps_barre_BH && delta_eps_barre_BH));
double delta_eps_equi = sqrt(2./3. * ( (delta_espBH && delta_espBH) - Sqr(delta_espBH.Trace()) /3. ));
save_resul.def_equi(1) = save_resul.def_equi_t(1) + delta_eps_equi ;
save_resul.def_equi(2) = sqrt(2./3. * (espBH && espBH)) ;
if (save_resul.def_equi(2) > save_resul.def_equi_t(3))
save_resul.def_equi(3) = save_resul.def_equi(2);
save_resul.def_equi(4) = delta_eps_equi;
};
// calcul de la fonction résidu de la résolution de l'équation constitutive: sig^ef(h/h0,b/b0) = 0
// avec "e,f" = 2 et 3
// h/h0, bb0 sont les inconnues du problème et sont l'élément d'entrée: <x(1),x(2)>
// l'argument test ramène
// . 1 si le calcul a été ok, -1 s'il y a eu un pb, mais on peut continuer, 0 s'il y a eu un pb
// fatal, qui invalide le calcul du résidu
Vecteur& LoiContraintesPlanesDouble::Residu_constitutif (const double & alpha,const Vecteur & x, int& test)
{
// récup du conteneur spécifique
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
test = 1; // init par défaut
// %%% on met à jour les conteneurs locaux en fonction des paramètres d'entrée %%%
// -- prise en compte de la variation de la déformation d'épaisseur
// pour nous ici on travaille avec giB_03D qui est normé et = 1 c-a-d h0 = 1
// donc à tdt il faut multiplier par h/h0 à tdt
double bsurb0 = x(1);double hsurh0 = x(2);
// vérif des grandeurs
bool erreur_sortie = false; // pour gestion de sortie d'infos supplémentaires
if (Limitation_h_b(hsurh0,bsurb0))
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{if (Permet_affichage() > 3)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
cout << "\n LoiContraintesPlanesDouble::Residu_constitutif(... " << endl;
erreur_sortie = true;
test = -1; // l'erreur n'est pas fatale mais elle signifie que ce n'est pas très normal !!
};
if (erreur_sortie && (Permet_affichage() > 4))
{cout << "\n giH_tdt: " << *(umat_cont_3D->giH_tdt);
cout << "\n defBB: ";eps_BB_3D.Ecriture(cout);
cout << "\n DepsBB: ";Deps_BB_3D.Ecriture(cout);
};
save_resul.hsurh0 = hsurh0; // sauvegarde de l'élongation d'épaisseur
save_resul.bsurb0 = bsurb0; // sauvegarde de l'élongation de largeur
CoordonneeB & giBtdt_2 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(2); // pour simplifier
CoordonneeB & giBtdt_3 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(3); // pour simplifier
giB_tdt_3D.CoordoB(3) = giBtdt_3 * hsurh0; giH_tdt_3D.CoordoH(3) = (giBtdt_3 / hsurh0).Bas_haut();
giB_tdt_3D.CoordoB(2) = giBtdt_2 * bsurb0; giH_tdt_3D.CoordoH(2) = (giBtdt_2 / bsurb0).Bas_haut();
// - et les tenseurs métriques à t et tdt, comme les vecteurs en épaisseurs
// et en largeurs sont normées initialement, seule l'intensité
// des normales est à mettre à jour compte tenue de la variation d'épaisseurs et de largeur
double hsurh0_2 = hsurh0 * hsurh0;
gijBB_tdt_3D.Coor(3,3) = hsurh0_2; gijHH_tdt_3D.Coor(3,3) = 1. / hsurh0_2;
double bsurb0_2 = bsurb0 * bsurb0;
gijBB_tdt_3D.Coor(2,2) = bsurb0_2; gijHH_tdt_3D.Coor(2,2) = 1. / bsurb0_2;
// mise à jour des informations liées à la déformation de 1 vers 3
Mise_a_jour_deformations_et_Jacobien_en_3D();
// calcul des invariants de déformation et de vitesse de déformation, et les def cumulées
// correspondant aux cas 3D
Calcul_invariants_et_def_cumul();
// initialisation du comportement tangent / au def
d_sigma_deps_3D.Inita(0.);
// appel du calcul de sig et dsig
bool en_base_orthonormee = false; // ici les tenseurs ne sont pas forcément en orthonormee
lois_interne->Calcul_dsigma_deps (en_base_orthonormee, sig_HH_t_3D,Deps_BB_3D,eps_BB_3D
,delta_eps_BB_3D,jacobien_0_3D,jacobien_tdt_3D
,sig_HH_3D,d_sigma_deps_3D
,save_resul.l_energ,save_resul.l_energ_t,module_compressibilite_3D
,module_cisaillement_3D,*umat_cont_3D);
// retour des grandeurs voulues
residu(1) = sig_HH_3D(2,2);
residu(2) = sig_HH_3D(3,3);
return residu;
};
// calcul de la matrice tangente de la résolution de l'équation constitutive: sig^ef(h/h0,b/b0) = 0
// avec "e,f" = 2 et 3
// h/h0, bb0 sont les inconnues du problème et sont l'élément d'entrée: <x(1),x(2)>
// l'argument test ramène
// . 1 si le calcul a été ok, -1 s'il y a eu un pb, mais on peut continuer, 0 s'il y a eu un pb
// fatal, qui invalide le calcul du résidu et de la dérivée
Mat_abstraite& LoiContraintesPlanesDouble::Mat_tangente_constitutif
(const double & alphap,const Vecteur & x, Vecteur& residu, int& test)
{
// récup du conteneur spécifique
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
test = 1; // init par défaut
// %%% on met à jour les conteneurs locaux en fonction des paramètres d'entrée %%%
// -- prise en compte de la variation de la déformation d'épaisseur et de largeur
// pour nous ici on travail avec giB_03D qui est normé et = 1 c-a-d h0 = b0 = 1
// donc à tdt il faut multiplier par h/h0 à tdt et idem pour b/b0
double bsurb0 = x(1);double hsurh0 = x(2);
// vérif des grandeurs
bool erreur_sortie = false; // pour gestion de sortie d'infos supplémentaires
if (Limitation_h_b(hsurh0,bsurb0))
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{if (Permet_affichage() > 3)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
cout << "\n LoiContraintesPlanesDouble::Mat_tangente_constitutif(... " << endl;
erreur_sortie = true;
test = -1; // l'erreur n'est pas fatale mais elle signifie que ce n'est pas très normal !!
};
if (erreur_sortie && (Permet_affichage() > 4))
{cout << "\n giH_tdt: " << *(umat_cont_3D->giH_tdt);
cout << "\n defBB: ";eps_BB_3D.Ecriture(cout);
cout << "\n DepsBB: ";Deps_BB_3D.Ecriture(cout);
};
save_resul.hsurh0 = hsurh0; // sauvegarde de l'élongation d'épaisseur
save_resul.bsurb0 = bsurb0; // sauvegarde de l'élongation de largeur
CoordonneeB & giBtdt_2 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(2); // pour simplifier
CoordonneeB & giBtdt_3 = (*(save_resul.meti_ref_tdt.giB_)).CoordoB(3); // pour simplifier
giB_tdt_3D.CoordoB(3) = giBtdt_3 * hsurh0;
giH_tdt_3D.CoordoH(3) = (giBtdt_3 / hsurh0).Bas_haut();
giB_tdt_3D.CoordoB(2) = giBtdt_2 * bsurb0;
giH_tdt_3D.CoordoH(2) = (giBtdt_2 / bsurb0).Bas_haut();
// - et les tenseurs métriques à t et tdt, comme les vecteurs en épaisseurs sont normées initialement, seule l'intensité
// des normales est à mettre à jour compte tenue de la variation d'épaisseurs
double hsurh0_2 = hsurh0 * hsurh0;
gijBB_tdt_3D.Coor(3,3) = hsurh0_2;
gijHH_tdt_3D.Coor(3,3) = 1. / hsurh0_2;
double bsurb0_2 = bsurb0 * bsurb0;
gijBB_tdt_3D.Coor(2,2) = bsurb0_2;
gijHH_tdt_3D.Coor(2,2) = 1. / bsurb0_2;
// mise à jour des informations liées à la déformation de 1 vers 3
Mise_a_jour_deformations_et_Jacobien_en_3D();
// calcul des invariants de déformation et de vitesse de déformation, et les def cumulées
// correspondant aux cas 3D
Calcul_invariants_et_def_cumul();
// initialisation du comportement tangent / au def
d_sigma_deps_3D.Inita(0.);
// appel du calcul de sig et dsig
bool en_base_orthonormee = false; // ici les tenseurs ne sont pas a priori en orthonormee
lois_interne->Calcul_dsigma_deps (en_base_orthonormee, sig_HH_t_3D,Deps_BB_3D,eps_BB_3D
,delta_eps_BB_3D,jacobien_0_3D,jacobien_tdt_3D
,sig_HH_3D,d_sigma_deps_3D
,save_resul.l_energ,save_resul.l_energ_t,module_compressibilite_3D
,module_cisaillement_3D,*umat_cont_3D);
// retour des grandeurs voulues
residu(1) = sig_HH_3D(2,2);
residu(2) = sig_HH_3D(3,3);
// la matrice tangente dépend du type de mesure de déformation, avec la relation
// en appelant hb/hb0 les deux ddl: h/h0 et b/b0, i,j,k et l variant de 2 à 3
// d_sig^ij(h/h0,b/b0) / d_ddl = d_sig^ij(h/h0,b/b0)/d_eps_kl * d_eps_kl/d_ddl
// comme les repères restent orthogonales à tous moments il n'y a pas de couplage entre les directions 2 et 3
Vecteur d_eps_kl_sur_d_ddl(2,0.); // init
// -- choix en fonction du type de déformation
switch (type_de_deformation)
{case DEFORMATION_STANDART : case DEFORMATION_POUTRE_PLAQUE_STANDART :
// cas d'une déformation d'Almansi
{ // epsBB33 = 1/2 * (1. - (h0/h)^2), en orthonormee
// dans le repère local: epsBB33 = 1/2 * (h^2 - 1.), or h0=1. donc : epsBB33 = 1/2 * ((h/h0)^2 - 1.)
// idem pour la direction 22
d_eps_kl_sur_d_ddl(1) = bsurb0;d_eps_kl_sur_d_ddl(2) = hsurh0;
};
break;
case DEFORMATION_LOGARITHMIQUE : case DEF_CUMUL_CORROTATIONNEL : case DEFORMATION_CUMU_LOGARITHMIQUE :
// cas d'une def logarithmique ou une approximation: eps_BB_3D(3,3) = log(hsurh0);
{ d_eps_kl_sur_d_ddl(1) = log(bsurb0);d_eps_kl_sur_d_ddl(2) = log(hsurh0);
};
break;
default :
cout << "\nErreur : type de deformation qui n'est pas actuellement pris en compte, type= "
<< Nom_type_deformation(type_de_deformation);
cout << "\n LoiContraintesPlanesDouble::Mat_tangente_constitutif \n";
Sortie(1);
};
// d'où la dérivée du résidu, en considérant qu'il n'y a pas de couplage entre les directions 2 et 3
// c-a-d d_eps33_sur_d_bSurb0 = 0 et d_eps22_sur_d_hSurh0 = 0
// non c'est faux, du à la partie sphérique des contraintes par exemple, il faut tout considérer
derResidu(1,1) = d_sigma_deps_3D(2,2,2,2) * d_eps_kl_sur_d_ddl(1);
derResidu(1,2) = d_sigma_deps_3D(2,2,3,3) * d_eps_kl_sur_d_ddl(2);
derResidu(2,1) = d_sigma_deps_3D(3,3,2,2) * d_eps_kl_sur_d_ddl(1);
derResidu(2,2) = d_sigma_deps_3D(3,3,3,3) * d_eps_kl_sur_d_ddl(2);
return derResidu;
};
//--- cas de la résolution de l'équation sigij(eps_meca)*V_j{.l}=0, l=2 et 3
// calcul de la fonction résidu de la résolution de l'équation constitutive
// l'argument test ramène
// . 1 si le calcul a été ok, -1 s'il y a eu un pb, mais on peut continuer, 0 s'il y a eu un pb
// fatal, qui invalide le calcul du résidu
Vecteur& LoiContraintesPlanesDouble::Residu_constitutif_3D(const double & alpha,const Vecteur & x, int& test)
{
// récup du conteneur spécifique
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
test = 1; // init par défaut
//récup de eps_mecaBB_t
Tenseur3BB& eps_mecaBB_t = *((Tenseur3BB*) (save_resul.eps_mecaBB_t));
// les inconnues sont les déformations mécaniques {epsilon'}_kl;
// (k,l) = (22) (33) (12) : on a donc 3 inconnues que l'on suppose rangées dans cette ordre
// c-a-d les def qui vont servir dans l'appel de la loi 3D
eps_P_BB.Coor(1,2) = x(3);eps_P_BB.Coor(2,2) = x(1);eps_P_BB.Coor(3,3) = x(2);
eps_P_BB.Coor(1,3) = 0.; eps_P_BB.Coor(2,3) = 0.;
if (Permet_affichage() > 4)
{cout << "\n eps_P_BB (dans Va) avant limitation : " ;eps_P_BB.Ecriture(cout);}
// vérif des grandeurs et calcul des variations d'épaisseur et de largeur
double& hsurh0 = save_resul.hsurh0;double& bsurb0 = save_resul.bsurb0;
// la méthode va mettre à jour les variations bsurb0 et hsurh0 et le jacobien
bool erreur_sortie = false; // pour gestion de sortie d'infos supplémentaires
int info_limit = Calcul_et_Limitation2_h_b(hsurh0,eps_P_BB,bsurb0,jacobien_tdt_3D,jacobien_0_3D);
if (info_limit)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{if (Permet_affichage() > 3)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
cout << "\n LoiContraintesPlanesDouble::Residu_constitutif_3D(... " << endl;
erreur_sortie = true;test = info_limit - 2;
};
if (erreur_sortie && (Permet_affichage() > 4))
{cout << "\n giH_tdt: " << *(umat_cont_3D->giH_tdt);
cout << "\n eps_BB_3D: ";eps_BB_3D.Ecriture(cout);
cout << "\n DepsBB: ";Deps_BB_3D.Ecriture(cout);
cout << "\n x(1)= "<< x(1) << " x(2)= "<< x(2) << " x(3)= "<< x(3) << flush;
};
// dans la direction V_1 c'est la déformation issue de la cinématique
// elle a déjà été calculée, et les autres composantes demeurent nulle
// donc eps_P_BB est au complet,
eps_BB_3D=eps_P_BB; // recopie dans eps_BB_3D: là on est dans le repère V_a
eps_BB_3D.ChBase(beta3D); // changement de base, là on est dans le repère g^i
// on s'occupe maintenant de l'incrément de déformation
delta_eps_BB_3D = eps_BB_3D - eps_mecaBB_t;
// au cas où, on impose (normalement c'est inutile) que les composantes 13 et 23 sont nulles
delta_eps_BB_3D.Coor(1,3) = 0.; delta_eps_BB_3D.Coor(2,3) = 0.;
// 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
#ifdef MISE_AU_POINT
if (unSurDeltat < 0)
{ cout << "\n le pas de temps est négatif !! "
<< "\n LoiContraintesPlanesDouble::Residu_constitutif_3D(...";
};
#endif
unSurDeltat = ConstMath::tresgrand;
};
// ici on considère que la déformation est de type Almansi
Deps_BB_3D = delta_eps_BB_3D * unSurDeltat;
// calcul des invariants de déformation et de vitesse de déformation, et les def cumulées
// correspondant aux cas 3D
Calcul_invariants_et_def_cumul();
/* // --- on s'occupe des jacobiens
// à tdt le jacobien est calculée à partir de la déformation supposée d'Almansi
// ici il s'agit du jacobien relatif au comportement mécanique
// |J| = sqrt(det(2 * eps_ij + g_ij(0)))
// calcul de la matrice correspondante à (2 * eps_ij + g_ij(0))
gij_meca_BB = (2. * eps_BB_3D + (*(umat_cont_3D->gijBB_0)));
gij_meca_BB.Matrice_composante(mat_inter);
// calcul du jacobien
jacobien_tdt_3D = sqrt(Dabs(mat_inter.Determinant()));
*/
// initialisation du comportement tangent / au def
d_sigma_deps_3D.Inita(0.);
// appel du calcul de sig et dsig
bool en_base_orthonormee = false; // ici les tenseurs ne sont pas en orthonormee
// car on travaille dans la base naturelle
#ifdef MISE_AU_POINT
if (Permet_affichage() > 8)
{cout << "\n LoiContraintesPlanesDouble::Residu_constitutif_3D(..."
<< "\n giH_tdt: " << *(umat_cont_3D->giH_tdt);
cout << "\n defBB: ";eps_BB_3D.Ecriture(cout);
cout << "\n DepsBB: ";Deps_BB_3D.Ecriture(cout);
cout << "\n x(1)= "<< x(1) << " x(2)= "<< x(2) << " x(3)= "<< x(3) ;
// en base orthonormee
Tenseur3BB epsBB_3D_inter;
eps_BB_3D.BaseAbsolue(epsBB_3D_inter,*(umat_cont_3D->giH_tdt));
cout << "\n eps_BB_3D en absolue 3D (dans Va): "<<epsBB_3D_inter;
cout << "\n jacobien_tdt_3D= "<< jacobien_tdt_3D
<< "\n jacobien_0= " << jacobien_0_3D
<< flush;
};
#endif
lois_interne->Calcul_dsigma_deps (en_base_orthonormee, sig_HH_t_3D,Deps_BB_3D,eps_BB_3D
,delta_eps_BB_3D,jacobien_0_3D,jacobien_tdt_3D
,sig_HH_3D,d_sigma_deps_3D
,save_resul.l_energ,save_resul.l_energ_t,module_compressibilite_3D
,module_cisaillement_3D,*umat_cont_3D);
#ifdef MISE_AU_POINT
if (Permet_affichage() > 8)
{cout << "\n sig_HH_3D: " << sig_HH_3D;
cout << "\n d_sigma_deps_3D: "<< d_sigma_deps_3D << flush;
};
#endif
//-------- calcul du residu -----------
// residu(mi) = V_{mj} * \sigma^{ij}(\epsilon'_{kl}) : m=2 et 3
CoordonneeH V_sig_H(sig_HH_3D * (*ViB_3D)(2));
// on considère que les densités d'effort interne dans le plan, la composante 23 devant
// être nulle compte tenu du tenseur des déformations
// residu_3D(1) = V_sig_H(1);residu_3D(3) = V_sig_H(2);
residu_3D(1) = V_sig_H * (*ViB_3D)(2); // composante 22
residu_3D(3) = (*ViB_3D)(1) * V_sig_H ; // composante 12
// essai avec le mixte
// residu_3D(3) = 0.5 * ((*ViB_3D)(1) * V_sig_H + V_sig_H * (*ViB_3D)(1)); // composante 12
#ifdef MISE_AU_POINT
if (Permet_affichage() > 8)
{cout << "\n en g_i: sig_HH * ViB(2)= " << V_sig_H
<< " df_ds_22= " << residu_3D(1) << " df_ds_21= " << residu_3D(3)
<< flush;
};
#endif
V_sig_H = sig_HH_3D * (*ViB_3D)(3);
// residu_3D(2) = V_sig_H(3); // seule la composante dans la direction 3 est utilisée
residu_3D(2) = V_sig_H * (*ViB_3D)(3); // composante 33
#ifdef MISE_AU_POINT
if (Permet_affichage() > 8)
{cout << "\n en g_i: sig_HH * ViB(3)= " << V_sig_H << " df_ds_33= " << residu_3D(2);
// pour info on calcul la projection symétrique
double f_12 = (*ViB_3D)(1)* sig_HH_3D * (*ViB_3D)(2);
cout << " df_ds_12= " << f_12 << flush;
};
#endif
return residu_3D;
};
//--- cas de la résolution de l'équation sigij(eps_meca)*V_j{.l}=0, l=2 et 3
// calcul de la matrice tangente de la résolution de l'équation constitutive
// l'argument test ramène
// . 1 si le calcul a été ok, -1 s'il y a eu un pb, mais on peut continuer, 0 s'il y a eu un pb
// fatal, qui invalide le calcul du résidu et de la dérivée
Mat_abstraite& LoiContraintesPlanesDouble::Mat_tangente_constitutif_3D
(const double & alphap,const Vecteur & x, Vecteur& residu_local, int& test)
{
// récup du conteneur spécifique
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
test = 1; // init par défaut
//récup de eps_mecaBB_t
Tenseur3BB& eps_mecaBB_t = *((Tenseur3BB*) (save_resul.eps_mecaBB_t));
// les inconnues sont les déformations mécaniques {epsilon'}_kl;
// (k,l) = (22) (33) (12) : on a donc 3 inconnues que l'on suppose rangées dans cette ordre
// c-a-d les def qui vont servir dans l'appel de la loi 3D
eps_P_BB.Coor(1,2) = x(3);eps_P_BB.Coor(2,2) = x(1);eps_P_BB.Coor(3,3) = x(2);
eps_P_BB.Coor(1,3) = 0.; eps_P_BB.Coor(2,3) = 0.;
// vérif des grandeurs et calcul des variations d'épaisseur et de largeur
double& hsurh0 = save_resul.hsurh0;double& bsurb0 = save_resul.bsurb0;
// la méthode va mettre à jour les valeurs de bsurb0 et hsurh0, et le jacobien
bool erreur_sortie = false; // pour gestion de sortie d'infos supplémentaires
if (Permet_affichage() > 4)
{cout << "\n eps_P_BB (dans Va) avant limitation : " ;eps_P_BB.Ecriture(cout);};
int info_limit = Calcul_et_Limitation2_h_b(hsurh0,eps_P_BB,bsurb0,jacobien_tdt_3D,jacobien_0_3D);
if (info_limit)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{if (Permet_affichage() > 3)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
cout << "\n LoiContraintesPlanesDouble::Mat_tangente_constitutif_3D(... " << endl;
erreur_sortie = true;test = info_limit - 2;
};
if (erreur_sortie && (Permet_affichage() > 4))
{cout << "\n giH_tdt: " << *(umat_cont_3D->giH_tdt);
cout << "\n defBB: ";eps_BB_3D.Ecriture(cout);
cout << "\n DepsBB: ";Deps_BB_3D.Ecriture(cout);
cout << "\n x(1)= "<< x(1) << " x(2)= "<< x(2) << " x(3)= "<< x(3) << flush;
};
// dans la direction V_1 c'est la déformation issue de la cinématique
// elle a déjà été calculée, et les autres composantes demeurent nulle
// donc eps_P_BB est au complet,
eps_BB_3D=eps_P_BB; // recopie dans eps_BB_3D: là on est dans le repère V_a
eps_BB_3D.ChBase(beta3D); // changement de base, là on est dans le repère g^i
// on s'occupe maintenant de l'incrément de déformation
delta_eps_BB_3D = eps_BB_3D - eps_mecaBB_t;
// au cas où, on impose (normalement c'est inutile) que les composantes 13 et 23 sont nulles
delta_eps_BB_3D.Coor(1,3) = 0.; delta_eps_BB_3D.Coor(2,3) = 0.;
// 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
#ifdef MISE_AU_POINT
if (unSurDeltat < 0)
{ cout << "\n le pas de temps est négatif !! "
<< "\n LoiContraintesPlanesDouble::Mat_tangente_constitutif_3D(...";
};
#endif
unSurDeltat = ConstMath::tresgrand;
};
// ici on considère que la déformation est de type Almansi
Deps_BB_3D = delta_eps_BB_3D * unSurDeltat;
// calcul des invariants de déformation et de vitesse de déformation, et les def cumulées
// correspondant aux cas 3D
Calcul_invariants_et_def_cumul();
/* // --- on s'occupe des jacobiens
// à t=0 on considère que c'est le même jacobien que transmis
jacobien_0_3D = *(umat_cont_3D->jacobien_0);
// à tdt le jacobien est calculée à partir de la déformation supposée d'Almansi
// ici il s'agit du jacobien relatif au comportement mécanique
// |J| = sqrt(det(2 * eps_ij + g_ij(0)))
// calcul de la matrice correspondante à (2 * eps_ij + g_ij(0))
gij_meca_BB = (2. * eps_BB_3D + (*(umat_cont_3D->gijBB_0)));
gij_meca_BB.Matrice_composante(mat_inter);
// calcul du jacobien
jacobien_tdt_3D = sqrt(Dabs(mat_inter.Determinant()));
*/
// initialisation du comportement tangent / au def
d_sigma_deps_3D.Inita(0.);
// appel du calcul de sig et dsig
bool en_base_orthonormee = false; // ici les tenseurs ne sont pas a priori en orthonormee
#ifdef MISE_AU_POINT
if (Permet_affichage() > 8)
{cout << "\n LoiContraintesPlanesDouble::Mat_tangente_constitutif_3D(..."
<< "\n giH_tdt: " << *(umat_cont_3D->giH_tdt);
cout << "\n defBB: ";eps_BB_3D.Ecriture(cout);
cout << "\n DepsBB: ";Deps_BB_3D.Ecriture(cout);
cout << "\n x(1)= "<< x(1) << " x(2)= "<< x(2) << " x(3)= "<< x(3) ;
// en base orthonormee
Tenseur3BB epsBB_3D_inter;
eps_BB_3D.BaseAbsolue(epsBB_3D_inter,*(umat_cont_3D->giH_tdt));
cout << "\n eps_BB_3D en absolue 3D (dans Va): "<<epsBB_3D_inter;
cout << "\n jacobien_tdt_3D= "<< jacobien_tdt_3D
<< "\n jacobien_0= " << jacobien_0_3D
<< flush;
};
#endif
lois_interne->Calcul_dsigma_deps (en_base_orthonormee, sig_HH_t_3D,Deps_BB_3D,eps_BB_3D
,delta_eps_BB_3D,jacobien_0_3D,jacobien_tdt_3D
,sig_HH_3D,d_sigma_deps_3D
,save_resul.l_energ,save_resul.l_energ_t,module_compressibilite_3D
,module_cisaillement_3D,*umat_cont_3D);
#ifdef MISE_AU_POINT
if (Permet_affichage() > 8)
{cout << "\n sig_HH_3D: " << sig_HH_3D;
cout << "\n d_sigma_deps_3D: "<< d_sigma_deps_3D << flush;
};
#endif
//-------- calcul du residu -----------
// residu(mi) = V_{mj} * \sigma^{ij}(\epsilon'_{kl}) : m=2 et 3
//residu(1) = V_2 . \sigma . V_2 = V_{2i} ~\sigma^{ij} ~V_{2j} = ( V_2 \otimes V_2) : \sigma
//residu(2) = V_3 . \sigma . V_3 = V_{3i} ~\sigma^{ij} ~V_{3j} = ( V_3 \otimes V_3) : \sigma
//residu(3) = V_1 . \sigma . V_2 = V_{1i} ~\sigma^{ij} ~V_{2j} = ( V_1 \otimes V_2) : \sigma
CoordonneeH V_sig_H(sig_HH_3D * (*ViB_3D)(2)); // = \sigma . V_2
// ici on utilise le vecteur passé en argument
residu_local(1) = V_sig_H * (*ViB_3D)(2); // composante 22 = V_2 . \sigma . V_2
// residu(3) = V_sig_H * (*ViB_3D)(1); // composante 21
residu_local(3) = (*ViB_3D)(1) * V_sig_H ; // composante 12 = V_1 . \sigma . V_2
// essai avec le mixte
// residu_local(3) = 0.5 * ((*ViB_3D)(1) * V_sig_H + V_sig_H * (*ViB_3D)(1)); // composante 12
#ifdef MISE_AU_POINT
if (Permet_affichage() > 8)
{cout << "\n en g_i: sig_HH * ViB(2)= " << V_sig_H
<< " df_ds_22= " << residu_3D(1) << " df_ds_21= " << residu_3D(3)
<< flush;
};
#endif
V_sig_H = sig_HH_3D * (*ViB_3D)(3); // = \sigma . V_3
residu_local(2) = V_sig_H * (*ViB_3D)(3); // composante 33 = V_3 . \sigma . V_3
#ifdef MISE_AU_POINT
if (Permet_affichage() > 8)
{cout << "\n en g_i: sig_HH * ViB(3)= " << V_sig_H << " df_ds_33= " << residu_3D(2);
// pour info on calcul la projection symétrique
double f_12 = (*ViB_3D)(1)* sig_HH_3D * (*ViB_3D)(2);
cout << " df_ds_12= " << f_12 << flush;
};
#endif
// calcul de la matrice tangente
// formule générique
//\frac{\partial \left( \sigma^{ij}(\epsilon'_{kl}) ~ V_{mj} \right )}{\partial \epsilon'_{rs}}
// = V_{mj}~\frac{\partial \sigma^{ij}(\epsilon_{kl}) }{\partial \epsilon_{op}}~V^r_{~.o}~ V^s_{~.p}
// ici on doit utiliser uniquement:
//\frac{\partial \left( \sigma^{3j}(\epsilon'_{kl}) ~ V_{3j} \right ) }{\partial \epsilon'_{rs}}
//\frac{\partial \left( \sigma^{\alpha j}(\epsilon'_{kl}) ~ V_{2j} \right ) }{\partial \epsilon'_{rs}}
// avec \alpha =1,2 ; (k,l) et (r,s) } = (2,2) ; (3,3) ; (1,2)
// d (V_{mi} ~\sigma^{ij} ~V_{nj} \right ) / d \epsilon'_{gh}
// = V_{mi}~V_{nj}~ d \sigma^{ij} / d \epsilon_{op} ~V^r_{~.o}~ V^s_{~.p}
// pour (m,n) et (g,h) = (2,2); (3,3); et (1,2)
derResidu_3D.Initialise(0.); // init
// for (int j=1;j<4;j++) // j représente la composante de l'axe de projection
// { for (int i=1;i<4;i++) // i représente le classement de (\epsilon'_{kl})
// // i= 1 pour (k,l) = (2,2), i=2 pour (k,l) = (3,3), i=3 pour (k,l) = (1,2)
// {// d (sig(2,j) V(2)(j)) / d (eps')(k,l)
// derResidu_3D(1,i) += (*ViB_3D)(2)(j) *
// (d_sigma_deps_3D(2,j,1,1) * (*ViH_3D)(indx(i))(1) * (*ViH_3D)(indy(i))(1)
// + d_sigma_deps_3D(2,j,1,2) * (*ViH_3D)(indx(i))(1) * (*ViH_3D)(indy(i))(2)
// + d_sigma_deps_3D(2,j,2,1) * (*ViH_3D)(indx(i))(2) * (*ViH_3D)(indy(i))(1)
// + d_sigma_deps_3D(2,j,3,3) * (*ViH_3D)(indx(i))(3) * (*ViH_3D)(indy(i))(3)
// );
// // d (sig(1,j) V(2)(j)) / d (eps')(k,l)
// derResidu_3D(3,i) += (*ViB_3D)(2)(j) *
// (d_sigma_deps_3D(1,j,1,1) * (*ViH_3D)(indx(i))(1) * (*ViH_3D)(indy(i))(1)
// + d_sigma_deps_3D(1,j,1,2) * (*ViH_3D)(indx(i))(1) * (*ViH_3D)(indy(i))(2)
// + d_sigma_deps_3D(1,j,2,1) * (*ViH_3D)(indx(i))(2) * (*ViH_3D)(indy(i))(1)
// + d_sigma_deps_3D(1,j,3,3) * (*ViH_3D)(indx(i))(3) * (*ViH_3D)(indy(i))(3)
// );
// // d (sig(3,j) V(2)(j)) / d (eps')(k,l)
// derResidu_3D(2,i) += (*ViB_3D)(3)(j) *
// (d_sigma_deps_3D(3,j,1,1) * (*ViH_3D)(indx(i))(1) * (*ViH_3D)(indy(i))(1)
// + d_sigma_deps_3D(3,j,1,2) * (*ViH_3D)(indx(i))(1) * (*ViH_3D)(indy(i))(2)
// + d_sigma_deps_3D(3,j,2,1) * (*ViH_3D)(indx(i))(2) * (*ViH_3D)(indy(i))(1)
// + d_sigma_deps_3D(3,j,3,3) * (*ViH_3D)(indx(i))(3) * (*ViH_3D)(indy(i))(3)
// );
//
// };
// };
// pour (m,n) et (g,h) = (2,2); (3,3); et (1,2)
// d (V_{mi} ~\sigma^{ij} ~V_{nj} ) / d\epsilon'_{gh} =
// (V_{mi} V_{nj}) (d \sigma^{ij} / d \epsilon_{op}) (V^g_{~.o} V^h_{~.p})
// pour 2,2 et 3,3 pas de pb de symétrie
// pour 1,2
for (int mn=1;mn<4;mn++) // mn représente la composante de l'axe de projection
{// d'abord les termes diagonaux
for (int gh=1;gh<3;gh++) // gh représente le classement de (\epsilon'_{gh})
derResidu_3D(mn,gh) = VhVg_BB(mn) && (d_sigma_deps_3D && VhVg_BB(gh));
// maintenant les termes non diagonaux à doubler
derResidu_3D(mn,3) = 2.* VhVg_BB(mn) && (d_sigma_deps_3D && VhVg_BB(3));
};
return derResidu_3D;
};
// limitation des variations d'épaisseurs et largeurs
// ramène true s'il y a eu une modif: sert pour la méthode 1
bool LoiContraintesPlanesDouble::Limitation_h_b(double& hsurh0, double& bsurb0)
{ bool retour = false;
// il ne faut pas que la nouvelle épaisseur soit inf au mini
if (hsurh0 < mini_hsurh0)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{if (Permet_affichage() > 3)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
{cout << "\n *** attention hsurh0 " << hsurh0 << ", est inferieur a la limite fixee "<<mini_hsurh0
<< " on le ramene a la limite ";
};
hsurh0 = mini_hsurh0;retour = true;
};
// il ne faut pas que la nouvelle largeur soit inf au mini
if (bsurb0 < mini_bsurb0)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{if (Permet_affichage() > 3)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
{cout << "\n *** attention bsurb0 " << bsurb0 << ", est inferieur a la limite fixee "<<mini_bsurb0
<< " on le ramene a la limite ";
};
bsurb0 = mini_bsurb0;retour = true;
};
// on regarde également pour les maxi
if (hsurh0 > maxi_hsurh0)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{if (Permet_affichage() > 3)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
{cout << "\n *** attention hsurh0 " << hsurh0 << ", est superieur a la limite fixee "<<maxi_hsurh0
<< " on le ramene a la limite ";
};
hsurh0 = maxi_hsurh0;retour = true;
};
if (bsurb0 > maxi_bsurb0)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{if (Permet_affichage() > 3)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
{cout << "\n *** attention bsurb0 " << bsurb0 << ", est superieur a la limite fixee "<<maxi_bsurb0
<< " on le ramene a la limite ";
};
bsurb0 = maxi_bsurb0;retour = true;
};
// retour
return retour;
};
// calcul de l'opérateur tangent (2ième méthode)
// d sigma^{ij} / d epsilon_{kl}
void LoiContraintesPlanesDouble::Calcul2_dsigma_deps(TenseurHHHH& d_sigma_deps)
{
#ifdef MISE_AU_POINT
if (Permet_affichage() > 8)
{ cout << "\n dans le repere locale d_sigma_deps_3D= \n";
int e=1;
for (int i=1;i<4;i++) for (int j=1;j<4;j++) for (int k=1;k<4;k++)for (int l=1;l<4;l++,e++)
{ cout << "("<<i<<","<<j<<","<<k<<","<<l<<")= "<<d_sigma_deps_3D(i,j,k,l) << " ; ";
if (e>6) {cout << "\n"; e=1;}
};
Tenseur3HHHH inter_HHHH;
d_sigma_deps_3D.ChangeBase(inter_HHHH,giB_tdt_3D);
cout << "\n dans le repere orthonormee d_sigma_deps_3D= \n";
e=1;
for (int i=1;i<4;i++) for (int j=1;j<4;j++) for (int k=1;k<4;k++)for (int l=1;l<4;l++,e++)
{ cout << "("<<i<<","<<j<<","<<k<<","<<l<<")= "<<inter_HHHH(i,j,k,l) << " ; ";
if (e>6) {cout << "\n"; e=1;}
};
};
#endif
// on calcule deux grandeurs intermédiaires
// a) d sigma^{ij} / d epsilon_{o'p'}(3D)~~V^1_{~.o'}~ V^1_{~.p'}
// on récupère que le pointeur, ce qui évite d'utiliser l'opérateur de copie du tenseur
// !!!!**** a vérifier -> oui c'est ok
// const Tenseur3BB & epsBB = *((Tenseur3BB*) &epsBB_tdt); // passage en dim 3
Tenseur3HH dSigDeps11P_HH(d_sigma_deps_3D && V1V1_BB);
#ifdef MISE_AU_POINT
if (Permet_affichage() > 7)
{ cout << "\n (*ViB_3D)(1)= "<<(*ViB_3D)(1);
cout << "\n V1V1_BB: "<< V1V1_BB
<< "\n dSigDeps11P_HH: "<< dSigDeps11P_HH
<< flush;
};
#endif
// b) d sigma^{ij} / d epsilon_{o'p'}(3D)~~V^g_{~.o'}~ V^h_{~.p'}
Tableau <Tenseur3HH* > dsigdeps_gh_P(3);
for (int gh=1;gh<4;gh++)
{dsigdeps_gh_P(gh) = ((Tenseur3HH*) &( d_sigma_deps_3D && VhVg_BB(gh)));
#ifdef MISE_AU_POINT
if (Permet_affichage() > 7)
{ cout << "\n VhVg_BB("<<gh<<")= "<< VhVg_BB(gh);
cout << "\n dsigdeps_gh_P(gh): "<< (*dsigdeps_gh_P(gh)) << flush;
}
#endif
};
// c) mat1(mn,gh) = [ V_{mi}~(d sigma^{ij} / d epsilon_{op}(3D)~~V^g_{~.o}~ V^h_{~.p}) ~~V_{nj}]
// SM(mn) = ( V_{mi}~~ d sigma^{ij} / d epsilon_{o'p'}(3D)~~V^1_{~.o'}~ V^1_{~.p'}~~V_{nj} )
Mat_pleine mat1(3,3);
Vecteur SM(3);
for (int mn=1;mn<4;mn++)
{SM(mn) = VhVg_BB(mn) && dSigDeps11P_HH;
for (int gh=1;gh<4;gh++)
mat1(mn,gh) = VhVg_BB(mn) && (*dsigdeps_gh_P(gh));
};
#ifdef MISE_AU_POINT
if (Permet_affichage() > 7)
{ cout << "\n mat1 "<< mat1;
cout << "\n SM: "<< SM << flush;
}
#endif
// d) calcul des variations des déformations transversales
// en fonction de la déformation longitudinale
// d epsilon'_gh / d epsilon'_11 =
// - [ V_{mi}~(d sigma^{ij} / d epsilon_{op}(3D)~~V^g_{~.o}~ V^h_{~.p}) ~~V_{nj}]^{-1}~
// * ( V_{mi}~~ d sigma^{ij} / d epsilon_{o'p'}(3D)~~V^1_{~.o'}~ V^1_{~.p'}~~V_{nj} )
Vecteur depsP_gh_depsP_11_BB(3);
depsP_gh_depsP_11_BB = - mat1.Inverse() * SM;
#ifdef MISE_AU_POINT
if (Permet_affichage() > 7)
{ cout << "\n depsP_gh_depsP_11_BB "<< depsP_gh_depsP_11_BB << flush;
}
#endif
// e) d sigma'^{11} / d epsilon'_11 (1D~CP) =
// = V_{1i}~V_{1j}~d sigma^{ij} / d epsilon_{op}(3D)~~V^1_{~.o}~ V^1_{~.p}
// + V_{1i}~V_{1j}~d sigma^{ij} / d epsilon_{op}(3D)~~~V^g_{~.o}~ V^h_{~.p}
// * d epsilon'_gh / d epsilon'_11
// ici on se sert du fait que l'on a déjà calculé:
// dSigDeps11P_HH = d sigma^{ij} / d epsilon_{op}(3D)~~V^1_{~.o}~ V^1_{~.p}
Tenseur1HHHH d_sigP11_depsP11_HHHH
( (V1V1_BB && dSigDeps11P_HH)
+ (V1V1_BB &&
( (*dsigdeps_gh_P(1))*depsP_gh_depsP_11_BB(1)
+ (*dsigdeps_gh_P(2))*depsP_gh_depsP_11_BB(2)
+ (*dsigdeps_gh_P(3))*depsP_gh_depsP_11_BB(3)
)
)
);
#ifdef MISE_AU_POINT
if (Permet_affichage() > 7)
cout << "\n d_sigP11_depsP11_HHHH "<< d_sigP11_depsP11_HHHH << flush;
#endif
// f) d_sigma^{ij} / d_epsilon_{kl} \approx
// V_1^{~.i}~V_1^{~.j}~d sigma'^{11} / d epsilon'_11 (1D~CP)~\gamma^{k}_{~.1} ~\gamma^{l}_{~.1}
//fonctions static définissant le produit tensoriel de deux vecteurs
// si les vecteurs sont égaux le tenseur est symétrique sinon il est non symétrique
TenseurHH & inter1_HH = Tenseur3HH::Prod_tensoriel((*ViH_3D)(1),(*ViH_3D)(1));
Coordonnee3H u_H(gamma3D.Colonne(1)); // qui correspond aux \gamma^{k}_{~.1}
TenseurHH & inter2_HH = Tenseur3HH::Prod_tensoriel(u_H,u_H);
d_sigma_deps = Tenseur3HHHH::Prod_tensoriel(inter1_HH,inter2_HH);
d_sigma_deps *= d_sigP11_depsP11_HHHH(1,1,1,1);
#ifdef MISE_AU_POINT
if (Permet_affichage() > 8)
{ Tenseur3HHHH inter_HHHH;
cout << "\n dans le repere locale d_sigma_deps= \n";
int e=1;
for (int i=1;i<4;i++) for (int j=1;j<4;j++) for (int k=1;k<4;k++)for (int l=1;l<4;l++,e++)
{ cout << "("<<i<<","<<j<<","<<k<<","<<l<<")= "<<d_sigma_deps(i,j,k,l) << " ; ";
if (e>6) {cout << "\n"; e=1;}
};
Tenseur3HHHH inter1_HHHH(d_sigma_deps);
inter1_HHHH.ChangeBase(inter_HHHH,giB_tdt_3D);
cout << "\n dans le repere orthonormee d_sigma_deps= \n";
e=1;
for (int i=1;i<4;i++) for (int j=1;j<4;j++) for (int k=1;k<4;k++)for (int l=1;l<4;l++,e++)
{ cout << "("<<i<<","<<j<<","<<k<<","<<l<<")= "<<inter_HHHH(i,j,k,l) << " ; ";
if (e>6) {cout << "\n"; e=1;}
};
};
#endif
};
// vérif des grandeurs et calcul des valeurs de l'épaisseur et de la largeur
// ainsi que calcul du jacobien final
// = 0 pas de modif
// = 1 modif due aux maxi permis de h et b, = 2 si def d'entrée incorrecte
// cas de la 2ième méthode
int LoiContraintesPlanesDouble::Calcul_et_Limitation2_h_b
(double& hsurh0,Tenseur3BB& eps_P_BB_3D, double& bsurb0,double& jacobien_tdt_3DD, const double & jacobien_0_3DD)
{ int retour = 0; // init par défaut
bool modif_hsurh0 = false; bool modif_bsurb0 = false; // init
// à partir de la valeur des déformations en cours on regarde les variations transversales correspondantes
// ici les déformations eps_BB_3D sont exprimées dans le repère V_a donc orthonormée
switch (type_de_deformation)
{case DEFORMATION_STANDART : case DEFORMATION_POUTRE_PLAQUE_STANDART :
// cas d'une déformation d'Almansi
{ // on est dans le repère V_a, V_3 est direction principale du tenseur de def -> formule classique
// epsBB33 = 1/2 * (1. -(h_0/h_tdt)**2) avec un maxi de o.5
if (eps_P_BB_3D(3,3) > 0.5)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{if (Permet_affichage() > 0)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
cout << "\n LoiContraintesPlanesDouble::Calcul_et_Limitation2_h_b(.."
<< "\n **erreur: la deformation d'epaisseur dans V_a "<<eps_P_BB_3D(3,3)
<< ") est superieur a 0.5, on ne peut pas s'en servir pour mettre a jour"
<< " la nouvelle epaisseur, on garde donc l'ancienne " << flush;
retour = 2; // là le pb est sérieux
}
else // sinon cas normal
{hsurh0 = 1./sqrt(1.-2.*eps_P_BB_3D(3,3));
// il ne faut pas que la nouvelle épaisseur soit inf au mini
if (hsurh0 < mini_hsurh0)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{if (Permet_affichage() > 3)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
{cout << "\n *** attention hsurh0 " << hsurh0 << ", est inferieur a la limite fixee "<<mini_hsurh0
<< " on le ramene a la limite ";
};
hsurh0 = mini_hsurh0; modif_hsurh0 = true;
retour = MaX(1,retour); // on signale qu'il y a eu un recadrage
};
// on regarde également pour les maxi
if (hsurh0 > maxi_hsurh0)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{if (Permet_affichage() > 3)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
{cout << "\n *** attention hsurh0 " << hsurh0 << ", est superieur a la limite fixee "<<maxi_hsurh0
<< " on le ramene a la limite ";
};
hsurh0 = maxi_hsurh0;modif_hsurh0 = true;
retour = MaX(1,retour); // on signale qu'il y a eu un recadrage
};
// si on a modifié via les extrèmes il faut mettre à jour la def
if (modif_hsurh0)
eps_BB_3D.Coor(3,3) = 0.5 * (1.-hsurh0 * hsurh0);
};
// -- suite seulement si retour est différent de 2 ---
// c-a-d que l'on n'a pas eu de pb sérieux sur h
if (retour < 2)
{// pour les deux autres def on n'est pas en directions principales mais c'est un plan principal
// on calcule les def principales via la technique de Mohr (possible car on est en repère ortho)
double epsP_c = 0.5 * (eps_P_BB_3D(1,1) + eps_P_BB_3D(2,2)); // le centre
double R = sqrt(Sqr(eps_P_BB_3D(2,2) - epsP_c) + Sqr(eps_P_BB_3D(1,2)));
double epsP_I = epsP_c + R; // la plus grande
double epsP_II = epsP_c - R; // la deuxième mais qui en toute rigueur pourrait être intercalé avec eps_P_BB_3D(3,3)
// d'où la variation de jacobien à condition de ne pas avoir de pb
bool pb = false;
// if ((epsP_I > 0.5) || (epsP_II > 0.5))
if (epsP_I > 0.5)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{if (Permet_affichage() > 0)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
cout << "\n LoiContraintesPlanesDouble::Calcul_et_Limitation2_h_b(.."
<< "\n **erreur: la deformation eps'_I "<<epsP_I
<< ") est superieur a 0.5, on ne peut pas s'en servir pour mettre a jour"
<< " la variation de volume " << flush;
epsP_I = 0.5-ConstMath::trespetit; // on le met au maxi
retour = 2; // là le pb est sérieux
pb = true;
};
if (epsP_II > 0.5)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{if (Permet_affichage() > 0)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
cout << "\n LoiContraintesPlanesDouble::Calcul_et_Limitation2_h_b(.."
<< "\n **erreur: la deformation eps'_II "<<epsP_I
<< ") est superieur a 0.5, on ne peut pas s'en servir pour mettre a jour"
<< " la variation de volume " << flush;
epsP_II = 0.5-ConstMath::trespetit; // on le met au maxi
retour = 2; // là le pb est sérieux
pb = true;
};
// on considère que la direction tranversale est celle correspondante à celle différente de
// eps_P_BB_3D(1,1) donc on en déduit une variation de b mais qui ne sera celle de V_2 que
// dans un cas isotrope, sinon c'est celle de la direction principale de def dans le plan !
if (retour < 2) // là encore on ne modifie que si retour est < 2
// car si retour == 2 on va obtenir des bsurb0 infinie
// eps_P_BB_3D(1,1), c'est la déformation cinématique imposée selon V1
// si elle est > eps_P_BB_3D(2,2), cela veut dire que epsP_I est de son côté
// par contre si c'est < , cela veut dire que c'est epsP_II qui est de son côté
// d'où le choix qui suit
{ if (eps_P_BB_3D(1,1) > eps_P_BB_3D(2,2))
{bsurb0 = 1./sqrt(1.-2.*epsP_II);}
else // cas d'une compression
{bsurb0 = 1./sqrt(1.-2.*epsP_I);};
};
// il ne faut pas que la nouvelle largeur soit inf au mini
if (bsurb0 < mini_bsurb0)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{if (Permet_affichage() > 3)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
{cout << "\n *** attention bsurb0 " << bsurb0 << ", est inferieur a la limite fixee "<<mini_bsurb0
<< " on le ramene a la limite ";
};
bsurb0 = mini_bsurb0; modif_bsurb0 = true;
retour = MaX(1,retour); // on signale qu'il y a eu un recadrage
};
if (bsurb0 > maxi_bsurb0)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{if (Permet_affichage() > 3)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
{cout << "\n *** attention bsurb0 " << bsurb0 << ", est superieur a la limite fixee "<<maxi_bsurb0
<< " on le ramene a la limite ";
};
bsurb0 = maxi_bsurb0;modif_bsurb0 = true;
retour = MaX(1,retour); // on signale qu'il y a eu un recadrage
};
// si on a modifié via les extrèmes il faut mettre à jour la def
// ici il s'agit uniquement de limiter eps_BB_3D(2,2) à des valeurs non absurdes donc on
// fait l'approximation (fausse en ortho mais cela ne devrait pas être un gros pb)
if (modif_bsurb0)
{eps_BB_3D.Coor(2,2) = 0.5 * (1.-bsurb0 * bsurb0);
// on est obligé de refaire une passe sur les def principales qui vont ^tre utilisées
// pour la variation de jacobien
epsP_c = 0.5 * (eps_P_BB_3D(1,1) + eps_P_BB_3D(2,2)); // le centre
R = sqrt(Sqr(eps_P_BB_3D(2,2) - epsP_c) + Sqr(eps_P_BB_3D(1,2)));
epsP_I = epsP_c + R; // la plus grande
epsP_II = epsP_c - R; // la deuxième mais qui en toute rigueur pourrait être intercalé avec eps_P_BB_3D(3,3)
// d'où la variation de jacobien à condition de ne pas avoir de pb
if (epsP_I > 0.5)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{if (Permet_affichage() > 0)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
cout << "\n LoiContraintesPlanesDouble::Calcul_et_Limitation2_h_b(.."
<< "\n **erreur: la deformation eps'_I (deuxieme passe!!) "<<epsP_I
<< ") est superieur a 0.5, on ne peut pas s'en servir pour mettre a jour"
<< " la variation de volume " << flush;
epsP_I = 0.5-ConstMath::trespetit; // on le met au maxi
retour = 2; // là le pb est sérieux
pb = true;
};
if (epsP_II > 0.5)
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{if (Permet_affichage() > 0)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
cout << "\n LoiContraintesPlanesDouble::Calcul_et_Limitation2_h_b(.."
<< "\n **erreur: la deformation eps'_II (deuxieme passe!!) "<<epsP_I
<< ") est superieur a 0.5, on ne peut pas s'en servir pour mettre a jour"
<< " la variation de volume " << flush;
epsP_II = 0.5-ConstMath::trespetit; // on le met au maxi
retour = 2; // là le pb est sérieux
pb = true;
};
};
if (!pb) // cas normal
{double var_vol = 1./sqrt((1-2.*eps_P_BB_3D(3,3))*(1-2.*epsP_I)*(1-2.*epsP_II));
// jacobien_0_3DD contient uniquement le jacobien_0_1D et on part par défaut
// avec ici h0 = b0 = 1, c-a-d ici dans la loi de comp: on ne calcule en fait que h/h0 et b/b0
// la véritable épaisseur et largeur est obtenue au niveau de l'élément fini
jacobien_tdt_3DD = jacobien_0_3DD * var_vol;
};
// sinon on laise en l'état
}; // fin du cas retour < 2
};
break;
// -- pour l'instant je laisse le log en attente
// case DEFORMATION_LOGARITHMIQUE : case DEF_CUMUL_CORROTATIONNEL : case DEFORMATION_CUMU_LOGARITHMIQUE :
// // cas d'une def logarithmique ou une approximation
// { hsurh0 = exp(eps_BB_3D(3,3));
// bsurb0 = exp(eps_BB_3D(2,2));
// };
// break;
default :
cout << "\nErreur : type de deformation qui n'est pas actuellement pris en compte, type= "
<< Nom_type_deformation(type_de_deformation);
cout << "\n LoiContraintesPlanesDouble::Limitation2_h_b(.. \n";
Sortie(1);
};
//---- debug
//jacobien_tdt_3DD = jacobien_0_3DD; // pour voir si c'est un pb de jacob ien
// retour
return retour;
};
// calcul de l'état final, dans le cas de la méthode de perturbation: version 2
void LoiContraintesPlanesDouble::Cal_avec_perturbation2()
{
// récup du conteneur spécifique
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
//récup de eps_mecaBB_t
// 1) dans le repère gi
Tenseur3BB& eps_mecaBB_t = *((Tenseur3BB*) (save_resul.eps_mecaBB_t));
// 2) dans le dernier repère V_a: a utiliser en sortie post, pas pour le calcul !!
Tenseur3BB& eps_P_mecaBB_t = *((Tenseur3BB*) (save_resul.eps_P_mecaBB_t));
Tenseur3BB& eps_P_mecaBB = *((Tenseur3BB*) (save_resul.eps_P_mecaBB));
// les inconnues sont les déformations mécaniques {epsilon'}_kl;
// (k,l) = (22) (33) (12) : on a donc 3 inconnues que l'on suppose rangées dans cette ordre
// c-a-d les def qui vont servir dans l'appel de la loi 3D
// contrairement au cas de la méthode de Newton, ici on utilise les grandeurs stockées à t
// sauf pour la cinématique
// --- initialisation des déformations mécaniques dans le repère V_a
eps_P_BB = eps_mecaBB_t; // là on est dans le repère g^i
// on passe dans le repère V_a qui est donc le nouveau repère, pour ce passage
eps_P_BB.ChBase(betaP_3D);
// on calcul la déformation cinématique dans la direction V_1
// \epsilon'_{11}(t+\Delta t) = \epsilon_{kl}(t+\Delta t) ~\gamma^k_{~.1} ~\gamma^l_{~.1}
// idem pour l'accroissement de déformation et la vitesse de déformation
eps_P_BB.Coor(1,1) = 0.; Deps_P_BB.Coor(1,1) = 0.; delta_eps_P_BB.Coor(1,1) = 0.;
eps_P_BB.Coor(1,3) = 0.; eps_P_BB.Coor(2,3) = 0.;
for (int k=1;k<4;k++) // calcul de la matrice gamma
for (int l=1;l<4;l++)
{eps_P_BB.Coor(1,1) += (*epsBB_tdt_cin)(k,l) * gamma3D(k,1) * gamma3D(l,1);
Deps_P_BB.Coor(1,1) += (*DepsBB_cin)(k,l) * gamma3D(k,1) * gamma3D(l,1);
delta_eps_P_BB.Coor(1,1) += (*delta_epsBB_cin)(k,l) * gamma3D(k,1) * gamma3D(l,1);
};
// vérif des grandeurs et calcul de la nouvelle épaisseur et largeur
double& hsurh0 = save_resul.hsurh0;double& bsurb0 = save_resul.bsurb0;
// la méthode va mettre à jour les valeurs bsurb0 et hsurh0 et du jacobien
bool erreur_sortie = false; // pour gestion de sortie d'infos supplémentaires
if (Calcul_et_Limitation2_h_b(hsurh0,eps_P_BB,bsurb0,jacobien_tdt_3D,jacobien_0_3D))
mise à jour 7.010
2023-05-03 17:23:49 +02:00
{if (Permet_affichage() > 3)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
cout << "\n LoiContraintesPlanesDouble::Cal_avec_perturbation2(... " << endl;
erreur_sortie = true;
};
if (erreur_sortie && (Permet_affichage() > 4))
{cout << "\n giH_tdt: " << *(umat_cont_3D->giH_tdt);
cout << "\n defBB: ";eps_BB_3D.Ecriture(cout);
cout << "\n DepsBB: ";Deps_BB_3D.Ecriture(cout);
};
// dans la direction V_1 c'est la déformation issue de la cinématique
// elle a déjà été calculée, et les autres composantes demeurent nulle
// donc eps_P_BB est au complet,
eps_BB_3D=eps_P_BB; // recopie dans eps_BB_3D: là on est dans le repère V_a
eps_BB_3D.ChBase(beta3D); // changement de base, là on est dans le repère g^i
// on s'occupe maintenant de l'incrément de déformation
delta_eps_BB_3D = eps_BB_3D - eps_mecaBB_t;
// au cas où, on impose (normalement c'est inutile) que les composantes 13 et 23 sont nulles
delta_eps_BB_3D.Coor(1,3) = 0.; delta_eps_BB_3D.Coor(2,3) = 0.;
// 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
#ifdef MISE_AU_POINT
if (unSurDeltat < 0)
{ cout << "\n le pas de temps est négatif !! "
<< "\n LoiContraintesPlanesDouble::Cal_avec_perturbation2(...";
};
#endif
unSurDeltat = ConstMath::tresgrand;
};
// ici on considère que la déformation est de type Almansi
Deps_BB_3D = delta_eps_BB_3D * unSurDeltat;
// calcul des invariants de déformation et de vitesse de déformation, et les def cumulées
// correspondant aux cas 3D
Calcul_invariants_et_def_cumul();
// --- on s'occupe des jacobiens
// à t=0 on considère que c'est le même jacobien que transmis
jacobien_0_3D = *(umat_cont_3D->jacobien_0);
// à tdt le jacobien est calculée à partir de la déformation supposée d'Almansi
// ici il s'agit du jacobien relatif au comportement mécanique
// |J| = sqrt(det(2 * eps_ij + g_ij(0)))
// calcul de la matrice correspondante à (2 * eps_ij + g_ij(0))
gij_meca_BB = (2. * eps_BB_3D + (*(umat_cont_3D->gijBB_0)));
gij_meca_BB.Matrice_composante(mat_inter);
// calcul du jacobien
jacobien_tdt_3D = sqrt(Dabs(mat_inter.Determinant()));
// initialisation du comportement tangent / au def
d_sigma_deps_3D.Inita(0.);
// appel du calcul de sig et dsig
bool en_base_orthonormee = false; // ici les tenseurs ne sont pas en orthonormee
// car on travaille dans la base naturelle
#ifdef MISE_AU_POINT
if (Permet_affichage() > 8)
{cout << "\n LoiContraintesPlanesDouble::Cal_avec_perturbation2(..."
<< "\n giH_tdt: " << *(umat_cont_3D->giH_tdt);
cout << "\n defBB: ";eps_BB_3D.Ecriture(cout);
cout << "\n DepsBB: ";Deps_BB_3D.Ecriture(cout);
};
#endif
lois_interne->Calcul_dsigma_deps (en_base_orthonormee, sig_HH_t_3D,Deps_BB_3D,eps_BB_3D
,delta_eps_BB_3D,jacobien_0_3D,jacobien_tdt_3D
,sig_HH_3D,d_sigma_deps_3D
,save_resul.l_energ,save_resul.l_energ_t,module_compressibilite_3D
,module_cisaillement_3D,*umat_cont_3D);
#ifdef MISE_AU_POINT
if (Permet_affichage() > 8)
{cout << "\n sig_HH_3D: " << sig_HH_3D;
cout << "\n d_sigma_deps_3D: "<< d_sigma_deps_3D << flush;
};
#endif
// sauvegarde des résultats
save_resul.def_P(1) = eps_P_BB(2,2); // récup de la solution dans le repère V_a
save_resul.def_P(2) = eps_P_BB(3,3); // "
save_resul.def_P(3) = eps_P_BB(1,2); // "
// on doit sauvegarder les def méca dans le repère g^i, ce sont directement
// les dernières composantes utilisées pour l'appel de la loi 3D
Tenseur3BB& eps_mecaBB = *((Tenseur3BB*) (save_resul.eps_mecaBB));
eps_mecaBB = eps_BB_3D; // en g^i
eps_P_mecaBB = eps_P_BB; // dans V_a pour le post traitement
};
// passage entre le calcul classique en 1D de Calcul_dsigma_deps et le calcul 3D
// méthode 2
void LoiContraintesPlanesDouble::Passage_Calcul_1D_3D_dsigma_deps
(bool en_base_orthonormee, TenseurHH & sigHH_t,TenseurBB& DepsBB
,TenseurBB & epsBB_tdt,TenseurBB & delta_epsBB,double& jacobien_0,double& jacobien
,TenseurHH& sigHH_tdt,TenseurHHHH& d_sigma_deps
,EnergieMeca & energ,const EnergieMeca &
,double& module_compressibilite,double& module_cisaillement
,const Met_abstraite::Umat_cont& ex
)
{
#ifdef MISE_AU_POINT
if (Permet_affichage() > 7)
cout << "\n LoiContraintesPlanesDouble::Passage_Calcul_1D_3D_dsigma_deps";
#endif
// récup du conteneur spécifique
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
module_compressibilite=module_cisaillement=0.; // init
energ.Inita(0.); // initialisation des énergies mises en jeux
// --- pour les contraintes passage en 3D
bool plusZero = true;
// on affecte les contraintes a t: en fait pour les lois incrémentales, les infos à t
// sont déjà stocké, donc cela ne sert à rien, mais pour une loi élastique par exemple ce n'est pas le cas
sig_HH_t_3D.Affectation_trans_dimension(*((Tenseur1HH*) &sigHH_t),plusZero);
// --- pour la cinématique
// passage des informations spécifique à la loi le_SaveResul
// initialisation du tenseurs contrainte
sig_HH_3D.Inita(0.);
sig_HH_3D.Affectation_trans_dimension(*((Tenseur1HH*) &sigHH_tdt),plusZero);
// passage des métriques de l'ordre 1 vers 3
Passage_metrique_ordre1_vers_3(ex);
// passage des informations liées à la déformation de 1 vers 3, et variation de volume
Tableau <TenseurBB *>* toto = NULL; // pour dire que ce n'est pas attribué
Vecteur* titi = NULL; // " " "
Passage_deformation_volume_ordre1_vers_3(DepsBB,epsBB_tdt,toto,delta_epsBB
,jacobien_0,jacobien,titi,*(ex.jacobien_t));
// -- construction de la base ortho associée à la direction de traction/compression
// ici il s'agit pour le premier vecteur de la direction 1,
// et compte tenu de Passage_metrique_ordre1_vers_3, le repère gi est orthogonale
// on va donc le normer
BaseB Vi_B(3); BaseH Vi_H(3);
BaseB * ViB = &Vi_B;BaseH * ViH = &Vi_H;
// tout d'abord on récupère la base gi actuelle et on la norme:
// on se sert de Vi_B comme stockage transitoire
// cout << "\n base giB(1) "<< (* ex.giB_tdt)(1)
// << "\n base giB(2) "<< (* ex.giB_tdt)(2);
#ifdef MISE_AU_POINT
if (Permet_affichage() > 8)
{ cout << "\n base giB_tdt_3D(1) "<< giB_tdt_3D(1)
<< "\n base giB_tdt_3D(2) "<< giB_tdt_3D(2)
<< "\n base giB_tdt_3D(2) "<< giB_tdt_3D(3) << flush;
};
#endif
Vi_B.CoordoB(1) = giB_tdt_3D(1);Vi_B.CoordoB(1).Normer();
Vi_B.CoordoB(2) = giB_tdt_3D(2);Vi_B.CoordoB(2).Normer();
Vi_B.CoordoB(3) = giB_tdt_3D(3);Vi_B.CoordoB(3).Normer();
// calcul des coordonnées locales de V_a en contravariants:
// V_a^{~i} = \vec V_a . \vec g^j
for (int a=1;a<4;a++) for (int i=1;i<4;i++)
Vi_H.CoordoH(a)(i) = Vi_B.CoordoB(a) * giH_tdt_3D(i);
// calcul des coordonnées contra
// Vi_B(a) = g_ij * Vi_H(a)
for (int a=1;a<4;a++)
Vi_B.CoordoB(a) = gijBB_t_3D * Vi_H(a);
// appel de la résolution 3D
// calcul des contraintes et ses variations par rapport aux déformations a t+dt
// ceci dans le cas où l'axe de traction simple est quelconque
// ici on suppose que:
// - tous les calculs sont en dim 3
// - l'axe 3 est normal aux deux autres (ceci pour Vi et pour gi)
//
// ViB et ViH : correspond à la base de traction telle que ViB(1) = ViH(1) = l'axe de traction
// La base Vi est supposée orthonormée
// 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
Calcul_dsigma_deps (en_base_orthonormee,ViB,ViH,sig_HH_t_3D,Deps_BB_3D
,eps_BB_3D,delta_eps_BB_3D,jacobien_0_3D,jacobien_tdt_3D
,sig_HH_3D,d_sigma_deps_3D,save_resul.l_energ,save_resul.l_energ_t
,module_compressibilite,module_cisaillement,*umat_cont_3D);
// passage des tenseurs résultats à l'ordre 1
((Tenseur1HH*) &sigHH_tdt)->Affectation_trans_dimension(sig_HH_3D,false);
// maintenant on renseigne le tenseur de sortie
bool pluszero = true;
d_sigma_deps.Affectation_trans_dimension(d_sigma_deps_3D, pluszero);
};
// passage entre le calcul classique en 1D de Calcul_DsigmaHH_tdt et le calcul 3D
// méthode 2
void LoiContraintesPlanesDouble::Passage_Calcul_1D_3D_dsigma_DsigmaHH_tdt
(TenseurHH& sigHH_t,TenseurBB& DepsBB,DdlElement & tab_ddl
,BaseB& giB_t,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& jacobien_0,double& jacobien
,Vecteur& d_jacobien_tdt,TenseurHH& sigHH_tdt,Tableau <TenseurHH *>& d_sigHH
,EnergieMeca & energ,const EnergieMeca &
,double& module_compressibilite,double& module_cisaillement
,const Met_abstraite::Impli& ex)
{
#ifdef MISE_AU_POINT
if (Permet_affichage() > 7)
cout << "\n LoiContraintesPlanesDouble::Passage_Calcul_1D_3D_dsigma_DsigmaHH_td";
#endif
// récup du conteneur spécifique
SaveResul_LoiContraintesPlanesDouble & save_resul = *((SaveResul_LoiContraintesPlanesDouble*) saveResul);
module_compressibilite=module_cisaillement=0.; // init
energ.Inita(0.); // initialisation des énergies mises en jeux
// --- pour les contraintes passage en 3D
bool plusZero = true;
// on affecte les contraintes a t: en fait pour les lois incrémentales, les infos à t
// sont déjà stocké, donc cela ne sert à rien, mais pour une loi élastique par exemple ce n'est pas le cas
sig_HH_t_3D.Affectation_trans_dimension(*((Tenseur1HH*) &sigHH_t),plusZero);
// --- pour la cinématique
// passage des informations spécifique à la loi le_SaveResul
// initialisation du tenseurs contrainte
sig_HH_3D.Inita(0.);
sig_HH_3D.Affectation_trans_dimension(*((Tenseur1HH*) &sigHH_tdt),plusZero);
// passage des métriques de l'ordre 1 vers 3
Passage_metrique_ordre1_vers_3(ex);
// passage des informations liées à la déformation de 1 vers 3, et variation de volume
Tableau <TenseurBB *>* toto = NULL; // pour dire que ce n'est pas attribué
Vecteur* titi = NULL; // " " "
Passage_deformation_volume_ordre1_vers_3(DepsBB,epsBB_tdt,toto,delta_epsBB
,jacobien_0,jacobien,titi,*(ex.jacobien_t));
// -- construction de la base ortho associée à la direction de traction/compression
// ici il s'agit pour le premier vecteur de la direction 1,
// et compte tenu de Passage_metrique_ordre1_vers_3, le repère gi est orthogonale
// on va donc le normer
BaseB Vi_B(3); BaseH Vi_H(3);
BaseB * ViB = &Vi_B;BaseH * ViH = &Vi_H;
// tout d'abord on récupère la base gi actuelle et on la norme:
// on se sert de Vi_B comme stockage transitoire
#ifdef MISE_AU_POINT
if (Permet_affichage() > 8)
{ cout << "\n base giB_tdt_3D(1) "<< giB_tdt_3D(1)
<< "\n base giB_tdt_3D(2) "<< giB_tdt_3D(2)
<< "\n base giB_tdt_3D(2) "<< giB_tdt_3D(3) << flush;
};
#endif
Vi_B.CoordoB(1) = giB_tdt_3D(1);Vi_B.CoordoB(1).Normer();
Vi_B.CoordoB(2) = giB_tdt_3D(2);Vi_B.CoordoB(2).Normer();
Vi_B.CoordoB(3) = giB_tdt_3D(3);Vi_B.CoordoB(3).Normer();
// calcul des coordonnées locales de V_a en contravariants:
// V_a^{~i} = \vec V_a . \vec g^j
for (int a=1;a<4;a++) for (int i=1;i<4;i++)
Vi_H.CoordoH(a)(i) = Vi_B.CoordoB(a) * giH_tdt_3D(i);
// calcul des coordonnées contra
// Vi_B(a) = g_ij * Vi_H(a)
for (int a=1;a<4;a++)
Vi_B.CoordoB(a) = gijBB_t_3D * Vi_H(a);
#ifdef MISE_AU_POINT
mise à jour 7.010
2023-05-03 17:23:49 +02:00
if (Permet_affichage() > 4)
ajout des lois de comportement
2021-09-23 11:21:15 +02:00
{cout << "\n LoiContraintesPlanesDouble::Passage_Calcul_1D_3D_dsigma_DsigmaHH_td( ";
cout << "\n base Vi_B(1) "<< Vi_B(1) << "\n base Vi_B(2) "<< Vi_B(2)
<< "\n base Vi_B(3) "<< Vi_B(3)
<< "\n base Vi_H(1) "<< Vi_H(1) << "\n base Vi_H(2) "<< Vi_H(2)
<< "\n base Vi_H(3) "<< Vi_H(3) << flush;
};
#endif
// appel de la résolution 3D
// calcul des contraintes et ses variations par rapport aux déformations a t+dt
// ceci dans le cas où l'axe de traction simple est quelconque
// ici on suppose que:
// - tous les calculs sont en dim 3
// - l'axe 3 est normal aux deux autres (ceci pour Vi et pour gi)
//
// ViB et ViH : correspond à la base de traction telle que ViB(1) = ViH(1) = l'axe de traction
// La base Vi est supposée orthonormée
// 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
bool en_base_orthonormee = false;
Calcul_dsigma_deps (en_base_orthonormee,ViB,ViH,sig_HH_t_3D,Deps_BB_3D
,eps_BB_3D,delta_eps_BB_3D,jacobien_0_3D,jacobien_tdt_3D
,sig_HH_3D,d_sigma_deps_3D,save_resul.l_energ,save_resul.l_energ_t
,module_compressibilite,module_cisaillement,*umat_cont_3D);
// passage des tenseurs résultats à l'ordre 1
((Tenseur1HH*) &sigHH_tdt)->Affectation_trans_dimension(sig_HH_3D,false);
// maintenant on renseigne le tenseur de sortie
bool pluszero = true;
// d_sigma / d_ddl = d_sigma / d_eps * d_eps / d_ddl
// ici l'opérateur tangent correspond à la sensibilité de d sigma'^11 / d eps'_11
// et la direction V_1 c'est la direction de g_1 (à la norme près donc
// d_sigma / d_eps normalement c'est == à d sigma'^11 / d eps'_11
// du coup ici on n'a pas à itérer sur toutes les déformations
int nb_ddlPlus1 = d_epsBB.Taille()+1;
for (int i_ddl=1;i_ddl<nb_ddlPlus1;i_ddl++)
d_sigHH(i_ddl)->Coor(1,1) = (d_sigma_deps_3D(1,1,1,1) * (*d_epsBB(i_ddl)))(1,1);
};
Reference in a new issue Copy permalink