NASA Technical Memorandum 4714

Artificial Boundary Conditions for

Computation of Oscillating External Flows

S. V. Tsynkov

August 1996

NASA Technical Memorandum 4714

Artificial Boundary Conditions for

Computation of Oscillating External Flows

S. V. Tsynkov

Langley Research Center * Hampton, Virginia

National Aeronautics and Space Administration

Langley Research Center • Hampton, Virginia 23681-0001

August 1996
Acknowledgments

I am very thankful to V. S. Ryaben'kii of the Keldysh Institute for

Applied Mathematics, Russian Academy of Sciences, whose ideas

constitute the foundation of our joint research in numerical methods

for external problems and who influenced much of the content of

this paper. S. Abarbanel of the School of Mathematical Sciences,

Tel-Aviv University, is acknowledged for his valuable help in formulating the problem and for numerous useful discussions. I would also

like to thank V. Matsaev and E. Shustin of the School of Mathematical Sciences, Tel-Aviv University, and D. Gottlieb of the Division of

Applied Mathematics, Brown University, for their helpful comments. The suggestions of H. Atkins, T. Roberts, and C. Swanson of

the NASA Langley Research Center have contributed much to the

improvement of this paper.

Available electronically at the following URL address: http://techreports.larc.nasa.gov/ltrs/Itrs.html

Printed copies available from the following:

NASA Center for AeroSpace Information

800 Eikridge Landing Road

Linthicum Heights, MD 21090-2934

(301) 621-0390

National Technical Information Service (NTIS)

5285 Port Royal Road

Springfield, VA 22161-2171

(703) 487-4650
Abstract

In this paper, we propose a new technique for the numerical treatment of external flow problems with oscillatory behavior of the solution in time. Specifically, we consider the case of unbounded

compressible viscous plane flow past a finite body (airfoil). Oscillations of the flow in time may be caused by the time-periodic injection of

fluid into the boundary layer, which in accordance with experimental

data, may essentially increase the performance of the airfoil. To conduct the actual computations, we have to somehow restrict the original

unbounded domain, that is, to introduce an artificial (external) boundary and to further consider only a finite computational domain. Consequently, we will need to formulate some artificial boundary conditions

(ABC's) at the introduced external boundary. The ABC's we are aiming to obtain must meet a fundamental requirement. One should be

able to uniquely complement the solution calculated inside the finite

computational domain to its infinite exterior so that the original problem is solved within the desired accuracy. Our construction of such

ABC's for oscillating flows is based on an essential assumption: the

Navier-Stokes equations can be linearized in the far field against the

free-stream background. To actually compute the ABC's, we represent

the far-field solution as a Fourier series in time and then apply the Difference Potentials Method (DPM) of V. S. Ryaben 'kii. This paper contains a general theoretical description of the algorithm for setting the

DPM-based ABC's for time-periodic external flows. Based on our

experience in implementing analogous ABC's for steady-state problems (a simpler case), we expect that these boundary conditions will

become an effective tool for constructing robust numerical methods to

calculate oscillatory flows.

1. Introduction

The numerical study of problems originally formulated on unbounded domains requires the implementation of special techniques for the "treatment of infinity" (which is necessitated by the restricted

facilities of modern computers). One of the corresponding techniques is based on an artificial truncation

of the original infinite domain, which implies that one must set special boundary conditions at the external (artificial) boundary of the newly formed finite computational domain. The aim of this paper is to

describe the theoretical foundations for constructing such artificial boundary conditions (ABC's) for the

computation of certain unsteady external flows.

Before proceeding to the actual description of the problem, let us first define the concept of exact

ABC's. Namely, exact ABC's are the boundary conditions that enable one to uniquely complement the

solution of the "truncated problem" to the unbounded exterior of the computational domain so that the

original problem is solved. The exact ABC's usually appear to be nonlocal for steady-state problems in

space and for time-dependent problems in both space and time.

Let us emphasize that our main objective in this paper is to construct special boundary conditions

that would model (and in the ideal case equivalently replace) the exterior part of the problem, i.e., the

part we eliminate by truncation. Many examples of such boundary conditions can be found in comprehensive reviews by Givoli. (See refs. 1 and 2.) This formulation differs from another well-known problem related to setting the boundary conditions for numerical algorithms, namely, to construct such

boundary conditions that would ensure well-posedness of the truncated problem and stability of the
integration process in time. In fact, these two formulations are not completely independent. For example, the issue of well-posedness for certain classes of (local) ABC's was thoroughly investigated by

Gustafsson. (See refs. 3-5.) On the other hand, a group of very delicate questions related to the issue of

long-time stability is studied by Carpenter, Gottlieb, and Abarbanel in reference 6 (for some specific

boundary-value problems). The issue of connections between the (highly accurate nonlocal) boundary

conditions that "model the infinity" and the boundary conditions that ensure the long-time stability will

be an interesting subject for a future investigation.

In this paper, we consider an unbounded compressible viscous flow past a finite body or configuration of bodies (e.g., single-element or multi-element airfoil). The behavior of the flow in time is

assumed to be oscillatory. We must emphasize that while talking about the oscillatory time behavior we

mean that some alternating (time-periodic) influence is exerted on the flow (e.g., see experimental work

by Seifert, et al. in ref. 7) and expect that those frequencies that are connected to this influence will

dominate in the solution. We expect that this circumstance will enable us to construct the ABC's without taking into account any other time-dependent effects. From a mathematical standpoint, this case fills

an intermediate position between the steady-state and true unsteady flows.

The steady-state case is relatively simple compared with time-dependent flows. In reference 8, we

have constructed the ABC's for calculating external viscous compressible steady-state flows. These

boundary conditions were based on the concept of far-field linearization and on the application of the

Difference Potentials Method (DPM) of Ryaben'kii. (See refs. 9 and 10.) The ABC's (ref. 8) differ only

slightly from the exact ABC's (a more rigorous formulation of the latter statement may be found in

ref. 8); therefore, the ABC's (ref. 8) turn out to be global in space. However, practical implementation

of these boundary conditions is fairly easy. (See refs. 11 and 12.) They were used along with the NavierStokes code by Jameson, Schmidt, Turkel, and Swanson (refs. 13-15) for computing different external

flows. Numerical experiments show that the global DPM-based ABC's (ref. 8) provide high accuracy of

computations, as well as fast convergence of the multigrid iteration procedure to a steady state. (See

refs. 11 and 12.) The computational cost of boundary conditions (refs. 8, 11, and 12) is not high in comparison with the total cost of the original procedure. (See refs. 13-15.) Generally, the numerical algorithm we used for integrating the Navier-Stokes equations became more robust (in comparison with the

standard procedure (refs. 13-15)) if supplemented by the DPM-based ABC's. (See ref. 8.)

Additionally, we would like to emphasize that the ABC's (ref. 8) were constructed specially for the

steady-state problem and on the basis of stationary governing equations, independent of any specific

technique for solving the stationary equations inside the computational domain. In practical computations, we use multigrid iterations (refs. 13-15) for calculating the steady-state solutions in references 11

and 12. In doing so, we set the ABC's (ref. 8) on each iteration on the upper time level. Of course, the

boundary data on the intermediate stage of the iteration procedure (i.e., until we achieve the steady

state) are not necessarily consistent with the formal "stationary" treatment of the far field. However,

treating the "time-intermediate" boundary data as if it were already steady has been found effective in

computational practice. (See refs. 11 and 12.) We are going to use a similar idea for the time-periodic

case studied in this paper.

True unsteady flows are much more complicated in terms of both theoretical analysis and practical

calculations. In general, the exact ABC's for unsteady problems will be nonlocal in both space and time.

Therefore, the corresponding computational cost may appear to be rather high. This is also true for the

global DPM-based boundary conditions which can be constructed as close to the exact ones as desired.

The corresponding general theory for unsteady problems is contained in work by Ryaben'kii. (See

ref. 16.)

However, an intermediate case of oscillatory time behavior must be less expensive in terms of

required computer resources since the global character of the ABC's in time will obviously be restricted

by the value of one period. Moreover, the theoretical analysis of this case based on the usage of

2
the Fourierrepresentationi time alsoappearsto belesscomplicatedthanthe generalonefrom

reference16sinceinouranalysisweactuallyreducethetime-dependent problemto afamilyof steadystateproblems.

Ontheotherhand,don'tassumethattheoscillatingflowisaparticularand,therefore,anunimportantcase.Forexample,experiments(ref.7)showthatthetime-periodicinjectionof fluid intotheturbulentboundarylayermayincreaseits resistanceto adversepressuregradientswithoutseparation.This

impliesanessentialimprovement of airfoilperformance, upto 60percentforhigh(poststall)anglesof

attack,accordingto reference7. Thephenomenonwasobservedon differentgeometries(original

NACA0015airfoil,thesameairfoilwiththedeflectedflap,andsomeothers),whichleadsusto believe

thatit maybeeffectivelyusedinaircraftdesign.Therefore,anaccuratenumericalinvestigationof the

phenomenonbecomesanimportantissue,andanaccurateprocedurefor settingtheABC'smustbeone

oftheprincipleelementsof anycomputationalgorithmusedfor suchaninvestigation.

Thepreviousexampleisprobablynotauniqueonewherethetime-periodictreatmentof flowinthe

farfieldmightberelevant.In general,fortheoscillatorycaseweproposethefollowingconstructionof

ABC's.First,linearizethegoverningequationsinthefarfield.Then,assumingthatthetimeperiodis

initiallyprescribed,applytheFouriertransformin timeandobtaina familyof steady-stateproblems

(wheretheunknownsareamplitudes).Thelatterproblemsarethentreatedbymeansof theDPM.(See

refs.9and10.)Thecentralideaof theDPM-basedapproachisto equivalentlyreplacetheproblemformulatedona domainbya certainoperatorequationformulatedonits boundary.Foreachoneof the

foregoingsteady-stateproblems(note,thefamilyof theseproblemsisparameterizedbythefrequency,

i.e.,bythedualFouriervariable),thisreplacement resultsin anoperatorequationformulatedat the

artificialboundaryof thecomputational domain(i.e.,connectingtheboundaryvaluesof thesolution).

Theoperatorinvolved(aprojection)is somewhatanalogousto theboundarypseudodifferential operatorsintroducedbyCalderon.(Seeref.17.)Becauseoftheequivalencetotheexteriorlinearproblem,the

previously-mentionedoperatorequation(moreprecisely,theentirefamilyof theseequations)canbe

consideredadesirablexactABC(limitedonlybytheaccuracyof far-fieldlinearization)fortheproblem solvedinsidethecomputational domain.In otherwords,thisoperatorequationadequatelytakes

intoaccounthestructureof thesolutionfromoutsidethecomputational domain,whichmightalsobe

calledthe exact transfer of boundary conditions from infinity. (See ref. 16.)

We actually develop the DPM-based ABC's for the already discrete formulation of the problem. In

doing so, the set of the frequencies involved is obviously finite. Therefore, we can actually compute the

corresponding boundary operator for each one of the steady-state problems arising after the Fourier

transform in time. Then, for reasons of numerical convenience, we represent the solution to the linearized exterior problem in the form of generalized potential. (See refs. 9 and 10.) The density of generalized potential serves as an unknown function in the previously-mentioned operator equation. By using

the generalized potential to set the ABC's we gain more generality from a geometric standpoint. Moreover, we can easily match the solutions of the interior nonlinear problem and the exterior linear problem

when conducting practical computations. (We need to actually calculate the generalized potential only

in some neighborhood of the computational domain, to be discussed later.) Finally, applying the inverse

Fourier transform, we obtain ABC's in a matrix form, which enables easy practical implementation. In

fact, the entire procedure may be thought of as solving the linearized problem outside the computational

domain and then using the obtained solution to close the "truncated system" that is solved inside the

computational domain.

To conclude this introduction, let us point out an analogy to the previously investigated steady-state

case. (See ref. 8.) Namely, we are looking here for a solution to the unsteady problem that is defined on

an initially prescribed time interval (one period) and that meets the periodicity condition in time (at least

in the far field). To develop the ABC's for this case, we solve a certain linear problem in the far field (by

means of the DPM). The latter problem is also formulated for the time interval of one period. The

ABC's for the time-periodic case are basically constructed independent of any specific technique for
integratingtheNavier-Stokesequationsinsidethecomputational domain (as the ABC's (ref. 8) were

constructed irrespective of any specific way for actual computation of the steady state). Based on the

assumption of periodicity in time, these boundary conditions simply close the system that is solved

inside the computational domain; the closure is constructed for the time interval of one period. In practice, however, achieving the true oscillatory regime may require long-time computational runs that

cover many periods. During this long-time integration, each moment we need to update the external

boundary data using the ABC's (i.e., each time step) we treat the flow as it were already time-periodic

(in some generalized sense, see section 3). In so doing, the boundary conditions should guarantee only

the desirable far-field behavior of the solution. This behavior is actually determined by the condition

that all perturbations vanish at infinity (as in refs. 8, 11, and 12 when we were treating the external

boundary data on each iteration as already steady and demanding that the ABC's ensure the decrease of

the solution to the linearized problem at infinity).

This paper is organized as follows. In section 2, we describe the basic formulations of the problems.

Specifically, in subsection 2.1 we describe a geometric setup typical for the numerical solution of external flow problems, i.e., configurations of the finite computational domain and its infinite exterior. In

this subsection, we also introduce the flow equations (parabolized Navier-Stokes) and linearize them in

the far field against the constant free-stream background. In so doing, we obtain a coupled problem,

which is nonlinear inside the finite computational domain and linear outside it. Then, assuming that the

period of oscillating motion is known, we Fourier transform the exterior linear system with respect to

time and obtain an equivalent family of steady-state systems. These steady-state systems must be solved

as a part of the solution to the aforementioned coupled problem. However, we do not solve them

directly since the corresponding domain is still infinite. Instead, we equivalently replace each of the

exterior linear steady-state systems by the generalized Calderon pseudodifferential equation formulated

at the external boundary of the computational domain. To calculate the pseudodifferential operation

(projection) we need a special auxiliary problem that is first formulated on the entire plane for the linearized thin-layer equations (after the Fourier transform in time) with a certain compactly supported

right-hand side. Solvability of this auxiliary problem (in the sense of tempered distributions) is studied

in subsection 2.2. Then, in subsection 2.3, we show how one can replace the original auxiliary problem

formulated on an unbounded domain (entire plane) by a new problem formulated on some rectangle so

that the solutions of the two problems are in a certain sense close to each other.

Section 3 of this paper is devoted to numerics. In subsection 3.1, we introduce a finite-difference

scheme that approximates the linearized thin-layer equations. Since we discretize the equations not only

in space but also in time, we now get a finite (discrete) series instead of the original infinite Fourier

series which implies that the family of steady-state systems to be solved outside the computational

domain becomes finite as well. In subsection 3.2, we construct a difference analogue to the auxiliary

problem on the rectangle, describe the numerical algorithm for its solution (referring to our previous

work for some details) and briefly address our somewhat non-standard concept of convergence for the

solutions of the difference auxiliary problem. Finally, in subsection 3.3 we show how one uses the

recently formulated difference auxiliary problem and obtains difference analogue to the Calderon

