r3.06.07 Diagonalisation de la matrice de masse thermique#
Résumé :
Pour améliorer la régularité de la solution dans les problèmes de thermique transitoire, une des solutions consiste à « lumper » (c’est-à-dire: condenser sur la diagonale) la matrice de masse thermique (matrice de capacité).
Cette possibilité est accessible par les modélisations PLAN_DIAG, AXIS_DIAG et 3D_DIAG pour le phénomène THERMIQUE. Elle est activée lors de l’appel aux commandes de calcul thermique THER_LINEAIRE et THER_NON_LINE.
Lorsqu’on utilise ces modélisations, seuls les éléments finis linéaires (2D et 3D) ont une matrice de masse lumpée. En effet, la diagonalisation directe ne donne pas de résultats satisfaisants pour les éléments finis quadratiques donc on ne le permet pas
Les résultats théoriques sont illustrés par le calcul thermo-mécanique d’un cylindre soumis à un choc thermique.
Principe du maximum#
Énoncé du principe pour le cas continu#
On donne ici un des énoncés possible du principe du maximum pour l’opérateur de la chaleur (en l’absence de termes de source, et en thermique linéaire homogène isotrope) [bib2] .
Soit \(\Omega\) un ouvert borné de \({ℝ}^{n}\) de frontière \(\Gamma\) , dont l’adhérence est notée \(\stackrel{ˉ}{\Omega}\) .
Soit \(u(\text{x,t})\) telle que :
\(\frac{\partial u}{\partial t}-\Delta u\le 0\text{sur}\Omega \times ]0,T[,(T>0)\)
de classe \({C}^{2}\) par rapport à \(x\) et \(u\) de classe \({C}^{1}\) par rapport à \(t\) sur \(\Omega '\times ]0,T[\)
Alors \(\underset{\overline{\Omega}\times \left[0,T\right]}{\text{Max}}u=\underset{P}{\text{Max}}u\) , où \(P=(\overline{\Omega}\times \left\lbrace 0\right\rbrace )\cup (\Gamma \times \left[0,T\right])\) est la frontière du cylindre \(\Omega \times ]0\text{,T}[\) .
Ce résultat assure donc que le maximum de \(u\) est forcément atteint soit lors des conditions initiales soit sur un bord du domaine au cours du transitoire.
Respect du principe du maximum au niveau discret#
On considère l’équation de la chaleur (conduction thermique):
\(div(\lambda .\nabla T)+s(x,t)={\mathrm{\rho C}}_{p}\frac{\partial T}{\partial t}\) + Conditions limites + condition initiale \(T({t}_{0}\text{,x})={T}^{0}(x)\)
avec |
\(T\) |
température |
\(s\) |
chaleur par unité de volume (sources internes) |
|
\(t\) |
variable de temps |
|
\(x\) |
variable d’espace |
|
\(\lambda\) |
coefficient de conductivité thermique |
|
\({\mathrm{\rho C}}_{p}\) |
chaleur volumique à pression constante |
Types de conditions limites (en linéaire):
Température imposée : condition limite de Dirichlet
\(T(x,t)={T}_{\text{imp}}(x,t)\text{sur}{\Gamma}_{\text{imp}}\)
Flux normal imposé : condition de Neumann définissant le flux entrant dans le domaine
\(-q(x,t).n=f(x,t)\text{sur}{\Gamma}_{\text{flux}}\)
Échange : condition limite de Fourier modélisant les échanges convectifs sur les bords du domaine
\(-q(x,t).n=h(x,t)({T}_{\text{ext}}(x,t)-T(x,t))\text{sur}{\Gamma}_{\text{échange}}\)
La formulation variationnelle du problème est la suivante : [bib3]
\(\underset{\Omega}{\int}{\mathrm{\rho C}}_{p}\frac{\partial T}{\partial t}.\mathrm{\nu d}\Omega +\underset{\Omega}{\int}\lambda \nabla T.\nabla \nu d\Omega +\underset{{\Gamma}_{\text{échange}}}{\int}\mathrm{hT}.\nu d\Gamma =\underset{\Omega}{\int}s.\nu d\Omega +\underset{{\Gamma}_{\text{flux}}}{\int}f.\nu d\Gamma +\underset{{\Gamma}_{\text{échange}}}{\int}{\mathrm{hT}}_{\text{ext}}.\nu d\Gamma\)
\(\forall v\) vérifiant \(v={T}_{\text{imp}}(x,t)\text{sur}{\Gamma}_{\text{imp}}\)
Après discrétisation en espace de cette équation, on obtient le système :
\(M\left\lbrace \frac{\partial T}{\partial t}(t)\right\rbrace +\text{KT}(t)=F(t)\) .
avec |
\(T(t)\) : vecteur des températures nodales |
|
\(M\) : matrice de masse thermique |
\(M=\sum_{e}\underset{{\Omega}_{e}}{\int}\mathrm{\rho C}N{N}^{T}\text{dV}\) |
|
\(K\) : matrice de rigidité thermique |
\(K=\sum_{e}(\underset{{\Omega}_{e}}{\int}\lambda \nabla N.\nabla {N}^{T}\text{dV}+\underset{{\Gamma}_{{\text{échange}}_{e}}}{\int}hN.{N}^{T}d\Gamma )\) |
|
\(F\) : vecteur du second membre |
\(F=\sum_{e}(\underset{{\Omega}_{e}}{\int}sNd\Omega +\underset{{\Gamma}_{{\text{flux}}_{e}}}{\int}fNd\Gamma +\underset{{\Gamma}_{{\text{échange}}_{e}}}{\int}{\mathrm{hT}}_{\text{ext}}Nd\Gamma )\) |
|
\(N\) : (fonctions de forme) |
||
Pour la discrétisation en temps, on applique une \(\theta\) -méthode \((\theta \in [0,1])\) , qui conduit à:
\((M+\theta \Delta tK){T}^{n+1}=(M+(\theta -1)\Delta tK){T}^{n}+\theta \langle {F}^{n+1}\rangle +(1-\theta )\langle {F}^{n}\rangle\)
où \({T}^{n},{T}^{n+1}\) sont les vecteurs des températures nodales aux instants \({t}_{n}{\text{,t}}_{\text{n+}1}\) .
Conditions suffisantes pour le respect du principe du maximum au niveau discret#
Une des caractéristiques de non respect du principe du maximum est l’apparition d’oscillations (temporelles ou spatiales) : si on observe la variation de la température en un nœud au cours du temps, on note que la solution oscille et dépasse les valeurs minimale et maximale déterminées par les conditions initiales et les conditions limites. Ou encore, à un instant donné, on observe des oscillations spatiales.
On cherche donc des conditions suffisantes sur \(\Delta t\) , \(K\) et \(M\) pour que la solution n’oscille pas au cours du temps ( [bib1] , [bib4] , [bib5] ). En effet, on ne sait pas obtenir de conditions nécessaires et suffisantes. On cherche donc des conditions de non oscillation de la solution au cours du temps. Si celles-ci sont vérifiées, on vérifiera que les oscillations spatiales ont également disparu, et alors le respect du principe du maximum est assuré.
Hypothèses :
Pour pouvoir exprimer ces conditions suffisantes de non oscillation, il faut ajouter deux hypothèses:
on se place au niveau élémentaire. Le respect des propriétés au niveau élémentaire suffit pour que les conditions de non oscillation soient vérifiées pour les matrices assemblées.
on considère que la matrice de rigidité \(K\) n’est formée que du terme volumique \({K}_{V}=\underset{{\Omega}_{e}}{\int}\lambda \nabla N.\nabla {N}^{T}\text{dV}\)
Cette hypothèse n’est pas valable pour toutes les conditions limites (voir paragraphe 3.1.3 ).
Les conditions suffisantes de non oscillation reviennent à exprimer des conditions sur le pas de temps et sur les termes diagonaux et extra-diagonaux de \(M\) et \(K\) pour que certaines propriétés de ces matrices soient vérifiées (basées sur la monotonie des matrices) [bib1] :
Dans le cas général, les termes extra-diagonaux peuvent être de signe quelconque. Un étude rapide permet de déterminer les conditions sur t en fonction de leurs signes pour que les équations précédentes soient vérifiées:
\({K}_{ij}\ge 0\) |
\({K}_{ij}\le 0\) |
|
\({M}_{ij}\ge 0\) |
\({M}_{ij}+\theta .\Delta {\text{tK}}_{ij}\le 0\text{}i\ne j\) [(2591)] inconditionnellement fausse sauf \({M}_{ij}={K}_{ij}=0\) |
\(\underset{i\ne j}{\max}(\frac{{M}_{ij}}{-\theta {K}_{ij}})\le \Delta t\le \underset{i}{\min}(\frac{{M}_{ii}}{(1-\theta ){K}_{ii}})\) |
\({M}_{ij}\le 0\) |
\({M}_{ij}+(\theta -1)\Delta {\text{tK}}_{ij}\ge 0\text{}i\ne j\) [(2592)] inconditionnellement fausse sauf \({M}_{ij}={K}_{ij}=0\) |
\(\underset{i\ne j}{\max}(\frac{{M}_{ij}}{(1-\theta ){K}_{ij}})\le \Delta t\le \underset{i}{\min}(\frac{{M}_{ii}}{(1-\theta ){K}_{ii}})\) |
Quels que soit \(\Delta t\) et la forme de \(M\) , il y a risque d’oscillations. |
Intervalle à respecter sur \(\Delta t\) . La diagonalisation de \(M\) permet de supprimer la borne inférieure. |
Les conditions suffisantes pour éviter les oscillations sont alors:
\({K}_{ij}\le 0\text{}i\ne j\)
\(\Delta {t}_{\min}\le \Delta t\le \Delta {t}_{\max}\)
avec:
\(\Delta {t}_{\max}=\underset{i}{\min}(\frac{{M}_{ii}}{(1-\theta ){K}_{ii}})\) et
\(\underset{}{\Delta {t}_{\min}=\underset{i\ne j}{\max}}(\frac{{M}_{ij}}{(1-\theta ){K}_{ij}},\frac{{M}_{ij}}{-\theta {K}_{ij}})\) ,
En conséquence, il faut d’abord que les matrices élémentaires vérifient \({K}_{ij}\le 0\) (c’est le cas des éléments finis linéaires étudiés plus loin).
En ce qui concerne l’intervalle sur le pas de temps:
Si les oscillations sont dues à un pas de temps trop grand \((\Delta t>\Delta {t}_{\max})\) , on peut conseiller:
soit de choisir un schéma d’intégration en temps de type implicite \((\theta =1)\) , pour éliminer la borne supérieure de l’intervalle.
soit de diminuer \(\Delta t\) . (En pratique, il est difficile de connaître un ordre de grandeur de \(\Delta {t}_{\max}\) ).
Assez souvent, le problème des oscillations se pose pour des pas de temps petits \((\Delta t<\Delta {t}_{\min})\) ; en effet, pour prendre en compte les variations de la solution (par exemple, lors d’un choc thermique), on est amené à choisir une discrétisation fine en temps. Dans ce cas, pour éviter les oscillations, on peut suggérer:
soit d’augmenter la valeur des pas de temps . En pratique, ceci n’est pas toujours possible car \(\mathrm{Dt}\) peut être imposé par la nature du problème (variation rapide du chargement). De plus, il est difficile d’avoir un ordre de grandeur de \(\Delta {t}_{\min}\) .
soit de diminuer la taille des mailles et donc augmenter le nombre d’éléments. En effet, la valeur de \(\Delta {t}_{\min}\) dépend directement de la discrétisation spatiale:
Les expressions des matrices élémentaires sont en effet:
\({M}_{ij}=\underset{{\Omega}_{e}}{\int}\rho {\text{CN}}_{i}{N}_{j}\text{dV}\)
\({K}_{ij}=\underset{{\Omega}_{e}}{\int}\lambda \nabla {N}_{i}\nabla {N}_{j}\text{dV}\)
Pour les éléments 2D, les termes de \(M\) sont donc de la forme \(\mathrm{\rho C}\times \mathrm{Surface}\) alors que ceux de \(K\) sont seulement fonction de \(\lambda\) . Cette solution reste la meilleure si on n’est pas limité par le coût calcul, car la solution thermique et surtout mécanique sera d’autant plus précise.
Soit de diagonaliser la matrice \(M\) , ce qui supprime la borne inférieure de l’intervalle.
C’est la solution proposée ici.
Dans la suite de l’étude, on s’intéresse seulement au problème des oscillations qui apparaissent pour les pas de temps trop petits: \(\Delta t<\Delta {t}_{\min}\) . On présente plus précisément la méthode de diagonalisation de la matrice \(M\) choisie, et les différents types d’éléments auxquels elle s’applique.
Méthode de diagonalisation retenue et types d’éléments#
Éléments et conditions limites tels que les termes extra-diagonaux de \(K\) soient négatifs#
On a vu que la diagonalisation de \(M\) n’est efficace que lorsque les termes extra-diagonaux de la matrice de rigidité \(K\) sont négatifs. Dans le cas contraire, une des conditions suffisantes de non‑oscillation est inconditionnellement fausse, quelle que soit la forme de \(M\) .
Pour chacun des éléments finis utilisés en thermique dans Code_Aster, on vérifie si la matrice élémentaire de rigidité de l’élément a des termes extra-diagonaux négatifs, en s’appuyant principalement sur la [bib11] , qui donne les expressions analytiques des matrices élémentaires pour les éléments finis classiques. On résume ici les observations faites dans la [bib11] .
Éléments linéaires#
Éléments TRIA3, TETRA4, PENTA6#
La matrice élémentaire \(K\) est fonction des cotangentes des angles. Si un des angles est obtus (\(\ge 90°\) ), certains termes extra-diagonaux de \(K\) sont positifs. Si tous les angles sont aigus, la propriété est vérifiée.
On a le même type de résultat en 3D pour le tétraèdre à 4 nœuds et le pentaèdre à 6 nœuds.
Éléments QUAD4 et HEXA8#
Certains termes extra-diagonaux de \(K\) peuvent être positifs si l’élément est trop allongé dans une direction. Sinon, la propriété est vérifiée.
On a le même type de résultat en 3D pour l’élément HEXA8.
Élément 3D pyramide à 5 nœuds#
Pour cet élément, les fonctions de forme ne sont plus des polynômes mais des fractions rationnelles en \(x,y,z\) . Pour ce type d’élément, on ne dispose pas de l’expression, même approchée de \(K\) .
Éléments quadratiques#
Élément TRIA6#
Dans \(K\) , certains termes extra-diagonaux sont nécessairement positifs.
Élément QUAD9#
De même, sur l’expression analytique des termes de \(K\) , on constate que certains des termes extra‑diagonaux sont nécessairement positifs.
Élément QUAD8#
Pour cet élément, on n’a pas l’expression complète de \(K\) pour l’élément réel. Mais pour l’élément de référence, on note que certains termes extra-diagonaux sont positifs.
Conditions limites#
La matrice \(K\) n’est pas toujours réduite au terme \({K}_{V}=\underset{{\Omega}_{e}}{\int}\lambda \nabla N.\nabla {N}^{T}\text{dV}\) . Suivant les conditions limites utilisées, un terme surfacique peut être ajouté (cas de l’échange ou du flux non-linéaire). Ce terme peut être positif.
Pour une condition d’échange par exemple, le terme \({K}_{\mathrm{échange}}=\underset{{\Gamma}_{{\text{échange}}_{e}}}{\int}hN.{N}^{T}d\Gamma\) s’ajoute à la matrice \(K\) . Il est toujours positif et ne vérifie donc pas les propriétés 2.3 énoncées ci-dessus (\({K}_{ij}\ge 0\) ). Le risque d’oscillation est donc toujours présent même en diagonalisant la matrice masse. Dans ce cas, raffiner le maillage là où la condition limite est appliquée, permet de réduire la contribution positive de la condition limite et souvent de supprimer les oscillations.
Conclusion sur les éléments : propriétés des matrices \(K\)#
Pour les éléments linéaires, si l’élément réel n’est pas trop irrégulier, les termes extra-diagonaux de \(K\) sont bien négatifs. Pour les éléments quadratiques (en 2D) et certaines conditions limites, certains termes extra-diagonaux de \(K\) sont positifs. Même en diagonalisant \(M\) , on ne peut pas assurer que la solution n’oscillera pas.
Dans Code_Aster , pour éliminer les problèmes d’oscillation et de dépassement du maximum, on diagonalise seulement les matrices de masse pour des calculs thermiques réalisés sur des éléments linéaires. Pour les éléments quadratiques, on a vu que l’on ne pouvait pas diagonaliser directement la matrice de masse. On découpe donc ces éléments en éléments linéaires qui sont eux-mêmes lumpés. Ceci est appliqué aux éléments quadratiques 2D dans Code_Aster , mais pas aux éléments quadratiques 3D, non pas pour des raisons de méthode mais parce que le découpage automatique est difficile à mettre en œuvre en 3D.
Méthode de diagonalisation : Intégration aux nœuds des éléments#
Si la matrice de masse élémentaire est calculée par intégration numérique, ses termes s’écrivent sous la forme [bib8]
\({M}_{ij}=\underset{{\Omega}_{c}}{\int}{\mathit{\rho CN}}_{i}{N}_{j}\text{dV}\mathrm{\simeq }\sum_{q=1}^{N}{W}_{q}({\mathit{\rho CN}}_{i}{N}_{j}{)}_{q}\)
où \({N}_{i}{\mathrm{\rho CN}}_{j}\) est évalué au \({q}^{\text{ième}}\) point d’intégration
et \({W}_{q}\) est le poids d’intégration associé à ce point.
Classiquement, les points d’intégration sont les points de Gauss ; la position des \(N\) points et leur poids sont définis de telle sorte que le schéma intègre exactement les polynômes de degré \(\mathrm{2N}+1\) . Si on choisit les points d’intégration aux nœuds de l’élément , on obtient: \({M}_{ij}=0\) pour \(i\ne j\) . Cette méthode d’intégration est aussi appelée méthode de Newton-Cotes.
Remarque1 : Problèmes axisymétriques :
Si les points d’intégration sont aux nœuds, on aura, pour tout type d’élément, des masses nulles sur l’axe de symétrie. En effet, \({M}_{ij}=\int\mathrm{\rho C}{N}_{i}{N}_{j}2\pi r\text{dr}\text{dz}\) \({M}_{ij}\simeq {\delta}_{ij}2\pi {\mathrm{\rho CW}}_{i}\text{Jac}r({x}_{i})\) Si le point d’intégration \(i\) est un nœud de l’axe, \(r({x}_{i})=0\) et la masse correspondante est nulle. Pour les éléments axisymétriques, la méthode d’intégration aux nœuds n’est donc pas adaptée près de l’axe . Dans ce cas, il faut intégrer aux points de Gauss les éléments qui touchent l’axe de révolution, en utilisant la modélisation habituelle ( AXIS ). |
Remarque 2: autres méthodes possibles de diagonalisation:
D’autres méthodes sont étudiées dans [bib1] , en particulier pour essayer de diagonaliser les éléments quadratiques. En pratique, elle ne sont pas retenues à l’heure actuelle dans le Code_Aster.
Mise à l’échelle des termes diagonaux ([bib9] , [bib10] ): Hinton suggère la mise à l’échelle des termes diagonaux de la matrice \(M\) consistante, de telle sorte que la masse totale de l’élément soit conservée. On note que les masses lumpées sont toujours positives, même pour les éléments quadrangles à 8 et 9 nœuds.
Sommation par ligne ([bib10] ): On somme les valeurs de \({M}_{ij}\) par ligne et on concentre le résultat sur la diagonale. Malheureusement, ce procédé peut conduire à des masses négatives, notamment pour le quadrangle à 8 nœuds.
Remarque 3:
Pour les éléments quadratiques, on constate dans [bib1] que, même en diagonalisant avec la méthode de mise à l’échelle des termes diagonaux, on obtient des oscillations. On ne peut donc pas utiliser ces éléments dans le cadre de la diagonalisation (c’est-à-dire pour un maillage relativement grossier vis à vis de la rapidité du transitoire thermique).
On peut bien sûr utiliser les éléments quadratiques en thermique, à condition d’adapter la finesse du maillage à la raideur du choc thermique.
Mise en œuvre dans Code_Aster#
Afin d’éliminer les oscillations de la température en espace et en temps, les modélisations AXIS_DIAG, PLAN_DIAG et 3D_DIAG effectuent la diagonalisation des matrices de masse lors de calcul thermique linéaire (THER_LINEAIRE) et non-linéaire (THER_NON_LINE). Pour en garantir l’efficacité, on a vu qu’il faut l’effectuer sur des éléments linéaires.
Dans le cas où le maillage est linéaire, on effectue simplement une diagonalisation des matrices de masse par intégration aux nœuds.
Les modélisations disponibles sont donc:
Modélisations 2D#
Modélisation |
PLAN_DIAG |
AXIS_DIAG |
Maille |
Élément |
Élément |
TRIA3 |
THPLTL3 |
THAXTL3 |
QUAD4 |
THPLQL4 |
THAXQL4 |
SEG2 |
THPLSL2 |
THAXSL2 |
Commentaires sur les calculs élémentaires 2D :
Pour les éléments linéaires : les termes de masse ( matrice au premier membre et vecteur au second membre) sont lumpés par intégration aux nœuds. Les nouveaux éléments ont des options de calculs élémentaires identiques aux éléments classiques. Les seules options élémentaires modifiées sont donc MASS_THER et CHAR_THER_EVOL.
En axisymétrique : si des éléments du maillage touchent l’axe, il ne faut pas intégrer aux nœuds qui se trouvent sur l’axe. Il faut donc isoler cette couche d’éléments et lui affecter la modélisation AXIS.
Modélisation 3D#
Modélisation |
3D_DIAG |
Maille |
Élément |
HEXA8 |
THER_HEXA8_D |
PENTA6 |
THER_PENTA6_D |
TETRA4 |
THER_TETRA4_D |
QUAD4 |
THER_FACE4_D |
TRIA3 |
THER_FACE3_D |
Commentaires sur les calculs élémentaires 3D :
Pour les éléments linéaires : comme en 2D, les termes de masse (matrice au 1er membre et vecteur au 2nd membre) sont lumpés par intégration aux nœuds (3eme famille de points de Gauss).
On ne diagonalise donc actuellement que les éléments 3D linéaires.
En ce qui concerne les pyramides à 5 nœuds, l’intégration aux nœuds a été essayée mais ne fonctionne pas bien. Cf. [r3.06.07-element-pyramide-noeuds] (on ne sait pas si tous les termes extra-diagonaux sont négatifs). La modélisation “3D_DIAG” n’existe donc pas pour les pyramides à 5 nœuds. De toutes façons ces éléments sont minoritaires dans un maillage 3D.
Calcul thermique d’un cylindre soumis à un choc froid#
On illustre sur un exemple numérique ce qui a été montré précédemment ; à savoir que la diagonalisation est efficace pour vérifier le principe du maximum.
On s’inspire de l’exemple industriel du refroidissement d’un coude moulé: on applique un choc thermique froid (\(289°C\) à \(20°C\) ) sur un coude fissuré. Pendant le transitoire du refroidissement, la température calculée en certains nœuds atteint \(310°\) sans diagonalisation des matrices de masse. Pour l’exemple traité ici, on se restreint à un cylindre creux de même dimension que le coude sur lequel on applique un choc thermique froid.
Données#
On étudie un cylindre creux supposé infini. Comme il n’y a pas de dépendance par rapport à \(z\) (cylindre infini), on limite l’étude à un calcul plan. Par raison de symétrie, on ne maille qu’une portion de la structure.
Coordonnées des points:
\(x\) (\(\mathrm{mm}\) ) |
\(y\) (\(\mathrm{mm}\) ) |
\(z\) (\(\mathrm{mm}\) ) |
|
\(\mathrm{M1}\) |
436.75 |
||
\(\mathrm{M2}\) |
\(436.75\cos45°\) |
\(436.75\sin45°\) |
Les calculs sont effectués sur un maillage linéaire (mailles TRIA3-QUAD4):
Caractéristiques du maillage:
Nombre de nœuds: \(90\)
Nombre et type de mailles: \(64\text{TRIA3},32\text{QUAD4}\)
Caractéristiques du matériau:
\(\lambda =19,97\text{W/m °C}\)
\(\rho \mathrm{Cp}=4.89488{10}^{6}{\text{J/m}}^{3}\text{°C}\)
Conditions limites et chargement:
Pour assurer l’invariance par rotation, on impose des conditions de flux de chaleur nul sur les faces \(\mathrm{AB}\) et \(\mathrm{CD}\) . La paroi externe est supposée parfaitement calorifugée. Sur la peau interne \(\mathrm{AD}\) , le transfert thermique entre le cylindre et le fluide est modélisé par un coefficient d’échange convectif élevé :
\(h=40000\text{W/m² °C}\) .
Le choc thermique froid appliqué sur les coudes moulés est représenté par une variation linéaire de la température du fluide circulant dans le tuyau : \(289°\to 20°\) en \(\mathrm{12s}\) . Afin d’accentuer le problème de dépassement du maximum et donc de mieux mettre en évidence l’influence du diagonalisation, on adopte un choc plus brutal : \(289°\to 20°\) en \(\mathrm{1s}\) .
On adopte la discrétisation en temps suivante:
de \(t=0s\) |
à \(t=10s\) , |
1 pas de temps |
de \(t=10s\) |
à \(t=11s\) , |
2 pas de temps |
de \(t=11s\) |
à \(t=25s\) , |
7 pas de temps |
de \(t=25s\) |
à \(t=60s\) , |
10 pas de temps |
Numériquement, la valeur retenue pour le paramètre de la discrétisation en temps est \(\theta =0,57\) .
Résultats#
Les figures suivantes montrent les profils de température dans l’épaisseur du cylindre à l’instant \(t=\mathrm{15s}\) (instant où les dépassements du maximum sont les plus grands) sans diagonalisation des matrices de masse.
On donne aussi l’évolution temporelle \(T(t)\) aux nœuds \({M}_{1}\) et \({M}_{2}\) situés à un quart de l’épaisseur de la peau interne.
Sans diagonalisation, on note que la température oscille en temps et en espace dépassent la valeur maximale de \(289°\) au début du transitoire.
Avec diagonalisation sur les éléments linéaires, on observe une solution régulière sans dépassement du maximum.
Une étude similaire a été menée sur des éléments linéaires 3D (tétraèdre à 4 nœuds, pentaèdre à 6 nœuds, hexaèdre à 8 nœuds). Les résultats conduisent aux mêmes conclusions : avec diagonalisation, les oscillations de la température disparaissent pour le calcul sur les éléments linéaires 3D.
Remarque complémentaire concernant le calcul thermo-mécanique:
Une autre étude a été effectuée dans la [bib1] pour estimer les conséquences de la diagonalisation thermique sur les résultats mécaniques. On constate que le calcul ISO-P2 (éléments quadratiques divisées en éléments linéaires, dont les matrices de masse sont lumpées) fournit des résultats satisfaisants. On élimine les oscillations spatiales de la température. Mais dans le cas étudié, avec un maillage relativement grossier, la solution mécanique reste peu précise. Bien que la solution thermique soit correcte, pour améliorer la solution en contraintes, il faut de toutes façons raffiner le maillage.
Pour les mailles TRIA3, la diagonalisation conduit à une solution régulière sans dépassement du maximum:
Pour les mailles QUAD4, la diagonalisation conduit à une solution régulière sans dépassement du maximum:
Conclusion#
Les modélisations AXIS_DIAG, PLAN_DIAG et 3D_DIAG sont proposées afin de résoudre les problèmes de dépassement du maximum avec oscillation de la solution en espace et en temps qui apparaissent lors de certains calculs thermiques transitoires avec brusque variation du chargement.
Au niveau discret, l’analyse conduit à une condition suffisante de non-oscillation sur le pas de la discrétisation en temps qui doit appartenir à un intervalle :
\(\Delta {t}_{\min}\le \Delta t\le \Delta {t}_{\max}\)
où les valeurs de \(\Delta {t}_{\min}\) et \(\Delta {t}_{\max}\) dépendent des coefficients de matrices de masse et de rigidité thermiques ainsi que du paramètre \(\theta\) de la discrétisation en temps.
En pratique, si les oscillations proviennent d’un pas de temps trop grand \((\Delta t>\Delta {t}_{\max})\) , on suggère le choix d’un schéma implicite en temps \((\theta =1)\) . Si les pas de temps sont trop petits, la diagonalisation de la matrice de masse peut permettre de supprimer les oscillations.
Pour les éléments linéaires, on montre que la diagonalisation permet effectivement d’éviter les oscillations de la solution. Pour les éléments quadratiques, une diagonalisation directe ne suffit pas à éviter les oscillations donc on interdit cette possibilité.
Bibliographie#
M.A. REDON, J.M. PROIX: Diagonalisation de la matrice de masse thermique. Note EDF/DER/HI-75/97/008/0
BREZIS : Analyse fonctionnelle, Théorie et applications. Masson (1983).
OUYANG, D. XIAO : Criteria for eliminating oscillation in analysis of heat-conduction equation by finite-element method. Communications in numerical methods in engineering. Vol.10, 453-460 (1994).
NITROSSO : Etude du principe du maximum discret - Septembre 1996. Communication interne EDF.
[bib6] O.C. ZIENKIEWICZ : La méthode des éléments finis (3ème édition). Mc Graw-Hill (1977).
[bib7] I. FRIED, D.S. MELKUS : Finite element mass matrix diagonalisation by numerical integration with no convergence rate loss. Int. Journal of solids and structures. Vol.11, 461-466 (1975).
COHEN, A. ELMKIES : Eléments finis triangulaires P2 avec condensation de masse pour l’équation des ondes. Rapport de recherche INRIA n°2418 - Septembre 1994.
HINTON, A. ROCK, O.C. ZIENKIEWICZ : A note on mass diagonalisation in related processes in the finite element method. Int. J. Earthquake Engineering and Structural Dynamics, Vol.4, 245-249 (1976).
Description des versions du document#
Version Aster |
Auteur(s) Organisme(s) |
Description des modifications |
5 |
J.M.PROIX-R&D/AMA |
Texte initial |