v4.22.002 TTNL02 - Transitoire thermique avec changement de phase#

Résumé:

Ce test élémentaire permet de traiter un problème unidirectionnel en thermique transitoire non-linéaire et de vérifier la prise en compte d’un changement de phase liquide/solide par Code_Aster en introduisant par l’intermédiaire de l’enthalpie volumique la chaleur latente de fusion. La solution est analytique et fait intervenir les fonctions d’erreur \(\mathrm{erf}\) et \(\mathrm{erfc}\) . Le problème est traité dans les cas plan et volumique.

Pour les modélisations présentées ici, les écarts des résultats obtenus par Code_Aster se situent entre 1 et 4% de la référence calculée analytiquement.

Solution de référence#

Méthode de calcul utilisée pour la solution de référence#

On dispose d’une solution semi-analytique faisant intervenir les fonctions d’erreurs :

\(\mathit{erf}(x)=\frac{2}{\sqrt{\pi}}{\int}_{0}^{x}{e}^{-{t}^{2}}\mathit{dt}\) et \(\mathrm{erfc}(x)=\frac{2}{\sqrt{\pi}}{\int}_{x}^{+\infty }{e}^{-{t}^{2}}\mathrm{dt}\)

Cette solution est valide pour un milieu semi-infini, elle ne pourra donc être utilisée que dans un domaine de variation limité de la variable de temps.

Soit

../../../../_images/Object_538.svg

la position de l’interface solide/liquide. Soient \({s}_{t}=\frac{L}{\sqrt{{t}_{\mathrm{total}}}}\) et \(\lambda =\frac{{s}_{t}}{2\sqrt{{d}_{s}}}\)\({d}_{s}\) et \({d}_{l}\) désignent la diffusivité des milieux solide et liquide \(({d}_{s}=\frac{{k}_{s}}{{c}_{s}},{d}_{l}=\frac{{k}_{l}}{{c}_{l}})\) . La solution de l’équation de la chaleur est de la forme :

\({T}_{s}(x,t)={T}_{0}+\frac{{T}_{m}-{T}_{0}}{\mathrm{erf}(\lambda )}\mathrm{erf}(\frac{x}{2\sqrt{{d}_{s}t}})\) si \(x\le {x}_{t}\)

\({T}_{l}(x,t)={T}_{i}+\frac{{T}_{m}-{T}_{i}}{\mathrm{erfc}(\lambda \sqrt{\frac{{d}_{s}}{{d}_{l}}})}\mathrm{erfc}(\frac{x}{2\sqrt{{d}_{l}t}})\) si \(x\ge {x}_{t}\)

La donnée de \({t}_{\mathrm{total}}\) suffit à définir la solution, on fixe donc \({t}_{\mathrm{total}}=420.\)

Résultats de référence#

TEMPS : Abscisse

0.5

1.0

1.5

2.0

2.5

3.0

.000

.005

682.43

661.33

647.50

638.74

632.69

628.20

.010

726.05

705.75

692.06

682.43

675.24

669.63

.015

738.11

728.70

718.44

709.60

702.23

696.06

.020

739.86

737.22

731.99

726.05

720.27

714.94

.025

739.50

737.56

734.47

730.81

727.00

.030

739.93

739.39

738.11

736.20

733.88

.035

739.99

739.88

739.45

738.61

737.40

.040

739.98

739.86

739.55

739.00

.045

739.97

739.87

739.65

.050

739.97

739.89

.055

739.97

.060

.065

.070

.075

.080

.085

.090

.095

.100

TEMPS : Abscisse

3.5

4.0

4.5

5.0

5.5

6.0

.000

.005

624.68

621.84

619.48

617.48

615.25

614.25

.010

665.09

661.33

657.43

653.65

650.37

647.49

.015

690.83

686.33

682.43

678.99

675.95

673.22

.020

710.11

705.75

701.81

698.25

709.92

692.06

.025

723.23

719.60

716.17

712.95

720.89

707.09

.030

731.34

728.70

726.05

723.43

728.48

718.44

.035

735.89

734.18

732.34

730.43

733.42

726.53

.040

738.21

737.22

736.07

734.79

736.44

731.99

.045

739.29

738.77

738.11

737.33

738.18

735.47

.050

739.74

739.50

