A reduced-order finite difference extrapolation algorithm based on POD technique for the non-stationary Navier–Stokes equations

Applied Mathematical Modelling 37 (2013) 5464–5473

Contents lists available at SciVerse ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

A reduced-order finite difference extrapolation algorithm based on POD technique for the non-stationary Navier–Stokes equations q Zhendong Luo a,⇑, Hong Li b,⇑, Ping Sun c, Junqiang Gao a a

School of Mathematics and Physics, North China Electric Power University, Beijing 102206, China School of Mathematical Sciences, Inner Mongolia University, Huhhot 010021, China c School of Mathematics and Computer Science, Guizhou Normal University, Guiyang 550001, China b

a r t i c l e

i n f o

Article history: Received 25 May 2012 Received in revised form 23 September 2012 Accepted 18 October 2012 Available online 7 November 2012 Keywords: Proper orthogonal decomposition technique Finite difference extrapolation algorithm The non-stationary Navier–Stokes equations Error analysis Numerical simulation

a b s t r a c t In this paper, a proper orthogonal decomposition (POD) technique is used to establish a reduced-order finite difference (FD) extrapolation algorithm with lower dimensions and sufficiently high accuracy for the non-stationary Navier–Stokes equations, and the error estimates between the reduced-order FD solutions and the classical FD solutions and the implementation for solving the reduced-order FD extrapolation algorithm are provided. Two numerical examples illustrate the fact that the results of numerical computation are consistent with theoretical conclusions. Moreover, it is shown that the reduced-order FD extrapolation algorithm based on POD method is feasible and efficient for solving the non-stationary Navier–Stokes equations. Ó 2012 Elsevier Inc. All rights reserved.

1. Introduction The finite difference (FD) method is one of most efficient and convenient methods for finding numerical solutions to non-stationary, i.e., time-dependent partial differential equations (TDPDEs) (see, e.g., [1,2]). However, some FD schemes for TDPDEs in practical science and engineering computations usually include many number of degrees of freedom (i.e., unknown quantities). Due to accumulation of truncation errors in the real applied computational process with many number of degrees of freedom, even if a very good FD scheme may also appear no convergence after some computing steps. Thus, an extremely meaningful and very important work is how to establish a reduce-order FD scheme with fewer degrees of freedom and sufficiently high accuracy so that it can alleviate the accumulation of truncation errors, save time required for calculations and memory resource demands in the actual computational process, and continuously simulate the development of fluid flow needed. The proper orthogonal decomposition (POD) method (see [9]) is an efficient means which can reduce degrees of freedom (i.e., unknown quantities) of FD schemes for TDPDEs and alleviate the accumulation of truncation errors in the computational process. The POD method has been widely and successfully applied to numerous fields, including signal analysis and pattern recognition (see [10]), statistics (see [11]), geophysical fluid dynamics and meteorology (also see [11]). The q This work is jointly supported by National Science Foundation of China (11271127, 11061009 and 11061021), Natural Science Foundation of Hebei Province (A2010001663), Natural Science Foundation of Inner Mongolia (2012MS0106), Science Research Program of Guizhou (QKJ[2011]2367), Science Research Program of Inner Mongolia Advanced Education (NJ10006), and Special Funds for Co-construction Project of Beijing and North China Electric Power University. ⇑ Corresponding authors. Tel./fax: +86 10 61772167. E-mail addresses: [email protected] (Z. Luo), [email protected] (H. Li).

0307-904X/$ - see front matter Ó 2012 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.apm.2012.10.051

Z. Luo et al. / Applied Mathematical Modelling 37 (2013) 5464–5473

5465

POD method essentially provides an orthogonal basis for representing the given data in a certain least squares optimal sense, that is, it provides a way to find optimal lower dimensional approximations of the given data. In early times, the POD method was mainly used to perform principal component analysis in computations of statistics and search the main behavior of dynamic systems (see [3–7]). The POD method has been widely applied to reduced-order models for fluid mechanics problems (see, e.g., [8–13]) since the method of snapshots was introduced by Sirovich (see [14]). Some Galerkin POD methods for parabolic problems and a general equation in fluid dynamics have been presented (see [21,22]), and the singular value decomposition (SVD) approach combined with POD technique has been used to treat the Burgers equation (see [23]) and the cavity flow problem (see [24]). Some reduced-order FD schemes for the upper tropical Pacific Ocean formulation, parabolic problems, Burgers equations, the non-stationary Stokes equations, the non-stationary conduction– convection equations, and chemical vapor deposit reaction equations based on POD method have been presented (see [19–25]). Though the long-term stability of POD reduced-order models (see, e.g., [26–29]) has been discussed, almost all existing reduced-order computational methods (see, e.g., [15–25]) based on POD technique are actually to test and verify only the comparison of solutions on the same time span ½0; T and belong to repeating computations. Especially, to the best of our knowledge, there are no published results addressing the case that any reduced-order FD extrapolation algorithm based on POD method for the non-stationary Navier–Stokes equations is established or providing the implementation for solving the reduced-order FD extrapolation algorithm. Therefore, in this paper, we thoroughly improve the existing methods, namely, we do only use the first fewer given numerical solutions on very short time span ½0; T 0  (T 0  T) as snapshots, which are also drawn from samples of the actual physical system trajectories, to construct POD basis and establish reduced-order FD extrapolation algorithm based on POD technique for finding the numerical solutions on total time span ½0; T. Thus, we sufficiently adopt the advantage of POD method, namely, sufficiently utilize the given data (on very short time span ½0; T 0  and T 0  T) to forecast (or infer) future physic phenomenons (on time span ½T 0 ; T), which is why we name this reduced-order method as the reduced-order FD extrapolation algorithm. We also give a criterion of renewing POD basis and the error estimates of solutions of the reduced-order extrapolation algorithm to guide us renewing POD basis, and provide the implementation for solving the reduced-order extrapolation algorithm. We also provide two numerical examples to illustrate that the reduced-order FD extrapolation algorithm based on POD method is feasible and efficient for finding numerical solutions for the non-stationary Navier–Stokes equations. The plan of this paper is organized as follows. Section 2 is to derive a classical FD scheme for the non-stationary Navier–Stokes equations, to generate snapshots from the first fewer transient solutions computed from the equation system derived by the classical FD scheme, to reconstruct optimal orthonormal POD bases from the elements of the snapshots with POD method, and to establish a reduced-order FD scheme with lower dimensions and time second-order accuracy based on POD technique for the non-stationary Navier–Stokes equations. In Section 3, a criterion of renewing POD basis, the error estimates of reduced-order FD extrapolation algorithm, and the implementation for solving the reduced-order FD extrapolation algorithm are provided. In Section 4, two numerical examples are presented illustrating that the results of numerical computation are consistent with theoretical conclusions and validating the feasibility and efficiency of the reduced-order FD extrapolation algorithm. Section 5 provides main conclusions. 2. Classical FD scheme and reduced-order FD extrapolation algorithm for non-stationary Navier–Stokes equations In this section, we establish a classical FD scheme for non-stationary Navier–Stokes equations, extract snapshots from the first fewer transient solutions computed from the equation system derived by the classical FD scheme, reconstruct orthonormal POD bases from the elements of the snapshots, and establish a reduced-order FD extrapolation algorithm with lower dimensions and sufficiently high accuracy based on POD technique for non-stationary Navier–Stokes equations. 2.1. A classical FD scheme for non-stationary Navier–Stokes equations Let X  R2 be a bounded and connected polygonal domain. Consider the following incompressible non-stationary Navier– Stokes equations. Problem I. Find U ¼ ðu; v ÞT and p such that, for T > 0,