boundary pseudodifferential projection. Using this difference boundary projection and also calculating

the generalized difference potential, we actually compute the nonlocal DPM-based ABC's. The ABC's

are first obtained in the Fourier variables and then, after implementing the inverse transform, in the

physical variables as well. Finally, section 4 contains some conclusions and possible generalizations.

2. Basic Formulations

2.1. Governing Equations and Geometric Setup

Let us start with the parabolized Navier-Stokes equations, which are the same as the thin-layer

equations for two dimensions (see ref. 18 by Anderson, Tannehill, and Pletcher):

4
_p _)pu _)pv ^

+ + :u

_u Ou _u _p l _ Ou

P_-t-+Pu_'-xx+Pv_-yy+ _--xx= R-e _"yy_'-yy

_v _v _v _p 1 4 _ _v

P_+PU_x+PV_y+_y : __3 _yl_

OE Oe ()e (Ou _v) 1 F (Ou'_ 2 4 (Ov) 2+_ _ _£q

(1)

Here, x and y denote the Cartesian coordinates, u and v denote the Cartesian velocity projections, p

denotes the density, p denotes the pressure, e denotes the internal energy, bt denotes the viscosity, and Y

denotes the ratio of specific heats. To derive the last of equations (1), we assume that the gas is perfect

and that the Prandtl number Pr = p.Cp/_: is constant (_: is the heat conduction coefficient). We denote the

free-stream parameters, specifically, uo, Vo, Po, Po, eo, and lao, by the subscript "0." We additionally

assume that v0 = 0 and u0 > 0, which does not imply any loss of generality. The system (1) is written in

dimensionless form. The following scales were used to obtain dimensionless quantities: u0 was used for

velocity; Po for density, p0uo 2 for pressure, Uo2 for internal energy, _q_for viscosity, characteristic size L

(typically, airfoil chord) for all distances, and L/u 0 for time. The factor l/Re that multiplies the viscous

P0u0 L

terms in equations (1) arises from the nondimensionalization. Here, Re - is the Reynolds

_0

number.

Note that in our previous work (refs. 8, 11, and 12) we used the full Navier-Stokes equations to construct the ABC's for steady-state problems. In this paper, we are going to use the thin-layer system

(eqs. (1)). This system appears to apply quite well to the description of certain viscous flows (ref. 18), in

particular, the far-field flows that we are studying hereafter. Moreover, for the thin-layer system

(eqs. (1)) we can justify some results on the solvability of its linearized counterpart on R 2, which is

important for the general justification of our construction of ABC' s. Finally, the usage of equations (1)

instead of the full Navier-Stokes equations may save an appreciable amount of computer resources, as

will be seen from further consideration.

Let us now assume that the actual values of u, v, p, p, e, and tx in the far field only slightly deviate

from the corresponding free-stream parameters. For dimensionless quantities, this means

p = l+p u : l+u v = _' {J. = 1+{_ p = (_'M02)-I +_]

£ = [('y'- 1)')'M2] -l +

(2)

where

<<1 _ <<I v <,<1 I.l <<1 _ ,<<(YM02)-I

_; <<[(I- 1)yM2] -I

Here, M 0 = Uo(Tpo/Po) -112 is the Mach number at infinity, which is always assumed to be less than

unity. By substituting expressions (2) into equations (1) and retaining only the first-order terms with

respect to small perturbations _, _, _, _, _, and I], we obtain the following system of linear partial

differential equations with constant coefficients:
where

C _u D_U+F_U+H _2u : 0

_t + _x Oy _y2

(3a)

U _

F _.

u

v C =

p

P

ilOl O0

O1

O0

i°°il O0 D=

10

0 1 -M_

[:oo 1 0 0

H = -Re 4/3 0

0 7Pr -1

IOOil Ol

10

0 l -M_

-pr-lM_

(3b)

The system (3) is the linearization of equations (1) against the free-stream background. We omit the

.tilde in equations (3) since we are going to deal only with linear equations in perturbations henceforth.

Additionally, we used the equation of state e- 1 p (more precisely, its linearization

T-lp

= p - _----:___p) to eliminate internal energy from equations (3).

We have mentioned that equations (3) will be used for the description of fluid motion in the far

field. Let us now define a general geometric setup for the problem under consideration. The original

Navier-Stokes equations are integrated on a grid (e.g., C-type) generated around the airfoil; this grid

covers the finite computational domain which is denoted D m hereafter. (See fig. 1.) We henceforth

assume that the linearization (eqs. (3)) is valid outside the computational domain Din, i.e., on its complement Dex. (See fig. 1.) This assumption is true for large computational domains, i.e., far enough from

the immersed body. As we approach the airfoil, the possibility of linearization in Dex can always be verified a posteriori by analyzing the corresponding computational results (as was done in refs. 11 and 12

for the steady-state case).

To integrate the Navier-Stokes equations on the grid inside Din, we use some finite-dimensional

approximation of these equations. The actual type of the resulting discrete operator (i.e., finitedifference, finite-element, etc.) is not that important from the standpoint of constructing the ABC's; for

definiteness we assume that the Navier-Stokes equations are integrated by means of a finite-difference

scheme. To begin with, we also suppose that this scheme is fully explicit in time. We may think that we

already know the solution for the time level t/ on the entire grid, in particular, 1 = 0 implies the initial

data. When we advance one time step, i.e., calculate the solution for the level t l+l by means of the

scheme, we cannot obtain this solution for the whole grid since some nodes located near the external

boundary of Din will be missing. The actual location of missing nodes depends on the specific structure

of the scheme stencil. For example, a typical central-difference second-order approximation to the spatial part of the Navier-Stokes operator on a structured grid requires a 3 x 3 stencil. Using such a spatial

approximation combined with an explicit integration procedure in time, we can obtain the solution on

the level t t+l at all nodes, except for those that belong to the outermost coordinate row of the grid (designated 1-"! in fig. 1). To advance the next time step (t/+2) we will have to somehow determine these

missing values of the solution on the level t l+l. This will be done by means of solving the linearized system in Dex (i.e., by representing its solution in the form of the generalized potential for each Fourier
Y/2

-Y/2

Y

1_- .................. _ l_ 1

/j/ ........ ii

Inflow !I', - -

X x

Figure 1. Configuration of domains.

mode). In other words, using the solution to equations (3) in Dex, we close the system of difference

equations inside the computational domain Din. The closure we obtain is actually the desirable ABC's.

Note that in the case of steady-state problems the ABC's (ref. 8) were also used to close the subdefinite system of difference equations inside Din. As previously mentioned, boundary conditions (ref. 8)

were implemented in references 11 and 12 together with the pseudo-time iteration procedure for achieving the steady state. (See refs. 13-15.) From an algorithmic standpoint, this approach is almost the same

as the true integration in time, so the ABC's (ref. 8) were applied on the upper time level for each iteration. However, there is an essential difference between the approach in reference 8 and the technique to

be described in this paper. Namely, the former is intended only for the treatment of steady-state problems and is based on the linearized stationary equations, and the latter will take into account the previous evolution of the solution in time.

Additionally, let us note that in the case of implicit schemes we also need ABC's that will complete

the system of difference equations inside Din. Indeed, while integrating the Navier-Stokes equations by

7
means of an implicit scheme one has to solve a certain discrete system on the upper time level (t/+l),

whereas the data from the lower time level(s) play the role of forcing terms. This system will obviously

be subdefinite unless we specify additional relations that connect the values of unknowns in the grid

nodes located near the external boundary. In particular, for the previously-mentioned example of a

structured grid and central differences on the 3 x 3 spatial stencil, these additional relations (i.e., the

ABC's) should connect the values of the solution at the penultimate (the curve F in fig. 1) and outermost

rows of grid nodes. (See also refs. 8, 11, and 12.) Including the missing relations provided by the ABC's

into the system solved on the upper time level, we close this system and then advance the next time step.

Let us now provide an exact formulation of the problem. First, we select those nodes of the grid

where the solution can no longer be determined by the scheme but must be obtained by means of special

additional relations (i.e., by means of the ABC's). We designate this set of nodes v 1. Second, we select

those nodes of the grid where we need to know the solution in order to obtain it on v 1 with the help of

the ABC's. The latter set is designated v. Both v and v I will depend on the structure of the specific stencil. In particular, for the 3 x 3 stencil on a structured grid, v and v I correspond to the penultimate and

outermost rows of grid nodes, respectively. (Also see refs. 8, 11, and 12.) Without loss of generality, we

assume that the artificial boundary F (see fig. 1) is formed by the penultimate row of nodes v, so that all

nodes v I that form the curve F 1 (see fig. 1) already belong to Dex (i.e., to the "linear zone").

Then, we designate the time period by T. Clearly, we can further consider our problem for the time

interval [0,T] without loss of generality. We will also need the following brief notations:

DeT = Dex X [0, T], D Tm= Din X [0, T] , F T = F x [0, T], and F_ = F 1 x [0, T] . The closure of the

finite-difference system in DiT, which we are looking for and which should be provided by the ABC's,

is actually a set of relations expressing u Fr in terms of some data specified on F T. As previously men-

| ....

tioned, these relations will be based on the solution to the hnearlzed system (3) m DTx. The latter system is supplemented (on DTex) by the periodicity condition in time,

u[t : 0 : ult : T ((x, y) _ Dex ) (4)

and the free-stream condition at infinity,

u ---_0 as x2 + y2 --->oo (t_ [0, T]) (5)

The choice of the data on F T that "drive" the ABC's is closely connected to the concept of clear trace,

delineated in references 9 and 10. The question of the possible proper constructions of clear traces for

equations (3) may require a special thorough investigation in addition to the general analysis from references 9 and 10; such an investigation is not a direct subject of this paper. Therefore, we will not comment on this question in our further discussion, we only point out the actual construction we use.

Namely, let us first represent the vector function u(x,y,t) in the form of a Fourier Series in time for any

space point (x,y),

n=oo 27t

int--

u(x, y, t) = _ fin(x, y)e T (6)

n _ _oo

where

T 2n

lin(x,y) = _. 1 _ U(x,y,t)e -int--r dt (n=0,+l,+2 .... ) (7)

0
Insteadof consideringequations(3a),(4),and(5) D r we henceforth consider Dex the family of "sta- ex'

tionary" systems,

• ^n+D_U n F3u n H_2a n

lC0nCU _x + _y + _)y2 = 0

(n = 0, +1, +2 .... ) (8)

parameterized by the frequency ton = 2nn/T, n = 0, +1, +2 ..... and supplemented at infinity by the

boundary conditions

fn(x,y)--)Oasx2+y2---)oo (n = 0,+1,+2 .... ) (9)

which directly follow from formula (5). The matrices C, D, F, and H in the system (8) are the same as in

formula (3b).

For Of_" ] (specified on F) the data that each frequency ton we consider the pair of functions fi_.,

as

/

"drive" the ABC's; here, _ is the normal to F. (Note that if the interior solution is already computed by

means of the scheme inside DT then ^n and can be easily calculated.)

zn ' UF

Our ultimate goal will be to provide a full classification of those and only those functions I_, _v-__ 1

(defined on F) that generate a solution fin(x, y) to system (8) (with boundary conditions (9) defined on

Dex and such that its trace on F coincides with the "source" function itself, i.e.,

r vr, )

(10)

As will be seen fr°m further c°nsiderati°n' the c°rresp°nding set °f functi°ns /^nvr,_') can be

described as an image of a certain boundary projection operator. In other words, the functions

^n Ova"]will satisfy some boundary with (The of this operator equation projection. equation type

vr, a_ )

was

mentioned in the introduction as the one equivalent to the linearized exterior problem.) Let us designate

the corresponding projection operator by P_ (we actually construct this operator in section 3). Then,

specifying any function I^nUr" 3f_'] b_ )

(from Din), we apply P_ projection

inside and consider the

Pr Ur, 0_ ) = vr, _ ) as the right-hand side in equality (10) for the problem (eqs. (8)-(10)).

After solving the problem (eqs. (8)-(10)) on Dex, we find the trace of its solution on F 1 (i.e., on vl),

which in turn enables us to obtain the missing boundary relations that close the system of difference

equations inside BiT. These relations (i.e., the ABC's) are derived using the inverse Fourier transform.

They can be symbolically written as
V

= Pnxop_oR fin, (t E [0, T]) (11)

where the operator R represents some (smooth) interpolation of the discrete functions along the curve F,

and the operator Penx involves the calculation of the generalized potential to solve the problem

(eqs. (8)-(10)). The specific structure of all operators from formula (11) will be delineated in section 3,

where we actually construct their discrete counterparts.

Let us make a few important remarks. First, to formally close the system solved in D.rIn _ we have to

obtain additional relations between the values of the unknowns on F T and on FT. Such relations would

provide ABC's that are completely independent of any specific numerical procedure employed inside

D T . However, to simplify our task and at the same time only slightly compromise the previouslyIll

mentioned independence, we take into account that we almost always integrate the Navier-Stokes equations step-by-step in time (explicitly or implicitly). Therefore, we do not have to construct such ABC's

that would connect the values of the solution at v and at v I for all time moments t e [0, T]. It suffices to

determine u, v,p, and p at v I only for t = T(i.e., at the upper time level) since for all previous moments

these values have been determined when calculating previous time steps. Moreover, the formulation of

the problem (eqs. (8)-(10)), where the right-hand side from equality (10) belongs to the projection

(

image, _, _ )e Im P_, assumes that these data are a result of operating by P_ on the Fourier transform ur, _ ) of some time-periodic function. However, in conducting the step-by-step integration in

(

time, the actual data _ur, 0_ ) may not be periodic until we achieve a true oscillatory regime. Therefore, as mentioned in the introduction, any time we use the ABC's we implement a certain generalized

treatment of the external flow as being already time-periodic. Namely, instead of the true boundary data

u r, 0_ ) at F T, we use the best approximation of this data by periodic functions in the sense of least

squares. This approach will be delineated in section 3, which is devoted to numerics.

Second, we are unable to directly solve the problem (eqs. (8)-(10)) on Dex since the domain is infinite. Handling of this problem will require the additional truncation. Recall that we have already truncated the original infinite domain and have obtained Din; now we also truncate Dex in order to get a new

linear problem formulated on a finite domain, and therefore, available for solution on the computer.

This issue is addressed in subsection 2.3.

Third, we certainly will not solve the problem (eqs. (8)-(10)) every time we need to obtain a closed

system inside Din (i.e., each time step). Instead, using the linearity of the problem, we will specify some

basis in the space of boundary data and solve the problem (eqs. (8)-(10)) one time for each basis function. This approach will enable us to obtain the ABC's in matrix form, which is very convenient for

practical computing. (Also see refs. 8, 11, and 12.)

Ultimately we will deal only with the finite-difference formulations and, consequently, with the

finite Fourier series (instead of the infinite series (6), see section 3). In so doing, the discretization in

time for the linearized exterior problem in DeT should not necessarily coincide with the one used for the

Navier-Stokes scheme inside D_. A more convenient method may be to use interpolation in time,

which was previously proposed in reference 16.

10
Finally, let us mention that since we need to know the solution on F for the whole period T to

