ENSIM VAC 5A 2019/2020

2016/2017
Éléments finis de frontière

Introduction aux
techniques de conditions aux limites absorbantes
pour la modélisation numérique
de problèmes acoustiques en domaines non bornés
---
Adrien Pelat
 Position du problème

Propagation dans un
Problème de diffraction guide semi infini

Le domaine physique d'intérêt est ouvert ou partiellement bornée

Le domaine numérique de calcul doit être de taille finie
Comment tronquer le domaine de calcul ?
- sans introduire de réflexions parasites
- sans modifier la solution dans le domaine
- en assurant les performances numériques (coût, stabilité, ...)
 Approches utilisées
pour ce type de problèmes

Méthodes numériques appropriées à des calculs en domaines ouverts

● Boundary Element Method (BEM) : formulation intégrale, discrétisation uniquement
des frontières du domaine
● Infinite Element Method : formulation intégrale, discrétisation du domaine et des
frontières, propriétés particulières des éléments, non traité ici
Frontières artificielles dotées de conditions absorbantes dans des schémas
numériques classiques (FEM, FD, …)

● Frontières absorbantes (NRBC) ● Couches absorbantes

* *

 

 

Dériver une condition transparente exacte ?

Approximer cette condition exacte dans un schéma numérique ?
 Un peu de chronologie

Source pour cette chronologie : [Givolli,2003]

1970 1980 1990 2000 2010

Zero order NRBC

Low order NRBC

[Givoli,1992] Exact non-local NRBC
[Engquist,1977] High order local NRBC
[Bayliss,1980]
[Keller,1989] PML
[Givoli,1990]

Dans ce cours : [Berenger,1994]

● Ecriture de la C.A. exacte

● High order local NRBC

● Exact non local NRBC (opérateur DtN)

[Collino,1993]
● PML
 Dérivation de la condition
transparente
Problème acoustique dans un domaine semi infini :

( 1 ∂2
)
Δ− 2 2 p(x 1 , x 2 , t)=0 dans D
c ∂t
x2 ∞
∂ p( x =x )=0 sur Γ
2 0 N
 P ∂ x2 D
pt=0= p 0 dans D N x1
∂
pt =0=u0 dans D
∂t

On défini la TF de p par rapport à x1 et t :

p  k 1 , x 2 , =∫
0
∞
∫
∞

−∞
p  x 1 , x 2 ,t e ik x
1
dx
1
1e−i  t
dt

Application de la TF à l'équation des ondes :

 ...
( ) ( )
d2 d
+i √ k −k 1 p̂ =0 ℑ [ √ k −k 1 ] < 0
2 2 2 2
+(k
2
−k
2
̂
1 ) p =0
∀
2
dx 2 dx 2
 Dérivation de la condition
transparente
Condition transparente :
TF−1
( ) ( )
d t,x ∂ + Η p=0 ̂
Η . f =i √ k 2−k 21 . f̂
+i √ k −k 1 p̂ =0 avec
2 2 1

dx 2 ∂ x2

Problème acoustique équivalent dans domaine borné :

 
2
1 ∂
 − 2 2 p  x 1 , x 2 , t=0 dans D
c ∂t x2
∂  CLA
p x 2=x 0 =0 sur  N
∂ x2

D

 P
 
∂ H p=0 sur 
∂ x2 CLA
N x1

pt=0= p 0 dans D
∂
pt =0=u0 dans D
∂t
 Dérivation de la condition
transparente

( )
∂ + Η p=0 avec Η
Analyse de la condition transparente : ̂ . f =i √ k 2−k 21 . f̂
∂ x2
● Équation différentielle d'ordre 1 : Eq. d'ondes « aller » uniquement
● Reproduit des conditions d'espace libre (pas de réflexion)

● Solution du problème (P) dans D = solution du problème équivalent (P ) dans D 

● L'opérateur H est non local en espace et en temps

Opérateur non local : Opérateur local :

● Pour appliquer l'opérateur en 1 point, p doit ● Pour appliquer l'opérateur, seule la connaissance
être connu en tous points de la frontière. de p au voisinage du point est nécessaire
● Tous les nœuds de la frontière sont couplés ● Un nœud n'est couplé qu'à ses voisins immédiats
● Prise en compte à l'échelle globale du pb. ● Prise en compte à l'échelle de l'élément (mat.
● Peu maniable numériquement (matrices de élémentaires)
maillages non creuses) ● Nécessité de définir des DDL non physiques
supplémentaires (fonctions auxiliaires)