! 8 @u u@u v @u @p 1 @2u @2u > > > þ f ; ðx; y; tÞ 2 X  ð0; T; þ þ þ þ ¼ > > @t @x @y @x Re @x2 @y2 > > > > > ! > > > @ v u@ v v @ v @p > 1 @2v @2v > > > < @t þ @x þ @y þ @y ¼ Re @x2 þ @y2 þ g; ðx; y; tÞ 2 X  ð0; T; > > > @u @ v > > þ ¼ 0; > > @x @y > > > > > > > uðx; y; tÞ ¼ u1 ðx; y; tÞ; v ðx; y; tÞ ¼ u2 ðx; y; tÞ; > > > : uðx; y; 0Þ ¼ u0 ðx; yÞ; v ðx; y; 0Þ ¼ v 0 ðx; yÞ;

ðx; y; tÞ 2 X  ð0; T; ðx; y; tÞ 2 @ X  ð0; T; ðx; yÞ 2 X;

ð1Þ

5466

Z. Luo et al. / Applied Mathematical Modelling 37 (2013) 5464–5473

where U ¼ ðu; v ÞT is the velocity vector, p the pressure, T the total time, Re the Reynolds number, f ðx; y; tÞ and gðx; y; tÞ are the given body force at x-direction and y-direction, respectively, and u1 ðx; y; tÞ; u2 ðx; y; tÞ; u0 ðx; yÞ; v 0 ðx; yÞ all are the given functions. Let Dx and Dy be the spatial step increments on x and y direction, respectively, Dt be the time step increment, and unjþ1;k ; v nj;kþ1 , and pnj;k denote the value of function u; v , and p at the point ðxjþ1 ; yk ; t n Þ, ðxj ; ykþ1 ; t n Þ, and ðxj ; yk ; t n Þ 2

2

2

2

ð0 6 j 6 J; 0 6 k 6 K; 0 6 n 6 N ¼ ½T=DtÞ, respectively. Then Problem I has the following classical FD scheme:

  unþ1  unþ1 Dy þ ðv nþ1  v nþ1 ÞDx ¼ 0; jþ1;k j1;k j;kþ1 j;k1 2

2

2

Dt n n  pnj;k Þ þ Dtf jþ1;k ; ðp 2 Dx jþ1;k

¼ unjþ1;k þ F njþ1;k  unþ1 jþ1;k

v nþ1 j;kþ

1 2

¼ v nj;kþ1 þ Gnj;kþ1  u

jþ1;k1 2

2

þ

v j;k1 2v j;kþ1 þv j;kþ3 Dy2

2

2

Dt where F njþ1;k ¼ Re 2

2

2

2

2

n

2

2u

þu 1 jþ1;k jþ ;kþ1 2 2

v

ð3Þ

Dt n  pnj;k Þ þ Dtg nj;kþ1 ; ðp 2 Dx j;kþ1

Dy2

 DDxt unj;kþ1 ð 2

ð2Þ

2

u

þ

n  jþ12;kþ12

2u 1 þu 3 j1;k jþ ;k jþ ;k 2 2 2 Dx2

v

n Þ DDyt j12;kþ12

v

n

ð4Þ

  Dt  DDxt unjþ1;k ðunjþ1;k unj;k Þ DDyt v njþ1;k unjþ1;kþ1 unjþ1;k1 ; Gnj;kþ1 ¼ Re 2

v

n ð j;kþ12

n j;kþ1 

2

v

n j;k Þ.

2

2

2

2

v

2

v

v

2 þ j1;kþ1 j;kþ1 jþ1;kþ1 2 2 2 Dx2

Inserting (3) and (4) into (2) could obtain the approximate

FD scheme of Poisson equation for p as follows

pnj1;k  2pnj;k þ pnjþ1;k

Dx2

þ

pnj;k1  2pnj;k þ pnj;kþ1

Dy 2

¼ R;

ð5Þ

where R ¼½F jþ1;k  F j1;k þ Dtðfjþ1;k  fj1;k Þn =ðDt DxÞ þ ½Gj;kþ1  Gj;k1 þ Dtðg j;kþ1  g j;k1 Þn =ðDtDyÞ. 2 2 2 2 2 2 2  2  If 0:25 junjþ1;k j þ jv nj;kþ1 j DtRe 6 1 and 4Dt 6 minfReDx2 ; ReDy2 g or Dt ¼ oðReDx2 ; ReDy2 Þ, the difference schemes (3)–(5) is 2

