Université Numérique Cheikh Hamidou Kane

(UN-CHK)

Pôle Sciences, Technologie et Numérique

Master Analyse Numérique et Modélisation

Mémoire de fin d’études

Pour l’obtention de grade de Master 2 en Analyse Numérique et Modélisation

Master délivré dans le cadre du centre national du calcul

scientifique (CNCS), et en partenariat avec ATOS
Résolution numérique des équations de Stokes par la méthode
de compressibilité artificielle en volumes finis en 2D

Réalisé par : Abdallah Abdoulaye MAHAMAT

Sous la supervision de : Dr Abdou Khadre DIOUF

Devant le jury composé de :

Président Pr Abdou SENE

Encadreur Dr Abdou Khadre DIOUF
Examinateur Dr Timack NGOM
Examinateur Dr Elhadj Ibrahima THIAM
Année universitaire : 2021-2022

i
Table des matières

Dédicace v

Remerciments vi

Résumé vii

Abstract viii

Introduction Générale 1

1 Généralités sur les équations aux dérivées partielles 5

1.1 Équations aux dérivées partielles du premier ordre . . . . . . . . . . . . . . . . . . . . 5
1.1.1 Équation d’advection . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.2 Équations aux dérivées partielles du second ordre . . . . . . . . . . . . . . . . . . . . 6
1.2.1 Équation de la chaleur . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.2.2 Équation des ondes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.3 Le système de Stokes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
1.4 Méthodes classiques de résolution des équations de Stokes . . . . . . . . . . . . . . . 8
1.4.1 Méthode semi-implicite pour les équations liées à la pression . . . . . . . . . . 8
1.4.2 Marqueur et cellule (Marker And Cell) . . . . . . . . . . . . . . . . . . . . . . 9
1.4.3 Compressibilité artificielle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9

2 Construction du schéma numérique par la méthode des volumes finis du système de Stokes 10
2.1 Méthode de compressibilité artificielle : . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2 Méthodes numériques de discrétisation . . . . . . . . . . . . . . . . . . . . . . . . . . 10
2.2.1 Méthode des différences finies . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2.2 Méthode des éléments finis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2.3 Méthode des volumes finis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.3 Présentation du maillage décalé . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.3.1 Maillage décalé . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.3.2 Gauss-Ostrogradski . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
2.3.3 Application de la méthode des volumes finis . . . . . . . . . . . . . . . . . . . 15
2.3.4 Discrétisation temporelle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.4 Convergence du schéma numérique issu de la MVF . . . . . . . . . . . . . . . . . . . 17
2.4.1 Consistance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.4.2 Etude de la stabilité par la méthode de Von neumann . . . . . . . . . . . . . . 18
2.4.3 Critère de convergence du schéma . . . . . . . . . . . . . . . . . . . . . . . . . 26

3 Essais numériques 27
3.1 Implémentation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.1.1 Formulation discrète des équations de Stokes . . . . . . . . . . . . . . . . . . . 27
3.1.2 Procédure de résolution numérique . . . . . . . . . . . . . . . . . . . . . . . . 27
3.1.3 Présentation du maillage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.2 Résultats numériques . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

ii
 3.2.1 Profils des composantes de la vitesse Ux et Vx . . . . . . . . . . . . . . . . . . 29
3.2.2 Tests de stabilité . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

Conclusion Générale 31

iii
Table des figures

2.1 Positions des composantes de vitesse U (rouge), V (bleu) et de la pression (noir) sur
la grille de calcul[4]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
2.2 Représentation complète des maillages décalés Mc , Mu et Mv : Positions des différents
jeux de coordonnées, de la pression (noir) et des composantes de vitesse U (rouge) et
V (bleu) ainsi que des volumes de contrôle associés[4]. . . . . . . . . . . . . . . . . . . 13
2.3 Positions des composantes de vitesse et de la pression sur le maillage. Composante
longitudinale de vitesse U (en rouge), composante normale de vitesse V (en bleu),
pression (en noir). Noeuds à l’intérieur du domaine de calcul (cercles pleins), noeuds
associés aux conditions aux limites (cercles vides)[4]. . . . . . . . . . . . . . . . . . . 14

3.1 Procédure de résolution numérique pour la méthode de compressibilité artificielle . . . 28

3.2 Génération du maillage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
3.3 Maillage en 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3.4 Le profil de chaque compossante de la vitesse suivant x . . . . . . . . . . . . . . . . . 29
3.5 Les composantes de la vitesse suivant x . . . . . . . . . . . . . . . . . . . . . . . . . . 30
3.6 Les composantes de la vitesse suivant x . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.7 Le temps de calcul . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

iv
 Dédicace
Je dédie ce mémoire à l’ensemble de ma famille et plus particulièrement à ma mère Soukaray
Abdarmane pour l’amour, la confiance, les conseils ainsi que le soutien inconditionnel qui m’a permis
de réaliser les études pour lesquelles je me destine.
Je dédie également ce mémoire à mon défunt père Abdallah Abdoulaye qui n’a pas pu voir ce travail.

v
 Remerciments
Tout d’abord, je tiens à remercier Dieu tout puissant, de m’avoir accordé une bonne santé ainsi
qu’un courage immense dans la réalisation de ce mémoire.
Je tiens à exprimer ma profonde gratitude envers le Recteur de l’Université Numérique Cheikh
Hamidou Kane le Professeur Moussa LO pour avoir créé un environnement académique stimulant
et propice à l’apprentissage. Sa vision et son engagement envers l’éducation ont été une source
d’inspiration constante tout au long de mon parcours de Master.
Mes remerciements vont également au Professeur Abdou SENE et au Professeur Mary Teuw
NIANE , les initiateurs de la formation ATOS et en particulier la formation Master Modélisation
Mathématique, Analyse et Simulation Numérique dont je suis ressortissant.
Je tiens à exprimer ma reconnaissance envers mon encadreur, Docteur Abdou Khadre DIOUF. Ses
conseils éclairés, sa patience infinie et son soutien constant ont été les piliers de mon parcours
académique. Sa capacité à guider mes recherches m’a permis d’explorer de nouvelles voies et m’a
permis de développer non seulement mes compétences en recherche, mais aussi ma confiance en moi.
Je tiens également à remercier l’ensemble du corps professoral du pôle sciences, technologies et
numérique de l’Université Numérique Cheikh Hamidou Kane en particulier notre responsable
Docteur Oumar DIOP (pour le suivi et les conseils), l’ensemble du corps professoral intervenant des
différentes universités du Sénégal et de l’étranger pour leurs enseignements enrichissants.
Mes camarades de classe méritent également une mention spéciale pour les échanges fructueux et
les discussions stimulantes que nous avons eus tout au long de cette formation. Leur soutien et leur
amitié ont rendu cette expérience encore plus mémorable.
Enfin, je n’oublierai jamais le soutien inconditionnel de ma famille et de mes amis, qui ont été mes
piliers tout au long de ce voyage académique.
Merci à tous ceux qui ont reçu de près ou de loin à la réalisation de ce mémoire. Votre impact sur
mon parcours restera gravé dans ma mémoire pour toujours.

vi
 Résumé
L’objectif de cette étude est de résoudre numériquement les équations de Stokes en utilisant la
technique de compressibilité artificielle, par la méthode des volumes finis dans un maillage décalée
en 2D et en particulier notre étude est orienté vers les études d’analyse de stabilité.

Le schéma d’Euler explicite à été appliqué pour approximer les dérivées partielles temporelles des
variables.
Pour étudier la stabilité on a utilisé la méthode de Von Neumann dont la méthodologie est basé sur
le calcul du module des valeurs propores de la matrice d’amplification si ceux derniers est inférieur
ou égal à un alors la stabilité est confirmé.

Aprés simulation les conditions d’analyse de stabilité ont été bien vérifié. Alors dans le cas de
notre étude, le schéma d’Euler explicite est conditionnellement stable.

Mots-clés : Equations de Stokes, Volumes finis, Compressibilité artificielle, maillage décalée en

2D, Analyse de stabilité au sens de Von Neumann.

vii
 Abstract
The aim of this study is to numerically solve the Stokes equations using the artificial compres-
sibility technique, by the finite volume method in a 2D offset mesh, and in particular our study is
oriented towards stability analysis studies.

The explicit Euler scheme has been applied to approximate the time partial derivatives of the
variables.
To study stability, we used the Von Neumann method, whose methodology is based on calculating
the modulus of the amplification matrix’s eigenvalues, if these are less than or equal to one, then
stability is confirmed.

After simulation, the stability analysis conditions have been verified. So in the case of our study,
the explicit Euler scheme is conditionally stable.

Key-words : Stokes equations, Finite volumes, Artificial compressibility, 2D offset mesh, Stability
analysis in the Von Neumann sense.

viii
Introduction Générale

Notre compréhension des phénomènes du monde réel et notre technologie sont aujourd’hui en
grande partie basées sur les équations aux dérivées partielles, qui sont notées en abrégé EDP[6].
C’est en effet grâce à la modélisation de ces phénomènes au travers d’EDP que l’on a pu com-
prendre le rôle de tel ou tel paramètre, et surtout obtenir des prévisions parfois extrêmement précises.

