Herezh++

remarque non bloquante :

il y a une faille dans la lecture des parametres "avec_parametres_de_reglage_" à partir de la ligne 566 du fichier Loi_ortho3D_entrainee.cc.

Si on met un seul paramètre et que ce paramètre est juste un mot-clé, il n'est pas pris en compte.
J'ai voulu mettre "seule_deviatorique". Ce qui se passe c'est que :
1) lecture d'un mot-clé ligne 584 (variable "nom" contient "seule_deviatorique"
2) passage dans le if ligne 595 car il y a une fin de ligne eof (variable "nom" prend la valeur "fin_parametres_reglage_" car je n'ai qu'un seul paramètre)
3) ne passe pas dans le if ligne 628
4) retour au while ligne 572
=> fin de boucle

c'est non bloquant car :
on peut activer "seule_deviatorique" si on met un autre paramètre (par exemple "sortie_post_ 0") ou bien un entier bidon après seule_deviatorique (par exemple : "seule_deviatorique 1")

test avec critère pli :
- ok en traction sur un élément

- ok aussi sur un gonflage de coussin carré
c'est ok au sens où la convergence est atteinte (avec un critère d'arrêt trs honorable de 5e-4).

mais avec cependant pas mal de messages de non convergence contrainte plane double en cours de calcul
extrait 1 :

