M2 MIR – Microsystèmes Multiphysiques

Utilisation de la méthode des éléments finis pour la simulation

numérique des équations multi-physiques

Vincent Laude & Mahmoud Addouche

Institut FEMTO-ST, Université de Bourgogne Franche-Comté et CNRS,

15 B avenue des Montboucons, F-25030 Besançon Cedex, France
[Link]@[Link]
[Link]

v2.2, septembre 2017

Contenu du cours

1 Principes, 2D, problèmes scalaires (électrostatique, acoustique)

2 Modèles 3D, éléments finis vectoriels (élasticité)

3 Couplages multi-physiques

Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 2 / 34
Organisation du cours
Objectifs :
Formuler des modèles simples mais multi-physiques
Comprendre et utiliser la méthode des éléments finis
Vincent Laude :
Modélisation, maillage, choix des méthodes (5h CM)
Acoustique, modes de vibration élastique (4h TP)
Électrostatique (4h TP)
Couplage déformation mécanique–électrostatique (4h TP)
Couplage acoustique–structure (4h TP)
Mahmoud Addouche :
Performance énergétique d’une façade d’habitation (4h TP)
Dissipateurs pour CPU. Couplage Thermoélectrique; Effet Joule (4h
TP)
Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 3 / 34
 Principes

Équations aux dérivées partielles

Toute équation du type

Lu = f
avec u(x, y ) une fonction définie sur un domaine Ω de l’espace physique,
L un opérateur différentiel contenant x, y , u, ∂u ∂u
∂x , ∂y ... et des fonctions de
ces quantités.
Exemples :
Équation de Laplace : −∇ · (∇u) = 0
Équation de Poisson : −∇ · (∇φ) = ρ
Équation de Helmholtz : −∇ · (c∇u) − k 2 u = 0

Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 4 / 34
 Principes

Méthode de Galerkin (for the dummies)

Lu = f est définie sur un certain domaine Ω. On choisit de représenter u
sur une base de fonctions wj , soit
n
X
u(x, y ) = aj wj (x, y )
j=1

où aj sont des coefficients (réels ou complexes).

On projette l’équation différentielle sur les fonctions wi
Z Z
wi Lu = wi f , ∀i = 1...n
Ω Ω

On obtient une équation linéaire

Aa = f

La solution formelle est a = A−1 f.

R R
avec Aij = Ω wi Lwj et fi = Ω wi f .
Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 5 / 34
 Principes

Il ne faut pas avoir peur des équations linéaires !

Il existe beaucoup de programmes informatiques pour résoudre les
équations linéaires. Exemples : Matlab, Comsol, etc.
Inversion de matrice
A−1 A = I
Problème linéaire : trouver x tel que

Ax = f

Problème aux valeurs propres : trouver (λ, x) tels que

Ax = λx

Problème aux valeurs propres généralisé : trouver (λ, x) tels que

Ax = λBx

Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 6 / 34
 Principes

Maillage et régions

On divise le domaine Ω en régions dans lesquelles les coefficients de

l’équation différentielle sont homogènes, ou dans lesquelles on sait
que la solution doit être continue.
Chacune des régions est maillée. La qualité du maillage est très
importante pour obtenir des solutions de qualité. Un maillage trop
serré conduit à des temps de calcul très longs ; un maillage trop lâche
conduit à des erreurs importantes.

Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 7 / 34
 Principes

Conditions aux limites

Le comportement de la solution aux frontières est décrite par des

conditions aux limites (CL). Il en existe 2 types principaux.
1 Dirichlet, ou conditions aux limites essentielles :

u = 0, ou u = u0 , sur σ
2 Neumann, ou conditions aux limites naturelles :
∂u
=0
∂n
où n est la normale à la frontière σ.
Les logiciels permettent aussi parfois des conditions de type mixte
(Neumann + Dirichlet), périodicité, symétrie, rayonnement, absorption,
etc.

Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 8 / 34
 Principes

