Anomalie #397
ouvertcontact : vecteur normale de norme 2 au lieu de 1
Description
Gérard,
l'exemple ci-joint est le fameux montage "joint axisymétrique écrasé entre chemise et fouloir". Le calcul converge jusqu'au bout. Mais j'ai repéré un potentiel souci au temps 0.418 :
La force de contact devient soudainement élevée à l'un des noeuds. Il se trouve que la normale de contact à cet instant et à ce noeud à une norme de 2 au lieu de 1.
On dirait que ça intervient juste quand le noeud est en train de passer d'un élément de contact à un autre. Donc au moment où il passe par une frontière "noeud" entre 2 quadrangles. Est-ce qu'il n'y aurait pas un pb de normalisation de la normale pour cette situation un peu particulière ?
FORCE_CONTACT :
NORMALE_CONTACT :
Fichiers
Mis à jour par Julien Troufflard il y a environ un mois
- Fichier temps_0p418_FORCE_CONTACT.png temps_0p418_FORCE_CONTACT.png ajouté
- Fichier temps_0p418_NORMALE_CONTACT.png temps_0p418_NORMALE_CONTACT.png ajouté
Les figures ne sont pas passées dans le message initial
FORCE_CONTACT :
NORMALE_CONTACT :
Mis à jour par Gérard Rio il y a 7 jours
- Statut changé de Nouveau à En cours
- % réalisé changé de 0 à 10
La raison est:
- la sortie de la normale est de norme 1 si le noeud fait partie d'un seul élément de contact et de n si le noeud fait partie de n éléments de contact.
Même moi je ne m'en rappelais plus, mais c'est marqué dans la doc.
En fait j'avais mis cela pour justement voir si le noeud était en contact multiple.
Bref ici il s'agit d'un noeud qui appartient à 2 éléments de contact. Mais ici ce n'est pas normal car il se trouve que les deux éléments de contact sont identiques ! et normalement il y a un filtre qui interdit cela ... il semble qu'il y ait quelque chose de pas normale !
Il faut que je trouve le pourquoi ...
à suivre
Mis à jour par Julien Troufflard il y a 6 jours
Le fait de compter le nombre de contacts dans la normale ne serait-il pas redondant avec NOEUD_PROJECTILE_EN_CONTACT ?
En tous les cas, il y a une force de contact très élevée à ce noeud. Exactement le double du noeud en dessous. On dirait que la force de contact est calculée pour une normale de norme 2. Donc ce ne serait pas uniquement un indicateur en sortie ?
Mis à jour par Gérard Rio il y a 5 jours
- % réalisé changé de 10 à 30
oui, tout à fait, mais c'était fait un peu à dessein.
Donc effectivement deux fois le même élément est comptabilisé ce qui double la force de contact.
Une des raisons de l'existence de ces 2 contacts simultanés est que lorsque l'on glisse d'une frontière à l'autre on peut avoir un contact sur 2 frontières, car le calcul prend en compte un léger débordement de la frontière pour le calcul de la projection et actuellement j'avais permis d'avoir un noeud en contact avec plusieurs frontières en même temps.
Du coup après réflexion voici ce que je propose:
- un noeud ne peut avoir qu'un seul contact par zone de contact.
Je pense que le fait de pouvoir définir plusieurs zones de contact permet quand d'avoir plusieurs contacts, mais ceux-ci sont alors bien identifiés. Par exemple on peut avoir sur un même noeud esclave un contact vertical + un contact horizontal
Et donc lorsque je détecte 2 contacts sur une même zone, je choisis celui dont la distance du noeud à tdt à sa projection est la plus faible.
Mis à jour par Julien Troufflard il y a 5 jours
je ne sais pas si c'est la bonne approche de ne garder qu'un seul contact par zone de contact. En tant qu'utilisateur, je ne pense pas que ce soit très pratique de devoir penser à séparer en deux les zones de contact pour penser à une partie verticale et une autre horizontale (tellement de gens ne vont pas y penser, y compris moi).
Pourquoi ne pas simplement diviser par la force de contact par le nombre de contact ? Cela permettrait de :
1) retrouver une force de contact correcte
2) continuer avec le fonctionnement actuel qui permet de converger plein de calcul. J'ai peur que ne garder qu'un seul contact conduise à tomber sur des cas particuliers avec risque de pénétrer dans la matière (dans un coin typiquement).
autre argument : imaginons une surface de contact très accidentée. Il serait très fastidieux de séparer cette zone de contact en plusieurs pour permettre à un noeud de tenir compte de 2 contacts avec cette surface.
ceci dit, je ne suis pas sûr de bien cerner le fonctionnement dans Herezh++ pour calculer la normale. J'imagine un noeud en contact avec un nombre quelconque N d'éléments. J'aurais tendance à imaginer que la normale devrait être égale à la somme des normales à ces surfaces, puis une normalisation de ce vecteur. Et ensuite, affecter un coef de pénalisation à cette normale résultante unitaire.
Mis à jour par Gérard Rio il y a 5 jours
je ne sais pas si c'est la bonne approche de ne garder qu'un seul contact par zone de contact. En tant qu'utilisateur, je ne pense pas que ce soit très pratique de devoir penser à séparer en deux les zones de contact pour penser à une partie verticale et une autre horizontale (tellement de gens ne vont pas y penser, y compris moi).
Je me suis mal expliqué: pour tout contact, l'effort est toujours suivant la normale donc il n'y a pas à penser décomposer entre verticale et horizontale.
Je voulais seulement évoquer le fait que l'on pouvait forcer un contact multiple, si par exemple on décrit les maîtres avec plusieurs figures géométriques superposées. Dans ce cas on peut avoir plusieurs contacts sur un même noeud. Mais c'est un cas particulier.
Pourquoi ne pas simplement diviser par la force de contact par le nombre de contacts ?
oui, a priori c'est possible et c'est une solution que j'ai envisagée. Cependant elle pose plusieurs pb dont un que je ne sais pas résoudre simplement (mais qui serait sans doute possible) :
1) il faut qu'au moment du calcul des résidus et raideurs de contact je prenne en compte le nombre de contacts. Dans Herezh, chaque élément de contact calcul ces grandeurs de manière autonome: ce sont des grandeurs locales indépendantes des autres éléments de contact. Une technique serait d'introduire une couche intermédiaire entre éléments de contact et le calcul global du résidu + raideur. Mais il y aurait d'autres méthodes possibles. En tout cas ce n'est pas un gros problème.
2) Dans le cas d'un glissement entre deux frontières, qui est en fait le cas qui introduit les pb de double (ou plus) de contact, ce pose le pb de l'historique du frottement. Soit au temps t un seul contact et au temps t+dt, 2 (ou plus ) contacts. Le (ou les) nouveau(x) contact n'ont pas d'historique du coup si le premier contact disparait il faut transmettre l'historique pour que le calcul de frottement soit correct, mais si j'ai plusieurs nouveaux contacts, à qui je transmets ?
Fondamentalement c'est le même noeud qui glisse donc il n'y a pas de raison d'avoir plusieurs contacts actifs pendant le glissement. L'apparition de plusieurs contacts est due à la méthode de calcul qui intègre des approximations.
Si je garde un seul contact, je n'ai pas de pb pour conserver l'historique.
Ceci étant, dans Herezh l'apparition momentanée de plusieurs détections de contact est réelle et il faut pouvoir la gérer.
Dans des premières versions, je partais du principe que si un contact actif d'un noeud esclave existait je ne cherchais pas d'autre contact possible. Bon, là je passais à côté de cas où les nouveaux contacts possibles étaient également acceptables.
Donc par la suite j'ai introduit la possibilité d'avoir plusieurs contacts avec les pb actuels
La méthode que j'ai maintenant implantée est de considérer les nouveaux contacts possibles, de choisir celui qui est le plus pertinent, et de transférer sa situation de contact sur le noeud historique qui lui continue de garder l'ensemble de l'historique.
Théoriquement j'ai l'impression que cela devrait fonctionner. Mais en fait il faut que je test si d'une manière pratique, cela fonctionne et ... là je vais peut-être déchanter et devoir revoir ma copie.
Je vais toujours tester ce que j'ai mis en place
À suivre !
Mis à jour par Gérard Rio il y a 4 jours
- Fichier deformee.png deformee.png ajouté
- % réalisé changé de 30 à 90
a priori cela semble assez bien fonctionner sur l'exemple de l'écrasement du joint.
- il n'y a plus de pb de normale
- la convergence est ok
Par contre au cours de mes tests j'ai relevé les points suivants:
1) comme la loi proposée est hyperélastique avec un module de compressibilité très grand par rapport au module de cisaillement, à partir d'un certain taux de compression, la convergence globale devient difficile. Cela n'a, apriori rien à voir avec le contact, c'est tout simplement que l'on observe un effet de blocage numérique venant du fait que la partie sphérique écrase complètement la partie déviatorique au niveau de la raideur, et cette dernière manque pour l'unicité de la solution.
Je fais le même essai avec un K 10 fois plus petit et là aucun pb pour écraser fortement: hauteur initiale 7.8 , déplacement imposé -5.6
donc il faut se méfier des non-convergences supposées dues au contact. Malgré le fait que le contact peut entraîner pas mal de pb sur la convergence.
2) toujours concernant la convergence due au contact, j'ai noté que lorsque l'on atteignait une forte compression, la géométrie finale avait tendance à osciller, ceci due à ce que j'ai expliqué au 1). Du coup certains noeuds en contact avaient tendance momentanément à sortir du maître. Or dans la mise en données on indique:
PENETRATION_BORNE_REGULARISATION 1e-5
ce qui veut dire que si le noeud sort de plus de 1.e-5, on ne tient plus compte de sa réaction. Et du coup sur la visualisation on voit le noeud qui s'écarte violemment du maître avec lequel il n'est donc plus accroché.
Pour éviter ce phénomène, il faut introduire un autre paramètre:
PENETRATION_BORNE_REGULARISATION_PLUS 0.1
Dans ce cas tant que le noeud ne s'écarte pas plus de 0.1, le noeud reste accroché.
Donc en utilisant ce paramètre, cela permet d'éviter des décollements intempestifs lorsque globalement la solution n'est pas parfaitement stable.
Néanmoins pour la loi considérée dans l'exemple, il arrive un moment où le calcul patine due au 1) et là il n'y a pas de solution en implicite. En tout cas je n'en ai pas trouvé. Peut-être en RD ?
La nouvelle version est la : V 7.061
NB: On peut également y tester les sorties individuelles pour chaque maillage
Mis à jour par Julien Troufflard il y a 3 jours
- Fichier calcul_Gmsh_v7.056.tar calcul_Gmsh_v7.056.tar ajouté
- Fichier calcul_Gmsh_v7.061.tar calcul_Gmsh_v7.061.tar ajouté
- Fichier test_v7.056.log test_v7.056.log ajouté
- Fichier test_v7.061.log test_v7.061.log ajouté
- Fichier test.tar test.tar ajouté
- Fichier cas_test.png cas_test.png ajouté
ok pour les pb de compressibilité.
J'ai construit un cas test en AXI assez intéressant. Au départ, c'était pour valider cette nouvelle version Herezh. Mais finalement, il montre des pb pour les 2 versions : celle que j'utilise (7.056) et la dernière 7.061.
ci-joint les fichiers de calcul dans test.tar, le résultat Gmsh par version et le log des 2 calcul.
Pour ce calcul, j'ai pris une loi Mooney-Rivlin assez compressible, équivalente en petite déf à une loi élastique log E=100, NU=0.35
Le principe est un carré (maillage plot_AXI.her) écrasé dans un logement en U (logement_AXI.her). La situation est la suivante :
Comme indiqué, le but est d'avoir une config qui teste :
- le cas d'un noeud esclave qui franchit un noeud maitre (ici 2 cas : un cas horizontal et cas vertical).
- un noeud dans le coin qui va subir un contact vertical + horizontal
- il y a un petit jeu initial selon Y de 1e-4 (pas de contact initial) et un gros jeu initial selon X pour décomposer le contact vertical/horizontal au niveau du coin
Autre point : le calcul utilise un suite_point_info pour décomposer en 2 phases
1) Il ne se passe rien entre t=0 et t=1 => le maillage du plot ne doit subir aucune déformation. Il ne doit y avoir aucune force de contact (c'est en référence au ticket https://herezh.irdl.fr/issues/391 car je testerai ensuite en 2D et ainsi qu'en 3D en faisant une extrusion Z de ces maillages)
2) L'écrasement a lieu entre t=1 et t=2 par déplacement imposé selon Y sur le bord du haut du plot => le contact vertical doit être détecté dès le premier pas de temps, puis interviendra le contact horizontal
Pour des raisons de nombre de fichiers pièce jointe, je vais montrer quelques résultats dans les messages suivants.
Mis à jour par Julien Troufflard il y a 3 jours
- Fichier v7.056_exemple1.png v7.056_exemple1.png ajouté
- Fichier v7.056_exemple2.png v7.056_exemple2.png ajouté
- Fichier v7.056_exemple3.png v7.056_exemple3.png ajouté
- Fichier v7.056_exemple4.png v7.056_exemple4.png ajouté
- Fichier v7.056_final.png v7.056_final.png ajouté
résultats pour l'ancienne version 7.056
pas pb pour phase 1) du calcul entre t=0 et t=1 => pas de déformation, pas de force de contact (le contact est quand même détecté car il y a une normale de contact en certains noeuds). Donc ok pour cette partie.
Ensuite, dans la phase d'écrasement, voici quelques exemples :
au temps 1.14 : le noeud du coin se prépare au contact horizontal. Il a bien une normale à 45°.
Problème : il a par contre une norme de normale égale à sqrt(2) au lieu de 1
au temps 1.21 : le plot est au contact des 2 parois avec des forces de contact non nulles selon X et Y
Problème : la norme de contact au coin est maintenant inférieure à 1
au temps 1.23 : il y a eu franchissement d'un noeud esclave en bas.
Problème : on retrouve le problème initial de ce ticket. Ce noeud a désormais une norme de contact de norme 2
au temps 1.9 : malgré les défauts de norme, le contact se passe correctement avant ce temps.
Problème : à ce moment du calcul, le coin perd le contact et pénètre dans le maitre. La normale de contact n'a plus ni la bonne norme, ni la bonne direction
au temps final : le calcul converge malgré la perte totale de contact du noeud coin.
Mis à jour par Julien Troufflard il y a 3 jours
- Fichier v7.061_exemple1.png v7.061_exemple1.png ajouté
- Fichier v7.061_exemple2.png v7.061_exemple2.png ajouté
- Fichier v7.061_final.png v7.061_final.png ajouté
ps : oubli de préciser pour la version 7.056 => il n'y a pas de pb de norme de normale pour le noeud esclave en haut à droite qui franchit un noeud maitre sur la paroi verticale. Elle reste à 1 contrairement au cas du noeud du noeud en bas à gauche sur la paroi horizontale. bizarre ??
résultats pour la nouvelle version 7.061
ok pour la phase 1) => pas de déformation ni forces de contact
pour la partie écrasement : en général, il n'y a plus de pb de norme de normale pour tous les noeuds. Le souci est désormais uniquement sur le noeud coin qui oscille entre une normale purement horizontale ou purement verticale avec une impossibilité de trouver le contact correctement.
Voici 2 exemples parmi tant d'autres :
au temps 1.21 : normale purement verticale et noeud a pénétré dans le maitre
au temps 1.26 : idem avec un normale purement horizontale
au temps final : on aboutit à peu près à la même situation qu'avec version 7.056 => calcul converge malgré un coin qui a perdu totalement le contact.
Mis à jour par Gérard Rio il y a 3 jours
Le cas que tu montres Julien, est celui que j'évoquais en parlant de décomposition de la zone de contact.
Je crois (?) de mémoire qu'il a déjà été évoqué dans un autre ticket dont je ne me rappelle plus le num.
bon, ceci étant je vois deux réponses:
1) pour un angle droit ou un angle aigu provenant d'une réelle géométrie du maître, je pense que la solution est de définir 2 zones de contact, une pour chaque face du contact. Ce qui permettra d'avoir 2 forces de contact, qui auront en plus leur propre historique. Par exemple dans le cas du joint compressé il y a en réalité 2 zones de contact.
2) Par contre s'il s'agit d'une modélisation d'une pièce réelle, par exemple un fond de moule. Une solution est de mailler plus finement le rayon de fond d'angle ce qui permettra d'avoir une évolution progressive du contact. Dans la réalité, il y a toujours un rayon de raccordement dû à la fabrication et le contact va dépendre logiquement de ce rayon.
Mis à jour par Julien Troufflard il y a 3 jours
- Fichier v7.056_1.png v7.056_1.png ajouté
- Fichier v7.056_2.png v7.056_2.png ajouté
- Fichier v7.061_1.png v7.061_1.png ajouté
- Fichier v7.061_2.png v7.061_2.png ajouté
oui tu as sûrement raison. Je crois même que ça fait parti des bonnes pratiques EF d'ajouter un petit rayon aux angles vifs, même fictifs (i.e pas issu d'une pièce réelle).
Donc, il faudra bel et bien penser à décomposer les zones de contact. J'y étais réticent, mais finalement, ça parait d'une part logique et d'autre part largement faisable dans la grande majorité des cas que l'on a traiter.
Si je décompose le logement en 2 zones de contact, ça conduit les 2 versions Herezh au bon résultat final (pas de pénétration).
avec néanmoins une normale étrange dans le cas de l'ancienne version 7.056
version 7.056 :
1) avant que le contact ne devienne effectif sur la partie verticale, la normale au coin est dans la bonne direction
2) ensuite, la normale n'est plus dans le bon sens. L'exemple ci-dessous est représentatif du reste du calcul avec une normale au coin rentrante dans la matière et une normale de norme 2 au noeud en bas à gauche du plot. Mais néanmoins, le calcul converge avec un contact correct. Bizarre pour que ça fonctionne vue la normale inversée??? (Mais bon inutile d'investiguer puisque ce ne sera pas le cas avec la version 7.061)
version 7.061 :
1) l'exemple ci-dessous est représentatif de la quasi intégralité des incréments du calcul => normale au coin dans la bonne direction
2) dans les derniers incréments, la normale au coin devient purement verticale. Mais ça n'empêche pas le caclul de converger et avec un contact correct sans pénétration.
Au final, on a une norme au coin qui peut être différente de 1 (somme des normales de contact). Je ne sais pas si c'est un réel problème. A priori non (chaque contact est calculé avec une norme 1 multiplié par un coef de pénalisation => la somme des 2 donne une résultante pour l'affichage Gmsh).
Par curiosité, je vais tenter d'écraser un peu plus (pour voir si la normale finale pose problème, car c'est quand même étrange de ne plus avoir une résultante horizontale+verticale).
Puis ensuite je vais faire l'ajout d'un rayon de courbure au coin et voir ce que ça donne avec une zone de contact unique. Et mettre le résultat ici.
Ensuite je vais donc passer au cas 2D, et donc sans doute basculer sur le ticket https://herezh.irdl.fr/issues/391
Mis à jour par Julien Troufflard il y a 3 jours
si tu passes le déplacement d'écrasement à 0.75 :
nom_mail= plot_AXI N_N 'UY= COURBE_CHARGE: courbe_dpiY ECHELLE: -0.75'
tu devrais obtenir un pb que je vois assez souvent en contact (segmentation fault) :
max des reactions = 409670.36243119 max du residu total = 409015.42953461 **** Norme Residu/Reaction --> 0.99840132 (assemblage 1) delta maxi: 0.16731475, ddl X1, noeud/mail 2/1 , var total:-0.89937342, sur:1.07062658 autres ddl R_X1 0.00000000 var 50.52351377; ============================================================================ ******* NON convergence des iterations d'equilibre ********* Segmentation fault
D'où vient ce pb qui est sans doute un pb de mémoire. Herezh essayerait-il d'utiliser une variable qui n'existe plus ?
Mis à jour par Julien Troufflard il y a 3 jours
- Fichier calcul.info calcul.info ajouté
- Fichier logement_AXI.her logement_AXI.her ajouté
- Fichier logement_AXI.lis logement_AXI.lis ajouté
j'ai oublié de te donner les nouvelles versions de maillage logement_AXI.her et fichier .info (avec la déclaration de 2 zones de contact et le déplacement -0.75)
Mis à jour par Julien Troufflard il y a 3 jours
je voudrais aussi signaler la sensibilité à ce paramètre :
para_calculs_geometriques POINT_INTERNE_PREC_THETAI_INTERNE 1.e-4
avec ce paramètre, on se retrouve avec une normale au coin purement verticale en fin de calcul (puis segmentation fault)
Si je le passe à 1-6 (la valeur par défaut) => le contact au coin n'est jamais correct (normale qui oscille entre sa partie horizontale ou verticale pure) et systématiquement une pénétration du noeud dans le maitre.
Si je le passe à 1e-3 (ou même 1e-2), la normale est correcte avec une composante horizontale+verticale, même en fin de calcul. Donc la bonne situation recherchée pour ce contact. Mais néanmoins toujours ce segmentation fault (le calcul n'atteint pas la fin).
D'où vient cette dépendance à ce paramètre ? est-ce une tolérance relative ou absolue ? Quelle est la bonne valeur et est-ce risqué de le mettre à 1e-3 ?