Anomalie #262
résultat erroné traction barre 1D à grand nombre d'éléments
Ajouté par Julien Troufflard il y a environ 4 ans. Mis à jour il y a environ 4 ans.
Description
Gérard,
je suis tombé sur un problème étrange en essayant de comparer le cas "fonction tabulée" par rapport au cas "fonction expression littérale". J'ai voulu mettre en place un exemple de barre 1D avec beaucoup d'éléments pour comparer les temps de calcul.
Le test est une barre d'éléments bielette en traction uniaxiale avec une loi ISOELAS. Je charge la barre par déplacement imposé. Le problème est donc très basique. Mais bizarrement la réussite du calcul dépend du nombre d'éléments. Dans l'exemple ci-joint, le résultat devient erroné (non homogène) à partir du maillage à 2000 éléments. Avec le maillage à 5000 éléments, on voit même qu'une grande partie de la part ne se déplace même pas.
Je soumets ce pb car je n'en vois pas l'explication (le cas 5000 éléments tourne sans pb sur d'autres codes par exemple) et en espérant contribuer à résoudre un pb qui pourrait être silencieux dans d'autres cas d'application avec des éléments 1D.
Fichiers
| pb_barre1D_nb_elts.tar (166 ko) pb_barre1D_nb_elts.tar | Julien Troufflard, 29/09/2020 09:40 | ||
| pb_barre1D_section_infini.tar (101 ko) pb_barre1D_section_infini.tar | Julien Troufflard, 29/09/2020 15:10 |
Mis à jour par Gérard Rio il y a environ 4 ans
- Statut changé de Nouveau à En cours
- % réalisé changé de 0 à 90
oui c'est normal avec la mise en données que tu est indiquée. La raison est la suivante:
1) plus le maillage est dense plus la taille de chaque élément devient petite
2) lors du premier incrément, la condition cinématique est appliquée telle quelle.
Du coup, par exemple si on a 5000 éléments pour une longueur de 100 -> un élément a une longueur de 100/5000 = 0.02
Dans ton exemple tu imposes un déplacement de 50mm avec un pas de 0.1 donc au premier incrément tu imposes un déplacement de 5mm.
L'élément qui contient le noeud imposé voit donc sa longueur augmentée de 5/0.02 = 250 fois. Avec une loi isoelas cela ne converge pas tel quel ... ça me paraît normal.
Donc en diminuant le deltat de départ par exemple en mettant
DELTAtMAXI 0.1
DELTAt 0.001
avec un maillage de 2000 et de 5000 éléments, il y a quelques divergences dans les premiers incréments puis c'est ok jusqu'à la fin, l'incrément de temps augmente progressivement jusqu'à 0.1 pour les derniers incréments.
Par contre le type de comportement du calcul est sans doute dépendant de la loi de comportement. Avec une loi hypo ou hyper on peut sans doute utiliser d'autres paramètres...
Mis à jour par Julien Troufflard il y a environ 4 ans
oui effectivement, il y a une dépendance au pas de temps. J'avais joué un peu avec. J'avais testé 1e-3 comme toi, mais sur le cas à 5000 éléments, il subsistait une erreur (résultat pas homogène d'après l'échelle Gmsh sur par exemple : SIG11, SECTION_MOY_FINALE). J'avais donc abandonné cette piste. Mais en fait il faut descendre à 1e-4 sur le cas 5000 éléments pour que cela donne réellement le bon résultat et homogène.
Je comprends ton explication, cela parait logique. En regardant de plus près le cas UX=0.5, DELTAt 0.1, maillage 5000 éléments, j'ai observé que l'avant dernier noeud (le n°5000 qui celui juste avant celui où le déplacement est imposé qui est le 5001) est projeté violemment plus loin que le noeud 5001. L'élément 5000 se retrouve en compression et pousse sur une extrémité de l'élément 4999 qui est lui en traction (comme le reste de la barre), ces 2 éléments ayant chacun une contrainte à peu près égale mais de signe opposé. Au final, cela se transforme en calcul de traction classique mais c'est l'élément 5000 qui transmet l'effort de traction au reste de la barre car il est en compression vu qu'il est retourné et que l'on pousse selon UX sur l'un de ses noeuds (le 5001). Tout ce bazar a l'air d'être en équilibre, étonnant comme solution.
Bon je trouve néanmoins ce comportement de convergence étrange par comparaison avec un code commercial bien connu ( ;) ) pour lequel quelque soit le pas de temps initial (même 1) le résultat final ne change pas d'un iota. De plus, quelque soit le pas de temps, le nombre d'itération est égal à 1. Autrement dit, le cas barre 1D en traction avec loi élastique devrait a priori converger direct.
Je n'ai pas d'idée sur la différence entre ces 2 codes : résolution de l'équilibre ? méthode d'application de la condition limite ?
Je ne sais pas s'il faut creuser plus que ça...
Mis à jour par Gérard Rio il y a environ 4 ans
- Statut changé de En cours à Résolu
- % réalisé changé de 90 à 100
Il faut bien noter que la loi de Hooke pour des déformations supérieures à 10% par exemple, ce n'est pas "du tout" adapté : il faut utiliser une loi hyperélastique.
D'un point de vue technique, plus l'élongation est grande moins la contrainte de Cauchy via la loi de Hooke et la mesure d'Almansi devient sensible. Du coup la convergence se dégrade complètement.
Si l'on utilise une loi de Mooney Rivlin 1D
elastomere MOONEY_RIVLIN_1D
# ----------------------------------------------------------------------
# |... loi de comportement hyper elastique 1D de type mooney rivlin ...|
# | .. deux coefficients C01 et C10 .. |
# ----------------------------------------------------------------------
#
C01= 0.0167 C10= 0.2
La convergence s'effectue avec 5000 éléments, 1 pas de temps, 4 itérations !
Maintenant si au lieu d'allonger la barre on la comprime (je change le sens de la sollicitation), avec une loi de Hooke 1D -> convergence avec un deltat de 1 et 6 itérations !
Cela ne veut pas dire que la loi de hooke est bonne en compression et pas bonne en traction dans les grandes transformations, c'est simplement due au fait que la contrainte de Cauchy calculée via la def d'Almansi est sensible en compression même pour des déformations très grandes, car la def d'Almansi varie de -l'infinie à +0.5
Si on utilise une def logarithmique, on a une sensibilité dans les deux sens. Mais encore une fois, même avec une def log, la loi de Hooke en grande élongation, ce n'est pas à utiliser.
Si dans la mise en données d'Herezh, j'indique une déformation logarithmique:
materiaux
MAT ISOELAS1D
250. 0.45
type_de_deformation DEFORMATION_LOGARITHMIQUE
j'ai une convergence en traction en 4 itérations en compression et en traction !
Mis à jour par Julien Troufflard il y a environ 4 ans
- Fichier pb_barre1D_section_infini.tar pb_barre1D_section_infini.tar ajouté
ah oui, ok pour la mesure de déf. L'autre code est effectivement en mesure log.
J'ai poursuivi mon but de comparer "loi tabulée" avec "loi expression littérale", et je retombe sur un cas bizarre malgré un pas de temps petit (même très petit). Toujours ma barre en traction, maillage 5000 éléments. Je mets une loi ISO_ELAS_ESPO1D avec une fonction bidon juste pour tester. Je voudrais comprendre pourquoi sur le cas en pièce jointe, l'intégralité de ma barre à une section "nan" à l'exception de l'élément en bout de barre ??? et ceci dès le premier incrément. Le calcul quant à lui converge sagement sans crier. Les 4999 éléments de section "nan" ne sont pas du tout distordus puisque leurs noeuds n'ont pas bougé.
rq 1 : le fichier gnu trace en fonction du temps la valeur de la section de l'élément 1, 2500 et 5000, soit début, milieu et fin de barre
rq 2 : j'obtiens le même résultat que ce soit avec la loi tabulée f_epsAl_tab ou la loi littérale f_epsAl_fct_deriv
Mis à jour par Gérard Rio il y a environ 4 ans
oui c'est également normal ...
La section est calculée (cf. théorie) à partir du module de compression (E/(3.*(1.-2.*nu)))qui lui-même est calculé avec le E = E_init*f(x) . Or celui-ci est nulle à l'origine car la fonction que tu données donne x c-a-d directement |eps_1^1| c-a-d 0 au début.
Si tu modifies la fonction en:
f_epsAl_fct_deriv COURBE_EXPRESSION_LITTERALE_AVEC_DERIVEE_1D
f(x)= (x+1.e-12)
f'(x)= 1.
f"(x)= 0.
fin_parametres_courbe_expression_litterale_
là il n'y a plus de nan
Par contre l'évolution du module d'Young résultant n'est pas (a priori ?) cohérente avec nu=0.3 du coup on obtient une augmentation de la section au lieu d'une diminution. Là aussi c'est sans due au fait que le nu utilisé dans les formules de calcul de K est un nu cohérent avec les petites def...
Si le phénomène persiste et que tu veux absolument avoir l'évolution du E= E_init*f(x) telle que tu l'écris, il faudra peut-être utiliser une loi plus adaptée aux grandes transformations...
Mais a priori (?) il n'y a pas de bug particulier
Mis à jour par Gérard Rio il y a environ 4 ans
par exemple si on ajoute uniquement "1" à la fonction f(x)
f_epsAl_fct_deriv COURBE_EXPRESSION_LITTERALE_AVEC_DERIVEE_1D
f(x)= (x+1.e-12)+1
f'(x)= 1.
f"(x)= 0.
fin_parametres_courbe_expression_litterale_là on obtient bien une section qui diminue...
Mis à jour par Julien Troufflard il y a environ 4 ans
ah oui, effectivement, j'ai choisi un très mauvais exemple pour tester la mise en données.
Bon, je vais pouvoir passer à la vraie fonction que je voulais tester.