2

stable (see [1,2]) and there hold the following error estimates:

    En unjþ1;k ; v nj;kþ1 ; pnj;k ¼ kðuðxjþ1 ; yk ; tn Þ; v ðxj ; ykþ1 ; tn Þ; pðxj ; yk ; t n ÞÞ  unjþ1;k ; v nj;kþ1 ; pnj;k k ¼ OðDt; Dx2 ; Dy2 Þ; 2

2

2

2

2

2

ð6Þ

where n ¼ 1; 2; . . . ; N and k  k denotes the usual normal of the vector. If Re the Reynolds number, body force f ðx; y; tÞ at x-direction and body force gðx; y; tÞ at y-direction, boundary value functions u1 ðx; y; tÞ and u2 ðx; y; tÞ, initial value functions u0 ðx; yÞ and v 0 ðx; yÞ, the time step increment Dt, and the spatial step increments Dx and Dy are given, by solving the FD scheme (3)–(5), we can obtain the FD approximate solutions unjþ1;k ; v nj;kþ1 , and pnj;k ð0 6 j 6 J  1; 0 6 k 6 K  1; 1 6 n 6 NÞ of Problem I. 2

2

Put uni ¼ unjþ1;k ; 2

v ni ¼ v nj;kþ , and pni ¼ pnj;k 1 2

ði ¼ kJ þ j þ 1, 1 6 i 6 m; m ¼ JK; 0 6 j 6 J  1; 0 6 k 6 K  1Þ, respectively. We

may chose the first L solutions to construct a set fuli ; v li ; pli gLl¼1 ð1 6 i 6 m; L  NÞ with L  m elements from the set funi ; v ni ; pni gNn¼1 ð1 6 i 6 mÞ with N  m elements, which are known as snapshots. 2.2. Construct POD bases and establish reduced-order FD extrapolation algorithm The sets of snapshots fuli ; v li ; pli gLl¼1 ð1 6 i 6 mÞ in Section 2.1 can constitute three m  L matrixes Au ¼ ðuli ÞmL ; Av ¼ ðv li ÞmL , and Au ¼ ðpli ÞmL as follows, respectively.

0

u11

B 1 B u2 B Au ¼ B . B .. @ u1m

u21



uL1

u22 .. .

 .. .

uL2 .. .

u2m

   uLm

1

C C C C; C A

0

v 11 v 21 B 1 B v 2 v 22 B



v 1m v 2m



Av ¼ B . B .. @

.. .

 .. .

1

v L1 C v L2 C C

; .. C . C A

v Lm

0

p11

B 1 B p2 B Ap ¼ B . B .. @ p1m

p21



pL1

p22 .. .

 .. .

pL2 .. .

p2m

   pLm

1

C C C C: C A

Since the order m for matrixes Au ATu ; Av ATv , and Ap ATp is far larger than the order L for matrixes ATu Au ; ATv Av , and ATp Ap , however, their positive eigenvalues are identical, therefore, we may first solve the eigenequation corresponding to matrixes e u ¼ rankAu ), kv 1 P kv 2 P    P k e ATu Au ; ATv Av , and ATp Ap to find the eigenvalues ku1 P ku2 P    P k e > 0 ( M e v > 0 (M v ¼ uM u vM e p ¼ rankAp ) and corresponding eigenvectors u ð1 6 j 6 M e u Þ; u ð1 6 j 6 M e v Þ, rankAv ), and kp1 P kp2 P    P k e > 0 ( M uj vj pMp e p Þ, and then by relationship and u ð1 6 j 6 M pj

pffiffiffiffiffiffi /uj ¼ Au uuj = kuj ;

pffiffiffiffiffiffi e u ; / ¼ Av u = kv j ; j ¼ 1; 2; . . . ; M vj vj

pffiffiffiffiffiffi e v ; / ¼ Ap u = kpj ; j ¼ 1; 2; . . . ; M pj pj

e p; j ¼ 1; 2; . . . ; M

e u ), / (j ¼ 1; 2; . . . ; M e v ), and / (j ¼ 1; 2; . . . ; M e p ) corresponding to the nonwe may obtain the eigenvectors /uj (j ¼ 1; 2; . . . ; M vj pj null eigenvalues for matrix Au ATu ; Av ATv , and Ap ATp , respectively.

5467

Z. Luo et al. / Applied Mathematical Modelling 37 (2013) 5464–5473

e u 6 L), M v (0 < M v 6 M e v 6 L), and M p (0 < M p 6 M e p 6 L) columns from three eigenmatrixTaking the first M u (0 < M u 6 M es U u ¼ ð/u1 ; /u2 ; . . ., / e Þ; U v ¼ ð/v 1 ; /v 2 ; . . . ; / e Þ, and U p ¼ ð/p1 ; /p2 , . . . ; / e Þ, respectively, we construct three orthonorv Mv

uM u

pM p

mal POD bases Uu ¼ ð/u1 ; /u2 ; . . . ; /uMu Þ, Uv ¼ ð/v 1 ; /v 2 ; . . ., /v Mv Þ, and Up ¼ ð/p1 ; /p2 ; . . . ; /pMp ÞðM u ; M v ; M p  LÞ (see [19]). Write

unm ¼ ðun1 ; un2 ; . . . ; unm ÞT ;

v nm ¼ ðv n1 ; v n2 ; . . . ; v nm ÞT ;

v ni ¼ v nj;kþ ,

where uni ¼ unjþ1;k ;

1 2

2

pnm ¼ ðpn1 ; pn2 ; . . . ; pnm ÞT ;

n ¼ 1; 2; . . . ; N;

ð7Þ

and pni ¼ pnj;k ði ¼ kJ þ j þ 1, 1 6 i 6 m; m ¼ JK; 0 6 j 6 J  1; 0 6 k 6 K  1Þ, respectively.

Then, there hold that (see [22] or [23])

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kuðMu þ1Þ ; kv lm  Uv UTv v lm k 6 kv ðMv þ1Þ ; kplm  Up UTp plm k 6 kpðMp þ1Þ ; 1 6 l 6 L:

kulm  Uu UTu ulm k 6

ð8Þ

Thus, (3)–(5) are written as the following vector scheme:

 where

nþ1 nþ1 unþ1 m ; v m ; pm

e F

is

T

 T ¼ unm ; v nm ; pnm þ e F ðunm ; v nm ; pnm Þ;

determined

from (3)–(5).

Let

n ¼ 0; 1; 2; . . . ; N  1;

ð9Þ

T  n n n T  um ; v m ; pm ¼ Uu anMu ; Uv bnMv ; Up cnMp ,

and

n n n T un m ¼ ðu1 ; u2 ; . . . ; um Þ ,

n n n T n n n n T v n m ¼ ðv 1 ; v 2 ; . . . ; v m Þ , and pm ¼ ðp1 ; p2 ; . . . ; pm Þ be three column vectors corresponding to u; v , and p, respectively. If n n n n n um ; v m ; pm in (9) are approximately replaced with un m ; v m ; pm (n ¼ 0; 1; 2; . . . ; N), by noting that three matrixes Uu ; Uv , and

Up are formed with the orthonormal vectors, we may obtain the reduced-order FD extrapolation algorithm which has only Mu þ M v þ M p ðM u ; M v ; M p  L  mÞ unknown values:

8 n n T T T < aMu ¼ Uu unm ; bMv ¼ Uv v nm ; cnMp ¼ Up pnm ; n ¼ 1; 2; . . . ; L;  T  T : anþ1 ; bnþ1 ; cnþ1 ¼ an ; bn ; cn e an ; bn ; cn Þ; þ Gð Mu Mp Mu Mp Mv Mv Mv Mu Mp

n ¼ L; L þ 1; . . . ; N  1;

ð10Þ

e an ; bn ; cn Þ = ðUT ; UT ; UT ÞT e F ðUu anMu ; Uv bnMv ; Up cMp Þ. where Gð u p Mv Mu Mp v If we obtain anMu ; bnMv , and cnMp from (10), the reduced-order FD solutions to Problem I are as follows n un m ¼ Uu a M u ;

n v n m ¼ Uv bM v ;

pn m ¼ Up cM p ;

Thus, we obtain that the component forms: 16i6m¼KJ).