LoiContraintesPlanesDouble::Calcul_et_Limitation2_h_b(..
 **erreur: la deformation eps'_I  0.56333369) est superieur a 0.5, on ne peut pas s'en servir pour mettre a jour la variation de volume 

extrait 2 :

**>> non convergence sur l'algo de la resolution de (sigma.V_2=0, sigma.V_3=0),
 mail: 1, ele= 30, pti=2 , nb_incr_total=6, nb_iter_total=36, val obtenues:
 eps_P(2,2)= 0.28840284,  eps_P(3,3)= -0.064309088,  eps_P(1,2)= -0.38861528
 precedentes: 0, 0, 0
 dernier_residu: 2249.8129, 2202.4613, -60.413879
 LoiContraintesPlanesDouble::Calcul_dsigma_deps (...: valeur final de  h/h_0(t+dt)= 1 et de  b/b_0(t+dt)= 1

J'ai aussi fait les mêmes modifs pour les lois:
Projection_anisotrope_3D , Loi_ortho2D_C_entrainee et Hypo_ortho3D_entrainee
qui utilise le même type de lecture, et qui donc devaient présenter le même type d'anomalie (je n'ai pas vérifié!)

concernant la non convergence CP/CP2 :

c'est surtout le fait d'avoir une déformation almansi supérieure à 0.5 qui me pose question. D'où vient ce dépassement qui normalement ne devrait jamais arriver cinématiquement parlant ?

Comment régler ce pb ?
- passer par une élongation ?
- borner avec une tangente hyperbolique ?

je pense que c'est normal dans un processus itératif de Newton.
- au début d'une nouvelle itération, les paramètres cibles (ici les def transversales ) sont updatés en fct de leurs précédentes valeurs + l'incrément. On peut très bien (en fonction de l'incrément) dépasser les seuils limites qui sont fixés par la physique.
En fait tout dépend du type de convergence. À chaque fois que l'on se trouve près d'une saturation, il y a toujours un risque de dépasser des seuils. C'est pour cela que j'ai mis le plus possible de paramètres de réglage et aussi de messages d'alerte.

Malheureusement je ne connais pas de méthodologie qui marche à tous les coups.
Par contre je suis bien conscient que ce n'est pas facile de se dépatouiller avec tous ces paramètres, en particulier il faut bien comprendre comment ils agissent pour espérer les manier de manière efficace.

j'ai l'impression que l'utilisation des def log, directement dans la loi de comp, devrait répondre à ton pb ?
a priori c'est ce que je prévois de faire à partir de mars.
D'ici là, peut-être que tu pourrais utiliser une loi hypo_ortho_3D ce qui te permettrait d'utiliser nativement des déformations cumulés log ?

Ce qui me parait étrange c'est que le module de compressibilité calculé dans la loi iso_ortho, est déterminé à l'aide d'une hypothèse log:
cf. source:

module_courant = untiers * sigBH.Trace() / (log_var_vol);

et ensuite la variation d'épaisseur est calculée à l'aide de ce module, en incrémental, en considérant une variation log !

houlà oui, tu as tout à fait raison pour la loi ortho_elas3D. C'est bien le bon module qui est remonté!! J'avais pas tracé eps33 via l'épaisseur pour cette loi et elle est conforme à DEF_EPAISSEUR.

Mais du coup, il faut appliquer ce même calcul à ISOELAS!! On voit bien que la variation d'épaisseur de cette loi est vraiment pas bonne.
Je pourrais en avoir besoin pour faire des lois non linéaires et compatible critère pli. Si ce n'est pas modifié, je peux toujours utiliser une loi orthoelas3D et la rendre isotrope, mais c'est un peu dommage.

Je ne peux pas utiliser les lois Hypo ou autre future loi Log.
Je suis obligé de passer par des lois Almansi pour avoir le critère pli.

dans ISOELAS3D, j'avais tenté un calcul de K via les valeurs propres de epsBH, puis conversion en def log. ça marchait pas mal. Mais j'ai remplacé tout ça par la méthode de orthoelas3D. ça donne pareil en plus court dans le code.
Pour info, j'ai à peu près trouvé la raison de la divergence de mon calcul en traction. En fait, pour se comporter comme une loi log, ça me conduit à avoir un nu "almansi" assez exotique. Assez tôt dans le calcul, j'ai des termes négatifs dans l'opérateur tangent d_sigma_depsHHHH. Herezh passe par le calcul de 2 coefs pour calculer le tenseur sigBH et l'opérateur tangent : coef1 et coef2.
Par exemple :
coef1 : -1025.3953
coef2 : 160.6144
d_sigma_depsHHHH 1111 : -4131.8385
d_sigma_depsHHHH 1122 : -15019.656
d_sigma_depsHHHH 1133 : -3659.5756
d_sigma_depsHHHH 2222 : -23263.207
d_sigma_depsHHHH 2233 : -6820.8096
d_sigma_depsHHHH 3333 : -1402.0724
d_sigma_depsHHHH 1212 : 723.92779
d_sigma_depsHHHH 1313 : 174.96021
d_sigma_depsHHHH 2323 : 576.88192

Le plus drôle, c'est que ça ne dérange pas trop la convergence. Le calcul avance et calcule comme il faut (ce qui donne le graphe précédent). Tout semble ok jusqu'à ce que le nu "almansi" se rapproche de 1 et que l'on ait peu à peu coef1 == coef2. Quand ces 2 coefs sont proches, il y a des termes assez petits sur la diagonale de d_sigma_depsHHHH.

C'est vraiment curieux comme histoire. bon un peu hors sujet du coup.

pour en revenir à orthoelas3D, j'ai regardé après le calcul du "module_compressibilite", comment était calculé le module de cisaillement.

dans le code il y a ce calcul de G moyen :

if ((cas_calcul == 0) || (cas_calcul == 1))
      {module_cisaillement = untiers * 2. * (G12+G13+G23);}

il y aurait pas un 2 en trop ?

en compression, log_var_vol est négatif donc toujours inférieur à "petit" et ne passe jamais dans ce "if".

Ah oui, je suis d'accord, j'ai modifié !!
- dans toutes les lois anisotropes, dans la loi Iso_elas_expo3D
- dans les lois de Maxwell j'ai fait une modif du même genre

version 7.039

Pour les lois isoelas classique type Hooke je conserve la compressibilité classique pour l'instant (c-a-d celle qui est définie en petite def) car isoelas linéaire avec Almansi est clairement une loi qui n'a pas pour vocation de modéliser les grandes déformations. Lorsque j'aurai la version log avec un opérateur tangent efficace (ce qui n'est pas le cas actuellement, mais ce sera le cas courant 2025 normalement) je pourrai mettre en log un calcul avec les volumes. C'est un peu incohérent avec les lois anisotropes, mais pour l'instant je laisse comme cela.

NB:J'avais introduit à un moment le calcul avec variation de volume et cela a posé des pb (je ne me rappelle plus lesquels) du coup je suis revenu à la def classique:
Pour mémoire j'ai écrit dans le code:

// -- calcul des modules
// 22 fev 2019 ** non, en fait on a tendance à utiliser le module sécant du coup, c'est le module de compressibilité
//  habituelle qui correspond, donc on revient à la forme générique cf. doc,
//    // on n'utilise plus la forme linéaire, mais à la place la variation relative de volume
//    // constaté, au sens d'une mesure logarithmique, sauf dans le cas où cette variation est trop petite
//    //module_compressibilite = E/(3.(1.-2.*nu));
//    //  calcul de la valeur de la variation relative de volume en log
//    double log_var_vol = log(()/((ex.jacobien_0)));
//    if (log_var_vol > ConstMath::petit)
//      {module_compressibilite = sigBH.Trace() * untiers / log_var_vol;}
//    else // si la variation de volume est trop faible on passe par la formule traditionnelle

sources version 7.039, mis à jour sur git

en fait ce calcul de K via ln(V) pose problème si ce K est utilisé pour autre chose que la variation d'épaisseur. Je pense au calcul du pas de temps critique par exemple, ou bien peut-être le contact (qui requiert un module K je crois ?), c'est-à-dire en général tout endroit du code qui demande une caractéristique mécanique réelle.
Il y a des cas de chargement pour lesquels le calcul K= trace(sig)/(3ln(V)) ne donne pas la caractéristique réelle du matériau. Par exemple, si sigBH a une trace nulle (cas du cisaillement simple), on obtient K = 0. Ce qui n'est pas la réalité du matériau, mais pour la variation d'épaisseur, c'est exactement ce qu'il faut renvoyer pour obtenir une non-variation d'épaisseur en cisaillement simple.

C'est un peu comme si ce calcul de K via ln(V) devait être stocké dans une variable séparée, bien distincte avec un nom qui va bien. Et ensuite, utiliser cette variable uniquement pour la variation d'épaisseur.

Pour ce qui est de laisser ISOELAS en l'état, pourquoi pas, même si je ne comprends pas ce choix.
Pas de problème puisque je peux utiliser ortho_elas3D pour faire une loi isotrope.

C'est un peu comme si ce calcul de K via ln(V) devait être stocké dans une variable séparée

oui, et il y a aussi le choix entre K tangent et/ou K sécant, et... le calcul d'un G matériau pour l'intégration dans le calcul du pas critique: c'est dans un coin de ma tête
bon... en attente

j'ai poursuivi l'écriture de ma loi ORTHO_ELAS3D avec la nouvelle version.

J'obtiens exactement ce que je veux en traction (contrainte non linéaire, coef poisson constant au sens log dans la largeur et l'épaisseur).
Mais ça ne passe pas en compression à cause de la variation d'épaisseur (module 3K inférieur à delta trace de sigma).

à propos de ce bout de code du .cc de la loi ortho_elas :

double log_var_vol = log((*(ex.jacobien_tdt))/(*(ex.jacobien_0)));
    // pour le module de compressibilité, choix entre les différents cas
    if ((cas_calcul == 0) || (cas_calcul == 2))
     {if (Dabs(log_var_vol) > ConstMath::petit)
        {module_compressibilite = untiers * sigBH.Trace() / (log_var_vol);}

il s'agit bien là du module sécant Ks ?

Tandis que la variation d'épaisseur requiert un module tangent Kt avec par exemple cet extrait du .cc membrane quadrangle (quadramemb) :

{if (atdt)
      {double epaisseur_moy_tdt = (epaisseur_moy_t * jacobien_moy_ini) / jacobien_moy_fin
                               * troisK_moy / (troisK_moy - traceSig_moy+traceSig_moy_ini);

je me demande :
1) pourquoi est-ce que ça marche quand même en traction ?
2) pourquoi ça passe pas en compression ?

je viens de comprendre pour mon pb en compression.

Il reste ma question qui concerne traction et compression :
pourquoi ça fonctionne avec un Ks alors que la méthode de variation d'épaisseur dans quadramemb demande un Kt ?