-

Calcul numérique des

structures
Différences finies

Hicham Meftah

BTP4 - 2015/2016

H. Meftah - ENSA d’Agadir - [Link]@[Link] 1 / 47

 Introduction - Généralités

Le calcul numérique est devenu l’un des éléments essen-

tiels de tout programme de recherche et de développement.

H. Meftah - ENSA d’Agadir - [Link]@[Link] 2 / 47

 Introduction - Généralités

Avantages majeurs du calcul numérique :

Réduction significative du temps de conception et de déve-
loppement
Reproduction de conditions non-faisables expérimentalement

Fourniture d’informations plus détaillées et plus précises

Meilleure rentabilité et plus économique que l’expérience

Des expériences de référence sont cependant nécessaires

à la validation initiale des calculs.

H. Meftah - ENSA d’Agadir - [Link]@[Link] 3 / 47

 Introduction - Généralités

Traditionnellement, pour concevoir un fuselage d’avion de

ligne, Boeing effectuait une vingtaine d’expérience en souf-
flerie.
A présent, avec l’utilisation du calcul numérique, seules deux
ou trois campagnes expérimentales sont nécessaires.
Dans les souffleries, des informations globales sont obte-
nues : la traînée, la portance, la distribution de pression en
certains points.
Le calcul numérique permet d’obtenir en plus le champ de
vitesse et de pression autour de la carlingue (ce que l’ex-
périence peut faire mais pour un coût énorme).

H. Meftah - ENSA d’Agadir - [Link]@[Link] 4 / 47

 Introduction - Structures des équations

Le principe du calcul numérique est de résoudre des équa-

tions différentielles décrivant l’évolution de variables phy-
siques dans l’espace et le temps.
Les difficultés de résolution du problème sont généralement
liées à deux éléments principaux :
La complexité de la géométrie étudiée.
La complexité des couplages entre les différents phénomènes
physiques

H. Meftah - ENSA d’Agadir - [Link]@[Link] 5 / 47

 Introduction - Exemples

Exemple Numéro 1 : Simulation météorologique.

Géométrie sphérique “simple” à grande échelle mais très

grande (et donc coûteuse en mémoire et en temps).
De plus, les interactions entre l’atmosphère, les océans, les
surfaces émergées, les échanges thermiques, les réactions
chimiques à échelles variables, les rapports d’échelle, etc.
Font de ce problème l’un des plus complexes qui soit.

H. Meftah - ENSA d’Agadir - [Link]@[Link] 6 / 47

 Introduction - Exemples

Exemple Numéro 2 : Voilier.

Géométrie complexe et particulièrement “fine”

Interactions fluide/structure : solveurs multiples avec cou-
plages entre les deux.
Calcul de la prise au vent de la voile, de la déformation du
mat et des efforts sur la coque.
Calcul de l’aérodynamique, de l’équilibrage de la quille, etc.
Les échelles de longueur sont plus à “taille humaine” mais
les calculs peuvent être très complexes.

H. Meftah - ENSA d’Agadir - [Link]@[Link] 7 / 47

 Introduction - Application

Application : Oscillateur Harmonique

Nous allons prendre, comme premier exemple, un modèle

mécanique pour un oscillateur harmonique. Par exemple,
une masse attachée à l’extrémité d’un ressort.

F V

k raideur du ressort masse m

Exercice : Etablire l’équation de mouvement de l’oscillateur

H. Meftah - ENSA d’Agadir - [Link]@[Link] 8 / 47

 Introduction - Application

Application : Oscillateur Harmonique

L’équation du mouvement est donnée par :

d2 x
m = −kx
dt2
Et la solution de cette équation est de la forme :

x = Acos(ωt + ǫ)

Avec A et ǫ despconstantes dépendantes des conditions aux

limites et ω = (k/m), la pulsation de l’oscillateur.

H. Meftah - ENSA d’Agadir - [Link]@[Link] 9 / 47

 Introduction - Application

Application : Oscillateur Harmonique

Ce modèle peut être affiné et plus proche de la réalité, en

