Herezh++

Je viens de regarder...
1) pour utiliser une loi anisotrope il est nécessaire de définir un repère d'anisotropie, ce qui n'est pas le cas dans les deux fichiers de calcul.
2) en implicite je pense qu'il y aura un pb avec les membranes, si elles ne sont pas stabilisées dans la direction normale.

"je n'ai pas le réflexe de relancer avec HZpp en cas de bug. C'est dans un coin de ma tête maintenant."
oui, il faut utiliser la version avec vérification tant qu'il y a des pb en mise en données, et même ensuite, c'est toujours mieux de démarrer une exécution avec vérification pour vérifier qu'il n'y a pas de message de pb.
De même s'il arrive un pb pendant l'exécution, c'est intéressant d'essayer un restart avec vérification, juste avant l'arrivée du pb car les .BI sont indépendant de la version avec ou sans vérification.

"pour l'instant le calcul RD avec loi PEL60 + anisotropie a l'air bien plus lent que la loi PEL60 isotrope."
Concernant la vitesse, normalement en anisotrope, c'est plus lent quand isotrope, car il faut calculer la convection du repère puis la dérivée de cette convection par rapport à la déformation ou directement les ddl (au choix). Donc tout ça représente pas mal de calcul (on peut s'en rendre compte via la documentation théorique).

"Le calcul implicite avec loi PEL60 3D isotrope était à peine plus rapide que la loi anisotrope, avec une convergence très similaire."
Il faut regarder le temps cpu spécifique à la loi de comp.

"Je me pose la question de la symétrie matrice dans le processus Newton de contrainte plane. Est-ce que le solveur dédié est forcément symétrique ou y a-t-il possibilité de le rendre non symétrique ?"
La loi CP avec Newton utilise une matrice 1x1, donc ici le pb de la symétrie ne se pose pas (cf. théorie)
Par contre le calcul de l'opérateur tangent final, lui, utilise certaines symétrie du comportement initial.

Finalement, je ne comprends pas bien pourquoi avec une loi élastique, tu obtiens un opérateur final orthotrope en 3D non symétrique ...
J'ai re-regarder les développements théoriques (cf. manuel théorique)et normalement on doit récupérer un opérateur symétrique ???
Peux-tu re-regarder ? (je parle avant les CP)

"""""
"Le calcul implicite avec loi PEL60 3D isotrope était à peine plus rapide que la loi anisotrope, avec une convergence très similaire."
Il faut regarder le temps cpu spécifique à la loi de comp.
"""""
Désolé je n'ai pas été assez précis. Je ne parle pas de temps cpu. Je parlais en terme de nombre d'itérations. Seulement quelques itérations de plus pour le cas anisotrope.

"""""
La loi CP avec Newton utilise une matrice 1x1, donc ici le pb de la symétrie ne se pose pas (cf. théorie)
"""""
ha, je ne savais pas. J'étais persuadé que cela utilisait le comportement tangent du matériau pour faire progresser la méthode de Newton. Donc c'est pas ça alors.

"""""
Finalement, je ne comprends pas bien pourquoi avec une loi élastique, tu obtiens un opérateur final orthotrope en 3D non symétrique ...
J'ai re-regarder les développements théoriques (cf. manuel théorique) et normalement on doit récupérer un opérateur symétrique ???
Peux-tu re-regarder ? (je parle avant les CP)
"""""
J'ai remarqué ça en effectuant les développements théoriques pour mettre en place la méthode d'identification. J'en ai déjà rédigé une partie sous CVS. Mais ce document ne montre pas la dissymétrie. J'ai constaté cela via le script que j'ai mis en place peu à peu pour étudier la forme et les valeurs de la matrice d'anistropie (expression du tenseur H sous forme matricielle, voir doc CVS). Quand je multiplie la matrice H par la matrice C du comportement élastique, j'obtiens une matrice non symétrique. Je pourrais te le transmettre si tu veux.

