- Email: [email protected]
Resolution of the Navier-Stokes equations in vorticity-velocity formulation via control techniques
PHYSlCA
Physica D 60 (1992) 185-193 North-Holland
Resolution of the Navier-Stokes equations in vorticity-velocity formulation via control techniques H a s s a n Manouzi D~partement de Math~matiques et de Statistique, Universit~ Laval, Quebec, Canada G I K 7P4
and L. Steven H o u t Department of Mathematics and Statistics, Simon Fraser University, Burnaby, B.C., Canada V5A 1S6
The velocity-vorticity formulation of the Navier-Stokes equations is turned into an optimal control problem subject to linear constraints which consist of Laplacian equations for the velocity and for the vorticity. Various optimization algorithms can be used for the numerical resolution of the optimal control problem in question; in particular, the use of a conjugate-gradient method is discussed.
onF
with
f ]g.ndF=0,
I. Introduction
u=g
(1.4)
We are concerned with an optimal control formulation, as well as its finite-element approximations, of the Navier-Stokes equations written in the vorticity-velocity variables (see e.g. refs. [2,3,16]):
n being the unit outward normal on F. Many engineering applications require the solution of the Navier-Stokes equations in two and three dimensions. Primitive variable formulation (i.e. the velocity-pressure formulation) is widely employed in numerical simulation of threedimensional flows, but difficulties occur because of the divergence free condition [5,8,17]. A popular alternative is the stream functionvorticity formulation, which is a variant of the vorticity-velocity formulation [2-4,16]. In this article, we will be mainly concerned with the vorticity-velocity formulation. This formulation is often used for the approximation of viscous incompressible plane or axially symmetric flows. The main advantage of the vorticity-velocity formulation is that the incompressibility constraint does not pose the same difficulty as in the primitive variable formulation. Furthermore, the vorticity field is resolved directly; the vorticity is an important physical quantity that characterizes
F
- v A t o + ~(to, u ) = c u r l f to=curlu divu=0
in 12,
inO, inO,
(1.1) (1.2) (1.3)
where @(to, u )
[ ( u . g r a d ) t o - (to.grad)u (u.grad)to i f d = 2 ,
if d = 3 ,
u stands for the velocity field, to the vorticity, and f a given body force. O C R d, d = 2 or 3, denotes a bounded open domain with a Lipschitz continuous boundary F. For boundary conditions, we choose i E-mail: [email protected].
0167-2789/92/$05.00 © 1992 - Elsevier Science Publishers B.V. All rights reserved
186
H. Manouzi, L.S. Hou / Navier-Stokes equations in vorticity-velocity formulation
certain futures of turbulence. However, the vorticity-velocity formulation also poses difficulty in the imposition of boundary conditions for the vorticity. In fact, the specification of the vorticity on the boundary requires knowledge of the velocity in the interior and on the boundary. This induces a coupling between the vorticity and velocity equations, and therefore, one is faced with the problem of solving a coupled system of partial differential equations. This imposes serious limitations for the rise of the vorticityvelocity formulation in calculating certain flows, e.g. flows in a multiply connected domain or flows whose mass flux rates are not known a priori. In practice one often employs a simple iterative scheme to handle the vorticity boundary conditions. Many different techniques for overcoming the difficulty in the imposition of vorticity boundary conditions have been employed in numerical methods. For example one can exploit the relationship between the vorticity and the stream function on the boundary (see, among others, refs. [9-11]) or use those schemes proposed by Quartapelle and Valz-Gris [12-14]. Thus in this article we attempt to devise a more efficient treatment of the boundary conditions that lead to decoupled equations of vorticity and of velocity. Our method is effected by the introduction of a control variable and the reformulation of the system of equations into an optimal control problem which involves only Laplacian equations for the vorticity and for the velocity; thus independent choices of approximating spaces for these variables are permitted. Consequently, we may use virtually all finite-element spaces that can be used for the approximation of the Laplacian equation. We introduce some notations that will be used in the sequel. Let HS(O) and H*(F) denote Sobolev spaces of order s defined on /2 and F, respectively; let HS(12) and H*(F) denote their vector-valued counter-parts. The inner product on H*(O) or H*(12) is denoted by ( . , . ) , ; that on H~(F) or H~(F) is denoted by ( ' , ' ) ~ . r . We use the convention: H ° = L 2, H ° = L 2, (. ," )0 =
(" , ' ) and ( - , . ) o . r = ( ' , ' ) r .
Also,
HI(o) = {vEHI(~Q): v =0
on F } ,
and
X(a) =
{r E Lz(~r]): div "r ~ L2(12)} ifd = 2 , (~r E (L2(a))d: div r E L2(~2)} ifd = 3.
( . , - ) denotes the duality pairing between a Banach space X and its dual space X*. The duality pairing between H1/2(F) and its dual H-1/2(F) is denoted by ( -,. }r. For details, see ref. [1]. We will make use of the following relation between the velocity and vorticity derived from (1.2): - A u = curl to
in O ,
(1.5)
which is a direct consequence of the vector identity curl(curl u) = grad(div u) - Au
(1.6)
and the continuity equation (1.3). The rest of the paper is organized as follows. In section 2, we will establish the equivalence of (1.1)-(1.4) and an optimal control problem subject to constraints that are derived from the vorticity-velocity formulation. In section 3, we describe our algorithm for solving the optimal control problem proposed in section 2. Note that in the sequel we will only consider the threedimensional problems.
2. An optimal control formulation of the Navier-Stokes equations
Our objective in this section is to reformulate (1.1)-(1.4) as an equivalent optimal control problem. Our idea is to introduce a new variable s = v grad to
187
H. Manouzi, L.S. Hou I Navier-Stokes equations in vorticity-velocity formulation
as a control and consider the optimization problem: seek u ~ H01(f/), m E H i ( f / ) and s E 2;(f/) that minimize ~ ( u , ca, s) = ½ f Is - v grad cal 2 dO O
+ ½ f Idiv s + curl f -
~(ca, u)l 2 dr/
/9
(2.1) subject to
-vAu=s T-s u=g
(2.2)
in 12,
on F ,
=s'n
(2.7)
and
(2.3) in f/
(2.4)
and 0ca
(1.4). Using (2.2)-(2.5) and the definition s = v grad ca, it is straightforward to show that ~T(u, ca, s) = 0, i.e. (u, ca, s) is an infimum of 3 in #Z,d. Second, suppose (u, to, s)EH~(12) x H i ( f / ) x ~7(f/) is an optimal solution of (2.6) which satisfy the constraint equations (2.2)(2.5), then 5"(u, ca, s) = 0 (since the first part of the proof guarantees that the infimum of f f in 0~,d is z e r o ) SO that s = v grad ca in 12
- v Aca + yca = - d i v s + 3' curl u
v On
Proof. First, let (u, ¢0) be a solution of (1.1)-
onF
(2.5)
where y is a positive constant parameter, By introducing the admissibility set
div s = @(ca, u) - curl f
in 12.
Thus (1.1) follows immediately from these two equations. (2.7) together with (2.4) implies (1.2). Taking the curl of (1.2) and using identity (1.6), we arrive at curl ca = - A u + grad(div u)
%~ = {(u, co, s) • n~(f/) x n ~ ( f / ) x Z(f/):
(2.8)
in 12.
(2.9)
(2.2) and (2.7) allow us to deduce
(u, ca, s) satisfies (2.2)-(2.5)}, vAu=vcurlca
in/2
the optimization problem (2,1)-(2.5) can be abbreviated as: seek (fi, &, g) E q/ad such that
so that by comparing with (2.9), we obtain
3 ( ~ , ,b, ~) -< ~ ( u , ca, s) V (u, ca, s) • %~.
grad(div u) = 0
in/2,
(2.6) or
We first establish the equivalence of the optimization problem (2.6) and the velocityvorticity formulation of the Navier-Stokes equations (1.1)-(1.4).
Proposition 2.1. If (u, ca) ~ H I ( O )
H1(12) is a solution of (1.1)-(1.4), then the triplet (u, ca, s) with s = v gradra constitutes a solution of (2.6) satisfying if(u, ca, s ) = 0 . Conversely, if (u, to, s) E H1(12) x Hi(f~) × Z(12) is an optimal solution of (2.6), then (u, ca) is necessarily a solution of (1.1)-(1.4). X
divu=C
inO,
where C is a constant. Now since J'a div u dO = fr u.ndr= fr g'ndF=0, we conclude C = 0, i.e. (1.2) holds. []
3. The Stokes problem and conjugate gradient methods Based on the optimal control formulation in-
H. Manouzi, L.S. Hou / Navier-Stokes equations in vorticity-velocity formulation
188
troduced in section 2, we may employ a gradient or conjugate gradient method to develop an iterative algorithm to solve this problem numerically. In this section, we will study the Stokes problem: - v Ato = curl f to = curl u divu=0 u=g
in 12,
(3.1)
in 12,
(3.2)
in 12,
(3.3)
on F .
(3.4)
Similar to section 2, we have the equivalence of the Stokes problem and this optimal control problem.
Proposition 3.1. (u, to, s) E H1(12) x Hi(O) x ~;(O) is a solution of the minimization problem (3.9) if and only if (u, to, s) ~ H1(12) X H1(12) X ,v(12) is a solution of the Stokes problem (3.1)(3.4) with s = v grad to.
Corollary 3.2. The minimization problem (3.9) has a unique solution (fi, tb, g ) ~ grad" Furthermore, we have
By dropping the nonlinearity, the functional (2.1) is now modified to
~c(~, ~,, ~) = 0.
5g(u, to, s) = ½f Is - v grad to[2 d12
Note that, for each given s, (3.5)-(3.8) determine a unique (u, to):= (u(s), to(s)). Now we set
D
+ -~ f Idiv s + curl f l 2 d12.
3C(s) = ~c(u(s), to(s), s) .
D
The constraint equations remain the same, i.e. -vAu=s T-s u =g
in 12,
on F ,
(3.6)
-vAto+yto=-divs+ycurlu
0to v~=s'n
(3.5)
in12,
onr.
(3.7) (3.8)
We now proceed as in section 2 to formulate an equivalent optimal control problem. We define the admissible set as follows:
Then the minimization problem (3.9) is equivalent to minimize :~(s) over 7(12). We intend to use a gradient method to solve this minimization problem. Since the functional ~(s) is quadratic and (u(s), to(s)) depends linearly on s, we immediately obtain the following results.
Proposition 3.3. The mapping s--~ ~(s) is twice continuously Frechet-differentiable with respect to s. The first Frechet-derivative of ~ is defined by: ( ~ ' ( s ) , t) = ; (s - v grad to) : t d12 D
% , = {(u, to, s) e H1(12) x x ' ( 1 2 ) x :r(12): (u, to, s) satisfies (3.5)-(3.8)}.
+ f (div s + curl f ) . div t dD D
Our optimization problem is to seek (ii, &, g) grad such that
-
~(¢i, tb, ~) -< 5g(u, to, s) '¢(u, to, s) E grad" (3.9)
+ ((t.n), tb)r V t ~ ( 1 2 ) ,
(div,). n
d12 + f (t
t)- y
D
(3.10)
H. Manouzi, L.S. Hou / Navier-Stokes equations in vorticity-velocity formulation
where y E H i ( / 2 ) and 4, E H i ( O ) are the solutions of
Ho(/2 ) x H i ( O ) are, respectively, the solution of -- u Az = qT -- q
- u A4, + 3'd} --- u(div s - v Am),
inO,
(3.11) - u A0 + 3'0 = - d i v q + 3' curl z
o4, = 0
onF,
On
189
(3.12)
00
V-~n = q . n
inO,
onF;
and and of - v Ay = --3" curl d}-
(3.13) --uAy=tT--t
in / 2 ,
:~'(s) can also be computed by: --
( ~ ' ( s ) , t)
;, A ~ + 3'~ = - d i v t + 3' curl y
in O ,
0¢, v -~-~n = t. n
onF.
= f (s - v grad ¢o):(t - z, grad ¢o,) d O Remark 1. (3.10) can be formally derived by taking variations in the Lagrangian functional
D
+ f (div s + f ) .div t d O
VtE ~(O),
O
(3.14)
where u, ~ Hol(/2) and oJ, ~ H 1(O) are the solutions of -uAu,=tx-t
~e(u, oJ, s, y, 4,) = ~r(u, o~, s) -
{ u(grad u, grad y) - (s r - s, y)}
-
{ v(grad to, grad 4,) + 3,(¢o, 4,) + (div s, 4,)
- y(curl u, 4,) - (s. n, 4,)r)
in/2,
- v A to, + y¢o, = - div t + 3' curl u ,
with respect s. (3.11)-(3.13) can be derived by taking variations with respect to u and to, respectively.
in/2
and
tgOJ,
v--=t.n
onF.
On
The second Frechet-derivative computed by:
~l"(s) may be
(~f"(s), ( q, t))
Remark 2. The calculation of ~ ' ( s ) is usually performed using (3.10). (3.14) will only be used to calculate ~"(s). Now we use the Polak-Ribiere version of the conjugate-gradient method to find the minimizer of ~r(. ) (see refs. [6,7]), which proceeds as follows. Step O: Initialization
= f (q - v grad O):(t - v grad O) d O D
sc°)~Z(O) +fdivq.divtd/2
Vq, t ~ ( / 2 ) ,
(3.15)
O
where
(y,O)~Hi(O)xHI(O)
and
(z, 0) E
given.
We then define r (°),d ( ° ) E ~ ( O )
by
190
H. Manouzi, L . S . H o u I Navier-Stokes equations in vorticity-velocity formulation
f b(ta, ,@) = v ] grad t~ : grad ~ d D
f r(°):td~+ f divr<°).divtdO n =
D
<~'C°)), t> vte ~;(a),
+
~,(,,,, 4,) v,,,, 4,
E .%(a).
d (o) = r (o) . Step O: Initialization Then for n >-0, assuming that s ("), r ('), d (n) are known, we obtain s ("+1), r ("+°, d ("÷1) by:
s (°) E ~ ( O ) given.
Step 1: Descent
Compute u~°) E V h and to (°) ~ vh:
Find % ~ R ,
a(u(O>, ,,~) = ((s(O))¢ _ s£O>, ,,,)
~ ( s (") s(n+l)
a . d (")) <- ~ ( s (") - ad(")),
= s(n) --
and(n)
Va E R,
V,,~ ~ V~o,
b(~(°) , Oh)= - ( d i v s (°) , Oh)
+ -/(curl u(~°), gJ~) VgJ, ~ v * ,
.
U(h° ) = g h on F . Step 2: Calculation o f the new descent direction
Compute y(O) E V~ and ~(°) E Vh:
Find r (n+l) E ~ ( D ) , such that
b(~(~ °), Oh) = - V(S(h°) -- V grad a)(h°), grad Oh) VOh E V h ,
f r("+x):t d a + f d i v r ( " + ~ ) ' d i v t d O O
a(y(h°), %) = y ( ~ (°), curl Zh) =(~'(s("+')),t)
VZh ~ V~.
VtE$(O), Compute r(h°) E I~ h:
(r(-+l), r (-+1) _ r(")) A. = (r(.), r(.) )
f r(h0): t h d f / + f div r(°). div t h dD Set
12
12
= f (S(h° ) - vgradt~(°)):th d a
d (n+l) = r (n+l) + A n d ( " ) .
/2
Step 3: Do n = n + 1, go to step 1.
- f ( d i v th)" *(h°) d ~ + f ( t T - h):y (°) dD D
We will employ a finite-element scheme to carry out the above algorithm. We choose finiteelement subspaces V h C H1(12), V h C H~(TI) and x h C,~(12). For each fixed h, the finiteelement analogue of the Polak-Ribiere algorithm is defined as follows. We first introduce the bilinear forms
D
+ f (div s(s°) + f ) . d i v t h d D
+ f(th. n ) . , 2 °) d r
Vth
Xh,
F
Set a(u, v) = v _f grad u :grad v d O
Vu, v ~ Hl(ff~)
d (o) = r {o) .
O
and
Then for n >-O, assuming that s <"), r ("), d (") arc
H . M a n o u z i , L . S . H o u / N a v i e r - S t o k e s equations in v o r t i c i t y - v e l o c i t y f o r m u l a t i o n
191
_(.+t) _(.+x) _ r(hn))
known, we obtain s (~+x), r (n+°, d (~+1) by:
rh
, rh
~n =
Step 1: Descent
(4n), 4 . ) )
Set d(n+l)
Choose a suitable a n E R, S ( n + l ) = S(n) ~ Otn d ( n ) .
_(n+l) = rh + A . d (~) .
Step 3: D o n = n + 1, go to step I.
Step 2: Calculation o f the new descent direction Compute u. ( -h+ a ) E
vh
and w. (hn + O E V h :
a(u~ "+'), Vh)
~,-("+'), v , )
= ((4"+") T b ( ('~ (n + 1) , * h ) =
-(div
Vv, ~V'o ,
ah'(n+l) ' * h ) '
v¢,,E
+ T(curl .".(n+l) ,¢,.) u(n+l) h
,
a(Y(h") , zh)
on F .
= gh
V h
We tested this (Stokes) algorithm with a 2D model problem on D = (0, 2) x (0, 1). The exact solution has the velocity u = ( - x , y)T and the Vorticity to = 0. Our numerical results conform very well with the exact solution. See figs. 1-3. For the Navier-Stokes equations, we may use the same scheme with the following replacements in corresponding steps: . (.) = ~ r. ,oh and O(h.) vh: Compute Yh
= ((S(hn + 1))T
Compute y~hn+l) E V h and 0 ~ "+1) E v h :
_(n+ 1)
-- a h
x
, Zh )
f
b(,~ "+'), oh) = _ v(s~n +1)
-[-
--
Yah ~ V h
,(n+l)
v grad u, a
] [(Zh ° grad)oJ(h") -- (oJ~") • grad)Zh] D
----.,1
, gl~tu 0h)
• [div Sth~) + f -
(U(h")" grad)co~ ")
,
+ (¢o~h")'grad)U(h~)]dD a(Y
~Zh ~
VZh Evho;
Vh" b(~(h") , Oh)
_(n+ 1)
Compute "h
~
~h
:
= --V((S(h"+1) -- U grad oa(h"+l)), grad Oh)
f - I( .h+ ~ ) , **t h d ~ + f div r~n+l). div g~
th
dg'l
+ f [(u(h")" grad)0h - (Oh" grad)u(h")l
D
=
f
n
. d~ (s~n+l) - v grad u,. (hn + l ) x):~a
• [div S(h") + f -
/-j
(U(h")" grad)CO~h")
+ (tOth")" grad)Uth " ) ] d D
- f (div th) • ~(n+t)v.h
dD
Compute rth") E ,~ h :
fr~h'):thd~+fdivrg').divt, da
+ f (t~ - th): y~hn÷" d a D
g2
+ f (div sl n+t) + f ) .div t h d/'~
F
Vt h
D
= f (S(h~) -- V g r a d n
n ~'h + f (t~.n).,~(n+~) dF
Vu h E V h .
,~h
-
f n
oJ(h~)) : (t.) da
(div t(~)) • d~(h") d a
192
H. Manouzi, L.S. Hou / Navier-Stokes equations in vorticity-velocity formulation
+ f (t x - t h): y(h~) dO o
+ f [div ~ ) + f - (.~"). grad)co~~) n
+ (~o~~) • grad)u~) 1 .div + f (t,. n). ,(~) dr
") d a
Vt(~~) e $h,
F
respectively. It is worth mentioning that these algorithms involve only positive-definite discrete problems and thus can be efficiently resolved by conjugate gradient methods. To summarize, our algorithm serves an approach to handle the difficulty of imposing boundary conditions in the vorticityvelocity formulation. It provides a systematic way of developing iterative algorithms for non-
Fig. 3. Approximate vorticity surface.
linear problems. It has the feature that all subproblems are well-conditioned positive definite (Laplacian type) problems for which efficient solvers have matured. When three-dimensional problems are considered, simple finite-element spaces can be used.
Acknowledgement The work of HM was supported by the Natural Science and Engineering Research Council of Canada under grant number OGP 0046545; that of LSH was supported by the Natural Science and Engineering Research Council of Canada under grant number OGP-0089763.
References Fig. 1. Approximate velocity vector field.
Fig. 2. Exact velocity vector field.
[1] R. Adams, Sobolev Spaces (Academic, New York, 1.975). [2] F.S,R. Dennis, D.B. Ingham and R.N. Cook, J. Cornput. Phys. 33 (1979) 325. [3] H.F. Fasel, Numerical Solution of the Complete Navier-Stokes Equations, Lecture Notes in Mathematics, Vol. 771 (Springer, Berlin, 1980). [4] M. Fortin and F. Thomasset, J. Comput. Phys. 31 (1979) 113. [5] V. Girault and P.-A. Raviart, Finite Element Methods for Navier-Stokes Equations (Springer, Berlin, 1986). [6] R. Glowinski, Numerical Methods for Nonlinear Variational Problems (Springer, New York, 1984). [7] R. Glowinski and O. Pironneau, SIAM Rev. 21 (1979) 167. [8] M. Gunzburger, Finite Element Methods for Incompressible Viscous Flows: A Guide to Theory, Practice and Algorithms (Academic, Boston, 1989).
H. Manouzi, L.S. Hou / Navier-Stokes equations in vorticity-velocity formulation [9] M. Gupta and R. Manohar, J. Comput. Phys. 31 (1979) 265. [10] M.D. Olson, Proc. Mc GilI-EIC Conf. on FEM in Civil Engg., Montreal, 1972. [11] S. Orszag and M. Israeli, Ann. Rev. Huid. Mech. 6 (1974) 281. [12] L. Quartapelle, J. Comput. Phys. 40 (1981) 453. [13] L. Quartapelle and M. Napolitano, Int. J. Num. Math. Huids 4 (1984) 109.
193
[14] L. Quartapelle and F. Valz-Gris, Int. J. Num. Math Huids 1 (1981) 129. [15] S.L. Smith and C.A. Brehbia, J. Comput. Phys. 17 (1975) 235. [16] C.G. Speziale, J. Comput. Phys. 73 (1987) 476. [17] R. Temam, Navier-Stokes Equations (North-Holland, Amsterdam, 1979).
Physica D 60 (1992) 185-193 North-Holland
Resolution of the Navier-Stokes equations in vorticity-velocity formulation via control techniques H a s s a n Manouzi D~partement de Math~matiques et de Statistique, Universit~ Laval, Quebec, Canada G I K 7P4
and L. Steven H o u t Department of Mathematics and Statistics, Simon Fraser University, Burnaby, B.C., Canada V5A 1S6
The velocity-vorticity formulation of the Navier-Stokes equations is turned into an optimal control problem subject to linear constraints which consist of Laplacian equations for the velocity and for the vorticity. Various optimization algorithms can be used for the numerical resolution of the optimal control problem in question; in particular, the use of a conjugate-gradient method is discussed.
onF
with
f ]g.ndF=0,
I. Introduction
u=g
(1.4)
We are concerned with an optimal control formulation, as well as its finite-element approximations, of the Navier-Stokes equations written in the vorticity-velocity variables (see e.g. refs. [2,3,16]):
n being the unit outward normal on F. Many engineering applications require the solution of the Navier-Stokes equations in two and three dimensions. Primitive variable formulation (i.e. the velocity-pressure formulation) is widely employed in numerical simulation of threedimensional flows, but difficulties occur because of the divergence free condition [5,8,17]. A popular alternative is the stream functionvorticity formulation, which is a variant of the vorticity-velocity formulation [2-4,16]. In this article, we will be mainly concerned with the vorticity-velocity formulation. This formulation is often used for the approximation of viscous incompressible plane or axially symmetric flows. The main advantage of the vorticity-velocity formulation is that the incompressibility constraint does not pose the same difficulty as in the primitive variable formulation. Furthermore, the vorticity field is resolved directly; the vorticity is an important physical quantity that characterizes
F
- v A t o + ~(to, u ) = c u r l f to=curlu divu=0
in 12,
inO, inO,
(1.1) (1.2) (1.3)
where @(to, u )
[ ( u . g r a d ) t o - (to.grad)u (u.grad)to i f d = 2 ,
if d = 3 ,
u stands for the velocity field, to the vorticity, and f a given body force. O C R d, d = 2 or 3, denotes a bounded open domain with a Lipschitz continuous boundary F. For boundary conditions, we choose i E-mail: [email protected].
0167-2789/92/$05.00 © 1992 - Elsevier Science Publishers B.V. All rights reserved
186
H. Manouzi, L.S. Hou / Navier-Stokes equations in vorticity-velocity formulation
certain futures of turbulence. However, the vorticity-velocity formulation also poses difficulty in the imposition of boundary conditions for the vorticity. In fact, the specification of the vorticity on the boundary requires knowledge of the velocity in the interior and on the boundary. This induces a coupling between the vorticity and velocity equations, and therefore, one is faced with the problem of solving a coupled system of partial differential equations. This imposes serious limitations for the rise of the vorticityvelocity formulation in calculating certain flows, e.g. flows in a multiply connected domain or flows whose mass flux rates are not known a priori. In practice one often employs a simple iterative scheme to handle the vorticity boundary conditions. Many different techniques for overcoming the difficulty in the imposition of vorticity boundary conditions have been employed in numerical methods. For example one can exploit the relationship between the vorticity and the stream function on the boundary (see, among others, refs. [9-11]) or use those schemes proposed by Quartapelle and Valz-Gris [12-14]. Thus in this article we attempt to devise a more efficient treatment of the boundary conditions that lead to decoupled equations of vorticity and of velocity. Our method is effected by the introduction of a control variable and the reformulation of the system of equations into an optimal control problem which involves only Laplacian equations for the vorticity and for the velocity; thus independent choices of approximating spaces for these variables are permitted. Consequently, we may use virtually all finite-element spaces that can be used for the approximation of the Laplacian equation. We introduce some notations that will be used in the sequel. Let HS(O) and H*(F) denote Sobolev spaces of order s defined on /2 and F, respectively; let HS(12) and H*(F) denote their vector-valued counter-parts. The inner product on H*(O) or H*(12) is denoted by ( . , . ) , ; that on H~(F) or H~(F) is denoted by ( ' , ' ) ~ . r . We use the convention: H ° = L 2, H ° = L 2, (. ," )0 =
(" , ' ) and ( - , . ) o . r = ( ' , ' ) r .
Also,
HI(o) = {vEHI(~Q): v =0
on F } ,
and
X(a) =
{r E Lz(~r]): div "r ~ L2(12)} ifd = 2 , (~r E (L2(a))d: div r E L2(~2)} ifd = 3.
( . , - ) denotes the duality pairing between a Banach space X and its dual space X*. The duality pairing between H1/2(F) and its dual H-1/2(F) is denoted by ( -,. }r. For details, see ref. [1]. We will make use of the following relation between the velocity and vorticity derived from (1.2): - A u = curl to
in O ,
(1.5)
which is a direct consequence of the vector identity curl(curl u) = grad(div u) - Au
(1.6)
and the continuity equation (1.3). The rest of the paper is organized as follows. In section 2, we will establish the equivalence of (1.1)-(1.4) and an optimal control problem subject to constraints that are derived from the vorticity-velocity formulation. In section 3, we describe our algorithm for solving the optimal control problem proposed in section 2. Note that in the sequel we will only consider the threedimensional problems.
2. An optimal control formulation of the Navier-Stokes equations
Our objective in this section is to reformulate (1.1)-(1.4) as an equivalent optimal control problem. Our idea is to introduce a new variable s = v grad to
187
H. Manouzi, L.S. Hou I Navier-Stokes equations in vorticity-velocity formulation
as a control and consider the optimization problem: seek u ~ H01(f/), m E H i ( f / ) and s E 2;(f/) that minimize ~ ( u , ca, s) = ½ f Is - v grad cal 2 dO O
+ ½ f Idiv s + curl f -
~(ca, u)l 2 dr/
/9
(2.1) subject to
-vAu=s T-s u=g
(2.2)
in 12,
on F ,
=s'n
(2.7)
and
(2.3) in f/
(2.4)
and 0ca
(1.4). Using (2.2)-(2.5) and the definition s = v grad ca, it is straightforward to show that ~T(u, ca, s) = 0, i.e. (u, ca, s) is an infimum of 3 in #Z,d. Second, suppose (u, to, s)EH~(12) x H i ( f / ) x ~7(f/) is an optimal solution of (2.6) which satisfy the constraint equations (2.2)(2.5), then 5"(u, ca, s) = 0 (since the first part of the proof guarantees that the infimum of f f in 0~,d is z e r o ) SO that s = v grad ca in 12
- v Aca + yca = - d i v s + 3' curl u
v On
Proof. First, let (u, ¢0) be a solution of (1.1)-
onF
(2.5)
where y is a positive constant parameter, By introducing the admissibility set
div s = @(ca, u) - curl f
in 12.
Thus (1.1) follows immediately from these two equations. (2.7) together with (2.4) implies (1.2). Taking the curl of (1.2) and using identity (1.6), we arrive at curl ca = - A u + grad(div u)
%~ = {(u, co, s) • n~(f/) x n ~ ( f / ) x Z(f/):
(2.8)
in 12.
(2.9)
(2.2) and (2.7) allow us to deduce
(u, ca, s) satisfies (2.2)-(2.5)}, vAu=vcurlca
in/2
the optimization problem (2,1)-(2.5) can be abbreviated as: seek (fi, &, g) E q/ad such that
so that by comparing with (2.9), we obtain
3 ( ~ , ,b, ~) -< ~ ( u , ca, s) V (u, ca, s) • %~.
grad(div u) = 0
in/2,
(2.6) or
We first establish the equivalence of the optimization problem (2.6) and the velocityvorticity formulation of the Navier-Stokes equations (1.1)-(1.4).
Proposition 2.1. If (u, ca) ~ H I ( O )
H1(12) is a solution of (1.1)-(1.4), then the triplet (u, ca, s) with s = v gradra constitutes a solution of (2.6) satisfying if(u, ca, s ) = 0 . Conversely, if (u, to, s) E H1(12) x Hi(f~) × Z(12) is an optimal solution of (2.6), then (u, ca) is necessarily a solution of (1.1)-(1.4). X
divu=C
inO,
where C is a constant. Now since J'a div u dO = fr u.ndr= fr g'ndF=0, we conclude C = 0, i.e. (1.2) holds. []
3. The Stokes problem and conjugate gradient methods Based on the optimal control formulation in-
H. Manouzi, L.S. Hou / Navier-Stokes equations in vorticity-velocity formulation
188
troduced in section 2, we may employ a gradient or conjugate gradient method to develop an iterative algorithm to solve this problem numerically. In this section, we will study the Stokes problem: - v Ato = curl f to = curl u divu=0 u=g
in 12,
(3.1)
in 12,
(3.2)
in 12,
(3.3)
on F .
(3.4)
Similar to section 2, we have the equivalence of the Stokes problem and this optimal control problem.
Proposition 3.1. (u, to, s) E H1(12) x Hi(O) x ~;(O) is a solution of the minimization problem (3.9) if and only if (u, to, s) ~ H1(12) X H1(12) X ,v(12) is a solution of the Stokes problem (3.1)(3.4) with s = v grad to.
Corollary 3.2. The minimization problem (3.9) has a unique solution (fi, tb, g ) ~ grad" Furthermore, we have
By dropping the nonlinearity, the functional (2.1) is now modified to
~c(~, ~,, ~) = 0.
5g(u, to, s) = ½f Is - v grad to[2 d12
Note that, for each given s, (3.5)-(3.8) determine a unique (u, to):= (u(s), to(s)). Now we set
D
+ -~ f Idiv s + curl f l 2 d12.
3C(s) = ~c(u(s), to(s), s) .
D
The constraint equations remain the same, i.e. -vAu=s T-s u =g
in 12,
on F ,
(3.6)
-vAto+yto=-divs+ycurlu
0to v~=s'n
(3.5)
in12,
onr.
(3.7) (3.8)
We now proceed as in section 2 to formulate an equivalent optimal control problem. We define the admissible set as follows:
Then the minimization problem (3.9) is equivalent to minimize :~(s) over 7(12). We intend to use a gradient method to solve this minimization problem. Since the functional ~(s) is quadratic and (u(s), to(s)) depends linearly on s, we immediately obtain the following results.
Proposition 3.3. The mapping s--~ ~(s) is twice continuously Frechet-differentiable with respect to s. The first Frechet-derivative of ~ is defined by: ( ~ ' ( s ) , t) = ; (s - v grad to) : t d12 D
% , = {(u, to, s) e H1(12) x x ' ( 1 2 ) x :r(12): (u, to, s) satisfies (3.5)-(3.8)}.
+ f (div s + curl f ) . div t dD D
Our optimization problem is to seek (ii, &, g) grad such that
-
~(¢i, tb, ~) -< 5g(u, to, s) '¢(u, to, s) E grad" (3.9)
+ ((t.n), tb)r V t ~ ( 1 2 ) ,
(div,). n
d12 + f (t
t)- y
D
(3.10)
H. Manouzi, L.S. Hou / Navier-Stokes equations in vorticity-velocity formulation
where y E H i ( / 2 ) and 4, E H i ( O ) are the solutions of
Ho(/2 ) x H i ( O ) are, respectively, the solution of -- u Az = qT -- q
- u A4, + 3'd} --- u(div s - v Am),
inO,
(3.11) - u A0 + 3'0 = - d i v q + 3' curl z
o4, = 0
onF,
On
189
(3.12)
00
V-~n = q . n
inO,
onF;
and and of - v Ay = --3" curl d}-
(3.13) --uAy=tT--t
in / 2 ,
:~'(s) can also be computed by: --
( ~ ' ( s ) , t)
;, A ~ + 3'~ = - d i v t + 3' curl y
in O ,
0¢, v -~-~n = t. n
onF.
= f (s - v grad ¢o):(t - z, grad ¢o,) d O Remark 1. (3.10) can be formally derived by taking variations in the Lagrangian functional
D
+ f (div s + f ) .div t d O
VtE ~(O),
O
(3.14)
where u, ~ Hol(/2) and oJ, ~ H 1(O) are the solutions of -uAu,=tx-t
~e(u, oJ, s, y, 4,) = ~r(u, o~, s) -
{ u(grad u, grad y) - (s r - s, y)}
-
{ v(grad to, grad 4,) + 3,(¢o, 4,) + (div s, 4,)
- y(curl u, 4,) - (s. n, 4,)r)
in/2,
- v A to, + y¢o, = - div t + 3' curl u ,
with respect s. (3.11)-(3.13) can be derived by taking variations with respect to u and to, respectively.
in/2
and
tgOJ,
v--=t.n
onF.
On
The second Frechet-derivative computed by:
~l"(s) may be
(~f"(s), ( q, t))
Remark 2. The calculation of ~ ' ( s ) is usually performed using (3.10). (3.14) will only be used to calculate ~"(s). Now we use the Polak-Ribiere version of the conjugate-gradient method to find the minimizer of ~r(. ) (see refs. [6,7]), which proceeds as follows. Step O: Initialization
= f (q - v grad O):(t - v grad O) d O D
sc°)~Z(O) +fdivq.divtd/2
Vq, t ~ ( / 2 ) ,
(3.15)
O
where
(y,O)~Hi(O)xHI(O)
and
(z, 0) E
given.
We then define r (°),d ( ° ) E ~ ( O )
by
190
H. Manouzi, L . S . H o u I Navier-Stokes equations in vorticity-velocity formulation
f b(ta, ,@) = v ] grad t~ : grad ~ d D
f r(°):td~+ f divr<°).divtdO n =
D
<~'C°)), t> vte ~;(a),
+
~,(,,,, 4,) v,,,, 4,
E .%(a).
d (o) = r (o) . Step O: Initialization Then for n >-0, assuming that s ("), r ('), d (n) are known, we obtain s ("+1), r ("+°, d ("÷1) by:
s (°) E ~ ( O ) given.
Step 1: Descent
Compute u~°) E V h and to (°) ~ vh:
Find % ~ R ,
a(u(O>, ,,~) = ((s(O))¢ _ s£O>, ,,,)
~ ( s (") s(n+l)
a . d (")) <- ~ ( s (") - ad(")),
= s(n) --
and(n)
Va E R,
V,,~ ~ V~o,
b(~(°) , Oh)= - ( d i v s (°) , Oh)
+ -/(curl u(~°), gJ~) VgJ, ~ v * ,
.
U(h° ) = g h on F . Step 2: Calculation o f the new descent direction
Compute y(O) E V~ and ~(°) E Vh:
Find r (n+l) E ~ ( D ) , such that
b(~(~ °), Oh) = - V(S(h°) -- V grad a)(h°), grad Oh) VOh E V h ,
f r("+x):t d a + f d i v r ( " + ~ ) ' d i v t d O O
a(y(h°), %) = y ( ~ (°), curl Zh) =(~'(s("+')),t)
VZh ~ V~.
VtE$(O), Compute r(h°) E I~ h:
(r(-+l), r (-+1) _ r(")) A. = (r(.), r(.) )
f r(h0): t h d f / + f div r(°). div t h dD Set
12
12
= f (S(h° ) - vgradt~(°)):th d a
d (n+l) = r (n+l) + A n d ( " ) .
/2
Step 3: Do n = n + 1, go to step 1.
- f ( d i v th)" *(h°) d ~ + f ( t T - h):y (°) dD D
We will employ a finite-element scheme to carry out the above algorithm. We choose finiteelement subspaces V h C H1(12), V h C H~(TI) and x h C,~(12). For each fixed h, the finiteelement analogue of the Polak-Ribiere algorithm is defined as follows. We first introduce the bilinear forms
D
+ f (div s(s°) + f ) . d i v t h d D
+ f(th. n ) . , 2 °) d r
Vth
Xh,
F
Set a(u, v) = v _f grad u :grad v d O
Vu, v ~ Hl(ff~)
d (o) = r {o) .
O
and
Then for n >-O, assuming that s <"), r ("), d (") arc
H . M a n o u z i , L . S . H o u / N a v i e r - S t o k e s equations in v o r t i c i t y - v e l o c i t y f o r m u l a t i o n
191
_(.+t) _(.+x) _ r(hn))
known, we obtain s (~+x), r (n+°, d (~+1) by:
rh
, rh
~n =
Step 1: Descent
(4n), 4 . ) )
Set d(n+l)
Choose a suitable a n E R, S ( n + l ) = S(n) ~ Otn d ( n ) .
_(n+l) = rh + A . d (~) .
Step 3: D o n = n + 1, go to step I.
Step 2: Calculation o f the new descent direction Compute u. ( -h+ a ) E
vh
and w. (hn + O E V h :
a(u~ "+'), Vh)
~,-("+'), v , )
= ((4"+") T b ( ('~ (n + 1) , * h ) =
-(div
Vv, ~V'o ,
ah'(n+l) ' * h ) '
v¢,,E
+ T(curl .".(n+l) ,¢,.) u(n+l) h
,
a(Y(h") , zh)
on F .
= gh
V h
We tested this (Stokes) algorithm with a 2D model problem on D = (0, 2) x (0, 1). The exact solution has the velocity u = ( - x , y)T and the Vorticity to = 0. Our numerical results conform very well with the exact solution. See figs. 1-3. For the Navier-Stokes equations, we may use the same scheme with the following replacements in corresponding steps: . (.) = ~ r. ,oh and O(h.) vh: Compute Yh
= ((S(hn + 1))T
Compute y~hn+l) E V h and 0 ~ "+1) E v h :
_(n+ 1)
-- a h
x
, Zh )
f
b(,~ "+'), oh) = _ v(s~n +1)
-[-
--
Yah ~ V h
,(n+l)
v grad u, a
] [(Zh ° grad)oJ(h") -- (oJ~") • grad)Zh] D
----.,1
, gl~tu 0h)
• [div Sth~) + f -
(U(h")" grad)co~ ")
,
+ (¢o~h")'grad)U(h~)]dD a(Y
~Zh ~
VZh Evho;
Vh" b(~(h") , Oh)
_(n+ 1)
Compute "h
~
~h
:
= --V((S(h"+1) -- U grad oa(h"+l)), grad Oh)
f - I( .h+ ~ ) , **t h d ~ + f div r~n+l). div g~
th
dg'l
+ f [(u(h")" grad)0h - (Oh" grad)u(h")l
D
=
f
n
. d~ (s~n+l) - v grad u,. (hn + l ) x):~a
• [div S(h") + f -
/-j
(U(h")" grad)CO~h")
+ (tOth")" grad)Uth " ) ] d D
- f (div th) • ~(n+t)v.h
dD
Compute rth") E ,~ h :
fr~h'):thd~+fdivrg').divt, da
+ f (t~ - th): y~hn÷" d a D
g2
+ f (div sl n+t) + f ) .div t h d/'~
F
Vt h
D
= f (S(h~) -- V g r a d n
n ~'h + f (t~.n).,~(n+~) dF
Vu h E V h .
,~h
-
f n
oJ(h~)) : (t.) da
(div t(~)) • d~(h") d a
192
H. Manouzi, L.S. Hou / Navier-Stokes equations in vorticity-velocity formulation
+ f (t x - t h): y(h~) dO o
+ f [div ~ ) + f - (.~"). grad)co~~) n
+ (~o~~) • grad)u~) 1 .div + f (t,. n). ,(~) dr
") d a
Vt(~~) e $h,
F
respectively. It is worth mentioning that these algorithms involve only positive-definite discrete problems and thus can be efficiently resolved by conjugate gradient methods. To summarize, our algorithm serves an approach to handle the difficulty of imposing boundary conditions in the vorticityvelocity formulation. It provides a systematic way of developing iterative algorithms for non-
Fig. 3. Approximate vorticity surface.
linear problems. It has the feature that all subproblems are well-conditioned positive definite (Laplacian type) problems for which efficient solvers have matured. When three-dimensional problems are considered, simple finite-element spaces can be used.
Acknowledgement The work of HM was supported by the Natural Science and Engineering Research Council of Canada under grant number OGP 0046545; that of LSH was supported by the Natural Science and Engineering Research Council of Canada under grant number OGP-0089763.
References Fig. 1. Approximate velocity vector field.
Fig. 2. Exact velocity vector field.
[1] R. Adams, Sobolev Spaces (Academic, New York, 1.975). [2] F.S,R. Dennis, D.B. Ingham and R.N. Cook, J. Cornput. Phys. 33 (1979) 325. [3] H.F. Fasel, Numerical Solution of the Complete Navier-Stokes Equations, Lecture Notes in Mathematics, Vol. 771 (Springer, Berlin, 1980). [4] M. Fortin and F. Thomasset, J. Comput. Phys. 31 (1979) 113. [5] V. Girault and P.-A. Raviart, Finite Element Methods for Navier-Stokes Equations (Springer, Berlin, 1986). [6] R. Glowinski, Numerical Methods for Nonlinear Variational Problems (Springer, New York, 1984). [7] R. Glowinski and O. Pironneau, SIAM Rev. 21 (1979) 167. [8] M. Gunzburger, Finite Element Methods for Incompressible Viscous Flows: A Guide to Theory, Practice and Algorithms (Academic, Boston, 1989).
H. Manouzi, L.S. Hou / Navier-Stokes equations in vorticity-velocity formulation [9] M. Gupta and R. Manohar, J. Comput. Phys. 31 (1979) 265. [10] M.D. Olson, Proc. Mc GilI-EIC Conf. on FEM in Civil Engg., Montreal, 1972. [11] S. Orszag and M. Israeli, Ann. Rev. Huid. Mech. 6 (1974) 281. [12] L. Quartapelle, J. Comput. Phys. 40 (1981) 453. [13] L. Quartapelle and M. Napolitano, Int. J. Num. Math. Huids 4 (1984) 109.
193
[14] L. Quartapelle and F. Valz-Gris, Int. J. Num. Math Huids 1 (1981) 129. [15] S.L. Smith and C.A. Brehbia, J. Comput. Phys. 17 (1975) 235. [16] C.G. Speziale, J. Comput. Phys. 73 (1987) 476. [17] R. Temam, Navier-Stokes Equations (North-Holland, Amsterdam, 1979).