ajoutant deux forces
1. Une force de contrôle de la pulsation ωc imposée par l’opé-
rateur. Elle est de la forme Fcos(ωc t)
2. Une force de résistanceau mouvement −f (dx/dt)
La première correspond à une oscillation de l’axe de l’expé-
rience et la seconde aux frottements appliqués à la masse.
L’équation du mouvement (qui a une solution analytique)
s’écrit alors :

d2 x dx
m +f + kx = Fcos(ωc t)
dt2 dt

H. Meftah - ENSA d’Agadir - [Link]@[Link] 10 / 47

 Introduction - Application

Application : Oscillateur Harmonique

Choisissons à présent un cas où la force de résistance n’est

plus liniéaire avec la vitesse mais à une puissance arbitraire
de la vitesse, soit −f |dx/dt|α
L’équation du mouvement est alors de la forme

d2 x dx
m 2
+ f | |α−1 + kx = Fcos(ωc t)
dt dt

Il n’existe pas de solution analytique, une résolution numé-

rique est alors nécessaire.

H. Meftah - ENSA d’Agadir - [Link]@[Link] 11 / 47

 Techniques de Discrétisation D.F.

Avec une résolution de type Différences Finies (DF), le do-

maine du problème (c.a.d l’espace de résolution) est dis-
crétisé. C’est à dire divisé en points.
Les fonctions inconnues seront résolues sur ces points dis-
crets.
f (x)

f (x)

H. Meftah - ENSA d’Agadir - [Link]@[Link] 12 / 47

 Techniques de Discrétisation D.F.

Les differentes dérivées sont approximées par un développ-

ment de Taylor adapté à la précision désirée.

La précision de la méthode ou du schéma numérique va dé-

péndre du nombre de points retenus pour approximer une
dérivée donnée.

L’utilisation des différences finies est simple mais une atten-

tion particulière doit être portée sur la précision des sché-
mas et la nature du maillage.

H. Meftah - ENSA d’Agadir - [Link]@[Link] 13 / 47

 Calcul de la dérivée 1ère

Soit une fonction continue f (x) ou discrète fi dont les valeurs

sont connues en tout point xi discret du maillage. Nous de-
vons déterminer la dérivée de cette fonction en tout point
de ce maillage.
Ecrivons l’expension en série de Taylor au voisinage de xi
de la fonction f (x)

(x − xi )2 ∂ 2 f
   
∂f
f (x) = f (xi ) + (x − xi ) +
∂x i 2! ∂x2 i
(x − xi )n ∂ n f
 
+ ... + +H
n! ∂xn i

H représente les termes d’ordre élevé.

H. Meftah - ENSA d’Agadir - [Link]@[Link] 14 / 47

 Calcul de la dérivée 1ère

Approchons la fonction en x = xi+1 (Equation 1)

(xi+1 − xi )2 ∂ 2 f
   
∂f
f (xi+1 ) = f (xi ) + (xi+1 − xi ) + + ...
∂x i 2! ∂x2 i
(1)
Approchons la fonction en x = xi−1 (Equation 2)

(xi−1 − xi )2 ∂ 2 f
   
∂f
f (xi−1 ) = f (xi ) + (xi−1 − xi ) + + ...
∂x i 2! ∂x2 i
(2)
Il est possible d’estimer la dérivée de f (x) en xi de trois
∂f
manières différentes : en isolant ∂x dans l’équation 1, ou
bien dans l’équation 2 ou encore dans l’équation résulante
de l’addition de 1 et 2.

H. Meftah - ENSA d’Agadir - [Link]@[Link] 15 / 47

 Calcul de la dérivée 1ère

Possibilité 1 (de l’équation 1)

(xi+1 − xi ) ∂ 2 f
   
∂f fi+1 − fi
= − + ...
∂x i xi+1 − xi 2! ∂x2 i

Possibilité 2 (de l’équation 2)

(xi − xi−1 ) ∂ 2 f
   
∂f fi − fi−1
= + + ...
∂x i xi − xi−1 2! ∂x2 i

Possibilité 3 (de l’équation 1+2)

fi+1 − fi−1 (xi+1 − xi )2 − (xi − xi−1 )2 ∂ 2 f

   
∂f
= − +...
∂x i xi+1 − xi−1 2(xi+1 − xi−1 ) ∂x2 i

H. Meftah - ENSA d’Agadir - [Link]@[Link] 16 / 47

 Calcul de la dérivée 1ère

