METHODES NUMERIQUES APPLIQUEES AUX … relative au transfert de chaleur par conduction et en régime non stationnaire sera l’exemple à résoudre durant toute la première partie

METHODES NUMERIQUES APPLIQUEES AUX CALCULS DES ECOULEMENTS ET DU TRANSFERT DE CHALEUR

(Version 1, Juin 2011)

Par : Pr. Abbès AZZI

Faculté de Génie-Mécanique USTO MB

BP.1505, El-Mnaouar, 31000, Oran, Algèrie.

Tel-fax:+213 (0) 41 416121

e-mail: [email protected] url : www.abbesazzi.com

T3028262422201816141210

21 x 21


www.abbesazzi.com Page 2

AVANT PROPOS

La pierre angulaire de la méthode des différences finies, est bel est bien le développement en série de Taylor. Brook Taylor, cet élève qui devint plus célèbre que ces professeurs, découvrit les séries appelées ‘développement de Taylor’. Par sa découverte, Taylor a mis entre nos mains le moyen de prédire la valeur d’une fonction en un point donné en fonction de sa valeur et la valeur de ces dérivées en un autre point tout proche du premier.

C’est bien à partir de cette série, qu’on peut obtenir les schémas algébriques pour remplacer les dérivées dans une équation de type EDP (Equation aux Dérivées Partielles). C’est la base même de la méthode des différences finies et des autres méthodes déduites de celle-ci. Tout le reste n’est qu’annexes servant à parler de stabilité, consistance, erreurs de troncature et autres.

Vous l’aurez compris, toute la philosophie de cette méthode est d’essayer de prédire ce qui se passerait dans un laps de temps sur la base de ce qui se passe à l’instant (valeur instantanée) et les tendances de changement actuelles (les dérivées successives). Ceci est vrai pour le temps mais aussi pour l’espace. Cette prédiction est d’autant plus juste que l’incrémentation est petite et/ou que les lois de changement et d’évolution sont connues.

Mon cours de différences finies, je le divise habituellement en trois grands chapitres classés par ordre de complexité. J’aime aussi construire mon cours autour d’exemples à résoudre ce qui permettra d’apprendre tout en appliquant.

Il est aussi important de dire que les équations de transport dont il est question en MDF, comportent essentiellement un terme non stationnaire, un terme de transport par convection, un terme de transport par diffusion et enfin un terme source.

La partie diffusion est la plus simple à traiter, puisqu’en générale le coefficient de diffusion est assimilé à une constante, d’où une équation linéaire plus simple à traiter. L’équation de Fourier, relative au transfert de chaleur par conduction et en régime non stationnaire sera l’exemple à résoudre durant toute la première partie du cours. Dans cette partie il est question d’introduire l’étudiant aux schémas numériques de base aussi bien pour l’espace que pour le temps. Les notions de précision (erreurs de troncature), de stabilité et de consistance compléteront cette première partie.

Dans un deuxième temps, la partie diffusion sera retraitée par l’approche des volumes finis. Les mêmes exemples seront repris et discutés sur la base de cette méthode. Tout comme pour les différences finies, la méthode des volumes finis repose sur un principe de base qui est le théorème de la divergence. Ce principe permet de substituer une intégrale de volume par une intégrale de surface. Cette partie du cours correspond à ce que je donne habituellement aux étudiants de graduation.



Les termes de convection sont non linéaires et par conséquent plus compliqués à traiter. Il s’agit là, d’un mouvement macroscopique de fluide, à qui on doit adapter les schémas de convection en fonction de la direction de l’écoulement. Cette partie sera traitée directement par la méthode des volumes finis et portera sur la dualité précision-stabilité. Les différents types de schéma et leurs propriétés seront étudiés à travers des exemples d’applications. En générale, je réserve cette partie pour les étudiants de post-graduation, mais n’empêche que des fois avec des étudiants studieux en graduation, on peut aborder une partie de ce chapitre.

La troisième partie du cours, concerne la résolution des systèmes d’équation (Navier-Stokes). A travers ce système d’équations quasi-non linéaires et couplées j’introduis les algorithmes de correction de pression utilisés pour les équations de fluides incompressibles. La partie compressible ne fait pas encore partie de ce cours.



SOMMAIRE

Les équations aux dérivées partielles, classification

PARTIE I :

1. Présentation de la méthode des différences finies 1. L’équation de conduction de la chaleur (Joseph Fourier) 2. Le problème stationnaire 3. Le problème non stationnaire 4. Schémas explicite et implicites 5. Le concept de stabilité (transformation de Fourier) 6. Schéma de Crank-Nicholson 7. Schéma de Duffort-Frankel 8. Le concept de consistance 9. Mini-projet (conduction thermique en 2D)

2. Présentation de la méthode des volumes finis 1. Application à la partie diffusion (1D) 2. Diffusion en 2D et 3D 3. Mini-projets (conduction thermique en 2D)

PARTIE II :