n ¼ 1; 2; . . . ; N:

un ¼un i ; jþ12;k

v

n ¼ j;kþ12

v

n i ,

ð11Þ and

pn ¼pn i j;kþ12

(06j6J 1; 06k6K 1; i¼kðJ þ1Þþjþ1;

Remark 1. Since the classical FD scheme (3)–(5) on each time level includes 3m unknown quantities, while the Eqs. (10) and (11) on each time level (when n > L) only include Mu þ M v þ M p unknown quantities (M u ; Mv ; M p  L  m, for first example in Section 4, L ¼ 20; Mu ¼ Mv ¼ Mp ¼ 6, but m ¼ 104 ), so the Eqs. (10) and (11) are a reduced-order FD extrapolation algorithm with very few dimensions. Especially, it is completely deferent from existing reduced-order modelings based on POD technique (see, e.g., [15–25]). Though we draw the snapshots from the first L approximate solutions of classical FD scheme in this paper, in fact, when one computes actual problems, one may obtain the ensemble of snapshots from physical system trajectories by drawing samples from experiments. Since the development and change of a large number of future nature phenomena (for example, weather change) are closely related to previous results, if the physical system of the nature phenomena is well-behaved in the sense that past dynamics are representative and encompassing of future dynamics, one can use the previous or existing results as the snapshots, reconstruct the POD basis, and simulate future developing laws for nature phenomena.

3. Error estimates and implementation of the reduced-order extrapolation algorithm In this section, we first give the error estimates of the reduced-order FD extrapolation algorithm and the criterion of renewing POD basis, and then provide an implementation for solving the reduced-order FD extrapolation algorithm. 3.1. Error estimates and criterion of renewing POD basis It is obvious that the second equation in (10) and (11) has also the following form like (2)–(4)

    nþ1 nþ1 ujþ v nþ1  v nþ1 Dx ¼ 0; 1;k  uj1;k Dy þ j;kþ1 j;k1 2

2

2

L  1 6 n 6 N þ 1;

ð12Þ

2

Dt n n  pn L  1 6 n 6 N þ 1; ðp j;k Þ þ Dtf jþ12;k ; Dx jþ1;k Dt n n ¼ v n þ Gn  pn L  1 6 n 6 N þ 1; ðp j;k Þ þ Dtg j;kþ1 ; j;kþ12  j;kþ12 2 Dx j;kþ1

n nþ1 n ujþ 1;k ¼ ujþ1;k þ F jþ1;k 

ð13Þ

v nþ1 j;kþ

ð14Þ

2

2

1 2

2

5468

Z. Luo et al. / Applied Mathematical Modelling 37 (2013) 5464–5473

Dt where F n jþ1;k ¼ Re

u

2

v j;k1 2v j;kþ1 þv j;kþ3

jþ1;k1 2

2u

jþ1;k 2 Dy2

n

u

þu

jþ1;kþ1 2

þ

j1;k 2

2u

jþ1;k 2 Dx2



þu

jþ3;k 2

n

v   1 2v j;kþ1 þv jþ1;kþ1 n n n Dt n n n Dt j1;kþ2 2 2 ; G DDxt un ðu u Þ v u u ¼ 1 1 1 1 1 1 1 j;kþ jþ1;k j;k Dy jþ ;k Re Dx2 jþ ;k jþ ;kþ jþ ;k 2

2



n 2 2 2 DDxt un v n v n ðv n DDyt v n j;kþ1  v j;k Þ. Dy2 j;kþ12 jþ12;kþ12 j12;kþ12 j;kþ12   0:25 jun jþjv n j DtRe61 and 4Dt6minfReDx2 ;ReDy2 g. jþ1;k j;kþ1

þ

2

The

2

2

2

stability

2

2

condition

of

(12)–(14)

is

2

We also write (10) and (11) as the following vector form.

 nþ1 nþ1 nþ1 T  n n n T n n um ; v m ; pm ¼ um ; v m ; pm þ e F ðun m ; v m ; pm Þ;

n ¼ L; L þ 1; L þ 2; . . . ; N:

ð15Þ

n n T Put en ¼ ðunm ; v nm ; pnm ÞT  ðun m ; v m ; pm Þ . Subtracting (15) from (9) and using (10) yield that

en ¼ ðun ; v n ; pn ÞT  ðUu UTu un ; Uv UTv v n ; Up UTp pn ÞT ; n n F ðunm ; v nm ; pnm Þ  e F ðun enþ1 ¼ en þ e m ; v m ; pm Þ;

n ¼ 1; 2; . . . ; L;

ð16Þ

n ¼ L; L þ 1; L þ 2; . . . ; N  1:

ð17Þ

By (8), we obtain that

ken k 6 Let d ¼ max

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kuðMu þ1Þ þ kv ðMv þ1Þ þ kpðMp þ1Þ ;

n

n ¼ 1; 2; . . . ; L:

ð18Þ

   o junjþ1;k j þ jv nj;kþ1 j DtRe; jun j þ jv n j DtRe; Dt= minfReDx2 ; ReDy2 g . Under the stability conditions of (2)–(4) jþ1;k j;kþ1 2

2

2

2

and (12)–(14), d 6 1. Therefore, there hold from (18), (2)–(4), and (12)–(14) that nL n n kenþ1 k 6 ken k þ k e F ðunm ; v nm ; pnm Þ  e F ðun keL k m ; v m ; pm Þk 6 ð1 þ dÞken k 6    6 ð1 þ dÞ hqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffii nL kuðMu þ1Þ þ kv ðMv þ1Þ þ kpðMp þ1Þ ; L 6 n 6 N  1: 6 ð1 þ dÞ

ð19Þ

Since the absolute of each component for vector is more than its norm, combining (6) with (18) and (19) yields the following results.     Theorem 1. Under the stability conditions 0:25 junjþ1;k j þ jv nj;kþ1 j DtRe 6 1; 0:25 jun j þ jv n j DtRe 6 1, and 0:25Dt 6 jþ1;k j;kþ1 2

2

2

2

minfReDx2 ; ReDy2 g of (2)–(4) and (12)–(14), there hold the following error estimates between the accuracy solution for the non-stationary Navier–Stokes equations and the solutions of the reduced-order FD extrapolation algorithm (10) and (11):

juðxjþ1 ; yk ; t n Þ  un j þ jv ðxj ; ykþ1 ; t n Þ  v n j þ jpðxj ; yk ; t n Þ  pn j;k j jþ12;k j;kþ12 2 2   hqffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffii ¼ O C n ðdÞ kuðMu þ1Þ þ kv ðMv þ1Þ þ kpðMp þ1Þ ; Dt; Dx2 ; Dy2 ; 1 6 n 6 N; where

d ¼ max

n

   o junjþ1;k j þ jv nj;kþ1 j DtRe; jun j þ jv n j DtRe; Dt= minfReDx2 ; ReDy2 g 6 1; C n ðdÞ ¼ 1ð1 6 n 6 LÞ jþ1;k j;kþ1 2

C n ðdÞ ¼ ð1 þ dÞ

nL

ð20Þ

2

2

and

2

ðL þ 1 6 n 6 NÞ.

Remark 2. The error estimates of Theorem 1 give a guide for choosing number of POD basis, namely, as long as taking pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Mu ; M v , and M p such that kuðMu þ1Þ þ kv ðMv þ1Þ þ kpðMp þ1Þ ¼ OðDt; Dx2 ; Dy2 Þ. C n ðdÞ ¼ ð1 þ dÞnL is a criterion for renewing pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi POD basis, namely, if C n ðdÞ½ kuðMu þ1Þ þ kv ðMv þ1Þ þ kpðMp þ1Þ  > max ðDt; Dx2 ; Dy2 Þ, the old POD basis is substituted with the

new POD basis regenerated from new snapshots. Moreover, if n    o n 2 2 max junjþ1;k j þ jv nj;kþ1 j DtRe; jun j þ j v j D tRe 6 D t= minfRe D x ; Re D y g, then d < 1. 1 1 jþ ;k j;kþ 2

2

2

Dt ¼ oðReDx2 ; ReDy2 Þ

and

2

3.2. An implementation of the reduced-order FD extrapolation algorithm Solving the reduced-order FD extrapolation algorithm (10) and (11) makes up of the following 5 steps. Step 1. For given the Reynolds number Re, body force f ðx; y; tÞ at x-direction and body force gðx; y; tÞ at y-direction, boundary value functions u1 ðx; y; tÞ and u2 ðx; y; tÞ, initial value functions u0 ðx; yÞ and v 0 ðx; yÞ, the time step increment Dt, and the spatial step increments Dx and Dy, solving the classical FD scheme (3)–(5) at the first fewer L steps (in usually, take L ¼ 20) yields the FD approximate solutions unjþ1;k ; v nj;kþ1 , and pnj;k ð0 6 j 6 J; 0 6 k 6 K; 1 6 n 6 LÞ and further they 2

2

construct a set of snapshots fuli ; v li ; pli gLl¼1 ð1 6 i 6 m) with L  m elements (for actual engineering problems, the ensemble of snapshots is obtained from physical system trajectories by drawing samples from experiments), where uni ¼ unjþ1;k ; v ni ¼ v nj;kþ1 , and pni ¼ pnj;k ði ¼ kJ þ j þ 1, 1 6 i 6 m; m ¼ JK; 0 6 j 6 J  1; 0 6 k 6 K  1Þ, respectively. 2

