Numerical integration of fluid dynamics problems by discrete dynamical systems

Chaos, Solitons and Fractals 11 (2000) 193±206

www.elsevier.nl/locate/chaos

Numerical integration of ¯uid dynamics problems by discrete dynamical systems F.N. Koumboulis a,*, M.G. Skarpetis b, B.G. Mertzios c b

a University of Thessaly, Department of Mechanical and Industrial Engineering, Pedion Areos, 383 34 Volos, Greece National Technical University of Athens, Department of Electrical and Computer Engineering, Division of Electroscience, Greece c Democritus University of Thrace, Department of Electrical and Computer Engineering, 67100 Xanthi, Greece

Abstract A discrete time algorithm for the numerical solution of the Navier±Stokes equations, namely the set of nonlinear partial di€erential equations (PDEs) characterizing the ¯ow of viscous and compressible ¯uids is presented. The good performance of the algorithm is guaranteed by preserving the passivity property of both the continuous time and discrete multidimensional dynamical descriptions of the Navier±Stokes equations. The present results are applied to compute the steady ¯ows between parallel plates (Couette shear ¯ows). Ó 1999 Elsevier Science Ltd. All rights reserved.

1. Introduction The problem of deriving the analytic solution of nonlinear (NL) partial di€erential equations (PDEs) is a rather dicult task. The problem has not as yet been solved for the majority of cases. For this reason a variety of numerical methods solving discretized versions of the NL PDEs (see for example [1±8] and the references therein). The solution of such a numerical problem appears to have many applications in ¯uid mechanics, electromagnetics and thermodynamics. One of the most important set of NL PDEs is that of the Navier±Stokes equations. The present paper is devoted towards the aim of deriving the numerical solution of this type of equations. Motivated by the results in [9,10] regarding the solution of NL PDEs (Maxwell equations, Euler equations) via multidimensional (MD) passive wave digital ®lters (WDF), the method in [11] is modi®ed to cover the solution dynamic problems for ¯uids with nonlinear viscosity. In particular, a discrete algorithm for the numerical solution of the Navier±Stokes equations, namely the set of nonlinear PDEs characterizing the ¯ow of viscous and compressible ¯uids is presented. The numerical solution is achieved by directly modeling the original continuous-domain physical system by mean of a discrete multidimensional passive (MD-passive) dynamical system, using principles of MD nonlinear digital ®ltering (see for example [11] and [12]). The resulting integration algorithm is highly robust and massively parallel, implying only local interconnections. The WDF characteristics of the discretization guarantee the robustness of the algorithm. In particular [15], in a WDF algorithm all rounding/truncation and over¯ow corrections can be carried out in such a way that passivity and even incremental passivity [16] hold even under ®nite-arithmetic conditions, thus ensuring full robustness [17] of the algorithm. The above results are applied to compute the steady ¯ows between parallel plates (Couette shear ¯ows). Here the viscosity coecient is assumed to be a nonlinear but monotone function of space as well as that the ¯uid is polytropic. Clearly, this case covers a

*

Corresponding author.

0960-0779/00/$ - see front matter Ó 1999 Elsevier Science Ltd. All rights reserved. PII: S 0 9 6 0 - 0 7 7 9 ( 9 8 ) 0 0 2 8 4 - 7

194

F.N. Koumboulis et al. / Chaos, Solitons and Fractals 11 (2000) 193±206

variety of dynamic ¯uid phenomena. It is important to mention that some ®rst results regarding the solution of the Navier±Stokes equations for polytropic ¯uids with constant viscosity, via MD WDFs, have been presented in [13], while the respective results for inviscous ¯uids in [11] and [18].

2. Equivalent expressions of the Navier±Stokes equations Consider a viscous compressible ¯uid with monotone viscosity. The change in the momentum of the ¯uid ¯owing through a control volume at any time instant, as well as the dynamics expressing the conservation of mass are governed by the following two partial di€erential equations [14]:   oq ‡ q
rq ˆ qf ÿ rp ÿ x0 ; …2:1a† q o~t oq ‡ r
qq ˆ 0; o~t

…2:1b†

T respectively. The vector ~tt ˆ ‰t1 t2 t3 ~tŠ is the vector having as elements the three spatial variables x, y, z and the time ~t…t1 ˆ x; t2 ˆ y; t3 ˆ z†. Based upon this de®nition the well known nabla operator r is expressed equivalently as follows:  T o o o : rˆ ot1 ot2 ot3 T

T