Espace d’éléments finis

Étant donné un domaine Ω et son maillage Th, on choisit de décrire

la solution par un nombre fini de degrés de liberté (dof ), par exemple
les valeurs nodales uje pour les j nœuds de l’élément e.
L’espace d’éléments finis Wh est l’ensemble des fonctions
représentables par ce choix. Propriété importante en pratique : c’est
un espace fonctionnel de dimension finie.
Représentation
P dans un élément fini
u (x, y ) = j Nj (x, y )uje où Nje (x, y ) sont les fonctions de base.
e e

Pour les éléments de Lagrange, Pn , il y a continuité de

u(x, y ) = e u e (x, y ) entre les éléments.
P

Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 9 / 34
 Principes

Formulation faible

Supposons que nous avons choisi un espace d’éléments finis Wh (pour

un domaine Ω et son maillage Th).
On remplace le problèmeRinitial Lu =
R f par une approximation :
Trouver u ∈ Wh tel que Ω w Lu = Ω wf pour toutes les fonctions
test w ∈ Wh.
u et w sont complètement déterminés par les valeurs aux nœuds
U = {uje } et W = {wje }.
Il existe donc une matrice K et un vecteur B tels que

WT K U = WT B, ∀W

Et donc l’équation linéaire : K U = B

Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 10 / 34
 Principes

Résolution pratique par le logiciel FEM ?

1 Dessiner (ou définir) le domaine et les régions

2 Mailler le domaine (discrétisation du domaine)
3 Choisir le type des éléments finis
4 Définir les conditions aux limites
5 Définir les coefficients des équations physiques, les constantes utiles
6 Choisir le solveur (et ses options si nécessaire)
7 Résoudre le problème
8 Post-traiter les résultats

Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 11 / 34
 Principes

Électrostatique 2D

Equation de Poisson :
−∇ · (∇φ) = ρ
avec  la constante diélectrique et ρ une densité de charge.

Exemple pratique : 2 cylindres métalliques portant des charges, séparés

par une courte distance
Définition du domaine ?
Conditions aux limites internes (potentiel électrique, charge de
surface) ?
Conditions aux limites externes ?
Utiliser le maillage par défaut et analyser le résultat
Comment améliorer le maillage ?

Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 12 / 34
 Principes

Cas 1D : Comment le logiciel FEM résout le problème ?

1 Projection de l’équation différentielle sur les fonctions test

Z 1   Z 1
∂ ∂φ(x)
− dxψ(x) (x) = dxψ(x)ρ(x)
0 ∂x ∂x 0

2 Intégration par parties

1
∂φ 1
Z   Z 1
∂ψ ∂φ
dx  − ψ = dxψρ
0 ∂x ∂x ∂x 0 0
h i1
3 Application des conditions aux limites : ψ ∂φ
∂x 0 est connu ; par
exemple :
Dirichlet : ψ(0) = 0
Neumann : ∂φ(1)
∂x = 0

Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 13 / 34
 Principes

Cas 2D : Comment le logiciel FEM résout le problème ?

1 Projection de l’équation différentielle sur les fonctions test

Z Z
− dr ψ(r )∇ · ((r )∇φ(r )) = dr ψ(r )ρ(r )
Ω Ω

2 Théorème de Gauss (vrai en 2D et 3D)

Z Z Z
∂φ(r )
dr ∇ψ · (∇φ) − dnψ = dr ψρ
Ω σ ∂n Ω

ψ ∂φ(r )
R
3 Application des conditions aux limites : σ ∂n est connu ; par
exemple :
Dirichlet : ψ(r ) = 0 sur une partie de σ
Neumann :  ∂φ(r )
∂n = 0 sur le reste de σ

Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 14 / 34
 Principes

Acoustique (harmonique)

Equation acoustique en régime harmonique, avec source :

∇ · (−1/ρ0 ∇p) − ω 2 /(ρ0 cS2 )p = f

