Herez++_test

Gérard,

ci-joint une version du calcul axi précédent sans axisymétrie (membrane QUADRANGLE LINEAIRE) avec les mêmes lois. Il y a un peu le même souci. Pas de masses nulles mais un gros problème de loi contrainte plane (qui cherche à établir une épaisseur négative). Donc déjà pas sûr que ça vienne de l'axisymétrie.

J'ai également ajouté une loi mooney-rivlin. Pas l'air de poser problème. Donc a priori pas de problème de fond sur les lois en invariant de B.

A remarquer que la loi MAHEO_HYPER et donc mon potentiel HYPER W est singulier au départ (paramètre de régularisation à 1.e-12). Est-ce que ça viendrait de là ? j'ai essayé plusieurs valeurs de regul à la hausse ou à la baisse => pas d'amélioration

ps : désolé, j'ai mis ce ticket dans Herez++_test sans faire attention.

a+
Julien

Julien,
avec la loi MAT__ j'obtiens le message suivant:

" ======================================================================
INCREMENT DE CHARGE : 1 intensite 0.01 t= 0.01 dt= 0.01 ======================================================================
recalcul de la pseudo-masse (type: 1)
Erreur : (J_r(1)-3.)= -4.44089e-16 est < 0. ! ce n'est pas possible avec les formules actuelles
Maheo_hyper::Calcul_dsigma_deps

passage dans la methode Sortie

erreur de loi de comportement element= 96 point d'integration= 3
"

Donc là je pense qu'il s'agit d'un pb près de 0.
Avec un facteur de régularisation de e-5 par exemple, c'est ok, mais ensuite vient un pb sur le calcul de l'épaisseur:

" *** pb dans le calcul de la nouvelle epaisseur: le coef de compressibilite moyen = 50880.045 est par rapport a delta la trace de sigma = 838630.29 10 fois inferieur !! ce n'est pas normal a priori on continue sans changer l'epaisseur ...
QuadraMemb::CalEpaisseurMoyenne_et_vol_pti(...)
"

bon... je suis pris cet après-midi,
je regarde demain matin

ok merci.

Au fait, je me disais pourquoi pas essayer avec une matrice masse physique classique. ça serait moins performant, mais juste pour voir si ça tourne. Mais :
1) pas moyen de trouver comment utilise les masses physiques en mode dynamique_relaxation_dynam

2) j'ai essayé en dynamique explicite avec une loi élastique => ça tourne mais pas d'amortissement cinétique

• erreur, pour l'instant le module d'Young n'est pas calcule !!
  Hyper_externe_W::Module_young_equivalent(...

4) et avec MAHEO_HYPER (MAT__), le calcul tourne mais la structure bouge pas beaucoup (pas de temps très petit à cause coef K).

ci-joint un .info en dynamique explicite si besoin

décidement...
Ma priorité reste un calcul dynamique_relaxation_dynam et si possible l'option qui va bien pour avoir les masses physiques.

a+
Julien

Gérard,

bon que des bonnes nouvelles :
1) les lois HYPER_EXTERNE_W ont l'air de tourner en relaxation dynamique si je supprime la partie volumique du potentiel (en mettant C__K3 = 0 dans le fichier de constantes utilisateur) et que l'additionne à une loi ISOHYPERBULK

conclusion : c'est bien la partie volumique de HYPER_EXTERNE_W qui pose problème

2) les calculs de dérivées sont bons dans la partie volumique dans mooney-rivlin du fichier source MooneyRivlin3D.cc ligne 1213
=> r.a.s de ce côté concernant les sources herezh et donc pour toutes les lois genre poly, hart-smith, etc...
(et j'avais une erreur de constante dans le tout 1er fichier potentiel de ce ticket => à écraser par la nouvelle version ci-jointe)

je joins ici des nouveaux fichiers pour pouvoir tester.
J'ai clarifié les noms de loi et de potentiel. Il n'y a plus besoin de changer de potentiel pour changer de loi. Cela se fait directement en changeant de nom de loi (cf la rubrique materiaux du .info ci-joint pour switcher de l'un à l'autre) sachant que :
1) j'ai mis "Natif" quand c'est une loi native Herezh, sinon c'est "HW" quand c'est une version HYPER_EXTERNE_W
2) quand le mot ISOHYPERBULK apparait, ça veut dire que c'est une loi de mélange avec un ISOHYPERBULK natif herezh
==> ATTENTION A BIEN METTRE C__K3 = 0 dans le fichier de constantes utilisateur dans le cas d'une loi additive ISOHYPERBULK + HYPER_EXTERNE_W (et inversement remettre C__K3 = 30000 si la loi n'est pas additive)