Toutes les expressions de la page précédentes sont exactes

si tous les termes d’ordre élevé sont pris en compte. Ce-
pendant, d’un point de vue pratique, il est nécessaire de
tronquer les expressions au niveau de précision désiré.
Au plus simple, après troncature des termes d’ordre 2 et
plus, nous obtenons :
 
∂f fi+1 − fi
= Dérivée avant (ordre 1)
∂x i xi+1 − xi
 
∂f fi − fi−1
= Dérivée arrière (ordre 1)
∂x i xi − xi−1
 
∂f fi+1 − fi−1
= Dérivée centrée (ordre 2)
∂x i xi+1 − xi−1

H. Meftah - ENSA d’Agadir - [Link]@[Link] 17 / 47

 Calcul de la dérivée 1ère

Schéma du calcul de la dérivée.

Dérivée entrée
Dérivée arrière Dérivée avant

fi

fi−1 fi+1

i−1 i i+1 x

Cal ul de la dérivée en i

H. Meftah - ENSA d’Agadir - [Link]@[Link] 18 / 47

 Calcul de la dérivée 2nde

Considérons un maillage régulier tel que xi+1 − xi = ∆x.

Définissons l, un nombre entier tel que xi+l = xi + l ∗ ∆x.

Ecrivons l’expension en série de Taylor au voisinage de xi

de la fonction f (xi+l )

(l∆x)2 ∂ 2 f
   
∂f
f (xi+l ) = f (xi ) + l∆x +
∂x i 2! ∂x2 i
(l∆x)n ∂ n f
 
+ ... + +H
n! ∂xn i

H. Meftah - ENSA d’Agadir - [Link]@[Link] 19 / 47

 Calcul de la dérivée 2nde

On approche la fonction en x = xi+1 :

(∆x)2 ∂ 2 f
   
∂f
f (xi+1 ) = f (xi ) + ∆x + + ... (3)
∂x i 2! ∂x2 i

et en x = xi+2 :

(2∆x)2 ∂2f
   
∂f
f (xi+2 ) = f (xi ) + 2∆x + + ... (4)
∂x i 2! ∂x2 i

Ce qui nous donne (Eq. 4 - 2*Eq.3)

 2 
∂ f f (xi+2 ) − 2f (xi+1 ) + f (xi )
2
=
∂x i (∆x)2

H. Meftah - ENSA d’Agadir - [Link]@[Link] 20 / 47

 Calcul de la dérivée 2nde

A nouveau, il est possible de définir des dérivées centrées

ou décentrées.
 2 
∂ f fi+2 − 2fi+1 + fi
2
= Dérivée avant (ordre 1)
∂x i (∆x)2

∂2f
 
fi − 2fi−1 + fi−2
= Dérivée arrière (ordre 1)
∂x2 i (∆x)2
∂2f
 
fi+1 − 2fi + fi−1
= Dérivée arrière (ordre 2)
∂x2 i (∆x)2
La dernière définition est très utilisée.

H. Meftah - ENSA d’Agadir - [Link]@[Link] 21 / 47

 Equations aux dérivées partielles - EDP

Définition : Une équation aux dérivées partielles (ou EDP)

est une relation entre une fonction inconnues de plusieurs
variables et ses dérivées partielles.

Définissons des variables indépendantes notées x, y, z, . . .

et les variables dépendantes définies par u, v, w, ...

Une EDP s’écrit sous la forme générale suivante :

∂u ∂ 2 u ∂ 2 u ∂ 2 u
F(x, y, z, u, , , , , ...) = 0
∂x ∂x2 ∂y2 ∂x∂y

H. Meftah - ENSA d’Agadir - [Link]@[Link] 22 / 47

 EDP - Ordre d’une EDP

L’ordre d’une EDP est défini par l’ordre maximal des déri-
vées présentes dans l’équation

∂u
Ordre 1 : ∂x
− c ∂u
∂y
=0
 3
∂2u ∂2u
Ordre 2 : ∂x2
+ ∂y2
=0

∂2u 4
Ordre 4 : ∂x2
+ c ∂∂xu4 = 0

Si des EDP interdépendantes sont impliquées, l’ordre n’est