The velocity vector is denoted by q ˆ ‰q1 q2 q3 Š , q is the density, f ˆ ‰f1 f2 f3 Š is the vector of external forces applied to the ¯uid inside the control volume and p is the pressure of the ¯uid in the control volume. T The vector x0 ˆ ‰x01 x02 x03 Š is the vector of forces due to viscosity. In particular for a viscous ¯uid ¯ow the vector x0 is l x0 ˆ ÿlr2 q ÿ r…r
q† ÿ Al ‰q1 q2 q3 ŠT ; …2:1c† 3 where l is the viscosity coecient (being a function of space and time) and where the operator Al is de®ned to be 24 3 l o ‡ l2 oto2 ‡ l3 oto3 ÿ 23 l1 oto2 ‡ l2 oto1 ÿ 23 l1 oto3 ‡ l3 oto1 3 1 ot1 6 7 2 o o 4 l o ‡ l3 oto3 ‡ l1 oto1 ÿ 23 l2 oto3 ‡ l3 oto2 7 Al ˆ 6 …2:2a† 3 2 ot2 4 ÿ 3 l2 ot1 ‡ l1 ot2 5; 4 ÿ 23 l3 oto1 ‡ l1 oto3 ÿ 23 l3 oto2 ‡ l2 oto3 l o ‡ l2 oto2 ‡ l1 oto1 3 3 ot3 l1 ˆ

ol ; ot1

l2 ˆ

ol ; ot2

l3 ˆ

ol : ot3

…2:2b†

In this section three equivalent expressions of the Navier±Stokes equations (2.1) will be presented. These equivalent expressions will facilitate the passive discretization algorithm that will be derived. It is mentioned that the steps required for the derivation of the equivalent forms are analogous to those in [10]. For presentation purposes the three equivalent expressions will be partitioned into three di€erent subsections. 2.1. Generalized time It is convenient for the discretization steps, that follow, to transform the time ~t to a new variable t4 having dimension of a spatial variable [11]. This new variable is de®ned by means of a new function q4 equal to the derivative q4 ˆ q4 …~t† ˆ dt4 =d~t and the new variable t4 . Similarly the vector ~tt is transformed to a T generalized time vector t ˆ ‰t1 t2 t3 t4 Š (including space and time). The time variable t4 can be considered as an arbitrary, for the time being, variable which will be speci®ed later in order to facilitate the proof passivity. Using this new vector and the notation

F.N. Koumboulis et al. / Chaos, Solitons and Fractals 11 (2000) 193±206 T

D ˆ ‰D1 D2 D3 D4 Š ;

Dk ˆ

o …k ˆ 1; 2; 3; 4†; otk

195

…2:3†

the equations in (2.1a) and (2.1b) can be rewritten equivalently as follows: rp ‡

4 X

qqk Dk q ˆ qf ÿ x0 ;

…2:4a†

kˆ1 3 X Dk …qqk † ‡ q4 D4 q ˆ 0:

…2:4b†

kˆ1

The new time variable introduced above does not in¯uence the Eq. (2.1c). The de®nition (2.3) could easily be used to reform, in a way analogous to (2.2), the Eq. (2.1c). 2.2. Normalization The goal of this transformation is to introduce new arbitrary parameters in the Navier±Stokes equations. These parameters will be useful in proving passivity. Consider the identity p ÿ p
1 …2:5a† xD y x ˆ xDy ‡ yDx 2 holding true for every real functions x, y with x scalar (x P 0) and y vector. The above identity is usually called normalization identity. De®ne the following normalized quantities p  T 1 t20 2q4 a a p; l^ ˆ l; …2:5b† q^ ˆ q^1 q^3 q^3 ˆ p q; q^ ˆ q; p^ ˆ p0 p0 p0 q4 a p  T 2 a f q; f^ ˆ f^1 f^3 f^3 ˆ p0

p  T 2 a ^ 02 x ^ 03 ˆ ^0 ˆ x ^ 01 x x x0 ; p0

…2:5c†

where t0 and p0 are positive constants and a ˆ a…~t† is an arbitrary function of time. Using (2.5), Eqs. (2.4a) and (2.4b) take on the form  r 4 r  X qk q4 qk q4 1 ^ ^ 0j ; j ˆ 1; 2; 3 fj ÿ x …2:6a† q^a 2 Dk q^j q^a 2 ‡ Dj p^ ˆ 2 t0 t0 kˆ1   3 X qk ˆ 0: Dk q^ q4 kˆ1

…2:6b†

^ 0 is now given by the following According to Eqs. (2.5c) and (2.1c), the transformed viscosity forces vector x formula l^ T ^ 0 ˆ ÿ^ lr2 q^ ÿ r…r
q^† ÿ Al^ ‰^ q1 q^2 q^3 Š ; x 3 where l^ is the normalized viscosity coecient and where 2 4 ÿ 23 l^1 D2 ‡ l^2 D1 l^ D1 ‡ l^2 D2 ‡ l^3 D3 63 1 2 4 l^ D ‡ l^3 D3 ‡ l^1 D1 Al^ ˆ 4 ÿ 3 l^2 D1 ‡ l^1 D2 3 2 2 ÿ 23 l^3 D1 ‡ l^1 D3 ÿ 23 l^3 D2 ‡ l^2 D3 l^1 ˆ

o^ l ; ox

l^2 ˆ

o^ l ; oy

l^3 ˆ

o^ l : oz