2

5469

Z. Luo et al. / Applied Mathematical Modelling 37 (2013) 5464–5473

Step 2. Make up of snapshot matrixes Au ¼ ðuli ÞmL ; Av ¼ ðv li ÞmL , and Ap ¼ ðpli ÞmL and solve the linear systems of equations ðATu Au  ku I L Þuu ¼ 0; ðATv Av  kv I L Þuv ¼ 0, and ðATp Ap  kp I L Þup ¼ 0 obtaining, respectively, the eigenvalues ku1 P ku2 e u ¼rankAu ), kv 1 P kv 2 P    P k e P    P k e > 0 (M e v > 0 ( M v ¼rankAv ), and kp1 P kp2 P    P ku M ep > 0 uM u vM e p ¼rankAp ) and corresponding eigenvectors u ðj ¼ 1; 2; . . . ; M e u Þ; u ðj ¼ 1; 2; . . . ; M e v Þ, and u ðj ¼ 1; 2; . . . ; M e p Þ. (M uj vj pj e u ), M v (M v 6 M e v ), and M p (M p 6 M e p ) of Step 3. For the error l ¼ OðDt; Dx2 ; Dy2 Þ needed, determine the numbers M u (M u 6 M pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi POD bases such that kuðMu þ1Þ þ kv ðMv þ1Þ þ kpðMp þ1Þ 6 l, and construct POD bases Uu ¼ ð/u1 ; /u2 ; . . . ; /uMu ) (where pffiffiffiffiffiffi pffiffiffiffiffiffi /uj ¼ Au uuj = kuj ; j ¼ 1; 2; . . . ; M u Þ; Uv ¼ ð/v 1 ; /v 2 ; . . . ; /v Mv ) (where /v j ¼ Av uv j = kv j ; j ¼ 1; 2; . . . ; M v ), and Up ¼ pffiffiffiffiffiffi ð/p1 ; /p2 ; . . . ; /pMp ) (where /pj ¼ Ap upj = kpj ; j ¼ 1; 2; . . . ; M p ). n n n n n n n Step 4. Solve (10) and (11) obtaining the reduced-order solution vectors un m ¼ ðu1 ; u2 ; . . . ; um Þ; v m ¼ ðv 1 ; v 2 ; . . . ; v m Þ, and n n n n n n n n pn ¼ ðp ; p ; . . . ; p Þ for Problem I, further, obtaining the component forms u ¼ u ; v ¼ v , and p ¼ pn m 1 2 m i i i jþ1;k j;kþ1 j;kþ1 2

