1
% - 2nd order scheme presentation: [FouPip-2004,BouMadMetRaz-1997]
2
% - Fig. compare ux and uy cuts [GhiGhiShi-1982]
3
% - PCG convergence versus Re/delta_t and meshes
4
% - converence in time to stationary solution
7
% - adapt criteria ? see freefem...
8
\section{The Navier-Stokes problem}
9
\label{sec-navier-stokes}
10
\pbindex{Navier-Stokes}%
14
\cindex{benchmark!driven cavity flow}%
20
\subsection*{Formulation}
21
This longer example combines most fonctionalities presented
22
in the previous examples.
24
Let us consider the Navier-Stokes problem for the driven cavity
25
in $\Omega = ]0,1[^d$, $d = 2,3$.
26
Let $Re>0$ be the Reynolds number,
27
and $T>0$ a final time.
30
\ \ \ $(NS)$: {\it find ${\bf u}=(u_0,\ldots,u_{d-1})$
31
and $p$ defined in $\Omega \times ]0,T[$ such that:}
34
Re\left( \Frac{\partial {\bf u}}{\partial t} + {\bf u}.\nabla {\bf u} \right)
35
&-\ {\bf div}(2 D({\bf u}))
36
&+& \bnabla p &=& 0 \ {\rm in}\ \Omega \times ]0,T[, \\
37
&-\ {\rm div}\,{\bf u} && &=& 0 \ {\rm in}\ \Omega \times ]0,T[, \\
38
&{\bf u}(t\!=\!0) && &=& 0 \ {\rm in}\ \Omega \times \{0,T\}, \\
39
&{\bf u} && &=& (1,0) \ {\rm on}\ \Gamma_{\rm top} \times ]0,T[, \\
41
\ {\rm on}\ (\Gamma_{\rm left} \cup \Gamma_{\rm right} \cup \Gamma_{\rm bottom})
43
&\Frac{\partial u_0}{\partial {\bf n}} &=&
44
\Frac{\partial u_1}{\partial {\bf n}} &=&
45
u_2 = 0 \ {\rm on}\ (\Gamma_{\rm back} \cup \Gamma_{\rm front})\times ]0,T[
49
where $D({\bf u}) = (\bnabla {\bf u} + \bnabla {\bf u}^T)/2$.
50
This nonlinear problem is the natural extension of the linear Stokes
51
problem, as presented in paragraph~\ref{sec-navier-stokes},
52
page~\pageref{sec-navier-stokes}.
53
The boundaries are represented on
54
Fig.~\ref{fig-square-cube}, page~\pageref{fig-square-cube}.
56
\subsection*{Time approximation}
57
% ------------------------------
58
\cindex{characteristic method}
61
Let us consider the following backward second
63
$\phi\in C^2([0,T])$~:
65
\Frac{{\rm d}\phi}{{\rm d}t}(t)
67
\Frac{3\phi(t) - 4\phi(t-\Delta t) + \phi(t-2\Delta t)}
70
\mathcal{O}(\Delta t^2)
72
The problem is approximated by the following
73
second-order implicit Euler scheme:
79
-4{\bf u}^{n} \circ X^n
80
+ {\bf u}^{n-1} \circ X^{n-1}
84
&-\ {\bf div}(2 D({\bf u}^{n+1}))
85
&+& \bnabla p^{n+1} &=& 0 \ {\rm in}\ \Omega, \\
86
&-\ {\rm div}\,{\bf u}^{n+1} && &=& 0 \ {\rm in}\ \Omega, \\
87
&{\bf u}^{n+1} && &=& (1,0) \ {\rm on}\ \Gamma_{\rm top}, \\
88
&{\bf u}^{n+1} && &=& 0
89
\ {\rm on}\ \Gamma_{\rm left} \cup \Gamma_{\rm right} \cup \Gamma_{\rm bottom}, \\
90
&\Frac{\partial u_0^{n+1}}{\partial {\bf n}} =
91
\Frac{\partial u_1^{n+1}}{\partial {\bf n}}
95
\ {\rm on}\ \Gamma_{\rm back} \cup \Gamma_{\rm front} \ \mbox{ when } d=3,
98
where, following~\cite{BouMadMetRaz-1997,FouPip-2004}:
100
X^n(x) &=& x - \Delta t\, {\bf u}^*(x) \\
101
X^{n-1}(x) &=& x - 2\Delta t\, {\bf u}^*(x) \\
102
{\bf u}^* &=& 2{\bf u}^n - {\bf u}^{n-1}
104
It is a second order extension of
105
the method previously introduced in
106
paragraph~\ref{characteristic-method}
107
page~\pageref{characteristic-method}.
108
The scheme defines a second order recurence
109
for the sequence $({\bf u}^{n})_{n\geq -1}$,
110
that starts with ${\bf u}^{-1}={\bf u}^{0}=0$.
112
\subsection*{Variationnal formulation}
113
% ------------------------------------
114
The variationnal formulation of this problem expresses:
116
\ \ \ $(NS)_{\Delta t}$: \ {\it find ${\bf u}^{n+1}\in {\bf V}(1)$ and $p^{n+1} \in L^2_0(\Omega)$ such that:}
118
\begin{array}{lcccll}
119
a({\bf u}^{n+1},{\bf v}) &+& b({\bf v}, p^{n+1})
120
&=& m({\bf f}^{n},{\bf v}),
121
& \forall {\bf v}\in {\bf V}(0), \\
123
b({\bf u}^{n+1},q) && &=& 0, & \forall q \in L^2_0(\Omega),
132
4\, {\bf u}^{n}\circ X^n
134
{\bf u}^{n-1}\circ X^n
139
a({\bf u},{\bf v}) &=& \Frac{3 Re}{2\Delta t} \int_\Omega {\bf u} . {\bf v} \, dx
140
+ \int_\Omega 2 D({\bf u}) : D({\bf v}) \, dx
142
and $b(.,.)$ and ${\bf V}(\alpha)$
143
was already introduced in paragraph~\ref{sec-stokes},
144
page~\pageref{sec-stokes}, while studying the Stokes problem.
146
\subsection*{Space approximation}
147
% ------------------------------------
148
The Taylor-Hood~\cite{hood-taylor-73} finite element approximation
149
of this generalised Stokes problem was also considered
150
in paragraph~\ref{sec-stokes}, page~\pageref{sec-stokes}.
151
We introduce a mesh ${\cal T}_h$ of $\Omega$
152
and the finite dimensional spaces
153
${\bf X}_h$, ${\bf V}_h(\alpha)$ and $Q_h$.
154
The approximate problem writes:
156
\ \ \ $(NS)_{\Delta t,h}$: {\it find ${\bf u}_h^{n+1} \in {\bf V}_h(1)$ and $p^{n+1}\in Q_h$ such that:}
158
\begin{array}{lcccll}
159
a({\bf u}_h^{n+1},{\bf v}) &+& b({\bf v}, p_h^{n+1}) &=&
160
m({\bf f}_h^{n},{\bf v}),
161
& \forall {\bf v}\in {\bf V}_h(0), \\
162
b({\bf u}_h^{n+1},q) && &=& 0, & \forall q \in Q_h.
164
\label{eq-fvh-navier-stokes}
172
4\, {\bf u}_h^{n}\circ X^n
174
{\bf u}_h^{n-1}\circ X^n
177
The problem reduces to a sequence resolution
178
of a generalized Stokes problems.
180
\subsection*{File \file{navier\_stokes\_solve.icc}}
181
\exindex{navier\_stokes\_solve.icc}
182
\myexample{navier_stokes_solve.icc}
184
\subsection*{Comments}
185
% ----------------------------------
186
The {\tt navier\_stokes\_solve} function is similar to the
187
\file{stokes\_cavity.cc}.
188
It solves here a generalised Stokes problem
189
and manages a right-hand side ${\bf f}_h$:
190
\begin{lstlisting}[numbers=none,frame=none]
191
characteristic X1 ( -delta_t*uh_star);
192
characteristic X2 (-2.0*delta_t*uh_star);
193
field l1h = riesz(Xh, compose(uh1,X1), qopt);
194
field l2h = riesz(Xh, compose(uh2,X2), qopt);
195
field lh = l0h + (Re/delta_t)*(2*l1h - 0.5*l2h);
197
\exindex{convect.cc}%
198
This last computation is similar to those done in the
199
\file{convect.cc} example.
200
\clindex{solver\_abtb}%
201
The generalized Stokes problem is solved by the \code{solver\_abtb} class.
202
The stopping criterion is related to the stationnary solution or the maximal
205
%\cindex{preconditioner!Cahouet-Chabart}%
206
% The generalized Stokes problem is solved by the \code{solver\_abtb} class.
207
% The preconditioner is here the Cahouet and Chabart
208
% one~\cite{CahCha-1988}.
209
% As showed in~\cite{KobOls-2000-precond-uzawa-stokes-gen},
210
% the number of iterations need by the conjugate gradient algorithm
211
% to reach a given precision is then independent of the mesh size.
212
% This preconditioner leads to the
213
% resolution of the follwoing subproblem:
215
% (M_h^{-1} + \lambda A_h^{-1})q_h = r_h
219
% where $\lambda=Re/\Delta t$,
220
% $r_h\in Q_h$ is given and $M$ and $A$ are respectively the
221
% the mass matrix and the discrete Poisson operator
222
% with Neumann boundary condition.
223
% The resolution of this subproblem
224
% hqs been previously developed in section~\ref{sec-neumann-laplace},
225
% page~\pageref{sec-neumann-laplace}.
227
%\subsection*{File \file{cahouet-chabart.h}}
228
% \exindex{cahouet-chabart.h}%
229
% \exindex{build-poisson-neumann.h}%
230
% \myexample{cahouet-chabart.h}
232
\subsection*{File \file{navier\_stokes\_cavity.cc}}
233
\exindex{navier\_stokes\_cavity.cc}%
234
\myexample{navier_stokes_cavity.cc}
237
\subsection*{File \file{navier\_stokes\_criterion.icc}}
238
\exindex{navier\_stokes\_criterion.icc}%
239
\myexample{navier_stokes_criterion.icc}
241
\subsection*{Comments}
242
% ----------------------------------
243
\cindex{mesh!adaptation!anisotropic}%
245
The code performs a computation by
246
using adaptive mesh refinement, in order to capture
248
The \code{adapt\_option\_type} declaration is used by \code{rheolef} to
249
send options to the mesh generator.
250
\exindex{cavity.icc}%
251
The code reuse the file \file{cavity.icc}
252
introduced page~\pageref{ex-cavity-icc}.
253
This file contains two functions that defines
254
boundary conditions associated to the cavity driven problem.
256
The {\tt criteria} function computes the
257
adaptive mesh refinement criteria:
259
c_h = (Re|{\bf u}_h|^2 + 2|D({\bf u}_h)|^2)^{1/2}
261
The {\tt criteria} function is similar to those presented
262
in the \file{embankment-adapt-2d.cc} example.
264
\subsection*{How to run the program}
265
% ----------------------------------
272
$Re=100$: 4804 elements, 2552 vertices &
273
$\psi_{max} = 9.5\times 10^{-6}, \psi_{min} = -0.103 $ \\
275
\includegraphics[width=7cm]{navier-stokes-Re=100-geo.pdf} &
276
\includegraphics[width=7cm]{navier-stokes-Re=100-psi.pdf} \\
278
$Re=400$: 5233 elements, 2768 vertices &
279
$\psi_{max} = 6.4\times 10^{-4}, \psi_{min} = -0.111 $ \\
281
\includegraphics[width=7cm]{navier-stokes-Re=400-geo.pdf} &
282
\includegraphics[width=7cm]{navier-stokes-Re=400-psi.pdf}
284
\caption{Meshes and stream functions associated to the solution of the
285
Navier-Stokes equations for $Re=100$ (top) and $Re=400$ (bottom).}
286
\label{fig-navier-stokes}
293
$Re=1000$: 5873 elements, 3106 vertices &
294
$\psi_{max} = 1.64\times 10^{-3}, \psi_{min} = -0.117 $ \\
296
\includegraphics[width=7cm]{navier-stokes-Re=1000-geo.pdf} &
297
\includegraphics[width=7cm]{navier-stokes-Re=1000-psi.pdf} \\
299
\caption{Meshes and stream functions associated to the solution of the
300
Navier-Stokes equations for $Re=1000$.}
301
\label{fig-navier-stokes-2}
303
The mesh loop adaptation is initiated from a \code{bamg} mesh
304
(see also appendix~\ref{sec-bamg}).
306
\exindex{square.bamgcad}%
307
\exindex{square.dmn}%
309
bamg -g square.bamgcad -o square.bamg
310
bamg2geo square.bamg square.dmn > square.geo
312
Then, compile and run the Navier-Stokes solver for the driven cavity
315
make navier_stokes_cavity
316
./navier_stokes_cavity square.geo 100
318
The program performs a computation with $Re=100$.
319
By default the time step is $\Delta t = 0.05$
320
and the computation loops for five mesh adaptations.
321
At each time step, the program prints an approximation of the time derivative,
322
and stops when a stationary solution is reached.
323
Then, we visualise the \file{square-5} adapted mesh and its associated solution:
324
\toindex{field!-scale}%
325
\toindex{field!-velocity}%
328
field square-5.field.gz -velocity -scale 4 -mayavi
330
Notice the \code{-scale} option that applies a multiplicative factor to
331
the arrow length when plotting.
332
\exindex{streamf\_cavity.cc}%
334
\toindex{field!-n-iso-negative}%
336
\cindex{stream function}%
337
The representation of the stream function writes:
340
zcat square-5.field.gz | ./streamf_cavity | field -bw -n-iso-negative 10 -
342
The programs \code{streamf\_cavity},
343
already introduced page~\pageref{ex-streamf-cavity-cc}, is here reused.
344
The last options of the \code{field} program draws isocontours
345
of the stream function using lines, as shown on Fig.~\ref{fig-navier-stokes}.
346
The zero isovalue separates the main flow from recirculations,
347
located in corners at the bottom of the cavity.
349
For $Re=400$ and $1000$ the computation writes:
351
./navier_stokes_cavity square.geo 400
352
./navier_stokes_cavity square.geo 1000
357
\hspace{-0cm} \includegraphics{navier-stokes-cut-u0.pdf}&
358
\hspace{-0cm} \includegraphics{navier-stokes-cut-u1.pdf}
360
\caption{Navier-Stokes: velocity profiles along lines passing
361
throught the center of the cavity, compared
362
with data from~\cite{GhiGhiShi-1982}:
363
(a) $u_0$ along the vertical line;
364
(b) $u_1$ along the horizontal line line.}
365
\label{fig-navier-stokes-cut}
367
The visualization of the cut of the horizontal velocity
368
along the vertical median line writes:
369
\toindex{field!-comp}%
370
\toindex{field!-cut}%
371
\toindex{field!-normal}%
372
\toindex{field!-origin}%
374
field square-5.field.gz -comp 0 -cut -normal -1 0 -origin 0.5 0
375
field square-5.field.gz -comp 1 -cut -normal 0 1 -origin 0 0.5
377
Fig.~\ref{fig-navier-stokes-cut} compare the cuts with data
378
from~\cite{GhiGhiShi-1982}, table~1 and~2
379
(see also~\cite{GupKal-2005}).
380
Observe that the solution is in good agreement with
381
these previous computations.
385
\begin{tabular}{lllllc}
387
$Re$ & & $x_c$ & $y_c$ & $-\psi_{\rm min}$ & $\psi_{\rm max}$ \\
390
& 0.613 & 0.738 & 0.103 & $9.5\times 10^{-6}$ \\
391
& Labeur and Wells \cite{LabWel-2007}
392
& 0.608 & 0.737 & 0.104 & - \\
393
& Donea and Huerta \cite{DonHue-2003}
394
& 0.62 & 0.74 & 0.103 & - \\
397
& 0.554 & 0.607 & 0.111 & $5.6 \times 10^{-4}$ \\
398
& Labeur and Wells \cite{LabWel-2007}
399
& 0.557 & 0.611 & 0.115 & - \\
400
& Donea and Huerta \cite{DonHue-2003}
401
& 0.568 & 0.606 & 0.110 & - \\
404
& 0.532 & 0.569 & 0.117 & $1.6\times 10^{-3}$ \\
405
& Labeur and Wells \cite{LabWel-2007}
406
& 0.524 & 0.560 & 0.121 & - \\
407
& Donea and Huerta \cite{DonHue-2003}
408
& 0.540 & 0.573 & 0.110 & - \\
412
\caption{Cavity flow: primary vortex position and stream
416
Finaly, table~\ref{tab-vortex} compares the
417
primary vortex position and its associated stream function value.
418
Notice also the good agreement with previous simulations.
419
The stream function extremal values are obtained by:
420
\toindex{field!-min}%
421
\toindex{field!-max}%
423
zcat square-5.field.gz | ./streamf_cavity | field -min -
424
zcat square-5.field.gz | ./streamf_cavity | field -max -
426
The maximal value has not yet been communicated to our knoledge and
427
is provided in table~\ref{tab-vortex} for cross validation purpose.
428
The small program that computes the primary vortex position is showed below.
431
zcat square-5.field.gz | ./streamf_cavity | ./vortex_position
434
\subsection*{File \file{vortex\_position.cc}}
435
\exindex{vortex\_position.cc}%
436
\myexample{vortex_position.cc}
438
For higher Reynolds number,
439
Shen~\cite{shen-1991} showed in 1991 that the flow converges
440
to a stationary state
441
for Reynolds numbers up to $10\,000$; for Reynolds numbers larger
442
than a critical value $10\,000<Re_1<10\,500$ and less than another critical
443
value $15\,000<Re_2<16\,000$, these authors
444
founded that the flow becomes periodic in time which
445
indicates a Hopf bifurcation; the flow loses time periodicity for
447
In 1998, Ould Salihi~\cite{Oul-1998} founded a loss of stationarity
448
between $10\,000$ and $20\,000$.
449
In 2002, Auteri et al.~\cite{AutParQua-2002} estimated
450
the critical value for the apparition of the first instability
451
to $Re_1\approx 8018$.
452
In 2005, Erturk et al.~\cite{ErtCorGok-2005} computed
453
steady driven cavity solutions up to $Re \leq 21\,000$.
454
Also in 2005, this result was infirmed by~\cite{GelLubOlsSta-2005}:
456
estimated $Re_1$ close to $8000$, in agreement with~\cite{AutParQua-2002}.
457
The 3D driven cavity has been investigated in~\cite{MinEth-1998}
458
by the method of characteristic (see also~\cite{MelLegDooWat-2011}
459
for 3D driven cavity computations).
460
In conclusion, the exploration of the driven cavity at large
461
Reynolds number is a fundamental challenge in computational
464
% ------------------------------------------------
465
% pb de stationarite discrete Re=1000
466
% ------------------------------------------------
467
% test d'arret |du/dt| < eps
468
% pour Re=1000 et les maillages adaptes 4 & 5,
469
% c'est periodique autour de 0.1 et ne % tends pas vers zero !
470
% => comment arreter les calculs autrement que avec n >= iter_max ?
472
% ----------------------------------------------------------------------
474
% ----------------------------------------------------------------------
475
% du/dt(n) = (11/6)u(n) - u(n-1) + (3/2)*u(n-2) - (1/3)*u(n-3)
476
% u* = 3*u(n-1) - 3*u(n-2) + u(n-3)
477
% REFS: BotForPasPeySab-2001
478
% O. Botella, M. Y. Forestier, R. Pasquetti, R. Peyret, C. Sabbah
479
% Chebyshev methods for the Navier-Stokes equations: algorithms and applications
480
% Nonlinear Analysis, Volume 47, Issue 6, August 2001, Pages 4157-4168
481
% ----------------------------------------------------------------------
482
% autres comparaisons:
483
% ----------------------------------------------------------------------
484
% The reader can compare the result shown on Fig.~\ref{fig-navier-stokes}
485
% to any available computation on this standard benchmark
486
% (e.g.~\cite[p.~526]{choe-kim-kim-kim-2001} for $Re=100$,
487
% \cite[p.~144]{pironneau-88} for $Re=500$).
488
% http://www.nag.co.uk/simulation/Fastflo/Documents/ToolBox/html/node22.htm
489
% Re=100-s05f02.gif (manque les recirculations)
492
% paper-30604.ps.gz : coupes u1 et u2
494
% http://cs.engr.uky.edu/~jzhang/colleague.html
499
% http://aimsciences.org/DCDS-B/Volume1-4/Choi.pdf
501
% page 32 : figures Re=1 et 100
503
% http://www.csc.fi/elmer/examples/cavity/
506
% http://www.math.temple.edu/~jxu/cfd.html
507
% Re=1000, 7500, 10000
508
% http://calvino.polito.it/~sberrone/Publications/Papers/Report_07_1999.ps.gz
512
% Hopf Bifurcation of the Unsteady Regularized Driven Cavity Flow, J. Comp. Phys. Vol. 95, 228-245 (1991)
514
% It is found that the flow
515
% converges to a stationary state
516
% for Reynolds numbers (Re) up to 10000; for Reynolds numbers larger
517
% than a critical value $10000<Re1<10500$ and less than another critical
518
% value $15000<Re2<16000$, the flow becomes periodic in time which
519
% indicates a Hopf bifurcation; the flow loses time periodicity for
524
% http://www.fem.gr.jp/fem/fluid/dcf/dcf7.html
525
% animation : Re=0..1000
added
removed
doc/pusrman/navier-stokes.tex