avec cS la vitesse acoustique (m/s), ρ0 la densité (kg/m3 ), f un terme

source (1/s2 )

Exemple pratique : haut-parleur dans une pièce ou dehors

Définition du domaine
Conditions aux limites (pression, accélération normale, soft and hard
boundary, radiation) ?
Appliquer une force : Dirichlet ou Neumann ? Quelle réponse ?
Qu’est-ce qu’un bon maillage pour une onde ?
Changer les éléments finis ; la solution change t-elle ?

Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 15 / 34
 Principes

Acoustique (valeurs propres)

Equation acoustique en régime harmonique, avec nombre d’onde axial :

∇ · (−1/ρ0 ∇p) + (kz2 /ρ0 − ω 2 /(ρ0 cS2 ))p = 0

avec cS la vitesse acoustique (m/s), ρ0 la densité (kg/m3 ) et kz un

nombre d’onde axial (1/m)

Exemple pratique : tuyau d’orgue

Définition du domaine
Conditions aux limites (soft and hard boundary) ?
Trouver les modes propres. Que représentent-ils ?
Varier kz : tuyau d’orgue ?

Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 16 / 34
 Principes

Problèmes aux valeurs propres (synthèse)

Soit une matrice carrée M de dimension n, à valeurs réelles ou

complexes. Une équation aux valeursP propres λ et aux vecteurs
propres u est de la forme Mu = λu ( ni=1 Mij uj = λui ). Les valeurs
propres sont racines du polynôme caractéristique |M − λI | = 0.
Il y a exactement n valeurs propres et au plus n vecteurs propres (a
priori complexes). Les vecteurs propres sont non nuls et donc
normalisables (on ne les connaı̂t qu’à une constante près).
Si M est réelle symétrique (ou Hermitienne), les valeurs propres sont
réelles et les vecteurs propres sont orthogonaux.
Problème aux valeurs propres généralisé : Au = λBu
En pratique, il existe des algorithmes très performants pour
déterminer valeurs et vecteurs propres.

Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 17 / 34
 Modèles 3D

Élasticité 3D : Déformations

Soit un point x de coordonnées (x1 , x2 , x3 ).

ui (x + dx) = ui (x) + ∂u
∂xj dxj au premier ordre.
i ∂ui
∂xj est le gradient des
déplacements.
On sépare ce gradient en partie symétrique (le tenseur des
déformations Sij ) et antisymétrique
 selon  
∂ui 1 ∂ui ∂uj 1 ∂ui ∂uj
∂xj = Sij + ASij , Sij = 2 ∂xj + ∂xi et ASij = 2 ∂xj − ∂xi
Seule la partie symétrique du gradient des déplacements mesure une
déformation locale (la partie antisymétrique mesure les rotations
locales).
La dilatation (variation locale du volume) est
S = S11 + S22 + S33 = ∇ · u.
Les termes S11 , S22 et S33 correspondent à des mouvements
longitudinaux, les termes Sij , i 6= j, à des mouvements de cisaillement.

Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 18 / 34
 Modèles 3D

Élasticité 3D : Contraintes

Trois forces indépendantes s’exercent sur

une surface : une contrainte de traction-
compression et deux contraintes de ci-
saillement.
Sur la face orthogonale à x1 d’un cube, la force par unité de surface
ou tension mécanique est T11 + T21 + T31 . Tij est un tenseur de rang
2 symétrique, le tenseur des contraintes. Pour une surface de normale
n, la tension mécanique est le vecteur Tij nj .
L’équation de la dynamique s’écrit (avec fi les forces internes)
∂Tij ∂ 2 ui
∂xj + fi = ρ ∂t 2

Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 19 / 34
 Modèles 3D

Élasticité 3D : Loi de Hooke

L’expérience montre que le comportement élastique de la plupart des

solides dans le cas de petites déformations obéit à la loi de Hooke :

Tij = cijkl Skl