ps : c'est un peu bizarre car du coup, la relaxation dynamique fonctionne avec la loi Maheo hyper en version HYPER_EXTERNE_W additionné à ISOHYPERBULK (nom de loi = MAT_HW_Maheo_modifie_ISOHYPERBULK) alors que la loi native MAT_Natif_Maheo_ISOHYPERBULK rencontre toujours un problème comme vu précédemment.

tant mieux pour moi mais donc il y a quelque chose à vérifier dans MAHEO_HYPER. Je vais regarder de mon côté les dérivées dans le code source. Mais j'ai vérifié toutes les lois en traction/compression : à paramètres identiques les lois de même type native, HW additif ou non sont strictement identiques sur une courbe contrainte-déformation. Donc j'y crois pas trop à une erreur de formules de dérivées. à suivre..

bon, déjà, ça fonctionne avec le potentiel HYPER_W, du coup ça devient moins pressant comme demande...

a+
Julien

Julien,
j'ai retravaillé les développements théoriques et j'ai trouvé des bugs (je ne me rappelle plus où) .
Je te joins les nouvelles formules. Il faut regarder les relations:
- 173 (n'ont pas changé)
- 190, 191, 192 et 193 (les deux premières et la dernière sont nouvelles )
J'ai aussi modifié le code en conséquence et j'ai remplacé les limites (Maxi) par des développements limités, car à la suite d'un petit test sur octave j'ai constaté que les maxi ne marchait pas bien pour V proche de 1. je joins la partie modifiée dans le texte:

Je suis intéressé que tu regardes ... si tu as un moment

Gérard,

j'ai regardé formules et code. Tout m'a l'air bon dans les dérivées.

Pour le module sécant, les formules 192 sont en accord avec 189 et le code aussi. Mais je ne suis pas sûr d'avoir vraiment compris cette notion de module sécant. Si je comprends bien, le module sécant est la relation entre une variation de pression et une variation de V (V=v/v0). Si non, alors le reste ci-dessous est à mettre à la poubelle :)

mais admettons, si ce module Ks est la relation entre delta(pression) et delta(V), alors ça revient en fait à dire que :
Ks = -d(p)/d(V)
Et du coup, j'ai du mal à comprendre pourquoi il faut postuler une forme pour Ks ( -p = Ks log(V) ) afin d'établir la formule 189.

par exemple, pour les potentiel w (volume actuel) :

si on prend directement la formule 18, ça donne pour un potentiel w :
Ks = -d(p)/d(V) = V.d^2(w)/d(V)^2 + 2*d(w)/d(w)

si on fait un calcul test sur un cube, cette formule pour Ks dédiée au potentiel w est exacte si :
1) on applique un chargement quelconque (traction, oedmétrique, un noeud sur lequel on applique UX,UY,UZ quelconque, etc...)
donc c'est a priori assez général.
2) on prend un loi qui fait la séparation exacte entre variation de volume et variation de forme (donc autre chose que Favier)

donc en gros, cette formule 18 n'est pas juste dédiée aux sollicitations purement sphériques. Du moment que l'on a une loi qui fait bien la séparation, ça marche tout le temps et le module sécant Ks que l'on peut en déduire correspond très bien au résultat donné par Herezh si on calcule la dérivée numérique : Ks = delta Spherique_sig / delta (VOLUME_PTI/volume initial du pti)

si on prend une loi type Favier, ça n'est plus exact car la partie déviatorique devrait entrer en jeu puisqu'elle génère de la variation de volume.

je ne sais pas si c'est utile comme remarque, si ça peut améliorer la convergence, ou si c'est pouillème par rapport à ce qui est déjà implanté.

Concernant les potentiels en invariants J3, je suis donc parti sur l'équation de base (19) => -p = d(W)/d(V)
et j'ai regardé les 4 types de potentiels.