restore the solution on v 1, the first few time steps (until the total time reaches T) will require some special treatment. It might be based on the usage of either a larger grid or some other external boundary

conditions for the initial stage of integration in time.

We now proceed to the actual construction of the operators involved in formula (11 ). This construction will be essentially the same for all wavenumbers n (n is contained as a parameter in the corresponding expressions hereafter).

As was mentioned before, the computation of the ABC's (eq. (11)) consists of two stages. First we

apply the projection P_ to provide the proper boundary data (right-hand side of equality (10) for the

problem (eqs. (8)-(10)). Then we find the solution to the problem (eqs. (8)-(10)) in the form of the generalized potential (operator P_x)' Both of these stages will require the application of the DPM. (See

refs. 9 and 10.) In particular, it appears that the computation of P_ and Pffx requires the solution of the

same auxiliary problem (AP) described in sections 2.2 and 2.3 for the continuous formulation and in

section 3.2 for the difference formulation. This AP is actually the main element of the DPM-based

approach. The Green operator of the AP plays in the theory of generalized potentials approximately the

same role as the Green function (or the fundamental solution) plays in classical potential theory. (See

refs. 9 and 10.) The AP is formulated on the entire plane (x,y) for the inhomogeneous counterpart of

system (8) with a certain compactly supported right-hand side i"_ = (.f_',)'_, _, )'g) (to be specified

later on). Namely, we will need to solve the following system,

i(0 Ca n +D _n F _n H _)2t_n

n 0x + ay + = i'n (12) 3y2

on R 2, suppfn(x, y) c Din, and we will require that the solution be unique in the class of functions vanishing at infinity. In other words, system (12) is supplemented by the following boundary condition,

un(x,y)-->O as x2+y2-->_ (13)

which is the same as boundary conditions (9).

Once we are able to solve the AP (eqs. (12) and (13)), we can construct the boundary operator P_,

properly formulate the problem (eqs. (8)-(10)), and finally obtain its solution in the form of a generalized potential. This is actually a very brief description of our DPM-based approach; it will be delineated

in section 3 for the discrete formulation of the problem. Now we will investigate the solvability of the

AP (eqs. (12) and (13)).

2.2. Solvability of Linearized Problem on Entire Plane

We will look for the solution to the AP (eqs. (12) and (13)) in the space of tempered distributions

G' (see ref. 19 by Hrrmander or ref. 20 by Vladimirov), which is a conjugate space to the space G of all

infinitely smooth functions defined on R" that decrease at infinity with all their derivatives faster than

any power of (x 2 + y2)-1/2. Take the Fourier transform

oo

1

un(_, 1]) = _ _ _ un( X, y)e-i_x-irlYdxdy (14a)

e_ t_o

1

_n(_, 1]) = _ _ _ _n(x ' y)e_i_x_iqYdxdy

--oo --oo

(14b)

11
of bothsidesof system(12)andrepresenttheresultintheformof amatrixequation,

Qu = f (15)

Note that in system (15) and henceforth in this subsection we drop the superscript n to simplify the notations. Then, the symbol Q (eq. (15)) is given by

Q = icoC+i_D+irlF-rl2H

it ill 0 i(to + _)

ri 2

f(to + _) + R'--e 0 it 0

4 1]2

iri 0 0 i(to+_)a 3 Re

0 0 i(t_ + _) _ Trl2 i

RePr MO

1"12

RePrM 2

(16)

We first show that system (15) is solvable in G'. For the time being, we do not need any restrictive

assumptions in regard to f; as previously mentioned, f is compactly supported (suppf c Din), and without loss of generality we may think that f is absolutely integrable on R2 (fe LI(R2)). Then, its Fourier

transform f is bounded and continuous on R2; consequently, if we formally write down the solution to

system (15) as

u = Q-I_ (17)

then the properties of the right-hand side in equality (17) are fully determined by the inverse symbol

Q-1. Indeed, it is well known (ref. 20) that if the right-hand side of equality (17) is locally absolutely

integrable on R 2 then it defines the tempered distribution, i.e., the generalized function from G'- The latter will coincide (in the sense of distributions) with the classical function Q-I (_, rl)_(_ ' 1]) everywhere

on R 2, except for the set of singularities of Q-I(_, rl)f(_, rl) (if any). Since in our case the function

(_, rl) is continuous and bounded on R 2, then it suffices to determine whether the function Q-I (_, 1])

belongs to L_oc(R2 ).

To do this, we have to find all singularities of Q-I(_, rl). Calculating the determinant of Q(_, rl),

we obtain

Q(_,rl ) = [_(_+_)4+__

(tO + _)2(_ 2 + 1.12 ) +

(to+_)21]4(4+_ p_ 1

4 _2/]4 116

3 M2prRe 2 M2prRe 2

+iI!to+_)3rl2(7 (to+_)_2r12(4 1 _ 4]((to+_)rl 6-

Re _,3+P_) M0_R e _,3+_r) (°)- 4- _)/]4(1 + L

(18)

Here _ and r I are the variables and to, M o, Re, Pr, and y are parameters. We emphasize that both variables _ and 1] are supposed to be real (see formulas (14)); however, the coefficients of Q(_, 1]) are, generally speaking, complex. Thus, to find singular points of the symbol Q (eq. (16)), one has to find the

real roots of Q(_, 1]) (eq. (18)), which actually implies to find common real roots of two polynomials,

9_Q(_, rl) and _3Q(_, 1"1).First, the point (_ = -to, r I = 0) is clearly one of such common roots. Then, we

12
notethat 3Q(_,rl) turnsinto zeroon thetwo entirestraightlines,_=-03 andr1= 0. Moreover,

_RQ(_,1])hasnootherrootsontheline_ = -03,exceptfor 1]= 0. Further,substitutingq = 0 intothe

equation_RQ(_,1])= 0(seeformula(18))andassumingthat_¢ -03,wefindthefollowingtworootsof

0)M0 o_M0

_RQ(_,1])= 0 thatbelongto thelineq = 0: _I - i - M 0 and _e - 1 + M o" We also observe that if

03 = 0 (which corresponds to the steady-state flows), then all three roots, (-03, 0), _,0)M° 0 ) , and

1 + M o' 0 , merge into one.

In an attempt to find other real roots (if any) of Q(_, q) (eq. (18)), we divide the equation

3Q(_, 1]) = 0 by (03 + _)1]2/Re. (This is possible since we have already proven that no other zeros exist

on the two lines _ = -03 and 11= 0, except for those already found.) The resulting equation,

(0)+_)2(_+ 7"_ _2_ 4 112

_rj- _02k,_ + _r)- _02(1 + _r)_ -g_, 3 p-Z-fie 1"142 - 0

(19)

is of fourth order, and taking into account that the equation _RQ(_, 1]) = 0 (see formula (1 8)) is of sixth

order, we conclude that the polynomial Q(_, 1]) may have not more than a finite number of isolated real

roots in total (three of which have already been found). We emphasize here that this property (finite

number of isolated real roots) presents an essential difference between the problem under investigation

and classical acoustics problems in which the viscous terms in the governing equations are usually

neglected. Namely, for the acoustics equations (i.e., linearized Euler equations) the singular points of

the symbol are no longer isolated. They usually form a curve on the plane R 2 which may cause noticeable difficulties with justification of the uniqueness of solution. These difficulties are similar to those

that arise in studying the Helmholtz equation, which may be referred to as describing acoustics in the

stationary medium. We do not deal with Helmholtz-like equations in this paper; we only note that contrary to the acoustics case, system (12) is presumably easier from this standpoint since the proof of

uniqueness appears to be elementary. (See proposition 4.)

Since equation (19) is of second order with respect to _ we can resolve it for each 1] and obtain

explicit function(s) _3 = _3(q)" Because we are interested only in real solutions, we have to consider a

few different cases.

First, assume that co ¢ 0. Then, rewrite equation (19) as

_2I_4 _/ ] (4 1 )1 +20)_(_ _'_+0)2( 7 _r) 112( + 1 ) 4. 1]4

+ + + - 1 - - 0 (20)

Pr M20 3 -P--r PrJ \3 _ -fir _7_-

and observe that, if M02 = +_r + , then equation (20) degenerates and therefore has a

unique real solution _(3°) = _°)(rl) for any 11. If M02 > + Fr , then we can easily make

sure that the discriminant

D = 40)2(_+_._r)2_4[_+_r _2(4 ___r)110)2 112 4 1"14 ]

(21)

13
is always positive, which means that equation (19) has two different real solutions, _:(1) ---- _(1)(TI) and

_(2)= _(2)(lq) for any rI. If Mo2 < (4_ + Pr 3 + Pr'Y, then the condition D _>0 (see formula (21))

imposes certain restrictions on rl. Namely, we have

3PrRe2I_ 1 1 1

- 8T h MO2( + _rr)+ _11 11/2-<'11 <(3PrRe2[ 1 1

PrRe2M2k,3+-_rX_+-p--r 3+l_r 12(4+_r

where

(22)

Therefore, in this case the real solutions to equation (19), _b = _1)(1] ) and "_Br(27= _2)(1] ), exist only

for rl within the above range. (See inequality (22)).

Now consider the case to = 0 (which corresponds to the steady-state problem). From equation (19),

we easily derive

_2[_+ P,- T Mo-12(4-3

+ 1)] = 1"14(p_) 1+ +4. r}4

(23)

Equation (23) has real solutions, _d) = _l)(rl) and _(2)_ = _2)(rl), only for

Mg > +Fr 3 +Fr . Otherwise, we conclude that the equation 3Q(_, rl)= 0 for to = 0 has no

other real roots except for ({ = 0, rl = 0) and therefore, the same is true for the equation Q(_, rl) = 0.

In practice, we have calculated explicit symbolic expressions for the functions {(0) {(30)(rl)

_l) = _l)(rl) ' and _2) = _(ff)(ri) using Mathematica. (See ref. 21 by Wolfram.)(These expressions

are not presented here because they are fairly cumbersome.) Then, substituting the functions

_0) = _0)(rl) ' _(31) = _O)(rl) ' and _2) = _(z)(rl) into the second equation, 9_Q(_, 1"!)=0, we obtain

the algebraic equations with respect to only one variable q. Clearly the above equations (which are different for the different solutions, _(0)3 = _0)(.q), _(31) = _(31)(rl), and _(2) = _(2)(I]) ) may have real

root(s) if and only if the original equation Q(_, rl) = 0 has other real zero(s) besides those that have

already been found, (-to, 0), _--_--_--_--_-_ 1o 0) 0, 0 , and ion0) 1 + M 0' 0 . Therefore, we finally have reduced the

question about the real zeros of the equation Q(_, 11)= 0 to the question about the real root(s) of certain

algebraic equations of one variable.

Regrettably, the resulting equations (after the substitution of {_0) = _0)(B), {_1) = _b(rl) ' and

{(if) = _(ff)(n) into 91Q({, q) = 0) appear to be too complicated for obtaining general expressions for

their real root(s). However, we may implement the following semi-numerical approach which provides

fairly convincing results.

First, note that the case to = 0 seems to be the simplest one. This case actually admits rigorous analysis without doing any simplifying assumptions. As previously mentioned, equation (23) has no real

14
solutionsfor M02< + Fr 3 + Fr (which implies that the determinant (eq. (18)) has no real roots);

for M02 = + _r 3 + equation (23) degenerates and any pair (_, I"1)of the kind _ is arbitrary,

r I = 0 is its root. Substituting this root into 9tQ(_, II)= 0 (see eq. (18)), we obtain _4(I/M I - 1) = 0,

which yields _ = 0. Therefore, we did not find any new real zero. For M02> + Fr 3 + _r , equation (23) has two different real solutions for any 1"1;moreover, _(1) = _(1)(11) = -_3_:(2) = __(2)(1]).

Since all powers of _ in 91Q({, 1]) are even, we do not need to separately consider {_l) = {_l)(rl) and

_(2) 3 = _2)(_). Substituting _0,2) 3 = _(l,2)(q ) into 9_Q(_, 11)= 0, we obtain the following eighthorder equation with respect to rl: a118 +/71] 6 + cTI 4 = 0, where the coefficients a, b, and c are obviously

real. The explicit expressions for a, b, and c were obtained by means of Mathematica (see ref. 21); we

do not present them here because they are cumbersome. However, using these expressions we can prove

that a > 0, b > 0, and c > 0. Then, it becomes clear that there are no other real roots except for the one

we have already found, 1"1= 0 (which also yields _ = 0). Indeed, the equation arl4+ brl2+ c = 0 for

a > 0, b > 0, and c > 0 may have only essentially complex roots 1"1.Therefore, we conclude that for

co = 0 the symbol Q (eq. (16)) has only one singular point (_ = 0, rl = 0).

Recall that all equations under study generally depend on five real parameters, to, M0, Re, Pr, and y.

To simplify our task, we fix the values of some of these parameters. Let us set y = 1.4 (two-atom gas)

and Pr = 0.72 (air). This choice of values for the ratio of specific heats and for the Prandtl number,

respectively, is most frequently used since it is closely related to numerous practical problems; we will

not consider any other numerical values for these two parameters. We now investigate another simple

case, toe0, M(] = +_r 3 +Fr .Then, wehave

COF112(l+p )+4. 114IF{7

Substituting this expression into 9_Q(_, rl), we obtain a sixteenth-order polynomial with respect to rI.

This polynomial contains only even degrees, namely, 0, 2, 4, 6, 8, 10, 12, 14, and 16. It is possible to

make sure (we always use Mathematica (ref. 21) to perform cumbersome transformations) that the coefficients of the above polynomial are positive for all to (to _ 0) and for all Re; consequently, the corresponding sixteenth-order equation has no real roots. Therefore, the determinant (eq. (18)) has no other

real zeros in this case as well.

We have finally come to the most complicated case, which so far allows us only approximate investigation. Let Mg < PrA3 + Fr . Then, we have to clarify whether the functions 91Q({_1)(11 ), 1"1)

and/or 91Q({_2)(r I), 11) turn into zero for rl within the range given in inequality (22). Both functions are

actually of a general algebraic type (they contain non-integer powers), which means we have only a

remote possibility of accurately (analytically) showing that they have no real roots, particularly because

these functions depend on many parameters. At least at this point we are unable to construct the corresponding rigorous proof; therefore, we use the following graphical approach.

15
To start,we select some representative discrete set of the parameters involved. The range for the

Mach number is known, so we simply choose a few points within this range. As for the Reynolds

number, the representative values for the graphical tests we are conducting may be chosen to be about a

few thousand. Indeed, we are not studying Stokes' flows that correspond to very low Re. As for typical

laminar solutions for the flows around an airfoil, they apparently cease to exist starting with Reynolds

numbers at around a few thousand. Moreover, for turbulent flows with true molecular Reynolds numbers of a few million, one can successfully model turbulence in the far field by introducing a new effective value of the Reynolds number, which also appears to be around a few thousand. (See ref. 12.)

Finally, recall that the periodicity of flow in time is caused by some external influence, and reference 7

reports that the maximum effect (i.e., response) of such an influence corresponds to nondimensional frequencies of about 1. Therefore, we will not consider frequencies much less than unity or much greater

than unity. The upper bound for the band of frequencies originates from the numerics since we are

going to pass from series (6) to the finite Fourier series while actually solving the problem on the computer. (See section 3.)

We also note that the limits for 11 (see inequality (22)) do not depend on

the sign of to. Moreover, since _1)(1],Mo ,Re,Pr, 7, [tol) = __2)(q, M0 ,Re,Pr,?,_lto[),

_(2)(q, MO ' Re, Pr, _[,It.0l) = -_(1)(1"1, M O,Re, Pr, 7, -Itol) (eq. (20)), and all powers of _ and (to + _) in

9_Q(_, 1]) are even, it suffices to investigate the behavior of only one of the above functions for both

positive and negative values of to. We do this by plotting the corresponding graphs for the following

specific values of the parameters involved: to = +0.5, _1, _10, and +_50; M 0 = 0.4 and 0.7; Re = 1000,

2000, and 5000; and T and Pr are still 1.4 and 0.72, respectively. The graphs drawn with the help of

Mathematica (ref. 21) in different scales show that neither of the above curves intersects the real axis.

(We do not present these plots here because they are not of interest except to show that the corresponding curve has no zeros). Relying on this approximate graphical investigation, we may expect that at least

within some range of the parameters involved the symbol Q (eq. (16)) has no other real singular points,

except for those that have already been found.

We use an analogous graphical approach for the case M 2 > + Fr 3 + Pr . We have no prescribed range for rI in this case. However, it is clear that the asymptotics of the functions

_Q(_I, 2)(1]), 11) for large 11is q8, so it suffices to study the behavior of the above functions only on