● L'opérateur H doit être approximé

 est discrétisé par un maillage (volume et frontières)
● Le domaine continu D

● Approcher l'opérateur par un autre opérateur pour le rendre plus maniable

 Zoologie des CLA

Performances numériques :
● Approximation non uniforme dans le plan (w , k1)
● Dimensionner la CLA pour une plage de valeurs (w , k1) : fréquence / angle d'incidence

Plusieurs type d'approximations pour des CLA locales :

● À l'ordre 1 (« zero order ») : analogue à une condition de Sommerfeld
● D'ordre faible : développement de Taylor (Engquist–Majda, ordre 2)
● D'ordre élevé : développement de Padé (Colino)

CLA exactes et non locales

● Opérateur Dirichlet to Neumann
 Approximation d'ordre faible   :
DL de la condition transparente

Réécriture de la condition transparente

● Réécriture valide pour les ondes aller
( ∂ + Η p=0
∂ x2 ) avec ̂
Η . f =i √ k 2−k 21 . f̂

√
2
k1 k1
ℑ [ √ k −k ] < 0
2 2
1 √ k −k = k 1−
2 2
1 2
= k √1−x
2
avec x=
k k

Développement de Taylor :
● Définition générale
n 
f ' ' x 0  2 f  x0 n n
f  x =f  x 0 f '  x 0  x−x 0   x−x 0  ...  x−x 0  o  x −x 0  
2! n!
● a
Application à la fonction (1+x) au voisinage de x=0 (position de la frontière artificielle):
a a a−1 2 a a−1a−2 3 a  a−1 a−2... a−n1 n
1x  =1ax x  x ... o  x 
2! 3! n!
● Pour nous, a=1/ 2 et x ≡−x 2
 Approximation à l'ordre 1

Réécriture de la condition transparente

● Réécriture valide pour les ondes aller
( ∂ + Η p=0
∂ x2 ) avec ̂
Η . f =i √ k 2−k 21 . f̂

√
2
k1 k1
ℑ [ √ k −k ] < 0
2 2
1 √ k −k = k 1−
2 2
1 2
= k √1−x
2
avec x=
k k

DL à l'ordre 1 :
●
 1− x ≃1 2

Expression de la CLA approchée :

TF−1
√ ( )
2
k 1 ̂
H.f ≃ik ̂f
t 1 ∂ ∂ + 1 ∂ p=0 sur Γ
H≃ik 1− ≃ik H.f ≃ CLA
k
2
c ∂t ∂ x 2 c ∂t
● Par cette approximation, l'opérateur est devenu différentiel
● Eq. diff. d'ordre 1 décrivant uniquement des ondes « aller »
● Eq. 1D dans la direction normale à la frontière

Condition totalement transparente pour des ondes en incidence normale à la frontière

La discrétisation des opérateurs différentiels en temps et en espace est aisée (FEM, FD)
 Approximation d'ordre faible

1 2
DL à l'ordre 2 [Enquist,1979] :  1− x ≃1− 2 x
2

● Expression de la CLA approchée :

TF−1
( )
2
1k ̂ t,x 1 ∂ 2 1 ∂2 1 ∂2
̂ ̂
H.f ≃ik f −ik
1
f
1
+ − p=0 sur ΓCLA
2k 2 2
c ∂t 2
c ∂t ∂ x2 2 2
∂ x1
● Prend en compte la propagation selon x1 d'ondes « aller » et « retour »
● Très bonne performances pour des ondes d'incidence non normale
→ Ex : R=0.03 pour un angle d'incidence de 45°

1 2 1 4
DL d'ordre supérieur (>2): √ 1−x ≃1− x + x
2