si on suppose Ks = -d(p)/d(V), alors => Ks = d^2(W)/d(V)^2

j'ai regardé le type 2 W = K/2 (V-1) => alors je ne sais pas si cela avait été remarqué avant mais il n'est pas du tout utilisable et je ne l'ai d'ailleurs jamais utilisé. Comme sa dérivée première donne K/2, donc une constante, on obtient une pression hydro constante égale à K/2, même au repos ====> à enlever d'Herezh, non ?

j'ai essayé le type 3 (W = K/2 log(V)**2) => je trouve Ks = d^2(W)/d(V)^2 = K*(1.-log(V))/V**2. Et cette formule correspond bien à la dérivée numérique calculée à partir du .maple de Herezh. Et ça marche pour n'importe quel type de sollicitation (excepté pour Favier). Par contre, le module sécant estimée par l'équation 192 (Ks_3 = K/V) donne un résultat différent.

j'ai essayé le type 1 (W = K*[1- (1 + log(V))/V]) => même remarque que pour type 3. J'ai Ks = K*(1.-2.*log(V))/V**3

et idem pour le type 4 (W = K/2 (V-1)**2 ) => j'ai Ks = K

Et donc pour finir, si je ne me trompe pas :
- dans le cas où la loi fait la séparation exacte volume/forme, l'équation (14) peut devenir générale et non plus dédiée aux seules sollicitations purement sphériques. En fait, cette équation (14) pourrait s'écrire autrement dans le cas général en ne tenant compte que de la partie sphérique de D. On peut alors dans ce cas se contenter d'utiliser les dérivées du potentiel w ou W en suivant l'équation 18 ou 19 pour obtenir le module sécant "exact" avec Ks = d(p)/d(V).
pour les lois type Favier/Orgéas, c'est une autre affaire.

j'ai mis en pièce jointe un répertoire pour tester. Voir Readme.

a+
Julien

Julien,
merci pour ta relecture, du coup cela valide ce que je pensais.
Concernant le module de compressibilité, c'est un sujet qui m'a donné du fil à retordre et je n'ai pas vraiment finit de conclure. Je vais essayer de faire une synthèse de ce que j'ai bidouillé.
Tout d'abord l'utilisation du module:
1) il est utilisé pour calculer la variation d'épaisseur pour les éléments plaques et coques et en contrainte plane et la variation de la section pour les éléments barres et en contrainte doublement plane.
2) il "devrait" être utilisé pour le calcul du pas de temps critique, mais dans la pratique c'est un pseudo module d'Young qui est implanté... mais je devrais changer en m'inspirant de ce qui est implanté dans Ls-dyna. Là j'ai déjà avancé mais je n'ai pas concrétisé (à faire dans le futur).

Donc pour l'instant seul 1) est opérationnel. Au niveau théorique pour les contraintes planes, cela correspond (doc: theorie_herezh++.pdf) aux 2 premières formules du chap. 16.1 (formules: 479,480), suivant que l'on considère un module sécant ou un module tangent. Ce qui donne les formules 482 ou 487 (suivant le type de module considéré).
Si le module tangent est constant, en intégrant 480 on obtient une formule de variation globale en log: 485. Mais si le module tangent n'est pas constant, la formule 485 définit un 3 module, le module sécant en log !!
Si le module sécant est constant, on a la formule classique : 479, qui est donc différente de 485.
Donc déjà on voit qu'il faut faire attention à la définition de base de ce que l'on veut appeler "compressibilité". En regardant la littérature dans le domaine, je me suis aperçu qu'il y avait effectivement plusieurs définitions.

Au niveau implémentation, ce n'est pas tout à fait cohérent ce que j'ai fait. Là aussi il faudrait que j'améliore.
a) au niveau du calcul de la variation d'épaisseur en CP, j'ai utilisé directement la formule de variation globale en log: 485 . Donc ce sera d'autant plus juste que le module sera constant et faux s'il varie beaucoup... Il faut noter que la variation d'épaisseur est celle relative à la loi, c'est une variation d'épaisseur que je qualifie de "mécanique".
b) par contre, la loi est utilisée avec des éléments plaques ou coques et là, je calcule une variation d'épaisseur que je qualifie de "géométrique", via la formule incrémentale 487 qui donc utilise normalement le module tangent. Il faut noter que j'utilise en fait une moyenne de module tangent pour tous les pti de l'élément pour les quadrangles et les triangles, et pour les SFE, le calcul est effectué pour chaque pti.