some finite interval of 1]. We used Mathematica (ref. 21) to plot the corresponding graphs for the same

values of to, Re, 7, and Pr as mentioned before and for M 0 = 0.8 and 0.9. The graphs (drawn in different

scales for different q-intervals, up to -105 < r1 < 105) show that neither of the curves has real zeros in

this case as well.

Summarizing, we conclude that at least for a certain reasonable range of the parameters involved,

M O,Re, Pr, T, and to, we have justified the following proposition.

Proposition 1: The symbol Q(_, 1"1)(eq. (16)) has only three real singular points on the (_, rl)-plane:

---, 0 . For to = O, these three points merge into one.

(-to, 0), 1_---2--_--_- o, 0 , and 1 + M o

To determine whether the inverse symbol Q-I(_, 1]) belongs to L_oc(R2 ), it suffices to investigate

the behavior (integrability) of this matrix function near the three singularities. This investigation actually means that we have to check integrability of each of the 16 elements of Q-I(_, 1])" These elements

are given by (Q-1)j, i = _i, j/Q, 1 _<i, j < 4, where _i,j are the corresponding cofactors.

16
iI

) ¢

\

\

\

\

\

\

\

\

\

\

\

) ",

\

\

\

\

\

\

\

\

) .,

0 1 k

Figure 2. Powers involved in the denominator QQ (black circles) and Newton's diagram (dashed line) for QQ ;

(_,= -¢o, rl = 0), ¢o;_O.

Let us first concentrate on the singularity (_ = -o, r I = 0) for to _: 0. We replace the above expressions for the elements of inverse symbol by their equivalents, (Q-I)j, i = (_i, jg-.)/(O_-.) (a means

complex conjugate), to make the denominator purely real. Since both the denominator QQ and the

numerator 8i, jQ are the sums of monomials of the type const. (to + _)Lrll_m (here const depends on M 0,

Re, Pr, T, and ¢o, and k, 1, m are nonnegative integers), then it would be sufficient to make sure that any

expression of the sort

[const°( to + _)krlt_[

Q?2 (24)

that originates from (8i, ja)/(OO.), l < i,j < 4, is integrable near (_ = -to, r I = 0). Since to _ 0, then the

factors _m do not contribute to the asymptotic behavior of expression (24) near (_ = -to, r I = 0), which is

an essential difference in comparison with the case to = 0. (See the following discussion.) Therefore, we

may investigate this asymptotic behavior by constructing Newton's diagram (see ref. 22 by Walker)

with respect to only two variables, to + _ and r I. Namely, we show in figure 2 the set of points (k, 1) that

correspond to all monomials const. (to + _)kl]l_m involved in QQ. The Newton diagram (ref. 22) is a

lower part of the convex hull of the above set. The diagram is shown by the dashed line in figure 2.

Those points (k, l) which belong to the Newton diagram determine the asymptotic behavior of Q_2 near

(_ = -to, _ = 0).

More precisely, the asymptotic behavior of QQ near the singularity is determined not only by the

lowest degree monomials (see Newton's diagram in fig. 2) but may also depend on some higher order

terms if the form

17
+ +

M o

16 _4 -8

9 Re4pS2M 4q +

_4 F(4 1 2

Re-'_M_ Ltk3+ pr) -81 (0") + _)2'q4

(25)

(which corresponds just to the previously mentioned lowest degree terms that constitute the Newton

diagram) degenerates under some conditions. However, in this specific case the form A_2C°)(eq. (25)) is

positive definite because + Pr - _ is positive for any Pr. Therefore, after some natural change

of variables (see the following text) the asymptotic behavior of the denominator QQ becomes uniform

with respect to the polar angle, which implies that while investigating the integrability of Q-I(_, q) one

may simply neglect all the higher-order terms (black circles above the dashed line on fig. 2) and consider the expression

Iconst. ( 03+ _)krlt_ml

(26)

instead of equation (24). Furthermore, we may only increase the ratio (eq. (26)) by neglecting the third

term (-(03 + _)2_a) in equation (25). Indeed, it is easy to see that in doing so we only decrease the

denominator but still preserve its positive definiteness. Finally, eliminate the factors _ra, for simplicity.

We have already mentioned that _m do not contribute to asymptotics near (_ = -03, rl = 0), 03 _ O. Therefore, to estimate the integrals, we may replace these factors by appropriate constants, e.g.,

tconst . (03+ _)kTiZlI_1_ax Iconst°(03+ _)krit_m[ <

(_4/M4)(03+_)n+(16/9)(_a/Renpr2M4)rl 8 (_4in/M4)(03+_)4+(16/9) (_min/Re 4 4Pr 2Mo)rl 4 8

where minimum (min) and maximum (max) are found on a sufficiently small neighborhood of

(_ = -03,q = o).

Thus, we have reduced the original question of integrability of (_i, jQ)/(QQ) to checking the integrability of the following function:

Iconst° (co + _)krll I

a(03 + _)4 + bq8

(a > 0; b > 0; k, l are nonnegative integers) (27)

on some neighborhood of (_ =-03, rl = 0). Because of the symmetry, it suffices to integrate function (27) only on one quadrant, for example, 03 + _ > 0 and rl > 0. Moreover, since we are studying local

integrability, we also introduce some upper limits for co + _ and for r I, e.g., 03 + _ < 1 and r I < 1. Let us

now change the variables 4ra(03 + _)2 = 4 and ,4rbTI 4 = _ and then proceed to the following integral:

11

1)/2X(/- 3)/4 co°s,II

oo

(28)

Further, make another change of variables, from Cartesian (4, X) to polar (p, 0) coordinates, and for

simplicity, truncate our rectangular domain, {0 < 4 < 1, 0 < Z < 1 } --->{4 2 + _2 < 1}, which obviously

does not influence the result (integrable or not integrable). Finally, we obtain, instead of integral (28),

1 n/2

const[ p 1 + (k - I)/2 + (l - 3)/4

d p2 _ (c°s0)(k- 1)/2(sinO)(t-3)/adOdp

0 0

(29)

18
o 1 k

Figure 3. Black circles represent powers of the monomials in cofactors for (_ = -¢o, r I = 0), ¢o _ 0. Gray area corresponds to

integrability conditions (30).

From formula (29) one can easily derive the conditions sufficient for the integral to exist. Namely, they

are

k-I 1-3

+ T > e (30a)

k-1

--5-- > e - 1 (30b)

1-3

T > e - 1 (30c)

where e is an arbitrarily small positive number.

We now have to make sure that all conditions (30) are satisfied for all cofactors _i,j, 1 < i,j < 4.

First, we note that since k and 1 are always nonnegative integers, then two conditions (eqs. (30b) and

(30c)) are met automatically. Then, to check the fulfillment of the third condition (eq. (30a)) one has to

accurately calculate all monomials involved in all cofactors _i,j, 1 < i,j < 4 and analyze the powers (k,/)

for (o) + _)krl/. This step was done with the help of Mathematica. (See ref. 21 .) In figure 3, we have collected all the relevant powers (k, l) for all cofactors _i,j, 1 < i,j < 4. We also show in figure 3 the range

of powers (k, l) which satisfies conditions (30) (gray area). Using figure 3, one can easil_y conclude that

all monomials involved satisfy all conditions (30). Therefore, the inverse symbol Q-J(_, rl) is absolutely integrable near the singular point (_ = -co, 1"1= 0) for c0 _ 0.

19
,I

w\

\

\

\

\

\

1

\

\

\

\

\

\

\

\

\

IK

\

\

\

\

\

\

\

\

\

\

\

\

\

1

0 1 k

Figure 4. Powers involved in the denominator Q(2 (black circles) and Newton's diagram (dashed line) for Q_2 ; (_ = O, 1"1= 0),

co_O,

The integrability of Q-I(_, rl) for to = 0 near the singular point (_ = 0, _q = 0) is investigated by the

same method. We only note that since _ and to + _ are now the same, both of them do contribute to the

asymptotic behavior of Q-I(_, rl) near (_ = 0, rl = 0). Therefore, the sets of monomials involved, as well

as the Newton diagram, for QQ will differ noticeably from those relevant to the case to _ 0. Indeed, the

asymptotic behavior of the denominator QQ near (_ = 0, Tl = 0) is now determined by the following

form (compare with expression (25)):

A_ ) = (1 + M41 M22 )_8+ M4pr2Re4+(_g _1._.__2 2 M2j_ 2 "_6T12+ _ _4TI4+M_ 1 I(!+ " _r)2-_21_2_ Pr_ 8

(31)

which corresponds to the Newton diagram presented in figure 4.

As in the case to _ 0, the form A_ ) (eq. (31)) also appears to be positive definite since all five coefficients in expression (31) are positive for all Re, Pr, and M 0 < 1. However, the Newton diagram shown

in figure 4 consists of two straight intervals, whereas the one in figure 2 contains only one interval. This

difference is essential because now each of the aforementioned two intervals (see the two-component

dashed polygonal line in fig. 4) will determine its own domain of integrability for the expressions

A_ )

(32)

2O
on the (k, /)-plane. Here k and I are the powers in the numerator of expression (32). Since the form A_ )

is not simply positive definite, but all powers involved are even, and each coefficient in formula (31) is

positive, we can find the corresponding domain of integrability on the (k, /)-plane independently for

each of the two parts of the Newton diagram. (See fig. 4.) To do this for either part of the diagram,

neglect those terms in the denominator which correspond to another part (in so doing, the denominator

may only decrease). Then, formally divide both the numerator and the denominator by the common factor _4. Using the changes of variables analogous to those previously implemented, we come to the following set of conditions sufficient for the integrability of function (32) near (_ = 0, rI = 0):

k-I 1-7

2 + T >e (33a)

k-1

T > e - 1 (33b)

1-7

> e - 1 (33c)

and

k-5 l-1

---_---+ T > e (34a)

k-5

--_ > e- 1 (34b)

l-1

> e - 1 (34c)

Note that three conditions (eqs. (33)) correspond to the upper part of the Newton diagram and three conditions (eqs. (34)) correspond to its lower part. (See fig. 4.)

In figure 5, we show (with black circles) all powers (k, 1) involved in all cofactors _i,j, 1 <_i,j < 4

for the case to = 0. Gray areas on this figure correspond to the range of those coefficients (k,/) which

satisfy integrability conditions (33) and (34). Note that conditions (33c) and (34b) impose some additional restrictions on l and k for the upper and lower components respectively of the Newton diagram in

figure 5. We did not have such restrictions in the case to _: 0. (See inequalities (30).) One can easily see

from figure 5 that all elements of Q-l, (_)i,j_)/(a__.), 1 <i,j<4, are absolutely integrable near

(_ = 0, q = 0) in the case to = 0 as well.

Finally, we only have to show that Q-I(_, 1]) is absolutely integrable on some neighborhood of each

of the singular points _--_0,0 1 and 1 +Mo '0 for to_e0. If we simply ensure that is

integrable on the same neighborhood, then the integrability of Q-l({, 1.1) follows. To do this, first note

that grad Q(_, rl) _e0 at either of these two points. Indeed, it is quite easy to see from equation (18) that

-- co 01ion0/ - ,0 for all M 0 < 1. Then, refer to reference 23, wherein

_ _ 0 at both _,1_---L-_0, 0 and 1 + M 0

Vainberg proves exactly the same statement we need, namely, the integrability of IQ -1 (_, 1"1)[on some

neighborhood of an isolated real zero of the polynomial Q(_, rl) when grad Q(_, 1]) _ 0 at this point.

21
0 1 k

Figure 5. Black circles represent powers of the monomials in cofactors for (_ = 0, rl = 0), co# 0. Light-gray area corresponds

to integrability conditions (33). Middle-gray area corresponds to integrability conditions (34). Dark-gray area is common to

both conditions (33) and (34).

Thus, we can finally formulate the following proposition.

Proposition 2: The inverse symbol Q-I (_ rl ) (eq. (16)) is absolutely integrable on any finite domain of

R 2, i.e., Q-I(_, 1.1) _ L]oc(R2).

In accordance with reference 20, proposition 2 immediately implies proposition 3.

Proposition 3 (existence): The system (15) is solvable in q" for any compactly supported f _ L I(R2);

in equation (15) is a Fourier transform of f ( eq. (14b ) ).

The solution to the AP (eqs. (12) and (13)) that we are looking for may generally be found by

means of the inverse Fourier transform (again, the superscript n is omitted below),

oo o_

1

ii(x, y) = _-_ _ _ u(_, TI)ei_x+irlYd_d_ (35)

--oo --oo

Using the brief notation, we may rewrite formula (35) as u = (u)V= (Q-If) v . However, in doing so

we still do not know whether the function u(x, y) (eq. (35)) satisfies boundary condition (13). Let us

first prove the following proposition.

^

Proposition 4 (uniqueness): lf the solution u of system (12) satisfies the boundary condition (13), then

it is unique in the class of distributions vanishing at infinity.

22
Indeed, any function u that solves system (12) is actually an inverse Fourier transform of some solution

to system (15), u = (u) v. In turn, any distribution u e q' that solves system (15) (see formula (17))

should coincide with the regular function Q-I (_, 11)f(_, rl ) everywhere on R 2 except at the three singular points of Q(_, 11) (since _(_, rl) has no singular points). Therefore, any other solution to system (15)

may differ from u only by a distribution with the support belonging to the three-point set

(-m, 0), _,1- M o' 0 I + M o' 0 . Such a distribution may be only a finite sum of g-functions and

their derivatives. (See ref. 20.) Therefore, if u = (u)v vanishes at infinity, then any other solution to

system (12) will differ from u by an inverse Fourier transform of a finite sum of k-functions and their

derivatives, and, consequently, it will not vanish at infinity since Fourier transforms of _i-functions and