calculé qu’après avoir combiné les équations en une seule.

H. Meftah - ENSA d’Agadir - [Link]@[Link] 23 / 47

 EDP - Linéarité

Une EDP est dite linéaire si l’ordre de la puissance appli-

quée à ses dérivées partielles est 1 :

∂u ∂u
∂x
− ∂y
= 0 est linéaire.

∂u 2 ∂2u


∂x
+ ∂y2
= 0 est non-linéaire.

H. Meftah - ENSA d’Agadir - [Link]@[Link] 24 / 47

 EDP - Linéarité

Considérons
∂u ∂u
A +B =C
∂x ∂y

Si on a A(x, y) et B(x, y) où A et B constants, l’équation est

linéaire.

Si on a A(x, y, u) et B(x, y, u), l’équation est quasi-linéaire.

Si on a A(x, y, u, ∂u ∂u ∂u ∂u
∂x , ∂y ) et B(x, y, u, ∂x , ∂y ), l’équation est non-
linéaire.
De manière générale, pour une EDP d’ordre n, si ses co-
efficients dépendent de dérivées d’ordre n, l’EDP est non-
linéaire. Si les coefficients dépendent de dérivées d’ordre
< n, l’équation est quasi-linéaire.

H. Meftah - ENSA d’Agadir - [Link]@[Link] 25 / 47

 EDP - Conditions aux limites

En général, il peut y avoir une infinité de solutions à une

EDP. Pour choisir celle qui représente la solution à un pro-
blème physique, il est nécessaire de préciser les conditions
spécifiques au problème : les conditions aux limites et les
conditions initiales.
u=0

∂u
u = u1 ∂n =0

u=0

H. Meftah - ENSA d’Agadir - [Link]@[Link] 26 / 47

 EDP - Conditions aux limites

Considérons le volume de résolution V. La frontière de cette

région est notée S . Une EDP définie par l’opérateur L tel
que L(u) = 0 peut être définie avec u la fonction inconnue
du problème. Il y a trois types de conditions aux frontières :

Les conditions de Dirichlet : u = g sur x ∈ S.

∂u
Les conditions de Neumann : ∂n = h sur x ∈ S.

Les conditions mixtes : ku + µ ∂u

∂n = k sur x ∈ S.

g, h, k, λ et µ sont des fonctions imposées sur S .

H. Meftah - ENSA d’Agadir - [Link]@[Link] 27 / 47

 EDP - Conditions initiales

Condition initiale : Il s’agit en fait de la condition limite lié à

la variable temporelle. La condition initiale donne la valeur
de la fonction inconnue à l’instant où e processus décrit par
l’EDP débute. Il peut s’agir d’une combinaison de la variable
indépendante et de ses dérivées temporelles.
Un problème est bien posé si :

une solution existe

la solution est unique

la solution dépend continument des données

de légères modifications des données correspondent à de

légères modifications de la solution.

H. Meftah - ENSA d’Agadir - [Link]@[Link] 28 / 47

 Méthode différence finie explicite

Equation de la chaleur

La loi de Fourier pour la conduction de la chaleur s’écrit :

∂T ∂2T
=k 2
∂t ∂x
Les trois points suivants doivent être pris en compte :

Evolution à la fois temporelle et spatiale.

La limite temporelle est "ouverte".

Attention aux problèmes de stabilité.

H. Meftah - ENSA d’Agadir - [Link]@[Link] 29 / 47

 Méthode différence finie explicite

Exemple en 2D + le temps
instant initial t = 0 t > 0 la température diuse dans le domaine

H. Meftah - ENSA d’Agadir - [Link]@[Link] 30 / 47

 Méthode différence finie explicite

Exemple en 1D + le temps
t, dire tion temporelle

n+1

n−1

t=0
i−1 i i+1

x dire tion spatiale

H. Meftah - ENSA d’Agadir - [Link]@[Link] 31 / 47

 Méthode différence finie explicite

Les dérivées sont déterminées par un développement de

Taylor et le domaine de résolution est divisé en pas constants
∆t pour le temps et ∆x pour l’espace.
Til = T(t = l ∗ ∆t, x = i ∗ ∆x)
Dérivées temporelles et spatiales.
n − 2T n + T n
∂T T n+1 − Tin ∂2T Ti−1 i i+1
= i et =
∂t ∆t ∂x2 ∆x2