739.15

738.71

739.12

737.56

.055

739.91

739.81

739.65

739.42

739.60

738.75

.060

739.97

739.93

739.86

739.75

739.83

739.39

.065

739.99

739.98

739.95

739.90

739.93

739.72

.070

739.99

739.98

739.96

739.97

739.88

.075

739.99

739.99

739.95

.080

739.98

.085

739.99

.090

.095

.100

(En \(°C\) , en fonction de l’abscisse en mètre et du temps en secondes).

Remarque :

On se limite aux variations pendant les 6 premières secondes, au-delà de 10 secondes la condition au limite à l’extrémité \(x=1\) n’est plus assurée.

Incertitude sur la solution#

Inconnue, due à l’évaluation des fonctions d’erreur.

Références bibliographiques#

    1. Necati Özisik - Heat Conduction - Chapter 10 : Phase-change problems example 10-3 - John Wiley & Sons.

Modélisation A#

Caractéristiques de la modélisation#

Modélisation 2D :

../../../../_images/100018E000002F3B000007C0CD6513C8CC904221.svg

Caractéristiques du maillage#

20 QUAD8

Remarque#

La chaleur latente de fusion est fournie par l’intermédiaire de l’enthalpie sur un intervalle de \(0.01°C\) .

../../../../_images/10000C4000000EAD000012B5DF895554DDF2B473.svg

Valeurs testées#

Les nœuds observés ont pour coordonnée \(y=0.0\)

Identification température

Référence

t = 0.5 s N6 (x = 0.005)

682.43

t = 1.0 s N6 (x = 0.005)

661.33

t = 3.0 s N6 (x = 0.005)

628.20

t = 6.0 s N6 (x = 0.005)

614.25

t = 0.5 s N11 (x = 0.010)

726.05

t = 1.0 s N11 (x = 0.010)

705.75

t = 3.0 s N11 (x = 0.010)

669.63

t = 6.0 s N11 (x = 0.010)

647.49

t = 0.5 s N16 (x = 0.015)

738.11

t = 1.0 s N16 (x = 0.015)

728.70

t = 3.0 s N16 (x = 0.015)

696.06

t = 6.0 s N16 (x = 0.015)

673.22

t = 0.5 s N21 (x = 0.020)

739.86

t = 1.0 s N21 (x = 0.020)

737.22

t = 3.0 s N21 (x = 0.020)

714.94

t = 6.0 s N21 (x = 0.020)

692.06

Le calcul par éléments finis nécessite une discrétisation en temps de \(\Delta t=5.{10}^{-4}s\) au moins pour les premiers pas. La condition au limite imposée à l’origine faisant passer brusquement la température de \(740.°C\) à \(580.°C\) . On observe au niveau des premiers pas de temps quelques oscillations qui se stabilisent ensuite assez rapidement, malgré tout la température maximum est dépassée, il n’y a pas respect du maximum discret. Ce phénomène est observé lors des chocs thermiques, seul un traitement numérique particulier au niveau de la matrice de masse peut remédier à ce dernier.

Modélisation B#

Caractéristiques de la modélisation#

Modélisation 3D :

../../../../_images/10001A060000352000000AA56E1CC11E010A9408.svg

Caractéristiques du maillage#

20 HEXA20

Valeurs testées#

Les nœuds observés ont pour coordonnées : \(x=y=0.005\)

Identification Température

Référence

t = 0.5 s No324 (z = 0.005)

682.43

t = 1.0 s No324 (z = 0.005)

661.33

t = 3.0 s No324 (z = 0.005)

628.20

t = 6.0 s No324 (z = 0.005)

614.25

t = 0.5 s No312 (z = 0.010)

726.05

t = 1.0 s No312 (z = 0.010)

705.75

t = 3.0 s No312 (z = 0.010)

669.63

t = 6.0 s No312 (z = 0.010)

647.49

t = 0.5 s No300 (z = 0.015)

738.11

t = 1.0 s No300 (z = 0.015)

728.70

t = 3.0 s No300 (z = 0.015)

696.06

t = 6.0 s No300 (z = 0.015)

673.22

t = 0.5 s No277 (z = 0.020)

739.86

t = 1.0 s No277 (z = 0.020)

737.22

t = 3.0 s No277 (z = 0.020)