…2:6c† 3 ÿ 23 l^1 D3 ‡ l^3 D1 7 ÿ 23 l^2 D3 ‡ l^3 D2 5; 4 l^ D ‡ l^2 D2 ‡ l^1 D1 3 3 3

…2:7a†

…2:7b†

196

F.N. Koumboulis et al. / Chaos, Solitons and Fractals 11 (2000) 193±206

2.3. Hadamard transformation Due to conservation of energy, physical systems are usually passive. The systems of PDEs describing a physical system does usually not re¯ect this passivity. The property of passivity is necessary for the discretization of the PDEs and the derivation of their numerical solution. This hold true since a passive MD circuit yield a causal MD WDF which is robust and stable [11]. In particular, in order to obtain a MD passive circuit the following Hadamard transformation is applied t0 ˆ H ÿ1 t; where t0 ˆ

…2:8†

T ‰t10 t20 t30 t40 Š

and where H is a suitable orthogonal matrix de®ned by [10] 3 1 ÿ1 ÿ1 1 1 6 ÿ1 1 ÿ 1 17 7: H ˆ HT ˆ 6 4 1 15 2 ÿ1 ÿ 1 1 1 1 1 2

…2:9†

The above transformation is a transformation of the coordinate system (including time). De®ning the operators D0k ˆ

o otk0

…k ˆ 1; 2; 3; 4†

…2:10†

and using relations (2.9) and (2.10), the following expression is derived D0 ˆ H T D:

…2:11†

From Eqs. (2.9) and (2.11) the operators D0i may further be expressed in terms, of new derivative operators as follows


1ÿ 1ÿ 1ÿ D1 ˆ D001 ÿ D004 ; D2 ˆ D002 ÿ D005 ; D3 ˆ D003 ÿ D006 ; …2:12† 2 2 2 D4 ˆ


1ÿ
1ÿ
1 ÿ 00 D1 ‡ D004 ˆ D002 ‡ D005 ˆ D003 ‡ D006 2 2 2

D001 ˆ D04 ‡ D01 ; D006 ˆ D01 ‡ D02 :

D002 ˆ D04 ‡ D02 ;

D003 ˆ D04 ‡ D03 ;

…2:13† D004 ˆ D02 ‡ D03 ;

D005 ˆ D01 ‡ D03 ; …2:14†

The above relations can be considered as a de®nition of the operators D00i …i ˆ 1; 2; 3; 4† in terms of the operators D0i …i ˆ 1; 2; 3; 4†. Using the above transformations and after some extensive manipulations Eqs. (2.6) and (2.7) may be rewritten as follows: 7 X ^ 0j ; uj;i ˆ f^j ÿ x

j ˆ 1; 2; 3

…2:15a†

iˆ1

10 X u4;i ˆ 0;

…2:15b†

iˆ1

where the variables ui;j , for j ˆ 1; 2; 3, are related to the normalized velocities as follows p p  uj;i ˆ Li D0i Li q^j ; i ˆ 1; . . . ; 5 D05 ˆ D4   uj;6 ˆ D00j q^j ‡ p^ ;

  uj;7 ˆ D00j‡3 q^j ÿ p^

F.N. Koumboulis et al. / Chaos, Solitons and Fractals 11 (2000) 193±206

for j ˆ 4, as follows p p  L0;i p^ ; u4;i ˆ L0;i D0i u4;i‡4 ˆ ui;6 ;

197

i ˆ 1; 2; 3; 4

u4;i‡7 ˆ ÿui;7 ;

i ˆ 1; 2; 3

and where Li ˆ q^a

u0;i q4 ÿ e; u20

i ˆ 1; 2; 3; 4

L5 ˆ 2…e ÿ 1† L0;i ˆ 2

q^t ÿ 1 K1 u^0;i ÿ 3 ‡ 2 ; p^ p^

where K1 and e arbitrary real constants with e P 1 and where u0;1 ˆ q1 ÿ q2 ÿ q3 ‡ q4 ; u0;2 ˆ ÿq1 ‡ q2 ÿ q3 ‡ q4 ; u0;3 ˆ ÿq1 ÿ q2 ‡ q3 ‡ q4 ; u0;4 ˆ q1 ‡ q2 ‡ q3 ‡ q4 u0;i u^0;i ˆ p ; a q4

1 tˆ p

Z

dp p0 ˆ 2 q^ t0 p

Z

dp ; q

where for the derivation of the above relations it has been assumed that the ¯uid is barotropic, i.e. that the pressure is a unique function of the density. Clearly this assumption covers the majority of ¯uid mechanics phenomena. The variables ui;j , being a nonlinear transformation of the ¯uid velocities are the unknowns in the nonlinear system of Eq. (2.15).