n+1

i−1 i i+1

H. Meftah - ENSA d’Agadir - [Link]@[Link] 32 / 47

 Méthode différence finie explicite

Après discrétisation, l’équation de la chaleur donne donc :

Tin+1 − Tin T n − 2Tin + Ti+1

n
= k i−1
∆t ∆x2
et finalement :

Tin+1 = r Ti−1
n n
+ (1 − 2r) Tin


+ Ti+1

avec r = k∆t/∆x2

H. Meftah - ENSA d’Agadir - [Link]@[Link] 33 / 47

 Méthode différence finie explicite

Consistance et convergence
Une méthode est convergente si, lorsque les pas de temps
∆t et d’espace ∆x tendent vers 0, la solution approchée
tend vers la solution exacte

(∆t −→ 0, ∆x −→ 0) =⇒ Ti = T(xi )

Une condition nécessaire à la convergence est la consis-

tance des équations. Lorsque les pas de temps et d’espace
tendent vers 0, les équations approchées doivent tendre
vers l’équation analytique.

H. Meftah - ENSA d’Agadir - [Link]@[Link] 34 / 47

 Méthode différence finie explicite

Consistance et convergence

On peut montrer que si (∆t → 0, ∆x → 0) et si r ≤ 1/2 ,

la solution numérique tend vers la solution analytique. La
méthode est donc consistante.

H. Meftah - ENSA d’Agadir - [Link]@[Link] 35 / 47

 Méthode différence finie explicite

stabilité

Un système est convergent si il est à la fois consistant et

stable.

La stabilité nous indique si l’erreur augmente ou non au

cours du calcul. Si elle est faible au départ puis dépasse
les limites raisonnables, le système est instable.

Au contraire, si l’erreur est maintenue dans un domaine suf-

fisamment petit, le système est stable.

La méthode explicite présentée est inconditionellement stable

si r ≤ 0.5.

H. Meftah - ENSA d’Agadir - [Link]@[Link] 36 / 47

 Méthode différence finie implicite

On utilise dans l’évaluation de la dérivée seconde des termes

inconnus au niveau n + 1

Le schéma le plus simple s’écrit :

Tin+1 − Tin T n+1 − 2Tin+1 + Ti+1

n+1
= k i−1
∆t ∆x2

qui peut se réécrire sous la forme :

n+1
−rTi−1 + (1 + 2r)Tin+1 − rTi+1
n+1
= Tin

avec r = k∆t/∆x2
Ce schéma est inconditionnellement stable.

H. Meftah - ENSA d’Agadir - [Link]@[Link] 37 / 47

 Méthode différence finie implicite

Il nous faut résoudre la matrice suivante :

 
1 + 2r −r  n+1
T1n
  
T1
 −r 1 + 2r −r 
   .   . 
 .. .. ..     

 . . .   .
∗
 
= . 

 .. .. ..   .   . 
 . . .     
   .   . 
 −r 1 + 2r −r 
Tmn+1 Tmn
−r 1 + 2r

T1n+1 et Tmn+1 sont connus (conditions limites).

Système tri-diagonal qui peut être résolu grâce à l’algo-
rithme de Thomas.

H. Meftah - ENSA d’Agadir - [Link]@[Link] 38 / 47

 Méthode différence finie implicite
Algorithme de Thomas

un Système tri-diagonal s’écrit sous la forme suivante :

ai xi−1 + bi xi + ci xi+1 = di

avec a1 = 0 et cn = 0, sous forme matricielle, on a :

 
b1 c1    
 a2 x1 d1
b2 c2 
   x2   d2 
 .. .. ..     

 . . .  
∗ .  
= . 

 .. .. ..   .   . 
 . ..     
   .   . 
 an−1 bn−1 cn−1 
xn dn
an bn

H. Meftah - ENSA d’Agadir - [Link]@[Link] 39 / 47

 Méthode différence finie implicite

Algorithme de Thomas

La première étape consiste à modifier les coefficients de la

manière suivante :

c1
; i=1


b1

c′i =
ci
; i = 2, 3, . . . , n − 1