quoiqu'il en soit, je pense que la dissymétrie doit être le cas quelque soit la loi, notamment PEL60. Et là, j'en suis sûr non pas par la théorie. Je le vois via les calculs Herezh. Sans l'option de matrice non symétrique, les calculs ne convergent pas avec une loi PEL60 anisotrope en 3D sur éléments HEXAEDRE. Avec l'option, ça converge très bien, voire aussi bien que le cas isotrope en terme de nombre d'itérations. Au niveau temps de calcul, a priori, la loi anisotrope demande par contre le double en temps cpu.

Je ne mets jamais la loi PEL60 sur redmine par une sorte de rigueur de confidentialité. C'est pourquoi je passe par un exemple élastique. Et au final, le problème de convergence pointé par ce ticket n'a pas l'air de dépendre de la loi, ni même des valeurs des coefficients d'anistropie.

Si on a tous les paramètres d'anisotropie égaux à 1, le calcul souffre malgré tout d'un problème de convergence. JE viens de le voir. Je ne l'avais pas remarqué avant que tu corriges le bug de lancement. Le calcul RD que je t'ai transmis précédemment devrait avoir la même convergence (en nombre d'itérations) que le cas isotrope que je t'envoie ici

Gérard,

suite à la réunion de ce matin, j'avais dit que A3 avait une légère influence. Mais en fait, A3 n'a aucune influence quand on est en contrainte plane. C'est ce que j'ai l'air d'obtenir dans mes simulations (traction, traction biaxiale). La légère influence que j'ai mentionnée ce matin était due à une erreur dans mes conditions limites de traction biaxiale conduisant à sig33 non nul.

je pense que ce résultat est intéressant à savoir pour tes tests sur membrane. Normalement, tu ne devrais pas voir d'influence de A3.

un peu de détail ci-dessous :

La seule preuve théorique que j'ai, c'est en élasticité linéaire via un logiciel de calcul formel pour état de contrainte quelconque ou un état de déformation quelconque.

Je raisonne en 3D avec un axe O_3 normal aux axes globaux I_1,I_2 du repère global et un angle entre l'axe O_1 et l'axe I_1. Mais pour une membrane, le repère d'anisotropie est toujours défini ainsi. Il est tel que O_3 est normal à la membrane et on peut définir un angle entre la direction O_1 et l'axe global I_1 (projeté sur la membrane, le repère ad hoc comme tu l'appelles dans Herezh).

J'ai ça pour un angle de 0° :
en déformation (epsilon = [H:C]^-1 . sigma) :
eps11 = -NU*sig33/(A3*E) - NU*sig22/(A2*E) + sig11/(A1*E)
eps22 = -NU*sig33/(A3*E) + sig22/(A2*E) - NU*sig11/(A1*E)
eps33 = sig33/(A3*E) - NU*sig22/(A2*E) - NU*sig11/(A1*E)
eps12 = sig12/(2*B12*G)
eps13 = sig13/(2*B13*G)
eps23 = sig23/(2*B23*G)

donc en ce qui concerne l'état de déformation, si sig33=0, alors epsij est indépendant de A3. Et en particulier, eps33 qui régit l'épaisseur.

pour les contraintes (sigma = H:C:epsilon):
sig11 = A1*E*(-NU*eps22 - NU*eps33 + eps11*(NU - 1))/(2*NU**2 + NU - 1)
sig22 = A2*E*(-NU*eps11 - NU*eps33 + eps22*(NU - 1))/(2*NU**2 + NU - 1)
sig33 = A3*E*(-NU*eps11 - NU*eps22 + eps33*(NU - 1))/(2*NU**2 + NU - 1)
sig12 = 2*B12*G*eps12
sig13 = 2*B13*G*eps13
sig23 = 2*B23*G*eps23

là encore, pas de A3 ailleurs que dans sig33 (qui sera nul en contrainte plane).

J'ai également regardé pour un angle non nul entre I_1 et O_1. Le résultat de calcul formel est barbare donc je le mets pas ici. Mais en pièce jointe, j'ai mis l'exemple d'un angle de 45°. En faisant une recherche de caractere sur A3, on voit bien que sig33=0 annulera l'influence de A3 sur eps11, eps22 et eps33.

J'ai essayé avec un angle moins classique (angle = 12° par exemple). Le résultat est très très barbare. Même avec une recherche de caractère, on n'y comprend rien.

Par la simulation en tout cas, je ne constate pas de différence entre A3=1 et A3=0.1, quelque soit l'angle <I_1, O_1>, quand sig33=0, pour une loi PEL60.

1) Concernant la convergence en RD, je viens de faire une vérification (version avec vérifications) sur un incrément de temps (0.1s) et j'obtiens exactement le même nombre d'itérations: 1016
Pour info, au niveau des temps de calcul, on a presque un rapport de 5 avec les paramètres CP que tu utilises, le codage en dur de la loi CP étant évidemment très performant !

Si on réduit la précision demandée pour les CP:

NEWTON_LOCAL  avec_parametres_de_reglage_
                 nb_iteration_maxi_  30
                 nb_dichotomie_maxi_  5
                 tolerance_residu_ 1.e-3#6
                 tolerance_residu_rel_ 1.e-2#5
              fin_parametres_reglage_Algo_Newton_

Là j'obtiens un rapport de 4.

2) Concernant la symétrie de l'opérateur tangent, en fait ma réponse n'était pas tout à fait claire:
- l'opérateur de projection conduit à un résultat symétrique (le tenseur contrainte), c'est de cet opérateur que je parlais en disant que théoriquement il était symétrique,
- l'opérateur tangent (classique) d sig(projeté)/d eps n'est effectivement pas symétrique, ce qui fait qu'en implicite la convergence est bien meilleurs avec un stockage non symétrique. En réfléchissant, je ne vois pas de contre-indication à ce fait: localement cela signifie que l'on n'a pas un potentiel associé convexe dans le cas d'une loi élastique, mais d'un point de vue matériau ce n'est pas obligatoire d'avoir une convexité.
Bon... à voir