2

2

(0 6 j 6 J  1; 0 6 k 6 K  1; i ¼ kðJ þ 1Þ þ j þ 1; 1 6 i 6 m ¼ KJ). n o  pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

j þ jv n j DtRe; minfReDDxt2 ;ReDy2 g . If ð1 þ dÞnL Step 5. Let d ¼ max jun kuðMu þ1Þ þ kv ðMv þ1Þ þ kpðMp þ1Þ 6 l (L þ 1 6 jþ1;k j;kþ1 2

2

n n n n n n n n n n n n 6 N), then un m ¼ ðu1 ; u2 ; . . . ; um Þ; v m ¼ ðv 1 ; v 2 ; . . . ; v m Þ, and pm ¼ ðp1 ; p2 ; . . . ; pm Þ ðn ¼ 1; 2; . . . ; NÞ are just solupffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

nL pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi tions satisfying accuracy needed. Else, namely, if ð1 þ dÞ kuðMu þ1Þ þ kv ðMv þ1Þ þ kpðMp þ1Þ > l (L þ 1 6 n 6 N), l l l l v lm Þ ¼ ðv l1 ; v l2 ;...; v lm Þ, and ðpl1 ;pl2 ;...;plm Þ ¼ ðpl1 ;pl2 ;...;plm Þ ðl ¼ nL; put ðul1 ;ul2 ;...;ulm Þ ¼ ðul 1 ;u2 ;...;um Þ;ðv 1 ; v 2 ;..., n n n n n n n n nL1;...;n1Þ, repeat Step 2 to Step 4 until the solutions un m ¼ ðu1 ;u2 ;...;um Þ; v m ¼ ðv 1 ; v 2 ;...; v m Þ; pm ¼ n n ;p ;...;p Þ ðn ¼ 1;2;...;L;Lþ1;...;NÞ satisfying accuracy needed are obtained. ðpn 1 2 m

Remark 3. Though d in Step 5 is different from that in Theorem 1, it is reasonable in actual computational process. The advantage above implementation of the reduced-order FD extrapolation algorithm consists in that when extrapolation constantly pushes forward process, if the numerical solutions dissatisfy accuracy request, the technique can automatically update POD basis and the reduced-order FD extrapolation algorithm so that it can continuously simulate down and find out the numerical solutions satisfying accuracy requirement. This is incomparable for the existing reduced-order numerical methods (see, e.g., [15–25]), especially, classical FD scheme. 4. Two numerical examples In this section, we present two numerical examples to validate the feasibility and efficiency of the reduced-order FD extrapolation algorithm based on POD technique.

4.1. Example of physical model of square cavity non-stationary flow In this subsection, the physical model is square cavity non-stationary flow. Thus, let the computational field  except that X ¼ ð0; 1Þ  ð0; 1Þ; Re ¼ 103 ; Dx ¼ Dy ¼ 0:01; Dt ¼ 0:0001; u0 ¼ v 0 ¼ g ¼ u1 ¼ u2 ¼ 0, and f ¼ 0 on X f ðx; yÞ ¼ 4yð1  yÞ (if x ¼ 1).

Fig. 1. Classical FD solution (left chart) and reduced-order FD extrapolation algorithm solution (right chart) of the velocity U at the time level t ¼ 1 and when Re ¼ 103 .

5470

Z. Luo et al. / Applied Mathematical Modelling 37 (2013) 5464–5473

Fig. 2. Classical FD solution (left chart) and reduced-order FD extrapolation algorithm solution (right chart) of the pressure field p at the time level t ¼ 1 and when Re ¼ 103 .

Fig. 3. Absolute error for Re ¼ 103 with different number of POD basis and at the time level t ¼ 1.

We find the numerical solutions ðunjþ1;k ; v nj;kþ1 Þ and pnj;kþ1 by classical FD scheme (3)–(5) when n ¼ 10000 (i.e., t ¼ 1), which 2