bi − c′i−1 ai


et 
d1
; i=1


b1

di′ = ′ a
di − di−1 i
; i = 2, 3, . . . , n


 ′
bi − ci−1 ai

H. Meftah - ENSA d’Agadir - [Link]@[Link] 40 / 47

 Méthode différence finie implicite

Algorithme de Thomas

La deuxième étape consiste à la résolution de la manière

suivante :
(
xn = dn′
xi = di′ − c′i xi+1 ; i = n − 1, n − 2, . . . , 1

H. Meftah - ENSA d’Agadir - [Link]@[Link] 41 / 47

 Méthode de Crank-Nicolson

Les dérivées spatiales sont à moitié évaluées au temps n et

à moitié au temps n + 1 :

Tin+1 − Tin
 
k 1 n n n 1 n+1 n+1 n+1
= (T − 2Ti + Ti−1 ) + (Ti+1 − 2Ti + Ti−1 )
∆t ∆x2 2 i+1 2

Ce schéma est inconditionnellement stable.

H. Meftah - ENSA d’Agadir - [Link]@[Link] 42 / 47

 Méthode de Crank-Nicolson

Généralisation de Crank-Nicolson

Une généralisation du schéma de Crank-Nicolson s’écrit

sous la forme suivante :

Tin+1 − Tin k  n
= (1 − θ)(Ti+1 − 2Tin + Ti−1
n
)
∆t ∆x2
n+1
− 2Tin+1 + Ti−1
n+1

+ θ(Ti+1 )

avec θ une constante comprise entre 0 et 1.

H. Meftah - ENSA d’Agadir - [Link]@[Link] 43 / 47

 Méthode de Crank-Nicolson

Généralisation de Crank-Nicolson

On peut constater que si θ = 0, on a un simple schéma

explicite ; si θ = 1/2, on retrouve Crank-Nicolson et enfin si
θ = 1 on a affaire à un schéma implicite

Ce schéma est stable si :

∆t 1 1

 0 ≤ k ∆x2 ≤ 2−4θ pour 0 ≤ θ < 2

1
inconditionnellement stable pour ≤θ≤1

2

H. Meftah - ENSA d’Agadir - [Link]@[Link] 44 / 47

 Méthode de Crank-Nicolson

Généralisation de Crank-Nicolson

On peut constater que si θ = 0, on a un simple schéma

explicite ; si θ = 1/2, on retrouve Crank-Nicolson et enfin si
θ = 1 on a affaire à un schéma implicite

Ce schéma est stable si :

∆t 1 1

 0 ≤ k ∆x2 ≤ 2−4θ pour 0 ≤ θ ≤ 2

1
inconditionnellement stable pour ≤θ≤1

2

H. Meftah - ENSA d’Agadir - [Link]@[Link] 45 / 47

 Méthode ADI

ADI signifie Alternating Direction Implicit (Peaceman, Rach-

ford et Douglas 1955).
Cette méthode rentre dans la catégorie des méthode du
pas fractionnaire ou méthode de splitting
C’est un schéma implicite à deux pas. Il est de type prédiction-
correction :

 n+1/2
 n+1/2 n+1/2 n+1/2

n Ti+1,j −2Ti,j +Ti−1,j n
Ti,j+1 n +T n
Ti,j −Ti,j −2Ti,j i,j−1
 Pas 1 : = k +


∆t/2 ∆x2 ∆y2
n+1 n+1/2
 n+1/2 n+1/2 n+1/2 n+1 n+1 n+1

Ti,j −Ti,j Ti+1,j −2Ti,j +Ti−1,j Ti,j+1 −2Ti,j +Ti,j−1
 Pas 2 : = k +


∆t/2 ∆x2 ∆y2

H. Meftah - ENSA d’Agadir - [Link]@[Link] 46 / 47

 Méthode ADI

On avance d’une moitié de pas à chaque fois. A chaque

demi pas, les terme implicites (n + 1/2 au pas 1 et n + 1 au
pas 2) ne portent que sur une des deux dimensions (ici x
au pas 1 et y au pas 2).
Avec la méthode ADI, on inverse deux fois une matrice tri-
diagonale.

H. Meftah - ENSA d’Agadir - [Link]@[Link] 47 / 47