1. Application de la méthode des volumes finis pour un problème de convection-diffusion

1. Les propriétés d’un schéma de convection 2. Schéma avant d’ordre un 3. Schéma centré d’ordre deux 4. Schéma hybride 5. Schémas à haute précision avec et sans limiteurs

PARTIE III :

1. Algorithme de couplage pression-vitesse 1. Relaxation 2. Maillage décalé 3. Interpolation de Rhie & Show (maillage colocatif) 4. Algorithmes : SIMPLE, SIMPLEC, SIMPLER et PISO



LES EQUATIONS AUX DERIVEES PARTIELLES (EDP)

Dans cette première partie il est question de proposer un classement des équations aux dérivées partielles de la mécanique des fluides et des conditions aux limites qui vont avec.

Classification :

Considérons la forme générale d’une Equation aux Dérivées Partielles (EDP) de second ordre suivant les deux variables indépendantes (x et y) :

02

22

2

2

GF

yE

xD

yC

yxB

xA

(1)

Une classification assez simple de cette équation peut être faite sur la base des coefficients associés aux dérivées d’ordre le plus élevé A, B et C. On calcule le déterminant définit par :

CAB 42 L’équation est dite de type

elliptique si 0 , parabolique si 0 , hyperbolique si 0 .

Dans le cas d’un système d’EDP, il faut écrire l’équation caractéristique du système pour trouver sa nature. La marche à suivre est illustrée par l’exemple suivant :

11111 EyVD

xVC

yUB

xUA

(2)

22222 EyVD

xVC

yUB

xUA

(3)

on écrit les déplacement :

dyy

Udxx

UdU

(4)

dyyVdx

xVdV

(5)

Les équations précédentes s’écrivent sous la forme compacte suivante :



dV

dU

E

E

yVxVy

Ux

U

dydx

dydx

DCBA

DCBA

2

1

2222

1111

00

00 (6)

Le déterminant :

02122112211221

21221 dxDBDBdydxCBCBDADAdyCACA (7)

On divise l’équation précédente par 2dx , et on définit dxdyf '

0'2' cfbfa (8)

cab 42 (9)

L’équation est dite de type elliptique si 0 , elle est parabolique si 0 , et hyperbolique si 0 . Une des utilités de cette classification est de prévoir le comportement de l’équation vis à vis des conditions aux limites. Si nous imaginons un écoulement de fluide de gauche vers la droite, une perturbation en un point donné n’a pas d’influence amont si l’équation est de type parabolique. Si par contre l’équation est de type elliptique une perturbation quelconque en un point quelconque aura une influence dans toutes les directions de l’espace. Une conséquence directe de cette caractéristique est qu’un problème de type parabolique peut être résolu par une marche avant, alors qu’une équation de type elliptique nécessite la prise en considération des conditions aux limites imposées sur toutes les frontières du domaine de calcul. Par exemple :

L’équation de Laplace 02

2

2

2

yx

elliptique

L’équation de diffusion 02

2

xt

parabolique

L’équation 02

2

2

2

yx

hyperbolique

L’EDP de nature parabolique : C’est le cas d’un problème de propagation associé à un mécanisme de dissipation tel que la conduction thermique non stationnaire.



L’équation 2

2

xt

liée aux conditions initiales : x sin et aux conditions aux

limites : 0,1,0 tt accepte la solution exacte suivante txtx 2expsin, C’est une équation linéaire d’ordre 2, parabolique par rapport à la variable du temps t. La propagation en avant dans le temps et la diffusion dans l’espace, font que la solution en un point P peut influencer n’importe qu’elle point pour itt . Cependant les points se situant

dans la zone itt ne sont pas influencés par la solution au point P. En d’autres termes on

dira que le passé influe sur le futur alors que l’inverse n’est pas vrai. La dissipation dans l’espace, fait que même si la distribution initiale pour 0t est discontinue, la solution devient continue pour 0t . L’EDP de nature elliptique : Cette catégorie d’EDP est associée aux problèmes de nature stationnaire ou d’équilibre tels que l’écoulement stationnaire d’un fluide visqueux, la répartition stationnaire du champ de température ou la distribution d’un potentiel.

L’équation de Laplace du type 02

2

2

2

yx

, associée aux conditions aux limites suivantes

xx sin0, , expsin1, xx et 0,1,0 yy accepte la solution exacte

suivante : yxyx expsin, La principale caractéristique de ce type d’équation elliptique est qu’une perturbation introduite en un point quelconque à l’intérieur du domaine de calcul influe sur la totalité du domaine. Ceci implique que pour résoudre un problème de type elliptique il est impératif de poser les conditions aux limites sur toutes les frontières du domaine. Ici aussi une discontinuité dans les conditions aux limites est rapidement effacer (lisser) à l’intérieur du domaine de calcul. L’EDP de nature hyperbolique : Cette catégorie d’EDP peut être considérée comme extension des équations elliptiques pour lesquels certaines valeurs critiques des paramètres doivent être déterminées en même temps que la distribution d’équilibre correspondante. La résonance de circuit électrique ou d’enceintes acoustiques ainsi que la détermination des fréquences propres des structures élastiques constituent des exemples de ce type d’équations.