Bonjour Julien
Le calcul du bulge test en membrane que tu proposes n'est pas trivial en implicite.
En simple élasticité, en implicite, j'ai réussi à le faire passer en utilisant des membranes stabilisées par de la flexion via des éléments coques sfe1 à 2 pti dans l'épaisseur.
Je te joins ma mise en données.
Concernant la projection anisotrope + contrainte plane, peut-être vaut-il mieux essayer dans une première étape de valider la loi dans un essai de chargement dans le plan de la membrane. Une fois que cela sera ok, on pourra passer à des déplacements hors plans.

ci-joint un cas test en cisaillement simple d'une plaque 100x100mm constituée de 9 éléments quadrangle.

à vrai dire, le calcul est en train de tourner et je ne vois pas pour l'instant d'évolution anormale de la convergence. Le premier incrément a convergé. C'est lent malgré le petit maillage mais ça semble converger

la pièce jointe précédente ne contient pas tout ce qu'il faut.

voici le cas test avec tous les fichiers nécessaires

syndrome d'oubli de pièce jointe. Voici le cas test bulge

Bonjour Julien,
en travaillant sur le cas du cisaillement, on utilisant des lois de complexités croissantes, j'ai vu qu'un des problèmes était le calcul de la compressibilité, lorsque la variation de volume est très faible. Or ce module est utilisé pour les loi CP ce qui induit des instabilités
Jusqu'à présent, j'avais utilisé une limite inférieure arbitraire à la variation log de volume pour statuer sur le type de calcul à employer pour la compressibilité c-a-d:
si var vol < val_min ==> utilisation de la compressibilité initiale sinon : calcul via pression et var vol
bon... j'ai essayer une autre technique, qui consiste à comparer la compressibilité nouvelle avec celle initiale et si diff trop grande, adoption de la compress initiale.
Avec un paramètre par défaut qui peut être modifié éventuellement par l'utilisateur:
là cela fonction beaucoup mieux .... avec une loi ortho3D projetée
@suivre !