●
2 8
Pour augmenter les performances de transparence (en angle d'incidence, ...)
● Abouti à un opérateur différentiel d'ordre 4 très complexe qui rend le pb très mal posé

Approximation par développement de Taylor :

● Opérateur différentiel et local → problème bien posé et numériquement maniable

● Limitations → doit rester d'ordre faible pour conserver ses avantages

● Pratique et urilisé lorsque les performances requises sont raisonables

 Approximations d'ordre élevé
par fractions rationnelles
Développement en fraction rationnelles : 2
Lx
Expression du développement :  1− x ≃ 1−∑ l=1  l
2
●
2
1−l x
● Il existe plusieurs choix possibles pour les coefficients bl et al .

Approximants de Padé :
● Expressions des coefficients :  l = cos 2  l
2L1  βl =
2
2 L+1
sin2 (
lπ
2 L+1 )
Application à la condition transparente (x2=0) [Collino,1993] :

( √ ) k 21
[ ( )]
2
d d L k
dx 2
+i 1− 2 ̂p =0 + ik 1−∑l=1 βl 2 1 2 ̂p = 0 sur ΓCLA
k dx 2 k −α l k 1

Introduction d'une fonction auxiliaire : 2

Définie uniquement sur la frontière k2
● ^
ϕ(k 1,k )=
1
.^
p (k ¿ 0, k )
k ² −αl k 2
1
 Approximations d'ordre élevé
par fractions rationnelles
Report de la fonction auxiliaire dans la CLA :
2
k1
ϕ(k
̂ 1, k )= . ̂p ( k 1, x2 =0, k )
( d
)
L
2
k −αl k
2
+i k −∑l=1 ik βl ϕ̂ l ̂p = 0
1 dx 2

[ ( )]
2
d L k
+ ik 1−∑l=1 βl 2 1 2 p = 0  −l k 1  l = k 1 p  x 2 =0
2 2 2
dx 2 k −α l k 1

Écriture de la CLA après TF −1

t ,x1
∂p ∂p L ∂ l
 −∑ l=1 l = 0
∂ x2 ∂ t ∂t
2 2
∂ l ∂ l 2
∂ p
2
−l 2
= 2
Ensemble de L+1 équations couplées : ∂t ∂ x1 ∂ x1
● 1 eq. de transport selon x
2
● L eq. d'ondes 1D selon x
1
Nécessité de définir des CL en x1 pour les L équations
Ordre des opérateurs différentiels £ 2
L termes correctifs : prise en compte sur la frontière des variations en x1 de
l'onde incidente (incidence quelconque)
Comparaison des coefficients de
réflexions des différentes C.L.A.
 Les couches Absorbantes   :
Perfectly Matched Layers
Contexte :
● Initialement développé dans le domaine temporel pour des ondes EM [Bérenger, 1994]
● Principe applicable au domaine fréquentiel pour des ondes scalaires
Formulation dans le cas général dans le domaine fréquentiel :
● Eq. d'ondes 2D pour le champ j :

*
● Prolongement complexe dans la PML (domaine *)


● Nouvelle eq. d'ondes :

 Les couches Absorbantes   :
Perfectly Matched Layers
Cas le plus simple de PML :
● Coefficients de PML constants : s0,x = s0,y = t
● Paramètre de forme : b = 1
● Eq. d'ondes 2D pour le champ j :
[ ( )
1 ∂ 1 ∂
+
1 ∂ 1 ∂
(
τ ∂x τ ∂x τ ∂ y τ ∂ y )
+ K
]
2
p = 0

Précautions :
● Affecter des CL à l'extérieur des PML
● Dimensionner le coefficient de réflexion de la PML
PML =exp( ℑ[ τ]k x 2d)
● R

● Compromis entre d et t

● Finesse du maillage : bien reconstruire l'exp.

● Efficace uniquement dans l'axe de la PML !
 Introduction to PML waveguides :
Green function of the 2D half space 16

2D free half space

z

zs x
+
z=-l
 Formalism in 2D
Introduction to PML waveguides : 1
2
Leaky modes
Modal-FE method Cross-section
Green function of the 2D half space 16
3
4
Applications
Conclusion
discontinuity
Transposition to 3D

2D free half space

z

zs x
+
z=-l

Exact solution :
 Formalism in 2D
Introduction to PML waveguides : 1
2
Leaky modes
Modal-FE method Cross-section
Green function of the 2D half space 16
3
4
Applications
Conclusion
discontinuity
Transposition to 3D

2D free half space

z

zs x
+
z=-l

Exact solution :
 Formalism in 2D
Introduction to PML waveguides : 1
2
Leaky modes
Modal-FE method Cross-section
Green function of the 2D half space 16
3
4
Applications
Conclusion
discontinuity
Transposition to 3D

2D free half space 2D PML waveguide

z z
z=h
PML
zs x x
+
z=-l z=-l

Exact solution :
 Formalism in 2D
Introduction to PML waveguides : 1
2
Leaky modes
Modal-FE method Cross-section
Green function of the 2D half space 16
3
4
Applications
Conclusion
discontinuity
Transposition to 3D

2D free half space 2D PML waveguide

z z
z=h
PML
zs x x
+
z=-l z=-l

Exact solution : Transverse problem :

 Formalism in 2D
Introduction to PML waveguides : 1
2
Leaky modes
Modal-FE method Cross-section
Green function of the 2D half space 16
3
4
Applications
Conclusion
discontinuity
Transposition to 3D

2D free half space 2D PML waveguide

z z
z=h
PML
zs x x
+
z=-l z=-l

Exact solution : Transverse problem :

Solution for the Green function

[Knockaert, De Zutter, IEEE, 2002]

 Formalism in 2D
Introduction to PML waveguides : 1
2
Leaky modes
Modal-FE method Cross-section
Green function of the 2D half space 16
3
4
Applications
Conclusion
discontinuity
Transposition to 3D

2D free half space 2D PML waveguide

z z
z=h
PML
zs x x
+
z=-l z=-l
 < 1% (10 modes)
Exact solution : Transverse problem :

Solution for the Green function

[Knockaert, De Zutter, IEEE, 2002]

 1 Leaky modes Formalism in 2D
Formalism of the Modal-FE method : 2 Modal-FE method Cross-section
The problem 3
4
Applications
Conclusion
discontinuity
17 Transposition to 3D

z
Computational domain : z=h
2D PML waveguide, rigid boundaries PML
x

z=-l
Wave equation : x=0 x=L
Introduction of PML :
 1 Leaky modes Formalism in 2D
Formalism of the Modal-FE method : 2 Modal-FE method Cross-section
The problem 3
4
Applications
Conclusion
discontinuity
17 Transposition to 3D

z
Computational domain : z=h
2D PML waveguide, rigid boundaries PML
x

z=-l
Wave equation : x=0 x=L
Introduction of PML :

Discretization of the transverse problem (FEM) :

N nodes grid + First order polynomials
Development of the solution :
 1 Leaky modes Formalism in 2D
Formalism of the Modal-FE method : 2 Modal-FE method Cross-section
The problem 3
4
Applications
Conclusion
discontinuity
17 Transposition to 3D

z
Computational domain : z=h
2D PML waveguide, rigid boundaries PML
x

z=-l
Wave equation : x=0 x=L
Introduction of PML :

Discretized wave equation :

 1 Leaky modes Formalism in 2D
Formalism of the Modal-FE method : 2 Modal-FE method Cross-section
General modal solution 3
4
Applications
Conclusion
discontinuity
18 Transposition to 3D

z
z=h
PML
x

z=-l
General modal solution x=0 x=L

Constants and
 1 Leaky modes Formalism in 2D
Formalism of the Modal-FE method : 2 Modal-FE method Cross-section
General modal solution 3
4
Applications
Conclusion
discontinuity
18 Transposition to 3D

z
z=h
PML
x

z=-l
General modal solution x=0 x=L

Eigensolutions of

Constants and
 1 Leaky modes Formalism in 2D
Formalism of the Modal-FE method : 2 Modal-FE method Cross-section
General modal solution 3
4
Applications
Conclusion
discontinuity
18 Transposition to 3D

z
z=h
PML
x

z=-l
General modal solution x=0 x=L

Source condition Radiation impedance

and
are functions of , ,
Constants and
 1 Leaky modes Formalism in 2D
Formalism of the Modal-FE method : 2 Modal-FE method Cross-section
General modal solution 3
4
Applications
Conclusion
discontinuity
18 Transposition to 3D

z
z=h
PML
x

z=-l
General modal solution x=0 x=L

Remark :
The radiation impedance can be
algebraically propagated back

Constants and
 1 Leaky modes Formalism in 2D
Solution for an initial condition problem : 2 Modal-FE method Cross-section
Incident Gaussian beam 3
4
Applications
Conclusion
discontinuity
20 Transposition to 3D

Incident Infinite
Gaussian beam waveguide
PML

∞
 1 Leaky modes Formalism in 2D
2 Modal-FE method Cross-section
Comparison with a full FEM 2D computation 3 Applications discontinuity
4 Conclusion
25 Transposition to 3D

FEM : Modal-FE :
N ≃ 13000 N ≃ 100
(10 nodes/) (10 nodes/ )
 1 Leaky modes Formalism in 2D
Illustration in 2D : 2 Modal-FE method Cross-section
Scattered wave by an irregular ground 3
4
Applications
Conclusion
discontinuity
26 Transposition to 3D

PML

1 discontinuity => N discontinuities

 1 Leaky modes Formalism in 2D
Transposition to 2 Modal-FE method Cross-section
3D open irregular waveguides 3
4
Applications
Conclusion
discontinuity
27 Transposition to 3D

Computational domain :
3D canyon with piecewise constant cross-section
Rigid walls

Wave equation :
 1 Leaky modes Formalism in 2D
An illustration in 3D : 2 Modal-FE method Cross-section
Street like geometry 3
4
Applications
Conclusion
discontinuity
29 Transposition to 3D

Typical competitiveTitle:/home/adrien/DATADRI/THESE/Multi
effects Title:/home/adrien/DATADRI/THES
Title:/home/adrien/DATADRI/THESE/Multi
guided waves ≠ radiated waves z = 0.1
Creator:MATLAB, z = 0.5
Creator:MATLAB,
The zThe
= 0.9
Mathworks, Inc. Mathworks,
Vers
Creator:MATLAB, The Mathworks,Title:/home/adrien/DATADRI/THESE/Mul
Inc. Vers
CreationDate:07/06/2009 CreationDate:07/06/2009
CreationDate:07/06/2009 09:07:38
09:06:47Creator:MATLAB, 09:07:58
The Mathworks, Inc. Ve
LanguageLevel:2 LanguageLevel:2 CreationDate:07/06/2009
LanguageLevel:2 09:07:58
Field partially confined LanguageLevel:2
wthin the street

+
Title:/home/adrien/DATADRI/THESE/MultiTitle:/home/adrien/DATADRI/THESE
Title:/home/adrien/DATADRI/THESE/Multi
Creator:MATLAB, zThe
Creator:MATLAB, zCreator:MATLAB,
= 1.1Mathworks, Inc.The
= Vers
1.5 The
Mathworks, zInc. Mathworks, In
= 2 Vers
Also partially radiated
CreationDate:07/06/2009 09:08:22 CreationDate:07/06/2009
CreationDate:07/06/2009 09:08:46 09:09:15
towards the sky
LanguageLevel:2 LanguageLevel:2LanguageLevel:2
 Bibliographie

CLA d'ordre 1 :

[Givoli,1992] Numerical Methods for Problems in Infinite Domains, Elsevier, Amsterdam, 1992.
[Givoli,1991] Non-reflecting boundary conditions: a review, J. Comput. Phys. 94 (1991) 1–29.
CLA d'ordre faible :
[Engquist,1979] Radiation boundary conditions for acoustic and elastic calculations, Commun. Pure
Appl. Math. 32 (1979) 313–357.
Opérateur DtN :
[Keller,1989] Exact non-reflecting boundary conditions, J. Comput. Phys. 82 (1989) 172–192.
[Givoli,1990] Non-reflecting boundary conditions for elastic waves, Wave Motion 12 (1990) 261–279.
CLA d'ordre élevé :
[Givoli,2004] High-order local non-reflecting boundary conditions: a review, Wave Motion 39 (2004)
319–326 (+ toute la biblio citée dans l'article)
[Colino,1992] Conditions absorbantes d'ordre élevé pour les modèles de propagation d'onde dans des
milieux rectangulaires, rapport de recherche de l'INRIA n°1790, nov. 1992
[Collino,1993] F. Collino, High order absorbing boundary conditions for wave propagation models.
Straight line boundary and corner cases, in: R. Kleinman, et al. (Eds.), Proceedings of the Second
International Conference on Mathematical and Numerical Aspects ofWave Propagation, SIAM,
Delaware, 1993, pp. 161–171..
PML
[Bérenger,1994] A perfectly matched layer for the absorption of electromagnetic waves, J. Comput. Phys.
114 (1994) 185–200.