2

2

are depicted graphically the left charts in Figs. 1 and 2, respectively. When we take L ¼ 20 and use 5 steps in Section 3.2 to solve the reduced-order FD extrapolation algorithm (11) and (12), pffiffiffiffiffiffiffi pffiffiffiffiffiffiffi pffiffiffiffiffiffiffi we obtain from computing that ku7 þ kv 7 þ kp7 6 4  104 . Thus, we only chose the first 6 POD bases, without renewing   POD basis, the reduced-order FD extrapolation algorithm is still convergent, whose solutions un ; v n at t ¼ 1 and pn jþ1;k j;kþ1 j;kþ1 2

2

2

(n ¼ 10000) are depicted graphically the right charts in Figs. 1 and 2. Every two figures in Figs. 1 and 2 are exhibiting quasiidentical similarity, respectively, but the reduced-order FD extrapolation algorithm solution is better and more stable than the classical FD scheme solution. Fig. 3 shows the errors between solutions obtained by reduced-order FD extrapolation algorithm with different number of POD basis and solutions obtained by classical FD scheme when t ¼ 1 and Re ¼ 1000. It also shows that the results for numerical example are consistent with those obtained for the theoretical case. Moreover, it is shown that the reduced-order FD extrapolation algorithm is computationally very effective for finding the approximate solutions of the non-stationary Navier–Stokes equations.

4.2. Example of two dimensional flow past cylinder problem  (30 long and 10 wide). An inlet boundary with a velocLet the cylinder have a radius of 1 in the computational domain X ity of 1 flows parallel to the domain length towards the right of the domain. The center of the cylinder is placed 5 from the inlet boundary. The Reynolds number (Re) is 100. The boundary conditions applied to the cylinder and both lateral sides are no normal flow and zero shear (slip) boundary conditions with a spin-up period of 8. The initial conditions are defined by running the full model from the ‘static’ state during the spin-up (the process of a model adjusting to its forcing) period . ½0; 8. And f ¼ 0 on X  (30 long and 10 wide) into 3000  1000 ¼ 3  106 small squares with side We first divide the computational domain X length Mx ¼ My ¼ 0:01, but refinement mesh is adopted nearly cylinder, and take Mt ¼ 0:001.

Z. Luo et al. / Applied Mathematical Modelling 37 (2013) 5464–5473



5471



We find the numerical solutions unjþ1;k ; v nj;kþ1 and pnj;kþ1 by classical FD scheme (3)–(5) when n ¼ 2000 (i.e., t ¼ 1), which 2

2

2

are depicted graphically the top charts in Figs. 4 and 5, respectively.

Fig. 4. Classical FD solution (top chart) and reduced-order FD extrapolation algorithm solution (bottom chart) of the velocity U at the time level t ¼ 2 and when Re ¼ 102 .

Fig. 5. Classical FD solution (top chart) and reduced-order FD extrapolation algorithm solution (bottom chart) of the pressure field p at the time level t ¼ 2 and when Re ¼ 102 .

5472

Z. Luo et al. / Applied Mathematical Modelling 37 (2013) 5464–5473

When we take L ¼ 20 and use 5 steps in Section 3.2 to solve the reduced-order FD extrapolation algorithm (11) and (12), pffiffiffiffiffiffiffi pffiffiffiffiffiffiffi pffiffiffiffiffiffiffi ku6 þ kv 6 þ kp6 6 3  103 . Thus, we only chose the first 5 POD bases. When n ¼ 1001,

we obtain from computing that

error exceeds 3  103 needing to renew one time POD basis, the reduced-order FD extrapolation algorithm can go on, whose   and pn ; v n at t ¼ 2 (n ¼ 2000) are depicted graphically the bottom charts in Figs. 4 and 5. Every two solutions un jþ1;k j;kþ1 j;kþ1 2

2

2

figures in Figs. 4 and 5 are also exhibiting quasi-identical similarity, respectively, but the reduced-order FD extrapolation algorithm solution is better and more stable than the classical FD scheme solution. The error between solutions obtained by reduced-order FD extrapolation algorithm with 5 POD basis and solutions obtained by classical FD scheme when t ¼ 2 and Re ¼ 100 does no exceed 2  103 . It also shows that the results for numerical example are consistent with those obtained for the theoretical case. Moreover, it is shown that the reduced-order FD extrapolation algorithm is also computationally very effective for finding the long-term approximate solutions of the non-stationary Navier–Stokes equations. 5. Conclusions and perspective In this paper, we have employed the SVD method and the POD technique to establish a reduced-order FD extrapolation algorithm for the non-stationary Navier–Stokes equations. We first compile ensembles of data from the first few L ðL  NÞ transient solutions computing an equation system derived with the classical FD scheme for the non-stationary Navier– Stokes equations, while in actual applications, one may obtain the ensemble of snapshots from physical system trajectories by drawing samples from experiments with interpolation (or data assimilation). Next, we employ the SVD method to deal with ensembles of data obtaining the POD basis. And then, the classical FD solution vectors are replaced with the linear combination of the most main POD basis to establish a reduced-order FD extrapolation algorithm for non-stationary Navier– Stokes equations. Finally, we provide the error estimates between the classical FD solutions and the reduced-order FD extrapolation algorithm solutions based on POD method and the implementation for solving the reduced-order FD extrapolation algorithm of non-stationary Navier–Stokes equations. Comparing with the theoretical errors and numerical computational error shows that the results for numerical examples are consistent with those obtained for the theoretical case. It is shown that our present method has improved and innovated the existing methods (see, e.g., [15–25]). The reduced-order FD extrapolation algorithm here does only use the first fewer given numerical solutions on very short time span ½0; T 0  (T 0  T) as snapshots to construct POD basis and establish reduced-order FD extrapolation algorithm based on POD technique for finding the numerical solutions on total time span ½0; T. Thus, the advantages of POD method are sufficiently adopted, namely, it sufficiently utilizes the given data (on very short time span ½0; T 0  and T 0  T) to forecast (or infer) future physic phenomenons (on time span ½T 0 ; T), which is an extremely meaningful and very useful work. In this paper, we do only present the reduced-order FD extrapolation algorithm for simpler staggered mesh Marker and Cell FD scheme (see [1,2]), but its ideas and technique can be easily extended other FD schemes and more complicated PDEs. Future work in this area will aim to extend the reduced-order FD extrapolation algorithm, implementing it for a realistic atmosphere quality forecast system and more complicated PDEs. Acknowledgments The authors acknowledge cordially the anonymous referees and Editor M. Cross for their valuable comments which led to the improvement of this paper. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15]