3. Generalized Kirchho€ circuit realization It can easily be observed that the transformed form of the Navier±Stokes equations, namely the four equations in (2.15) describe a generalized Kirchho€ circuit, [10], involving four main loops and 27 supplementary loops, appropriately interconnected. 3.1. Main loop The circuit involving the four main loops is shown in Fig. 1. The circuit in Fig. 1 is an obvious realization of the four equations in (2.15). Each equation corresponds to one Kirchho€ loop. The generalized voltages ^ 0j are the outputs of the supplementary circuits that will be realized later. Clearly, the overall generalized x Kirchho€ circuit is passive if the above circuit (Fig. 1) involving the four main loops is passive as well as if every supplementary circuit is also passive. The passivity of the circuit of Fig. 1, is not guaranteed. For the four main loops of the circuit of Fig. 1, passivity is guaranteed if the following conditions are satis®ed [12] L1 ; . . . ; L5 P 0;

L01 ; . . . ; L05 P 0:

…3:1†

To analyze the above conditions consider the case of a polytropic ¯uid, i.e. a ¯uid with pressure depending upon the density by a relation of the form: p ˆ Cqn where C; n are positive constant real numbers with n > 1. Using the assumption that the ¯uid is polytropic and after some algebraic manipulations, the conditions in Eq. (3.1) can be expressed in form of denormalized variables as in the following Lemma [10].

198

F.N. Koumboulis et al. / Chaos, Solitons and Fractals 11 (2000) 193±206

Lemma 3.1. If the following inequalities are satisfied p 3…n ÿ 1†pmax p0 q p 6 6 min …q4 ÿ 3†qmax aq4 e 2q4 ÿ 2 3qmax

…3:2†

where pmax ˆ maxp…t†; t1 ;t2 ;t3

qmax ˆ maxq…t†; t1 ;t2 ;t3

qmax ˆ maxq…t† t1 ;t2 ;t3

^ 0j are passive, then the overall generalized Kirchhoff circuit in and the circuits corresponding to the voltages x Fig. 1 is passive. Lemma 3.2. The inequality in Eq. (3.2) is satisfied if the generalized derivative of time is chosen to be s 3…n ÿ 1†pmax p …3:3† ‡ 3qmax ; q4 P 2qmin while the arbitrary parameters a and p0 are chosen to be r p0 3…n ÿ 1†pmax ˆ q4 qmin pmax : a 2e

…3:4†

Proof. To prove the above relations, substitute Eqs. (3.3) and (3.4) into Eq. (3.2), to yield inequalities identically true.  ^ 0j 3.2. Circuit realizations for the generalized voltages x Till the end of this section the 27 supplementary generalized Kirchho€ circuits will be derived, while their passivity will be certi®ed. Note that the realization of these circuits will be made on the basis of four characteristics circuits given in Figs. 2±5 that follows. These circuits involve the nonlinear viscosity terms thus realizing the terms of the Navier±Stokes equation not included in the Euler equations. To realize these circuits new type of loops, involving gyrators as standard type of approximations, will be introduced.

Fig. 1. Generalized Kirchho€ circuit (Main loops). ^ 0j2 ; x ^ 0j3 . ^ 0j1 ; x Fig. 2. MD passive Kirchho€ circuit for the description of x

F.N. Koumboulis et al. / Chaos, Solitons and Fractals 11 (2000) 193±206

199

^ 00j . Fig. 3. MD passive Kirchho€ circuit for the description of x Fig. 4. MD passive Kirchho€ circuit for the description of x11 .

^ 02 x ^ 03 ŠT can be considered as the sum of three di€erent ^ 01 x According to Eq. (2.6c) the viscosity forces ‰x terms, i.e. the sum of the terms: ÿ^ lr2 q^;

l^ ÿ r…r
q^†; 3

ÿAl^ ‰^ q1

q^2

T q^3 Š :

…3:5†

Each of the ®rst two terms can not be realized by an MD-passive circuit and thus it has to be modi®ed. Following the method in [13] ten new circuits can be derived to approximate describe the above two terms. These circuits can be considered as modi®cations of the standard Kirchho€ realizations guaranteeing passivity. The modi®ed circuits are achieved by deleting the inductance L5 from the reference circuit in T ^j ˆ x ^ 0j (for ^ 0j ˆ ‰ÿ^ l=3†r…r
q^† ÿ Al^ ‰q1 ; q2 ; q3 Š Šj (j ˆ 1; 2; 3), by x lr2 q^ ÿ …^ Fig. 1, replacing, each x 2 symbolic simpli®cation) (where ‰ÿ^ lr q^ ÿ …^ l=3†r…r
q^† ÿ Al^ ‰q1 ; q2 ; q3 ŠT Šj is the jth element of ÿ^ lr2 q^ ÿ …^ l=3†r…r
q^† ÿ Al^ ‰q1 ; q2 ; q3 ŠT ) and using the following lemma. Lemma 3.3. The first two terms in Eq. (2.6c), namely the expression ‰ÿ^ lr2 q^ ÿ …^ l=3†r…r
q^†Šj can be ap0 0 0 00 ^ j2 ; x ^ j3 ; x ^ j as follows ^ j1 ; x proximated by the voltages x # " l^ 2 ^ 0j1 ‡ x ^ 0j2 ‡ x ^ 0j3 ‡ x ^ 00j ; …3:6† ÿ l^r q^ ÿ r…r
q^† ˆ x 3 j