c’est-à-dire que les contraintes sont une fonction linéaire des

déformations.
cijkl est le tenseur des rigidités, de rang 4. Il a a priori 34 = 81
composantes. Mais la symétrie de Tij et Skl impose que cjikl = cijkl et
cijkl = cijlk . Il n’y a donc que 36 composantes indépendantes au plus.
Avec la symétrie cijkl = cklij on tombe à 21 composantes
indépendantes au plus.

Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 20 / 34
 Modèles 3D

Élasticité 3D : Notation contractée (ou matricielle)

(11) −→ 1; (22) −→ 2; (33) −→ 3

(23) = (32) −→ 4; (13) = (31) −→ 5; (12) = (21) −→ 6
TI = Tij ; cIJ = cijkl ; TI = cIJ SJ
S1 = S11 ; S2 = S22 ; S3 = S33 ; S4 = 2S23 ; S5 = 2S13 ; S6 = 2S12

Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 21 / 34
 Modèles 3D

Maillage 3D

Comme en 2D on peut mailler tout domaine connexe par des

triangles...
En 3D on peut utiliser des tétraèdres

Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 22 / 34
 Modèles 3D

Éléments finis vectoriels

On considère le vecteur (u1 , u2 , u3 ).

Étant donné un domaine Ω et son maillage Th, on décrit la solution
par un nombre fini de degrés de liberté (dof ), les valeurs nodales uije
pour les j nœuds de l’élément e.
L’espace d’éléments finis Wh est l’ensemble (fini) des fonctions
représentables par ce choix.
Représentation
P dans un élément fini
ui (x, y ) = j Nj (x, y )uije où Nje (x, y ) sont les fonctions de base
e e

(comme dans le cas scalaire).

Pour les éléments
P e de Lagrange, Pn , il y a continuité de
ui (x, y ) = e ui (x, y ) entre les éléments.

Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 23 / 34
 Modèles 3D

Formulation faible de l’élasticité 3D (1)

1 Équation de l’élastodynamique

∂2u
−∇T + ρ =f
∂t 2
2 Projection sur les fonctions test v

∂2u
Z Z Z
− v · ∇T + v·ρ 2 = v·f
Ω Ω ∂t Ω

3 Théorème de Gauss
∂2u
Z Z Z Z
∇vT − v · Tn + v·ρ 2 = v·f
Ω σ Ω ∂t Ω

Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 24 / 34
 Modèles 3D

Formulation faible de l’élasticité 3D (2)

1 Loi de Hooke : TI = cIJ SJ avec S = ∇u

∂2u
Z Z Z Z
S(v)I cIJ S(u)J − v · Tn + v·ρ 2 = v·f
Ω σ Ω ∂t Ω
R
2 Conditions aux limites : σ v · Tn est connu
Dirichlet : v = 0 sur une partie de σ (blocage)
Neumann : Tn = 0 sur le reste de σ (traction nulle)
3 Cas monochromatique (ou excitation harmonique)
Z Z Z Z
S(v)I cIJ S(u)J − v · Tn − ω 2 v · ρu = v·f
Ω σ Ω Ω

Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 25 / 34
 Modèles 3D

Déformation d’une poutre sous l’effet d’une force

On se donne une poutre en silicium, ancrée d’un côté.

Elle fléchit sous sa propre masse, en statique. De combien ?
Visualiser la déformation.

Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 26 / 34
 Modèles 3D

Modes propres de vibration d’une poutre

On considère de nouveau la poutre en silicium, ancrée d’un côté, des

deux, ou simplement “flottant” dans l’espace.
On néglige la gravité, mais on veut estimer son comportement
dynamique par ses modes propres de déformation.
Obtenir ces modes propres.
Que change la condition d’ancrage ou de traction nulle ?

Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 27 / 34
 Couplages multi-physiques

Couplage des équations physiques

Beaucoup de matériaux induisent un couplage de grandeurs physiques

