Pour la simulation de problèmes multi   échelles comme les écoulements turbulents, il est possible de combiner deux algorithmes de Runge-Kutta explicites et implicites afin de traiter implicitement l'intégration temporelle seulement dans la direction normale à la paroi, et ainsi de relaxer la contrainte sur le pas de temps. On préserve aussi le caractère peu diffusif et peu dispersif du schéma, et on limite par ailleurs la taille du système à inverser. On présentera l'algorithme et les premiers résultats obtenus en aéroacoustique

BAILLY, Christophe

BOGEY, Christophe

KREMER, François

I-Revues

20e`me Congre`s Franc¸ais de Me´canique Besanc¸on, 29 aouˆt au 2 septembre 2011
E´tude de sche´mas de Runge-Kutta semi-implicites pour
la simulation des e´coulements de paroi
F. Kremer, C. Bogey, C. Bailly
Laboratoire de Me´canique des Fluides et d’Acoustique
UMR CNRS 5509, Ecole Centrale de Lyon, 69134 Ecully, France
Re´sume´ :
Pour la simulation de proble`mes impliquant un raffinement de maillage, deux algorithmes de Runge-
Kutta semi-implicites sont de´veloppe´s afin de traiter de manie`re implicite l’inte´gration temporelle dans
la direction raffine´e et de manie`re explicite dans les autres directions, et ainsi relaxer la contrainte sur
le pas de temps. On impose e´galement aux sche´mas d’eˆtre peu diffusifs et peu dispersifs. Ces sche´mas
peuvent ainsi eˆtre utilise´s pour la simulation de proble`mes multi-e´chelles comme les e´coulements de
couche limite turbulente compressible. Dans ce papier, l’algorithme et les premiers re´sultats obtenus
pour un cas test acoustique sont pre´sente´s.
Abstract :
For the simulation of problems involving mesh refinement, two semi-implicit Runge-Kutta algorithms
are developped in order to treat in an implicit way the time integration in refined direction and in
an explicit way in other directions, so that constraint on the timestep can be relaxed. Low-diffusive
and low-dispersive properties of the scheme are also imposed. Then, these schemes can be used for
simulation of multiscale problems such as compressible turbulent boudary layer flows. Algorithm and
first results for an acoustic test case are presented.
Mots clefs : Runge-Kutta semi-implicite ; faible dissipation et faible dispersion
1 Introduction
Les phe´nome`nes physiques reproduits par la simulation ae´roacoustique directe couvrent une large
gamme d’e´chelles spatiales et temporelles [7], pour lesquelles des me´thodes nume´riques spe´cifiques ont
e´te´ de´veloppe´es au cours des vingt dernie`res anne´es [3]. En particulier, plusieurs me´thodes d’inte´gration
temporelle ont e´te´ propose´es par Hu et al. [2], Stanescu & Habashi [5] et Berland et al. [1]. Ces me´thodes
de type Runge-Kutta (RK) explicite ont la particularite´ de pre´senter une faible dissipation et une faible
dispersion jusqu’a` des fre´quences proches de la fre´quence de coupure impose´e par le pas de temps.
Toutefois, dans un calcul impliquant une paroi, l’apparition de tre`s petites e´chelles spatiales dont la
stabilite´ doit eˆtre assure´e implique une se´ve`re re´duction du pas de temps. Une solution a` ce proble`me
est d’utiliser un sche´ma implicite, qui assure la stabilite´ des petites e´chelles tout en gardant un pas de
temps important. Cependant, ce gain sur le pas de temps risque d’eˆtre annule´ par la ne´cessite´ d’inverser
des syste`mes de tre`s grande taille, comprenant les trois dimensions de l’espace, ce qui implique un
temps de calcul excessivement long.
Il existe une solution interme´diaire qui consiste a` combiner deux sche´mas, l’un explicite et l’autre
implicite, ce dernier n’e´tant utilise´ que pour l’inte´gration temporelle dans la direction ou` le maillage est
le plus fin. Dans un calcul de couche limite, par exemple, on peut ne traiter implicitement que la de´rive´e
des flux normale a` la paroi. Le syste`me est alors beaucoup plus rapide a` inverser car il n’implique
qu’une dimension de l’espace. Un de ces sche´mas dits semi-implicites de´veloppe´ par Zhong [8] a e´te´
utilise´ pour la simulation aux grandes e´chelles d’une couche limite turbulente compressible par Suh et
al. [6].
1
20e`me Congre`s Franc¸ais de Me´canique Besanc¸on, 29 aouˆt au 2 septembre 2011
Dans ce papier, deux me´thodes semi-implicites adapte´es a` l’ae´roacoustique sont de´veloppe´es. Chacune
de ces me´thodes combine le sche´ma explicite RK46-Ber d’ordre 4 a` 6 sous-e´tapes de Berland et al. [1]
avec un sche´ma implicite a` 6 sous-e´tapes. La me´thode suivie pour le de´veloppement de ces sche´mas,
ainsi que les caracte´ristiques obtenues en termes de stabilite´ et de pre´cision sont pre´sente´es dans
la section 2. Ces caracte´ristiques sont ensuite e´tudie´es sur des cas test acoustiques fournis dans la
section 3. La section 4 conclut l’e´tude mene´e sur ces sche´mas semi-implicites.
2 De´veloppement de me´thodes semi-implicites
La strate´gie de de´veloppement des sche´mas s’inspire des travaux de Kennedy & Carpenter [4]. Ceux-ci
expriment un sche´ma semi-implicite a` s sous-e´tapes de la manie`re suivante :
ui = un +∆t
∑i
j=1 a
[I]
ij f(uj) + ∆t
∑i
j=1 a
[E]
ij g(uj)
un+1 = un +∆t
∑s
i=1 b
[I]
i f(ui) + ∆t
∑s
i=1 b
[E]
i g(ui)
(1)
On peut remarquer qu’un sous-sche´ma spe´cifique est de´die´ a` chacun des termes f et g. Dans nos
travaux, les coefficients (a
[E]
ij , b
[E]
i ) du sous-sche´ma explicite sont choisis e´gaux a` ceux du sche´ma
RK46-Ber d’ordre 4 a` 6 sous-e´tapes de Berland et al. [1]. Un sous-sche´ma implicite d’ordre 4 est
alors recherche´, tel que l’erreur de couplage entre les deux sous-sche´mas soit d’ordre 3. Les coefficients
(a
[I]
ij , b
[I]
i ) sont donc soumis aux contraintes suivantes [4] (par souci de clarte´, les contraintes (2) portant
sur l’ordre 4 du sous-sche´ma implicite sont e´crites sans l’exposant [I]) :
(O1)
∑s
i=1 bi = 1 (O2)
∑s
i=1 bici =
1
2
(O3)a 12
∑s
i=1 bic
2
i =
1
3! (O3)b
∑s
i,j=1 biaijcj =
1
3!
(O4)a 16
∑s
i=1 bic
3
i =
1
4! (O4)b
∑s
i,j=1 biciaijcj =
3
4!
(O4)c 12
∑s
i,j=1 biaijc
2
j =
1
4! (O4)d
∑s
i,j,k=1 biaijajkck =
1
4!
(2)
b
[I]
i = b
[E]
i
(C3)a 12
∑s
i=1 b
[E]
i c
[E]
i c
[I]
i −
1
3! = 0 (C3)b
∑s
i,j=1 b
[E]
i a
[E]
ij c
[I]
j −
1
3! = 0
(C3)c
∑s
i,j=1 b
[E]
i a
[I]
ij c
[E]
j −
1
3! = 0
(3)
avec ci =
∑i
j=1 aij et s = 6. Les jeux de coefficients (a
[E]
ij , a
[I]
ij , bi) (avec bi = b
[E]
i = b
[I]
i ) obtenus
constituent donc des sche´mas semi-implicites a` 6 sous-e´tapes d’ordre 3, note´s SIRK63 (SIRK signifiant
Semi-Implicit Runge-Kutta). L’ordre 3 est lie´ au couplage entre les deux sous-sche´mas.
Dans ce travail, deux sche´mas de type SIRK63 ont e´te´ construits. Le premier sche´ma posse`de des
proprie´te´s de stabilite´ assez fortes. Il est donc note´ SIRK63-S (S pour Stability). Le second sche´ma, note´
SIRK63-A (A pour Accuracy), a e´te´ de´termine´ avec des contraintes supple´mentaires pour ame´liorer
la pre´cision :
c
[I]
i = c
[E]
i
(C4)
∑s
i,j,k=1 bia
[I]
ij a
[E]
jk ck −
1
4! = 0
(4)
Le second sche´ma est toujours formellement d’ordre 3, mais les contraintes supple´mentaires ont pour
effet de diminuer l’amplitude des erreurs induites.
Dans les descriptions qui suivent, on fera re´fe´rence aux sous-sche´mas explicites et implicites d’un
sche´ma X par les notations X(E) et X(I) respectivement. L’erreur de dissipation 1 − |Gimp| en fonc-
tion de la pulsation adimensionne´e ωy∆t est montre´e sur la figure 1(a) pour les sous-sche´mas im-
plicites SIRK63-A(I) et SIRK63-S(I) de´veloppe´s dans ce papier, ainsi que pour le sous-sche´ma im-
plicite ASIRK-3C(I) de Zhong [8]. L’erreur de dissipation du sche´ma Runge-Kutta standard d’ordre 4
est e´galement fournie comme re´fe´rence.
Le sous-sche´ma implicite de Zhong est plus dissipatif que les sche´mas propose´s : sa dissipation est
supe´rieure d’un ordre de grandeur de celle de RK4 pour ωy∆t < pi/4, tandis que dans cette meˆme
gamme de pulsation adimensionne´e, le premier sous-sche´ma de´veloppe´, SIRK63-S(I), a une dissipation
2
20e`me Congre`s Franc¸ais de Me´canique Besanc¸on, 29 aouˆt au 2 septembre 2011
proche de celle de RK4. Il devient meˆme le´ge`rement plus pre´cis pour ωy∆t > pi/4. Le second sous-
sche´ma SIRK63-A(I), qui a e´te´ de´veloppe´ spe´cifiquement pour la pre´cision, montre une dissipation
infe´rieure de deux ordres de grandeur de celle de SIRK63-S(I), pour ωy∆t < pi.
L’erreur de dispersion |ω∗imp∆t − ωy∆t|/pi en fonction de ωy∆t, montre´e sur la figure 1(b), fait ap-
paraˆıtre les meˆmes tendances que celles observe´es sur l’erreur de dissipation. Le sous-sche´ma implicite
de Zhong est plus dispersif que les sche´mas propose´s. Parmi ces derniers, le sous-sche´ma SIRK63-A(I)
est de nouveau plus pre´cis que le SIRK63-S(I), pour ωy∆t < pi/4. Pour des pulsations addimentionne´es
supe´rieures, les deux sous-sche´mas de´veloppe´s pre´sentent une dispersion e´quivalente.
(a) (b)
pi/8 pi/4 pi/2 pi
10−6
10−4
10−2
100
ωy∆t
1−
|G
im
p|
pi/8 pi/4 pi/2 pi
10−6
10−4
10−2
100
ωy∆t
|ω*
im
p∆
t−
ω
y∆
t| /
 pi
Fig. 1 – (a) Dissipation, et (b) dispersion des sous-sche´mas implicites. · · · · ·× RK4 standard d’ordre 4,
* ASIRK-3C(I) d’ordre 3 de Zhong, + SIRK63-A(I) d’ordre 4, 3 SIRK63-S(I) d’ordre 4
L’erreur de dissipation des diffe´rents sous-sche´mas explicites en fonction de la pulsation adimensionne´e
ωx∆t est donne´e sur la figure 2(a). Le sous-sche´ma explicite de Zhong pre´sente encore l’erreur la plus
importante. Les deux sche´mas propose´s ont le meˆme sous-sche´ma explicite : RK46-Ber. Ce dernier
a une dissipation infe´rieure d’au moins deux ordres de grandeur de celle de RK4 pour ωx∆t < pi.
L’erreur de dispersion fournie sur la figure 2(b) fait apparaˆıtre globalement les meˆmes tendances.
En de´finitive, l’e´tude mene´e dans cette section permet de montrer que la pre´cision des sche´mas semi-
implicites fournis dans ce papier est plus grande que celle du sche´ma semi-implicite de Zhong, et au
moins e´quivalente a` celle de RK4.
(a) (b)
pi/8 pi/4 pi/2 pi
10−6
10−4
10−2
100
ω
x
∆t
1−
|G
e
xp
|
pi/8 pi/4 pi/2 pi
10−6
10−4
10−2
100
ω
x
∆t
|ω*
e
xp
∆t
−ω
x∆
t| /
 pi