^ 0j1 ; x ^ 0j2 ; x ^ 0j3 and x ^ 00j are the outputs of the circuits in Figs. 2 and 3, respectively, where where the voltages x 0 0 2 00 c ˆ l^=2^ s , d ˆ l^=…^ s † , a^ ˆ l^=2^ s , b ˆ 3^ l=^ s00 , s^0 ˆ u0 s0 , s^00 ˆ u0 s00 and where s0 is a time constant. Proof. Solve the circuits in Figs. 2 and 3 with respect to their outputs. Substituting the resulting relations in Eq. (3.6) an approximate equality is derived.  Note that the respective result for the special case of constant normalized viscosity coecient l^ has been reported in [13]. T In order to implement the third term in Eq. (2.6c), namely the term ÿAl^j ‰q^1 q^2 q^3 Š , i.e. the jth element T of ÿAl^ ‰q^1 q^2 q^3 Š , can be expressed as the sum three variables h iT …3:7† ÿAl^;j q^1 q^2 q^3 ˆ ÿxj1 ÿ xj2 ÿ xj3 ;

200

F.N. Koumboulis et al. / Chaos, Solitons and Fractals 11 (2000) 193±206

where

8 > < 43 l^1 D1 ‡ l^2 D2 ‡ l^3 D3 xj1 ˆ ÿ 23 l^2 D1 ‡ l^1 D2 > : ÿ 2 l^3 D1 ‡ l^1 D3 3

if j ˆ 1 if j ˆ 2 ; if j ˆ 3

8 > < ÿ 23 l^1 D3 ‡ l^3 D1 xj3 ˆ ÿ 23 l^2 D3 ‡ l^3 D2 > : 4 l^3 D3 ‡ l^2 D2 ‡ l^3 D3 3

if j ˆ 1 if j ˆ 2 ; if j ˆ 3

8 > < ÿ 23 l^1 D2 ‡ l^2 D1 xj2 ˆ 43 l^2 D2 ‡ l^1 D1 ‡ l^3 D3 > : ÿ 2 l^3 D2 ‡ l^2 D3 3

if j ˆ 1 if j ˆ 2 if j ˆ 3

…3:8†

For j ˆ 1 the ®rst term, i.e. the term x11 can be described via an MD passive circuit as in Fig. 4. The second term, i.e. the term x12 can be described with an MD passive circuit (inductance, capacitor, gyrator) as in Fig. 5. Clearly, a circuit analogous to that in Fig. 5 can easily be constructed for x13 . In order to guarantee passivity, in all three circuits for x11 ; x12 and x13 , the respective derivative viscosity coecients l^i …i ˆ 1; 2; 3† must be positive. This can easily be translated to the property of increasing monotone normalized viscosity coecient l^ with respect to space ‰x; y; zŠ. Using Eqs. (2.5b) and (3.4), it can readily be observed that the original viscosity coecient l must be increasing monotone with respect to space ‰x; y; zŠ. In the introduction it has already noted that the ¯uid has been assumed to have monotone viscosity. Since the positioning and orientation of the space coordinate system can be chosen by the analyst as well as since the circuits in Figs. 4 and 5 can be easily be re-realized to have opposite sign, the case of decreasing monotone viscosity is considered to be covered by the case of increasing monotone viscosity. The circuits in Figs. 4 and 5 are designed in order to describe the variables x11 ; x12 and x13 . Analogous circuitry can readily be designed for x21 ; x22 , x23 and x31 ; x32 , x33 , thus completing the overall generalized Kirchho€ circuit. We are now in position to establish the following proposition. Proposition 3.1. The viscosity forces term in the Navier±Stokes equations can be approximated by the gen^ 0j2 ; x ^ j3 ; x ^ 00j ; xj1 ; xj2 ; xj3 (see Figs. 2±5), ^ 0j1 ; x eralized Kirchhoff circuits, realizing the generalized voltages x as follows ^ 0j1 ‡ x ^ 0j2 ‡ x ^ 0j3 ‡ x ^ 00j ÿ xj1 ÿ xj2 ÿ xj3 ^j ˆ x x

…3:9†

Combining the realization in Fig. 1 with the results of Proposition 3.1, the following theorem is established. Theorem 3.1. The Navier±Stokes equations, for a monotone viscosity fluid, presented in (2.1) can approximately be realized by a passive generalized Kirchhoff circuit. The circuit is given in Fig. 1 where the subcircuits ^ 0j are of the form presented in Figs. 2±5, according to Proposition 3.1. having output voltages x

4. Discretization algorithm In order to derive a robust algorithm it is necessary to adopt power waves instead of the voltage waves usually preferred in wave digital ®ltering [12]. Thus for a port of voltage u, current i and non-constant port resistance r P 0, the forward wave a and the backward wave b are de®ned as follows [12] aˆ

u ‡ ri p ; 2 r

bˆ

u ÿ ri p : 2 r