L’équation de propagation d’une onde suivante 2

2

2

2

xt

représente un très bon exemple

pour l’équation de type hyperbolique. Cette équation associée aux conditions initiales xx sin0, , 00, xt et aux conditions aux limites 0,1,0 tt accepte la

solution suivante : txtx cossin,



Enfin, la figure 1 représente schématiquement l’influence d’une perturbation au point P sur l’ensemble du domaine de calcul pour les trois types d’équations.

Figure 1 : Nature des équations et conditions aux limites,

(a) Hyperbolique, (b) Parabolique et (c) Elliptique.

Les conditions aux limites Soit un problème définit dans un domaine R, limité par la frontière R . Les conditions aux limites peuvent être de trois natures : Dirichlet : Dans ce type de conditions la valeur de la variable dépendante est imposée sur la frontière du domaine de calcul

Rsurf (10) Newman : La variable dépendante n’est pas connue sur la frontière mais sa dérivée est bien définit

Rsurqs

oufn

(11)

Mixte : Une combinaison linéaire des deux premières conditions est imposée sur la frontière

Rsurkfkn

0,

(12)

Un problème de transfert de chaleur ou d’écoulement est dit bien posé si en résolvant les équations du problème liées aux conditions aux limites et initiales

La solution numérique existe. La solution numérique est unique. La solution numérique dépend de façon continue de la variation des conditions

aux limites.



PRESENTATION DE LA METHODE DES DIFFERENCES FINIES

TAYLOR BROOK (1685-1731) Sir Brook Taylor est un homme de sciences anglais, né le 18 août 1685 à Edmonton (Angleterre). Il est mort à l’âge de 46 ans, le 29 décembre 1731 à Londres. Son domaine d’intérêt inclus les mathématiques, la musique, la peinture et la philosophie. Son amour pour les mathématiques, lui a été transmis par ces professeurs John Machin et John Keill. Il compléta ses études à l'université de Cambridge et devint célèbre pour ses contributions au développement du calcul infinitésimal. Le 3 avril 1712 (à l’âge de 27 ans), Taylor fut admis à la Royal Society de Londres (l'équivalent de l'Académie des sciences de Paris), non sur la base de ces publications scientifiques mais sur recommandation et expertise de Machin, Keill et autres. Environ deux années après il fut élu secrétaire de la Royal Society, et il y resta du 14 janvier 1714 au 21 octobre 1718, lorsqu'il dut se résigner pour raisons de santé d'une part, d'autre part par manque de motivation. La période où il fut secrétaire de la Royal

Society de Londres fut celle de sa vie où il fut le plus productif en mathématiques. Il publia deux ouvrages en 1715, qui sont extrêmement important pour l'histoire des mathématiques. Dans son ouvrage, Methodus incrementorum directa et inversa, Taylor ajouta aux mathématiques supérieures une nouvelle branche appelée ‘calcul de différences finies’, inventa l'intégration par parties, et découvrit les séries appelées ‘développement de Taylor’. En fait, la première mention par Taylor de ce qui est appelé aujourd'hui théorème de Taylor apparaît dans une lettre que ce dernier écrivit à Machin le 26 juillet 1712. L'importance du théorème de Taylor ne fut pas perçue avant 1772 quand Lagrange proclama que c'était le principe de base du calcul différentiel. Le terme ‘série de Taylor’ semble avoir été utilisé pour la première fois par L'Huilier en 1786. Taylor présenta aussi les principes de base de la perspective dans Linear Prospect (1715). La seconde édition fut appelée New principles of linear perspective. Enfin, Taylor fit de nombreux séjours en France. C'était d'une part suite à des problèmes de santé et d'autre part pour garder le contact avec ces amis mathématiciens. Actuellement, la pierre angulaire de la méthode des différences finies n’est autre que le développement des séries de Taylor. (Wikipédia Encyclopédie)

La méthode des différences finies : Cette méthode est basée sur la technique du

développement en séries de Taylor qui permet d’approximer la valeur d’une fonction en un point donné si on connaît la valeur de la dite fonction ainsi que toute ces dérivées en un point voisin en espace ou en temps. Cette technique permet de développer des schémas pour remplacer les dérivées premières et secondes des EDP pour pouvoir envisager une solution numérique par calculateur.

Pour obtenir une solution numérique il faut tout d’abord définir un domaine numérique constitué par un ensemble de points discrets appelé grille de calcul. Les valeurs instantanées et locales des variables dépendantes du problème sont définit sur l’ensemble des points de la grille de calcul. La différence entre cette vue numérique à travers un certain nombre de points et la distribution continue exacte représente l’erreur commise par la méthode numérique. Il est tout à fait logique de penser que plus le nombre de points est important plus la visualisation est claire, un peu comme les pixels d’une photo numérique. La Figure 1 représente des exemples de grilles de calcul.