Fig. 2 – (a) Dissipation, et (b) dispersion des sous-sche´mas explicites. · · · · ·× RK4 standard d’ordre 4,
* ASIRK-3C(E) d’ordre 3 de Zhong, SIRK63-S(E) et SIRK63-S(E) : RK46-Ber
3
20e`me Congre`s Franc¸ais de Me´canique Besanc¸on, 29 aouˆt au 2 septembre 2011
3 Application a` des cas test acoustiques
Les proprie´te´s des sche´mas semi-implicites sont e´tudie´es sur un cas de propagation d’impulsion acous-
tique sans e´coulement moyen. Les e´quations d’Euler 2D sont re´solues sur un maillage carte´sien et
anisotrope. On de´finit le rapport d’aspect AR = ∆x/∆y (AR pour Aspect Ratio) ou` ∆x et ∆y sont
les pas de discre´tisation spatiale dans les directions x et y respectivement. Plus AR est grand, et plus
le maillage est raffine´ dans la direction y. L’impulsion initiale a une demi-largeur e´gale a` 3∆x et une
amplitude de 10 Pa. Pour comparer les performances des diffe´rents sche´mas d’inte´gration temporelle,
l’erreur de discre´tisation spatiale doit eˆtre ne´gligeable devant l’erreur de discre´tisation temporelle.
On utilise donc un sche´ma aux diffe´rences finies centre´es sur 21 points, d’ordre 20. Pour les tests de
stabilite´, des conditions de pe´riodicite´ sont applique´es aux bords du domaine, ce qui permet d’utiliser
le sche´ma centre´ sur l’ensemble du maillage, et d’e´viter ainsi d’e´ventuelles instabilite´es cause´es par
l’utilisation de sche´mas de´centre´s.
Les simulations sont re´alise´es avec les sche´mas d’inte´gration temporelle suivants : SIRK63-S, SIRK63-A,
ASIRK-3C et RK4. Chacun est teste´ pour des rapports d’aspect de 2, 4, 8, 16 et 32, avec ∆x = 1 fixe´
et en faisant varier ∆y. Le nombre de CFL est de´fini par CFL = c∆t/∆x.
3.1 Stabilite´
Des tests de stabilite´ sont effectue´s pour trouver le nombre de CFL maximum CFLmax avant diver-
gence. Les re´sultats sont fournis sur la figure 3. Les axes sont en e´chelle logarithmique. Le sche´ma
RK4 e´tant explicite, son nombre de CFL maximum diminue progressivement, comme pre´vu, lorsque le
rapport d’aspect augmente, c’est a` dire lorsque ∆y diminue. Le sche´ma de Zhong pre´sente un nombre
de CFL maximum constant, e´gal a` 0.2. Le nombre de CFL maximum de SIRK63-A diminue lorsque
AR augmente, jusqu’a` se stabiliser a` CFLmax = 0.25 pour AR ≥ 16. Le sche´ma SIRK63-S, enfin, est
le plus stable, avec CFLmax = 1.2 quelquesoit la valeur de AR.
1 2 4 8 16 32
10−2
10−1
100
101
AR
CF
L m
a
x
Fig. 3 – Nombre de CFL maximum en fonction de AR, pour les sche´mas : ·····× RK4 standard, * ASIRK-3C
de Zhong, + SIRK63-A, 3 SIRK63-S
3.2 Pre´cision
Pour les tests de pre´cision, on re´alise des calculs aux meˆmes rapports d’aspect que pre´ce´demment. La
simulation est mene´e jusqu’a` t = 32∆x/c ou` c est la vitesse du son, pour des nombres de CFL variant
entre 0.1 et 1. Les re´sultats sont compare´s a` des calculs de re´fe´rence. Les maillages e´tant diffe´rents
pour chaque rapport d’aspect, un calcul de re´fe´rence est re´alise´ sur chaque maillage, avec un nombre
de CFL tre`s faible. Il est par exemple e´gal a` 10−2 pour AR = 1, et a` 10−4 pour AR = 32. La solution
donne´e par un des calculs de re´fe´rence est repre´sente´e sur la figure 4.
On de´finit alors l’erreur globale de la manie`re suivante :
Eglob =
√∫∫
S
(p− pref )2 ds√∫∫
S
(pref − p0)2 ds
(5)
4
20e`me Congre`s Franc¸ais de Me´canique Besanc¸on, 29 aouˆt au 2 septembre 2011
Fig. 4 – Solution de re´fe´rence a` t = 32∆x/c. Repre´sentation de p− p0 pour une amplitude ∆p = ±1 Pa
ou` pref est issu du calcul de re´fe´rence. La figure 5 pre´sente l’erreur obtenue pour les diffe´rents sche´mas,
pour des rapports d’aspect de 1 et 32. Les axes sont en e´chelle logarithmique. La comparaison des
figures 5(a) et 5(b) montre que les courbes sont sensiblement les meˆmes aux deux rapports d’aspect.
Seule une diffe´rence est note´e pour le nombre de CFL maximum de SIRK63-A et de RK4. On peut donc
se limiter a` l’analyse des re´sultats obtenus avec AR = 1, donne´s sur la figure 5(a). Les algorithmes semi-
implicites pre´sentent tous la meˆme pente d’ordre 3. Le sche´ma de Zhong a une pre´cision relativement
faible : son erreur est de plus d’un ordre de grandeur supe´rieure de celle de RK4. Les sche´mas propose´s
sont en revanche plus pre´cis : celui de´veloppe´ pour sa stabilite´ (SIRK63-S) a une erreur similaire a`
celle de RK4, et celui de´veloppe´ pour sa pre´cision (SIRK63-A), a une erreur d’un ordre de grandeur
infe´rieure a` celle de RK4.
(a) (b)
0.1 0.2 0.4 0.8 1
10−8
10−6
10−4
10−2
100
CFL
E g
lo
ba
le
0.1 0.2 0.4 0.8 1
10−8
10−6
10−4
10−2
100
CFL
E g
lo
ba
le
Fig. 5 – Erreur globale en fonction du nombre de CFL : (a) AR=1, (b) AR=32. · · · · ·× RK4 standard d’ordre 4,
* ASIRK-3C d’ordre 3 de Zhong, + SIRK63-A d’ordre 3, 3 SIRK63-S d’ordre 3
4 Conclusions
Deux sche´mas semi-implicites d’ordre 3, SIRK63-A et SIRK63-S, sont propose´s dans ce papier. Ces
deux sche´mas utilisent le meˆme sous-sche´ma explicite : le sche´ma RK46-Ber de Berland et al. [1].
Leurs sous-sche´mas implicites respectifs sont inconditionnellement stables et permettent d’inte´grer les
e´quations dans la direction ou` les tailles de maille sont tre`s petites, sans diminuer le pas de temps.
L’erreur de dissipation et l’erreur de dispersion des sche´mas SIRK63-A et SIRK63-S sont infe´rieures
ou e´gales a` celles du sche´ma de Runge-Kutta d’ordre 4 (RK4) standard. Des tests effectue´s sur un
cas de propagation d’impulsion acoustique en deux dimensions permettent de retrouver ces proprie´te´s.
Des mesures de temps CPU ont e´te´ effectue´es sur des calculs non pre´sente´s ici. Ces mesures montrent
que le sche´ma SIRK63-S permet d’obtenir un temps de calcul jusqu’a` dix fois plus petit qu’avec un
sche´ma explicite, pour un couˆt en me´moire raisonnable.
5
20e`me Congre`s Franc¸ais de Me´canique Besanc¸on, 29 aouˆt au 2 septembre 2011
Le sche´ma SIRK63-S pourrait donc eˆtre utilise´ pour des calculs tels que la simulation de couche limite
turbulente compressible, ou` le raffinement du maillage pe´nalise fortement le pas de temps des sche´mas
explicites. Le sche´ma SIRK63-A se re´ve`le e´galement inte´ressant, notamment pour sa grande pre´cision.
Il ne´cessite toutefois un traitement ame´liore´ de la re´solution du syste`me implicite, pour le rendre plus
rapide.
Re´fe´rences
[1] Berland J., Bogey C., Bailly C. Optimized explicit schemes : Matching and boundary schemes and
4th-order Runge-Kutta algorithm In 10th AIAA/CEAS Aeroacoustics Conference, 10-12 May 2004,
Manchester, UK AIAA Paper 2004-2814
[2] Hu F.Q., Hussaini M.Y., Manthey J.L. Low-dissipation and low-dispersion Runge-Kutta schemes
for computational acoustics Journal of Computational Physics 124 177-191, 1996
[3] Colonius T., Lele S.K. Computational aeroacoustics : Progress on nonlinear problems of sound
generation Progress in Aerospace Sciences 40(6) 345-416, 2004
[4] Kennedy C.A., Carpenter M.H. Additive Runge-Kutta schemes for convection-diffusion-reaction
equations Applied Numerical Mathematics 44 139-181, 2003
[5] Stanescu D., Habashi W.G. 2N-Storage Low Dissipation and Dispersion Runge-Kutta Schemes for
Computational Acoustics Journal of Computational Physics 143 674-681, 1998
[6] Suh J., Frankel S.H., Mongeau L., Plesniak M.W. Compressible large eddy simulations of wall-
bounded turbulent flows using a semi-implicit numerical scheme for low Mach number aeroacoustics
Journal of Computational Physics 215 526-551, 2006
[7] Tam C.K.W. Computational aeroacoustics : issues and methods. AIAA Journal 33(10) 1788-1797,
1995
[8] Zhong X. Additive semi-implicit Runge-Kutta methods for computing high-speed nonequilibrium
reactive flows Journal of Computational Physics 128 19-31, 1996
6


Etude de schémas de Runge-Kutta semi-implicites pour la simulation des écoulements de paroi

http://documents.irevues.inist.fr/bitstream/2042/46369/1/cfm2011_936.pdf

Etude de schémas de Runge-Kutta semi-implicites pour la simulation des écoulements de paroi

Abstract

Similar works

Full text

Available Versions

I-Revues