Une fois modéliser le phénomène, par fois, il est difficile de trouver de solution explicite. C’est
pourquoi avec l’apparition des ordinateurs extrêmement puissants permet néanmoins aujourd’hui
d’obtenir des solutions approchées. C’est ce qui s’est passé par exemple lorsque nous regardons
les prévisions météorologiques ou bien lorsque nous voyons des images animées d’une simulation
d’écoulement d’air sur l’ail d’un avion. Le rôle des mathématiciens est de construire des schémas
d’approximation, et de démontrer la pertinence des simulation en établissant des estimations a priori
sur les erreurs commises.

Les EDP ont été probablement formulées pour la première fois lors de la naissance de la méca-
nique relationnelle au cours du 17e siècle (Newton Leibniz...). En suite, le ≪catalogue≫ des EDP
s’est enrichi au fur et à mesure le développement des sciences et en particulier de la physique. S’il
ne faut retenir quelques noms, on se doit de cité celui de Euler, puis ceux de Navier et Stokes, pour
les équations de la mécanique des fluides, ceux de Fourier pour l’équation de la chaleur, de Maxwell
pour celles de l’éléctromagnétisme, Schrödinger et Heisenberg pour les équations de la mécanique
quantique, et bien sûr de Einstein pour les EDP de la théorie de la relativité.

En mécanique des fluides, les équations de Navier-Stokes sont des équations aux dérivées partielles
non linéaires qui décrivent le mouvement des fluides newtoniens (donc des gaz et de la majeure partie
des liquides). L’équation de Stokes est une forme simplifiée de l’équation de quantité de mouvement
contenue dans les équations de Navier-Stokes. L’équation de Stokes, décrit l’écoulement d’un fluide
newtonien incompressible et à faible nombre de Reynolds.

La résolution numérique des équations de Stokes en utilisant la méthode de compressibilité arti-

ficielle en volumes finis découle de la nécessité de surmonter les difficultés associées à la simulation
des écoulements visqueux incompressibles, qui se caractérisent par leur complexité et leurs diverses
conditions aux limites. Les équations de Stokes, cruciales pour comprendre les comportements des
fluides visqueux à faible vitesse et les phénomènes géophysiques, se prêtent rarement à une solution
analytique. En conséquence, la résolution numérique joue un rôle critique pour obtenir des solutions
précises et stables.

1
 L’analyse de stabilité de Von Neumann joue un rôle central dans cette approche, permettant d’éva-
luer la convergence et la stabilité des schémas numériques employés pour discrétiser les équations de
Stokes. En utilisant la méthode de Fourier, cette analyse examine comment les erreurs numériques
se propagent au fil des itérations, orientant ainsi la sélection de schémas appropriés et le réglage des
paramètres de simulation. La combinaison de la méthode de compressibilité artificielle en volumes
finis avec une analyse de stabilité rigoureuse permet d’atteindre des solutions numériques exactes et
fiables pour les équations de Stokes en 2D. Cette approche trouve des applications variées, allant de la
recherche scientifique à l’ingénierie, contribuant ainsi à une meilleure compréhension des écoulements
fluides et à des progrès technologiques significatifs.

Dans le cas de notre travail, les équations de Stokes ont été discrétisées sur une grille décalée.
Sur ce, la méthode des volumes finis est appliquée sur un maillage rectangulaire 2D et la technique
de compressibilité artificielle a été appliquer pour faciliter la résolution numérique des équations
discrétisées.
L’un des principaux objectifs du présent travail est de trouver les paramètres qui satisferont la
stabilité du schéma numérique. Par conséquent, ce travail de recherche peut être résumé comme
suit : Le chapitre I présent les généralités sur les équations aux dérivées partielles. Le chapitre II
présent la construction du schéma numérique par la méthode des volumes finie du système de Stokes.
Il comprend la technique de discrétisation utilisée pour les grilles décalées. Le chapitre III présent les
essais numériques. Et bien évidemment, notre travail sera clôturé par une conclusion générale et de
perspective.

2
 Revues Bibliographiques
Beaucoup des travaux publiés ont été élaborés dans le cadre de la résolution numérique des équa-
tions de Navier Stokes par la méthode de compressibilité artificielle, parmi ces travaux, on peut citer
ceux de :

P. Louda, K. Kozel, J. Přı́hoda[7] ont étudié numériquement les équations de Navier Stokes par
la méthode de compressibilité artificielle en volumes finis appliqués aux écoulements stables et ins-
tationnaires d’un fluide incompressible newtonien. Les cas instationnaires comprennent l’instabilité
auto-induite et l’instabilité forcée. La simulation en régime permanent d’écoulements turbulents 2D
et 3D sur une marche arrière confinée montre que, contrairement à la viscosité de tourbillon, le mo-
dèle 3D EARSM donne des résultats acceptables en raison de sa capacité à capturer les écoulements
secondaires. Dans la méthode à temps unique, la magnitude du paramètre de compressibilité artifi-
cielle est limitée par le haut pour des raisons de calcul et ne permet pas de prédire de manière fiable
les écoulements secondaires. En revanche, l’extension à temps double fonctionne bien pour les cas
d’essais laminaires ainsi que pour la simulation de jet synthétique libre turbulent.

[Link], [Link] [11] ont étudié numériquement les équations de Navier Stokes par la mé-
thode de compressibilité artificielle, tridimensionnelles et instables sur des grilles non échelonnées.
L’approche modifie le schéma d’itération standard à double temps et à compressibilité artificielle (AC)
en incorporant des idées provenant de formulations fractionnées basées sur la pression. La méthode
hybride FSAC (fractional-step/artificial-compressibility) qui en résulte est précise au second ordre et
fait progresser les équations de Navier-Stokes dans le temps via une procédure en deux étapes. Dans
la première étape, qui est identique à l’étape de convection-diffusion dans les méthodes FS basées
sur la pression, un première champ de vitesse préliminaire est calculé, qui n’est pas sans divergence.
Dans la deuxième étape, cependant, au lieu de dériver une équation de pression-Poisson comme dans
les méthodes FS, la projection du champ de vitesse dans l’espace vectoriel solénoı̈dal est mise en
œuvre en utilisant une formulation AC à double pas de temps. Contrairement aux formulations AC
standard à double temps, où les itérations à double temps sont effectuées avec l’ensemble du système
non linéaire, dans le schéma FSAC, les termes convectifs et visqueux ne sont calculés qu’une fois par
pas de temps physique. Les expériences numériques montrent que la méthode proposée fournit des
solutions précises du second ordre et nécessite beaucoup moins de temps de calcul que la formulation
AC standard largement utilisée.

[Link], [Link], [Link], A.F. Di Rienzo [1] ont étudié numériquement les équations
de Navier Stokes par la méthode de compressibilité artificielle en différence finie sur un maillage car-
tésien régulier, en analogie avec la méthode de Boltzmann en treillis (LBM). Les efforts préliminaires
vers des implémentations optimales ont montré que LW-ACM est capable d’une vitesse de calcul
similaire à celle du LBM optimisé (BGK). De plus, la demande en mémoire est significativement
nettement inférieure à celle de (BGK-) LBM. Des repères bidimensionnels et tridimensionnels sont
étudiés, et une étude comparative approfondie entre l’approche actuelle et l’état de l’art est réalisée.
Les preuves numériques suggèrent que LW-ACM représente une excellente alternative en termes de
simplicité, de stabilité et de précision.

3
 T. Ohwada, P. Asinari, D. Yabusaki [10] ont étudié numériquement les équations de Navier
Stokes par la méthode de compressibilité artificielle et la méthode de Boltzmann en treillis qui
donnent les solutions suivantes dans la limite de la disparition du nombre de Mach. L’inclusion de la
viscosité apparente est considérée comme l’une des raisons du succès de la méthode de Boltzmann en
treillis, du moins dans le cas de la méthode 2D. Dans le présent article, la robustesse de la méthode
de compressibilité artificielle est améliorée en introduisant un nouveau terme de dissipation, qui
rend possible le calcul d’un nombre de cellules de Reynolds élevé. L’augmentation de la stabilité
est également confirmée dans l’analyse de stabilité linéaire ; la magnitude des valeurs propres est
considérablement réduite pour une faible résolution. drastiquement réduite pour une faible résolution.
Des comparaisons sont faites avec la méthode de Boltzmann en treillis. Il est confirmé que l’ACM
fortifié est plus robuste et plus précis que la méthode de Boltzmann en treillis.

4
Chapitre 1

Généralités sur les équations aux dérivées

partielles

Une équation dans laquelle figurent une fonction u de plusieurs variables dépendantes x1 , .., xn et
des dérivées partielles de u par rapport à ces variables, c’est-à-dire, une équation de la forme :

∂u ∂u ∂ 2u ∂ mu
F (x1 , .., xn , u, 2 , ..., 2 , , ..., 2 m ) = 0, (1.1)
∂ x1 ∂x1 ∂x1 ∂x2 ∂ xn

est une équation aux dérivées partielles.

Une telle équation est dite d’ordre m quand elle contient au moins une dérivée d’ordre m sans
en contenir d’autres d’ordre supérieur. Toute fonction u(x1 , .., xn ) qui satisfait identiquement à cette
équation est une solution de celle-ci. L’équation ci-dessus est dite linéaire lorsque F est une combi-
naison linéaire de u et ses dérivées.

1.1 Équations aux dérivées partielles du premier ordre

Définition 1.1.1. Une équation dans laquelle figurent une fonction u de plusieurs variables dépen-
dantes x1 , .., xn et des dérivées partielles du premier ordre de u par rapport à ces variables, c’est-à-
dire, une équation de la forme :