…4:1†

In order to transform the circuit of Fig. 1 to an equivalent reference WDF circuit, we apply the general trapezoidal rule (see for example [15]) to the inductances appearing in the circuit. With regard to the trapezoidal rule it is important to point out that it is one of the most accurate discretization rule. Consider

F.N. Koumboulis et al. / Chaos, Solitons and Fractals 11 (2000) 193±206

201

the case of an one dimension inductance u ˆ D…Li† where D ˆ d=d~t, ~t is the original time variable, u is the voltage and i is the current. According to the trapezoidal rule, the discretization takes on the form u…~t† ‡ u…~t ÿ T~0 † ˆ

2L ‰i…~t† ÿ i…~t ÿ T~0 †Š; T~0

…4:2a†

where T~0 is the sampling period. The operator mapping the discrete i to the discrete u is denoted by rD…T~0 †, i.e. u…t† ˆ rD…T †i;

rˆ

2L : T~0

…4:2b†

Since the inductanceÕs appearing in Fig. 1, are strongly nonlinearly dependent on various other quantities in the circuit the usual de®nition u ˆ D…Li† is not acceptable from the passivity point of view. Where D denotes a di€erential operator de®ned upon one of the many variables t10 ; t20 ; t30 ; t40 . An appropriate de®nition guaranteeing passivity is p p u ˆ LD… Li†: …4:3† Applying to the above generalized inductance the generalized trapezoidal rule (referring to many variables) the following discretized relation is derived     p  i u u 2 hp  0 p …t0 † ‡ p …t0 ÿ T 0 † ˆ Li …t † ÿ Li …t0 ÿ T 0 † ; …4:4† T0 L L T

where t0 ˆ ‰t10 ; t20 ; t30 ; t40 Š is the generalized time (involving space and time in spatial dimensions) and T 0 is a vector of shift of the sampling having the form  T …4:5† T 0 ˆ T10 T20 T30 T40 ˆ aT0 ; T0 > 0; a ˆ ‰ a1 a2 a3 a4 ŠT ; ai P 0: The above sampling vector can be interpreted as follows: T0 is the basic sampling period (scalar and constant). T 0 ˆ ‰T10 T20 T30 T40 ŠT is the vector of sampling period for each transformed space and time coordinate. The coecients ai are the weight coecients for every sampling period, respectively. It is clear that as ai ! 0 the discretization becomes perfectly accurate. It is important to note that, in many cases, since all T elements of t0 ˆ ‰t10 ; t20 ; t30 ; t40 Š are in spatial dimensions the weighting coecients can be considered to be 0 0 0 0 equal ai ˆ a ; Ti ˆ a T0 ˆ T0 . Applying Eq. (4.4) into Eq. (4.5) the relation between the multidimensional forward and the backward wave takes on the form b…t0 † ˆ ÿa…t0 ÿ T 0 †:

…4:6†

Using Eqs. (4.2b) and (4.5) the reference circuit yielding a WDF arrangement corresponding to the physical system (1.1) (see Fig. 1) is shown in Fig. 6 with ri ˆ

2Li T00

 T10 ˆ T00 T50 ˆ

…i ˆ 1; 2; 3; 4†;

0

1 0 T 2 0

0 T00

0

T

T00

 T100 ˆ T00

0

 T500 ˆ T00

0 T00

0

0

T

2L0i T00

…i ˆ 1; 2; 3; 4†

 T20 ˆ 0 T00

;

T00

T00

r0i ˆ

;

T

;

0

0

T

 T30 ˆ 0

;

…4:7†

0

T00

0

T

 T40 ˆ 0

;

0

0 T00

T

T

;

…4:8†  T200 ˆ 0

 T600 ˆ T00

T00

T00

T00

0

0

T

0

;

T

:

 T300 ˆ 0 0

T00

T00

T

;

 T400 ˆ 0

T00

T00

0

T

;

…4:9†

202

F.N. Koumboulis et al. / Chaos, Solitons and Fractals 11 (2000) 193±206

Fig. 5. MD passive Kirchho€ circuit for the description of x12 . Fig. 6. Reference circuit.

The operators D…Ti0
† and D…Ti00
† in Fig. 6 are analogous to that in Eq. (4.2b) with the di€erence they are related to the multidimensional shifts of the sampling Ti0
and Ti00
de®ned in relations (4.8) and (4.9), respectively. The circuit in Fig. 6, can serve as a reference circuit and it yields the WDF arrangement which satis®es all the necessary requirements. The corresponding WDF arrangements of Figs. 2±5 are shown in Figs. 7±10 respectively. The respective reference circuit of Figs. 2±5 can easily be realized on the bases of Fig. 6. This is why they are not presented here. The overall WDF arrangement (involving sub circuits of the form of Figs. 6±10) is shown in Fig. 11. The adaptors in Figs. 7±9 are appropriately chosen in order to ®t [15] to the continuous time circuits of Figs. 2±4. An adaptor is analytically presented in Fig. 11.