their derivatives are polynomials. (See ref. 20.) Thus, proposition 4 is justified.

Let us now select a finite ball B where

BcR2, B-ED (-0,), 0), (_--_0, 1 0 I+Mo '0

(E>0)

and construct a partition of unity, 1 = gz + g_, where the functions gB and g_ are infinitely smooth

and bounded on R 2. The function g8 is identically zero outside the ball B + e, therefore, g_ is identically zero inside the ball B- e. Note that such functions always exist. (See, e.g., ref. 20.) Obviously,

= (Q-I_) V^ = 1 ,_v ^

v fi (Q- gzf ) + (Q-lg_). We will separately analyze each term on the right-hand side

of the above sum. First, it is clear that Q-lgB[• LI(R 2) because _ is bounded and Q-l• L/oc(R2).

,_V ^V

Therefore, (Q-lgBf) --_ 0 while _x2+ y2 _ oo. For the second term (Q-lg_i') we cannot yet construct a general proof of its decay at infinity. The difficulties here arise from the fact that

Q-1 • Lloc(R2 ) but Q-1 _ LI(R2), i.e., it is not absolutely integrable near infinity. Therefore, a general

proof may require an appropriate regularization of the corresponding oscillatory integral. However, we

retain this question for a future investigation. For the time being, we can formulate the following two

statements. Each will address the vanishing of the solution at infinity for some particular case (or in a

weaker formulation).

First, assume that f _ L2(R2), which is actually not restrictive for our purposes. Then, f • L2(R 2)

(we may treat the Fourier transform here in the sense of Plancherel, ref. 24). As mentioned before,

Q-1 ¢ LI(R2); however, Q-lg_ can be shown to be bounded on R 2. Therefore, Q-lg_f• L2(R2),

^V

which immediately yields (Q-lgb_) • L2(R2). Thus, in this case the solution u to system (12) is represented as a sum of two terms, fi(_ ) + u (2) , where u (l) _ 0 while _ y2 ___oo (true vanishing in the

sense of boundary conditions (13) and u(2)• L2(R2), which may be treated as a "generalized decay".

23
Wealsonoteherethatthestatementonuniquenessproveninproposition4 alsoappliestothefunctions

fromL2(R 2) since the polynomials obviously do not belong to L2(R2).

Second, if we impose some additional restrictions on f, namely, if we require that f be sufficiently

smooth on R 2 so that fe LI(R2), then we obtain a true decay for the second term as well,

(Q-lg_f) v _ 0 while _ + y2 _ oo. Therefore, for a more particular class of the right-hand sides we

may affirm that the problem (eqs. (12) and (13)) is uniquely solvable in q'. We note that for some other

cases (see ref. 9) such a restriction of the class of admissible right-hand sides does not influence the construction of a DPM-based numerical algorithm. We will not rigorously formulate and prove this statement for the specific case currently under study. However, we expect that this property does take place.

These expectations are based on the numerical experience. (See refs. 8, 11, and 12.)

2.3. Truncation of Linearized Problem

As mentioned in section 2.1, we are not going to directly solve the problem (eqs. (8)-(10)). Instead,

we will implement some additional truncation and further solve only a new finite substitute for the linearized problem. Since the problem (eqs. (8)-(10)) will be solved by means of the DPM, we must construct an equivalent finite substitute for the auxiliary problem (eqs. (12) and (13)). Moreover, the same

finite substitute for the AP (eqs. (12) and (13)) will be used for calculating the operator P_, which provides boundary data for the problem (eqs. (8)-(10)). (See section 2.1 .) In this section, we construct the

finite substitute for the AP introducing some additional assumptions in regard to both the smoothness of

the solution we are looking for as well as the rate of its decrease at infinity. This is done in order to simplify the presentation and to avoid unnecessary complications that are not essential for the purpose of

constructing the numerical algorithm. We hope to provide a more rigorous analysis of the approach

described here in a forthcoming paper.

For reasons of numerical convenience and effectiveness, we will use a different method for calculating the solution of the AP, rather than the one from section 2.2. Using this new solution technique, we

will equivalently reformulate the AP on a new finite domain. Namely, let us again take the Fourier

transform of both sides of system (12); however, now we do so only in one Cartesian direction, y (compare with eqs. (14)),

oo

^ 1

/l(X, 1"1) = _ .[ II(X, y)e-irlYdy

--oo

(36a)

oo

1

f(x, 11) = _ S f(x, y)e-iqYdy

-oo

(36b)

(Again, we drop the subscript n hereafter in this section to simplify the notation. Moreover, we retain

here the same notation, u and f, as in section 2.2; however, the left-hand sides of expressions (14) and

(36) are obviously not the same.) Then, we obtain the following family of systems of ordinary differential equations (ODE's):

du(x, rl)

dx

+ Q(rl)u(x, rl) = f(x, 11) (37)

24
where

Q(1]) = D-I

0 iq 0 i¢o

i_+ _2 0 0 0

Re

0 /co+ 4 7]2 /n o

0 0 im+ TTI2 im _2

RePr M 2 RePrM2 0

and

f(x, 13) = D-I}(x, 1])

(the matrix D is defined in formula (36)). The family (eq. (37)) is parameterized by the continuous variable r I, _oo < 11< oo, and x is an independent variable. Recall that the solution u(x, y) we are going to

calculate should vanish at infinity. (See boundary conditions (13).) Consequently, we will generally

impose the following boundary condition on the solution of system (37):

2

ux, 1] _ 0 as Ixl _ _o (38a)

However, in particular cases (see the following discussion and ref. 8 for details) the condition

(eq. (38a)) may appear too restrictive--namely, the cases when Q(rl) has purely imaginary (or zero)

eigenvalues. Therefore, for some selected values of co and r I we will only require

I ,x,nl<- const as Ix[ --_ oo (38b)

Note that we do not consider solutions that grow polynomially, the latter solutions correspond to the

case when Q(rl) has multiple purely imaginary eigenvalues and does not have a basis composed of

eigenvectors.

Once we are able to find (for every rl) a solution to system (37) that would satisfy boundary condition (38) at infinity, then the solution to the AP (eqs. (12) and (13)) can be restored by means of a onedimensional inverse Fourier transform,

oo

1 ^

u(x, y) - 4r_ I u(x, rl)ei_lYdrl (39)

Let us designate the inverse operator for the one-dimensional problem (eqs. (37) and (38)) by Gx(rl).

That is, the solution u(x, rl) to this problem is given by

u(x, rl) = Gx(1])f(x , 1]) (40)

The operator Gx(rl) is obviously linear. Combining formulas (36), (39), and (40), we obtain the following formula for the solution of the AP (eqs. (12) and (13)):

e_ e_

1