∂u ∂u
F (x1 , .., xn , u, , ..., ) = 0, (1.2)
∂x1 ∂xn

est dite une équation aux dérivées partielles (en abrégé : EDP ) du premier ordre.
Toute fonction u(x1 , .., xn ) qui satisfait identiquement à cette équation est une solution de celle-ci.
Dans le cas de deux variables x, y, on a

∂u ∂u
F (x, y, u, , ) = 0, (1.3)
∂x ∂y

5
1.1.1 Équation d’advection
On note Ω un ouvert borné de Rd , d=1,2 ou 3 T un réel, T > 0, v la vitesse d’advection et c(x,t)
la concentration au point x au temps t avec x ∈ Ω. L’équation d’advection est donnée par :

 ∂c ⃗
(x, t) + v.∇c(x, t) = 0 x ∈ Ω, t ∈]0, T [
 ∂t


c(x, 0) = c0 (x) x∈Ω (1.4)


x ∈ ∂Ω, t ∈ [0, T [

c(x, t) = 0

L’interprétation physique de l’équation d’advection est la suivante : on considère un canal ou un

fleuve représenté par Ω et qui s’écoule à la vitesse v. La quantité c(x,t) représente la concentration
au point x et à l’instant t d’un produit polluant qui a été déversée dans le canal à t=0 selon le profil
c0 . Le polluant est transporté par l’écoulement et à un temps t fixé, la fonction x → c(x, t) modélise
la répartition de polluant le long du canal.

1.2 Équations aux dérivées partielles du second ordre

Définition 1.2.1. Soit u une fonction de deux variables x et y. On appelle équation aux dérivées
partielles du second ordre, une relation de la forme

∂u ∂u ∂ 2 u ∂ 2 u ∂ 2 u
F (x, y, u, , , , , ) = 0, (1.5)
∂x ∂y ∂x2 ∂x∂y ∂y 2

faisant intervenir u et ses dérivées partielles jusqu’à l’ordre 2.

La forme générale d’une équation aux dérivées partielles linéaire d’ordre 2 à coiefficients constants
est :

∂ 2u ∂ 2u ∂ 2u ∂u ∂u
α 2
+ β + γ 2
+δ +ϵ + ζu = f, (1.6)
∂x ∂x∂y ∂y ∂x ∂y

Où α,β,γ,δ,ϵ,ζ sont des constantes et f : D → R est une application appelée le second membre de
l’équation(1.6).
On dit que l’équation (1.6) est :
− elliptique si 4αλ − β 2 > 0,
− parabolique si 4αλ − β 2 = 0,
− hyperbolique si 4αλ − β 2 < 0,
−− L’équation de Laplace(ou Poisson) ∆u = 0(∆u = f ) est elliptique.
2
−− L’équation de la chaleur ∂u ∂t
− c ∂∂xu2 = 0 est parabolique (c>0).
2 2
−− L’équation de des ondes ∂∂t2u − c2 ∂∂xu2 = 0 est hyperbolique.

1.2.1 Équation de la chaleur

L’équation de la chaleur en une dimension d’espace est donnée par l’équation aux dérivées par-
tielles suivante :

6
 
∂u 2
 (x, t) = c ∂∂xu2 (x, t), x ∈ R, t > 0
 ∂t


u(x, 0) = u0 , x∈R (1.7)



u(0, t) = u(t), t>0

où c>0 est une constante donnée, u est une fonction inconnue réelle de deux variables réelles x et
t. Cette fonction u=u(x,t) représente la température dans un conducteur de dimension un. La valeur
de u(x,t) dépend du temps t ≥ 0 et de la position x, u0 est la condition initiale et u(t) est la condition
aux bords.
L’interprétation physique de l’équation de la chaleur est la suivante : Ces équations décrivent au
cours du temps t de la température u d’un milieu continu homogène soumis à une source de chaleur
nulle lorsque les constantes physiques sont les prises égales à c, la température étant fixé égale u(t)
au cours du temps sur la frontière du milieu continu et la température initiale étant donné égale à
u0 .

1.2.2 Équation des ondes

L’équation des ondes modélise des phénomènes de propagation ou de vibration. Par exemple, en
deux dimensions d’espace, elle est un modèle pour étudier les vibrations d’une membrane élastique
tendue (comme la peau d’un tambour). En une dimension d’espace, elle est aussi appelée équation
des cordes vibrantes. Au repos, la membrane occupe un domaine plan Ω. Sous l’action d’une force
normale à ce plan d’intensité f, elle se déforme et son déplacement normal est noté u. On suppose
qu’elle est fixée sur son bord, ce qui donne une condition de Dirichlet. L’équation des ondes dont u
est solution est donnée par :

∂2u


 ∂t2
− ∆u = f dans Ω × R+
∗


∂Ω × R+

 u=0 sur ∗
(1.8)


 u(t = 0) = u0 dans Ω


 ∂u (t = 0) = u

dans Ω
∂t 1

1.3 Le système de Stokes

Le système de Stokes modélise l’écoulement d’un fluide visqueux incompressible à petite vitesse.
On suppose que le fluide occupe un domaine Ω et qu’il adhère à la paroi de celui-ci, c’est-à-dire que
sa vitesse est nulle sur la paroi (ce qui conduit à une condition aux limites de Dirichlet). Sous l’action
d’une force f(x) (une fonction de Ω dans RN ), la vitesse u(x) (un vecteur) et la pression (scalaire)
sont solutions de



 ∇p − µ∆u = f dans Ω

∇.u = 0 dans Ω (1.9)



u = 0 sur ∂Ω

7
 Où µ >0 est la viscosité du fluide. Remarquons qu’en plus des N équations ∇p−µ∆u = f (corres-
pondant à la conservation de la quandité de mouvement), il y a une autre équation divu=0 appelée
condition d’incompressibilité (qui correspond à la conservation de la masse). Si la dimension d’espace
est N=1, le système de Stokes est sans intérêt car on voit facilement que la vitesse est nulle et que la
pression est une primitive de la force. Par contre en dimension N ≥ 2, le système de Stokes a bien
un sens : en particulier, il existe des champs de vitesses incompressibiles non triviaux (prendre, par
exemple, un rotationnel).

Un écoulement de Stokes (ou écoulement rampant) caractérise un fluide visqueux qui s’écoule
lentement en un lieu étroit ou autour d’un petit objet, dont les effets visqueux dominent alors sur les
effets inertiels. Il est en effet régi par une version simplifiée de l’équation de Navier-Stokes, l’équation
de Stokes, dans laquelle les termes inertiels sont absents. Le nombre de Reynolds mesure le poids
relatif des termes visqueux et inertiel dans l’équation de Navier-Stokes et est faible (beaucoup plus
petit que 1) dans un écoulement de Stokes.
L’équation de Stokes permet en particulier de décrire la décantation de très petites particules
dans les liquides voire les gaz (cas des cristaux de glace décantant dans la haute atmosphère ou du
brouillard commun), ainsi que les écoulements de liquide dans les dispositifs microfluidiques. Les
écoulements de Couette et de Poiseuille sont aussi décrits par cette équation.

L’une des propriétés intéressantes des équations incompressibles de Stokes est que la variable
primitive p n’apparaı̂t pas dans l’équation dite de continuité. Cela implique que l’information de
pression ne peut pas être incorporée dans l’équation de continuité, ce qui pose un sérieux problème.
Par conséquent, tout schéma de discrétisation doit intégrer une sorte de stratégie pour résoudre ce
problème.

Cependant, les méthodes numériques pour les écoulements incompressibles ont connu un déve-
loppement rapide au cours des dernières décennies. Nous présentons ci-dessous quelques-unes des
méthodes qui ont acquis une bonne réputation pour la résolution des équations incompressibles de
Navier-Stokes :

1.4 Méthodes classiques de résolution des équations de Stokes

1.4.1 Méthode semi-implicite pour les équations liées à la pression
Il s’agit de l’un des algorithmes utilisés pour résoudre les équations incompressibles de Navier-
Stokes. Il s’agit de l’acronyme de Semi-Implicit Method for Pressure Linked Equations(Méthode semi-
implicite pour les équations liées à la pression) [2]. Cette méthode utilise une équation de correction
de pression pour mettre à jour le champ de pression. En utilisant les corrections de pression, l’équation
de continuité est résolue pour la correction des vitesses. Ces corrections sont utilisées pour mettre à
jour le champ de pression et de vitesse jusqu’à ce que l’on obtienne une solution convergente où les
corrections sont négligeables.

8
1.4.2 Marqueur et cellule (Marker And Cell)
Cette méthode tire ses origines dans (Harlow and Welch 1965 et Welch and al 1966). Elle utilise
un maillage décallé. La pression, la vitesse suivant X et la vitesse suivant Y ne sont pas calculées aux
mêmes points. La méthode (MAC) est une technique numérique utilisée pour résoudre des équations
aux dérivées partielles (EDP) en dynamique des fluides numérique (CFD). L’idée de base de la
méthode MAC est de discrétiser le domaine de calcul en une grille de cellules, puis de calculer les
vitesses aux faces et la pression du fluide au centre cellulaire [8], ou ≪ points marqueurs ≫.
La méthode MAC est couramment utilisée pour résoudre les équations de Navier-Stokes, qui
régissent le comportement des écoulements de fluides. En général, l’avantage de la grille décalée est
qu’il n’y a pas de couplage entre vitesse et pression. De plus, on n’a pas besoin d’une condition de
pression sur la frontière.

1.4.3 Compressibilité artificielle

Développée par Chorin [3] à la fin des années 1960, l’une des méthodes les plus faciles à mettre
en œuvre pour les équations incompressibles de Navier-Stokes, qui introduit un terme dérivé de la
pression dans l’équation de continuité permettant de résoudre un champ de pression actualisé. D’une
certaine manière, le système d’équations modifié, appelé équations de compressibilité artificielle,
ressemble souvent aux équations gouvernantes d’un cas d’écoulement compressible. La méthode de
compressibilité artificielle (ACM) offre plusieurs avantages pour la résolution numérique des équa-
tions de Stokes.

9
Chapitre 2

Construction du schéma numérique par la

méthode des volumes finis du système de
Stokes

2.1 Méthode de compressibilité artificielle :

Soit Ω un domaine régulier comme mentionné dans le chapitre précédent, la méthode de com-
pressibilité artificielle (ACM) sera utilisée dans le présent travail pour résoudre les équations in-
compressibles de Stokes. Le système de Stokes incompressible et adimentionnel par la méthode de
compressibilité artificielle peut s’écrire comme suit :

∂U


 ∂t
+ ∇P − ∆U = f dans Ω


 1 ∂P + ∇.U = 0

dans Ω
β ∂t
(2.1)


CI dans Ω



CL de Dirichlet priodique sur ∂Ω

Où β est appelé le paramètre de compressibilité artificielle.

Pour résoudre numériquement ce système, il existe plusieurs méthodes numériques pour résoudre ce
système, dans notre cas, on s’intéresse en particulier à la méthode des volumes finis.

2.2 Méthodes numériques de discrétisation

Les méthodes numériques de discrétisation qui sont la méthode des différences finies, la méthode
des éléments finis et la méthode des volumes finis.
Ces méthodes consiste à ramener la résolution du système d’équations aux dérivées partielles dans
le domaine d’étude, avec des conditions aux limites appropriées, à celles d’un système d’équations
algébriques.

10
2.2.1 Méthode des différences finies
La méthode des différences finies, présente une technique de résolution des équations aux dérivées
partielles, par l’approximation de dérivées par des différences finies.
Cette méthode consiste à subdiviser le domaine d’étude en un nombre déterminé de nœuds et
à représenter la fonction recherchée en chacun des nœuds du domaine par un développement limité
en série de Taylor [5]. Ainsi, l’équation différentielle est transformée en équation algébrique pour
chaque nœud. La résolution du système d’équations algébriques permet d’obtenir la distribution de
la fonction étudiée dans le domaine d’étude.
La méthode de différence finie ne permet pas la prise en compte des conditions de passage d’un milieu
physique à un autre et des non-linéarités, cela nécessite un traitement spécifique. D’autre part, elle
s’adapte mal aux objets de la géométrie complexe à cause de la rigidité du maillage.

2.2.2 Méthode des éléments finis

La méthode des éléments finis est de toutes les méthodes de discrétisation la plus utilisée car
elle peut traiter des problèmes de géométrie complexe, elle couvre de nombreux domaines de la phy-
sique, les moyens informatiques actuels la rende facile de mise en œuvre[8]. Cette méthode consiste à
transformer les équations différentielles sous des formes intégrale basé sur le concept de minimisation
d’une quantité (comme l’énergie. . .), conduisant à la solution exacte. Le principe fondamental de la
méthode des éléments finis consiste à : définir une partition du domaine d’étude c’est-à-dire subdiviser
le domaine d’étude en régions élémentaires (Éléments Finis) ; Représenter la fonction inconnue sur
chacun de ces éléments par une approximation polynomiale ; Construire les formes intégrales ; Mini-
miser l’intégrale ; Une organisation matricielle des calculs ; Une résolution du système algébrique ; La
MEF est une méthode très puissante pour la résolution des équations aux dérivées partielles surtout
dans les géométries complexes. Sa mise en œuvre par contre, est assez compliquée et nécessite une
place mémoire assez importante.

2.2.3 Méthode des volumes finis

La méthode des volumes finis est basée sur l’intégration de l’équation différentielle partielle pour
un volume de contrôle[9]. Cette méthode consiste à subdiviser le domaine d’étude en volumes de
contrôle de telle manière que chaque volume entoure un nœud principal P. Chaque volume est délimité
par quatre interfaces (e,w,s,n). Chaque nœud principal P et limité par des nœuds voisins (E,W,S,N).
L’équation différentielle est intégrée dans chaque volume de contrôle. Pour calculer l’intégrale dans
ce volume de contrôle l’inconnue est représentée à l’aide d’une fonction d’approximation (linéaire, pa-
rabolique, exponentielle, de puissance etc.) entre deux nœuds consécutifs. Ensuite, la forme intégrale
est discrétisée dans le domaine d’étude. Le résultat de discrétisation donne une équation algébrique
qui est constituée par les valeurs nodales.

2.3 Présentation du maillage décalé

Les méthodes numériques présentées sont basées sur une approche volumes finis d’ordre un en
temps et deux en espace. La discrétisation spatiale des équations est réalisée sur un maillage décalé

11
en se basant sur les schémas numériques. La discrétisation temporelle repose sur la méthode d’Euler
explicite.

2.3.1 Maillage décalé

La représentation discrète des champs de vitesse et de pression est effectuée sur maillages décalés.
Chaque grandeur (pression et composantes de vitesse) est définie sur son propre maillage et les
maillages se trouvent décalés les uns par rapport aux autres d’une demi-maille suivant une direction
particulière . Autrement dit, si on définit la pression P au centre de la maille, les composantes de
vitesse U et V sont localisées au centre des interfaces de la maille (figure 2.1).

Figure 2.1 – Positions des composantes de vitesse U (rouge), V (bleu) et de la pression (noir) sur
la grille de calcul[4].

Considérons un domaine bidimentionnel rectangulaire D, de longueur Lx et de hauteur Ly . La

position de chaque point du domaine est définie dans un repère cartésien orthonormé dont l’origine
est placée dans le coin inférieur gauche par les coordonnées x ; y. Dans un souci de simplification, le
maillage est considéré cartésien, orthogonal et la discrétisation du domaine conduit à Nx × Ny mailles
suivant les deux directions. Le maillage associé à la pression P, également commun à l’ensemble des
grandeurs scalaires, que l’on notera Mc , admet les caractéristiques suivantes (figure 2) :
– les coordonnées discrètes xc (i) et yc (j) définissent la position du centre de la maille Mc (i, j), où
sont localisées les grandeurs scalaires. (Le maillage étant cartésien, x ne dépend que de i et y que de
j)
– les coordonnées xi (i) et yi (j) définissent la position des interfaces supérieures de la maille (i ; j) et
serviront à repérer les composantes de vitesse.

Les relations entre les coordonnées des interfaces et centrées sont :

xi (i − 1) − xi (i)
xc (i) =
2

yi (j − 1) − yi (j)
yc (j) =
2
– les dimensions de chaque maille Mc (i, j) sont définies par :

12
 ∆Xci = xi (i) − xi (i − 1)

∆Ycj = yi (j) − yi (j − 1)

Le volume de contrôle de la maille (la surface en 2D) est défini comme le produit des tailles de maille
dans chacune des directions :
∆Vci,j = ∆Xci · ∆Ycj

A partir des caractéristiques de Mc , il est possible de définir les deux maillages décalés Mu et Mv
associés aux composantes de vitesse U et V , respectivement. Les caractéristiques de chacune des
grilles sont les suivantes (figure 2.2) :
– Pour Mu
Coordonnées de Ui,j : (xi (i), yc (j))
Dimensions de maille : ∆Xui = xc (i + 1) − xc (i) et ∆Yui = yc (j) − yc (j − 1) = ∆Ycj
Volume de contrôle : ∆Vui,j = ∆Xui × ∆Yuj
– Pour Mv
Coordonnées de Vi,j : (xc (i), yi (j))
Dimensions de maille : ∆Xvi = x(i + 1) − x(i) = ∆Xci et ∆Yvi = yc (j) − yc (j − 1)
Volume de contrôle : ∆Vvi,j = ∆Xvi × ∆Yvj

Figure 2.2 – Représentation complète des maillages décalés Mc , Mu et Mv : Positions des différents
jeux de coordonnées, de la pression (noir) et des composantes de vitesse U (rouge) et V (bleu) ainsi
que des volumes de contrôle associés[4].

Les mailles périphériques de la grille de calcul sont réservées à la gestion des conditions aux li-
mites, on les nomme mailles fictives ou fantômes (ghost cells). Si on considère que les trois grilles de
calcul ont les mêmes dimensions Nx × Ny (pour des raisons pratiques de programmation liés à la
taille des tableaux), alors les points discrets à l’intérieur du domaine de calcul sont :

13
Pour Mc : [2,Nx - 1] × [2,Ny - 1]
Pour Mu : [2,Nx - 2] × [2,Ny - 1]
Pour Mv : [2,Nx - 1] × [2,Ny - 2]

On remarquera que les maillages décalés associés aux composantes de vitesse sont définies avec
une maille de moins suivant la direction de la composante de vitesse concernée .
Dans ce cas, les points associés aux conditions aux limites sont en coı̈ncident avec les bords du
domaine de calcul. Sinon, ils se trouvent à l’extérieur du domaine de calcul et appartiennent aux
mailles fictives (ghost cells) (figure 2.3).

Figure 2.3 – Positions des composantes de vitesse et de la pression sur le maillage. Composante
longitudinale de vitesse U (en rouge), composante normale de vitesse V (en bleu), pression (en noir).
Noeuds à l’intérieur du domaine de calcul (cercles pleins), noeuds associés aux conditions aux limites
(cercles vides)[4].

2.3.2 Gauss-Ostrogradski
Théorème 2.3.1. Le théorème de divergence indique que le flux sortant d’un champ vectoriel à travers
une surface fermée est égal à l’intégrale en volume de la divergence sur la région à l’intérieur de la
surface.
Intuitivement, il indique que la somme de toutes les sources moins la somme de tous les puits donne
le débit net à travers une région.
ZZZ I I
⃗ dV =
∇.U ⃗ .ndΓ
U (2.2)
v

14
2.3.3 Application de la méthode des volumes finis
Considérons l’équation de Stokes en coordonnées cartésiennes 2D pour les écoulement incompres-
sible écrite sous la forme suivante :


∂U 2 2


 ∂t
+ ∂P
∂x
− ( ∂∂xU2 + ∂∂yU2 ) = fx

∂V 2 2
∂t
+ ∂P
∂y
− ( ∂∂xV2 + ∂∂yV2 ) = fy (2.3)


 1 ∂P + ∂U + ∂V = 0

β ∂t ∂x ∂y

Où f=(fx ,fy ) représente le terme source

Pour de simplification des calculs considérons fx et fy est égale zéro.

∂U 2 2
 + ∂P − ( ∂∂xU2 + ∂∂yU2 ) = 0
 ∂t ∂x


∂V 2 2
∂t
+ ∂P
∂y
− ( ∂∂xV2 + ∂∂yV2 ) = 0 (2.4)


 1 ∂P + ∂U + ∂V = 0

β ∂t ∂x ∂y

Ces équations peuvent s’écrire sous forment des opérateurs comme suit :

∂U
 + ∇P − ∆U = 0 sur Ω
 ∂t


1 ∂P (2.5)
β ∂t
+ divU = 0 sur Ω



U⃗ = (U, V )

Nous présentons dans cette section la discrétisation de ces opérateurs pour les deux composantes
de vitesse U et V . Les propriétés physiques de l’écoulement sont supposées uniformes. La discrétisa-
tion est formulée à partir de schémas centrés du second ordre suivant une approche volumes finis, en
se référant aux éléments de la métrique associés aux différents maillages décrits précédemment (les
maillage étant supposés réguliers).
L’intégrale appliquée au système (2.4) sur le volume fini correspondant au nœud P, nous donne :
R R
∂U ∂P 2 ∂2U
( ∂∂xU2 +
RR RR

∂t
dxdy + ∆V ∂x
dxdy − ∂y 2
)dxdy =0
R R∆V ∆V


∂V ∂P 2 ∂2V
( ∂∂xV2 +
RR RR
∆V ∂t
dxdy + ∆V ∂y
dxdy − ∆V ∂y 2
)dxdy = 0 (2.6)


1 ∂P
R R ( ∂U ∂V
 RR
∆V β ∂t
dxdy + ∆V ∂x
+ ∂y
)dxdy = 0

Z Z
∂U ∂U ∂U
dxdy = ∆Vui,j = ∆xi ∆yj (2.7)
∆V ∂t ∂t ∂t

Z Z
∂P
dxdy = (Pi,j+1 − Pi,j )∆yj (2.8)
∆V ∂x

∂ 2U Ui,j+1 − 2Ui,j + Ui,j−1

Z Z
2
dxdy = ∆yj (2.9)
∆V ∂x ∆xi

15
 ∂ 2U Ui+1,j − 2Ui,j + Ui−1,j
Z Z
2
dxdy = ∆xi (2.10)
∆V ∂y ∆yj

– Flux visqueux

ϕi,j = Dϕi,j + ax,p .ϕi+1,j + ax,m .ϕi−1,j + ay,p .ϕi,j+1 + ax,m .ϕi,j−1 (2.11)

ϕ = U, V (2.12)

Où
Les indices p et m désignent respectivement les interfaces de maille supérieures et inférieures où
sont calculés les flux pour chacune des directions x et y. D et a sont les coefficients diagonaux et
extra-diagonaux résultant du calcul de bilan de flux.

Pour U
1
ax,p =
∆x2

1
ax,m =
∆x2

1
ay,p =
∆y 2

1
ay,m =
∆y 2
Pour V
1
ax,p =
∆x2

1
ax,m =
∆x2

1
ay,p =
∆y 2

1
ay,m =
∆y 2

D = −(ax,p + ax,m + ay,p + ay,m ) (2.13)

- Pour l’équation de continuité suivant P

Z Z
1 ∂P 1 ∂P
dxdy = ∆x∆y (2.14)
∆V β ∂t β ∂t

16
 - Pour l’équation de continuité suivant U
Z Z
∂U
dxdy = (Ui,j − Ui,j−1 )∆y (2.15)
∆V ∂x

- Pour l’équation de continuité suivant V

Z Z
∂V
dxdy = (Vi−1,j − Vi,j )∆x (2.16)
∆V ∂y

2.3.4 Discrétisation temporelle

Le schéma explicite d’Euler

Le schéma explicite d’Euler est une méthode numérique consiste à discrétiser les équations diffé-
rentielles dans le temps en utilisant une méthode de différences finies, c’est-à-dire en approximant les
dérivées temporelles. Cependant, comme pour toute méthode numérique, le schéma explicite d’Euler
peut présenter des limitations en termes de stabilité et de précision. Des méthodes alternatives telles
que les schémas implicites ou semi-implicites peuvent être utilisées pour surmonter ces limitations
dans certains cas.
Dans le cas de notre travail le terme ∂U ∂t
on l’approxime en temps par une différence finie et en
considérant le schéma explicite d’Euler et ce qui nous donne le terme suivant :
n+1 n
∂U Ui,j − Ui,j
(tn , xi , yj ) ≃ (2.17)
∂t ∆t

n+1 2∆t 2∆t n ∆t n ∆t n ∆t n ∆t n ∆t n n

Ui,j = (1 − − )Ui,j + Ui,j+1 + Ui,j−1 + Ui+1,j + Ui−1,j − (P − Pi,j )
∆x 2 ∆y 2 ∆x 2 ∆x 2 ∆y 2 ∆y 2 ∆x i,j+1
(2.18)

n+1 2∆t 2∆t n ∆t n ∆t n ∆t n ∆t n ∆t n n

Vi,j = (1 − − )Vi,j + V i,j+1 + Vi,j−1 + Vi+1,j + Vi−1,j − (P − Pi+1,j )
∆x2 ∆y 2 ∆x2 ∆x2 ∆y 2 ∆y 2 ∆y i,j
(2.19)

n+1 n ∆tβ n n ∆tβ n n

Pi,j = Pi,j − (Ui,j − Ui,j−1 )− (Vi−1,j − Vi,j ) (2.20)
∆x ∆y

2.4 Convergence du schéma numérique issu de la MVF

2.4.1 Consistance
L’étude de consistance est une étape importante dans la validation des méthodes numériques, car
elle permet d’évaluer l’erreur de troncature, c’est-à-dire l’erreur introduite par la discrétisation des
équations continues en équations discrètes.

Définition 2.4.1. On appelle erreur de discrétisation en tn , Xi , Yj , la différence entre la solution

exacte en t, X, Y et les n-ème, i-ème et j-ème composantes de la solution donnée par le schéma

17
numérique :

⃗ (t, X, Y ) − U
⃗ei,j = U ⃗ (tn , Xi , Yj ) (2.21)

- L’erreur de consistance associée au flux l’expression :

n n n n
Rc = (Ui,j − Ui,j−1 ) + (Vi−1,j − Vi,j ) (2.22)

On dit que le flux numérique est consistant si

lim max|Rc | = 0 (2.23)

∆t→0,∆X→0,∆Y →0 κ∈τ

Où ∆t,représente le pas de temps ∆X, ∆Y représentent les pas des espaces, τ représente le
maillage admissible sur Ω et κ est son maille.

2.4.2 Etude de la stabilité par la méthode de Von neumann

La notion de stabilité s’applique à des schémas, pour lesquels on calcule des solutions de façon
itérative. Les calculs s’effectuent sur des ordinateurs avec une précision finie, et donc sont sujet à
des erreurs d’arrondis. Lors d’un calcul itératif, ces erreurs peuvent être amplifiées par le schéma
numérique. Le but de l’étude de stabilité est donc de déterminer quelle est l’amplification des erreurs
(ou perturbations) par le schéma.
L’analyse de stabilité au sens de Von Neumann consiste à étudier l’impact du schéma sur un mode
de Fourier. Pour que le mode de Fourier en question soit solution du problème continu, on remplace
les conditions de Dirichlet homogènes du problème par des conditions périodiques, et pour alléger les
notations, on considère les intervalles ]0, 2π[ × ]0, 2π[ comme intervalles d’étude en espaces.
Un schéma itératif est dit stable, si les perturbations de la solution numérique ne sont pas ampli-
fiées au cours des itérations.

Définition 2.4.2. Un schéma numérique est dit stable si

|G(k∆x, l∆y, β)| ≤ 1 (2.24)

Où G(k∆x, l∆y, β) est la matrice d’amplification

On pose

U (tn , x, y, ) = eikx · eily (2.25)

En interpolant sur les grilles spatiales des pas ∆x et ∆y, on obtient ainsi :

n n ik∆xi
Ui,j = Ûk,l e · eil∆yj (2.26)

18
 Et qui sera remplaçer dans les expressions suivantes :

n+1 2∆t 2∆t n ∆t n ∆t n ∆t n ∆t n ∆t n n

Ui,j = (1 − 2
− 2
)Ui,j + 2
Ui,j+1 + 2
U i,j−1 + 2
Ui+1,j + 2
Ui−1,j − (Pi,j+1 − Pi,j )
∆x ∆y ∆x ∆x ∆y ∆y ∆x

n+1 2∆t 2∆t n ∆t n ∆t n ∆t n ∆t n ∆t n n

Vi,j = (1 − − )Vi,j + Vi,j+1 + Vi,j−1 + Vi+1,j + Vi−1,j − (P − Pi+1,j )
∆x 2 ∆y 2 ∆x 2 ∆x 2 ∆y 2 ∆y 2 ∆y i,j

n+1 n ∆tβ n n ∆tβ n n

Pi,j = Pi,j − (Ui,j − Ui,j−1 )− (Vi−1,j − Vi,j )
∆x ∆y
Pour U.

n+1 n+1 ik∆xi

Ui,j = Ûk,l e · eil∆yj

n n ik∆x(i+1)
Ui+1,j = Ûk,l e · eil∆yj

n n ik∆x(i−1)
Ui−1,j = Ûk,l e · eil∆yj

n n ik∆xi
Ui,j+1 = Ûk,l e · eil∆y(j+1)

n n ik∆xi
Ui,j−1 = Ûk,l e · eil∆y(j−1)

Pour V

n+1 n+1 ik∆xi

Vi,j = V̂k,l e · eil∆yj

n n ik∆x(i+1)
Vi+1,j = V̂k,l e · eil∆yj

n n ik∆x(i−1)
Vi−1,j = V̂k,l e · eil∆yj

n n ik∆xi
Vi,j+1 = V̂k,l e · eil∆y(j+1)

n n ik∆xi
Vi,j−1 = V̂k,l e · eil∆y(j−1)

Pour P

n+1 n+1 ik∆xi

Pi,j = P̂k,l e · eil∆yj

n n ik∆x(i+1)
Pi+1,j = P̂k,l e · eil∆yj

19
 n n ik∆xi
Pi,j+1 = P̂k,l e · eil∆y(j+1)
2∆t 2∆t ∆t ∆t ∆t ∆t
On pose A = (1 − ∆x2
− ∆y 2
), B= ∆x2
,C = ∆y 2
, E= ∆x
et F = ∆y

n+1 ik∆xi il∆yj n ik∆xi il∆yj n ik∆xi il∆y(j+1) n ik∆xi il∆y(j−1) n ik∆x(i+1) il∆yj
Ûk,l e ·e = AÛk,l e ·e +B Ûk,l e ·e +B Ûk,l e ·e +C Ûk,l e ·e +

n ik∆x(i−1)
C Ûk,l e · eil∆yj − E(P̂k,l
n ik∆xi
e · eil∆y(j+1) − P̂k,l
n ik∆xi
e · eil∆yj )

n+1 ik∆xi il∆yj n ik∆xi il∆yj n ik∆xi il∆y(j+1) n ik∆xi il∆y(j−1) n ik∆x(i+1) il∆yj
V̂k,l e ·e = AV̂k,l e ·e +B V̂k,l e ·e +B V̂k,l e ·e +C V̂k,l e ·e +

n ik∆x(i−1)
C V̂k,l e · eil∆yj − F (P̂k,l
n ik∆xi
e · eil∆yj − P̂k,l
n ik∆x(i+1)
e · eil∆yj )

n+1 ik∆xi
P̂k,l e n ik∆xi
· eil∆yj = P̂k,l e n ik∆x
· eil∆yj − Eβ(Ûk,l e n ik∆xi
· eil∆yj − Ûk,l e · eil∆y(j−1) )−

n ik∆x(i−1)
F β(V̂k,l e · eil∆y − V̂k,l
n ik∆xi
e · eil∆Y j )

En simplifiant par eik∆Xi · eil∆Y j , on aurra :

n+1 n n il∆y n −il∆y n ik∆x n −ik∆x n il∆y n

Ûk,l = AÛk,l + B Ûk,l e + BUk,l e + C Ûk,l e + C Ûk,l e − E(P̂k,l e − P̂k,l )

n+1 n n il∆y n −il∆y n ik∆x n −ik∆x n n ik∆x

V̂k,l = AV̂k,l + B V̂k,l e + BVk,l e + C V̂k,l e + C V̂k,l e − F (P̂k,l − P̂k,l e )

n+1 n n n −il∆y n −ik∆x n

P̂k,l = P̂k,l − Eβ(Ûk,l − Ûk,l e ) − F β(V̂k,l e − V̂k,l )

En possant
eik∆x + e−ik∆x = 2 cos(k∆x)

et
eil∆y + e−il∆y = 2 cos(l∆y)

et en transformant le système sous la forme matricielle on a :

 n+1
   n 
Ûk,l A + 2C cos(k∆x) + 2B cos(l∆y) 0 −E(eil∆y − 1) Ûk,l
 n+1   ik∆x  n 
V̂k,l  =  0 A + 2C cos(k∆x) + 2B cos(l∆y) −F (1 − e )  V̂k,l 
n+1 −il∆y −ik∆x n
P̂k,l −Eβ(1 − e ) −F β(e − 1) 1 P̂k,l

20
 On pose la matrice d’amplification :
 
A + 2C cos(k∆x) + 2B cos(l∆y) 0 −E(eil∆y − 1)
G(k∆x, l∆y, β) =  0 A + 2C cos(k∆x) + 2B cos(l∆y) −F (1 − eik∆x )
 

−Eβ(1 − e−il∆y ) −F β(e−ik∆x − 1) 1

 
1 − 2C + 2C cos(k∆x) − 2B + 2B cos(l∆y) 0 −E(eil∆y − 1)
= 0 1 − 2C + 2C cos(k∆x) − 2B + 2B cos(l∆y) −F (1 − eik∆x )
 

−Eβ(1 − e−il∆y ) −F β(e−ik∆x − 1) 1

En passant par la transformation trigonométrique

x
1 − cos(x) = 2 sin2 ( )
2

1 − eix = 1 − cos(x) − i sin(x)

1 − e−ix = 1 − cos(x) + i sin(x)

e−ix − 1 = cos(x) − i sin(x) − 1

x
(1−eix )(e−ix −1) = (1−cos(x)−i sin(x))(cos(x)−i sin(x)−1) = −2+2 cos(x) = −2(1−cos(x)) = −4 sin2 ( )
2

x
(eix −1)(1−e−ix ) = (cos(x)+i sin(x)−1)(1−cos(x)+i sin(x)) = −2+2 cos(x) = −2(1−cos(x)) = −4 sin2 ( )
2

On cherche à calculer le polynome caratéristique de la matice G.

|G(k∆x, l∆y, β)| =

1 − 4C sin2 ( k∆x
2
) − 4B sin2 ( l∆y
2
)−λ 0 −E(eil∆y − 1)
= 0 1 − 4C sin2 ( k∆x
2
)) − 4B sin2 ( l∆y
2
) − λ −F (1 − eik∆x )
−Eβ(1 − e−il∆y ) −F β(e−ik∆x − 1) 1−λ

a−λ 0 −Eb1
= 0 a − λ −F c1
−Eβb2 −F βc2 1 − λ
et en calculant le polynome caratéristique de la matice on obtient l’expression :

PG (λ) = (λ − a)((λ − a)(1 − λ) − F 2 βc1 c2 )) − E 2 βb1 b2 (λ − a) (2.27)

= (λ − a)(λ2 − (1 + a)λ + a − F 2 βc1 c2 − E 2 βb1 b2 )

21
Avec
k∆x
c1 c2 = −4 sin2 ()
2
l∆y
b1 b2 = −4 sin2 ( )
2
Dire que PG (λ) = 0, =⇒ λ = a ou λ2 − (1 + a)λ + a − F 2 βc1 c2 − E 2 βb1 b2 = 0
Pour λ2 − (1 + a)λ + a − F 2 βc1 c2 − E 2 βb1 b2 = 0

∆ = (1 + a)2 − 4a + 4F 2 βc1 c2 + 4E 2 βb1 b2

= 1 − 2a + a2 + 4F 2 βc1 c2 + 4E 2 βb1 b2

= (a − 1)2 + 4F 2 βc1 c2 + 4E 2 βb1 b2

k∆x l∆y
= (a − 1)2 − 16β(F 2 sin2 ( ) + E 2 sin2 ( ))
2 2
Si PG (λ) ∈ R, alors

Premier cas :

Si 16β(F 2 sin2 ( k∆x

2
) + E 2 sin2 ( l∆y
2
)) = (a − 1)2 alors il existe une seule solution double λ1 =λ2 =
1+a
2
k∆x l∆y
λ1 = λ2 = 1 − 2C sin2 (
) − 4B sin2 ( )
2 2
On remarque que dans ce cas on a l’absence de β, ce que ne nous avantage pas sur l’analyse de
stabilité.

Deuxième cas :

Si 16β(F 2 sin2 ( k∆x

2
) + E 2 sin2 ( l∆y
2
)) > (a − 1)2 alors l’équation n’admet pas de solution dans R.

Troisième cas :

Si 16β(F 2 sin2 ( k∆x

2
) + E 2 sin2 ( l∆y
2
)) < (a − 1)2 alors il existe deux valeurs propres distincts λ1 et
λ2 dans R.

k∆x l∆y k∆x l∆y 2 k∆x l∆y

16β(F 2 sin2 ( )+E 2 sin2 ( )) < (a−1)2 = (−4C sin2 ( )−4B sin2 ( )) = 16(C sin2 ( )+B sin2 (
2 2 2 2 2 2

k∆x l∆y k∆x l∆y 2

β(F 2 sin2 ( ) + E 2 sin2 ( )) < (C sin2 ( ) + B sin2 ( ))
2 2 2 2

(C sin2 ( k∆x
2
) + B sin2 ( l∆y
2
))2
β<
F 2 sin2 ( k∆x
2
) + E 2 sin2 ( l∆y
2
)
q
1 + a − (a − 1)2 − 16β(F 2 sin2 ( k∆x
2
) + E 2 sin2 ( l∆y
2
))
λ1 =
2

22
 q
1 + a + (a − 1)2 − 16β(F 2 sin2 ( k∆x
2
) + E 2 sin2 ( l∆y
2
))
λ1 =
2

q
2 − 4C sin2 ( k∆x
2
) − 4B sin2 ( l∆y
2
) − (−4C sin2 ( k∆x
2
) − 4B sin2 ( l∆y
2
))2 − 16β(F 2 sin2 ( k∆x
2
) + E 2 sin2 ( l∆y
2
))
λ1 =
2

Alors les trois valeurs propres sont données par les expressions suivantes :

k∆x l∆y
λ0 = 1 − 4C sin2 ( ) − 4B sin2 ( )
2 2

r
k∆x2 2 l∆y k∆x l∆y 2 k∆x l∆y
λ1 = 1−2C sin ( )−2B sin ( )−2 (C sin2 ( ) + B sin2 ( )) − β(F 2 sin2 ( ) + E 2 sin2 ( ))
2 2 2 2 2 2

r
k∆x2 2 l∆y k∆x l∆y 2 k∆x l∆y
λ2 = 1−2C sin ( )−2B sin ( )+2 (C sin2 ( ) + B sin2 ( )) − β(F 2 sin2 ( ) + E 2 sin2 ( ))
2 2 2 2 2 2

Si
(C sin2 ( k∆x
2
) + B sin2 ( l∆y
2
))2
β≤
F 2 sin2 ( k∆x
2
) + E 2 sin2 ( l∆y
2
)
Pour
∆x = ∆y = h

∆t kh lh
λ0 = 1 − 4 2
(sin2 ( ) + sin2 ( ))
h 2 2

s
2∆t kh lh 2∆t (sin2 ( kh ) + sin2 ( lh2 ))2 kh lh
λ1 = 1 − 2 (sin2 ( ) + sin2 ( )) − 2
2
− β(sin2 ( ) + sin2 ( ))
h 2 2 h h 2 2

s
2∆t kh lh 2∆t (sin2 ( kh ) + sin2 ( lh2 ))2 kh lh
λ2 = 1 − 2
(sin2 ( ) + sin2 ( )) + 2
2
− β(sin2 ( ) + sin2 ( ))
h 2 2 h h 2 2

λ1 et λ2 sont réels si :
sin2 ( kh
2
) + sin2 ( lh2 )
β≤
h2
On soit que, ∀ kh, lh ∈ ]0, 2π[ × ]0, 2π[

kh
0 ≤ sin2 ( )≤1
2

23
, alors
kh lh
0 ≤ sin2 ( ) + sin2 ( ) ≤ 2
2 2

4∆t kh lh 8∆t
0≤ 2
(sin2 ( ) + sin2 ( )) ≤ 2
h 2 2 h
Alors ∀ kh, lh ∈ ]0, 2π[ × ]0, 2π[ on a |λ0 | est strictement inférieur à 1 si et seulement si 0< 8∆t
h2
≤ 1,
h2
i.e 0<∆t ≤ 8 .
Pour |λ1 | est strictement inférieur à 1 si et seulement si
s
2∆t 2 kh 2 lh 2∆t (sin2 ( kh
2
) + sin2 ( lh2 ))2 2 kh 2 lh
| (sin ( ) + sin ( )) − − β(sin ( ) + sin ( ))| ≤ 1
h2 2 2 h h2 2 2

sin2 ( kh
2
) + sin2 ( lh2 ) 2
0≤ 2
≤ 2
h h

2
0≤β≤ (2.28)
h2

Or ∀ β s
2∆t (sin2 ( kh ) + sin2 ( lh2 ))2 kh lh
0≤ 2
2
− β(sin2 ( ) + sin2 ( ))
h h 2 2
et s
2∆t (sin2 ( kh
2
) + sin2 ( lh2 ))2 2 kh 2 lh
− − β(sin ( ) + sin ( )) ≤ 0
h h2 2 2
et on a aussi
kh lh
0 ≤ sin2 ( ) + sin2 ( ) ≤ 2
2 2

2∆t kh lh 4∆t
0≤ 2
(sin2 ( ) + sin2 ( )) ≤ 2
h 2 2 h

s
2∆t kh lh 2∆t (sin2 ( kh ) + sin2 ( lh2 ))2 kh lh 4∆t
0 ≤ 2 (sin2 ( ) + sin2 ( )) − 2
2
− β(sin2 ( ) + sin2 ( )) ≤ 2
h 2 2 h h 2 2 h

Alors ∀ kh, lh ∈ ]0, 2π[ × ]0, 2π[ on a |λ1 | est strictement inférieur à 1 si et seulement si 0< 4∆t
h2
≤ 1,
i.e

h2
∆t ≤ (2.29)
4

Pour |λ2 | est strictement inférieur à 1 si et seulement si

s
2∆t kh lh 2∆t (sin2 ( kh ) + sin2 ( lh2 ))2 kh lh
| 2 (sin2 ( ) + sin2 ( )) + 2
2
− β(sin2 ( ) + sin2 ( ))| ≤ 1
h 2 2 h h 2 2

24
 sin2 ( kh
2
) + sin2 ( lh2 ) 2
0≤ 2
≤ 2
h h

2
0≤β≤
h2
Or ∀ β s
2∆t (sin2 ( kh
2
) + sin2 ( lh2 ))2 2 kh 2 lh
0≤ − β(sin ( ) + sin ( ))
h h2 2 2
et on a
kh lh
0 ≤ sin2 ( ) + sin2 ( ) ≤ 2
2 2

2∆t kh lh 4∆t
0≤ 2
(sin2 ( ) + sin2 ( )) ≤ 2
h 2 2 h

s
2∆t 2 kh 2 lh 2∆t (sin2 ( kh
2
) + sin2 ( lh2 ))2 2 kh 2 lh 4∆t
0≤ (sin ( ) + sin ( )) + − β(sin ( ) + sin ( )) ≤
h2 2 2 h h2 2 2 h2

Alors ∀ kh, lh ∈ ]0, 2π[ × ]0, 2π[ on a |λ1 | est strictement inférieur à 1 si et seulement si 0< 4∆t
h2
≤ 1,
i.e

h2
∆t ≤ (2.30)
4

D’où le schéma est conditionnellement stable sous β et ∆t

Si PG (λ) ∈ C, alors on obtient λ1 et λ2 comme suit :

s
2∆t kh lh 2∆t (sin2 ( kh ) + sin2 ( lh2 ))2 kh lh
λ1 = 1 − 2 (sin2 ( ) + sin2 ( )) − i − 2
2
+ β(sin2 ( ) + sin2 ( ))
h 2 2 h h 2 2

s
2∆t kh lh 2∆t (sin2 ( kh ) + sin2 ( lh2 ))2 kh lh
λ2 = 1 − 2
(sin2 ( ) + sin2 ( )) + i − 2
2
+ β(sin2 ( ) + sin2 ( ))
h 2 2 h h 2 2

|λ1 |2 =|λ2 |2 , alors

2∆t 2 kh 2 lh 4∆t2 2 kh 2 lh 2 4∆t2 kh lh

|λ1 |2 = (1 − 2
(sin ( ) + sin ( )))2
− 4
(sin ( ) + sin ( )) + 2
β(sin2 ( ) + sin2 ( ))
h 2 2 h 2 2 h 2 2

4∆t 2 kh 2 lh 4∆t2 kh lh
|λ1 |2 = 1 −2
(sin ( ) + sin ( )) + 2
β(sin2 ( ) + sin2 ( ))
h 2 2 h 2 2
On sait que ∀ kh, lh ∈ ]0, 2π[ × ]0, 2π[ on a :

4∆t 2 kh 2 lh 4∆t2 2 kh 2 lh 8∆t 8∆t2

0≤ (sin ( ) + sin ( )) + β(sin ( ) + sin ( )) ≤ + β
h2 2 2 h2 2 2 h2 h2

25
La condition de stabilité est vérifier si

8∆t 8∆t2
+ β ≤1
h2 h2

i.e

1
β≤ (h2 − 8∆t) (2.31)
8∆t2

et comme β est positive alors

h2
∆t ≤ (2.32)
8
2
Alors on peut conclure que le schéma explcite est stable si ∆t ≤ h8 et β ≤ h22 dans le cas le
2
polynôme caractéristique est réel et si ∆t ≤ h8 et β ≤ 8∆t
1 2
2 (h − 8∆t) dans le cas le polynôme

caractéristique est complexe.

2.4.3 Critère de convergence du schéma

La convergence d’un schéma numérique est une propriété théorique globale assurant que l’écart
(au sens d’une norme) entre la solution approchée et la solution exacte tend vers 0 lorsque le pas de
discrétisation tend vers 0 (ou lorsque chacun des pas globaux associés aux différentes directions tend
vers 0).
La solution approchée d’un schéma numérique reste peu crédible tant que sa convergence n’a pas
été montrée.
La notion de convergence d’un schéma exprime sa propriété à tendre vers la solution exacte du
problème traité pour des paramètres de maillage, tel que le nombre de nœuds tendant vers l’infini. Il
est possible de montrer que celle-ci est assurée sous deux conditions d’après la définition suivante :

Définition 2.4.3. Un schéma numérique est dit convergent s’il est à la fois stable et consistant.
Convergence = stabilité + consistance

26
Chapitre 3

Essais numériques

Dans ce dernier chapitre nous allons présenter dans un premier temps la partie implémentation
sur lequel nous allons présenter la formulation discrète, la procédure de résolution numérique et la
présentation du maillage. Dans un second temps nous allons exposer les résultats numériques sur
lesquels nous allons montrer les profils de la vitesse et effectuer des testes de stabilité.

3.1 Implémentation
3.1.1 Formulation discrète des équations de Stokes

n+1 2∆t 2∆t n ∆t n ∆t n ∆t n ∆t n ∆t n n

Ui,j = (1 − − )Ui,j + Ui,j+1 + Ui,j−1 + Ui+1,j + Ui−1,j − (P − Pi,j )
∆x 2 ∆y 2 ∆x 2 ∆x 2 ∆y 2 ∆y 2 ∆x i,j+1
(3.1)

n+1 2∆t 2∆t n ∆t n ∆t n ∆t n ∆t n ∆t n n

Vi,j = (1 − 2
− 2
)Vi,j + 2
Vi,j+1 + 2
Vi,j−1 + 2
Vi+1,j + 2
Vi−1,j − (Pi,j − Pi+1,j )
∆x ∆y ∆x ∆x ∆y ∆y ∆y
(3.2)

n+1 n ∆tβ n n ∆tβ n n

Pi,j = Pi,j − (Ui,j − Ui,j−1 )− (Vi−1,j − Vi,j ) (3.3)
∆x ∆y

3.1.2 Procédure de résolution numérique

27
 Figure 3.1 – Procédure de résolution numérique pour la méthode de compressibilité artificielle

3.1.3 Présentation du maillage

Figure 3.2 – Génération du maillage

28
 Figure 3.3 – Maillage en 2D

3.2 Résultats numériques

Cette partie à pour objectif de présenter les résultats numériques. On va d’abord présenter le
profil de chaque compossante de la vitesse avant de terminer par les tests de stabilité.

3.2.1 Profils des composantes de la vitesse Ux et Vx

Figure 3.4 – Le profil de chaque compossante de la vitesse suivant x

29
 On vient de traçer ci-dessous la figure 3.4 les profils des composantes de la vitesse Ux et Vx suivant
l’axe x pour un maillage de taille 11 × 11.

Ces deux composantes de la vitesse sont obtenus en prénant les paramétres suivants :

— β = 10
— ∆t = 0.001
— le critère de convergence 10−6
Les compossantes de la vitesse on convergé au bout de 5500 itération (figure 3.4) où on voit le résidu
décroit.

3.2.2 Tests de stabilité

Dans cette partie, on va aborder les tests de stabilité. On présente deux cas. Le premier cas
correspond à un test où le critère de stabilité est vérifié. Quant au second cas, le critère n’est pas
vérifié.

Cas 1 : stabilité de Ux et Vx

Figure 3.5 – Les composantes de la vitesse suivant x

En respectant le critère de stabilité, on constate au niveau de cette figure que le principe du

maximum est bien vérifié.

30
Ces deux profils sont obtenus pour les paramètres :

— β = 10
— ∆t = 0.001
— e = 10−6

Cas 2 : instabilité de Ux et Vx

Figure 3.6 – Les composantes de la vitesse suivant x

Cette discontinuité est due au fait que les paramères choisi ne respectent pas le critère de stabilité.
Ces deux profils sont obtenus pour les paramètres :

— β = 1000
— ∆t = 0.001
— e = 10−6

Remarque 3.2.1. moins le paramètre de compreessibilté artificielle est choisi pétit dans l’intervalle
considéré plus la stabilité est efficace.

Remarque 3.2.2. plus le nombre des points augmentent moins l’intervalle de pas de temps dimuni et
l’intervalle de paramètre de compreessibilté artificielle augmente.

31
Conclusion Générale

Dans cette étude nous avons commencé par la bibliographie sur lesquels nous avons énuméré quels
que anciennes travaux sur la résolution numérique par la méthode de compressibilité artificielle. Dans
le chapitre I nous avons présenté quels que équations aux dérivées partielles des différents phéno-
mènes sur lequel nous avons parlé de système de stokes. Dans le chapitre II nous avons étudié le
critère de stabilité au sens de Von Neumann par la méthode de Fourier, sur ce, nous avons obtenus
des conditions de stabilité sur β et ∆t.
Pour la discrétisation en temps la méthode d’Euler explicite a été appliquer et pour la discrétisation
en espace la méthode des volumes finis a été appliquer. Dans le chapitre III on a présenter les essaies
numériques pour valider les conditions de stabilité.
les résultats obtenus après simulation :
-Pour un nombre de points est égal à 11, la stabilité est confirmé pour un pas de temps compris
dans l’intervalle de [0 ; 001, 0 ; 00125] et pour un paramètre de compressibilité artificielle est égale à
]0, 250] ; -Pour un nombre de points est égal à 21 ; la stabilité est confirmée pour un pas de temps
compris dans l’intervalle de [0.0001, 0 ; 00031] et pour un paramètre de compressibilité artificielle
est égal à]0, 800] ; -Pour un nombre de points est égal à 31 ; la stabilité est confirmé pour un pas
de temps compris dans l’intervalle de [0 ; 0001, 0 ; 00014] et pour un paramètre de compressibilité
artificielle est égal à ]0, 1875].