714.94

t = 6.0 s No277 (z = 0.020)

692.06

Modélisation C#

Caractéristiques de la modélisation#

Modélisation PLAN_DIAG :

../../../../_images/100018E000002F3B000007C0CD6513C8CC904221.svg

Caractéristiques du maillage#

20 QUAD8

Remarque#

La chaleur latente de fusion est fournie par l’intermédiaire de l’enthalpie sur un intervalle de \(0.01°C\) .

../../../../_images/10000C4000000EAD000012B5DF895554DDF2B473.svg

La matrice de masse est diagonalisée pour limiter les oscillations

Valeurs testées#

Les nœuds observés ont pour coordonnée \(y=0.0\)

Identification température

Référence

t = 0.5 s N6 (x = 0.005)

682.43

t = 1.0 s N6 (x = 0.005)

661.33

t = 3.0 s N6 (x = 0.005)

628.20

t = 6.0 s N6 (x = 0.005)

614.25

t = 0.5 s N11 (x = 0.010)

726.05

t = 1.0 s N11 (x = 0.010)

705.75

t = 3.0 s N11 (x = 0.010)

669.63

t = 6.0 s N11 (x = 0.010)

647.49

t = 0.5 s N16 (x = 0.015)

738.11

t = 1.0 s N16 (x = 0.015)

728.70

t = 3.0 s N16 (x = 0.015)

696.06

t = 6.0 s N16 (x = 0.015)

673.22

t = 0.5 s N21 (x = 0.020)

739.86

t = 1.0 s N21 (x = 0.020)

737.22

t = 3.0 s N21 (x = 0.020)

714.94

t = 6.0 s N21 (x = 0.020)

692.06

Le calcul par éléments finis nécessite une discrétisation en temps de \(\Delta t=5.{10}^{-4}s\) au moins pour les premiers pas. La condition au limite imposée à l’origine faisant passer brusquement la température de \(740.°C\) à \(580.°C\) . On observe au niveau des premiers pas de temps quelques oscillations qui se stabilisent ensuite assez rapidement, malgré tout la température maximum est dépassée, il n’y a pas respect du maximum discret. Ce phénomène est observé lors des chocs thermiques, seul un traitement numérique particulier au niveau de la matrice de masse peut remédier à ce dernier.

Modélisation D#

Caractéristiques de la modélisation#

Modélisation D_PLAN_HHO linéaire :

../../../../_images/100018E000002F3B000007C0CD6513C8CC904221.svg

Caractéristiques du maillage#

20 QUAD8

Valeurs testées#

Les nœuds observés ont pour coordonnées : \(x=y=0.005\)

Identification température

Référence

t = 0.5 s N6 (x = 0.005)

682.43

t = 1.0 s N6 (x = 0.005)

661.33

t = 3.0 s N6 (x = 0.005)

628.20

t = 6.0 s N6 (x = 0.005)

614.25

t = 0.5 s N11 (x = 0.010)

726.05

t = 1.0 s N11 (x = 0.010)

705.75

t = 3.0 s N11 (x = 0.010)

669.63

t = 6.0 s N11 (x = 0.010)

647.49

t = 0.5 s N16 (x = 0.015)

738.11

t = 1.0 s N16 (x = 0.015)

728.70

t = 3.0 s N16 (x = 0.015)

696.06

t = 6.0 s N16 (x = 0.015)

673.22

t = 0.5 s N21 (x = 0.020)

739.86

t = 1.0 s N21 (x = 0.020)

737.22

t = 3.0 s N21 (x = 0.020)

714.94

t = 6.0 s N21 (x = 0.020)

692.06

Le calcul par éléments finis nécessite une discrétisation en temps de \(\Delta t=5.{10}^{-4}s\) au moins pour les premiers pas. La condition au limite imposée à l’origine faisant passer brusquement la température de \(740.°C\) à \(580.°C\).

Modélisation E#

Caractéristiques de la modélisation#

Modélisation 3D_HHO linéaire :

../../../../_images/10001A060000352000000AA56E1CC11E010A9408.svg

Caractéristiques du maillage#

20 HEXA20

Valeurs testées#

Les nœuds observés ont pour coordonnées : \(x=y=0.005\)

Identification Température

Référence

