Anomalie #326
dilatation thermique LOI_ISO_THERMO et type de loi méca
Description
Bonjour Gérard,
je suis en train de regarder le fonctionnement de la loi thermique LOI_ISO_THERMO. En gros le contexte : j'ai voulu regarder la sensibilité d'une sphère sous pression à la dilatation thermique. Je me suis rendu compte que le résultat final dépend du type de loi méca que l'on choisit.
je suis donc reparti de la base sur un cube en 3D pour comprendre. Avec bien sûr, un fort sentiment que c'était une histoire de type de mesure de déformation, la loi thermo étant définie en mesure Almansi, tandis que la loi méca ça dépend selon qu'on prenne ISOELAS, HYPOELAS, ISOHYPERBULK, etc...
Le cas test est un cube en dilatation libre :
- pas de chargement méca
- des conditions de symétrie
- une grosse variation de température de +100
- un fort coefficient de dilatation pour amplifier les différences entre loi => alphaT = 1000e-6 (si on met vraiment trop, ça ne converge pas)
remarque : j'ai mis une grosse précision (PRECISION 1.e-9) pour être bien sûr d'avoir établi un équilibre correct
dans le fichier .maple, j'ai regardé les grandeurs :
- Almansi_totale11 : déformation almansi totale 11 (colonne 16)
- EPS11 : déformation almansi 11 => partie mécanique (colonne 13)
- SIG11 : contrainte 11 (colonne 22)
J'ai testé divers loi méca. Pour toutes les lois, j'ai mis ce qui correspondrait à peu près à un module de compressibilité de 250 MPa. Mais, au final, les résultats ne dépendent pas de cette valeur (heureusement!), et ceci pour aucune des lois testées.
Dans ce cas test, on s'attend à avoir une partie méca EPS11 nulle et une contrainte nulle. Et donc une déformation totale égale à la déformation thermique (1000e-6*100 => eps_thermique = 0.1)
Concernant la loi méca, voici ce que j'ai testé et ce que j'obtiens :
- ISOELAS
=> pas de problème. J'obtiens bien : Almansi_totale11 = def thermique (EPS11=0, SIG11=0)
(remarque : une certaine difficulté à converger avec alphaT=1000e-6. Pas de pb de convergence pour des alphaT plus classiques autour de 100e-6)
- ISOHYPERBULK3
=> la partie méca est non nulle. J'obtiens :
Almansi_totale11 = 8.023350651496e-02 , EPS11 = -1.976649348504e-02
on a bien : Almansi_totale11 - EPS11 = 0.1 (soit la def thermique prévue)
malgré EPS11 non nulle, il y a une contrainte nulle : SIG11 = 3.416070845000e-13
- HYPO_ELAS3D
=> la partie méca est non nulle. J'obtiens :
Almansi_totale11 = 9.146359655596e-02 , EPS11 = -8.536403444036e-03
comme précédemment : Almansi_totale11 - EPS11 = 0.1 (soit la def thermique prévue)
comme précédemment, il y a une contrainte nulle : SIG11 = -7.899711840594e-11
- MOONEY_RIVLIN_3D avec potentiel vol=4 : K/2(V-1)**2
=> le calcul tourne mais le cube ne se déforme pas du tout
Almansi_totale11 = 0
EPS11 = -0.1
Je comprends bien le résultat pour ISOELAS. Je comprends également qu'il doit y avoir un bug quelque part pour MOONEY_RIVLIN_3D.
Mais pour ISOHYPERBULK3 et HYPO_ELAS3D, pas évident de comprendre. Le fonctionnement n'est pas conforme à ce qu'il faudrait théoriquement obtenir (pas de déformation mécanique).
Pour ISOHYPERBULK3, j'ai essayé d'en prédire le résultat depuis l'équation théorique hyperélastique -p = = 0 = V(dw/dV) + w.
Et idem avec HYPO_ELAS3D en partant de trace(sigma) = 0 = K_c.trace(eps_tot - eps_thermique).
Mais au final, je ne trouve jamais ce que trouve Herezh.
Quel calcul est fait par Herezh sur ce cas test de cube ?
N'y a-t-il pas une modif à faire quelque part sur le fonctionnement de la dilatation thermique pour éviter une déformation méca non physique ? (rendre tout ça indépendant de la loi méca choisie)
question supplémentaire : dans le cas de ISOHYPERBULK3, la variable associée est V=v/v0. Quelle est la mesure de déformation que l'on pourrait associer à cette loi ? Peut-on imaginer qu'il y ait une partie de V "mécanique" et une partie "thermique" et donc un V "total" ????
a+
Julien
Fichiers
Mis à jour par Gérard Rio il y a presque 2 ans
- Statut changé de Nouveau à En cours
jolie travail Julien !!
bon, je vais creuser et te tiens au courant.
Mis à jour par Gérard Rio il y a presque 2 ans
- % réalisé changé de 0 à 20
je me suis tout d'abord intéressé au cas de la loi hypo.
- A priori, le fonctionnement d'Herezh est conforme à ce que j'avais pensé, je ne détecte pas d'erreur informatique.
Du coup cela relève des choix théoriques.
1) La dilatation thermique est implantée selon: eps_thermique = alpha * (T - T0) et delta_eps_thermique = alpha * delta T
la déformation mécanique vaut: eps_totale - eps_thermique = eps_meca et delta_eps_meca = delta_eps_totale - delta_eps_thermique
avec le choix: eps_totale est une def d'Almansi, idem pour le delta_eps_totale
2) dans le cas de la loi hypo: la contrainte est calculée à partir du delta_eps_meca (ou de D_eps_meca) ce qui fait qu'à chaque incrément la condition sigma 0 , qui correspond dans notre cas à delta sigma 0 (car la contrainte du pas précédent est nulle) conduit à un delta_eps_meca nulle. Donc au finale on a la somme des delta_eps_meca nulle, ce qui dans un cas sans thermique correspond à une def log approchée nulle.
Compte tenu de ces hypothèses et des conditions du test on a à la fin du calcul:
. eps_thermique = 0.1 car c'est une valeur imposée
. une déformation totale qui par défaut est calculée en Almansi:
0.5*(1. - 1./(1.1**2)) = 0.086777
. la déformation méca (1,1) calculée par Herezh: -8.536403444036e-03
Si je calcule la def totale en log cela donne : log(1.1) = 0.095310
d'où une différence entre les deux def: 0.095310 - 0.086777 = 8.5330e-03 ,
c'est-à-dire très proche de la def méca (1,1) donnée par Herezh
Donc je pense que le résultat est conforme aux choix fait dans l'implantation informatique.
Je suppose que si on utilise pour le 1) le choix: eps_totale est une def de log, on aura une def_log_meca(1,1) nulle.
Par contre ça ne marchera pas avec une loi de Hooke qui utilise Almansi. Il faut que toutes les lois utilisent du log.
Donc avec les choix actuels on observe pour des grandes def thermique (10% ici) une deff finale qui est de l'ordre de la diff entre mesure d'Almansi et mesure Log.
a suivre...
Mis à jour par Julien Troufflard il y a presque 2 ans
- % réalisé changé de 20 à 0
merci Gérard. (ps : je n'ai pas encore lu ton message précédent, j'étais en train d'écrire).
si j'ai bien compris, il y aurait 2 pistes pour assurer un bon comportement de la dilatation thermique dans Herezh++.
Le point de départ est d'origine théorique (physique) : la dilatation doit produire une variation de volume en fonction d'une variation de température. Ensuite, il faut que ça suive en terme de mesure de déformation, peu importe le choix. Et il faut que ça suive également en terme de séparation "partie méca" / "partie thermique" dans les tenseur déformation.
possibilités :
1) on dit que la dilatation crée une variation de volume d'origine thermique V_total = V_meca + V_thermique
mais je ne vois pas trop l'intérêt, car il faudra convertir ça en déformation (calcul de EPS pour ISOELAS, calcul de Q_eps pour Favier, etc...)
et il faut que tout soit cohérent : si on a un certain Q_eps (pour Favier), il faut que ça soit cohérent avec V_meca et EPS (qui la déf méca).
2) donc autant appliquer basiquement la définition de la dilatation thermique avec coefficient linéaire. La dilatation est une relation linéaire donnant la variation de longueur dL en fonction de la dimension initiale L0, du coef de dilatation linéaire alpha et de la variation de température dT :
dL = alpha*L0*dT
=> autrement dit, la dilatation thermique génère une déformation de type ingénieur : dL/L0 = alpha*dT
et donc, en 3D, une variation de volume égale à dL_1*dL_2*dL_3 = dV
d'où la généralisation 3D donnant le tenseur de def ingénieur d'origine thermique eps_ing_th en fonction du tenseur identité I :
eps_ing_th = alpha*dT*I
ce tenseur serait défini dans le repère global. Et j'imagine que cela est raccord avec ce que l'on appelle justement le tenseur de dilatation (Cauchy-Green gauche ???).
Le travail à faire serait alors de convertir correctement cette déformation selon les lois et repère : en déf Almansi pour ISOELAS ou calcul Q_eps pour Favier.
Mais globalement, tout ça est compliqué. Par exemple avec ce choix de eps_ing_th, je parie que ISOELAS va générer une EPS méca non nulle. ça ne parait pas évident de rendre tout ça cohérent.
Mis à jour par Julien Troufflard il y a presque 2 ans
merci pour l'explication concernant Hypo. C'était évidemment une histoire de différence de mesure de déf Almansi/log. Je comprends pas comment je suis passé à côté, j'ai fait également des conversions qui ne donnait pas ces résultats ??? breeeeef.
Mis à jour par Julien Troufflard il y a presque 2 ans
au final, si je reviens sur ton explication pour la loi Hypo :
tu sembles dire que Herezh++ considère des déformations calculées pour une élongation de 1.1 (cf calcul de déf almansi 0.086 et déf log 0.095). Donc Herezh++ considèrerait déjà que la déf thermique est de type ingénieur ??? enfin, dans certains cas (ce n'est pas le cas pour loi ISOELAS).
Mis à jour par Gérard Rio il y a presque 2 ans
- % réalisé changé de 0 à 20
En fait d'une manière pratique dans Herezh la déformation thermique est du même type que celle qui est utilisée pour la def totale et pour la def méca. Cela est du à la formule: eps_totale - eps_thermique = eps_meca
Je m'aperçois que je n'ai pas écrit grand chose dans la doc théorique et la doc utilisateur, concernant la dilatation.
Effectivement, on peut se poser légitimement la question de : quoi t'est-ce cette alpha ...d'un point de vue théorique, en particulier en fonction de la thermo. J'ai quand même qq idées mais je vais un peu revoir avant de développer: d'ailleurs je vais en profiter pour remplir la doc théorique histoire de ne pas ce reposer les mêmes questions dans qq temps !
Je te tiens au courant
Mis à jour par Julien Troufflard il y a presque 2 ans
ok.
pas sûr que ce soit le bon choix. La déf thermique devrait être la même pour tout le monde. Si on considère une déf thermique de 0.1 (ingénieur) :
=> la loi Hypo devrait la voir comme étant log(1.1)
=> la loi ISOELAS, Favier, etc... comme étant 0.5*(1 - 1.1**-2)
=> loi de type Maheo, mooney-rivlin, etc... comme étant 1.1**2 (ou un truc du genre car c'est B barre qu'il faut considérer et non B tout court)
=> bien répercuter ça également dans les lois dépendent de V comme ISOHYPERBULK
car au final, quelle que soit la loi méca, la loi de dilatation thermique doit générer une variation de volume linéaire de pente alpha.
et après tout ça, il faudra également faire attention au cas des lois contraintes planes (la déformation 33 thermique existe).
et à voir aussi si une déformation thermique est sensée intervenir dans l'épaisseur pour une loi méca déformation plane. ??? c'est encore autre chose ça. Je n'ai pas la réponse et pas testé dans Herezh.
Mis à jour par Julien Troufflard il y a presque 2 ans
juste pour ajout à mon message précédent :
imaginons une loi additive qui combine ISOELAS, ISOHYPERBULK, FAVIER, HYPO. Dans une telle loi, ses composantes ne verraient pas toutes la même déformation thermique globale ??? On se retrouve souvent ce genre de situation. Par exemple : une loi additive ISOHYPERBULK + FAVIER + MAHEO. ou bien encore ISOHYPERBULK + FAVIER + HYSTERESIS. Et c'est d'ailleurs en comparant ce genre de loi sur une sphère sous pression que j'ai constaté des différences et démarré ce ticket.
Mis à jour par Julien Troufflard il y a presque 2 ans
2 remarques (dont une qui vaut le coup) :
1) petite erreur sur un message precedent. La dilatation thermique ne genere pas une variation de volume dL_1*dL_2*dL_3 = dV.
bon, c'est pas important. L'important c'était surtout de l'imaginer au sens d'une déf ingénieur.
2) surtout : il semble que dans Ab,aqus, un choix un peu analogue à Herezh a été fait. C'est-à-dire que la mesure de déf thermique dépend de la loi méca. Par exemple, une loi élastique (mesure log) et une loi polynomiale (mesure B) ne conduit pas à une même dilatation. D'après les déplacements du cube, dans le cas élastique, la déf thermique sera de type log. Et dans le cas polynomial, il a l'air de considérer la déf thermique comme une déf ingénieur pour le calcul méca.
Et donc, dans le cas polynomial, il y a incohérence. Par exemple, si on a une déf thermique de 0.1 (log) => on obtient une déf log totale de 0.09531 (conversion de déf ing 0.1 en log). Donc si d'un côté il y a une déf thermique log de 0.1 et en résultat une déf log totale de 0.09531, ça implique qu'il y aurait une partie méca non nulle (et pourtant des contraintes nulles).
je vais refaire le test en mettant 2 cubes en même temps dans le même calcul avec chacun une loi différente.
Mis à jour par Julien Troufflard il y a presque 2 ans
et donc oui, c'est bien ça.
=> en mettant 2 cubes en même temps dans le même calcul qui subissent la même évolution thermique, même conditions limites (dilatation libre). Pour les 2 cubes, la déformation thermique est de 0.1. En résultat, l'un des cubes finit avec une déformation totale log de 0.1 (loi élastique) et l'autre une déformation totale log de 0.09531 (loi polynomiale).
Mis à jour par Gérard Rio il y a presque 2 ans
- % réalisé changé de 20 à 30
merci pour tes retours, on voit que ce n'est pas tout simple !
J'ai commencé à rédiger une partie théorique sur la dilatation thermique, que j'ai déposé sur le site: cf. chap. 24
Cela m'amène aux réflexions suivantes:
- comment est mesuré expérimentalement le coefficient ?
- Une fois ceci posé, on peut transposer ce coefficient au cas de mesure utilisé pour le calcul à savoir dans Herezh la mesure d'Almansi.
Si par exemple le coefficient est mesuré expérimentalement via une mesure log, il faut le traduire pour une utilisation avec Almansi.
Dans Herezh on peut utiliser une fonction pour définir alpha donc a priori tout type de traduction est possible.
Je crois que c'est peut-être un début de réponse.
a suivre
Mis à jour par Julien Troufflard il y a presque 2 ans
tout à fait. Il faut voir la loi thermique comme n'importe quelle loi, c'est-à-dire qu'elle est définie dans une certaine mesure de déformation dans le code de calcul et l'utilisateur, lui, il utilise ce qu'il veut pour identifier les coefs via les conversions qu'il veut.
Ce qui pose problème, c'est quand le code change la mesure de déf de la loi thermique selon les cas sans qu'on ait la main dessus.
Mais alors du coup, il me semble que Herezh ne soit pas dans ce cas. La déformation thermique est toujours almansi malgré ce que laissait suggérer mes messages précédents.
J'ai refait un calcul Hypo mais cette fois en un seul pas de temps. Et j'obtiens bien une déf totale almansi de 0.1 pour une déf thermique almansi de 0.1. J'ai alors tenté en divers pas, et du coup, ça donne ça :
1 pas : tot= 9.999999998986e-02 meca= -1.014471839892e-11
2 pas : tot= 9.499999999791e-02 meca= -5.000000002088e-03
4 pas : tot= 9.274687499193e-02 meca= -7.253125008070e-03
8 pas : tot= 9.167409816761e-02 meca= -8.325901832386e-03
20 pas : tot= 9.104653120118e-02 meca= -8.953468798821e-03
50 pas : tot= 9.079877466174e-02 meca= -9.201225338260e-03
100 pas : tot= 9.071659765251e-02 meca= -9.283402347490e-03
1000 pas : tot= 9.064281177681e-02 meca= -9.357188223187e-03
étonnant, on dirait donc une erreur due à l'intégration. Mais contrairement à d'habitude, plus le pas de temps est petit, plus l'erreur est grande. L'erreur évolue toutefois de manière asymptotique, heureusement elle n'explose pas.
=> mais ça reste un point à regarder. Sur un cas aussi simple (irrotationnel), normalement il ne devrait pas y avoir de problème lié à l'intégration de D.
j'ai essayé avec la loi ISOHYPERBULK :
1 pas : tot= 8.023350651512e-02 meca= -1.976649348488e-02
10 pas : tot= 8.023350651496e-02 meca= -1.976649348504e-02
50 pas : tot= 8.023350651497e-02 meca= -1.976649348503e-02
100 pas : tot= 8.023350651710e-02 meca= -1.976649348290e-02
1000 pas : tot= 8.023350651631e-02 meca= -1.976649348369e-02
=> aucune dépendance au pas de temps. Mais le résultat est largement erroné. Je me demande quel calcul est fait pour cette loi qui dépend de V=v/v0. Peut-être que Herezh++ utilise la trace de déf almansi pour évaluer un v/v0 thermique ? si oui, ça expliquerait l'erreur. Il faut vraiment modifier ça.
La conclusion importante, c'est quand même que Herezh++ a l'air de considérer toujours une déformation almansi thermique, quelque soit la loi méca contrairement à ce que je pensais au départ. Mais alors il faut comprendre au cas par cas l'origine des erreurs selon les lois et améliorer.
Mis à jour par Gérard Rio il y a presque 2 ans
je continue de rédiger (et d'approfondir) la partie théorique.
Je me suis aperçu que j'ai négligé un terme dans l'opérateur tangent de la contrainte: dépendance à la dilatation via la métrique, cf. la documentation théorique que j'ai mise à jour.
Cela pourrait expliquer pourquoi tu as certaine fois de pb de convergence.
J'ai réfléchi pour y remédier... ça demande un peu de murir mais je crois que c'est faisable
Pour les calculs que tu présentes, peut-être que c'est normal:
- avec un seul pas de temps, c'est comme si tu faisais du hooke
- avec plein de pas, c'est comme si la déformation utilisée pour la méca était une def cumulée, car dans ce cas c'est l'accroissement de def thermique qui est utilisé et non la def totale ... cf. les formules que j'utilise. Là aussi il faut que je développe au niveau théorique.
Je n'oublie pas les pb avec les autres lois !!!
Mis à jour par Julien Troufflard il y a presque 2 ans
Ok.
Pour la loi isohyperbulk, c'est bien la trace de eps thermique (almansi) qui est utilisée et retranché au V total.
Verifié avec la formule -p = V.dw/dV + w, je trouve environ -p=1e-11 comme herezh si je calcule w en retranchant 3*alpha*delta(T) à VOLUME_ELEMENT (v0=1 car cube unitaire)Mis à jour par Gérard Rio il y a presque 2 ans
effectivement pour les lois hyper type grenobles, je modifie la trace pour prendre en compte la variation de volume.
Mais.... j'utilise aussi le tenseur de déformation mécanique c-a-d sans la def thermique !! Dommage ! je crois que c'est pour cela que les résultats ne sont pas correcte.
Du coup je me suis plongé dans la biblio de (Guelin, Favier Pégon...) car il y a des choses que je veux éclaircir et en particulier revoir les interactions avec la thermo...
Bref, un chantier !! mais j'avance