Fig. 7. Realization of circuits of Fig. 2 in the WDF domain. Fig. 8. Realization of circuit of Fig. 3 in the WDF domain.

F.N. Koumboulis et al. / Chaos, Solitons and Fractals 11 (2000) 193±206

203

Fig. 9. Realization of circuits of Fig. 4 in the WDF domain. Fig. 10. Realization of circuits of Fig. 5 in the WDF domain.

The arrangement of Fig. 11 contains four circuits …0; 1; 2; 3†. Each circuit contains one series adaptor corresponding to the main loops in Fig. 6. In Fig. 11 there are also additional adaptors. Some of them are denoted by N 0 …i† and N 00 …i† …i ˆ 1; 2; 3†. The required multiplier coecients in these adaptors can be

Fig. 11. Analytic WDF arrangement.

204

F.N. Koumboulis et al. / Chaos, Solitons and Fractals 11 (2000) 193±206

determined if the corresponding port resistance are determined. According to Fig. 6, and the results presented in [15] the port resistances can be computed to be 2 4 Ri;k ˆ 0 Lk ; Ri;6 ˆ 0 ; k ˆ 1; 2; 3; 4; 5; i ˆ 1; 2 T0 T0 R0;k ˆ

Ri;0 ˆ

2 L0;k ; T00 6 X

Ri;k ;

R0;4‡i ˆ

4 ; T00

k ˆ 1; 2; 3;

i ˆ 1; 2; 3

i ˆ 1; 2; 3:

kˆ1

Using simple algebraic manipulations all other port resistances of the WDF can easily be computed. According to Eq. (4.8) and Fig. 11 the forward waves are related to the backward waves at previous time instants. Some indicative relations, governing the adaptors, are the following ai;1 …n1 ; n2 ; n3 ; n4 † ˆ ÿbi;1 …n1 ÿ 1; n2 ; n3 ; n4 †;

i ˆ 0; 1; 2; 3

ai;2 …n1 ; n2 ; n3 ; n4 † ˆ ÿbi;2 …n1 ; n2 ÿ 1; n3 ; n4 †;

i ˆ 0; 1; 2; 3

ai;3 …n1 ; n2 ; n3 ; n4 † ˆ ÿbi;3 …n1 ; n2 ; n3 ÿ 1; n4 †;

i ˆ 0; 1; 2; 3

ai;4 …n1 ; n2 ; n3 ; n4 † ˆ ÿbi;4 …n1 ; n2 ; n3 ; n4 ÿ 1†;

i ˆ 0; 1; 2; 3

It is important to mention that the forward and backward waves T

a…t0 † ˆ a…‰t10 ; t20 ; t30 ; t40 Š †;

T

b…t0 † ˆ b…‰t10 ; t20 ; t30 ; t40 Š †

are expressed in the discrete domain as T

a…tn0 † ˆ a…‰…t10 †n1 ; …t20 †n2 ; …t30 †n3 ; …t40 †n4 Š †;

T

b…t0 † ˆ b…‰…t10 †n1 ; …t20 †n2 ; … t30 †n3 ; …t40 †n4 Š †;

T

where n ˆ ‰n1 ; n2 ; n3 ; n4 Š ; tni ˆ ni T00 …i ˆ 1; 2; 3; 4†, or equivalently as a…n† ˆ a…n1 ; n2 ; n3 ; n4 †;

b…n† ˆ b…n1 ; n2 ; n3 ; n4 †:

The equations de®ning N 0 …i† and N 00 …i† …i ˆ 1; 2; 3† are [13] a0;4‡i ˆ a00;4‡i ‡ a0i;6 ;

ai;6 ˆ ÿa00;4‡i ‡ a0i;6

…4:10†

with 00 †; a0i;6 …tn † ˆ b0i;6 …tn ÿ T3‡i

a00;4‡i …tn † ˆ b00;4‡i …tn ÿ Ti00 †;

…4:11†

where 1 b0i;6 ˆ ÿ …b0;4‡i ‡ bi;6 †; 2

1 b00;4‡i ˆ …b0;4‡i ÿ bi;6 †: 2

The port resistances of the adaptors in Fig. 11 at a sampling point are unknown but they can be computed since the incident waves at the adaptors are known. Expressing for each port, the port voltage in terms of the corresponding incident wave, the corresponding port resistance, and the loop current, and writing that the sum of the loop voltages is zero, the following equations are derived 2

12 X kˆ0

2

ak;k

12 X p Rk;k ÿ q^k Rk;k ˆ 0;

7 7 X X p a0;k R0;k ÿ p^ R0;k ˆ 0: kˆ1

k ˆ 1; 2; 3

…4:12a†

kˆ0

kˆ1

…4:12b†

F.N. Koumboulis et al. / Chaos, Solitons and Fractals 11 (2000) 193±206

205

Fig. 12. Couette ¯ow between parallel plates. Fig. 13. Solution of the x-axis velocity u with respect to y (dotted line is for the numerical solution and solid line is for the exact solution) both lines have no visual di€erences.