t = 0.5 s No324 (z = 0.005)

682.43

t = 1.0 s No324 (z = 0.005)

661.33

t = 3.0 s No324 (z = 0.005)

628.20

t = 6.0 s No324 (z = 0.005)

614.25

t = 0.5 s No312 (z = 0.010)

726.05

t = 1.0 s No312 (z = 0.010)

705.75

t = 3.0 s No312 (z = 0.010)

669.63

t = 6.0 s No312 (z = 0.010)

647.49

t = 0.5 s No300 (z = 0.015)

738.11

t = 1.0 s No300 (z = 0.015)

728.70

t = 3.0 s No300 (z = 0.015)

696.06

t = 6.0 s No300 (z = 0.015)

673.22

t = 0.5 s No277 (z = 0.020)

739.86

t = 1.0 s No277 (z = 0.020)

737.22

t = 3.0 s No277 (z = 0.020)

714.94

t = 6.0 s No277 (z = 0.020)

692.06

Modélisation F#

Caractéristiques de la modélisation#

Modélisation 3D_HHO constante :

../../../../_images/10001A060000352000000AA56E1CC11E010A9408.svg

Caractéristiques du maillage#

20 HEXA20

Valeurs testées#

Les nœuds observés ont pour coordonnées : \(x=y=0.005\)

Identification Température

Référence

t = 0.5 s No324 (z = 0.005)

682.43

t = 1.0 s No324 (z = 0.005)

661.33

t = 3.0 s No324 (z = 0.005)

628.20

t = 6.0 s No324 (z = 0.005)

614.25

t = 0.5 s No312 (z = 0.010)

726.05

t = 1.0 s No312 (z = 0.010)

705.75

t = 3.0 s No312 (z = 0.010)

669.63

t = 6.0 s No312 (z = 0.010)

647.49

t = 0.5 s No300 (z = 0.015)

738.11

t = 1.0 s No300 (z = 0.015)

728.70

t = 3.0 s No300 (z = 0.015)

696.06

t = 6.0 s No300 (z = 0.015)

673.22

t = 0.5 s No277 (z = 0.020)

739.86

t = 1.0 s No277 (z = 0.020)

737.22

t = 3.0 s No277 (z = 0.020)

714.94

t = 6.0 s No277 (z = 0.020)

692.06

Modélisation G#

Caractéristiques de la modélisation#

Modélisation D_PLAN_HHO constante :

../../../../_images/100018E000002F3B000007C0CD6513C8CC904221.svg

Caractéristiques du maillage#

20 QUAD8

Valeurs testées#

Les nœuds observés ont pour coordonnées : \(x=y=0.005\)

Identification température

Référence

t = 0.5 s N6 (x = 0.005)

682.43

t = 1.0 s N6 (x = 0.005)

661.33

t = 3.0 s N6 (x = 0.005)

628.20

t = 6.0 s N6 (x = 0.005)

614.25

t = 0.5 s N11 (x = 0.010)

726.05

t = 1.0 s N11 (x = 0.010)

705.75

t = 3.0 s N11 (x = 0.010)

669.63

t = 6.0 s N11 (x = 0.010)

647.49

t = 0.5 s N16 (x = 0.015)

738.11

t = 1.0 s N16 (x = 0.015)

728.70

t = 3.0 s N16 (x = 0.015)

696.06

t = 6.0 s N16 (x = 0.015)

673.22

t = 0.5 s N21 (x = 0.020)

739.86

t = 1.0 s N21 (x = 0.020)

737.22

t = 3.0 s N21 (x = 0.020)

714.94

t = 6.0 s N21 (x = 0.020)

692.06

Le calcul par éléments finis nécessite une discrétisation en temps de \(\Delta t=5.{10}^{-4}s\) au moins pour les premiers pas. La condition au limite imposée à l’origine faisant passer brusquement la température de \(740.°C\) à \(580.°C\).

Synthèses des résultats#

L’erreur obtenue par rapport à la solution analytique reste raisonnable pour les points d’observation listés dans les tableaux. Signalons toutefois que le choc thermique imposé au début du transitoire provoque des oscillations (lorsqu’on observe la variation de la température en un point au cours du temps) qui s’amortissent rapidement et qui ont disparues au temps \(t=0.5s\) .