iI(x, y) = _ I Gx(rl) I f(x's)e-i(s-y)qdsdTI

--e,o -o_

(41a)

25
Now, we will show how one can pass from the AP (eqs. (12) and (13)) to the new AP formulated on

the strip {-oo<x<oo} x{-Y-<y- <Y} and periodic in the y direction, with Y being the value of the

period. In doing so, we expect that when the period Y grows, Y --_ oo, the solution to the new periodic

AP will uniformly converge to the solution of the original AP (eqs. (12) and (13)) on any strip

{_oo < x < oo} x {-_ _<y < _} where _ isfixed and always less than YI2. We note that the same approach

was used in reference 8 for the steady-state problems.

Hereafter, we assume that all functions involved are defined on the infinite strip

{-_ < x < _, } x - < y < _ . The width of the strip g is initially supposed to be greater than the

diameter of suppf. (Later we will consider the limit g _ _.) We assume periodicity of the solution to

the new AP in the y direction. Then the solution that vanishes as [xl --_ _ is given by

k = o0 Y/2 .2nk

IIy(x,_ y) _ _ {2_k_l - -- -

= .__, Gx_--_)-_ _ f(x,s)e 'Y (s Y)ds (41b)

k = -_ -Y/2

In formula (41b), we use the Fourier series of a periodic function instead of the Fourier integral used in

formula (41a). Our aim is to estimate [u(x, y) - ur(x, y)l fromabove on a finite (fixed) interval (-#, #),

< Y/2, uniformly with respect to x. Let us introduce a uniform mesh in rl, where hrl = 2rdY is the

mesh size, and designate rlk = kh W k = 0, +1, +_2..... Let us then fix some interval (-A, A); we will

always choose hrl (and consequently Y) so that A =hq(K + 1/2), K being an integer. Then,

lu(x,y)-uy(x,y) oo 00 k = oo Y/2 - Y)ds

: _/ooGx<q)/f(x,s)e-i(s-Y'1]dsdl]- _ hrlGx(1]k) _ f(x,s) e-iTIk(s

k = _oo -Y/2

A k=K

<2-_ ""- "'" + " = 2re

k = -K A

Let us separately estimate each of the two terms (the first one corresponds to the finite interval, and the

second one corresponds to the complementary infinite interval).

l/2)hq 0_ Y/2 s)e _lqk(s Y)ds

Gx(TI) _ f(x,s)e-i(s-Y)rldsd_-hrlGx(_k) _ f(x, -" -

1/2)h n -._ -Y/2

k=K[ (K+

- 2_

k =-r Ll(r_

1/2)hrl i oo Y)ds Gx('q) - f(x, s)e-i(s- y)rldsdrl - hrlGx(rlk ) f f(x, s)e -illk(sl/2)hn ....

+ hT1Gx(rlk) Isl>_Y/2f(x, s)e-iTl*(s-Y)ds ]

26
Clearly, the right-hand side of this inequality is actually the sum of errors of the quadrature formula of

^ oo ^

rectangles for the function u(x,'q)e iqy = Gx("q)f_oof(x , s)e-irl(s-y)ds (see formula (40)) on elementary segments of the kind [(k - l/2)hq, (k + 1/2)hrl], k = -K ..... K. Indeed, for each k, k = -K ..... K, the

third term that corresponds to the integration over Is[ > Y/2 turns into zero for sufficiently large Y's

since i'(x, s) is compactly supported. Therefore, one can obtain the following estimate:

-"27t "[l <c°nst'h A xER max u(x,q)e_rly l

_ (-A, A)

<const.h2A max _21_(X' I])+21y[ +yZlu(x, q)[

11e (-a, A)

= (c l+c2[yr+c3y2)h2A (Cl, C2, C3>0)

Note that if we initially assume that the solution u(x, y) decreases at infinity sufficiently fast, then the

^

differentiability of its Fourier transform u(x, n) (see the right-hand side of the previous inequality) follows directly.

For the second expression, we obtain

_1"12-< _-_ Gx(rl) f(x, s)e -i_(s

JnJ A _oo

+

k >K hrl Gx(rlk) Yf/2 f(x,s)e-in(s-y)ds ]

-Y/2

Y/2

Let us replace the integration limits .[ in the second term on the right-hand side of this inequality by

,,_ -Y/2

, as was done when estimating J'Jl" Then,

Additionally we assume that the solution we are looking for has two absolutely integrable derivatives.

Then its Fourier transform decreases faster than Inl-2 and the previous inequality straightly implies

1 c4

_f'12---- _ (Ca>0)

Combining the two obtained estimates, one easily gets

^ c 4

lu(x, y)-uvCx, y)l <-coh_A+-- A

where c o ae__f max (C1 + c2[Yl + c3y2), Co > 0. Clearly, all constants involved in the foregoing estiyE (-_, _)

mates depend, generally speaking, on the specific nonperiodic function u(x, y) that we approximate by

the periodic functions ur(x, y).

27
Nowlet e be an arbitrary positive number. We will choose sufficiently large Ye (i.e., sufficiently

small hnE) so that the following inequality

c 4

coh _ A + -_ < e (42)

is satisfied for all Y > YE" In other words, we require that for a prescribed e and for any h n < hrl _ inequality (42) has real positive solutions A of the special kind, A = (K + 1/2)hrl (K being an integer). The

latter requirement is always met if, e.g., the distance between the real roots of the quadratic equation

cob _ A 2- £A + c 4 = 0 is greater than hrr This, in turn, yields the inequality e 2- 4c0c4h2 _ cohr 12 6 > 0

for hW This inequality is obviously satisfied for any 0 < hrl < h_, where h_ e R is a unique positive

root of the equation e 2 - 4CoC4h2 - coh2 _6 = 0. Since the fulfillment of inequality (42) is sufficient for

the estimate

[u(x, y) - ur(x, Y)l < e (43)

to be true, then we have shown that for all e > 0 one can always find a sufficiently large period Ye so

that for any Y > Ye the absolute value of the discrepancy between the nonperiodic solution u(x, y)

and its periodic approximation ur(x, y), lu(x, y)- Ur(X, Y)I, does not exceed e for all x and for all

-__<y_<_.

Thus, we have reduced the original AP (eqs. (12) and (13)) to the new AP formulated on the strip

{-_ < x < _} x {-Y _ y _ Y }. In section 3, we show that we will only need to know the solution of the

AP in some neighborhood of suppf, therefore, the approximation of the nonperiodic function u(x, y) by

a periodic one, ur(x, y), only on afinite interval -_ <_y < _ is sufficient for our purposes. Let us now

show how to pass from the domain {-oo < x < _} x {-Y < y < 2}, which is still infinite, to a truly finite

domain for the new AP.

,nstea of t now onsi er rcOn u, omain

D o = (0, X) x (-Y/2, Y/2). (See fig. 1.) This new domain D o should completely contain F and F I.

We will reformulate the new AP so that its solution will be determined only on this finite domain D O

and will coincide there with the corresponding fragment of the solution found on

{-_ <x< _o} x{-Y_<y< 2} before the reformulation. As previously mentioned, we only need to calculate the solution to the AP in some neighborhood of Din. Therefore, we are always able to choose an

appropriate X and Y so that this neighborhood belongs to D O , and consequently we only need to construct special boundary conditions at the lines x = 0 and x = X so that the reformulated new AP being

solvedon D 0 is equivalent to the periodic APon the strip {-_<x<_}x{-Y_<y_<Y 1

described earlier in this section. These boundary conditions at x = 0 and x = X will be set separately for each

28
wavenumber k, k = 0, ±1, +2 ..... (see formula (41b)) involved in the Fourier representation of the function ur(x, y). Namely, for each k, k = 0, +1, +__2..... we require that the corresponding Fourier mode,

2

u(x, rlk) -=u(x, 2nk/Y), meets boundary condition (38) at infinity. To exactly transfer condition (38)

from infinity to the finite boundaries x = 0 and x = X, we use the following consideration. Since system (37) consists of ODE's with constant coefficients and since it is homogeneous outside (0, X) (recall

that suppt'(x, y) c Din, and consequently suppf(x, rl) c (0, X) for all _), then it obviously has four linearly independent eigensolutions (in the region of homogeneity). Depending on the structure of the set

of eigenvalues of the matrix Q(rl), these eigensolutions may either increase or decrease (eq. (38a))

exponentially, or they may oscillate (eq. (38b)) while x --->,,_ and while x --->-_. As previously mentioned, we do not consider the last possible case when Q(rl) has multiple purely imaginary (or zero)

eigenvalues and does not have a basis composed of eigenvectors, which leads to polynomially growing

solutions. Sometimes one can analytically make sure that this case really does not take place. For example, we do this in section 3 in the discrete formulation for some particular values of co and q. In other situations, this question may require some additional numerical investigation as in reference 8. At any rate,

to satisfy boundary condition (38), we must prohibit at x = 0 all solutions that do not decrease to the left

(i.e., as x --->-_ ), and prohibit at x = X all solutions that increase to the right (i.e., as x _ _ ). The reason for this asymmetry was mentioned before: once we have purely imaginary (or zero) eigenvalues of

Q(rl) and, consequently, oscillating or constant-in-space solutions (see formula (38b)), then we cannot

always prohibit at both ends of the interval (0, X) all modes that do not decrease in the corresponding

direction. However, it should not affect the result since the final solution we are looking for (u(x, y))

decreases at infinity. (See subsection 2.2.) Moreover, we have proven in reference 8 that once we have a

selected nondecreasing mode in Fourier representation of the solution, then after the inverse Fourier

transform the entire solution will nevertheless decrease. Therefore, we can take into account selected

nondecreasing modes (if any) by simply admitting them at one of the two boundaries, x = 0 or x = X. (In

case we do not do this, the problem may appear overdetermined.) It seems more natural to admit the

nondecreasing Fourier modes (if any) at the downstream boundary x = X. (See ref. 8.)

Now, we calculate the eigenvalues _'r(qk), r = 1..... 4, for the matrix Q(rlk). Those eigensolutions

that increase to the right correspond to eigenvalues _'r < 0, and those eigensolutions that do not

decrease to the left correspond to eigenvalues 9_, r > 0. Therefore, the following boundary conditions at

x = 0 and x = X may be considered to provide an exact transfer of boundary condition (38) from infinity:

?,

S-(rlk)u(0, rlk) = 0 (k=0,+l,+2 .... ) (44a)

+

S ('qk)U(X,'flk) -- 0 (k=0,+l,+2 .... ) (44b)

Here S-(rlk) and S+(qk) are the special rank-deficient 4 x 4 matrices that depend on Q(rlk), with their

ranks equal to the numbers of eigenvalues _,r(rlk) with nonnegative and negative real parts, respectively.

These matrices are given by

S-(rlk) = H (Q(rl k) - _.r(lqk)I) (45a)

9_Lr(qk) <0

S+(TIk) = I"I (Q(rlk) - _'r(qk )I)

91_.r(rlk) ->0

(45b)

29
HereI isanidentitymatrixandproductsinformulas(45)arecalculatedin accordancewiththemultiplicitiesof theeigenvalues. Analogousconditionswill beusedinsection3whiledealingwiththefinitedifferenceformulationof theAP.

Thus,theformulationof the newfiniteAP is nowcomplete.Namely,we haveto solveequations(12)for thecompactlysupportedright-handsidef (suppfc Din ) on the domain D O (see fig. 1)

with the periodicity boundary conditions in the y direction (Y being the value of the period) and with

boundary conditions (44a) and (45a) at x = 0 and (44b) and (45b) at x = X. In the next section, we proceed to the finite-difference formulation of the problem and describe the numerical algorithm for setting

the global DPM-based ABC's.

3. Numerical Method

3.1. Finite-Difference Scheme

Let us introduce a uniform Cartesian grid in D O x [0, T], with _, hy, and x being the sizes of the

grid in x, y, and t directions, respectively. We designate this grid N O ,

N OT = {(Xm, y j, tl) -- (mh x, jhy - Y/2, lx)lh x, hy, x > 0;

m = O, 1..... M, M = X/hx; j = O, 1..... 2J + 1, 2J + 1 = Y/hy;

1=0,1 ..... 2L + l,2L + l = T/x} (46)

We will construct a second-order finite-difference approximation of the system (eq. (3a)) on the

grid N or (see formula (46)) using the stencil shown in figure 6.

(m,j+l,l+l)/_ .......... -_l (m+l,j+l,l+l)

I I

I I

I I

(mj l+l) t ", t.... l,_,l.lj

I 1

• I I

(m,j-l,l+l)( +l.j- ,_/+1) ........... 7/0 (m+l,j+l,l)

(m,j-I,l)l_" ............ t_ (m+ l,j-l,l )

IX

Figure 6. Stencil.

3O
Namely, we use the first-order differences in the x and t directions and second-order central differences in the y direction, and we center the scheme with respect to the point (m + 1/2,j, 1 + 1/2), which

yields

l t . :ul+ ut+j ut uI )

-Um, j Um+l,j-Um+l,j_+l_DI. m+l,j- m,j m+l,j- m,j

C m,j__ + -% J 2[, h x + hx

• /+1 _UI+I UI+I _ul+l "_

l _U/m+ t,j_ 1 + +

+ _1F(u/m, j+ 1 - U/m,j -1 am+ l,j+ 1 urn, j+ 1 m,j-I m+ 1, j+l -m+ l,j- 1

4_ 2hy + 2hy 2hy 2by J

(u t -2Utm +u t u l _2u t +u t

+ H m,j+l ,j m,j-I m+l,j+l m+l,j m+l,j-1

2 + h 2

hy Y

+

UI+I _2ul+I +UI+I

m,j+l m,j m,j-1

h 2 2

y hy

lal+l _2ul+1 ul+I

m+l,j+l m+l,j + m+l,j-I

+ = 0 (47)

The finite-difference scheme (47) is written for the nodes (m, j, l), m = 0, 1..... M - 1, j = 0, 1..... 2J,

l = 0, 1..... 2L with the assumption that we later impose periodicity boundary conditions in time as well

as in the y direction. Note that the stability of the scheme of type (47) was examined for the model scalar

equation

_u _)u _)u 1 _2u

Ot + -_x + _y Re 2 _y2 (48)

It turns out that the corresponding finite-difference scheme for equation (48) is unconditionally stable in

the yon Neumann sense.

Then, using the periodicity conditions (compare with formula (4))

tim,jO = tiZLm, j+ 1 (m = 0, 1, ..., M; j = 0, 1.... , 2J + 1) (49)

we implement a discrete Fourier transform in time (compare with formulas (7) and (6)),

2L _inlx2rt

^n _ 1 Eu I je T (m = 0, 1, M;j=0, 1, 2J+l;n=-L, L) (50) urn, j 2L + 1 m, •...........

l=0

L inlx2rt

ul m, j" = E Um, ^n je T

n =-L

and instead of system (47) obtain

(m = 0, 1..... M; j = 0, 1..... 2J+ 1; l = 0, 1..... 2L) (51)

2 (Um'j+Um+l'j)+Cn hx re, j+__ F .um, J+2hym, j-1 + tim+ l,j+l-u#l+l,j-I 2hy

c. I;,.,7,,j+j a" +_-H -2urn'J+2 m,j-I

hy ^n 2tin+n- j+^n )

Urn+ l,j+ 1 -- 1, Urh+ 1,j- 1

+ ,_) = 0

Y

(52)

_.. (1 2rt'_, (1 2_'_

Here s = zts,n[,,_n'_T)/'t" , c n = cos(5. y);

j = 0,1 ..... 2J.

n = -L, ...,L; m = 0, 1..... M - 1 ; and

31
The finite-difference system (52) is a discrete analogue of the continuous system (8) on the twodimensional grid N O,

No = {(Xm' Yi) = (mhx' jhy - Y/2)lh x, hy > 0;

re=O, 1..... M,M = X/hx; J=O, 1..... 2J + l,2J + l = Y/hy} (53)

In formula (53) hx and hy are the same as in formula (46).

We also note that once x _ 0, then s n --_ i2rcn/T = iton (see section 2) and cn _ 1.

3.2. Difference Auxiliary Problem

Let us construct another Cartesian grid in D ° ,

M0 = {(Xm+ I/2' Yi)- ((m + l/2)h x, jhy- Y/2)lh x, hy > 0;

m=0,1 ..... M-l,M=X/hx;j=0,1 ..... 2J+l,2J+l = Y/hy} (54)

The grid sizes h x and hy are the same as before. The difference AP is formulated for the inhomogeneous

counterpart of system (52) with a certain compactly supported right-hand side. The unknowns for the

difference AP are defined on the grid N O (see formula (53)), and the right-hand side is defined on the

grid M 0. (See formula (54).) In doing so, we obviously have the second order of approximation. We

will define the specific right-hand side for the AP, f_ + l/2, j, m = 0, 1 ..... M - l, j = 0, 1..... 2J, in

section 3.3. As for now, assuming that this right-hand side is already known (suppi',_ + 1/2, j c Din), we

provide an exact formulation of the difference AP and describe an effective algorithm for its numerical

solution.

In accordance with the results of section 2.3, we impose the periodicity boundary conditions in the

y direction,

^ n u n (m = 0, 1, M) (55a) am, 0 = m, 2J + 1 "",

^n ^n (m 0, 1, M) am,-1 = Um, 2J = .-.,

Then, we implement a discrete Fourier transform (compare with formulas (36)),

2J . . 2n

_n

I1 n je -'kjhy-'_ (m = O, 1, M; k = -J, J)

1

Um'k - 2J+ 1 m, -", "",

j=O

(55b)

(56a)

2J _ikjhr _

^ 1 __ ^n fn+l/2'k - 2J+l fm+l/2'je (m=0,1 ..... M-1;k=-J ..... J)

j=0

and obtain, instead of the inhomogeneous counterpart to system (52),

AkUm+l,k+ n n B_ un, k = tn + 1/2, k (m =0, 1, ......M-l'k=-J .... J)

(56b)

(57)

where the 4 × 4 matrices A t and BE are given by

Sn Cn Cnr k Cnt k

A t = _-C+_D+-_--F+-_-H

(58a)

s n c n Cnr k Cnt k

B_ = _-C-_xD+--_-F+-_--H

(58b)

32
fl 97r'_ _

Here

rk = isin(khy=_)/hy, tk = -4sin2_2kh_yyhy, ,, and C, D, F, and H are defined in formula (3b). For each wavenumber k, k = -J .... J, system (57) is composed of ordinary difference equations, and it is a discrete analogue of system (37). To find a solution to the difference AP, we will have

to solve system (57) for all k, k = -J .... J. However, the formulation of the difference AP is still incomplete. To complete it, we have to set some boundary conditions at m = 0 and m = M (as was done at

x = 0 and x = X for the continuous case in section 2.3). These boundary conditions should guarantee the

desirable far-field behavior of the solution (i.e., decay at infinity). They will be formulated separately

for each wavenumber k, k = -J .... J, i.e., the system (eq. (57)) will be supplemented for each k by some

boundary conditions at m = 0 and m = M. The idea for constructing these boundary conditions in the

discrete case is analogous to the one implemented in constructing boundary conditions (44) and (45) for

continuous system (37). Namely, when formally considered on an infinite one-dimensional mesh,

---oo< m < oo, system (57) obviously becomes homogeneous at least for m > M and m < 0. The homogeneous system has four linearly independent eigensolutions: those that correspond to [a_'(k)f <1

decrease to the right (i.e., as m _ co); those that correspond to [am(k) > 1 decrease to the left (i.e., as

m _-oo); and those that correspond to [lan(k) = 1 have either constant or oscillatory behavior. Here,

[am(k), r = 1..... 4, are the eigenvalues of the matrix Q_ d--ef(A_)-lB_. Let us note that while calculating the eigenvalues [ar(k) for the stationary Navier-Stokes equations (ref. 8) (the eigenvalues are calculated numerically using standard NAG subroutines), we have found that for all specific sets of the

parameters involved (i.e., grid sizes hx and hy and hydrodynamic parameters Mo, Re, Pr, and 7) the absolute values of eigenvalues were never equal to unity except for the case of zero wavenumber, k = 0. For

k = 0, we have obtained a multiple eigenvalue I[a(0)l = 1. (See ref. 8.) However, even in this case the

system matrix still has a basis composed of eigenvectors, which provides us with the reason for not considering the polynomially growing solutions in reference 8. For system (57), we also have a particular

case when the eigenvalues of the system matrix become equal to unity in absolute value. Namely, it is

easy to see from formula (58) that Q0 = (AO)-IBO = -I (identity matrix). Obviously, Qo has four linearly independent eigenvectors; therefore, we do not have polynomially growing solutions in this case

as well. As for other values of k and n, a numerical check (as was done in ref. 8) will always be necessary to determine whether the eigenvalues I_trn(k)l = 1 exist. If such eigenvalues do exist, a check is

also necessary to determine what their multiplicities are and if there is a basis composed of eigenvectors. Relying on our previous experience (ref. 8), we assume that while solving system (57), we can

restrict ourselves by considering only these two cases: la_(k) _ 1 and [am(k) = 1 with the full system

ofeigenvectors. Nontrivial Jordan blocks (of order more than l) for [am(k) = l are excluded from consideration. Note, if the basis composed of eigenvectors does exist for [a_n(k) = 1, then system (57) will

be treated exactly in the same way as in the case Jan(k) ¢ 1 (the only difference is that the stability constant becomes proportional to M).

Returning to the question of setting the boundary conditions for system (57) at m = 0 and m = M, we

require that, analogous to the continuous case (see section 2.3), boundary conditions at m = 0 should

prohibit all modes that do not decrease to the left (i.e., as m -_ -_,,) and boundary conditions at m = M

should prohibit all modes that increase to the right (i.e., as m --> _). Therefore, we may represent the

desirable boundary conditions in the form of matrix relations (compare with formulas (44)),

33
^

sn-(k)u_, k = 0 (k = -J ..... J) (59a)

n+ '_n

S (k)UM, k = 0 (k = -J ..... J) (59b)

where

sn-(k) = H (Q_- _trn(k)I)

l..,"(k)[> 1

(60a)

Sn+(k) = H (60b)

(compare with formulas (45)).

Thus, the final formulation of the difference AP is the following. One should solve the inhomogeneous counterpart to system (52) in D O on the grid N 0 (see formula (53)), where the right-hand side

_n+ 1/2, j (suppfn+l/2, j CDin) is specified on the grid M 0 (see formula (54)), with periodicity

boundary conditions (55) in the y direction and boundary conditions (59a) and (60a) at the line m = 0

and (59b) and (60b) at the line m = M.

To solve the difference AP, we implement the following numerical procedure. First, apply discrete

Fourier transform (56) to both sides of the finite-difference system, then solve the system of ordinary

difference equations (57) with boundary conditions (59) for each wavenumber k, k =-J .... J, and

finally restore the solution by means of the inverse Fourier transform,

k = J ^ ikjhy2_y_

Um, ^n j : Z ^nUm, ke (m = 0,1 ..... M;j = 0,1 , ..., 2J) (61)

k = -J

The type of boundary conditions (59) (which are imposed separately for each wavenumber k) makes the

choice of numerical method most relevant. An effective algorithm for solving one-dimensional problems (eqs. (57) and (59)) is delineated in our work (ref. 25). We do not reproduce the corresponding

results here, we only note that this algorithm may be thought of as a version of the well-known successive substitution technique but without its "inverse" or "resolving" part. The computational cost of the

numerical procedure in reference 25 as applied to solving the problem (eqs. (57) and (59)) is O (M)

operations (for each k, k = -J, ... J).

Let us now briefly describe the concept of convergence for the solutions Of the difference AP.

According to section 2.3, we approximate the nonperiodic solution by a periodic one on a finite interval

-y < y <__ when the period Y grows, Y--> _. In its own turn, an approximate solution to the periodic

problem is found by a finite-difference method on the grid with sizes hx and hy. Therefore, we will consider (uniform) convergence of the periodic difference solution (i.e., solution of the difference AP) to

the nonperiodic continuous solution (i.e., to the solution of the original continuous AP) only on a finite

rectangle (0, X) × (-y, _) (this rectangle should be large enough to contain at least F1) rather than on the

whole domain of the difference AP. Moreover, we will consider this convergence not only when the

grid size vanishes but also when the period Y synchronously increases, i.e., as (h x, hy, Y) ---->(0, O, _).

Of course, the rate of decrease for the grid sizes hx and hy and the rate of increase for the period Y are

not independent; some estimates connecting these rates can be found in reference 8. Furthermore, some

numerical experiments from reference 8 show that the presented construction of the difference AP does

ensure the convergence of its solution to the solution of the continuous AP in the sense just described.

34
3.3. Computation of ABC's

(

In accordance with section 2.1, to set the ABC's we need to know the following data: _[^nuv' "uv|04)"

Here, _ is the normal to F. When integrating the Navier-Stokes equations step-by-step in time, we

/, \

assume that ./^nuv,04 )]is provided from inside Din; then we use these data to restore u nvl,which enables

us to advance the next time step. However, as we carry out our analysis in the Fourier space, we cannot

 onsi rl Uv, uv, )' 04 ) as the actual values obtained inside the computational domain. To get ^n

we first have to calculate the Fourier transform of the function uv, 04 )" Without loss of generality, we

may always think that the latter is specified at the following points:

vx{xlll = 0 ..... 2L+ 1, x(2L + 1) = T} (62)

Of course, actual discretization in time for the Navier-Stokes equations inside DiT should not necessarily coincide with the one used for the solution of the exterior linearized problem. (See formula (46)).

However, we may always use some interpolation in time to obtain the boundary data on a uniform mesh

with respect to t (eq. (62)), which is convenient for further consideration. Hereafter, we simply assume

that this interpolation (which is one-dimensional in time and of sufficiently high order) has already been

implemented for each node v, if necessary.

Another important issue related to the step-by-step integration in time is that the function u v, _4 J'

which provides the boundary data, is not necessarily time-periodic until we achieve a true oscillatory

regime. However, for the purpose of constructing the ABC's, we will propose some generalized treatment of the boundary data as being already periodic. Namely, let us formally calculate the Fourier coefficients of uv, _4 )'

^n _u_'/ 1 2L+I Z (U/v,_.__).- _U/_-in lx_ uv, 04 ) - 2£+1

/=1

(n = -L ..... L) (63)

Then, it is well-known (see ref. 24 by Kolmogorov and Fomin) that the time-periodic function

V lv ' -_ ) Z _l n 0 _1_ in T

n = -L

(1=0, 1..... 2L+I) (64)

minimizes the following functional

35
(11"112is a usual Euclidean norm) on the class of periodic functions; i.e., 1%1 vv, _4 ) from formula (64) is

the best periodic least squares approximation of u v, 24 )" Relying on this property, we will further use

Fourier coefficients (63) as the boundary data that "drive" the ABC's (which may be referred to as the

generalized treatment of the boundary data as being time-periodic). Clearly, as we integrate the NavierStokes equations in time and approach the true oscillatory regime, the "source" function u v, _4 ) and

its Fourier series v,,, 04 ) (eq. (64)) also approach each other.

We now implement the DPM (refs. 9 and 10) to actually calculate the ABC's. We note that the

boundary data are specified on the curve F, which is positioned arbitrarily with respect to LUv, "-_-_j

coordinate lines of the grid N O. (See formula (53).) Moreover, we do not impose any restrictions on the

shape of F itself. In our opinion, the DPM (refs. 9 and 10) provides an ideal tool for treating such geometrically complicated problems.

Let us introduce the following discrete sets. We consider a six-node two-dimensional stencil

Stra+ I/2, j = {(Xra' Yj)' (Xm' Yj+ 1)' (Xm' Yj- !)' (Xm+ 1' Yj)' (Xm+ 1' Yj+ 1)' (Xm+ 1' Yj- 1)}

This stencil is actually a projection of the one from figure 6 onto the plane t = const. Obviously, the discretization (52) was obtained using Strn + 1/2,j" Then, we define

Min = ff_Of'_Din M= PdO\Min N = _.) Stm+l/2, j Nin = _-) Stm+l/2, j It = Nf'_Ni n

(x_+i/2, yj)_ M (xs_+i/2, yi) E 9_j.

Clearly,thesetof gridnodes I'islocatedneartheartificial boundary F. We willcallthissetthegrid

boundary (byanalogy).The setsMin,_ Nin, N, and _,areshown infigure7.

Further,we willneed tointerpolategridfunctionsfrom N O tothepointsv Ic F l.Let us selectall

thosenodes _cc N O thatshouldbe takenintoaccountonce constructinglocalinterpolationformulasof

sufficientlyhigh (e.g.,second)order.All the nodes _carelocatednot farfrom F I.Without lossof

generality, we always may assume that _¢c N- We denote the operation of local interpolation from the

Cartesian grid Nto v 1 by Rvl N .

Let us also introduce the set of collocation points a c F and the space of eight-component vector

functions Wo^ n _ wa^n defined on the set o. The elements of v&_ will be used as unknowns for the boundary operator equation, which will replace the exterior linear difference problem. Henceforth, we will

^n treat wa as vectors containing the values of fin, _n, /3n, and On and the values of the derivatives Ofin

24 ' 04 ' and ,-,_ at the points o; here, _ is the (outward) normal to F. Note that the functions wa are

/

bu_ _

the discrete approximations of l Uv, ^n _) from section 2.1.

Generally, the sizes h x and hy of the grid N Oand the size ha of the one-dimensional collocation grid

on the curve 1-"are not independent; they should be correlated to each other in a certain way. This

requirement is a consequence of convergence conditions for the DPM algorithm. (See ref. 9.) The

theoretical questions concerning the correlation between the sizes of grids N O and o are delineated in

36
I

I

b

I\

r _

p b

i r

i i

J i i

_ r r

; i i

i

I

I i I

r

F

i

i

i

i

J i

i

I

i

I

I I

I I

I

I

I

r

I

I I r

I I I

i I i !

I I I

I I ,

I I I I

I I q r

I I

I i ,

' II P r 'I

I , I

, i I , 1

I

_ r

I

I

I I

I ,

I

I

I

I

I

I

I

I

I

I

I

r

I

I

I

0

Figure 7. Grid sets. Grid 9_ --continuous thin vertical lines; MO--continuous thin horizontal lines & dashed thin vertical

lines; continuous boundary F--thick dark line; Min--gray boxes; _---gray circles; 5_--white boxes; 9_--white circles;

_[ = :7( _ Nin --big white circles & boxes.

reference 9. As concerns practical applications, the final choice of grids is always done taking into

account some previous computational results. In particular, it seems useful to conduct the computations

(see refs. 11 and 12) for the set of collocation points, which is more concentrated at the outflow part of

the external boundary in the wake region and uniformly spaced at the inflow part of the external boundary. Moreover, sometimes the relation Iol -h'l l/2 appears to be proper. At any rate, for each specific

class of problems (determined both by the geometry of computational domain and by the parameters of

fluid at infinity) one always can make an appropriate choice of the grids N O and c relying on general

theory (ref. 9) and on the numerical experience.

Art ^/1

Let us now specify some wo e Wo and implement the following procedure. First, we smoothly

^ ^

interpolate wg along F (i.e., along the smooth components of l") and get the function Rwg. Here, R is

an interpolation operator. Then, we drop the normals from the points T to F and find the values of Rwg

at the foots of these normals. Since w_n (and consequently R6,g ) contains the values of both _, _n _b',

and 15rtand their normal derivatives and since the distance between any node T and the curve F is small

(of order h), we may approximately find _n _n,/Sn, and On at the nodes _, using the two first terms of

the Taylor expansion. We will designate the entire operation of continuation of the boundary data from

37
A n All

to ), as r_, n.¢owa = u_,. Note that the preceding algorithm of continuation generally applies to the

smooth parts ofF (where the normal exists). In practice, however, the curve F is usually not smooth (see

fig. 1), and it is impossible to construct an appropriate normal when the node ), is located in some neighborhood of the breaking point of the curve. The construction of the operator r_ in this case is based on

the existence of two linearly independent directions along the curve, which enables us to obtain the

desirable continuation anyway.

Now, using the calculated continuation of the boundary data, ^II _ II u,t = 7r_,awa , we construct the following grid function:

UNO ^ll m,=

J

(65)

which is defined already on the entire grid N 0. (See formula (53).) Then, we substitute the function ^II UN0

from formula (65) into the left-hand side of system (52). Generally speaking, UN0 ^ n does not satisfy equations (52). Therefore, we generate some nonzero right-hand side, which we designate L0UNO. II II Here, L_

is the linear operator defined by the left-hand side of system (52). This operator takes the functions

defined on the grid N O (eq. (53)) as an input and generates the functions defined on the grid M 0

(eq. (54)) as a result. Finally, we truncate the function LouNo n II to the set M/n, which yields

f_o m+ l/2, j

L.,,rl"l^l'l

= --0U:%0 m + 1/2, j

0

if (Xm+ 1/2' Yj) _ Min

if (Xm+ 1/2' Yj) q_Min

(66)

We will use the function fecP^ll --f_+^n 1/2, j from formula (66) as the right-hand side for the difference AP;

by definition, suppf,_ + 1/2, j c Dil l . Once we solve the difference AP with the right-hand side ^ll fM 0 (eq.

n^ll (66)), we get the function G0fM0. Here, G_ is the Green (i.e., inverse) operator of the difference AP.

ll^n The function GofM0 is defined on the grid N °. As it is considered only on the sub-grid N c N °, it is

called the difference potential with the density u,^n,t P_cz/u_,n ^ II = G0fM0n^ll N" (See ref. 9.) Clearly, the differ-

^n

ence potential satisfies equations (52) since f_0 = 0 on _ moreover, it satisfies the boundary condiII "II is a discrete realization of the generalized tions of the difference AP. The difference potential P_tu_,

potential mentioned in the introduction. Later, we will find an approximate (i.e., difference) solution to

the problem (eqs. (8)-(10)) in the form of a difference potential and then use this solution to construct

the ABC's, i.e., to obtain the missing relations between the unknowns at v c F and at v 1 c F 1. Therefore, we will need to know the difference potential only on the two subsets of N located closely to F

and F 1 on _, _ N and on _cc N, respectively.

II

Indeed, once we calculate the difference potential on ](, we can then construct the operator P"t as

the trace of the potential, II^ll defnn ^n "t "

P,tuv = rN_,u, t This operator will be the key element of the boundary

operator equation of the DPM. Actually, P_ is a difference boundary projection operator (ref. 9), which

substitutes P_ (see section 2.1, eq. (11)) in practical computations.

38
Oncewecalculatethedifferencepotentialon n, we can interpolate it and find the trace of the solution to the linear problem on v l, Uv_ ^n = pv13,u3,n ^ n = RvlNPa@/uy.n ^n Thereby, we obtain the desirable relations between the unknowns at F and F 1 through the solution of the linearized exterior problem. Let us

now recall that we have replaced the original infinite-in-space problem by the periodic problem formulated on the strip {-_ < x < _ } x - < y < _ , claiming that we will need to know the solution only on

some neighborhood of Din, and therefore the convergence on a fixed interval, -_ _<y _<_, is sufficient

for our purposes (see subsection 2.3). Since we do not need to calculate the difference potential anywhere except at 7 and _, the previous statement is now justified.

We now formulate the main result of the DPM theory. (See ref. 9.) Consider the entire space of grid

functions u_, ^ n defined on 7. Those and only those elements of this space which satisfy the equation

I1^_ ^n

P_,uv = uv (67)

can be complemented to Nso that the complement solves system (52) (with boundary conditions (55)

and (59)). As previously mentioned, the operator P_ is a projection; it is a discrete analogue of the Cal-

^/1 the result n ^/1 (as well as deron boundary pseudodifferential operators. (See ref. 17.) For any u_,, Pvu_,

/1 ^/I

Pv_,u'_ ) can be explicitly calculated in accordance with the aforementioned procedure,

^/'1 ^n tlAn An /'/^/1 I1 At/ /1^/1 J1 A/1

u_t _ uNo _ LouNo _ f_vto_ Gofe¢p _ PN_,ul, ==>Pl, ul,, Pv:/U_/•

( In section 2.1, we have declared

our goal as to characterize those and only those (v_, _ ) that

would solve equations (8) with boundary conditions (9) on Dex and coincide with the trace of the solution on F. Instead, we have provided an analogous classification (see eq. (67)) for the discrete rather

than for the continuous formulation of the problem. Therefore, we have equivalently replaced linear

system (52) on N, along with the boundary conditions (55) and (59), by the boundary equation with

projection (eq. (67)). Consequently, we can now specify the proper boundary data (see equality (10)) for

the discrete counterpart of the problem (eqs. (8)-(10)). Namely, let us take any (u_, _u_] 0_ )

(provided

from inside Din ) and interpolate it along F to the set of collocation points o, w_ = Roy uv, 0_ ).

Then, continue A/1 '_ WO using the operator r_o, and finally apply P_, (which requires solving the AP). In

accordance with the previously formulated main result, the grid function

v_,^n= p_,n_,oRov Uv, _ )

(68)

admits the complement to N that solves the problem (eqs. (52), (55), and (59)).

Thus, we have completed the first stage of constructing the ABC' s (section 2.1 ) and now proceed to

the second one. Instead of the problem (eqs. (8)-(10)), we will consider its discrete counterpart: to solve

system (52) on N with external boundary conditions (55) and (59), and with boundary condition

39
^n ^n (69) Uy = V T

^n in equality (69) comes from formula (68). The solvability of the problem (eqs. (52), (55), (59), at 7;, wt

and (69)) is guaranteed by the special type of boundary data provided in formula (68).

To actually find the solution to the pffoblem (eqs. (52), (55), (59), and (69)), we calculate the differ-

. tl ^rg /1

ence potential P_v_, with the density vr from formula (68). Then, we interpolate the potential from N

to v 1, which yields

finVl = RvlNP_(yp._,_R_v nn lAn Uv, /d"An/An / _ ) = T Uv, _)_ )

(70)

Equality (70) provides the missing relations between the unknowns at v and at v I in the Fourier space;

these relations are based on the solution to the linearized exterior problem. Therefore, equality (70) is

almost the desirable ABC, and the only remaining step is the inverse Fourier transform in time. Before

implementing the inverse transform, let us note that the entire algorithm becomes most convenient from

a practical standpoint if we calculate the matrix representation of the operator "i_n from formula (70). To

do that, we choose some basis in W_, ^ n e.g., the simplest one, composed of the vectors like (0, .... 0, 1,0,

.... 0), and implement the entire proceeding procedure. More precisely, we calculate

Uvl^n = RvlNP_(_p.f_._owo,n n ^n for each basis vector w_. In so doing, we obtain the matrix of

^n

RvlNP_6,P._n_, o (each column will be the response to a specific basic function wa ) and then, multiplying the above matrix from the right by the interpolation matrix Row, we finally obtain the matrix representation of ,_n (Note that we do not start from basis functions (^ n bu_ 1 since the number of nodes c

• _Uv'

is usually much less than the number of nodes v.) In fact, it is possible to show (see ref. 9) that for any

wo, ^ n p_Ca_p_,nraw on n ^ n _=p_(_n,towon^ n and therefore we need to calculate the matrix RvlNP_n_, _n . Clearly,

the computation of each column of the matrix Rv_ NP_L_n_to requires one solution of the difference AP,

which, in turn, involves the direct (eq. (56b)) and inverse (eq. (61)) Fourier transforms and the solution

of the problem (eqs. (57) and (59)) for each wavenumber k, k = -J ..... J. Either Fourier transform will

require only O(M • J) rather than O(M • j2) operations. (For definitions of M and J, see formula (46).)

Indeed, the support of the right-hand side ^nfM 0 is actually concentrated near F since UN0 ^n differs from

zero only on Tand the operator L_ is local. Therefore, while calculating direct Fourier transform (56b)

for each m, m = 0, 1 ..... M - 1, only a few values f,_ + l/v, j differ from zero, and consequently the total

cost of this computation is O(M. J) operations. Analogously, while calculating the inverse Fourier

transform (61) for each m, m = 0, 1..... M, we need to know urn, ^_ j only for a few selected values ofj

since all other (xm, yj) do not belong to K. Therefore, the total cost of this computation is O(M •J)

operations as well. Finally, the solution of the problem (eqs. (57) and (59)) for each wavenumber k,

k =-J ..... J costs O(M) operations. (See ref. 24.) Adding all these quantities, we obtain a total of

O(M. J) operations for the computation of each colunm of the matrix RvlNP_D/n_, a. We see that

although the entire algorithm requires repeated solution of the difference AP, the solution may be

obtained by means of an efficient procedure, which should make the total expense for calculating the

4O
ABC's quite acceptable. Note that in our previous work (see refs. 8, 11, and 12) we have used a different version of the algorithm. We expect that the total cost per one Fourier mode in time will decrease for

the current version because of using the thin-layer rather than the full Navier-Stokes equations. (Indeed,

the matrices A_ and B_ in system (57) are 4 × 4 and the matrices for the full Navier-Stokes equations

are 8 x 8.) (See ref. 8.)

Recall, our final goal is to express the values of physical variables at Vl' Uvl - (Uv: Vv_' Pv l' Pv_ ),

in terms of u v, _ ). Choosing the same discretization in time as in the formula (eq. (62)) and implementing inverse Fourier transform (64), we obtain from formula (70),

_U inlx-_

Ulv = Z'r n ,

I

n : -L

(I=0,1 ..... 2L+l)

Then, substituting expression (63) into formula (71) and changing the summation order, we get,

(71)

= _- s

U/v' 2L +-------_ v, --_--_f

n=-L = 1

n = L 2L + I _( _. s'_ . 2re

1 Z Z Tneinl" Uv _-msx-T

-- S

2L+l tUv' -_-_ ?

n=-Ls= 1

= 2L__I(2L_...._.__lln=L_._ein(l-s)x_'_^n'_sT J_'luv , _Uv/--ff_-_ j

s-1 " n=-L

Finally, designating

T l, s = 1 s)'r-- ^

2L + 1 Tn

rl _ -L

(l=0, 1..... 2L+l)

we obtain

Ulvl = Tl, S Uv,-_ )

s=l

(1=0,1 ..... 2L+l) (72)

Equality (72), which is a specification of equality (11 ), provides the missing boundary relations between

the values of the unknowns at F T and at F_r (in the discrete formulation). Therefore, equality (72) provides the ABC's we were aiming to obtain. Additionally, we note that the ABC's (eq. (72)) can be simplified for the case of integrating the Navier-Stokes equations step-by-step in time inside Din. In doing

so, we only need to know u v on the upper time level, i.e., for t = T, which corresponds to 1 = 2L + 1 or

• [ . , ,

to 1 = 0 because of the peno&oty. Substituting l = 0 into equality (72), we obtain

Uvl, = Z T0'S v'--_ J (73)

=T s=l

Equality (73) is a desirable global ABC for implementation together with the step-by-step integration

procedure in time. Indeed, formula (73) expresses the values of u, v, p, and p (perturbations) at the

41
outermost coordinate row vl on the upper time level t = T as a function of the prescribed data Iu v, _-_v /

through the time-periodic solution of the linearized thin-layer equations with the free-stream boundary

condition at infinity. We note that the matrices of operators T o, s are calculated explicitly and therefore,

the practical implementation of the ABC's (eq. (73)) is reduced to a few matrix-vector multiplications.

We also note that this practical implementation may preliminarily require some interpolation in time at

the nodes v, which may be represented in a matrix form as well. If we use an explicit scheme for integrating the Navier-Stokes equations inside Din, then we directly implement formula (73) at each time

step for determining the missing values of the unknowns at the outermost coordinate row of the grid on

the upper time level. If the scheme inside Din is implicit, then we include the relations (eq. (73)) into the

l

system of equations we solve on the upper time level treating all T O,s_Uv' s --10_ fors<2L+l as forcing

terms.

4. Concluding Remarks

We have constructed the DPM-based nonlocal artificial boundary conditions for computation of

oscillating external flows, specifically, compressible viscous fluid flows past finite bodies. To develop

the ABC's, we used linearization of the governing equations against the constant free-stream background in the far field. To justify the constructions of difference potentials used for computation of the

ABC's, we provided some results on solvability of the linearized thin-layer equations. The nonlocal

nature of the proposed ABC's arises from their closeness to the exact boundary conditions. In spite of

this nonlocal nature, our ABC's apply to artificial boundaries of irregular shape with equal ease, which

is very important for applications. We expect that these boundary conditions may become an effective

numerical tool in practical computing. The numerical results on the implementation of the described

ABC's will be presented in future work.

Note that we describe the algorithm for calculating the ABC's only for a particular class of methods

used for integrating the Navier-Stokes equations inside Din, namely, for such methods that the knowledge of missing relations between only two external coordinate rows of the grid (v and Vl) is sufficient

for closing the discrete system inside the computational domain. Obviously, once the method used

inside Din is of higher (than the second) order, the consideration of only two curves, F and F I, might be

insufficient. However, we always can assume that the "linear region" Dex contains more than one curve,

e.g., 1-"1 and F2 instead of only F 1, and can treat this case in the same way as described in this paper.

Moreover, one can use higher order schemes for solving the linearized exterior problem as well. Such

modifications may extend the possible range of applications for the described technique by including,

for example, some computational problems of aeroacoustics.

NASA Langley Research Center

Hampton, VA 23681-0001

April 1I, 1996

References

1. Givoli, Dan: Non-Reflecting Boundary Conditions. J. Comput. Phys., vol. 94, May 1991, pp. 1-29.

2. Givoli, Dan: Numerical Methods for Problems in Infinite Domains. Elsevier, 1992.

3. Gustafsson, Bertii: Far Field Boundary Conditions for Time-Dependent Hyperbolic Systems. SIAM J. Sci. Star. Comput.,

vol. 9, no. 5, Sept. 1988, pp. 812-828.

4. Gustafsson, B.: The Choice of Numerical Boundary Conditions for Hyperbolic Systems. J. Comput. Phys., vol. 48,

Nov. 1982, pp. 270-283.

42
5. Gustafsson, Bertil: Inhomogeneous Conditions at Open Boundaries for Wave Propagation Problems. Appl. Numer. Math.,

vol. 4, Mar. 1988, pp. 3-19.

6. Carpenter, Mark H.; Gottlieb, David; and Abarbanel, Saul: Time-Stable Boundary Conditions for Finite-Difference

Schemes Solving Hyperbolic Systems: Methodology and Application to High-Order Compact Schemes. J. Comput. Phys.,

vol. 111, no. 2, Apr. 1994, pp. 220-236.

7. Seifen, A.; Daraby, A.; Nishri, B.; and Wygnanski, I.: The Effects of Forced Oscillations on the Performance of Airfoils.

AIAA-93-3264, July 1993.

8. Ryaben'kii, V. S.; and Tsynkov, S. V.: Artificial Boundary Conditions for the Numerical Solution of External Viscous Flow

Problems. SIAMJ. Numer. Anal., vol. 32, no. 5, 1995, pp. 1355-1389.

9. Ryaben'kii, Viktor Solomonivich: The Difference Potential Method for Some Problems in Continuum Mechanics.

lzdatel 'stvo Nauka, (Moscow), 1987.

10. Ryaben'kii, V. S.: Boundary Equations With Projections. Russian Math. Surv., vol, 40, no. 2, 1985, p. 121-183.

11. Tsynkov, S. V.: An Application of Nonlocal External Conditions to Viscous Flow Computations. J. Comp. Phys., vol. 116,

no. 2, 1995, p. 212-225.

12. Tsynkov, S. V.; Turkel, E.; and Abarbanel, S.: External Flow Computations Using Global Boundary Conditions. AIAA J.,

vol. 34, no. 4, Apr. 1996, pp. 70(0706.

13. Jameson, A.; Schmidt, Wolfgang; and Turkel, Eli: Numerical Solution of the Euler Equations by Finite Volume Methods

Using Runge Kutta Time Stepping Schemes. AIAA-81-1259, June 1981.

14. Swanson, R. C.; and Turkel, Eli: A Multistage Time-Stepping Scheme for the Navier-Stokes Equations. AIAA-85-0035,

Jan. 1985.

15. Swanson, R. C.; and Turkel, Eli: Artificial Dissipation and Central Difference Schemes for the Euler and Navier-Stokes

Equations. AIAA-87-1107, June 1987.

16. Ryaben'kii, V. S.: Exact Transfer of Boundary Conditions. Comp. Mech. Solids, Issue 1, 1990, pp. 129-145.

17. Calderon, A. P.: Boundary-Value Problems for Elliptic Equations. Proceedings of the Joint Soviet-American Symposium on

Partial Differential Equations, Novosibirsk, Aug. 1963, p. 303-304.

18. Anderson, Dale A.; Tannehill, John C.; and Pletcher, Richard H.: Computational Fluid Mechanics and Heat Transfer.

Hemisphere Publ. Corp., 1984.

19. H6rmander, L.: Linear Partial Differential Operators. Springer, 1963.

20. Vladimirov, V. S.: Equations of Mathematical Physics. M. Dekker, 1971.

21. Wolfram, Stephen: MathematicaTM--A System for Doing Mathematics by Computer. Addison-Wesley Publ. Co., Inc.,

1991.

22. Walker, Robert John: Algebraic Curves. Princeton Univ. Press, 1950.

23. Vainberg, B. R.: Some Problems for Hypo-Elliptic Equations Which are Well-Posed in the Entire Plane. Mat. Sb, vol. 62,

no. 104, 1963, p. 186-248.

24. Kolmogorov, A. N.; and Fomin, S. V.: Elements of the Theory of Functions and Functional Analysis, Nauka, (Moscow),

1981.

25. Ryaben'kii, V. S.; and Tsynkov, S. V.: An Effective Numerical Technique for Solving a Special Class of Ordinary Difference Equations. Appl. Numer. Math., vol. 18, no. 4, Oct. 1995, pp. 489-501.

43
!

REPORT DOCUMENTATION PAGE I Fo,_A_oro_,d

I

OMB No. 0704-0188

Pul_ic reporting burden for this collection of ;n_m,=t;on is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sc_rc_,

gathering and maintaining the data needed, and completing and reviewing the collection of information, Send comments regarding this burden estimate or any other aspect of thl

collection of information, including suggestions for reducing this burden, to Washington Headquarters Servmes, Directorate for Information Operations and Reports, 1215 Jefferso

Dev_s Highway, Suite 1204, Aflington, VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Prolect (0704-0188), Washington, DC 20503.

1. AGENCY USE ONLY (Leavebtank) !2, REPORT DATE

August 1996

4. TITLE AND SUBTITLE

Artificial Boundary Conditions

Flows

for Computation

6. AUTHOR(S)

S. V. Tsynkov

7. PERFORM|¢;GORGANIZATIONAME(S)ANDADDRESS(ES)

NASA Langley Research Center

Hampton, VA 23681-0001

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES)

National Aeronautics and Space Administration

Washington, DC 20546-0001

11. SUPPLEMENTARY NO =¢:S

3. Technical REPORT TYPEMemorandum AND DATES COVERED

5. FUNDING NUMBERS

of Oscillating External

WU 505-59-53-01

8. PERFORMING ORGANIZATION

REPORT NUMBER

L-17513

10. SPONSORINGIMONITORING

AGENCY REPORT NUMBER

NASA TM-4714

Tsynkov: NRC-NASA Resident Research Associate, Langley Research Center, Hampton, VA.

12a. DiSTrIBUTION/AVAILABILITY STATEME¢.'T

Unclassified-Unlimited

Subject Category 02

Availability: NASA CASI (301) 621-0390

12b. DISTRIBUTION CODE

13. ABSTRACT (Maximum 200 words)

In this paper, we propose a new technique for the numerical treatment of external flow problems with oscillatory

behavior of the solution in time. Specifically, we consider the case of unbounded compressible viscous plane flow

past a finite body (airfoil). Oscillations of the flow in time may be caused by the time-periodic injection of fluid into

the boundary layer, which in accordance with experimental data, may essentially increase the performance of the

airfoil. To conduct the actual computations, we have to somehow restrict the original unbounded domain, that is, to

introduce an artificial (external) boundary and to further consider only a finite computational domain. Consequently, we will need to formulate some artificial boundary conditions (ABC's) at the introduced external boundary. The ABC's we are aiming to obtain must meet a fundamental requirement. One should be able to uniquely

complement the solution calculated inside the finite computational domain to its infinite exterior so that the original

problem is solved within the desired accuracy. Our construction of such ABC's for oscillating flows is based on an

essential assumption: the Navier-Stokes equations can be linearized in the far field against the free-stream background. To actually compute the ABC's, we represent the far-field solution as a Fourier series in time and then

apply the Difference Potentials Method (DPM) of V. S. Ryaben'kii. This paper contains a general theoretical

description of the algorithm for setting the DPM-based ABC's for time-periodic external flows. Based on our experience in implementing analogous ABC's for steady-state problems (a simpler case), we expect that these boundary

conditions will become an effective tool for constructing robust numerical methods to calculate oscillatory flows.

14. SUBJECT TERMS

Time-periodic flows; Artificial boundary conditions; Boundary projection operators;

Difference potentials method

18. SECURITY CLASSIFICATION

OF THIS PAGE

Unclassified

17. SECURITY CLASSIFICATION

OF REPORT

Unclassified

NSN 7540-01-280-5500

19. SECURITY CLASSIFICATION

OF ABSTRACT

Unclassified

15. NUMBER OF PAGES

44

16. PRICE CODE

A03

20. LIMITATION

OF ABSTRACT

Standard Form 298 (Rev. 2-88)

Prescribed by ANSI Std. Z39-t8

298-102