autrement considérées séparément
Électro-optique, magnéto-optique
Piézomagnétique, piézoélectrique
Des relations constitutives relient les grandeurs physiques (constantes
matériaux)
Il suffit de considérer un espace d’éléments finis vectoriel dont les
inconnues sont les différentes grandeurs physiques.

Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 28 / 34
 Couplages multi-physiques

Exemple : piézoélectricité

Effet piézoélectrique direct : sous l’action d’une contrainte ou d’une

déformation, une polarisation électrique apparaı̂t
Effet piézoélectrique inverse : un champ électrique appliqué provoque
une déformation de la maille cristalline ou une contrainte. L’effet
piézoélectrique n’apparaı̂t que pour les structures cristallines non
centro-symétriques.
Équations couplées (ekij : constantes piézoélectriques)

Tij = cijkl Skl − ekij Ek ; Di = eikl Skl + ij Ej

∂φ
Le champ électrique dérive d’un potentiel : Ei = − ∂xi
En notation contractée

TI = cIJ SJ − ekI φ,k ; Di = eiJ SJ + ij φ,j

Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 29 / 34
 Couplages multi-physiques

Couplage par les conditions aux limites

Dans de nombreux cas, le système considéré est composé de

différentes régions, chacune gouvernée par un ensemble de lois
physiques
Couplage fluide - structure
Couplage acoustique - élastique (transducteurs acoustiques)
Couplage électrostatique - piézoélectrique
A l’interface entre les régions, les conditions aux limites permettent
de coupler les équations d’éléments finis :
Milieu 1 : A1 u1 + (intégrale de frontière sur u2 ) = f1
Milieu 2 : A2 u2 + (intégrale de frontière sur u1 ) = f2
Système complet
    
A1 C12 u1 f1
=
C21 A2 u2 f2

Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 30 / 34
 Couplages multi-physiques

Exemple : couplage acoustique - élastique

Comment une membrane qui vibre communique t’elle son mouvement

à l’air (ou l’eau) ? Et le problème inverse (le micro, ou le
transducteur) ?
Modèle élastique pour la membrane
Modèle acoustique pour l’air ou l’eau
A l’interface :
L’accélération normale de la membrane crée une variation de pression

∂ 2 un 1 ∂p
2
=
∂t ρF ∂n

La pression du fluide exerce une force (de surface) sur la membrane

Tij nj = −pni

Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 31 / 34
 Couplages multi-physiques

Couplage indirect de modèles physiques

On peut aussi modéliser des situations variées induisant des couplages

multiphysiques
Un système produit une action à distance, par exemple une force, sur
un autre système...
Une grandeur physique produit une déformation du système, le
maillage est donc déformé et change la solution pour une autre
grandeur qui était pourtant découplée...
Il n’y a pas de règle simple pour savoir comment résoudre le
problème : il faut improviser un modèle de couplage !

Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 32 / 34
 Couplages multi-physiques

Exemple : couplage électrostatique - déformations

Soit une membrane diélectrique métallisée sur une face, séparée par
un espace vide (un gap) d’une contre-électrode sur un substrat
La membrane est ancrée à un cadre support, ce qui définit les
conditions aux limites mécaniques
On résout le problème électrostatique, qui dépend de la distance entre
membrane et substrat
La force électrostatique déforme la membrane... ce qui change le
champ électrique... et ainsi de suite !
(La solution doit donc être obtenue de façon itérative)

Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 33 / 34
 Couplages multi-physiques

Quelques théorèmes intégraux utiles

H
V un volume, S une surface. désigne une intégrale sur une surface ou
un contour fermés.
Théorème du gradient
Z I
dV ∇f = nfdS
V S

Théorème de Gauss
Z I
dV ∇ · f = dSn · f
V S

Théorème du rotationnel
Z I
dV ∇ × f = dSn × f
V S

Vincent Laude (UFC, CNRS) M2 MIR – Microsystèmes Multiphysiques v2.2, septembre 2017 34 / 34