C'est la même idée concernant les contraintes doubles, mais là j'utilise directement eps33 et eps22 pour les déformations dites "mécaniques" et la formule incrémentale 493 pour la variation de section.

Donc en résumé, on a l'impression que c'est le module tangent qu'il faudrait toujours calculer.
Et bien c'est le module sécant que je calcule dans les lois !! (en général) mais pas toujours.
Par exemple dans les lois hypo, c'est évidemment le module tangent qui est naturellement calculé.

En fait ce que j'ai prévu, et commencé à implanter, c'est de calculer les deux: tangent et sécant...

à suivre ... (c'est pas fini !)

Julien,
juste concernant ta remarque concernant le développement limité: oui on pourrait mettre la limite, mais je pense que c'est plus satisfaisant d'avoir une meilleur continuité à l'interface du "petit" ...
Ceci étant, si ConstMath::petit est suffisamment petit le résultat doit être quasiment identique.
Pour info voici les constantes qui sont implantées et qui sont utilisées dans tout Herezh++:

Concernant la définition de Kt que j'ai introduite en 480, c'est bien une division par v et non v0 que j'utilise donc c'est exactement celle que tu as trouvé dans wikipedia.
Pour autant le pb n'est pas tant le choix d'une définition, mais plutôt: à quoi cela va servir. Car en fait dans mes développements théoriques, la compressibilité est une grandeur intermédiaire qui sert à plusieurs choses. L'important, je crois c'est que l'utilisation soit cohérente avec la définition... ce qui n'est pas tout à fait le cas dans Herezh pour l'instant.

Suite du #9 en réponse à #8
Je suis d'accord avec ta remarque ( -p = Ks log(V) ) en fait défini un module sécant "pour une variation en log" donc c'est différent du Ks du 479 et c'est également différent du Kt du 480 "sauf si Kt est constant" (cf. 485).
Je vais clarifier la chose dans la doc théorique et je vais l'appeler : -p = K_slog log(V)
NB: on retrouve finalement le même type de différence que sur la notion de déformation.

Si je reviens sur l'utilisation du module de compressibilité, K_slog est uniquement utilisé en CP, pour calculer la variation d'épaisseur dite mécanique.
Je pourrais me débrouiller pour éviter, et par exemple utiliser eps33 en tenant compte des différents types de def utilisés... (à voir)

Du coup, on pourrait peut-être tout ramener au module tangent.
Si je reviens sur le cas des lois hyper-élastique, les formules de départ sont: 18 ou 19 (suivant que l'on calcule le potentiel sur le volume initial ou final).
Une remarque concernant le chap. 2.1 de hyper-elasticite.pdf:
Je me place dans un cas purement sphérique de vitesse de déformation pour exhiber la partie pression qui concerne la partie purement sphérique du tenseur de contrainte relativement à une partie purement sphérique du tenseur D.
Si on considère un comportement isotrope, a priori il y aura toujours ce type de découpage: D sphérique vs pression, et D déviatorique versus contrainte déviatorique.
Mais on peut étendre les formule.

Donc oui, les formules 18 et 19 sont correctes pour un matériau hyper-élastique isotrope.
Je vais généraliser dans un cas quelconque, je vais mettre à jour la doc... à suivre

Pour les potentiels utilisant les invariants de Mooney-Rivlin on a directement 188, et effectivement, si on s'attache uniquement à calculer le module tangent on peut directement dérivée 2 fois la formule 18, et donc pour les potentiels implantés on peut ne retenir que la partie volumique du potentiel qui est clairement différencié. Du coup on obtiendra directement K_t.
Les formules que j'ai calculées sont donc relatives à K_slog.
Je vais compléter avec les formules en K_t

Pour les potentiels en V, Qeps et phase que j'appelle de "Grenoble" car Favier et Orgéas sont de Grenoble (c'est historique), je vais préciser la chose, mais on peut également étendre les formules

à suivre ...

concernant : -P = Kt (log(vol) - log(vol_0) ) = K_t log(vol/vol_0)
=> ok, mais du coup, ça serait pas plutôt Ks ? ça ressemble à la définition de Ks. Un peu la même remarque sur équation (480), il y a tantôt Kt, tantôt Ks.

concernant loi Grenoble :
=> alors là j'ai du mal à comprendre. J'ai repris mon répertoire de test de Kt pièce jointe précédente. J'ai calculé une loi Favier (MAT_2), sa dérivée numérique (script perl) et regardé le résultat avec fichier gnu_petit_w. Je n'ai pas fait d'erreur dans les formules. Et pourtant, pour une sollicitation quelconque, effectivement je suis d'accord avec toi, la formule (18) fonctionne pour calculer la pression (pièce jointe pression.png).
Par contre, j'ai des oscillations en calculant Kt (pièce jointe Kt.png). Le comportement tangent par dérivation du calcul Herezh ne donne pas pareil exactement pareil que la formule théorique : Kt = - d(p)/d(V) avec V=v/v0. Je vois bien que l'échelle du graphe est très zoomé, donc le résultat est grosso modo pas mal. Mais je n'avais pas de telles d'oscillations avec les modèles en invariant de B.

Et malgré tout, ce graphe a l'air d'aller dans ton sens (formule (18) vraie aussi pour Almansi dans le cas général, et donc ses dérivées => Kt). Au final, comme Favier/Orgéas est lié à un tenseur de trace nulle, il ne génère pas de pression hydro. Et pourtant, je sais qu'Almansi ne fait pas la décomposition volume-forme. ça me laisse un peu perplexe. Je ne suis pas encore bien persuadé. L'oscillation de la courbe Kt.png me donne l'impression qu'il y a un problème, même si léger.

oui,
-P = Kt (log(vol) - log(vol_0) ) = K_t log(vol/vol_0)
correspond à Kt si Kt constant,
mais de manière générale c'est effectivement un Ks, que j'ai appelé K_s,log, car il s'agit d'un K sécant pour une variation log.
J'ai mis à jour la doc théorique, voir le chap 16.1

Pour les lois de Grenoble, c'est effectivement assez dur à comprendre, c'est un peu pour cela que je voulais échanger de vive voix,
à suivre,
bonne réunion

je suis dispo si besoin pour appel.

et sinon, pour revenir sur mes graphes précédents :
j'ai modifié mon graphe "pression.png" pour mieux comprendre. ci-joint même graphe mais avec ajout de l'écart relatif en % entre simulation et théorie équation (18). Et donc, là, je comprend mieux. Il n'y a pas égalité entre pression simulée et pression théorique(18). Il y a un écart même si faible. Le fait d'utiliser la déf d'Almansi conduit bien à une influence sur le volume même si on utilise que sa partie déviatorique.
Et on le voit bien sur un essai de cisaillement simple. Si tu prends un cube 1x1x1, que tu lui mets les conditions X Y pour faire du cisaillement dans le plan XY, et que tu bloques une des faces selon Z, ça donne bien :
- si loi en invariant almansi => variation de l'épaisseur du cube le long de Z => volume diminue
- si loi en invariant de B => pas de variation (volume reste à 1)

Et j'ai donc fait l'erreur d'associer volume et trace nulle du tenseur de déf. La partie déviatorique Almansi ne génère pas directement de termes sur la trace de sigma mais elle influencera le volume donc la pression.

ça veut dire la formule (18), dans le cas général pour Almansi, devrait introduire des termes liés à la dérivation par rapport à Q_eps.
Parce que la comparaison à la simulation montre que :
1) pour une valeur de V donnée, il y a un écart entre théorique(18) et simulation. Donc il y a quelque chose qui manque dans la dérivée première du potentiel selon V.
2) pas exact non plus Kt. Donc => manque aussi quelque chose dans les dérivées secondes

modif théorie,
- précisions apportées aux modules de compressibilités : K_s devient K_{s1} et K_{s2} avec une remarque lorsque K_{s2} est constant (chap. 16.1)
- pour l'élasticité de Hooke: calcul de l'approximation du module tangent dans le cas de faible déformation d'Almansi et dans le cas de def log