The terms ak;k for k P 7 correspond to Figs. 7±10. The above nonlinear equations can easily numerically be solved using one of the standard numerical methods (Interval Halving, Fixed Point Iteration, NewtonÕs method, Secant method, etc.). Note these equations lead to a unique solution after expressing the boundary and initial conditions of (1.1) in terms of WDF, i.e. in terms of ak;k and q^k . The above mentioned results and observations are summarized in the following theorem. Theorem 4.1. The approximate discretization of the Navier±Stokes equations (2.1) is the set of equations in (4.12). The numerical solution of the Navier±Stokes equations is produced by solving the equations in (2.12) with appropriate initial and boundary conditions.

5. Numerical example In order to present a ®rst check of the numerical discretization algorithm presented in Sections 2±4, the rather simpli®ed case of Couette ¯ow between parallel plates [14] is presented. In Fig. 12 the two in®nite dimension plates are presented. For the above ¯uid ¯ow the Navier±Stokes equation take on the form ou ˆ 0; ox

0ˆÿ

op^ o2 u ‡l 2; ox oy

where u is the x-axis velocity of the ¯uid, p^ is the total hydrostatic pressure (including the gravitation). The pressure is assumed to vary linearly along the x-axis, i.e. to hold p^ ˆ ÿBx ‡ c; where B ˆ 0:1 and c ˆ 1 the viscosity l is considered to be 0.002. Applying the numerical discretization algorithm presented in Sections 2±4, the x-axis velocity u with respect to y is given in Fig. 13. 6. Conclusions In this paper the problem of numerically integrating the Navier±Stokes equations has been studied. A numerically stable and numerically robust discretization has been presented, via appropriate passive circuit realizations. The present results appear to contribute to the numerical solution of many ¯uid mechanics problems. Based on the present work results regarding analysis of turbulence is under investigation. It is further noted that the present discretization results appear to be necessary step to control distributed parameter systems, namely systems described by PDEs, and particularly for control applications with ¯uid dynamic characteristics.

206

F.N. Koumboulis et al. / Chaos, Solitons and Fractals 11 (2000) 193±206

References [1] R.D. Richtmyer, K.W. Morton, Di€erence method for initial-value problems, Interscience, New York, 1967. [2] R.C. LeBail, Use of fast Fourier Transform for solving partial di€erential equations, J. Comput. Phys. 9 (1972) 440±465. [3] K. Gustafson, K. Halasi, Cavity ¯ow dynamics at higher Reynolds number and higher aspect ratios, J. Comput. Phys. 70 (1979) 271±283. [4] S.F. McCormick, Multilevel adaptive methods for partial di€erential equations, in: SIAM Studies in Applied Mathematics, 1989. [5] S.F. McCormick, Multilevel projection methods for partial di€erential equations, in: SIAM Studies in Applied Mathematics, 1992. [6] R.E. Bank, PLTMG: a software package for solving elliptic partial di€erential equations, in: SIAM Studies in Applied Mathematics, 1994. [7] G.R. Sell, Global attractors for the three dimensional Navier±Stokes equations, J. Dynamics and Di€erential Equations 8 (1996) 1±31. [8] M. Prestin, L. Shtilman, A parallel Navier±Stokes solver: the Meikp implementation, J. Supercomputing 9 (1997) 347±364. [9] A. Fettweis, Multidimensional wave digital ®lters for discrete-time modeling of MaxwellÕs equations, Int. J. Num. Modeling 5 (1992) 1±29.  46 (1992) 209±218. [10] A. Fettweis, Discrete modeling of losseless ¯uid dynamic systems, AEU [11] A. Fettweis, G. Nitsche, Numerical integration of partial di€erential equations using principles of multidimensional wave digital ®lters, Journal of VLSI Signal Proc. 3 (1991) 7±24. [12] A. Fettweis, Wave digital ®lters: theory and practice, Proc. IEEE 74 (1986) 270±327. [13] A. Fettweis, Discrete passive modeling of viscous ¯uids, Proceedings of the IEEE International Symposium on Circuits and Systems, San Diego, CA, WSA, May 1992, pp. 1640±1643. [14] F.M. White, Viscous Fluid Flow, McGraw-Hill, New York, 1974. [15] A. Fettweis, G. Nitsche, Transformation approach to numerically integrating of PDEs by mean of WDF principles, Multidimensional Systems and Signal Proc. 2 (1991) 127±159. [16] K. Meerk oter, Incremental passivity of wave digital ®lters, Proceedings of the European Signal Processing Conference, Lausanne, Switzerland, 1980, pp. 27±31. [17] A. Fettweis, On assessing robustness of recursive digital ®lters, European Transactions on Telecommunications 1 (1990) 103±109. [18] F.N. Koumboulis, M.G. Skarpetis, B.G. Mertzios, On the derivation of the nonlinear discrete equations numerically intergrating the euler PDEs, Discrete Dynamics in Nature and Society (DDNS), 2 (1998) 41±51.