- Email: [email protected]
Reduction of modes for the control of viscous flows
International Journal of Engineering Science 39 (2001) 177±200
www.elsevier.com/locate/ijengsci
Reduction of modes for the control of viscous ¯ows H.M. Park *, O.Y. Kim Department of Chemical Engineering, Sogang University, Shinsoo-Dong, Mapo-Gu, Seoul, South Korea Received 28 July 1999; received in revised form 17 September 1999; accepted 17 January 2000
Abstract This paper applies the Karhunen±Loeve Galerkin procedure to the solution of optimal control problems of viscous ¯uid ¯ow through a grooved channel. Through the Karhunen±Loeve Galerkin procedure, one can reduce the Navier±Stokes equation to a small number of ordinary dierential equations by limiting the solution space to the smallest linear subspace that is sucient to describe the observed phenomena. This technique is well suited for the problem of control or optimization, where one has to solve the governing equation repeatedly but one can also estimate the approximate solution space based on the range of control variables. In our previous work [H.M. Park, Y.D. Jang, Int. J. Engrg. Sci. 38 (2000) 785±805], we showed that the Karhunen±Loeve Galerkin procedure could solve an optimal control problem of one-dimensional Burgers equation very eciently. The present paper extends the applicability of the Karhunen±Loeve Galerkin procedure to the solution of optimal control problems of multi-dimensional Navier±Stokes equations. Ó 2000 Elsevier Science Ltd. All rights reserved. Keywords: Control of viscous ¯ows; Karhunen±Loeve Galerkin procedure
1. Introduction Flow control problems have been addressed by many scientists, mathematicians and engineers due to its frequent occurrence in various systems in science and engineering. The control of ¯uid ¯ow is crucial to achieve desired design objectives or to optimize the performance of the systems that exploit ¯uid motion. The development of a technique of ¯ow control, which is ecient as well as rigorous, is one of the central issues to be resolved before achieving optimal performance and safe operation of various systems with ¯uid ¯ow. But the quest for successful control of viscous ¯uid remains one of the most arduous challenges in the ¯uid and control sciences today. The *
Corresponding author. E-mail address: [email protected] (H.M. Park).
0020-7225/01/$ - see front matter Ó 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 0 - 7 2 2 5 ( 0 0 ) 0 0 0 2 4 - 0
178
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
Navier±Stokes equation, which is the governing equation for the viscous ¯uid ¯ow, has almost in®nite degrees of freedom (d.f.) with strong nonlinearity, thus causing dicult computational and theoretical problems. With the optimal control problems in ¯uid ¯ows, candidate optimal states and controls may belong to in®nite-dimensional function spaces and one may have to minimize a nonlinear functional of the state and control variables subject to nonlinear constraints that take the form of a system of partial dierential equations [2]. As a consequence the optimization problems are nonconvex, and this creates some technical diculties in studying the existence and uniqueness of an optimal control. Perhaps one of the most important obstacles to the application of rigorous control theories to viscous ¯uid ¯ows is the tremendously large requirement of computer time and memory. The actual implementation of control schemes requires repeated computation of the governing and related equations, and one cannot expect on-line control without reducing the number of nonlinear equations in the system. Therefore, one of the most crucial factors that determines the success of the numerical algorithm for the optimal control problems of ¯uid ¯ow is the employment of very ecient numerical method to solve the system equations. The traditional methods of the discretization of partial dierential equations usually involve functions that have very little to do with the dierential equations. For example, piecewise polynomials are used in ®nite element methods, grid functions are used in ®nite dierence methods, and trigonometric functions or Chebyshev polynomials are used in spectral methods. As a result they yield large systems of equations whose inversion consumes most of the computation time. Recently, several research groups [3,4] have used the Karhunen±Loeve decomposition as a procedure for computation of coherent structures (empirical eigenfunctions) in turbulence from experimental or numerical data of ¯ow ®elds. Usually these empirical eigenfunctions have been used to interpret the statistical characteristics of the turbulent ¯ow. Furthermore, employing these empirical eigenfunctions and introducing additional drastic approximations to the Navier±Stokes equations they could derive low-dimensional dynamic models for turbulent ¯ows [5] or transitional ¯ows [6]. These low-dimensional models [5,6] do not simulate the ¯ow ®eld exactly when the control variables like the Reynolds number varies, but are found to approximately reproduce some interesting characteristics of turbulence or transitional ¯ows such as intermittency and bifurcation sequences. Also there have been some attempts to apply the dynamic system methodology to these approximate low-dimensional models for the purpose of ¯ow control [7]. Because the original Karhunen±Loeve decomposition technique [4±6] is applicable only to stationary stochastic ®elds such as turbulence or oscillatory transitional ¯ows, it could not be employed in the rigorous control schemes of ¯uid ¯ows where the control variables vary. But recent works [1,8,9] have extended the applicability of the Karhunen±Loeve decomposition to the analysis of nonstationary, nonhomogeneous deterministic as well as stochastic ®elds and allowed the derivation of rigorous low-dimensional dynamic models that simulate the given systems almost exactly even when the control variables vary arbitrarily within certain bounds. This extension of the original Karhunen±Loeve decomposition is called the Karhunen±Loeve Galerkin procedure [1,8,9] and can be adopted in rigorous implementation of modern control schemes for viscous ¯uid ¯ows. In our previous work [1], we showed that the Karhunen±Loeve Galerkin procedure could solve an optimal control problem of one-dimensional Burgers equation very eciently. The present paper applies the Karhunen±Loeve Galerkin procedure to the solution of optimal control
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
179
problems of multi-dimensional Navier±Stokes equations. Since the reduction of modes in the multi-dimensional case is much larger than that in the one-dimensional case, the present technique is found to be more powerful when applied to the control problems of the Navier±Stokes equation than those of the Burgers equation. The details of the Karhunen±Loeve Galerkin procedure is well documented in our previous works [1,8,9] and reference therein. 2. System and governing equations We consider a viscous incompressible ¯ow through a grooved channel as shown in Fig. 1(b). The ¯ow is induced by the time-varying inlet velocity U t. The unsteady ¯ow ®eld of the system is obtained by solving the following incompressible Navier±Stokes equation with relevant initial and boundary conditions: ov v rv ÿrP lr2 v; 1 q ot r v 0;
2
where v u; w is the velocity vector, P the pressure, q the density ( 1) and l is the viscosity ( 1). The initial and boundary conditions are: t 0;
v v0 x;
x 0 inlet :
3
u U t a t 1 ÿ h z2 ;
w 0;
Fig. 1. (a) The trial control; (b) the target velocity ®eld.
4
180
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
x 1:0 outlet :
ou 0; ox
at all other boundaries :
w 0;
5
u w 0;
6
where h z z ÿ 0:875=0:125 for 0:75 6 z 6 1:0 and a t is the maximum of the inlet velocity of parabolic pro®le. This set of Eqs. (1)±(6) is solved accurately by using the ®nite volume method [10]. The problem at hand is that of controlling the system above to produce a state variable v x; T at the ®nal time T that is as close as possible to a target ¯ow ®eld vT x. We aim at achieving this goal by adjusting the maximum of the inlet velocity a t. Moreover, we try to achieve the goal at the minimum cost of the control a t. The problem may be speci®ed as the following optimization problem: min J a;
7
where 1 J a 2
Z
v x; T ÿ vT x dX 2 X 2
Z 0
T
2
a t ÿ a0 dt;
8
a0 the initial value of a and is a positive weighting factor which represents the expense of applying the control. The magnitude of is small if the control is cheap and large if it is expensive. To ensure that the target vT x is a reachable state, we determine vT x by integrating Eqs. (1)±(6) with a given function a t and setting vT x v x; t T . The trial control a t employed to obtain vT x is given by 80t 1 0 6 t 6 0:05; a t 9 ÿ80 t ÿ 0:1 1 0:05 6 t 6 0:1: The trial control (9) and the corresponding target ¯ow ®eld vT x are depicted in Fig. 1(a) and (b), respectively. 3. Low-dimensional dynamic model In this section, we employ the Karhunen±Loeve Galerkin procedure to derive a low-dimensional dynamic model which shall be used in the solution of the optimal control problems. The low-dimensional dynamic model obtained through the Karhunen±Loeve Galerkin procedure must predict the system behavior accurately for various trajectories of a t, which is the maximum of the inlet velocity, that include not only the optimal trajectory but also other trajectories appearing during the iterative minimization of the objective function equation (8) by means of the conjugate gradient method. The prerequisite to this requirement is that the empirical eigenfunctions employed in the Karhunen±Loeve Galerkin procedure must span the admissible solution space of the Navier±Stokes equation for these various trajectories of a t. According to the
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
181
Schmidt±Hilbert theory [1,11], the empirical eigenfunction /i x can be expressed linearly in terms of N snapshots: /i x
N X i1
ai j vj x;
10
where vj x is the jth snapshot. Therefore, once we obtain an ensemble of snapshots fvn xg which encompasses the admissible solution space of the Navier±Stokes equation for the above trajectories of a t, the resulting empirical eigenfunctions obtained through the Karhunen±Loeve decomposition shall span the same solution space, and consequently the low-dimensional model shall simulate the system accurately for those various trajectories of a t. In the present work, it is assumed that the admissible range of a t is between 1 and 5. Then for any arbitrary trajectory a t, we ®nd that the variation of a t during the time increment Dt is Da a t Dt ÿ a t, which can be positive or negative but its absolute value must be less than 4. This simple observation suggests the strategy to obtain snapshots. The snapshots of the ¯ow ®eld are obtained by solving Eqs. (1)±(6) when the maximum of the inlet velocity a t changes abruptly from 1 to 5 at t 0. The initial condition is taken as the steady ¯ow ®eld when the inlet velocity is 1, and the transient ¯ow ®elds are recorded at an appropriate time interval until a new steady state is reached. These recordings are served as snapshots to be used in the Karhunen±Loeve decomposition. In this way, we obtain the ®rst set of 500 snapshots. The second set of 500 snapshots are obtained by solving the governing equations (1)±(6) with the abrupt change of a t from 5 to 1 at t 0 when the initial condition is the steady ¯ow ®eld of a 5:0. In these computations, the ®nite volume method [10] with 80 80 grids is employed, which is found to resolve the ¯ow ®elds accurately. As explained by Park and Cho [8], it is necessary to make these snapshots satisfy homogeneous boundary conditions. For this purpose, we obtain the reference steady ¯ow ®eld vr ur ; wr by solving the following set of equations: qvr rvr ÿrPr lr2 vr ;
11
r vr 0
12
with the following boundary conditions: 2 inlet : ur ar 1 ÿ h z ; wr 0;
outlet :
our 0; ox
wr 0;
at all other boundaries :
ur wr 0;
13
14 15
where ar 3. We make the ®rst set of 500 snapshots satisfy homogeneous boundary condition by the following arrangement:
182
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
5 v x; t ÿ vr x ; 3 n
n 1; 2; . . . ; 500:
Similarly, the second set is transformed to satisfy homogeneous boundary condition as 1 n v x; t ÿ vr x ; n 1; 2; . . . ; 500: 3
16
17
To these 1000 snapshots satisfying homogeneous boundary conditions, we apply the Karhunen±Loeve decomposition technique to get empirical eigenfunctions in the order of their importance in characterizing the system. Fig. 2(a)±(d) shows the ®rst, second, third and fourth
Fig. 2. Some dominant empirical eigenfunctions with large eigenvalues: (a) ®rst eigenfunction k1 0:9413; (b) second eigenfunction k2 4:943 10ÿ2 ; (c) third eigenfunction k3 7:042 10ÿ3 ; (d) fourth eigenfunction k4 1:565 10ÿ3 .
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
183
empirical eigenfunctions with the corresponding normalized eigenvalues k1 0:9413, k2 4:943 10ÿ2 , k3 7:042 10ÿ3 , k4 1:565 10ÿ3 , respectively. Also shown in Fig. 3(a)± (d) are typical empirical eigenfunctions with smaller eigenvalues, i.e., the 5th±8th eigenfunctions with the corresponding normalized eigenvalues, k5 3:595 10ÿ4 , k6 1:935 10ÿ4 , k7 6:809 10ÿ5 , k8 1:835 10ÿ5 : In these ®gures, solid lines denote positive-valued streamlines and dotted lines denote negative-valued streamlines. As previously [1,8], the empirical eigenfunctions with large eigenvalues represent the large-scale structures of the velocity ®eld v. On the contrary, Fig. 3(a)±(d) reveals that empirical eigenfunctions with small eigenvalues represent the small-scale structures caused by system boundaries and nonlinearity of the governing equation.
Fig. 3. Empirical eigenfunctions with small eigenvalues: (a) ®fth eigenfunction k5 3:595 10ÿ4 ; (b) sixth eigenfunction k6 1:935 10ÿ4 ; (c) seventh eigenfunction k7 6:809 10ÿ5 ; (d) eighth eigenfunction k8 1:835 10ÿ5 .
184
as
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
The velocity ®eld v x; t is now represented as a linear combination of empirical eigenfunctions v
n X i1
ai t/i x b tvr ;
18
where /i is the ith empirical eigenfunction, ai t the corresponding spectral coecient, b t a t=ar where a t is the maximum of the inlet velocity and ar 3, and N is the total number of empirical eigenfunctions employed in the Karhunen±Loeve Galerkin procedure. Applying the Galerkin principle which enforces the residual R to be orthogonal to each of the basis functions, Z X
R /i dX 0
i 1; 2; . . . ; N;
19
where R
ov P l v rv r ÿ r2 v; ot q q
20
Eqs. (1)±(6) are reduced to the following set of ordinary dierential equations: Mk
N N X N N X X dak lX al Hkl ÿ al am Qklm ÿ b al Rkl ÿ b b ÿ 1Ck ÿ Nk ÿ q l1 dt l1 m1 l1
db dt
k 1; 2; . . . ; N
21
with R ak t 0
X
v x; t 0 ÿ b 0vr /k dX R : / /k dX X k
22
In Eq. (21), Z Mk
/k /k dX;
23
/k vr dX;
24
X
Z Nk
X
Z Hkl
T
X
r/k : r/l dX;
25
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
Z Qklm
X
/k /l r/m dX;
185
26
Z Rkl
X
/k /l rvr vr r/l dX;
27
/k vr rvr dX;
28
Z Ck
X
where X is the system domain. The set of ordinary dierential equations (21) is solved by using the Adams±Bashforth method for the nonlinear terms and the Crank±Nicolson method for the linear terms. The number of empirical eigenfunctions employed, N, is determined by comparing the result of the low-dimensional dynamic model with that of the numerical solution of the original partial dierential equation for given a t. Usually, the error of the low-dimensional model decreases as the number of empirical eigenfunctions employed increases up to the optimal number. But further increase of number of empirical eigenfunctions beyond the optimal number does not always improve the accuracy because the empirical eigenfunctions with very small eigenvalues are contaminated with round-o errors. The optimal number of eigenfunctions for the set (21) is found to be 8. With this number of empirical eigenfunctions, the relative error of the low-di-
Fig. 4. Various control functions a t considered to investigate the accuracy of the low-dimensional dynamic model.
186
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
mensional dynamic model with respect to the ®nite volume solution is about 1%. In the sequel, all the results of the Karhunen±Loeve Galerkin procedure are based on the low-dimensional model constructed using 8 empirical eigenfunctions. To demonstrate the robustness of the low-dimensional dynamic model with respect to the variation of the inlet velocity, the low-dimensional dynamic model has been solved using various shapes of control a t and the results are compared with those obtained by ®nite volume method. Fig. 4(a)±(d) plot various control functions a t and Fig. 5 shows the corresponding relative error
Fig. 5. Relative errors of the low-dimensional dynamic model with respect to the ®nite volume solution for various controls a t displayed in Fig. 4(a)±(d).
Fig. 6. The temporal variation of the velocity v, as computed by the low-dimensional dynamic model and by the ®nite volume method, respectively, at selected points indicated when the control a t is given by Fig. 4(d): (a) x-component; (b) z-component.
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
187
of the low-dimensional model with respect to the ®nite volume solution of the Navier±Stokes equation. These results demonstrate that the set of empirical eigenfunctions obtained in this section spans the admissible solution space of the Navier±Stokes equation for various a t considered. Fig. 6(a) and (b) shows the comparison of velocity components obtained by the Karhunen±Loeve Galerkin procedure with those by the ®nite volume solution of the Navier±Stokes equation at several locations indicated when a t varies according to Fig. 4(d). It is shown that both results are almost the same. 4. Solution of optimal control problems employing the Navier±Stokes equation In this and the following sections, we solve the optimal control problems posed in Section 2 employing the Navier±Stokes equation and the low-dimensional dynamic model, respectively, by means of the same conjugate gradient technique so that the accuracy and eciency of the Karhunen±Loeve Galerkin procedure be assessed. As the ®rst step, we transform the continuous control function a t into discrete variables as a t
M X m1
bm wm t;
29
where M is the number of shape functions used to discretize a t, which is 31 in the present work, and wm t is the mth shape function depicted in Fig. 7. Among 31 discrete variables, b1 equals the value of a t at t 0. Therefore, we set b1 a t 0 and the problem at hand is to select the vector b of the discrete control variables bT b2 ; b3 ; . . . ; bM
30
Fig. 7. De®nition of shape functions.
188
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
such that the objective function, Eq. (8), shall be minimized. The conjugate gradient method based on the Fletcher±Reeves algorithm [12] is employed in the iterative minimization of the objective function. The variation of the objective function with respect to the vector b is determined by the following gradient vector:
oJ ob
T
oJ oJ oJ ; ;...; ob2 ob3 obM
;
31
where oJ obm
Z
ov v x; T ÿ vT x dX obm X
Z
T 0
a t ÿ a0 wm t dt
32
and
ov obm
T
ou ow : ; obm obm
33
For the brevity of notation, we denote: s m s x; y; t; m
ou x; y; t ; obm
34
h m h x; y; t; m
ow x; y; t ; obm
35
p m p x; y; t; m
1 op : q obm
36
and
On the other hand, the variation of the velocity ®eld with respect to the vector b is
ov ob
T
ov ov ov : ; ;...; ob2 ob3 obM
37
The governing equations for the elements of the sensitivity matrix ov=ob are obtained by partially dierentiating Eqs. (1)±(6) with respect to the parameter vector b. For values of m 2; 3; . . . ; M we have the following set of equations:
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
o o ou o ou s m u s m s m w s m h m ot ox ox oz oz 2 2 o o o ÿ p m m s m 2 s m ; ox ox2 oz o o ow o ow h m u h m s m w h m h m ot ox ox oz oz 2 o o o2 ÿ p m m h m 2 h m ; oz ox2 oz o o s m h m 0 ox oz
189
38
39
40
with the initial conditions s m h m 0
41
and the relevant boundary conditions: 2 s m Wm t 1 ÿ h z ;
x 0 inlet :
x 1:0 outlet :
o s m 0; ox
at all other boundaries :
h m 0;
42
h m 0;
43
s m h m 0:
44
These sets of equations (38)±(44), can be solved using the ®nite volume method [10]. The conjugate direction or the search direction d is calculated and renewed by the Fletcher±Reeves method: d
i1
i1 oJ ÿ ud i ob
i P 1;
45
where PM u
m2
PM
m2
oJ obm
i1 2
oJ obm
i 2 :
46
190
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
On the other hand, the search direction d at the ®rst step is computed by d 0 ÿ
oJ : ob
47
This renewed d i1 is used for the search direction at the next iterative state. The parameter vector b i1 at the i 1th iterative state is obtained from b i through moving in the conjugate direction d i: b i1 b i ad i :
48
The optimal step length a is determined by minimizing J b ad with respect to a. Exploiting the formula, v b ad v b
M X ov adm ; obm m2
49
where dm is the mth component of the conjugate direction vector d, i.e., d T d2 ; d3 ; . . . ; dM ;
50
we ®nd the expression for the optimal step length. B aÿ ; A
51
where Z A
M X X
!2 s mdm
m2
M X
!2 h mdm
Z dX
m2
0
T
M X m2
!2 dm wm t
dt
52
and Z B
X
u x; t ÿ uT x
! s mdm
m2
Z
M X
T 0
M X m1
bm Wm t ÿ a0
!
v x; t ÿ vT x
M X m2
! h mdm
dX
m2
! dm wm
M X
dt:
53
The iterative process of the conjugate gradient method is summarized as follows: 1. Assume b. Then a t is determined by Eq. (29). Obtain u x; t, w x; t, s x; t; m, h x; t; m by solving Eqs. (1)±(6) and Eqs. (38)±(44), respectively.
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
191
2. Calculate d 0 ÿ oJ =ob by Eq. (32). 3. Determine the optimal step length by Eq. (51). 4. Calculate b i1 by Eq. (48). 5. Obtain v b i1 , s b i1 ; m and h b i1 ; m by solving Eqs. (1)±(6) and Eqs. (38)±(44) , respectively. 6. Calculate u using Eq. (46). 7. Calculate d i1 by Eq. (45). If jd i1 j < , then stop. 8. Set i i 1 and go to step 3. 5. Solution of optimal control problems employing the low-dimensional dynamic model If the conjugate gradient method of Section 4 is applied to the low-dimensional dynamic model, the minimization of the objective function can be done very eciently. The d.f. of Eq. (21) is only 8, whereas that of the original partial dierential equation, which is equivalent to the grid number in the ®nite volume approximation multiplied by 3 (i.e., u, w, p), is about 2 104 . The drastic reduction of the d.f. through the Karhunen±Loeve Galerkin procedure shall correspondingly reduce the computer time in the solution of optimal control problem. It is convenient to reformulate the equations of the low-dimensional dynamic model by regarding b t as an additional state variable, aN 1 : # " N N X N N X X X dak 1 al Hkl ÿ al am Qklm ÿ aN1 al Rkl ÿ aN1 aN 1 ÿ 1Ck ÿ Nk c ÿm Mk dt l1 l1 m1 l1 k 1; 2; . . . ; N ; daN1 c: dt
54
55
The control variable in this formulation is c t. The objective function, Eq. (8), is rewritten in terms of the spectral coecients, ai , as J b
N N ÿ ÿ ÿ 2 X 1X ai T ÿ aTi Mi aN 1 T ÿ aTN1 ai T ÿ aTi Ni 2 i1 i1 Z Z T ÿ 1 2 aN 1 T ÿ aTN1 vr vr dX ar aN 1 ÿ a0 2 dt; 2 2 0 X
56
where Mi and Ni are given by Eqs. (23) and (24) and aTi i 1; 2; . . . ; N 1 are determined by exploiting vT x
N X i1
aTi /i x aTN1 vr x:
57
192
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
At the inlet (x 1), /i x is zero and Eq. (57) is reduced to: "
#
uT 1; y wT 1; y 0
" aTN1
ar 1 ÿ h z2 0
# :
58
Since uT 1; y aT 1 ÿ h z2 ;
59
we ®nd aTN1
aT : ar
60
The remaining coecients aTi i 1; 2; . . . ; N are given by R ÿ aTi
X
vT ÿ aTN 1 vr /i dX R : / /i dX X i
61
As previously, the continuous control variable c t is discretized as c t
M X m1
bm
1 dwm : ar dt
62
This follows from Eq. (29) and the fact that b t a t=ar . Since b1 is ®xed by the initial value of control, a 0, the minimization is performed over the parameter vector bT b2 ; b3 ; . . . ; bM . The gradient of the objective function with respect to b is N N X oai T ÿ ÿ ÿ oJ oaN 1 T X ai T ÿ aTi Mi ai T ÿ aTi Ni aN 1 T ÿ aTN1 obm obm obm i1 i1 Z N X oai T oaN 1 T ÿ T Ni aN 1 T ÿ aN1 vr vr dX obm obm X i1 Z T oaN 1 dt: ar aN 1 t ÿ a0 ar obm 0
The sensitivity equations that determine oak =obm are:
63
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
193
"
N N X N N X N X X X d oak 1 oal oal oai al ÿm Hkl ÿ ai Qkli ÿ Qkli dt obm Mk obm obm obm l1 l1 i1 l1 i1 N N X oaN1 X oal oaN 1 oaN 1 Rkl ÿ 2aN 1 ÿ Ck ÿ al Rkl ÿ aN1 obm obm obm obm l1 l1 # dwm 1 k 1; 2; . . . ; N ; m 2; 3; . . . ; M; 64 ÿ Nk dt ar d dt
oaN1 obm
1 dwm ar dt
m 2; 3; . . . ; M:
65
The initial conditions for the sensitivity equations are: oak t 0 0; obm
66
oaN1 t 0 0 obm
k 1; 2; . . . ; N; m 2; 3; . . . ; M:
67
The conjugate direction at the i 1th iteration is determined by d
i1
i1 oJ ÿ ud i ; ob
68
where PM u
m2
PM
m2
oJ obm
i1 2
oJ obm
i 2
69
and d T d2 ; d3 ; . . . ; dM . The parameter vector is updated in the conjugate direction according to the following equation: b i1 b i ad i :
70
The optimal step length a is determined by partially dierentiating J b ad with respect to a and setting the resulting equation equal to zero: B aÿ ; A
71
194
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
where " #2 ! ! N M M N M X X X X X oai T oaN1 T oai T A dm Mi 2 dm d m Ni obm obm obm i1 m2 m2 i1 m2 !2 Z !2 Z T M M X X oaN1 T oa t N 1 dm vr vr dX a2r dm dt ob obm 0 X m m2 m2
72
and " ! # N M M N X X X X ÿ ÿ oa T oa T i N 1 B ai T ÿ aTi dm Mi dm ai T ÿ aTi Ni obm obm i1 m2 m2 i1 ! N M M ÿ ÿ X X oai T X oaN 1 T T T aN 1 T ÿ aN 1 dm Ni aN 1 T ÿ aN1 dm ob obm m i1 m2 m2 ! Z Z T M X oaN1 t vr vr dX ar ar aN1 ÿ a0 dm dt: obm 0 X m2
73
The following relation has been exploited to derive Eqs. (71)±(73): ai b ad ai a
M X oai dm : obm m2
74
The algorithm of the optimal control problem employing the low-dimensional dynamic model follows the following procedure: 1. Assume b b2 ; b3 ; . . . ; bM and calculate ai t i 1; 2; . . . ; N and oai =obm t i 1; 2; . . . ; N ; m 2; 3; . . . ; M using Eqs. (54) and (55) and Eqs. (64) and (65), respectively. 2. Calculate d 0 ÿoJ =ob0 by Eq. (63). 3. Determine the optimal step length a by Eqs. (71)±(73). 4. Update b i1 by Eq. (70). 5. Calculate ai t i 1; 2; . . . ; N and oai =obm t i 1; 2; . . . ; N ; m 2; 3; . . . ; M using Eqs. (54) and (55) and Eqs. (64) and (65), respectively. 6. Calculate u using Eq. (69). 7. Calculate d i1 by Eq. (68). If jd i1 j < , then stop. 8. Set i i 1 and go to step 3. 6. Results The optimal control problem of the incompressible viscous ¯ow through the grooved channel, posed in Section 2, is solved by the traditional method explained in Section 4 and by the Karhunen±Loeve Galerkin procedure of Section 5, respectively, and the accuracy and eciency of
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
195
these methods are compared. This is an extension of the previous work done for the Burgers equation [1] to the case of the multi-dimensional Navier±Stokes equations. As before, we call the method employing the original partial dierential equation the SM±CG, while the method employing the low-dimensional model is called the KLG±CG. The optimal control problem is to determine the control parameter a t, the maximum value of the parabolic inlet velocity, that yields the target velocity ®eld, depicted in Fig. 1(b), at the ®nal time with the minimum value of the objective function J (cf. Eq. (8)). For 1 10ÿ3 in Eq. (8), the process of the minimization of J with respect to the iteration number is depicted in Fig. 8 for both the FDM±CG and the KLG±CG, when the initial approximation of the control is the trial control, shown in Fig. 1(a), employed to obtain the target velocity ®eld, Fig. 1(b). Fig. 8 shows that the convergence rate of the KLG±CG is faster than that of the FDM±CG at the initial stage, but the minimum value of the objective function is attained only after 200 iterations with both methods. The converged pro®les of optimal control obtained by the FDM±CG and the KLG±CG are plotted in Fig. 9. Also shown in the same ®gure is the trial control a t, given by Eq. (9), which has been used to obtain the target velocity ®eld, for comparison. Fig. 9 reveals that the optimal controls aFDM and aKLG obtained by the FDM±CG and the KLG±CG, respectively, are almost the same, the cost of both aFDM and aKLG being much less than that of the trial control atrial . Contrary to the trial control, the optimal control is almost zero initially but increases to a high value at the ®nal stage to ®t the target velocity ®eld. To assure that the KLG±CG and the FDM±CG yield correct velocity ®eld at the ®nal time, the velocity ®elds obtained by the KLG±CG and the FDM± CG at the ®nal time are compared with the target velocity ®eld at a cross-section of the domain. Fig. 10 shows that the states at the ®nal time, obtained by the FDM±CG (Fig. 10(c) and (d)) and the KLG±CG (Fig. 10(e) and (f)) are indistinguishable from the target velocity ®eld (Fig. 10(a) and (b)). Fig. 11 compares the transient ¯ow ®elds induced by the trial control, Eq. (9), and the optimal control. It is observed that the ¯ow ®eld caused by the trial control experiences larger
Fig. 8. Convergence rate of the conjugate gradient method (the FDM±CG and the KLG±CG).
196
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
Fig. 9. The pro®les of the optimal control obtained either by the FDM±CG or by the KLG±CG for the target velocity ®eld of Fig. 1(b). The trial control, (9), is also displayed for comparison.
excursion than that caused by the optimal control before settling down to the target velocity ®eld, implying that the former incurs larger control cost than the latter. As expected the consumption of CPU time for the case of the KLG±CG is much less than that for the case of FDM±CG since the d.f. for the former is only about 1/2000 of the latter. In fact, when the ultrasparc workstation is used, the CPU time requirement per one iteration of the conjugate gradient method for the FDM±CG and KLG±CG are about 22 h and 1.5 min, respectively. Since the dierence in the d.f. between the low-dimensional dynamic models and the original partial dierential equations shall become much larger as the dimensionality of the problem increases, the reduction of CPU time with the use of the KLG±CG instead of the FDM± CG shall be much more signi®cant in the case of three-dimensional problems. The CPU time requirement to construct the low-dimensional dynamic model is not signi®cant compared with that needed for the FDM±CG. To construct the low-dimensional dynamic model, it is necessary to solve the Navier±Stokes equation with two dierent control function a t to generate snapshots. This step requires 1.2 h. It takes 2 h to obtain the reference velocity ®eld vr . The remaining steps in the Karhunen±Loeve Galerkin procedure consume about 2 h. All in all the construction of the low-dimensional dynamic model requires about 5.2 h, which is much less than the CPU time requirement of just one iteration of the FDM±CG. Moreover, once the low-dimensional dynamic model is constructed, it can be used repeatedly for similar problems. Because the optimization problem under consideration is nonconvex, the proof for existence and uniqueness of the solution is not available in general. But the existence of the optimal control can be assured since the target velocity ®eld has been obtained by solving the Navier±Stokes equation with a given control a t. The uniqueness of the optimal control is con®rmed experimentally by solving the optimization problem using various initial approximation of a t as shown in Fig. 12. Fig. 12(a) shows the convergence rate of the KLG±CG when using dierent initial
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
197
Fig. 10. The comparison of the ®nal state, v x; T , obtained either by the FDM±CG or by the KLG±CG, with the target velocity ®eld vT x. These three velocity pro®les are virtually indistinguishable.
198
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
Fig. 11. The comparison of the transient ¯ow ®elds when a t is given by Eq. (9) and when a t is the optimal control.
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
199
Fig. 12. Convergence of the KLG±CG for various initial approximations of a t: (a) the convergence rate; (b) the converged pro®les of optimal control a t for various initial approximations.
approximations of a t, and the corresponding converged pro®les of optimal control are plotted in Fig. 12(b). This ®gure reveals that the ®nal converged pro®les of the optimal control a t are the same regardless of the initial approximation within numerical error.
7. Conclusion The Karhunen±Loeve Galerkin procedure is applied to the optimal control problem of viscous ¯uid ¯ow through a grooved channel. The Karhunen±Loeve Galerkin procedure, which is a Galerkin method employing the empirical eigenfunctions of the Karhunen±Loeve decomposition as trial functions, limits a priori the function space to the smallest linear subspace that is sucient to describe the observed solutions, and consequently reduces the Navier±Stokes equation to a minimal set of ordinary dierential equations. The low-dimensional dynamic model derived in this
200
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
way simulates the system accurately for a given space of control functions a t, and can be employed in the solution of optimal control problems eciently. This technique, called the KLG± CG, is compared with the traditional technique of employing the Navier±Stokes equation (called the FDM±CG), and is found to yield accurate results at the drastically reduced computational cost. This is an extension of our previous work where the same technique has been used to solve the optimal control problems of the one-dimensional Burgers equation [1]. Since the reduction of modes through the Karhunen±Loeve Galerkin procedure is much larger in the cases of multidimensional Navier±Stokes equations than in the case of one-dimensional Burgers equation, the KLG±CG is found to be more powerful as regards the computational eciency when applied to the control problems of the Navier±Stokes equation than those of the Burgers equation. Although the importance of control of ¯uid ¯ow has been recognized for a long time, the development of rigorous as well as practical methodology of ¯ow control has been retarded due to the inherent mathematical and computational complexities of the Navier±Stokes equations. In this regard, the Karhunen±Loeve Galerkin procedure is suggested as one of the ecient methods of solving problems of ¯ow control, which is practical enough to be implemented in industrial processes. References [1] H.M. Park, Y.D. Jang, Control of Burgers equation by means of mode reduction, Int. J. Engrg. Sci. 38 (2000) 785± 805. [2] M.D. Gunzburger, L. Hou, T.P. Svobodny, Boundary velocity control of incompressible ¯ow with an application to viscous drag reduction, SIAM J. Control Optim. 30 (1991) 167±181. [3] P. Moin, R.D. Moser, Characteristic-eddy decomposition of turbulence in a channel, J. Fluid Mech. 200 (1989) 471±506. [4] H.M. Park, L. Sirovich, Turbulent thermal convection in a ®nite domain: numerical results, Phys. Fluids A 2 (9) (1990) 1659±1668. [5] N. Aubry, P. Holmes, J.L. Lumley, E. Stone, The dynamics of coherent structures in the wall region of a turbulent boundary layer, J. Fluid Mech. 192 (1988) 115±173. [6] A.E. Deane, I.G. Kevrekidis, G.E. Karniadakis, S.A. Orszag, Low-dimensional models for complex geometry ¯ows: application to grooved channels and circular cylinders, Phys. Fluids 3 (1991) 23±27. [7] P. Moin, T. Bewley, Feedback control of turbulence, Appl. Mech. Rev. 47 (6), part 2, (1994) s3±s13. [8] H.M. Park, D.H. Cho, The use of the Karhunen±Loeve decomposition for the modeling of distributed parameter systems, Chem. Eng. Sci. 51 (1) (1996) 81±98. [9] H.M. Park, M.W. Lee, An ecient method of solving the Navier±Stokes equation for the ¯ow control, Int. J. Numer. Meth. Engrg. 41 (1998) 1131±1151. [10] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, McGraw-Hill, New York, 1980. [11] R. Courant, D. Hilbert, Methods of Mathematical Physics, vol. 1, Interscience, New York, 1953. [12] R. Fletcher, R.M. Reeves, Function minimization by conjugate gradients, The Computer J. 7, 149±154.
www.elsevier.com/locate/ijengsci
Reduction of modes for the control of viscous ¯ows H.M. Park *, O.Y. Kim Department of Chemical Engineering, Sogang University, Shinsoo-Dong, Mapo-Gu, Seoul, South Korea Received 28 July 1999; received in revised form 17 September 1999; accepted 17 January 2000
Abstract This paper applies the Karhunen±Loeve Galerkin procedure to the solution of optimal control problems of viscous ¯uid ¯ow through a grooved channel. Through the Karhunen±Loeve Galerkin procedure, one can reduce the Navier±Stokes equation to a small number of ordinary dierential equations by limiting the solution space to the smallest linear subspace that is sucient to describe the observed phenomena. This technique is well suited for the problem of control or optimization, where one has to solve the governing equation repeatedly but one can also estimate the approximate solution space based on the range of control variables. In our previous work [H.M. Park, Y.D. Jang, Int. J. Engrg. Sci. 38 (2000) 785±805], we showed that the Karhunen±Loeve Galerkin procedure could solve an optimal control problem of one-dimensional Burgers equation very eciently. The present paper extends the applicability of the Karhunen±Loeve Galerkin procedure to the solution of optimal control problems of multi-dimensional Navier±Stokes equations. Ó 2000 Elsevier Science Ltd. All rights reserved. Keywords: Control of viscous ¯ows; Karhunen±Loeve Galerkin procedure
1. Introduction Flow control problems have been addressed by many scientists, mathematicians and engineers due to its frequent occurrence in various systems in science and engineering. The control of ¯uid ¯ow is crucial to achieve desired design objectives or to optimize the performance of the systems that exploit ¯uid motion. The development of a technique of ¯ow control, which is ecient as well as rigorous, is one of the central issues to be resolved before achieving optimal performance and safe operation of various systems with ¯uid ¯ow. But the quest for successful control of viscous ¯uid remains one of the most arduous challenges in the ¯uid and control sciences today. The *
Corresponding author. E-mail address: [email protected] (H.M. Park).
0020-7225/01/$ - see front matter Ó 2000 Elsevier Science Ltd. All rights reserved. PII: S 0 0 2 0 - 7 2 2 5 ( 0 0 ) 0 0 0 2 4 - 0
178
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
Navier±Stokes equation, which is the governing equation for the viscous ¯uid ¯ow, has almost in®nite degrees of freedom (d.f.) with strong nonlinearity, thus causing dicult computational and theoretical problems. With the optimal control problems in ¯uid ¯ows, candidate optimal states and controls may belong to in®nite-dimensional function spaces and one may have to minimize a nonlinear functional of the state and control variables subject to nonlinear constraints that take the form of a system of partial dierential equations [2]. As a consequence the optimization problems are nonconvex, and this creates some technical diculties in studying the existence and uniqueness of an optimal control. Perhaps one of the most important obstacles to the application of rigorous control theories to viscous ¯uid ¯ows is the tremendously large requirement of computer time and memory. The actual implementation of control schemes requires repeated computation of the governing and related equations, and one cannot expect on-line control without reducing the number of nonlinear equations in the system. Therefore, one of the most crucial factors that determines the success of the numerical algorithm for the optimal control problems of ¯uid ¯ow is the employment of very ecient numerical method to solve the system equations. The traditional methods of the discretization of partial dierential equations usually involve functions that have very little to do with the dierential equations. For example, piecewise polynomials are used in ®nite element methods, grid functions are used in ®nite dierence methods, and trigonometric functions or Chebyshev polynomials are used in spectral methods. As a result they yield large systems of equations whose inversion consumes most of the computation time. Recently, several research groups [3,4] have used the Karhunen±Loeve decomposition as a procedure for computation of coherent structures (empirical eigenfunctions) in turbulence from experimental or numerical data of ¯ow ®elds. Usually these empirical eigenfunctions have been used to interpret the statistical characteristics of the turbulent ¯ow. Furthermore, employing these empirical eigenfunctions and introducing additional drastic approximations to the Navier±Stokes equations they could derive low-dimensional dynamic models for turbulent ¯ows [5] or transitional ¯ows [6]. These low-dimensional models [5,6] do not simulate the ¯ow ®eld exactly when the control variables like the Reynolds number varies, but are found to approximately reproduce some interesting characteristics of turbulence or transitional ¯ows such as intermittency and bifurcation sequences. Also there have been some attempts to apply the dynamic system methodology to these approximate low-dimensional models for the purpose of ¯ow control [7]. Because the original Karhunen±Loeve decomposition technique [4±6] is applicable only to stationary stochastic ®elds such as turbulence or oscillatory transitional ¯ows, it could not be employed in the rigorous control schemes of ¯uid ¯ows where the control variables vary. But recent works [1,8,9] have extended the applicability of the Karhunen±Loeve decomposition to the analysis of nonstationary, nonhomogeneous deterministic as well as stochastic ®elds and allowed the derivation of rigorous low-dimensional dynamic models that simulate the given systems almost exactly even when the control variables vary arbitrarily within certain bounds. This extension of the original Karhunen±Loeve decomposition is called the Karhunen±Loeve Galerkin procedure [1,8,9] and can be adopted in rigorous implementation of modern control schemes for viscous ¯uid ¯ows. In our previous work [1], we showed that the Karhunen±Loeve Galerkin procedure could solve an optimal control problem of one-dimensional Burgers equation very eciently. The present paper applies the Karhunen±Loeve Galerkin procedure to the solution of optimal control
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
179
problems of multi-dimensional Navier±Stokes equations. Since the reduction of modes in the multi-dimensional case is much larger than that in the one-dimensional case, the present technique is found to be more powerful when applied to the control problems of the Navier±Stokes equation than those of the Burgers equation. The details of the Karhunen±Loeve Galerkin procedure is well documented in our previous works [1,8,9] and reference therein. 2. System and governing equations We consider a viscous incompressible ¯ow through a grooved channel as shown in Fig. 1(b). The ¯ow is induced by the time-varying inlet velocity U t. The unsteady ¯ow ®eld of the system is obtained by solving the following incompressible Navier±Stokes equation with relevant initial and boundary conditions: ov v rv ÿrP lr2 v; 1 q ot r v 0;
2
where v u; w is the velocity vector, P the pressure, q the density ( 1) and l is the viscosity ( 1). The initial and boundary conditions are: t 0;
v v0 x;
x 0 inlet :
3
u U t a t 1 ÿ h z2 ;
w 0;
Fig. 1. (a) The trial control; (b) the target velocity ®eld.
4
180
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
x 1:0 outlet :
ou 0; ox
at all other boundaries :
w 0;
5
u w 0;
6
where h z z ÿ 0:875=0:125 for 0:75 6 z 6 1:0 and a t is the maximum of the inlet velocity of parabolic pro®le. This set of Eqs. (1)±(6) is solved accurately by using the ®nite volume method [10]. The problem at hand is that of controlling the system above to produce a state variable v x; T at the ®nal time T that is as close as possible to a target ¯ow ®eld vT x. We aim at achieving this goal by adjusting the maximum of the inlet velocity a t. Moreover, we try to achieve the goal at the minimum cost of the control a t. The problem may be speci®ed as the following optimization problem: min J a;
7
where 1 J a 2
Z
v x; T ÿ vT x dX 2 X 2
Z 0
T
2
a t ÿ a0 dt;
8
a0 the initial value of a and is a positive weighting factor which represents the expense of applying the control. The magnitude of is small if the control is cheap and large if it is expensive. To ensure that the target vT x is a reachable state, we determine vT x by integrating Eqs. (1)±(6) with a given function a t and setting vT x v x; t T . The trial control a t employed to obtain vT x is given by 80t 1 0 6 t 6 0:05; a t 9 ÿ80 t ÿ 0:1 1 0:05 6 t 6 0:1: The trial control (9) and the corresponding target ¯ow ®eld vT x are depicted in Fig. 1(a) and (b), respectively. 3. Low-dimensional dynamic model In this section, we employ the Karhunen±Loeve Galerkin procedure to derive a low-dimensional dynamic model which shall be used in the solution of the optimal control problems. The low-dimensional dynamic model obtained through the Karhunen±Loeve Galerkin procedure must predict the system behavior accurately for various trajectories of a t, which is the maximum of the inlet velocity, that include not only the optimal trajectory but also other trajectories appearing during the iterative minimization of the objective function equation (8) by means of the conjugate gradient method. The prerequisite to this requirement is that the empirical eigenfunctions employed in the Karhunen±Loeve Galerkin procedure must span the admissible solution space of the Navier±Stokes equation for these various trajectories of a t. According to the
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
181
Schmidt±Hilbert theory [1,11], the empirical eigenfunction /i x can be expressed linearly in terms of N snapshots: /i x
N X i1
ai j vj x;
10
where vj x is the jth snapshot. Therefore, once we obtain an ensemble of snapshots fvn xg which encompasses the admissible solution space of the Navier±Stokes equation for the above trajectories of a t, the resulting empirical eigenfunctions obtained through the Karhunen±Loeve decomposition shall span the same solution space, and consequently the low-dimensional model shall simulate the system accurately for those various trajectories of a t. In the present work, it is assumed that the admissible range of a t is between 1 and 5. Then for any arbitrary trajectory a t, we ®nd that the variation of a t during the time increment Dt is Da a t Dt ÿ a t, which can be positive or negative but its absolute value must be less than 4. This simple observation suggests the strategy to obtain snapshots. The snapshots of the ¯ow ®eld are obtained by solving Eqs. (1)±(6) when the maximum of the inlet velocity a t changes abruptly from 1 to 5 at t 0. The initial condition is taken as the steady ¯ow ®eld when the inlet velocity is 1, and the transient ¯ow ®elds are recorded at an appropriate time interval until a new steady state is reached. These recordings are served as snapshots to be used in the Karhunen±Loeve decomposition. In this way, we obtain the ®rst set of 500 snapshots. The second set of 500 snapshots are obtained by solving the governing equations (1)±(6) with the abrupt change of a t from 5 to 1 at t 0 when the initial condition is the steady ¯ow ®eld of a 5:0. In these computations, the ®nite volume method [10] with 80 80 grids is employed, which is found to resolve the ¯ow ®elds accurately. As explained by Park and Cho [8], it is necessary to make these snapshots satisfy homogeneous boundary conditions. For this purpose, we obtain the reference steady ¯ow ®eld vr ur ; wr by solving the following set of equations: qvr rvr ÿrPr lr2 vr ;
11
r vr 0
12
with the following boundary conditions: 2 inlet : ur ar 1 ÿ h z ; wr 0;
outlet :
our 0; ox
wr 0;
at all other boundaries :
ur wr 0;
13
14 15
where ar 3. We make the ®rst set of 500 snapshots satisfy homogeneous boundary condition by the following arrangement:
182
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
5 v x; t ÿ vr x ; 3 n
n 1; 2; . . . ; 500:
Similarly, the second set is transformed to satisfy homogeneous boundary condition as 1 n v x; t ÿ vr x ; n 1; 2; . . . ; 500: 3
16
17
To these 1000 snapshots satisfying homogeneous boundary conditions, we apply the Karhunen±Loeve decomposition technique to get empirical eigenfunctions in the order of their importance in characterizing the system. Fig. 2(a)±(d) shows the ®rst, second, third and fourth
Fig. 2. Some dominant empirical eigenfunctions with large eigenvalues: (a) ®rst eigenfunction k1 0:9413; (b) second eigenfunction k2 4:943 10ÿ2 ; (c) third eigenfunction k3 7:042 10ÿ3 ; (d) fourth eigenfunction k4 1:565 10ÿ3 .
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
183
empirical eigenfunctions with the corresponding normalized eigenvalues k1 0:9413, k2 4:943 10ÿ2 , k3 7:042 10ÿ3 , k4 1:565 10ÿ3 , respectively. Also shown in Fig. 3(a)± (d) are typical empirical eigenfunctions with smaller eigenvalues, i.e., the 5th±8th eigenfunctions with the corresponding normalized eigenvalues, k5 3:595 10ÿ4 , k6 1:935 10ÿ4 , k7 6:809 10ÿ5 , k8 1:835 10ÿ5 : In these ®gures, solid lines denote positive-valued streamlines and dotted lines denote negative-valued streamlines. As previously [1,8], the empirical eigenfunctions with large eigenvalues represent the large-scale structures of the velocity ®eld v. On the contrary, Fig. 3(a)±(d) reveals that empirical eigenfunctions with small eigenvalues represent the small-scale structures caused by system boundaries and nonlinearity of the governing equation.
Fig. 3. Empirical eigenfunctions with small eigenvalues: (a) ®fth eigenfunction k5 3:595 10ÿ4 ; (b) sixth eigenfunction k6 1:935 10ÿ4 ; (c) seventh eigenfunction k7 6:809 10ÿ5 ; (d) eighth eigenfunction k8 1:835 10ÿ5 .
184
as
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
The velocity ®eld v x; t is now represented as a linear combination of empirical eigenfunctions v
n X i1
ai t/i x b tvr ;
18
where /i is the ith empirical eigenfunction, ai t the corresponding spectral coecient, b t a t=ar where a t is the maximum of the inlet velocity and ar 3, and N is the total number of empirical eigenfunctions employed in the Karhunen±Loeve Galerkin procedure. Applying the Galerkin principle which enforces the residual R to be orthogonal to each of the basis functions, Z X
R /i dX 0
i 1; 2; . . . ; N;
19
where R
ov P l v rv r ÿ r2 v; ot q q
20
Eqs. (1)±(6) are reduced to the following set of ordinary dierential equations: Mk
N N X N N X X dak lX al Hkl ÿ al am Qklm ÿ b al Rkl ÿ b b ÿ 1Ck ÿ Nk ÿ q l1 dt l1 m1 l1
db dt
k 1; 2; . . . ; N
21
with R ak t 0
X
v x; t 0 ÿ b 0vr /k dX R : / /k dX X k
22
In Eq. (21), Z Mk
/k /k dX;
23
/k vr dX;
24
X
Z Nk
X
Z Hkl
T
X
r/k : r/l dX;
25
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
Z Qklm
X
/k /l r/m dX;
185
26
Z Rkl
X
/k /l rvr vr r/l dX;
27
/k vr rvr dX;
28
Z Ck
X
where X is the system domain. The set of ordinary dierential equations (21) is solved by using the Adams±Bashforth method for the nonlinear terms and the Crank±Nicolson method for the linear terms. The number of empirical eigenfunctions employed, N, is determined by comparing the result of the low-dimensional dynamic model with that of the numerical solution of the original partial dierential equation for given a t. Usually, the error of the low-dimensional model decreases as the number of empirical eigenfunctions employed increases up to the optimal number. But further increase of number of empirical eigenfunctions beyond the optimal number does not always improve the accuracy because the empirical eigenfunctions with very small eigenvalues are contaminated with round-o errors. The optimal number of eigenfunctions for the set (21) is found to be 8. With this number of empirical eigenfunctions, the relative error of the low-di-
Fig. 4. Various control functions a t considered to investigate the accuracy of the low-dimensional dynamic model.
186
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
mensional dynamic model with respect to the ®nite volume solution is about 1%. In the sequel, all the results of the Karhunen±Loeve Galerkin procedure are based on the low-dimensional model constructed using 8 empirical eigenfunctions. To demonstrate the robustness of the low-dimensional dynamic model with respect to the variation of the inlet velocity, the low-dimensional dynamic model has been solved using various shapes of control a t and the results are compared with those obtained by ®nite volume method. Fig. 4(a)±(d) plot various control functions a t and Fig. 5 shows the corresponding relative error
Fig. 5. Relative errors of the low-dimensional dynamic model with respect to the ®nite volume solution for various controls a t displayed in Fig. 4(a)±(d).
Fig. 6. The temporal variation of the velocity v, as computed by the low-dimensional dynamic model and by the ®nite volume method, respectively, at selected points indicated when the control a t is given by Fig. 4(d): (a) x-component; (b) z-component.
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
187
of the low-dimensional model with respect to the ®nite volume solution of the Navier±Stokes equation. These results demonstrate that the set of empirical eigenfunctions obtained in this section spans the admissible solution space of the Navier±Stokes equation for various a t considered. Fig. 6(a) and (b) shows the comparison of velocity components obtained by the Karhunen±Loeve Galerkin procedure with those by the ®nite volume solution of the Navier±Stokes equation at several locations indicated when a t varies according to Fig. 4(d). It is shown that both results are almost the same. 4. Solution of optimal control problems employing the Navier±Stokes equation In this and the following sections, we solve the optimal control problems posed in Section 2 employing the Navier±Stokes equation and the low-dimensional dynamic model, respectively, by means of the same conjugate gradient technique so that the accuracy and eciency of the Karhunen±Loeve Galerkin procedure be assessed. As the ®rst step, we transform the continuous control function a t into discrete variables as a t
M X m1
bm wm t;
29
where M is the number of shape functions used to discretize a t, which is 31 in the present work, and wm t is the mth shape function depicted in Fig. 7. Among 31 discrete variables, b1 equals the value of a t at t 0. Therefore, we set b1 a t 0 and the problem at hand is to select the vector b of the discrete control variables bT b2 ; b3 ; . . . ; bM
30
Fig. 7. De®nition of shape functions.
188
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
such that the objective function, Eq. (8), shall be minimized. The conjugate gradient method based on the Fletcher±Reeves algorithm [12] is employed in the iterative minimization of the objective function. The variation of the objective function with respect to the vector b is determined by the following gradient vector:
oJ ob
T
oJ oJ oJ ; ;...; ob2 ob3 obM
;
31
where oJ obm
Z
ov v x; T ÿ vT x dX obm X
Z
T 0
a t ÿ a0 wm t dt
32
and
ov obm
T
ou ow : ; obm obm
33
For the brevity of notation, we denote: s m s x; y; t; m
ou x; y; t ; obm
34
h m h x; y; t; m
ow x; y; t ; obm
35
p m p x; y; t; m
1 op : q obm
36
and
On the other hand, the variation of the velocity ®eld with respect to the vector b is
ov ob
T
ov ov ov : ; ;...; ob2 ob3 obM
37
The governing equations for the elements of the sensitivity matrix ov=ob are obtained by partially dierentiating Eqs. (1)±(6) with respect to the parameter vector b. For values of m 2; 3; . . . ; M we have the following set of equations:
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
o o ou o ou s m u s m s m w s m h m ot ox ox oz oz 2 2 o o o ÿ p m m s m 2 s m ; ox ox2 oz o o ow o ow h m u h m s m w h m h m ot ox ox oz oz 2 o o o2 ÿ p m m h m 2 h m ; oz ox2 oz o o s m h m 0 ox oz
189
38
39
40
with the initial conditions s m h m 0
41
and the relevant boundary conditions: 2 s m Wm t 1 ÿ h z ;
x 0 inlet :
x 1:0 outlet :
o s m 0; ox
at all other boundaries :
h m 0;
42
h m 0;
43
s m h m 0:
44
These sets of equations (38)±(44), can be solved using the ®nite volume method [10]. The conjugate direction or the search direction d is calculated and renewed by the Fletcher±Reeves method: d
i1
i1 oJ ÿ ud i ob
i P 1;
45
where PM u
m2
PM
m2
oJ obm
i1 2
oJ obm
i 2 :
46
190
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
On the other hand, the search direction d at the ®rst step is computed by d 0 ÿ
oJ : ob
47
This renewed d i1 is used for the search direction at the next iterative state. The parameter vector b i1 at the i 1th iterative state is obtained from b i through moving in the conjugate direction d i: b i1 b i ad i :
48
The optimal step length a is determined by minimizing J b ad with respect to a. Exploiting the formula, v b ad v b
M X ov adm ; obm m2
49
where dm is the mth component of the conjugate direction vector d, i.e., d T d2 ; d3 ; . . . ; dM ;
50
we ®nd the expression for the optimal step length. B aÿ ; A
51
where Z A
M X X
!2 s mdm
m2
M X
!2 h mdm
Z dX
m2
0
T
M X m2
!2 dm wm t
dt
52
and Z B
X
u x; t ÿ uT x
! s mdm
m2
Z
M X
T 0
M X m1
bm Wm t ÿ a0
!
v x; t ÿ vT x
M X m2
! h mdm
dX
m2
! dm wm
M X
dt:
53
The iterative process of the conjugate gradient method is summarized as follows: 1. Assume b. Then a t is determined by Eq. (29). Obtain u x; t, w x; t, s x; t; m, h x; t; m by solving Eqs. (1)±(6) and Eqs. (38)±(44), respectively.
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
191
2. Calculate d 0 ÿ oJ =ob by Eq. (32). 3. Determine the optimal step length by Eq. (51). 4. Calculate b i1 by Eq. (48). 5. Obtain v b i1 , s b i1 ; m and h b i1 ; m by solving Eqs. (1)±(6) and Eqs. (38)±(44) , respectively. 6. Calculate u using Eq. (46). 7. Calculate d i1 by Eq. (45). If jd i1 j < , then stop. 8. Set i i 1 and go to step 3. 5. Solution of optimal control problems employing the low-dimensional dynamic model If the conjugate gradient method of Section 4 is applied to the low-dimensional dynamic model, the minimization of the objective function can be done very eciently. The d.f. of Eq. (21) is only 8, whereas that of the original partial dierential equation, which is equivalent to the grid number in the ®nite volume approximation multiplied by 3 (i.e., u, w, p), is about 2 104 . The drastic reduction of the d.f. through the Karhunen±Loeve Galerkin procedure shall correspondingly reduce the computer time in the solution of optimal control problem. It is convenient to reformulate the equations of the low-dimensional dynamic model by regarding b t as an additional state variable, aN 1 : # " N N X N N X X X dak 1 al Hkl ÿ al am Qklm ÿ aN1 al Rkl ÿ aN1 aN 1 ÿ 1Ck ÿ Nk c ÿm Mk dt l1 l1 m1 l1 k 1; 2; . . . ; N ; daN1 c: dt
54
55
The control variable in this formulation is c t. The objective function, Eq. (8), is rewritten in terms of the spectral coecients, ai , as J b
N N ÿ ÿ ÿ 2 X 1X ai T ÿ aTi Mi aN 1 T ÿ aTN1 ai T ÿ aTi Ni 2 i1 i1 Z Z T ÿ 1 2 aN 1 T ÿ aTN1 vr vr dX ar aN 1 ÿ a0 2 dt; 2 2 0 X
56
where Mi and Ni are given by Eqs. (23) and (24) and aTi i 1; 2; . . . ; N 1 are determined by exploiting vT x
N X i1
aTi /i x aTN1 vr x:
57
192
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
At the inlet (x 1), /i x is zero and Eq. (57) is reduced to: "
#
uT 1; y wT 1; y 0
" aTN1
ar 1 ÿ h z2 0
# :
58
Since uT 1; y aT 1 ÿ h z2 ;
59
we ®nd aTN1
aT : ar
60
The remaining coecients aTi i 1; 2; . . . ; N are given by R ÿ aTi
X
vT ÿ aTN 1 vr /i dX R : / /i dX X i
61
As previously, the continuous control variable c t is discretized as c t
M X m1
bm
1 dwm : ar dt
62
This follows from Eq. (29) and the fact that b t a t=ar . Since b1 is ®xed by the initial value of control, a 0, the minimization is performed over the parameter vector bT b2 ; b3 ; . . . ; bM . The gradient of the objective function with respect to b is N N X oai T ÿ ÿ ÿ oJ oaN 1 T X ai T ÿ aTi Mi ai T ÿ aTi Ni aN 1 T ÿ aTN1 obm obm obm i1 i1 Z N X oai T oaN 1 T ÿ T Ni aN 1 T ÿ aN1 vr vr dX obm obm X i1 Z T oaN 1 dt: ar aN 1 t ÿ a0 ar obm 0
The sensitivity equations that determine oak =obm are:
63
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
193
"
N N X N N X N X X X d oak 1 oal oal oai al ÿm Hkl ÿ ai Qkli ÿ Qkli dt obm Mk obm obm obm l1 l1 i1 l1 i1 N N X oaN1 X oal oaN 1 oaN 1 Rkl ÿ 2aN 1 ÿ Ck ÿ al Rkl ÿ aN1 obm obm obm obm l1 l1 # dwm 1 k 1; 2; . . . ; N ; m 2; 3; . . . ; M; 64 ÿ Nk dt ar d dt
oaN1 obm
1 dwm ar dt
m 2; 3; . . . ; M:
65
The initial conditions for the sensitivity equations are: oak t 0 0; obm
66
oaN1 t 0 0 obm
k 1; 2; . . . ; N; m 2; 3; . . . ; M:
67
The conjugate direction at the i 1th iteration is determined by d
i1
i1 oJ ÿ ud i ; ob
68
where PM u
m2
PM
m2
oJ obm
i1 2
oJ obm
i 2
69
and d T d2 ; d3 ; . . . ; dM . The parameter vector is updated in the conjugate direction according to the following equation: b i1 b i ad i :
70
The optimal step length a is determined by partially dierentiating J b ad with respect to a and setting the resulting equation equal to zero: B aÿ ; A
71
194
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
where " #2 ! ! N M M N M X X X X X oai T oaN1 T oai T A dm Mi 2 dm d m Ni obm obm obm i1 m2 m2 i1 m2 !2 Z !2 Z T M M X X oaN1 T oa t N 1 dm vr vr dX a2r dm dt ob obm 0 X m m2 m2
72
and " ! # N M M N X X X X ÿ ÿ oa T oa T i N 1 B ai T ÿ aTi dm Mi dm ai T ÿ aTi Ni obm obm i1 m2 m2 i1 ! N M M ÿ ÿ X X oai T X oaN 1 T T T aN 1 T ÿ aN 1 dm Ni aN 1 T ÿ aN1 dm ob obm m i1 m2 m2 ! Z Z T M X oaN1 t vr vr dX ar ar aN1 ÿ a0 dm dt: obm 0 X m2
73
The following relation has been exploited to derive Eqs. (71)±(73): ai b ad ai a
M X oai dm : obm m2
74
The algorithm of the optimal control problem employing the low-dimensional dynamic model follows the following procedure: 1. Assume b b2 ; b3 ; . . . ; bM and calculate ai t i 1; 2; . . . ; N and oai =obm t i 1; 2; . . . ; N ; m 2; 3; . . . ; M using Eqs. (54) and (55) and Eqs. (64) and (65), respectively. 2. Calculate d 0 ÿoJ =ob0 by Eq. (63). 3. Determine the optimal step length a by Eqs. (71)±(73). 4. Update b i1 by Eq. (70). 5. Calculate ai t i 1; 2; . . . ; N and oai =obm t i 1; 2; . . . ; N ; m 2; 3; . . . ; M using Eqs. (54) and (55) and Eqs. (64) and (65), respectively. 6. Calculate u using Eq. (69). 7. Calculate d i1 by Eq. (68). If jd i1 j < , then stop. 8. Set i i 1 and go to step 3. 6. Results The optimal control problem of the incompressible viscous ¯ow through the grooved channel, posed in Section 2, is solved by the traditional method explained in Section 4 and by the Karhunen±Loeve Galerkin procedure of Section 5, respectively, and the accuracy and eciency of
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
195
these methods are compared. This is an extension of the previous work done for the Burgers equation [1] to the case of the multi-dimensional Navier±Stokes equations. As before, we call the method employing the original partial dierential equation the SM±CG, while the method employing the low-dimensional model is called the KLG±CG. The optimal control problem is to determine the control parameter a t, the maximum value of the parabolic inlet velocity, that yields the target velocity ®eld, depicted in Fig. 1(b), at the ®nal time with the minimum value of the objective function J (cf. Eq. (8)). For 1 10ÿ3 in Eq. (8), the process of the minimization of J with respect to the iteration number is depicted in Fig. 8 for both the FDM±CG and the KLG±CG, when the initial approximation of the control is the trial control, shown in Fig. 1(a), employed to obtain the target velocity ®eld, Fig. 1(b). Fig. 8 shows that the convergence rate of the KLG±CG is faster than that of the FDM±CG at the initial stage, but the minimum value of the objective function is attained only after 200 iterations with both methods. The converged pro®les of optimal control obtained by the FDM±CG and the KLG±CG are plotted in Fig. 9. Also shown in the same ®gure is the trial control a t, given by Eq. (9), which has been used to obtain the target velocity ®eld, for comparison. Fig. 9 reveals that the optimal controls aFDM and aKLG obtained by the FDM±CG and the KLG±CG, respectively, are almost the same, the cost of both aFDM and aKLG being much less than that of the trial control atrial . Contrary to the trial control, the optimal control is almost zero initially but increases to a high value at the ®nal stage to ®t the target velocity ®eld. To assure that the KLG±CG and the FDM±CG yield correct velocity ®eld at the ®nal time, the velocity ®elds obtained by the KLG±CG and the FDM± CG at the ®nal time are compared with the target velocity ®eld at a cross-section of the domain. Fig. 10 shows that the states at the ®nal time, obtained by the FDM±CG (Fig. 10(c) and (d)) and the KLG±CG (Fig. 10(e) and (f)) are indistinguishable from the target velocity ®eld (Fig. 10(a) and (b)). Fig. 11 compares the transient ¯ow ®elds induced by the trial control, Eq. (9), and the optimal control. It is observed that the ¯ow ®eld caused by the trial control experiences larger
Fig. 8. Convergence rate of the conjugate gradient method (the FDM±CG and the KLG±CG).
196
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
Fig. 9. The pro®les of the optimal control obtained either by the FDM±CG or by the KLG±CG for the target velocity ®eld of Fig. 1(b). The trial control, (9), is also displayed for comparison.
excursion than that caused by the optimal control before settling down to the target velocity ®eld, implying that the former incurs larger control cost than the latter. As expected the consumption of CPU time for the case of the KLG±CG is much less than that for the case of FDM±CG since the d.f. for the former is only about 1/2000 of the latter. In fact, when the ultrasparc workstation is used, the CPU time requirement per one iteration of the conjugate gradient method for the FDM±CG and KLG±CG are about 22 h and 1.5 min, respectively. Since the dierence in the d.f. between the low-dimensional dynamic models and the original partial dierential equations shall become much larger as the dimensionality of the problem increases, the reduction of CPU time with the use of the KLG±CG instead of the FDM± CG shall be much more signi®cant in the case of three-dimensional problems. The CPU time requirement to construct the low-dimensional dynamic model is not signi®cant compared with that needed for the FDM±CG. To construct the low-dimensional dynamic model, it is necessary to solve the Navier±Stokes equation with two dierent control function a t to generate snapshots. This step requires 1.2 h. It takes 2 h to obtain the reference velocity ®eld vr . The remaining steps in the Karhunen±Loeve Galerkin procedure consume about 2 h. All in all the construction of the low-dimensional dynamic model requires about 5.2 h, which is much less than the CPU time requirement of just one iteration of the FDM±CG. Moreover, once the low-dimensional dynamic model is constructed, it can be used repeatedly for similar problems. Because the optimization problem under consideration is nonconvex, the proof for existence and uniqueness of the solution is not available in general. But the existence of the optimal control can be assured since the target velocity ®eld has been obtained by solving the Navier±Stokes equation with a given control a t. The uniqueness of the optimal control is con®rmed experimentally by solving the optimization problem using various initial approximation of a t as shown in Fig. 12. Fig. 12(a) shows the convergence rate of the KLG±CG when using dierent initial
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
197
Fig. 10. The comparison of the ®nal state, v x; T , obtained either by the FDM±CG or by the KLG±CG, with the target velocity ®eld vT x. These three velocity pro®les are virtually indistinguishable.
198
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
Fig. 11. The comparison of the transient ¯ow ®elds when a t is given by Eq. (9) and when a t is the optimal control.
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
199
Fig. 12. Convergence of the KLG±CG for various initial approximations of a t: (a) the convergence rate; (b) the converged pro®les of optimal control a t for various initial approximations.
approximations of a t, and the corresponding converged pro®les of optimal control are plotted in Fig. 12(b). This ®gure reveals that the ®nal converged pro®les of the optimal control a t are the same regardless of the initial approximation within numerical error.
7. Conclusion The Karhunen±Loeve Galerkin procedure is applied to the optimal control problem of viscous ¯uid ¯ow through a grooved channel. The Karhunen±Loeve Galerkin procedure, which is a Galerkin method employing the empirical eigenfunctions of the Karhunen±Loeve decomposition as trial functions, limits a priori the function space to the smallest linear subspace that is sucient to describe the observed solutions, and consequently reduces the Navier±Stokes equation to a minimal set of ordinary dierential equations. The low-dimensional dynamic model derived in this
200
H.M. Park, O.Y. Kim / International Journal of Engineering Science 39 (2001) 177±200
way simulates the system accurately for a given space of control functions a t, and can be employed in the solution of optimal control problems eciently. This technique, called the KLG± CG, is compared with the traditional technique of employing the Navier±Stokes equation (called the FDM±CG), and is found to yield accurate results at the drastically reduced computational cost. This is an extension of our previous work where the same technique has been used to solve the optimal control problems of the one-dimensional Burgers equation [1]. Since the reduction of modes through the Karhunen±Loeve Galerkin procedure is much larger in the cases of multidimensional Navier±Stokes equations than in the case of one-dimensional Burgers equation, the KLG±CG is found to be more powerful as regards the computational eciency when applied to the control problems of the Navier±Stokes equation than those of the Burgers equation. Although the importance of control of ¯uid ¯ow has been recognized for a long time, the development of rigorous as well as practical methodology of ¯ow control has been retarded due to the inherent mathematical and computational complexities of the Navier±Stokes equations. In this regard, the Karhunen±Loeve Galerkin procedure is suggested as one of the ecient methods of solving problems of ¯ow control, which is practical enough to be implemented in industrial processes. References [1] H.M. Park, Y.D. Jang, Control of Burgers equation by means of mode reduction, Int. J. Engrg. Sci. 38 (2000) 785± 805. [2] M.D. Gunzburger, L. Hou, T.P. Svobodny, Boundary velocity control of incompressible ¯ow with an application to viscous drag reduction, SIAM J. Control Optim. 30 (1991) 167±181. [3] P. Moin, R.D. Moser, Characteristic-eddy decomposition of turbulence in a channel, J. Fluid Mech. 200 (1989) 471±506. [4] H.M. Park, L. Sirovich, Turbulent thermal convection in a ®nite domain: numerical results, Phys. Fluids A 2 (9) (1990) 1659±1668. [5] N. Aubry, P. Holmes, J.L. Lumley, E. Stone, The dynamics of coherent structures in the wall region of a turbulent boundary layer, J. Fluid Mech. 192 (1988) 115±173. [6] A.E. Deane, I.G. Kevrekidis, G.E. Karniadakis, S.A. Orszag, Low-dimensional models for complex geometry ¯ows: application to grooved channels and circular cylinders, Phys. Fluids 3 (1991) 23±27. [7] P. Moin, T. Bewley, Feedback control of turbulence, Appl. Mech. Rev. 47 (6), part 2, (1994) s3±s13. [8] H.M. Park, D.H. Cho, The use of the Karhunen±Loeve decomposition for the modeling of distributed parameter systems, Chem. Eng. Sci. 51 (1) (1996) 81±98. [9] H.M. Park, M.W. Lee, An ecient method of solving the Navier±Stokes equation for the ¯ow control, Int. J. Numer. Meth. Engrg. 41 (1998) 1131±1151. [10] S.V. Patankar, Numerical Heat Transfer and Fluid Flow, McGraw-Hill, New York, 1980. [11] R. Courant, D. Hilbert, Methods of Mathematical Physics, vol. 1, Interscience, New York, 1953. [12] R. Fletcher, R.M. Reeves, Function minimization by conjugate gradients, The Computer J. 7, 149±154.