Figure 1. Exemples de grilles de calcul

L’étape suivante consiste à approximer ou remplacer toutes les dérivées partielles par des schémas discrets (différence finies). L’EDP sera transformée en équation algébrique. Cette équation algébrique est ensuite appliquée sur l’ensemble des nœuds de la grille de calcul. Le résultat sera un système d’équation comportant autant d’équations que d’inconnues (nœuds). Ce système sera ensuite résolu par une méthode appropriée. Le résultat sera une distribution discrète de la solution sur l’ensemble des points du domaine de calcul.



Grille de calcul :

Figure 2 ; Grille de calcul structurée 2D.

Avant de commencer, il faut trouver un moyen qui nous permettra de localiser spatialement et temporellement tous les points de la solution numérique. C’est ce qu’on va appeler création de la grille de calcul. Dans la suite, on va résonner sur un espace plan (2D) et l’extension pour le 3D sera faite de manière intuitive. La Figure 2 représente la manière la plus directe pour repérer les points suivant la procédure structurée. C’est un peu comme une matrice, chaque point sera affecté de deux indexes (i,j) qui le positionneront par rapport à ces voisins. Soit U, la variable à calculer. Sa valeur aux différents points de la grille s’écrit de la manière suivante :

),( 00,1 yxxUU ji (1)

),( 00,1 yxxUU ji (2)

),( 001, yyxUU ji (3)

),( 001, yyxUU ji (4)

Maillage non-structuré : L’autre façon de mailler un domaine de calcul est de définir un

nuage de points, pas nécessairement structuré. Dans ce cas-là, il faudra numéroter les points de calcul un par un. Chaque point aura ces coordonnées x et y. En plus il faudra relier ces points entre eux de façon à créer des éléments (généralement des triangles, voir Figure 1). Le fichier de la grille de calcul sera compléter par une liste des éléments (eux-mêmes numéroter) et les points composants chaque élément.



Le développement en série de Taylor ;

n

n

n

n

Rnx

xUx

xUx

xUyxUyxxU

!

...!2

,,0

2

02

2

00000 (5)

n

n

n

n

Rnx

xUx

xUx

xUyxUyxxU

!

...!2

,,0

2

02

2

00000 (6)

Une autre écriture de l’équation (5), on oubli temporairement la deuxième dimension.

nn

iii

n

iii

iiiii Rxxn

xUxxxUxxxUxUxU 12

1

''

1'

1 !...

!2

Le terme Rn, représente les termes omis d’ordre (n+1 à l’infini). Théoriquement, on aura besoin d’un nombre infini de termes pour pouvoir calculer la valeur de U(xi+1). En pratique, on se limite à un nombre fini de terme et tout le reste sera considéré en tant que l’erreur de l’approximation (erreur de troncature).

Construction des schémas pour la dérivée d’ordre un et deux :

En arrangeant l’équation (5), on obtient le schéma aux différences avant:

xx

yxUyxxUx

U

0000

0

,, (7)

L’équation (6), donne le schéma aux différences arrière :

xx

yxxUyxUx

U

0000

0

,, (8)

Le schéma aux différences centrées s’obtient en soustrayant l’équation (6) de l’équation (5) :

20000

0 2,, x

xyxxUyxxU

xU

(9)

La dérivée seconde est obtenue en additionnant l’équation (5) à l’équation (6) :

22

000000

02

2 ,,2, xx

yxxUyxUyxxUxU

(10)

Les schémas ci-dessus s’écrivent sous forme indicielle :

xx

UUx

U jiji

ji

,,1

,

(11)



xxUU

xU jiji

ji

,1,

,

(12)

2,1,1

, 2x

xUU

xU jiji

ji

(13)

22

,1,,1

,2

2 2x

xUUU

xU jijiji

ji

(14)

Application 1: A titre d’exercice, construire un schéma pour approximer la dérivée croisée.

Application 2 : En utilisant un schéma d’ordre un, calculer la première dérivée de la fonction suivante pour x=0.5 et pour deux incrémentations h=0.5 et h=0.25

2.235.015.025.0 23 xxxxf

Répéter l’opération avec un schéma d’ordre deux. Comparez avec la solution exacte.

Le fait de dire qu’un schéma est d’ordre deux veut dire qu’il est plus précis que celui d’ordre un. L’erreur de troncature est proportionnelle à h2 au lieu de h (pour le schéma d’ordre 1). De ce fait un schéma d’ordre deux est toujours préféré en CFD. Pour la dérivée par rapport au temps, il est d’usage d’utiliser un schéma avant d’ordre un. C’est un peu par rapport à la nature de la variable temps.

Erreur de troncature : C’est l’erreur qui résulte de l’utilisation d’une approximation (schéma) à la place de la solution exacte (dérivée). Erreur d’arrondi : C’est l’erreur engendrée lorsqu’on se limite le nombre de décimales pris en compte après la virgule.



La principale remarque est que le schéma centré est d’ordre 2 est plus précis que les deux autres. Malheureusement ce schéma ne peut être utilisé pour les nœuds de frontières où le domaine de calcul est définit seulement d’un seul côté du nœud de calcul.

La formule d’un schéma d’ordre 2 applicable aux nœuds des frontières peut être construite en utilisant trois points au lieu de deux. La procédure est la suivante :

2,2,1,

,

xx

UcUbUax

U jijiji

ji

(15)

...6!2

3

,3

32

,2

2

,,,1

xxUx

xUx

xUUU

jijijijiji (16)

...6

2!2

223

,3

32

,2

2

,,,2

xxUx

xUx

xUUU

jijijijiji (17)

En multiplie l’équation (16) par b et l’équation (17) par c ;

3

,2

22

,,

,2,1,

42

2 xxUbcx

xUbcxUcba

UcUbUa

jijiji

jijiji

(18)

L’identification de l’équation (18) à l’équation (15), donne :

04120

bcbc

cba (19)

La résolution de ce système d’équation, donne l’expression suivante pour un schéma de second ordre utilisant trois points pour la dérivée première.

2,2,1,

, 243

xx

UUUx

U jijiji

ji

(20)

Application 3: Construire un schéma d’ordre 2 utilisant les points, i, i+1 et i+2 pour approximer la première dérivée.



1. L’équation de conduction de la chaleur (Joseph Fourier)

L’équation de Fourrier traduisant le transfert de chaleur par conduction sera utilisée dans la suite du cours comme exemple de base pour illustrer l’application de la méthode des différences finies.

QTtTc

(21)

où : tyxT ,, : La température, fonction de l’espace et du temps.

c : La chaleur spécifique.

: La masse volumique.

Q : Source de chaleur par unité de temps et de volume.

: Le coefficient de conductivité thermique.

t : Le temps.

Bien que la conductivité thermique, la chaleur spécifique et la masse volumique peuvent varier en fonction de la température, elles seront considérées constantes dans la suite du cours.

Notre première approche du problème sera d’appliquer cette équation pour un cas assez simple tel que le transfert de chaleur en 1D. Soit un fil métallique de section droite très petite par rapport à sa longueur de façon à ce que le flux de chaleur existe seulement suivant la longueur du fil. Si en plus la source de chaleur est absente, l’équation précédente prend la forme suivante :

2

2

xTa

tT

(22)

Où

ca , représente la diffusivité thermique.

Si les températures maximale et minimale du processus sont connues, la température sera adimensionalisée comme suit :

minmax

min

TTTT

(23)

et en introduisant la variable d’espace adimensionnelle, Lxx /' , où L est la longueur du fil, l’équation précédente s’écrit :



2

2

xa

t

(24)

où x’ a été remplacée par x pour simplifier l’écriture.

LE PROBLEME STATIONNAIRE

Si en plus le problème est stationnaire, l’équation devient :

02

2

x

(25)

Le problème sera complété par la pose des conditions aux limites.

L Longueur du fil.

NI = 6 Nombre de nœuds du maillage.

Les conditions aux limites seront du type Dirichlet :

11 , 0NI (26)

On calcul x par l’expression suivante : 11 NIx (27)

et on génère la grille de calcul par la portion de programme :

x(1) = 0.0

Do I=2,NI

x(i) = x(i-1)+x

enddo

L’équation (25) sera discrétisée par un schéma centré de second ordre :

022

11

xiii

(28)

x

L

I = 1 I = 2 I = 3 I = 4 I = 5 I = 6



Le nombre de nœuds global étant 6 dont deux sont réservés pour les conditions aux limites et quatre sont à calculés par la méthode des différences finies.

L’application de l’équation algébrique (28) aux quatre nœuds donne le système suivant :

I=2 02 321 soit 12 32 (29)

I=3 02 432 soit 02 432 (30)

I=4 02 543 soit 02 543 (31)

I=5 02 654 soit 02 54 (32)

Mathématiquement parlant, on dispose d’un système de quatre équations à quatre inconnus :

0001

2100121001210012

5

4

3

2

(33)

Ce type de matrice est appelée, matrice tri diagonal et elle est facilement résolu par la méthode du pivot (triangulation).

Solution :

0011

2100121002300012

5

4

3

2

0111

2100340002300012

5

4

3

2

1111

5000340002300012

5

4

3

2

15 5 2.05

134 54 4.04

123 43 6.03

12 32 8.02



On a aussi : 11 et 06

Il est clair que la solution est une droite en parfaite concordance avec la conduction thermique uni directionnelle qui possède un caractère linéaire.

Remarque : La solution de ce type de problème est possible analytiquement (deux intégrations successives) et la solution et celle d’une ligne droite.

LE PROBLEME NON-STATIONNAIRE

On reprend l’équation (24)

2

2

xua

tu

Dans ce genre de problème, en plus des conditions aux limites on a besoin des conditions initiales. C’est à dire une distribution initiale de la solution pour le temps zéro. Les variables auront deux indices : le premier se rapportant au temps et le deuxième à l’espace.

xitnUxtU ,., sera représentée par niU .

x

L

I = 1 I = 2 I = 3 I = 4 I = 5 I = 6

n = 1

n = 2

n = 3

t



Schéma explicit

L’équation précédente sera approximer par le schéma suivant :

22

111 2 x

xUUUat

tUU n

ini

ni

ni

ni

(34)

On remarque qu’on a utilisé un schéma avant d’ordre un pour la dérivée par rapport au temps et un schéma centré d’ordre deux pour la dérivée par rapport à l’espace.

Lors de cette discrétisation nous avons choisi de prendre les termes de droites au temps n. ce schéma s’appelle un schéma explicite, puisqu’il permet de formuler l’expression de la variable au point i et à l’instant n+1 explicitement en fonction de la solution déjà calculée au temps n. Ce schéma est représenté par la molécule suivante.

L’équation (34) sera arrangée comme suit :

ni

ni

ni

ni UUUU 11

1 21 (35)

avec 2xta

(36)

L’équation (35) sera appliqué aux nœuds d’une même rangé (c.a.d. n = cste).

Reprenons le problème de conduction de la température précèdent 2

2

xt

et posons

les conditions aux limites suivantes ( 0.10, t , 0.01, t ) et les conditions initiales (

0.0,0 x pour 10 x )

Si on reprend le même nombre de nœuds que précédemment (NI=6) le pas d’espace sera 2.0x

Cas 1 : 1.0 t ( 5.2 )

x .0000 .2000 .4000 .6000 .8000 1.0000 1 1.0000 .0000 .0000 .0000 .0000 .0000 2 1.0000 2.5000 .0000 .0000 .0000 .0000



3 1.0000 -7.5000 6.2500 .0000 .0000 .0000

A ce niveau, on peut arrêter les calculs puisqu’on remarque que les résultats numériques de la prédiction ne peuvent être acceptés physiquement. En l’absence de source de chaleur les valeurs de la température doivent être bornées par les conditions aux limites, pire encore on voit apparaître des valeurs négatives de la température adimensionnelle. On conclue que le schéma numérique n’est pas stable puisqu’il amplifie les erreurs introduites par les conditions initiales.

Cas 2 : 01.0 t ( 25.0 )

x .0000 .2000 .4000 .6000 .8000 1.0000

1 1.0000 .0000 .0000 .0000 .0000 .0000 2 1.0000 .2500 .0000 .0000 .0000 .0000 3 1.0000 .3750 .0625 .0000 .0000 .0000 4 1.0000 .4531 .1250 .0156 .0000 .0000 5 1.0000 .5078 .1797 .0391 .0039 .0000 6 1.0000 .5488 .2266 .0654 .0117 .0000 7 1.0000 .5811 .2668 .0923 .0222 .0000 8 1.0000 .6072 .3018 .1184 .0342 .0000 9 1.0000 .6291 .3323 .1432 .0467 .0000 10 1.0000 .6476 .3592 .1663 .0591 .0000 D’après les résultats ci-dessus, on remarque que la première variante avec 1.0 t est

instable. Elle ne peut pas aboutir à une solution raisonnable. Alors qu’avec 01.0 t le

processus est stable. Conclusion : la stabilité d’un schéma explicite n’est pas toujours assurée.

Concept de stabilité d’un schéma :

Un schéma est dit stable s’il amorti les erreurs provenant des C.I., des C.L. et de l’approximation utilisée. S’il amplifie les erreurs, le schéma sera instable et ne pourra pas converger vers une solution réaliste.

Pour introduire le concept de stabilité nous allons utiliser le schéma de l’équation (35)

ni

ni

ni

ni UUUU 11

1 21

Soit nu la solution exacte (en minuscule) et nU la solution numérique à l’instant n. ces deux quantités seront liées par :

ni

ni

ni uuU (37)



Où niu est l’erreur introduite dans le calcul par l’approximation du schéma (erreur de

troncature).

Remplaçons l’équation (37) dans (34), nous obtenons :

22

111

,2 xt

xuuu

tuu n

ini

ni

ni

ni

(38)

Ou

211

1 ,21 xttuuuu ni

ni

ni

ni (39)

Cette dernière équation décrit l’évolution de l’erreur en fonction du temps. Comme il est dit précédemment, un schéma numérique stable ne doit pas amplifier les erreurs. Cette

conditions est bien vérifiée si 021 , puisque 2xt est toujours positif.

211

1 21 xttuuuu ni

ni

ni

ni (40)

2maxmax

1 xttuu ni

ni

En d’autres termes l’erreur introduite par un pas de temps t ne peut être supérieur à

2xtt

ANALYSE DE LA STABILITE PAR LA TRANSFORMATION DE FOURIER

Contrairement à l’erreur de troncature qui peut être estimer pour n’importe quel problème (aussi complexe soit-il), il est pratiquement très difficile d’analyser la stabilité d’un schéma donnée. Il est même impossible d’étudier la stabilité d’un algorithme pour des équations non linéaires. Une méthode d’analyse de la stabilité basée sur la transformation de Fourier peut être appliquée au schéma précèdent (35) :

ni

ni

ni

ni UUUU 11

1 21

La solution d'un tel problème peut s’écrire sous la forme suivante :

xijni etnU (41)

L’injection de cette solution dans l’équation (35) donne :

xijxijxijxij eeetnetn 11 211 (42)

Qui peut aussi s’écrire :



tnGtn 1 (43)

Où :

xjxj eeG 21 (44)

est appelé facteur d'amplification.

Pour qu’un schéma soit stable il faut que 1G .

Si on applique les égalités trigonométriques suivantes

jj ee cos2 ; 2

sin2cos1 2 (45)

à l’expression de G, nous obtenons :

xG cos221 (46)

2sin41 2 xG

(47)

Enfin : 1G sera vérifier si

12

sin411 2

x (48)

21

2sin 2

x (49)

Cette conditions est vérifié si :

210 (50)

En conclusion nous dirons que le schéma explicite étudié précédemment est stable pour la condition (50).

En analysant l’exemple cité précédemment nous constatons que l’algorithme est instable pour un 1.0 t qui correspond à 5.2 , et que nous avons stabilisé le calcul en

adoptant une valeur plus petite du pas du temps ; 01.0 t ( 25.0 ).

Les conclusions seront :

Pour un 2.0x la valeur maximale du pas du temps pour un calcul stable sera08.0t .



Si nous voulons augmenter la précision du schéma en adoptant par exemple un 1.0x , on doit aussi vérifier 02.0t

C’est à dire, plus la précision spatiale est grande plus le calcul sera plus long, puisque le pas du temps exigé pour la stabilité du schéma explicit sera plus petit. Du point de vue capacité de stockage en mémoire, ce schéma exige un espace double pour la distribution de la solution numérique ( n et n+1).

Schéma implicite

Reprenons le problème de la conduction thermique non stationnaire et re écrivons l’équation discrète (34) comme suit (les termes de droite sont au temps n+1)

22

11

111

1 2 xx

UUUatt

UU ni

ni

ni

ni

ni

(51)

Après groupement et arrangement :

ni

ni

ni

ni UUUU

1

111

1 21 (52)

Cette équation présente trois inconnus en même temps, ce qui ne permet pas de la résoudre directement comme c’était le cas pour le schéma explicite. Cette forme de discrétisation est appelée schéma implicite. Pour trouver la solution il faut écrire l’ensemble des équations issues de l’application de (52) sur tous les nœuds de la même ligne et ensuite résoudre le système tout entier.

Si nous reprenons l’exemple précédent composé de six nœuds, le système s’écrira :

121

31

2212 UUUUi nnn

nnnn UUUUi 31

41

31

2 213

nnnn UUUUi 41

51

41

3 214

651

51

4 215 UUUUi nnn

1U et 6U sont connues et représentent les conditions aux limites.



On dispose maintenant d’un système de quatre équations à quatre inconnus.

*5

*4

*3

*2

5

4

3

2

2100210

0210021

U

UU

U

UUUU

Les variables de type *iU représentent la solution numérique à l’itération précédente. La

solution de ce système donne directement la solution de l’équation. On constate que l’adoption de n’importe qu’elle valeur du paramètre aboutit à une solution numérique stable. On conclue que le schéma implicite est inconditionnellement stable.

Application 3: Utiliser l’analyse de fourrier comme précédemment pour montrer que le schéma implicite est inconditionnellement stable.

Schéma de Crank-Nickolson :

Suivant ce schéma l’équation (24) s’écrira de la manière suivante :

211

2

11

111

1 2212

21

xUUU

xUUU

atUU n

ini

ni

ni

ni

ni

ni

ni (53)

Un tel schéma prend une moitié en explicite et l’autre moitié en implicite. Une façon plus généralisée de discrétiser l’équation (24) est :

211

2

11

111

1 212x

UUUx

UUUatUU n

ini

ni

ni

ni

ni

ni

ni (54)

Pour 0 le schéma est explicite, pour 1 il est implicite et pour 5.0 il devient Crank-

Nicholson.



Schéma de Duffort Frankel

C’est un schéma explicite et inconditionnellement stable

21

111

11

2 xUUUU

atUU n

ini

ni

ni

ni

ni

(55)

Concept de consistance d’un schéma

Un schéma est dit consistant si et seulement si l’erreur de troncature tend vers zéro quand

tous les pas ix et t tendent vers zéro. En d’autres termes : plus on raffine le maillage de

calcul plus le résultat doit être précis. Le schéma implicite et explicit introduits

précédemment sont consistants puisque l’erreur de troncature 2, xt tend vers zéro

quand x et t tendent vers zéro.

Examinons le schéma de Duffort Frankel de l’équation (55).

ni

ni

ni

ni

ni

ni UUUU

xa

tUU

111

12

11

2

L’erreur de troncature a la forme suivante :

...61

122

,3

32

,2

22

,4

4

t

tU

xt

tUx

xU

jnjnjn

(56)

Tout va pour le mieux si 0lim

xt

quand 0t et 0x .

Par contre si t et x tendent vers zéro avec le même taux telle que

xt

, alors ce

schéma ne sera plus consistant.

Mini-Projets : (L’énoncé des applications ci-dessous est inspiré du cours de Lars Davidson,

Chalmers Tekniska Hogskola, Termo- och Fluiddynamik, thanks to Dr. Lars Davidson)

Le projet consiste à résoudre le problème de conduction thermique (diffusion) dans un domaine rectangulaire (2D) en appliquant des conditions aux limites de type Dirichlet et Newman.

L’équation de Fourrier :

0STGradDiv



02

2

2

2

S

yT

xT

Figure 1 : Domaine de calcul et conditions aux limites.

L =1 et H = 0.5 Lx10015

cas Sud Est Nord Ouest S

1 10 Hysin2010 10 0 xT -1.5

2 15 HyHy sin151510 10 0 xT -1.5

3 15 Hy2cos15 15 0 xT -1.5

4 10 HyHy sin10510 15 0 xT -1.5

5 15 HyHy 2cos155 10 0 xT -1.5

Tableau 1 : Les conditions aux limites du groupe 1

L =1.5 et H = 0.5 autrementyetxpour 204.03.01.17.001.0


1 10 Hysin2010 0 yT 10 0

2 10 Hysin2010 0 yT 30 0

3 10 HyHy cos1515 0 yT 15 0

y = 0

y = H

x = 0 x = L Sud

Nord

Ouest Est



4 10 HyHy cos1515 0 yT 30 0

5 10 Hysin2010 0 yT 10 0


L =1 et H = 1 20


1 10 20 Lx20 0 xT 0

2 10 20 Lx2110 0 xT 0

3 10 20 Lx515 0 xT 0

4 10 20 Lx155 0 xT 0

5 10 20 Lx5135

0 xT 0


S’inspirer de l’exemple du cas 1 ci-dessous pour adapter le programme à votre cas et présenter le rapport de votre mini-projet. Le rapport doit comporter la formulation du problème, les conditions aux limites, la discrétisation, la méthodologie utilisée, l’étude de sensibilité de la solution par rapport à la taille de la grille de calcul, les figures des résultats

(isothermes et le vecteur flux de chaleur définis par xTqx

et

yTq y

) et les

discutions.

Solution par la méthode des différences finies

SyT

xT

tT

2

2

2

2

On utilise un schéma explicit, avant pour le temps et centré pour l’espace. L’équation précédente prend la forme suivante :

nji

nji

nji

nji

nji

nji aTfTeTdcTbTT ,1,1,,1,1

1,



avec

2**xtkcb , 2**ytked , tSf * et

2***22***21 ytkxtka

Appliquer cette équation aux nœuds de la grille de calcul et obtenir un système d’équations qu’il faut résoudre par la méthode de Gauss-Seidel.

Programme Fortran à télécharger ici, Diffusion 2D en différences finies

http://www.abbesazzi.com/wp-content/uploads/2011/06/Diffusion-2D-DF.rar

Cas 1 :

Sensibilité de la solution à la taille de la grille de calcul : les calculs ont été conduits pour trois grilles ayants 10 x 10, 20 x 20 et 40 x 40 points et nommées G1, G2 et G3 respectivement. La distribution de la température pour (y = H / 2) et le long du rectangle est représenté sur la figure 2. Pour assurer la stabilité du schéma explicit il faut que le pas du temps vérifie la condition suivante :

2221

yxt

Figure 2: Sensibilité de la solution numérique vis-à-vis de la taille de la grille de calcul.

La distribution de la température sous forme de lignes isothermes est représentée sur la figure 3:

Pour illustrer le flux de chaleur, On trace les vecteurs du flux définis

par xTqx

et

yTq y

0,0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1,05

10

15

20

25

30

35

5

1 0

1 5

2 0

2 5

3 0

3 5

T

X

11x11 21x21 41x41



Figure 3: Isothermes, cas 1, 21x21 nœuds.

Figure 4: Distribution du flux thermique.

La figure 4 illustre la direction du flux thermique de conduction à l’intérieur du domaine de calcul.

Variante : Reprendre les mêmes cas en utilisant la méthode ADI pour résoudre l’équation stationnaire.

T3028262422201816141210

21 x 21

T3028262422201816141210

21 x 21



Documents

METHODES NUMERIQUES APPLIQUEES AUX … relative au transfert de chaleur par conduction et en régime non stationnaire sera l’exemple à résoudre durant toute la première partie