T. Chung, Computational Fluid Dynamics, Cambridge University Press, Cambridge, 2002. X.K. Xin, R.X. Liu, B.C. Jiang, Computational Fluid Dynamics, National Defense Science Technology Press, Changsha, 1989 (in Chinese). P. Holmes, J.L. Lumley, G. Berkooz, Turbulence, Coherent Structures, Dynamical Systems and Symmetry, Cambridge University Press, Cambridge, 1996. K. Fukunaga, Introduction to Statistical Recognition, Academic Press, New York, 1990. I.T. Jolliffe, Principal Component Analysis, Springer-Verlag, Berlin, 2002. 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 Dynam. 192 (1988) 115–173. J.L. Lumley, Coherent structures in turbulence, in: R.E. Meyer (Ed.), Transition and Turbulence, Academic Press, New York, 1981, pp. 215–242. H.V. Ly, H.T. Tran, Proper orthogonal decomposition for flow calculations and optimal control in a horizontal CVD reactor, Quart. Appl. Math. 60 (2002) 631–656. A.J. Majda, I. Timofeyev, E. Vanden–Eijnden, Systematic strategies for stochastic mode reduction in climate, J. Atmos. Sci. 60 (2003) 1705–1723. P. Moin, R.D. Moser, Characteristic-eddy decomposition of turbulence in channel, J. Fluid Mech. 200 (1989) 417–509. M. Rajaee, S.K.F. Karlsson, L. Sirovich, Low dimensional description of free shear flow coherent structures and their dynamical behavior, J. Fluid Mech. 258 (1994) 1401–1402. R.D. Roslin, M.D. Gunzburger, R.A. Nicolaides, G. Erlebacher, M.Y. Hussaini, A self-contained automated methodology for optimal flow control validated for transition delay, AIAA J. 35 (1997) 816–824. F. Selten, Baroclinic empirical orthogonal functions as basis functions in an atmospheric model, J. Atmos. Sci. 54 (1997) 2100–2114. L. Sirovich, Turbulence and the dynamics of coherent structures: Part I–III, Quart. Appl. Math. 45 (3) (1987) 561–590. K. Kunisch, S. Volkwein, Galerkin proper orthogonal decomposition methods for parabolic problems, Numer. Math. 90 (2001) 117–148.

Z. Luo et al. / Applied Mathematical Modelling 37 (2013) 5464–5473

5473

[16] K. Kunisch, S. Volkwein, Galerkin proper orthogonal decomposition methods for a general equation in fluid dynamics, SIAM J. Numer. Anal. 40 (2002) 492–515. [17] K. Kunisch, S. Volkwein, Control of Burgers equation by a reduced order approach using proper orthogonal decomposition, J. Optim. Theory Appl. 102 (1999) 345–371. [18] D. Ahlman, F. Södelundon, J. Jackson, A. Kurdila, W. Shyy, Proper orthogonal decomposition for time-dependent lid-driven cavity flows, Numer. Heat Trans. Part B – Fund. 42 (2002) 285–306. [19] Z.D. Luo, J. Chen, J. Zhu, R.W. Wang, I.M. Navon, An optimizing reduced order FDS for the tropical Pacific Ocean reduced gravity model, Int. J. Numer. Meth. Fluids 55 (2) (2007) 143–161. [20] P. Sun, Z.D. Luo, Y.J. Zhou, Some reduced finite difference schemes based on a proper orthogonal decomposition technique for parabolic equations, Appl. Numer. Math. 60 (2010) 155–164. [21] Z.D. Luo, X.Z. Yang, Y.J. Zhou, A reduced finite difference scheme based on singular value decomposition and proper orthogonal decomposition for Burgers equation, J. Comput. Appl. Math. 229 (1) (2009) 97–107. [22] Z.D. Luo, Q.L. Ou, Z.H. Xie, A reduced finite difference scheme and error estimates based on POD method for the non-stationary Stokes equation, Appl. Math. Mech. 32 (7) (2011) 847–858. [23] Z.D. Luo, R.W. Wang, J. Zhu, Finite difference scheme based on proper orthogonal decomposition for the non-stationary Navier-Stokes equations, Sci. China Ser. A: Math. 50 (8) (2007) 1186–1196. [24] P. Sun, Z.D. Luo, Y.J. Zhou, A reduced finite difference scheme based on POD for the non-stationary conduction–convection equations, Math. Numer. Sin. 31 (3) (2009) 234–323. [25] J. Du, J. Zhu, Z.D. Luo, I.M. Navon, An optimizing finite difference scheme based on proper orthogonal decomposition for CVD equations, Int. J. Numer. Meth. Biomed. Eng. 27 (2011) 78–94. [26] A.E. Deane, I.G. Kevrekidis, G.E. Karniadakis, S.A. Orsag, Low-dimensional models for complex geometry flows: Application to grooved channels and circular cylinder, Phys. Fluids A 3 (10) (1991) 2337–2354. [27] X. Ma, G.E. Karniadakis, A low-dimensional model for simulating three- dimensional cylinder flow, J. Fluid Mech. 458 (2002) 181–190. [28] M. Bergmann, C.H. Bruneau, A. Iollo, Enablers for robust POD models, J. Comput. Phys. 228 (2) (2009) 516–538. [29] Z. Wang, I. Akhtar, J. Borggaard, T. Iliescu, Two-level discretizations of non-linear closure models for proper orthogonal decomposition, J. Comput. Phys. 230 (2011) 126–146.