Figure 3.7 – Le temps de calcul

Plus le nombre de mailles augmente plus le temps de calcul devient important.

Pour traiter le problème industriel il est nécessaire d’optimiser le temps de calcul. Par exemple
la méthode de calcul parallèle peut-être appliqué pour accélérer le temps de calcul.

32
Bibliographie

[1] Pietro Asinari, Taku Ohwada, Eliodoro Chiavazzo, and Antonio F Di Rienzo. Link-wise artificial
compressibility method. Journal of Computational Physics, 231(15) :5109–5143, 2012.
[2] BR Baliga and SV Patankar. A new finite-element formulation for convection-diffusion problems.
Numerical Heat Transfer, 3(4) :393–409, 1980.
[3] Alexandre Joel Chorin. A numerical method for solving incompressible viscous flow problems.
Journal of computational physics, 135(2) :118–125, 1997.
[4] Yann Fraigneau. Principes de base des méthodes numériques utilisées dans filecode sunfluidh
pour la simulation des écoulements incompressibles eta faible nombre de mach. Notes et docu-
ments. LIMSI, 2013.
[5] Djazia Hadji, Nabila Mhamdi, Meriem Khenig, and Halima Chelik. Résolution numérique des
équations aux dérivées partielles par la méthode des différences finies. PhD thesis, Université
Ibn Khaldoun-Tiaret-, 2018.
[6] Bernard Helffer. Introduction aux équations aux dérivées partielles. Université de Parie-sud,
mémoire de master, Version de janvier-mai, 2007.
[7] P Louda, K Kozel, and J Přı́hoda. Numerical solution of 2d and 3d viscous incompressible steady
and unsteady flows using artificial compressibility method. International journal for numerical
methods in fluids, 56(8) :1399–1407, 2008.
[8] Alexei Lozinski, Zoubida Mghazli, and Khallih Ould Ahmed Ould Blal. Méthode des éléments
finis multi-échelles pour le problème de stokes. Comptes Rendus Mathematique, 351(7-8) :271–
275, 2013.
[9] Obaid Ullah Mehmood, Norzieha Mustapha, and Sharidan Shafie. Unsteady two-dimensional
blood flow in porous artery with multi-irregular stenoses. Transport in porous media, 92 :259–
275, 2012.
[10] Taku Ohwada, Pietro Asinari, and Daisuke Yabusaki. Artificial compressibility method and
lattice boltzmann method : similarities and differences. Computers & Mathematics with Appli-
cations, 61(12) :3461–3474, 2011.
[11] HS Tang and Fotis Sotiropoulos. Fractional step artificial compressibility schemes for the uns-
teady incompressible navier–stokes equations. Computers & fluids, 36(5) :974–986, 2007.

33