.. - . 

NASA-CR-166565-PT -1 

1984001~435 

A Reproduced Copy 

OF 

-TX -- - -.-- - :::--

Reproduced for NASA 

by the 

NASA Scientific and Technical"" Information Facility 

\ \11\\\1\ 1\11 \\\\ 11\\1 \\111 \\\1\ \1\1\ 1\\\ \11\ 

FFNo 672 Aug 6S NF02385 

https://ntrs.nasa.gov/search.jsp?R=19840019435 2018-11-29T21:14:32+00:00Z
., 

! 

.. 

• 

, 

.. ·'l ~ 

~!-~j .... .,..'Io:-'.J 

-1.9 1"--.1 

~ ..... ,.,. f~ 

-.,. . ..-l!" ~:::,H"" 

NASA CONTRACTOR REPORT 166555 _ .o.Jl·' 

(UASA-Cn-1bb!>oS-Pt-1) SOLU'l ie ~ OF 

J-Oltll::li$lO:lAL TI1U.-UEHliIH.N'I Vl!i~O(J$ FLlJilS. 

~AllT 1: IhV£STIGAI1Cb Cf CANCILATE ALGO~l~US (3C1unt1ii~ ~e~€dcch A~socidles, 

lnc.) Uo p He AUS/Af A~l C~~~ 12& GJ/bij 

Solution of Thrcc-Di~ensional Time-Dependent Viscous Flows 

Part 1: Investigation of Candidate Algoritl~s 

Bernard C. Weinberg 

lIenr>' ~lcDonald 

CONTRACT NAS2-10016 

Apr 11 1979 

N8~-lI5CJ 

Uu~lds 

1Jubu 

~~ ~ ~t·~ -.... ---.------.-----.. --.--.----.--..-. dr!i.-:.2 7,52l,,;5, 'E. _".'."' __ " ._. ____ . ____ . __ .~ __ 4 .....-.~ ......... I·~ .... -.... • ---- - -.- .. ~- ... -........... _-- _ .. 
i. 

,. 

;~ 

• 

\ 

• 

,~ .... ... 

~ 

.... 

.. .. .~, 

NASA CONTRACTOR REPORT 166565 

Solution of Thr~e-Dirnensional Time-Dependent Viscous Flows 

Part 1: Investigation of Candidate Algorithms 

Bernard C. Weinberg 

itcnry HcDonald 

Scientific Research Associates. Inc. 

Glastonbury, Connecticut 

Prepared for 

Ames Research C~nter 

under Contract NAS2-10016 

Natlonat AeronautiCS and 

Space Administration 

Ames Research Center 

Moffett Field. California 94035 

.-._-.-.." ... _. -- . 
' .. 

TABLE OF CONTENTS 

SUMMARY. . . . . . . INTRODUCTION 

LIST OF SYMBOLS • 

ANALYSIS 

Background • 

Temporal S~heme8 • 

Consistently Split Block Implicit Scheme 

Predictor-Corrector Method 

The Hopscotch Algorithm • 

Spatial Difference Approximations 

Implicit Tridiagonal Finite Differences • 

Basis Functions • 

RESULTS • 

Numerical Experiments 

CONCLUSIONS • 

APPENDICES 

Appendix A 

Appendix B 

REFERENCES 

TABLES 

FIGURES 

Linearization Technique 

Discretization of Model Equations 

------------.. - .... _---- .--------- ~---.. ~ - --_._--- - • _ . ..", ...... 1· ...",..;. ___. . .,-.. • __ ~ 

Page 

1 

2 

4 

7 

7 

13 

14 

17 

19 

21 

21 

32, 

46 

46 

53 

54 

54 

59 

64 

68 

78 
Stw.M.\RY 

The flow over a helicopter rotor blade in forward flight is an important 

example of three-dimensional time-dependent flow. The boundary layers on the 

rotor blade set loss levels and control retreating blade stall. As a consequence there is considerable interest in developing a numerical scheme for 

so'ving the time-dependent viscous compressible three-dimensional flow to aid 

in the design of helicopter rotors. In the present report candidate numerical 

algorithms are examined to determine their overall suitability for the 

efficient and routine solution of an appropriate system of partial differential 

equations. It is concluded that a consistently split time-linearized block 

implicit scheme using either quintic B-spline collocation or the generalized 

operator compact impliCit approach to generate a fourth order accurate 

algorithm is particularly well suited for use on the present problem. High 

cell Reynolds number behavior leads to favoring the generalized operator 

compact implicit approach over the quintic B-spline collocation method. 

1 

-_··· ____ 9.--_~_ ... ___ .,. _________ . ____ ... --.. - ___ "-s. .. ~ _.~ --l. ____ '..... • __ _ 
~, 

.. .;a 

INTRODUCTION 

The behavior of boundary layers on wings and bodies has long been of 

interest to aerodynamicists. In both Rteady and unsteady flows the boundary 

layers are known to govern a major portion of the losses and to significantly 

influence the vehicle lift and moment coefficients. When the flow is steady, 

boundary layer prediction schemes based on numerical solution to the governing 

partial differential equations of motion have reached a high level of sophistication and predictive accuracy, even in three space dimensions. In 

unsteady flows, such as are commonly encountered in rotary winged aircraft, 

some progress has been made in two space dimensions but little to date has 

appeared on unsteady three-dimensional boundary layers. 

Two particular problems arise with time-dependent three-dimensional 

boundary layers relative to the steady case. The first of these is the 

rather obvious one of time integration with its added requirements of transient 

accuracy coupled with an increase in the computational labor. The second of 

these is the so-called negative cross flow problem, which to some extent has 

troubled the steady boundary layer prediction schemes. Kendall et al., 

(Ref. 1) discuss the negative cross flow problem for steady three-dimensional 

boundary layers in a very illuminating fashion. This particular problem 

arises when the spanwise component of velocity changes sign and will be 

discussed in detail subsequently. Because of the interest by external 

aerodynamicists in swept wing boundary layers ~here the negative cross flow 

problem (in this case flow" from tip to root) is not usually encountered, the 

negative cross flow problem has not received a great deal of attention to 

date. However in transien"t flows, particularly those encountered on rotor 

blades in forward flight, negative cross flows are frequently encountered. 

For instance, the advancing rotor blade has cross flows of one sign during 

the first ninety degrees of rotation Gnd these can change sign over part of 

the blade during the second ninety degrees. 

Thus to be of practical value. time-dependent three-dimensional boundary 

layer prediction schemes require high computational efficiency and transient 

accuracy coupled to the ability to treat arbitrary cross flow profiles. 

These attributes are not available in any existing available computer code 

and hence in view of the potential use for a code of this type its development 

2 

- --------------"---:"., 
is appropriate and timely. In this report the initial phase in the development 

of an efficient time-dependent three-dimensional boundary layer prediction 

procedure is investigated, namely, the choice of the computational algorithm 

and spatial differencing technique. 

3 

-_____ 1L .. 'U.JIII 
". 

~ .. 

~,. 

::~ ~ : ,_ ..... 

: ~ ~ . ,~ 

, 

:'..: 

a 

A1 -A6 

b 

B1- B6 

Bi 

Bki 

c 

d 

Dy 

Dz 

D+,D_,DO 

:.b x 

:.by 

:.b z 

F 

Fj 

G 

h 

H 

I 

J 

J~!AX 

K 

L 

L(u) 

1 

n 

N 

Ni 

P 

- c + qj,qj,qj 

Q 

- c + rj,rj,rj 

R 

LIST OF SYMBOLS 

Coefficient in Eq. (6) 

Coefficients in Eq. (54) 

Coefficient in Eq. (6) 

Coefficients in Eq. (55) 

itlt B-spl1ne 

ith B-spline of order k 

Coefficient in Eq. (6) 

Coefficient in Eq. (6) 

Spatial operator - auyy + cuy 

Spatial operator - buzz + duz 

First order difference operators 

Spatial operator in x direction 

Spatial operator in y direction 

Spatial operator in z direction 

Function defined in Appendix A 

First derivative approximation at point j (Eq. 12) 

Function defined in Appendix A 

Spatial step size 

Function defined in Appendix A 

Number of intervals in x direction 

Number of intervals in y direction 

Maximum number of intervals 

Number of intervals in z direction 

Spatial operator 

Spatial operator 

Operator defined in Eq. (A7c) 

nth time level 

Number of equations 

Nonlinear spatial operator 

Pressure 

Difference weights in Eq. (15) 

Operator defined in Eq. (17) 

Difference weights in Eq. (15) 

Operator defined in Eq. (17) 

4 

'-, ---~.----'-'--'- ' -- - -- .- ,~ .: / \ 
.; 

, 

I 

" 

. 

, , 

" 

'1 

w 

y 

z 

a,Il,y,IS,£: 

ai 

aO,al 

Il 

Il i 

6t 

).m 

). .. z 6t/6z2 

II 

v 

v,a 

p 

T 

~ 

111 

LIST OF SYMBOLS (CONT'D) 

Cell Reynolds r~er 

Spline approximation to u 

Second derivative approximation at point j (Eq. 13) 

Source term 

Time 

Truncation error terms, Eq. (23) 

Streamwise velocity component 

Velocity components 

Normal velocity component 

Spanwise velocity component 

Weighting function 

Streamwise direction 

Surface normal direction 

Spanwise direction 

Coefficients in model equations (Eq. 56 and Eq. 57) 

Coefficient of B-spline 

Boundary condition values 

Parameter to allow centering of time step 

Coefficient of B-spline 

Time step 

Parameter in Eq. (25) 

Viscous stability parameter 

Roots of characteristic equation (cf. Eq. (28» 

Exponent defined in Eq. (23) 

Coefficients in Burgers eQuation 

Density 

Truncation error 

Vector of unknowns with components ~l' ~2' ~3 

~ _ ~n _ increltent value 

5 

- ---_._----_ .. _---.-. - .. ,~-------.-------.---------.~~--.---- ---------.. 
• 

~"~. 

;J 

,. , 

, 

-:-v 

1 ;. 

.. 

'~~. 

.~ ". 

,'.~ 

~-

· . 

· '.-. -! ~. 

.- . 

~ 

.: 

· : 

.J 

.--. ) 

~~ 

.. ::-.... ~ 

, 

Superscripts 

(1) 

m 

n 

n+1 

-,c,+ 

Subscripts 

j 

m 

LIST OF SYMBOLS (CONT'D) 

Order of derivative 

Defined in Eq. (25) 

nth time step, t a n6t 

(n+l)st time step, t a (n+l)6t 

Indicates where rand q coefficients are evaluated 

Primes indicate derivatives 

Terms in Eq. (9) not divided by a(x) 

Index - designates grid pOint 

Defined in Eq. (39) 

6 
.; ',:., 

'.~] ".: '>. .• 

... 

~. 

.. <0 • 

ANALYSIS 

Background 

In this section we discuss the requirements of a three-dimensional 

unsteady boundary laye.r code and demonstrate how the physics of the flow 

influences the choice of a computational technique. 

Three-dimensional boundary layers occur on the wings and fuselages of 

both conventional and rotary wing aircraft. In both types of vehicles, the 

boundary layers are importnnt in settin~ loss levels and determining useful 

operating ranges. As is well known, boundary layers are sensitive to pressure 

gradients. In time-dependent flow the temporal acceleration terms appear in 

the momentum equation in a form very similar to the conventional imposed 

pressure gradient and so for qualitative evaluation purposes can be regarded 

as 'pseudo' or 'auxiliary' pressure gradients. Viewed in this manner the 

temporal acceleration terms are likely to be able to influence quantities of 

practical importance such as skin friction, displacement thickness and the 

onset of separation. At the range of frequencies typically encountered in 

rotary wing aircraft aerodynamic problems, it is clear, for instance, from 

the very thorough review of ~!cCroskey (Ref. 2), that very significant transient 

boundary layer effects can be observed • 

In examining the flow problems of practical interest such as loss levels 

or the onset of separation it is evident that all three space dimensions must 

be considered. In conventional aircraft the sweep effect is of interest and 

inherently three-dimensional. In rotary wing aircraft in forward flight 

clearly very substantial transient changes occur in what might be termed the 

local sweep angle. Howeve·r generally speaking, the boundary layers remain 

·thin unless catastrophic flow separation occurs or the flow at the wing or 

rotor tip is considered. As a consequence it might be supposed that the usual 

three-dimensional thin boundary sheet approximations (Nash and Patel, Ref. 3) 

could be used to produce a valid set of governing equations. Fortunately some 

improvements in thin boundary sheet approximations are possible as a result 

of having to treat the negative cross flow problem mentioned earlier. 

The negative cross flow problem is best explained in a somewhat intuitive 

manner, and for steady boundary layers a very good physical description of the 

problem is given by Kendall et al., (Ref. 1). Looking at the suction surface of 

a conventional swept back wing the boundary layer cross flow, w, is usually 

7 
.~. 

>.1 

'. ··1 

outward in the z positive direction along the span from root to tip. Thus 

conventional steady boundary layer integration schemes have developed by 

forward marching the streamwise velocity u in the streamwise x direction and 

simultaneously marching out along the span in the z positive direction. In 

view of the physics of the problem, the spanwise marching scheme does not 

normally encounter negative w, i.e., spanwise inflow. This is very fortunate 

because it is difficult, indeed it could be argued impossible, to structure a 

physically satisfactory unconditionally stable scheme which permits forward 

marching in the spanwise direction with a negative w cross flow. At least 

j.ntuitively the problem of negative cross flow implies information being 

transferred upstream against the spanwise marching direction. Conventional 

stability analyses confirm the inability to forward march into regions of 

significant negative w. From experience with attempts to march the twodimensional boundary layer equations into a region of separated flow and its 

obvious relationship to the negative cross flow problem, it is not surprising 

that spanwise marching into a negative cross flow region is not accomplished 

without special treatment, for instance the Krause "zig-zag" scheme (Ref. 4). 

Recently conventional boundary layer developers have been turning to a spanwise as well as normal implicit formulation to remove the restriction of only 

positive cross flows (Kendall et al., Ref. 1). With a spanwise impliCit 

formulation spanwise diffusion is allowed, and the resulting implicit system 

of equations can be treated by direct elimination (Ref. I), by a predictorcorrector iterative approach (Ref. 5), or by the process of matrix splitting which 

reduces the matrix elimination labor (Refs. 6 and 7). Lin and Rubin (Ref. 5) 

in their predictor-corrector boundary region solutions for flow over a yawed 

cone at moderate incidence showed that allowing diffusion in the spanwise 

direction not only eliminates the problems associated with negative cross flow, 

but improves upon the solutions obtained by three-dimensional boundary layer 

techniques. Again intuitively a spanwise implicit construction permits 

information transfer in either direction. Boundary conditions applied at the 

tip can influence the flow inboard, if required by the physics of the flow. 

For these reasons the implicit spanwise construction has been a feature of 

the three-dimensional duct flow analysis of Briley (Ref. 6) and McDonald and 

Briley (Ref. 7). Based on the experience in Refs. 6 and 7, the additional 

computational effort resulting from a spanwisc implicit formulation could be 

8 
,-1 

-J 

... 

as little as a 2D% increase relative to the explicit spanwise marching approach. 

The extension of the steady three-dimensional boundary layer equations to allow 

spanwise diffusion is easily accomplished, and in view of the improved physical 

representation which thus follows, it is recommended. As a matter of course it 

has been assumed that normal to the wall an implicit formulation would be 

structured. In recent years for boundary layer type problems there has been 

little dispute as to the efficiency gains to be had from an implicit formulation normal to the wall (Ref. 8). However in the streamwise direction for 

steady two-dimensional flow, the equations are normally forward carched and 

the implicit stability obtained entirely from being implicit in th~ normal to 

the wall direction. 

For unst~ady boundary layers the problem can again be formulated in 

either an explicit or implicit manner. As with spatial marching of steady 

boundary layers for usual aerodynamic applications, the locally refined 

spatial mesh required to define a (turbulent) boundary layer influenced by a 

(transiently) varying pressure distribution, when an explicit (stability 

restricted) scheme is employed, results in a maximum time step that is much 

less than the time scale of the physical processes of interest. Thus for 

solving unsteady boundary layers of the type usually encountered in rotary 

winged aircraft an implicit formulation is desirable. Since in time-dependent 

flow diffusion in the streamwise direction is normally negligible due to the 

usual boundary layer approximations, it is possible to formulate an implicit 

time-dependent scheme that retains the implicit structure in the spanwise 

and normal directions (which was found desirable for the steady boundary 

layer) and march th~ solution in the streamwise direction. 

As mentioned earlier the streamwise marching sweep would probably require 

less computational effort by about 20% than a fully .implicit formulation and 

of course less storage. However since the solution is being time marched, the 

opportunity to use a streamwise implicit formulation at roughly the same cost 

as the streamwise marching sweep does arise. If one does perform a streamwise 

marching sweep, then the linearization of nonlinear terms is performed about 

the known spatial marching level. If a fully implicit structure is adopted, 

then full time linearization can be utilized. That is the linearization of the 

nonlinear terms is performed about the known time level. As is pointed out in 

Ref. 7, it is easier to obtain a consistent spatial-temporal order accurate 

9 
, .. 

linearization by marching in time than in space (in time the nonlinear marching 

derivatives have the form PUi whereas in space marching they have the form 

puiuj ). Further by structuring implicitly in the space marching direction, 

(small) regions of axial reverse flow would be permitted. As a result of 

these combined benefits of linearization and separation, a fully implicit 

structure is advocated. 

Transient calculations mean that, in essence, a full three-dimensional 

spatial integration is carried out at each time step. Thus, spatial accuracy 

is very important to minimize the spatial grid point density for efficiency 

since many time steps are contemplated in a given cycle. In order t~ get the 

most out of a given spatial difference formula, the errors from repLesenting 

nonlinear terms by linear combinations of terms should be less than or equal to 

the spatial discretization errors. If the linearization introduces a greater 

error than the spatial differencing, then either a coarser spatial mesh could 

be used, or iteration, or some form of linearization improvement is called 

for. Iteration across a time step is not recommended since this only reduces 

the linearization error and computationally costs as much as a complete time 

step. Cutting back the time step would be preferable to iterating to preserve 

the linearization error at some acceptable level, since cutting back on the 

time step would improve both' the transient error and the linearization error. 

To obtain a linearization, which introduces errors of at most the same as the 

spatial difference formulae, a Taylor series expansion about the known time 

level can be performed. This process clearly demands a formal block, i.p-., 

coupled, treatment of the system of equations. For instance in the streamwise 

momentum equation a typical t~rm is linearized: 

n+1 and clearly one cannot lag w at the old time level n without introducing a 

n+l first order time error in order to get an uncoupled system, i.e., w not 

appearing in the streamwise momentum equation., Thus forcal linearization and 

consideration of the resulting errors indicate the coupled system ough~ to be 

treated from the accuracy point of view. This is further reinforced when it 

is realized that direct elimination of block, i.e., coupled, banded systems 

are not computationally expensive compared to the iterative solution of an 

10 
, 

'k 

1. 

1 

-l 

.. j 

.~ -:-; 

'."> '. ~,;: 

"::~ -;t 

"·'~~L 

,~ . ~,. ,".-

':?,I .;';..~,"": 

::';'\' 

· -".,.-

.. ~ .. f'" .... 

.. --:'~ 

,'~ .. ~. '. ~;r'(~ n :> _-.• -1' 

.. '.~.II'" 

· ~ .~~;,~. 

-:'~~ . ' .. ~" 

· '~ :.~·.l~ -'; ,':~~h-. 

t·jti 

'I .. 

uncoupled system. The linearization technique is described in detail in 

Appendix A. 

Additionally a second type of approximation arioes unconnected Yith 

linearization but arising from baoic coupling terms in the original equations 

and if indeed some terms in an equation arc time lagged in order to uncouple 

the equation system and these terms are of equal importarce to tbe terms 

retained, then again an iterative updating is called for in order to achieve 

stability, accuracy and consistency. (This could be termed ad hoc equation 

uncoupling). Blottner (ReC. 8) has shown that trulny iterZltior.:I around the 

ad ~ uncoupled 3et (>10) io sometimes required in order to achieve an overall 

solution accuracy commensurate Yith the local diffc:rence molecule accuracy. 

Tho~pson and liacDonald (Ref. 9) found in a three-dimensional momentum integral 

procedure that a lagged sequential iterative type of cnlculation would not. 

in a number of instances, even converge. The weight of opinion definitely 

favors the block coupled approach. 

As a general observation, care is required' to obtain acceptable transient 

accuracy for long tim~ integration with conventional finite difference schemes. 

A Crank-Nicolson cent~red time implicit scheme for instance, although second 

order in time, shows quite a dispersion problem (relative to other schemes) on 

the simple pure convection problem. However the problem of transient accuracy 

is significantly reduced in the typical boundary layer problem since the time 

dependency io continuously input through initial and boundary conditions and 

relatively the concern is Yith 'short' time integrations. The computational 

problem is more of what the phase lag of the wall shear is. relative to the 

prescribed free streao disturbance, than concern over tbe convection velocity 

of a wave in a shear after a long propagation time. The interest is in forced 

oscillations with a minimum scale of the boundary layer thickness over a few 

cyclcs of thc !!lotion, just enough to obtain rcpctitioncycl1cally. It is 

thercfore expected that a significant dispersion problea vill not arise with 

a conventional impliCit scheme. 

The equation system which will be considered is forcally of block size 

four, consisting of the continuity and two moaentua cquations and an 'energy' 

equation for p, u, v, 1.1 (p is specified everywhere). If constant stagnation 

temperature is assumed, p, u and ware related by an algebraic equation and 

the problea can be reduced to a block-three syste~ rather than block-four upon 

11 
... 

option, with significant reduction in computer time. Matrix splitting 

techniques such as those described in Refs. 6 and 7 have considerable potential 

to reduce the computational labor of solving the block implicit system of 

(linear) algebraic equations which result from discrctizing the governing 

equations. ·Schemes of this general type are termed aplit linearized block 

implicit or for brevity split LBI schemes and arc reviewed in detail by 

Briley and McDonald (Ref. 10). With a careful ordering of the sweeps with 

a split LBI scheme it is possible to use a block-two with the third equation 

uncoupled o~ two of the sweeps so that this would be a potential major 

advantage of a split LBI approach. 

The ultimate goal of the "optimum" scheme is to diminish both storage 

requirements and running times in order to achieve a desired accuracy level. 

Although schemes can be constructed to satisfy either one or both of these 

goals, the robustness of a method can only be verified by considering its 

applicability to a general class of problems. Our concern here is with an 

approximate form of the three-di~ensional unsteady compressible Navier-Stokes 

equations, so that the method chosen will by necessity be required to treat a 

coupled system of nonlinear partial differ~ntial equations. At the outset the 

following observations can be ~de concerning the characteristics of the method. 

1. Implicit methods (preferably in all three spatial directions) are desired 

in order to eliminate stability restrictions and permit solutions with both 

positive and neRative streamwir.e and spanwlse velocities. 

2. The nature of the nonlinear coupling of the variables in the governing 

equations require that the equations be solved coupled. 

3. Iteration Ghould be a~oided and "ticp. linearization" procedures employed 

(a discussion of this point is given in the following section). 

4. The method should allow for general boundary conditions. 

5. The method should allow the flexibility of incorporating higher order 

spatial differencing methods. 

Although Item 5 pertains to npatial differencing, in the next section we 

will demonstrate that it also plays a crucisl role in choosing the type of 

temporal scheme. It would thus appear that first an overall temporal discretization procedure must be chonen and only then an efficient spatial schece that 

could be incorporated with it. 

The arguments above, linearization, stability considerations and physics 

lead us to suggest that the current state of the art dictates that the reco=ended 

scheme Ghould be in the framework of a time linearized block impliCit method. 

12 
• 

. 

'\ 

" 

,~ 

~ 

.~:>:1 ""',. .... t'~ 

With the foregoing as background indicating the recocmended overall 

approach to the problem of developing an efficient accurate computer code to 

predict time-dependent viscous flows, attention is turned to the specific problem dealt with in this report. This report concerns itself with first which 

overall numerical algorithm to adopt and secondly which spatial differencing 

scheme to employ or what basis of approximating solution functions to use. 

In the next section we will disCUGS three methods for solving general 

systems of cultidimensional parabolic equations, the consistently split block 

implicit scheme, Rubin'a Predictor-Corrector technique and the Hopscotch 

algorithm. The following two sections will describe two spntial differencing 

methods, IoIhich loIe classify as Q-R operator schemes and bas10 functiclll schemes 

(which deal ~1inly with B-splinea). Applications of these methods to model 

two point boundary value prob1e=~ and to n coupled cyctc= of two nonlincar oncdimensional parabolic equations possessing n "three-dimensional boundary 

layer-like behavior" arc then given. 

Temporal Schemes 

Consider a system of three-dimensional nonlinear parabolic differential 

equntions 

(1) 

where ¢ is a vector of unknowns (9 1, ¢2' ¢3). With the equations appropriately 

linearized, at each time step, a sy~tem of N - 3(1-1)(J-l)(K-l) linear equations 

result, where I, J, K are the nu:nber of intervals in the x, y, z directions, 

respectively. Direct inversion of the system is not practical in three space 

dimensions except for extremely coarse meshes since the operation count is 

3 * proportional to N /3 Although higher order spatial schemes allows one to 

reduce the number of grid points, resultG of model problems indicate that for 

the range of accuracy desired, one can expect at best a reduction of a factor 

of four in grid points in each coordinate direction. As significant as this may 

seem, it may not appreciably affect the overall computation unless the original 

problem is reduced to a more tractable form. 

~For I - J - K • 10 

For I • J - K • 50 

N a 2187 and the matrix inversion operation count is 

proportional to 1010. 

N - 352947 and the matrix inverRion operation count is 

proportional to 1016 • 

13 
' . 

. , 

• 

., . 'i . 

" 

The usual procedure io to tranaform the original problem into a sequence 

of simpler problems or into a sequence of one-dimensional problems. We 

consider the linearized version of Eq. (1) (see Appendix A for a description of 

the time linearization procedure) with the differential operators identified 

with their coordinate direction, namely 

(2) 

where superscripts indicate the time level, i.e., t

n - ~t, Sn is a source 

term and 0 ~ S ~ 1 is a parameter allowing one to center the time step, i.e., 

B - 0 corresponds to a forward difference, 8 • 1/2 to Crank-NicolGon and 

8 • 1 to a backlJard difference. 

Consistently Split Block Implicit Scheme 

We will first consider the consistently split block implicit scheme. 

Solution of Eq. (2) is accomplished by application of a generalization to 

systems of PDE's of an alternating-direction implicit (ADI) technique for 

parabolic-hyperbolic equations. The original ADI method was introduced by 

Peaceman and Rachford (Ref.'ll) and Douglas (Ref. 12); however, the 

alternating-direction concept has since been expanded and generalized. A 

discussion of various alternating-direction techniques is given by Mitchell 

(Ref. 13), Yanenko (Ref. 14) and more recently by Briley and McDonald (Ref. 10). 

The present technique Is simply an application of a generalization of 

the procedure developed by Douglas and Gunn (Ref. 15) for generating conaistently split ADI schemes as perturbations of fundamental implicit difference 

schemes such as the backward-difference or Crank-Nicolson schemes in its 

natural extension to systems of partial differential. equations. In this 

context a consistent scheme is one where the intermediate levels represent 

a discrete approximation to the governing equations whoae truncation error can 

be made to vanish as the time or spatial mesh io arbitrarily reduced. Consistency 

in this sense is a very valuable property as it can greatly simplify the 

accurate implementation of boundary conditions (Ref. 10). 

For the present, it will be assumed that ~(¢) contains derivatives of first 

and second order with respect to the coordinate direction, but no mixed derivatives. Mixed d~rivatives arc allowable within the formal framework but unless 

they are i~portant they are best treated explicitly (lagging) or by extrapolation. 

14 
The Douglas-Gunn representation of Eq. (2) can be wri~ten a9 the following 

three-step solution procedure: 

(3) 

In increment form Eq. (3) reduces to the algorithm given by Briley and 

McDonald (Refs. 10 and 17) for solving the compressible time-dependent thrcedimensional Navier-Stokes equations, viz., 

(4) 

(I _ f3llt~)( ~*** _ cp"> • cp** _ ~n 

where ~ * and ¢ ** are intermediate solutions. Each of Eqs. (4) can be written 

in narrow block-banded matrix form and solved by efficient block-elimination 

methods. If ¢ * and ¢ ** are eliminated, Eqs. (4) become 

(I - {36f:Dx HI - (36f:Dy )( I - {3M:DzHtP n+1 - tPn> • llt [ (:Dx + ~y + :Dz> cpo + SO] 

(5) 

15 
.. 

··.~1< ... 

~. " 

I 

I 

\ 

t 

t 

i 

~ 

·1 

I 

... 

If the multiplication on the left-hand side of Eq. (5) is performed, it 

becomes apparent that Eq. (5) approximates Eq. (2) to order (6t)2. Although 

the stability of Eqs. (4) has not been established in circumstances sufficiently general to encompass the Navier-Stokes equations, it is often suggeste~ 

(e.g., Richtmyer and Morton, Ref. 16, p. 215) that the scheme is stable and 

accurate under conditions more ge~eral than those for which rigorous proofs 

are available. This latter notion was adopted her~ as a working hypothesis 

supported by favorable results obtained in actual computations (e.g., Refs. 

17, 18 and 19). 

Several observations can be made concerning Eqs. (3) and (4). 

1. The system of three-dimensional equations has been reduced to three 

systeus of one-dimensional equations. 

2. The inversion of the total system is now at most approximately 108 

(IJ + JK + IK) operations compared to approximately (3IJK)3 operations for the 

direct inversion problem. 

3. The first step i.nvolves at least 40% - 50i. of the operations (Ref. 7). 

4. The method does not have a eFL stability condition and for S ~ 1/2 is 

von Neumann stable. 

5. The method is applicable to rectangular domains. Although it does not 

necessarily preserve symmetry along diagonals of rectangular domains, the 

procedure can be corrected if so desired (Ref. 20). 

A major attraction of the Douglas-Gunn scheme is that the intermedIate 

* ** +1 solutions ¢ and ~ are consistent approximations to ~n • Furthermore, 

n * ** n+l for steady solutions. ¢ a ¢ a ¢ a ¢ independent of 6t. Thus, physical 

bOlJndary conditions for 4>n:-l can be used in the intermediate steps without a 

serious loss in accura:y and with no loss for steady solutions. In this 

respect. the Douglas-Gunn scheme appears to have an advantage over locally onedimensional (LOD) or "splitting" schemes. and other schemes whose intermediate 

steps do not satisfy the consistency condition. The lack of con3istency in 

the intermediate steps complicates the treatment of boundary conditions and. 

according to Yanenko (Ref. 14. p. 33). does not permit the use of asymptotically large time steps. 

It is worth noting that the operator ~ can be split into any number of 

components which need not be associated with a particular coordinate direction. 

As pointed out by Douglas and Gunn (Ref. 15). the criterion for identifying 

16 
·., 

l 

.. 

" 

.. 

"-

'J 

j 

-i 

H 

Cj 

./ .' ..- ,,/ .~ 

.- ,~ 

.>1 

fi 

j 

.. ~ 

.~ 

'" .~. 

" Ii 

".i 

~, 

~l 

'!i 

" 'i 

; i 

f 

~: 

h 

r 

L ." . 

I 

, 

t 

I 

~ 

.. 

lr 

f 

;,. 

3. 

'J 

;'1 

,:! 

" ., 

Q 

., 

f: 

~. 

, 

;'1 

':1 ,-I 

;j 

,. 

:r1 

;1 j 

,. 

J 

. CII • 

sub-operators is that the associated matrices be "easily solved" (i.e., narrowbanded). Thus, mixed derivatives can be treated implicitly within such a 

framework, although this would increase the number of intermediate oteps 

and thereby complicate the solution procedure. Finally, only minor changes 

are introduced 1f, in the foregoing development of the numerical method, ~, 

and S are functions of the spatial coordinates and time, as well as ~. 

Alternative techniques to the consistently split linearized block implicit 

scheme presented here have been proposed as general algorithms. We discuss 

two such methods, a predictor-corrector scheme due to Rubin (which we will 

refer to as pIC) and the Hopscotch algorithm of Gourlay and his co-workers. 

Although these authors have had success with these methods, it will be shown 

that in order to meet the requirements of the problem under consideration, 

they are not as versatile or as efficient as the consistently split linearized 

block implicit scheme previously discussed • 

Predictor-Corrector Method 

A predictor-corrector method has been successfully employed by Rubin and 

Lin for three-dimensional viscous flows in which diffusion is important in 

two directions (Refs. 5 and 21). Their objective in developing a compromise 

between explicit techniques-and implicit methods (AUl) was to eliminate viscous 

stability restrictions and to minimize CFL stability limitations common to 

explicit methods and to reduce the total work per time step by eliminating one 

of the block tridiagonal inversions required in the usual ADl procedure. In 

addition they desired a "symmetric" method that would easily treat symmetry 

conditions and other derivative boundary conditions. 

A comparison of ADl with Rubin's predictor-corrector method can be 

obtained by conSidering the model two-dimensional linear parabolic equation • 

The pIc method reduces to the following expression 

{I - {3ll1D y + 2'X. z{3b} [u** - un] • ll1{ Dy + Dz}U" 

+ .B 6tD z[u*-u"] + 2'X.z.Bb(u*-u") 

17 

(6) 

(7) 
'J,j 

11 

il 

.. ! 

!l·1 !I 

!~ 

·:1 F·"", 

V d 

. ~ 

~j 

1j 

'j! 

;1 

'~.j 

I! 

U .; 

~i !'I 

~: 1;; 

ji 

~. 

Jt: 

'I 

r 

'~ 

;. 

~, 

.y, 

. ~ 

/ 

. . 1 . 

• 

_ ." __ 'r. _ __ ~ __ .~_ .. _ ~~~~ _ .. _v~_ .... ~/ .~--. --~"':;"'i ... ,,.,,..-----. 

* ** where u denotes u at a previous iteration, u denotes u at the latest 

iteration,Dy

• au

yy 

+ cu

y , Dz a bUzz + dUz and Az • At/Az2. As a first 

guess one of the following extrapolations is used; 

linear replacement 

O(At2) Taylor series 

u* • un 

u* - 2un _ un- l 

The use of linear replacement achieveD consistency only with many iterations 

or small time steps. 

The sequence of steps for solving Eq. (7) by the pIc method are: 

1. In Eq. (7) compute At(Dy + Dz)un• 

2. Compute BAtDz(u* - un) and 2A Bb(u* - un) using u* from a previous z 

iteration. 

( ** n 3. Compute the coefricient matrix of u - u ). 

4. Obtain (u** - un) (requires one tridiagonal sweep). 

5. Repeat steps 2-4 until convergence. usually two to three iterations. 

Note that even for linear problems. iteration Is necessary to obtain the 

desired accuracy. Thus at a minimum two to three tridiagonal sweeps and two or 

three explicit evaluations of the terms in step two are required • 

By comparison the Douglas-Gunn ADI procedure gives for the combined two 

steps 

* ** where u corresponds to an intermediate solution and u to the solution at the 

(n+l)st time step. Note that Eq. (8) differs from Eq. (7) in appearance only 

in' the underlined terms in the former equation. 

The sequence of steps in the Douglas-Gunn procedure for the two-Jimensional 

problLm is as follows: 

1. Compute the right-hand side to be used in th~ first sweep At[Dy + Dzlun• 

2 • Compute (u* - un). which requires one tridiagonal inversion, i.e., 

3. 

4 • 

The 

(I - BAtDz)(u* - un) n At(Dy + Dz)un • 

Compute (u** - un) which requires one tridiagonal inversion. i.e., 

(I - BAtDy)(u ** - u n) * n a (u - u ). 

Evaluate un+l = (u** - un) + un. 

major effort is expended by both methods in evaluating the term 

18 
.' ... . -, .-, 

... -~----- r «FS;:,ilR44 P41· ~. 

6t(Dy + Dz)un • Neglecting the evaluation of the coefficient matrix of 

(u*· - un), two passes of pic is equivalent to one double sweep of ADI. Hence 

by comparing the Gcquence of steps for both Pic and ADI we see there is no 

apparent advantage to the predictor-corrector method. Furthermore, due to the 

"explicit" nature of the Pic method in the "z" direction it possesses a CFL 

stability condition while ADI coes not. 

If we were to consider higher order spatial approximations in both 

coordinate directions, each operator evaluation could involve a matrix inversion, (see section on spatial approximations). Hence the Pic ~ethod could be 

more costly than an ADI procedure. Finally the pic methcd, to the authors' 

knowledge, has not been formulated for a full three-dimensional problem (with 

diffusion in all three directions, such as is required for the present problem) 

so that its applicability under such conditions is unknown. We see no 

advantage to the pic method and therefore do not recommend it for the present 

problem. 

The HopscotCh Algorithm 

The so-called family of Hopscotch algorithms have been advocated by 

Gourlay and h:l.s coworkers (Refs. 22, 23 and 24) for the solution of multidimensional parabolic equations. Several variants of Hopscotch exist; for twodimensional problems there is the fully explicit odd-even Hopscotch which 

resembles the DuFort-Frankel algorithm and the partially implicit line and ADI 

Hopscotch procedures which resemble the Peaceman-Rachford ADI method. These 

methods have lower operation counts than ADI schemes due to their partial 

explicit nature that, depending on the scheme, eliminates some or all of the 

matrix inversions. For instance, line Hopscotch requires only half the number 

of matrix inversions of a comparable ADI computation. However, the Hopscotch 

methods have stability restrictions and are only first order accurate in time. 

Numerical results of model linear scalar parabolic equations (Ref. 23) 

when the stability conditions are not violated, confirm the above conclusions, 

i.e., Hopscotch is more efficient than ADI. However, recent results for the 

driven cavity problem (Ref. 25), which requires the solution of a coupled 

system of nonlinear equations, lead to contrary conclusions. The Poisson 

equation for the stream function was solved separately by a direct method 

19 
" ':;-, " .. ",.: .. : ..... ..- ~.... "., . . d'::tW, . • 

. ,. Q ;A . $(( $. ; $,,;:4' • > r;C;W<4J444" f' , ... ii. ;." tt.Akp¢tDlibt-¥S;:,¥,.*UU .• , 

} 'l:': ~ •. _ .• _. 

) 

1 

(Ref. 26). while the vorticity transport equation was solved by both ADI and 

odd-even Hopscotch. The ADI solutions were 20% faster than Hopscotch to 

obtain a converged steady state solution. Although Hopscotch requires fewer 

operations per time step than ADI. the CFL stability condition necessitate~ 

the use of a smaller time step which results in the noted increased runni~& 

times. 

Hore recently. Greenberg (Ref. 27) has added a number of new members to 

the Hopscotch family for three-dimensional parabolic problems. As with the 

two-dimensional procedures. these new members also have restrictive stability 

conditions. To date there has not been widespread use of Hopscotch-type 

schemes. in particular for coupled nonlinear parabolic equations. and their 

viability under such circumstances is still an open question. 

In general. Hopscotch owes its favorable characteristic to what Gourlay 

terms E-operators. of which the standard tridiagonal finite differences is a 

member. For other higher order spatial differencing. where the operators are 

tridiagonal but not E-operators. e.g •• evaluations of spatial operators of the 

form auxx + bux which even for an explicit temporal scheme involves a matrix 

inversion, the method loses most of its desirable features (see section on 

spatial approxications). In addition the ability to handle coupled implicit 

boundary conditions is also not as flexible with Hopscotch. Finally. there 

have been claims that one of the attributes of Hopscotch is its ease of 

programming. While this may be true for model problems, for more complex 

problems, i.e., coupled systems of three-dimensional nonlinear parabolic 

equations. the basic logic for setting up the block inversions. whether in one 

direction as in Hopscotch,.or in three directions as in split LBI. is 

comparable. Hence. in this case. programming considerations should not greatly 

influence one's choice of method. 

The lack of versatility of the method. i.e •• stability restriction. 

inability to incorporate higher order spatial methods, and results of more 

realistic problems (Ref. 25) lead us not to recommend the Hopscotch algorithm. 

20 
'. 

~ . 

Spatial Difference Approximations 

Implicit Tridiagonal Finite Differences 

Q-R Operator Notation 

In this section implicit tridiagonal finite difference approximations to 

the first and second derivatives and to the spatial differential operator 

will be considered. The very versatile Q-R operator notation will be introduced. which allows as special cases a variety of schemes such as standard 

second order finite differences. first order upwind differences. fourth order 

operator compact implicit (OCI). fourth order generalized OCI and exponential 

type methods. Since all these schemes are of the same form. a single subroutine which defines the difference weights is all that is required to 

identify the method. while leaving the basic structure of the program 

unaltered. Subsequently. the results of numerical experiments for a number 

of these schemes will be presented. 

The Q-R formulation allows for ADI methods and permits the treatment of 

systems of coupled equations. i.e •• LBI methods. Although variable mesh 

schemes can be employed within the Q-R framework. it is believed preferable 

to use analytic transformations to obtain a uniform computational mesh. hence 

attention is restricted to uniform mesh formulations. 

The general concepts and notation for two point boundary value problems 

will be introduced and then the methodology entended to more general linear 

and nonlinear parabolic partial differential equations in one dimension. The 

extension to multidimensional problems will also be indicated. 

Consider the two point boundary value problem 

L(U) cC(X)Uu + b(X)Ux + c(X)u -'f(x) (9) 

with u(O) and u(l) prescribed. Derivative boundary conditions. although not 

treated here. can easily be incorporated into the framework of the Q-R 

operator notation. Let the domain be discretized so that Xj - (j-1)h. 

j - i. 2 ••••• J + 1. and Uj ~ u(xj ). Fj ~ ux(xj ). Sj ~ uxx(xj ) and h = l/J 

1s the meah width. The numbering convention was chosen here to be compatible 

with FORTRAN coding. 

21 
Without loss in generality for a(x) ~ 0, Eq. (9) Can be divided by a(x) 

so that we may treat instead the following equation 

L(u) • uxx + b(x)ux + c(x)u • f(X) (10) 

where 

b(x) • b(x) /c(x), c(x) • c(xl/o(x) Orid f(x). f(x)/o(x) 

The spatial differential operator is identified as 

(11) 

Substituting the finite difference approximations to the first and second 

derivatives 

(12) 

(13) 

into Eq. (10) and rearranging, we obtain 

or 

(14) 

where RCj = hbj is the cell Reynolds number. 

22 
_ " 1':.-

Equation (14) can be generalized by introducing operator format, i.e., 

(15) 

where the superscripts (-) minus, (c) center, and (+) plus indicate the 

difference weight that multiplies the variable evaluated at the (j-l), (j) and 

(j+l) grid points, respectively, and where the rj's and qj'S for grid point j 

are functions of h, bj _l , bj , bj+l , c

j

_l ' c

j and c

j

+l • Comparing Eqs. (14) and 

(15) we can identify the rj's and qj'S, viz., 

r j a I - RCj 12 qj • 0 

C 

• h2Cj - 2 C • I rj qj 

+ + rj al+RCj/2 qj -0 

We now define the tridiagonal difference operators Q and R 

R[Uj ] arjUj_1 +rjuj+rjuJ+1 

a[ fj] a qjfj_1 + q~fj + qjfj+1 

Noting that L(u) = f and substituting Eq. (17) into Eq. (15) we obtain 

. -1 Alternatively by employing the inverse operator Q an exp~ession for L(u)j 

can be obtained 

23 

(16) 

(17) 

(18) 

(19) 
.·.r· .. ··~ "<;:B~i~'4:;'~;~:;~~~'~:""'" .. 

. ~., , 

. ,. 

c . 

;-

:? 

f~ 

r 

f' 

, I 

~. 

-1 For standard central finite differences Q ~ Q = I, the identity matrix, 

(the spatial nperator is given explicitly in terms of Uj_l , Uj and Uj+l ) so 

that nothing was gained in obtaining Eq. (19). However, in general, for higher 

order methods Q is tridiagonal and Q-l is a full matrix. Hence Eq. (19) gives 

us a means of expressing the spatial operator for a wider class of difference 

'approximations. The formalism in Eq. (19) is also applicable for first and 

second derivatives appearing alone (cf. Ref. 28). It must be pointed out however 

that Eq. (19) is not the most general formulation since the compact implicit 

formulas cannot be combined to yield a single scalar equation relating the 

spatial operator to the function values (Ref. 28). 

In the next section a method due to Berger et al., (Ref. 29) is described 

that enables one to construct fourth order tridiagonal methods with certain 

desirable properties, i.e., evaluate the qj and rj coefficients. 

Generalized Operator Compact Implicit Schemes 

Given 

L(U) = Uxx + b(x)ux + c(x)u 

an expression relating L(u)' and u is sought in the form 

(20) 

where T

j is the truncation error and Q and R are tridiagonal displacement 

4 operators. The maximu~ accuracy attainable is fourth order, i.e., T

j ~ O(h ). 

Expanding Eq. (20) in terms of q-'c,+ and r-'c,+ we obtain, j j 

I ( 1[_ c +] Tj = h2 RUj - a LU)j • ~ fj UJ

-1 + rj Uj + rJ Uj +1 

- [qj(LU)j_1 + qj(LU)J + q;(LU)j+l] 

(21) 

24 
.j 

',j 

'.'1 

J 

1 

! 

j 

J 

'.'1 

., 

J 

1 

, 

J 

-'-,',-,,',1 

'1 

oJ 

A Taylor series expansion yields for tj 

·" 

ORIGINI\L p:,~::: 1:-

OF POOR QUi\LIr{ 

T) • TJUeX)} +T)IU(I)(lJ) + Tfu'l!)(x J, + TJ3u(~)(xJ' 

+T"U(4)(X )+T~u(!5)( x, + T6u(6'(X} + O(h~) j J 1 ) ) ) . 

(22) 

where superscripts in parentheses denote derivatives with respect to x. and 

where 

o I [( _ e +, 2 - C • >1 TJ • V rJ + rl + rJ + h (QJ CJ-I+ QJCJ + QJc j +1 J 

TJ

I • + [crt - rj>- h (QibJ-1 + Qj b J + Qjbj+I' - hZ(Qj'J+I- QjCJ-I'] 

(23) 

. T v V-Zl' + v _ h ( + v-I - J • h 7! (rJ + (-I) r,) - (v-I)! Qj bJ+1 + (-t) Q) b)-I' 

- IV~2)! I Q j + I _Ov Qj) + h'lqjCj" + I-lIvQjcl_,) I V' 3,4,5,6 

For second order central finite differences we net TO • Tl • T2 • O. This 

c - + yields, when qj • land qj - qj a 0, the following relations 

e (- ., hZ( _ c +, r

J • - r) +r) + QJ'j-1 +qjC) + qJCJ+1 

rt -rJ

- • hbJ • RCJ 

rt+rj-2 

25 
~- -

: I 

'. 

I , 

• I' 

\ 

I, 

" I. 

i; 

,. 

, 

L 

~ 

t 

• 

~ 

i 

] 

j 

which recovers Eq. (16), i.e., 

ORIGIN!',:' r-,~,:,: " ~~ 

OF POOH Q'J, .!.Il·{ 

r

C • -2 + hZc J j 

r

J

- • 1- RCj 12 

rt · 1 + RCJ/2 

To obtain the fourth order operator compact implicit scheme we again set 

o 1 2 -c+ T • T • T DO to obtain three expressions for r

j " in terms of the 

-c+ - + c qj' , , (note that q ,q ; 0 and q is not necessarily unity), i.e., 

(24a) 

(24b) 

3 4 - c + 

now T and T must be set to evaluate qj' , • The standard Swartz OCI 

3 4 method requires T • TaO, 

I (+ _) I [ + _] [+ _] "6 r - r -"2 RCj+IQj + RCj_IQj - Qj - Qj • 0 (24d) 

- h2 CQjCJ+1- QrCj-I) 

1 C + - ) 1 [ + _] I [ + 24- r +r - 6" RCj~IQj - RCj-IQj -"2 Qj + QJ-] • 0 

- h2(Q{Cj+1 + QTcj - l ) (24e) 

26 
--

- (5) - (6) 4 and results in a leading truncation error term of (Au + Bu )h. Substi-

+ - tuting (24b) and (24c) into (24d) and (24e) rj and rj can be elioinated and a 

- c + system of two £quations in qj and qj with qj as a parameter is obtained. The 

parameter q~ is proportional to the determinant of the system. The values of 

-c+ -c+ 

q • • and r •• are presented in Table I. j j 

As shown in Ref. 28. a cell Reynolds number stability condition exists 

for the Swartz OCI scheme. i.e •• for Rc ~ If:Z nonrealistic or oscillatory 

solutions will be obtained. 

relax the conditions T3 • T4 

to be of O(h4). 

In order to eliminate this restriction one can 

• O. and allow the coefficients of u(3) and u(4) 

di -,c,+ iii R By expan ng qj n a ser es n c 

(25) 

. - c + 12 parameters. ~m" m - O. 1. 2. 3. arc introduced. The equatio~s for 

T3 • O(h4) and T4. O(h4) yield 5 linear relations, leaving for disposal 6 

"free" parameters plus a factor. 

These parameters can be set according to some criteria that would yield 

certain desirable properties for the difference equations. The following 

constraints arc prescribed 

and h is sufficiently small so that 

IObJ-bl_l-bj+1 > 0 and hCj+l/bl

+1 < 2 

for j-2, .. ·,Jond cISO 

27 

(26) 

\ 
·~ t. 

I. 

... 

• 

~ . , ~. 

~ (~ 

r 

I

:· . . / .. , . .... " 

These conditions assure that R is diagonally dominant and Q is invertible for 

all Re' Further details arc given in Ref. 29. The significance of this 

approach is that one can construct 8 scheme possessing certain desired properties employing a set of preassigned rules. This is contrary to usual practice, 

in which a scheme is chosen, and then its properties are determined. It is 

important to note that the derivation of the q and r coefficients is not a 

trivial task. In computational effort it is also not cheap. This point will 

be discussed in greater detail in 8 later section. 

The q and r coefficients for the generalized OCI scheme described in 

Ref. 29 are given in Table II. In this report comparisons are made with another 

generalized OCI sche~e. whose coefficients are given in Table III. Numerical 

experiments indicate that these two schemes are comparable. differing only in 

the magnitude of the truncation error • 

Exponential Type Schemes 

Another family of schemes that can be expressed in Q-R operator notation 

are the so-called exponential methods. The idea. originally due to Allen 

(Ref. 30) (independently derived by Il'in (Ref. 31) and McDonald (Ref. 32» 

and employed by Dennis (Ref. 33). is to set the difference weights so that the 

numerical solution is equated to the analytic solution for the locally frozen 

constant coefficient equation. Allen (Ref. 30) and Il'in (Ref. 31) considered 

the homogeneous constant coefficient equation. 

(27) 

so that the difference approximation was set identically equal to the analytic 

solution. 

The analytic solution of Eq. (27) has the form 

-by 

u • A + Be 

where A and B are determined from the boundary conditions. 

The Q-R operator formulation of Eq. (27), 

28 

(28) 
\. 

, 

! 

, : . 

r 

@ 

rr 

f

; 

. , 

, 

": 

; 

- I 

, 

. 

I 

I 

~ 

.. 

I 

i 

I 

! 

i 

I , 

.1 

i 

I 

"; 

I 

• l 

possesses a fundamental solution of the form ~j. Solving for p and employing 

Eq. (24a) we obtain 

With the aid of Eq. (24b) and setting ~j to the analytic solution e-bjx' we 

obtain 

jL' • e-bjX • e- bjhj • e- Rc, • rj Ir,-

+ - (Tl With - c,+ the condition r - r a Rc • 0 in Eq. 23), we can define r' as j j j j j 

r-. Rc e-RCj I (I _ e-RCj) 

I 

rj • RCj I (I - e -RCj) (29) 

- + c and q • q D 0 and q • 1 where we have allowed b

j to vary. Equation (29) can 

be rearranged to yield an alternate form 

where D+, D_, Do are the forward, backward and central first difference operators, respectively. 

This method is second order accurate for Rc ~ 0(1) and becomes first order 

accurate as Rc + ~ where the scheme reverts to first order upwind differencing 

- + c (in Eq. (28) rj + 0, rj + RCj , rj +-RCj for bj > 0 to give RCj(Uj+l - Uj ) - 0). 

Another exponential scheme which is uniformly second order accurate was 

developed by El-Mistikawy and Werle (Refs. 34 and 35). The "exponential bex 

scheme" which is incorporated in their Dolution of the boundary layer equations 

with strong blowing, 1s based on a spatial operator of the form given in 

Eq. (10). Berger et al., (Ref. 36) derived the counterpart for an operator of 

29 
;t 

J '! 

the form given in Eq. (27). The q and r coefficients are presented in 

Table IV. Although this scheme reverts to second order upwind differences as 

Rc + -. it does not possess a maximum principle analogous to the ordinary 

differential equation it is approximating as does the exponential scheme of 

Allen (Ref. 30). 

The Allen exponential scheme can be applied to differential equations of 

the type given in Eq. (10). Substituting the Q-R representation for L(u). we 

obtain 

(31) 

Multiplying through by Q. and combining terms. the difference approximation to 

Eq. (10) becomes 

or 

(32) 

Note that the only difference between Eq. (32) and Eq. (18) (whe:-e c

j :: 0) 1:~ 

the coefficient ~~trix multiplying Uj • Hence the methodology is unaltered. 

Application of Q-R Operator Schemes to l-D Parabolic Equations 

Consider the one-dimensional linear parabolic equation * 

U t ~ auxx + bux + cu + d (33) 

with appropriate boundary conditicns and initial conditions. where a. b. c. d 

could be functions of x and t. Dividing by a ; O. identifying L(u) • u + bu • 

xx x 

* If Eq. (33) were nonlinear, assume it was linearized by the method described 

in Appendix B. 

30 

, 
I 

I 

I 

I 

I 

I 

I 

r 

-! 

.0 ~ _ 0. .. '. 

and employing a two point temporal difference scheme. we obtain 

where B • 1 reduces to a backward difference and B • 1/2 to the Crank-Nicolson 

approximation. 

Letting 

n+f3 I [ n+f3]-f [ n+f3] L • h2 a R 

and with some algebra we obtain. 

(34) 

where A • 6t/h2

• 

Again, Eq. (34) is general and permits second order finite differences. 

upwind differencing, exponential type schemes. and OCI schemes. Aside from the 

evaluation of the Q and R operators. the problem is no more complicated than 

standard second order finite differences. Results of model problems can be 

found in Ref. 28. 

Application to Coupled Nonlinear Parabolic Equations 

Given a system of m nonlinear parabolic equations in m unknowns. 

.r ml 

1=1 

j·I.2 ..... J+1 

31 
.. : . 

III • 

n+6 where Ni is a nonlinear nonconservative spatial operator. the Q-R formalism 

carries directly over provided that for any equation only one independent 

variable is operated upon by the differential operator. For example. 

o(U,w,v) Ut • Uxx + b(U,v,w)ux + c(u,V,w) 

is allowed since x derivatives of u only appear, while 

( ) Ut • Uu + b(u,v,w)ux + c(u,v,w) + d(u,v,w)w a U,VI,V x 

is not allowed since x derivatives of both u and w appear. The modified 

unsteady Navier-Stokes equations for the three-dimensional time-dependent 

boundary layer. when written in quasilinear form. fall within the class of 

allowable differential operators. Thus for the problem being addressed 

in the present study the OCI schemes are applicable. 

Multidimensional problems and/or more general equation forms can usually 

be accommodated by a splitting procedure, which reduces the differential 

operator to a sequence of one-dimensional problems which have the appropriate 

allowable form. However, as with standard finite differences. special procedures must be applied to cross derivative terms. e.g •• extrapolation or lagging 

at the previous time step or increasing the nuu.ber of intermediate steps in the 

splitting. 

In Appendix B an example of the above procedure is presented for a 

coupled system of one-dime~sional nonlinear parabolic equations possessing 

"boundary layer-like reverse flow" behavior. 

Basis Functions 

A convenient approach for developing numerical procedures for the solution 

of partial differential equations is founded on a basis function representation 

of the dependent variable, i.e •• it. one dimension, 

u ( X) • I a· B· (x) I I I 

32 
1 

" .. 

where u(x) is represented by a linear combination of suitably chosen basis 

functions, ni(x), and the ai are determined by some set of constraint relationships. The advantage of the basis function representation is that once the ai 

are determined, one can obtain at very little additional cost the function 

and its derivatives anywhere in the domain. 

The major concern here will be with the B-spline basis functions. However 

by way of comparison some of the characteristics of piecewise Lagrange polynomials and piecewise Hermite polynomials will also be given. In the following 

basis functions are sought, polynomials of order k (degree k-l), that possess 

certain smoothness properties, and the computational efficiency of these 

functions for the solution of differential equations are investigated. 

The Simplest functions are the piecewise Lagrange polynomials, which can 

be computed by the cardinal basis functions. Consider a grid numbered from 

j • 1 to J + 1 (to be consistent with FORTRAN coding) so that one is considering J + 1 grid points and J intervals. Over each interval co~sisting of k 

knots the basis functions are polynomials of degree k-l which are equal to 

unity at one particular knot and zero at the other knots. The dimension of 

the piecewise Lagrange polynomials (the number of basis functions) is J + 1 

(independent of the degree of the polynomial) and thus J + 1 constraints 

must be satisfied to determine all the a i • However, in order to evaluate a 

function and its derivatives at a particular location, only k basis functions 

distributed over k adjacent grid points are required. In addition Lagrange 

polynomials are only CO at the end points ~f each sub-interval, again 

independent of the degree of the polynomial and thus allows for j~ps in the 

first derivative there. Due to the purely interpolatory character of piecewise Lagrange polynomials the a i are nothing more than the functions evaluated 

at the appropriate knots. 

The interpolation problem reduces to an explicit procedure as does the 

evaluation of the derivatives at the internal knots. When incorporated with 

the method of collocation (the analog of interpolation for the solution of 

ordinary differential equations), the piecewise Lagrange basis functions used 

locally on a uniform mesh (with k odd) recovers the standard centered finite 

difference approximations. Thus a quadratic polynomial leads to a three point 

formula and the inversion of a tridiagonal matrix of order J + 1 and a quartic 

33 
,,4)$C $ ( 

------.. _ ..... --~ .. ~-_. \ J? "'W.,...,-r<;:~ --:'~"'!T:;j -','-

,.)..0. , 

-;\ .. 

~ 

:~ 

...... , 

''1. 

1 

1· 

s ..\ 

.: 

0, 

r 

] 

• 

I 

t 

) 

• ! 

~ ? 

? 

.. 

l+ 

,T 

\ :)-

->. :-..., 

~~ :'1 

t'?j 

polynomial leads to a five point formula and the invcraion of a pentadlagonal 

matrix of order J + 1. 

In order to obtain smoother representations additional constraints must be 

applied on the derivatives of the approximating functions. This leads to 

more complicated basis functions. For example, the family of Hermite polynomials requires that for a (2m - 1) degree polynomial em ~ 2) on two adjac~nt 

knots, m interpolatory constraints must be satisfied for the function and its 

em-I) derivatives. The approximation is still local, i.e., given a function 

and an appropriate number of derivatives on two adjacent knots, the function 

and its derivatives can be computed implicitly in the interior of that domain. 

However, the solution of a differential equation by the method of collocation 

involves the inversion of a matrix of order m(J + 1) of bandwidth (3m - 1) 

(- n + m where n = degree of polynomial). Thus in the process of obtaining a 

Cn- 2 representation a substantial increase in labor in comparison to Lagrange 

polynomials has been incurred. Further details can be found in Ref. 37. 

B-Splines 

Another family of functions are the polynomial splines, i.e., polynomials 

of degree k-l that arc Ck- 2• One would expect that the additional smoothness 

of the spline approximation would translate into a better behaved solution -

perhaps even more accurate - for "smooth" functions, than a Lagrange or 

Hermite polynomial of the same degree. Numerical experiments have shown this 

to be the case in many instances. However, the additional smoothness 

constraints that must be specified not only change the nature of the approximation, but may add considerably to the computational effort, depending on the 

manner in which they are applied. 

There are several approaches for representing polynomial splines. The 

first method involves specifying the polynomial and the smoothness conditions .. separately, and then solving the imposed constraint equations simultaneously 

for the undetermined coefficients. This approach has succesafully been 

.. The constraint equations may either be interpolatory conditions or a 

differential equation. 

34 

- -
i 

I 

./ ~ I 

/ 

.r 

/ 

! 

employed by Kendall and Bartlett (Ref. 38) in their solution of the chemically 

reacting boundary layer equations, and Murphy et al., (Ref. 39) for the solution of the boundary layer equations. With the aid of a skillful partitioning 

of the resulting banded block matrix, the above authors were able to solve 

for the function values directly. However. by so doing they lost ready access 

to the interpolatory polynomial representation, which is one of the acknowledged benefits of a spline technique. An alternate approach described in 

Ahlberg, Nilson and Walsh (Ref. 40) and adapted by Rubin and Khosla (Refs. 41 

and 42) in their solution of the boundary layer equations also employs a 

polynomial formulation. However, the additional continuity and smoothness 

conditions are used to transform their equations into expressions relating the 

function values and their derivatives. Their "spline equations" are solved 

simultaneously with the appropriate constraint equations, which entails for a 

single scalar differential equation the inversion of 2 x 2 block tridiagonal 

matrix to obtain fourth order accuracy and a 3 x 3 block tridiagonal matrix 

to obtain sixth order accuracy. 

In order to avoid inverting block banded matrices that arise even in 

the solution of scalar differential equations a third alternative is 

considered, the B-spline representation, which has the appropriate smoothness 

built into the functions themselves. For the solution of the coupled NavierStokes equations this procedure has a definite advantage in that the order of 

the block submatrices, typically three or four, can be reduced by one. A 

detailed discussion of B-splines, their construction and their mathematic~! 

properties are given in Refs. 37 and 43 and will not be repeated here. 

Instead, the properties of.n-splines that make them attractive for the solution 

of partial differential equations will be emphasized. In the following 

discussion, without loss in generality, a uniform mesh will be considered. 

The normalized B-spline (the sum of the basis functions at any knot is 

unity) of order k (degree k-l) is bell-shaped and spans k+ 1 knots. Within 

the interval the B-spline is positive while outside the interval, including 

the end points, it is zero. For knots of multiplicity unity the first k - 2 

derivatives of the B-spline are zero with a jump possible in the (k - l)st 

derivative. In general, for a knot of multiplicity k - v the vth derivative 

of the B-spline is discontinuous at the point with all lower order derivati'Tes 

continuous. Due to the compact support property of the B-splines at most k 

35 
,,+. 

./ 

/ 

~.' 

0< AiP!JWWOi.,. ',4'" ' .. 4~'i'-4W4 ···t. P;' (~.e=t='A" t· i .. A;,,~" '''-, t •. . 4iZ .J paft i ) y: .S sq;W;'fel1ij#kiR!'!ifiJA:'J'!CtWPti-, 

B-splines are nonzero over any interval. while at a knot k - 1 n-splines are 

nonzero. This property enables one to reduce the bandwidth of the matrices 

that are obtained in the solution of differential equations. while retaining 

the scalar structure. 

DeBoor (Ref. 44) has shown that for B-splines of order k the dimension 

of the basis is J + 1 + k - 2. where it is assumed that the knots are of 

multiplicity unity and there are J + 1 nodal points (J intervals) in the 

domain. Thus for cubics (k - 4) there are J + 3 independent basis functions 

and for quintics (k a 6) there are J + 5 basis functions. Hence J + 1 + k - 2 

conditions are required to fully specify the B-spline representation. In view 

of the above. for a B-spline of order k. in order to determine all the a i • the 

resulting matrix is scalar. is of order J + 1 + k - 2 and has a bandwidth of k. 

After the a i are det~rmined. the evaluation of the function and its k - 2 

derivatives each require at most k additional multiplications. 

For the pure interpolation problem 

(1.)( ) • f(e) U XI I t· 1 .... ,K-3 

U(X J

) • f

J -I, ... ,J+I 

u(t.)(X ) _ f(£) 

J+I .1+1 R -I, ... , K-3 

where superscript ~ designates a derivative of order 1. and k even. Prenter 

(Ref. 37) shows that a unique function exists. the polynomial spline of order 

k. What is apparent from the statement of this problem is that the B-spline 

interpolating polynomial is global in character. being "tied together" by 

the derivative constraints at the end points. In contrast to piecewise 

Lagrange and Hermite polynomials where coefficients were determined for each 

subinterval by the appropriate interpolatory constraints. the evaluation of 

B-splines requires the inversion of a matrix. In addition the ai's that are 

determined are not equal to the values of the functions and its derivatives 

at the knots 85 was the case for the Lagrange and Hermite representations. 

It is in part due to this global character that B-splines derive their 

smoothness prop~rties. 

36 
-; .. , : .... :-

DeBoor (Re!. 44) has presented a computer program to determine the 

normalized B-splines of order k < 20, which allows for nonuniform meshes and 

for discontinuities in lower order derivatives at the knots (knots ~f multiplicity> 1). He has also conditioned the B-splines near the boundaries to 

simplify the setting of boundary conditions. Alternatively the procedure 

suggested by Prenter (Ref. 34) could be employed (for a uniform mesh with 

simple knots), but that would require additional algebraic manipulation to 

modify the coefficient matrix so that one can obtain the appropriate band 

structure. Since both methods are equivalent (the uniqueness of the B-splinc 

representation) the choice is a matter of convenience. Prenter's approach 

employs the same basis functions throughout the domain, which reduces the 

memory and co~puter logic requirements but increases the algebraic manipulation, in particular for higher order splines. The more general DeBoor 

approach requires more computer memory, but reduces the preprocessing required 

of the user. In the subsequent work DeBoor's approach was used. 

The accuracy of the B-spline representation is now considered. For a 

function u that is sufficiently smooth, with derivatives of order k that are 

k continuous, u c C (a,b), the B-spline representation of u, s (of order k, 

k even), Prenter gives the following error estimates 

Hence, to solve a second order ODE by collocation using cubic splines (k = 4) 

2 

we can expect O(h ) accuracy but require only a simple tridiagonal matrix 

elimination. However, to obtain greater accuracy larger bandwidth matrices 

are required. Consideration of this point is given in the following section. 

Application of B-Splines to the Solution of Differential Equations 

Since the B-spline representation will be employed as a spatial approximation, its properties, order of accuracy, efficiency and spatial stability 

37 
.', 

. ' ... ~~'('a." .. 

behavior can be investigated by considering its applicaticn to the solution 

of two point boundary value problems. 

In this section B-spline basis functions of order k are applied to 

the solution of ordinary differential equations of the form 

L(U) • Uxx + bux + cu • f 

with boundary conditions 

(35) 

Several different approaches can be employed for the solution of Eq. (35), 

depending on the degree of accuracy desired and the level of complexity one is 

willing to accept. The techniques of collocation, Galerkin, subdomain, 

least squares, etc., can be viewed as special cases of the method of weighted 

residuals (of Crandall (Ref. 45) for instance). More recently Murphy (Ref. 46) 

has shown the relationship of these methods to orthogonalization processes 

and has characterized them as generalized Galerkin techniques. Herp. we will 

review the relationship of the different approaches to the method of weighted 

residuals and indicate how. they translate into computational effort. 

The basic idea of the method of weighted residuals is to choose a trial 

function for the independent variable 

u(x)· ~aj8ki(x) 

I 

. where Bki are the basis functions (here B-splines of order k) and a i are 

coefficients to be determined. The derivatives that are given by 

u' ( X) • 2: a j B ~i ( X) i 

u"(x)· 2:aiB~i(x) 

i 

are substituted into Eq. (35) to yield 

38 

(36) 

(37) 
'. 

,,;-

.... ~ 

.. 

L (u) • L t a I Bkl ( x ). t a I L Bkl ( x ) 

L(u)· tal {B~I(X) + b( x) B~I(X) + C(X)Bkl(X)} • f(x) 

(38) 

We now require that the reoidual error in some sense vanishes over the domain. 

i.e •• 

(39) 

where Wm 1s a sUitably chosen weighting function. The simplcBt method 1s to 

force the residual to vanilih at each of the nodal points 1n the clomain. This 

method of collocation is analogous to the interpolation problem which requires 

the approximating function nnd/or its derivatives to be equal to the true 

function at the nodal points. 

In this case W is the Dirac delta function and we recover the system m 

of equations 

or 

(40) 

For B-splines of order k there are J + 1 + k - 2 basis functions 

(dimension of the B-spline basis) so that in addition to the collocation 

relationships at each of the nodal points. including the boundary points. k - 2 

p'Jpplemental conditions arc -required. For cubic B-spline (k • 4) the two 

boundary conditions arc sufficient to clone the system. lienee a linear system 

of equations of order J + 3 and bandwidth k - 1 • 3 (tridiagonal) is obtained • 

Note that the problem is no cere difficult than the standard second order 

finite difference case, except that J + 3 equations arc solved instead of at 

39 
• 

most J + 1 equations • . k-2 The overall accuracy of the method is O(h ) which 

corresponds to the accuracy of spline approximation of the second derivative. 

Once the ai's are deter.nined, u, u' and u" can be obtained from Eq. l.17). 

For higher order splines, e.g., k • 6 (0(h4) accuracy) an additional k - 2 .. la 

conditions are required to close the system, which reduces to specifying two 

Clore conditions than nre available, i.e., the two boundary condition!>. If one 

has information concerning the behavior of the differential equation at the 

boundaries (values of derivatives thereh then that can be used. However, 

in general this is not nlways the case, so it was decided to collocate at two 

points that were nonco1ncldent with the knots. These locations were set at 

x a h/2 and at x • I - h/2 solely for ease of computation (modifying the 

matt"ix), but insight into the behavior of the differential equation could be 

used to choose the collocation pointD, i.e., regions of rapid change. It is 

important to note that thesa collocation p~ints are located in a subregion 

between two knots. There is no need to add grid points (knots). 

The bandwidth of the resulting matrix has nOll been increased to 5 in 

order to obtain an 0(h4) solution. If the coefficients, b, c and f arc not 

known at x - h/2 and x - h/2 then .these values must be obtained by interpolation. 

It can bp shown that more nccurate results cnn be obtained (one order 

better) if derivatives arc introduced as unkno_~s, i.e., 

A system of two equations in 01 and a

i are thus obtained 

and 

40 

(41) 

(42) 

(43) 
...... -_ .. 

1 

! 

~, 

., 

However, a 2 x 2 block matrix is now required to be inverted and boundary 

conditions on u' must also be supplied. 

Calerkin Method 

If we choose the test function Wm to be from the same family as the trial 

function, then the weighted residual becomes analogouG to an orthogonality 

condition. TIlis leads to the following relationship, 

(44) 

and we obtain (J + 1 + k - 2) equations which is exactly equal to the dicension 

of the basis. However, the bandwidth of the coefficient matrix for the ai's 

is increased to 2k - 1 (7 for cubics and 11 for quintics). This result can be 

obtained by focusing our attention at the B-spline test function centered at 

knot j - 11:. It will contribute to the integral in the interval (m - k/2, 

m + k/2) where it is nonzero. However, the B-spl:nes centered in the interval 

[m - (k - 1), m + (k - I)} will also contribute to the integral. Hence the 

Bki's will span (k - 1) + (k - 1) + 1 • 2k - 1 knots and will contain 2k - 1 

entries. 

The Galerkin scheme can be shown to be O(h4) accurate for cubic B-splines. 

However, as compared to the method of collocation the bandwidth has increased 

from k - 1 to 2k - 1 so that the procedures would be much more costly. Furthermore, the Galerkin scheme requires four nontrivial integrals which necessitates 

an integration scheme consistent with the order of the method, e.g., for cubics 

a fourth order Simpson's rule. Since the B-splines are not orthogonal in the 

sense that integrals of products of basis functions are not zero, no simplifications exist. 

Since in general one is not solving a problem that stems from a variational 

formulation, one ~ould be justified in searching for test functions from a 

different family which could simplify the integrations and reduce the bandwidth 

of the reSUlting matrix. One such technique, the subdomain method (Ref. 45) 

or what Murphy (Ref. 46) terms the generalized Calerkin procedure, employs a 

unit step function as the test function. This formulation yields the following 

system of equations 
.---.'. 

· 

, 

I' 

J • 1,2, ... ,J (45) 

Since only J equations are derived for the (J + 1 + k - 2)ais. k - 3 additional 

relationships are required. For cubics. this leaves. aside froo the two 

boundary conditions. one as yet unspecified condition. Murphy. (Ref. 37) and 

Bartlett and Kendall (Ref. 38) in their viscous flow solutiona specified an 

additional derivative boundary condition at the edge of the viscous layer. 

Although their condition was exact mathematically for the problem they 

considered. for oore general equations such conditions could overconstrain 

the solution. Hence it is felt that collocating at some point. i.e •• one of 

the boundaries. would be helpful and would not deteriorate the order of 

accuracy of the solution (this point is further discussed in the section on 

nuoerical results). 

The bandwidth of the resulting matrix has now been reduced to k as 

compared to 2k - 1 for the standard Galerkin scheme while still retaining 

4 O(h ) accuracy. However. the integration scheme (now over two adjacent grid 

points) needs to be appropriately handled in order to achieve the desired 

accuracy of the method. In order to employ Simpson's rule for cubic B-splines. 

values of the coefficients b. c and f are required at the intermediate points 

Xj + 1/2' If they are not known analytically. interpolations would be 

required which would increase the total comput3tion time. 

The advantages of the cubic B-spline generalized Galerkin procedure are 

its decreased blndwidth of four and its applicability to treat equations in 

conservative form. However. it has the disadvantage of being sensitive to 

the type of boundary condition used to satisfy the extraneous condition 

(cf •• section on numerical results). 

The advantages of the quintic B-spline collocation procedure are that 

boundary conditions are easy to apply and it possesses the lovest truncation 

error for a given grid spacing. However. it is less efficient than the 

generalized Galerkin method due to its larger bandwidth of five. At the 

expense of increasing computer memory requirements. quintic B-spline collocation 

can be made more efficient. 

Both methods have the advantage over nonbasis function schemes of allowing 

42 
\ 

" I 

I' 

t 

\ 

, 

! ,. 

t 

. 

,. 

\, 

. ' ,. 

. ~ 

;.. 

l' • ,. 

~ .. 

• 

. 

i· 

I 

I. 

, 

-:,- ---- . 

one to easily obtain derivatives of the solution vector as well as the solution 

at any point in the domain. Their main shortcoming is in the cell Reynolds 

number stability condition (cf., section on numerical results) which manifests 

itself as oscillations in the solution. In view of the above a blanket 

recommendation cannot be given. Only after an analyais of the problem under 

consideration is made and the goals of the solution (accuracy and efficiency) 

have been stated can one recommend either of the schemes. 

Time was not available to progr~ either Rubin's method or the MurphyKendall-Bartl~tt technique and thus a rigorous comparison with the B-spline 

formulation presented here cannot be given. However, it would appear that 

both B-spline procedures are competitive with these other spline variants, 

and they do have the advantage of reducing the order of the blocks for coupled 

systems of equations • 

A comparison of the B-spline algorithms with the Q-R operator techniques 

is given in the section on numerical experiments. 

Application of B-splines to 1-0 Parabolic Equations 

Consider the one-dimensional linear parabolic equation * 

U; • OUxx + bux + cu + d (46) 

with appropriate boundary conditions and initial conditions, where a, b, c, d 

could be functions of x and t. Employing a two point temporal difference 

scheme, we obtain, 

n /3 n+1 n /3 n+1 n OJ + (f3uj' + (I-f3)uJ

' ) + bt+ (f3uj + (I-(3)uj ) 

+ ct+/3(f3 ur+1 + C1-(3)ut> + dt+/3 (47) 

* If Eq. (46) were nonlinear, assume it vas linearized by the method described 

in the appendix. 

43 
.--

I 

I 

f 

~ • I· 

I 

t 

I: 

t 

.. 

U 

L 

~ 

~ 

t~ 

• • 

;.. 

---::1 .a · . 

t:: 

J.1 

• .. 

~ 

t3 

' ... " 

where B • 1 reduces to a backward difference and B • 1/2 to the Crank-Nicolson 

approximation. 

The B-spline spatial approximation can be combined directly with a 

temporal discretization scheme to solve Eq. (46). Here collocation and the 

generalized Galerkin procedures are considered. Although no numerical 

experiments were run with these methods, they are presented here to indicate 

how they can be applied. 

1. Collocation 

Substituting Eqs. (36) and (37) into (47) and rearranging one obtains 

j-I,2,.",J+1 

where 

is known from the previous time step and 

and similarly for 

2. 

b n+f3 n+/3 and d n+/3 J ' Cj j 

Generalized Galerkin (cubic B-splines) 

Integrating Eq. (48) over (j, j + 1) one obtains 

44 

(4S) 
,I 

" : 

By the use of Simpson's rule the integrals can be evaluated in terms of the 

function values at xj • xj+1/

2 and xj

+l • If the values at the "hali points". 

n+B n+B n+B n+B 

aj +l / 2, bj +l / 2• cj +l / 2 and dj+l / 2 are not known analytically they must be 

obtained by an interpolation routine. 

45 

~---.,..---' .-- -, . 
. I 

.. 

,. 

i 

r 

I 

~ 

i 

p ··--~-"'r"-· 

RESULTS 

Numerical Experiments 

Spatial schemes can be analyzed by considering two point boundary value 

ordinary differential equations. The properties of these schemes, accuracy 

versus grid spacing, running time for a given accuracy, ease of programming. 

storage requirements and cell Reynolds number effects can easily be obtained 

without the encumbrance of temporal discretizations. 

Here the following differential equation is considered 

with appropriate boundary conditions specified. 

The first type of problem investigated was the homogeneous constant 

coefficient case 

which has as its exact solution 

u(O) - 0, u(I)-1 

1- e -bx 

Uexact • I - e -6 

(49) 

(50) 

(51) 

The aim of these numerical experiments was to verify the convergence rates 

of the various methods and the effect of cell Reynolds number on obtaining 

physically meaningful solutions. The convergence rates were verified by 

considering the case of b a 1 for 10. 20 and 40 intervals, so that effects of 

cell Reynolds number. which ranged from .025 to .100 \10uld be insignificant. 

In Table V these results are presented for the generalized OCI schemes (OCI-Gl, 

(cf., Table Ill), OCI-G2, (cf., Table II», cubic B-spline collocation, (CBS-COL), 

quintic B-spline collocation, (QBS-COL), and cubic n-spline generalized 

Galerkin, (CnS-GAL). All the schemes except for CBS-COL have fourth order 

convergence rates, with QBS-COL possessing the lowest truncation error. 

46 
I 

I 

i 

I 

I 

I 

I 

I 

," 

/ 

/ 

~~'-~~~~r~:~f~~~;--~'.~1.4~~-~-'--~~--~----~-~ _________ ~,_.~,--------_?,~,,~----~,-~'~"~,~,~'~~4~~~~b*~M~·'~Ai~.P~~'~.~.~~ 

I 

I 

i 

The behavior of QBS-COL is due to the different orders of approximation of the 

derivatives, i.e., 0(h4) for the second derivative and 0(h5) for the first 

derivative, which reduces to 0(h4) for the overall approximation to the 

spatial operator. It is also interesting to note that CBS-COL gives 

identical results as the standard second order centered finite difference 

method, (CFD), since the resultin~ matrices of the two methods are linearly 

related. 

It has been observed (for instance in R~f. 47) that for the solution of 

the boundary layer equations, near the outer edge (at large values of the 

normal coordinate), where the flow is nearly uniform, oscillations and/or 

overshoots in the velocity profile may occur. This nonphysical behavior, 

which can be traced to the violation of a cell Reynolds number stability 

condition could deteriorate the entire solution. Therefore an understanding 

of the behavior of the spatial approximations with respect to cell Reynolds 

number is desirable. Considered here is the constant coefficient equation 

once more with b set equal to 80 for 10, 20 and 40 int~rvals which correspond 

to cell Reynolds numbers of 8, 4 and 2, respectively. In Table VI these 

results are presented. Solution profiles for Rc = 10 are shown in Table VII. 

The generalized OCI schemes give uniformly monotonic solutions (they are 

constructed to do just that), while the B-spline schemes give oscillatory 

solutions for cell Reynolds numbers approximately greater than 2. 

The behavior of the OCI schemes in the range of Reynolds number 2 to 4 

is a property of such an unrestricted (with Rc) sche~e (cf., Ref. 29). As is 

true for all schemes, only in the limit as h ~ 0 docs one obtain the convergence 

rates predicted by the the~ry. Error estimates obtained by Berger et al •• 

(Ref. 29) for the generalized OCI schemes indicate there are several overlapping regions that are dependent on the magnitude of Rc which varies as 

h-P with h fixed. For p positive and large, Rc ~ 0 and one recovers an 0(h4) 

convergence rate. while for.p negative and large, Rc ~ m and one obtains an 

2 O(h ) convergence rate. The transition, which is automatic. result a in a 

second order upwind diff~rencing formula when Rc is large. A more detailed 

discussion for the entire range of p values is presented in Ref. 29. 

By way of comparison, the exponential scheme of Il'in (Ref. 31) has been 

5ho'~ to be uniformly first order accurate (Ref. 48) while the exponential box 

47 
, 

.' 

I 

/ 

, 

j 

i 

/ 

~ ,. 

.. 

z 

:;;. 

51 ." 

;;, / ~ .. 

.:~ ... 

) 

... • 

.' 'J , , 

" _ E 

VQZHf)Q4(C 0: .. 

scheme due to El-Mistikawy and Werle (Ref. 35) is uniformly second order 

accurate (Ref. 36). For the problem under consideration (constant coefficient) 

both exponential schemes give the exact solution as expected. The above 

discussion points out one of the major advantages of the Q-R operator format. 

Once a program has been written within such a framework, various schemes can 

be implemented easily. Even hybrid type schemes that are evaluated pointwise 

are allowed. For example, schemes can be chosen by considerations of accuracy, 

cell Reynolds number behavior, or running time. 

For actual viscous flow problems, where moderate to high cell Reynolds 

numbers appear in uniform regions (ncar the outer edge of the viscous layer), 

the generalized OCI schemes work well in practice (Ref. 49). 

It is important to note that the cubic B-spline generalized Galerkin 

procedure is sensitive to the type of boundary condition set at x = O. When 

collocation at x • 0 was employed, the results were nonphysical while setting 

the first derivative to its exact value at that point gave results comparable 

to QBS-COL. The boundary condition was applied in a region of steep gradients 

so that the second order collocation approximation is not sufficient to 

prevent the erratic behavior. However, at the low cell Reynolds number range, 

the collocation boundary condition worked well as the results in Table V 

indicate. 

The second case considered is the linearized Burgers equation 

(52) 

which has the exact solution 

(53) 

and 

u( -CD) ... I, U{CD) - 0 

The coefficient a was set to 1/2 while v was varied; small v corresponds 

to a shock near x ~ O. Calculations were carried out in the regions 

48 

". 
... ./ 

i 

.:', I.~ .. . ~ 

-5 ~ x ~ 0 and/or 0 ~ y ~ 5, with the boundary cond1tiqns set to their 

exact function values. The Burgers equation mimics a true boundary layer in 

that regions of sharp gradients (near x - 0) corresponds to low cell Reynolds 

numbers while in the "uniform flow" region (I xl .... CD) the cell Reynolds number 

reaches a maximum. 

Results are presented in Table VIII for v - 1/16 and for mesh intervals 

50, 100 and 150 for the following schemes; OCI-GI, OCI-G2, OCI-Swartz, CBSGAL, QBS-COL and Allen's exponential scheme. The QBS-COL method again has 

the lowest truncation error for a given mesh distribution. However, the 

computation times (on a CDC 7600 machine) as presented in Tabl~ VIII indicate 

that in order to attain a given accuracy, OCI-G2 is the most efficient scheme 

for this problem while there is lIttle difference between OCI-Gl and CBS-GAL. 

Although QBS-COL does not fare as well, its performance c~~ b~ imp~c~ed by 

storing the values of the basis functions, Bi , at the nodal pointP an~ not 

computing them as needed. In Table IX OCI-Gl, OCI-G2, CBS-GAL, QBS-COL, the 

Allen exponential and the EI-Mistikawy-Werle exponential schemes are compared 

for the case v = 1/24 as the cell Reynolds number is increased beyond 2. 

Again, the OCI schemes and the exponential schemes have monotonic behavior 

while the B-Spline techniques lead to oscillatory solutions. 

The results indicate that with respect to accuracy there is little to 

choose from among the generalized OCI schemes and the higher order B-spline 

procedures. However, the B-spline procedures lead to larger banded matrices 

and have cell Reynolds number stability restrictions. Note that it is 

certainly not ruled out here that at some point in the future the B-splines 

C9uld be modified to eliminate the oscillatory behavior at high cell Reynolds 

number. Although the B-spline schemes do have the advantage of being able to 

evaluate derivatives and treat derivative boundary conditions more easily than 

the OCI methods, at present they do not seem to be as versatile as OCI at 

least on equation systems that OCI can be applied to. It is therefore 

recommended that for the problem under consideration here for the spatial 

scheme, the family of Q-R operator schemes, in particular the generalized 

OCI schemes be adopted. 

As a model problem for investigating the chosen spatial and temporal 

schemes, a coupled system of one-dimensional, nonlinear, parabolic equations 

49 
· ~·~~~~;i~~~t;~:f.~~ii~.l!l. 

was constructed which exhibit typical boundary layer-like behavior, and have the 

form 

where 

Ut • Uyy + A2Uy + A3U +A4UW + A~W + As 

Wt • Wyy + B2Wy + B3w + B4 uW + B~U + BS 

A~ • -a( E + slnwf), 

and 

-ay 2( . ) A • wcoswt( 1- e ) + a ,,+ Stnwl 6 

B~ • -aosinwf - ye- aY(2y +,Bsinw1), 

{ -CJ.y -ay } B6 .. wcosw1 y,Bye + o(J-e ) + 

- 2ye-ay {1-a(2y+,Bsinwtl} 

with boundary conrlitions 

U(O, t) .. w<O,t) .. 0 

U(a> ,t) .. " + sinw1 " > 1 w(co, t) .. osinw1 8>0 

The initial conditions were chosen to be the exact solution values at the 

initial time. 

The exact solutions are 

-ay U • (E + sinwtH I - e ) 

W • 

50 

(54) 

(55) 

(56) 

(57) 

I 
. . 

-.: 

" 

The exact solutions possess the following properties: 

1. periodicity in time, 

2. exponential decay in space, 

3. nonlinearities and equation coupling, 

and w exhibits a form of "flow reversal". 

The Q-R operator spatial schemes were examined, including second order 

central finite differences, second order Il'in (exponential) and the fourth 

order generalized operator compact implicit method. As noted previously, all 

the above schemes require the inversion of a scalar tridiagonal matrix for a 

single equation in one unknown. 

In the computer code, Eqs. (54) .and (55) were solved simultaneously, 

resulting in a 2 x 2 block tridiagonal system. Both first order fully implicit 

and second order Crank-Nicolson discretizations were ~ncluded in the program. 

Appendix B presents the appropriate linearizations and describes the implementation of the method. 

The functions u and w were chosen, as previously mentioned, to mimic the 

velocity components of a time-dependent boundary layer flow. The IOU velocity 

component" is alw..!ys positive, with its maximum value (which is reached at the 

outer edge) varying between e: + land c - 1. The lOW cross flow velocity 

component", as shown in Fig. I, has a maximum value in the interior, and 

exhibits a cross flow reversal in that w changes sign through a time cycle. 

The numerical results to be discussed are for the case shown, i.e., a a 1, 

a a 4.0, y a -2.0, 0 a 0.40, and e: a 1.5. 

A comparison of Newton iteration with the second order noniterative time 

linearization procedure wa~ made. The results shown in Figs. 2 and 3 verify the 

contention that the effort spent at iterating can be more effectively used in 

decreasing the step size, thus reducing both the temporal truncation and 

linearization errors. The spatial step size for the calculation shown in 

Figs. 2 and 3 was chosen (at h = .2), and the method employed was OCI-GI, so 

that the predominant error would be due to temporal effects, e.g., linearization· 

and discretization. As an auxiliary benefit, time linearization requires one 

less time level of storage, and is easier to program (cf., Appendix B). 

Figures 4 and 5 present a comparison of the generalized OCI scheme with 

second order central differences. The benefits in accuracy of the higher order 

OCI scheme. is evident in these figures. In fact, for the J c 40 case, the 

51 
'.. 

-4 temporal truncation error, which can be shown to be of order 10 , dominated 

the spatial truncation error so that the OCl scheme did not attain its 

theoretical convergence rate. 

52 
. , 

..... -. . , 

.' 

CONCLUSIONS 

In this report several temporal schemes have been investigated. consistently split block icplJ.cit. Rubin's predictor-corrector and liopscotch. 

It has been shown that. in order to meet tho requirements of three-dimensional 

unsteady boundary layer flows. the consistently split linearized block implicit 

scheme 1s the most versatile and efficient of the three. and was thus 

recommended. 

Tvo approaches to spatial approxication were described. the Q-R operator 

formulution and the B-spllne basis function technique. Results of numerical 

experim!nts indicate that both quintic B-spline collocation and generalized 

OCI schcmes performed well and that the Q-R operator fort:lUlation, in particular 

the generalized OCI scheces, arc at present particularly well suited for the 

problem of time-depcndent boundary layers in regard to efficiency, cell 

Reynolds number stability restrictions, and flexibility. 

The linearized block implicit temporal scheme, in conjunction with a 

generalized OCI spatial scheme, was employed to solve a system of two coupled 

nonlinear parabolic equations that exhibit "three-dimensional unsteady boundary 

layer" behavior. The results of n\lccrical experitllents indicate that the 

generalized ocr approach 1s viable and can be applied to viscous flow problecs. 

53 
APPENDIX A 

Linearization Techniquc 

A number of techniques have been used for implicit solution of the 

following first-order nonlinear scalar equation in one dependent variable 

~(x,t): 

Special cases of Eq. (AI) include the conservation fo~ if F(~) • 1, and 

quasi-linear flow if G(~) .~. Previous implicit methoda for Eq. (AI) 

(AI) 

which employ nonlinear difference equationo and aloo methods based on twostep predictor-corrector schemes arc discussed by Ames (Ref. 50, p. 82) and 

von Rosenburg (Ref. 51), p. 56). One such method is to difference nonlinear 

terms directly at the impliCit time level to obtain nonlinear implicit 

difference equations; these are then solved iteratively by a procedure such 

as Newton's method. Although otherwise attractive, there may be difficulty 

with convergence in the iterative solution of th~ nonlinear difference 

equations, and some efficiency is sacrificed by the need for iteration. An 

implicit predictor-corrector te~hnique has been dev~sed by Douglas and Jones 

(Ref. 52) which is applicable to the quasI linear Cdse (G •• ) o[ Eq. (AI). 

The first step of their procedure is to lineari1.c the equation by evaluatin& 

n n+l/2 the nonlinear coefficient as F(¢ ) and to predict values of • using either 

n+l the backward difference or the Crank-Nicolson scheme. Values [or ~ are 

then computed in a similar.manner using F(¢n+l/2) and the Crank-Nicolson scheme. 

Gourlay and Horris (Ref. 53) have also proposed implicit predictor-corrector 

techniques which cnn be applied to Eq. (AI). In the conservative case (F - I), A A 

their technique is to define G(¢) by the relation Get) • ?G(.) when such a A n+l n+l definition exists, nnd to evaluate GCt ) using values for ~ computed by 

an explicit predictor scheme. With G thereby known at the implicit time level, 

h i b d Ii d d I f n+1 d t e equat on can e treate as near an correcte va ues o. are compute 

by the Crank-Nicolson 6che~e. 

A technique is described here for deriving linear implicit difference 

approximations for nonlinear differential equations. TIle technique is based 

51, 

. 

• 
on an expansion of nonlinear implicit teres about the solution at the known 

time level, tn, and leads to a one-step, two-level scheme which, being linear 

in unknown (implicit) quantitics, can be solved cfficiently without iteration. 

This idea was applied by Richtmyer and Horton (Ref. 16, p. 203) to a scalar 

nonlinear diffusion equation. Here, the technique is develop2d for problems 

governed by t nonlinear equations in t dependent variables which are functions 

of time and space coordinates. The technique will be described for the threedimensional, unsteady equations. 

The solution domain is discretized by 

1 2 in the comp~tational coordinates, 6y , hy 

grid points having equal spacings 

3 I 2 3 and hy in the y ,y and y 

directions, respectively, and an arbitrary time step, ht. 

k and superscript n are grid point indices associated with 

The subscripts i, j, 

:,1 2 3 d : , Y ,y nn t, 

n I 2 3 n respectively, and thus ~i,j,k denotes ~(Yi' Yj' Yk' t). It is assumed that 

the solution is known at the n level, tn, and is desired at the (n+1) level, 

n+l t • At the risk of an occasional ambiguity, one or more of the subscripts 

is frequently omitted, so that ¢n is equivalent to ¢~,j,k' 

The numerical method employed is quite general and is formally derived for 

systems of governing equations which have the following fore: 

where ¢ is a column vector containing t dependent variables, It and 5 are 

column vector functions of 9, and ~ is a column vector whose elements are 

(A2) 

spatial differential operators which may be multidimensional. The generality 

of Eq. (A2) allows the method to be developed conCisely and permits various 

extensions and modifications (e.g., noncartesian coordinate systems, turbulence 

modelR) to be made more or less routinely. It should be emphasized, however, 

that the Jacobian all/a¢ must usually be nonsingular if the ADI techniques as 

applied to Eq. (A2) are to be valid. A necessary condition is that each 

dependent variable appear in one or more of the governing equations as a time 

derivative. An exception would occur if for instance, a variable having no 

time derivative also appeared in only one equation, so that this equation could 

be decouplcd from the remaining equations and solved ~ posteriori by an alternate method. 

55 

• I 
The linearized difference approximation is derived from the following 

implicit time-difference replacement of Eq. (A2): 

. 

(Hn+l-H")/6t =.8[.2> (4)"+') +5"+'] +(.- .Bl[~ (cpn)+ SO] (A3) 

where, for example, Hn+l : H(Qn+l). The form of ~ and the spatial differencing arc as yet unspecified. A parameter B(O ~ B ~ 1) has been introduced 

so as to permit a variable centering of the scheme in time. Equation (A3) 

produces a backward difference formulation for B • 1 and a Crank-Nicolson 

formulation for B a 1/2. 

The linearization is performed by a two-step process of expansion about 

the known time level t

n and subsequent approximation of the quantity 

( I n n+l . n il¢ 3t) lit, which arises from chain rule differentiation, by (~ - ¢ ). The 

result is 

(A4a) 

(Mb) 

(Me) 

The matrices ilH/il¢ and aS/3¢ are standard Jacobians whose elements are defined, 

for example, by (aU/3¢) = 3H 13¢. TIle operator elements of the matrix 

qr q r 

3:D13¢ arc similarly ordered, i.e., (3~/3¢) = 3~ 13¢ ; however, the 

qr q r 

intended meaning of the operator elements requires some clarification. For 

the qth row, the operation (3:D la?)n(?n+l - ¢n) is understood to mean that ~ n q n {a/3t-u [¢(x,y,z.t»)l 6t 15 computed and that all occurrences of (a¢r/3t) 

q n+l n arising from chain rule differentiation arc replaced by (¢ - ¢ )/llt. r r 

After linearization as in Eqs. (M). Eq. (AJ) becomes the following linear 

implicit time-differenced scheme: 

56 

... 

iii 
, 

, 

y---- .. I 

• 

,d 

Si 

pi 

, . 

,.: 

.... 

~ 

-

p 

U 

, 

, 

n+1 Although H is linearized to second order in Eq. (A4). the division by At 

in Eq. (Al) introduces an error term of order At. A technique for maintaining 

formal second-order accuracy in the presence of nonlinear tiDe derivatives is 

discussed by McDonald and Briley (Ref. 7). however. a three-level scheme 

results. Second-order temporal accuracy can also be obtained (for 8 • 1/2) by 

a change in dependent variable to ~ = H(~). provided this is convenient. since 

the nonlinear tiDe derivative is then eliminated. The temporal accuracy 

is independent of the spatial accuracy. 

On examination. it can be seen that Eq. (A5) is linear in the quantity 

n·H n (~ - ~ ) and that all other quantities are either known or evaluated at 

the n level. Computationally. it is convenient to solve Eq. (AS) for 

n+l n n+l (~ - ~ ) rather than ~ • This both simplifies Eq. (AS) and reduces 

roundoff errors, since it is presumably better to compute a small 0(6t) change 

in an 0(1) quantity than the quantity itself. To simplify the notation. a 

new dependent variable ~ defined by 

is introduced, and thus ~n+l 9n+l _ ~n. and ~n • O. It is also convenient 

to rewrite Eq. (AS) in the following simplified form: 

(A+At .!)lJI n+1 = At [~(c/-In)+Sn] (A7a) 

where the following symbols have been introduced to simplify the notation: 

(A7b) 

.I~ - /3(iJ ~ liJc/-I) (A7c) 

It is noted that .1(~) is a linear transformation and thus .1(0) • O. Furthermore if .1(9) is linear. ~hen 1(~) - -13))(~). 

57 
• 

I, 

" 

.. 

Spatial differencing of Eq. (A7a) is accomplished simply by replacing 

d i er vat i ve operators suc h as "/"yi. Q ,,2/">yi"yi by di fi i 0 Q 0 0 correspon ng n te 

difference operators. 0i' O~. Henceforth. it is assumed that ~ and ~ have 

been discretized in this manner. unless otherwise noted. 

Before proceeding. some general observations seem appropriate. The 

foregoing linearization technique assumes only Taylor expandability. an assumption already implicit in the use of a finite difference method. The governing 

equations and boundary conditions are addressed directly as a system of coupled 

nonlinear equations which collectively determine the solution. The approach 

thus seems more natural than that of making ad hoc linearization and decoupling 

approxim.~tlons. as is often done in applying implicit schemes to coupled 

and/or nonlinear partial differential equations. With the present approach. 

it is not necessary to associate each governing equation and boundary condition 

with a particular dependent variable and then to identify various "nonlinear 

coefficients" and "coupling terms" which must then be treated by lagging. 

predictor-corrector techniques. or iteration. The Taylor expansion procedure 

is analogous to that used in the generalized Newton-Raphson or quasilinearization methods for iterative solution of nonlinear systems by expansion 

about a known current guess at the solution (e.g., Bellman & Kalaba, Ref. 54). 

However, the concept of expanding about the previous time level apparently 

had not been employed to produce a noniterative implicit time-dependent sche~e 

for coupled equations, wherein nonlinear terms are approximated to a level of 

accuracy commensurate with that of the time differencing. The linearization 

technique also permits the implicit treatment of coupled nonlinear boundary 

conditions, such as stagnation pressure and enthalpy at subsonic inlet 

boundaries, and in practica, this latter feature was found to be crucial to 

the stability of the overall method (Ref. 17). 

58 
• • 

, 

/ 

/ 

/ 

i, 

APPENDIX D 

I 

.• -~ .• -" - -. ""':'-:."'..,..---...., .. "'-1" ••. , 

Discretization of Model Equations 

P 5 

Herein we shall consider the system of model equations discussed in the 

section on numerical experimento namely, 

where 

and 

where 

A~ • -a( E + sinwt), 

B .. • a, 

B" • -a8sinwt - ye- aY(2y + ,Bsinwt), 

86 • wcos~t { y,Bye-ay + 8C1-e-ay>} "+ a2 8sinwt 

- 2ye -ay { 1- a( 2y + ,Bslnwn} 

Since nonlinear terms wu ,uw and uw appear in Eqs. (Bl) and (B2). y y 

(Dl) 

(B2) 

linearization is required. Two types of linearization procedures are described, 

time linearization and Newton iteration (quasilinearization). The application 

of these methods within the Q-R operator framework 1s also demonstrated. 

59 

·-~I 
.-~-.~-

.. 

• 

I 

, I 

~-

~ 

I ; 

i 

! 

I· 

/ 

/ 

.j 

/ 

. 

.' 

I 

• 

t 

~ ... 

~ 

;.. 

L·" . ~' 

L. 

.., 

-. 

· · I 

f ,-

r: 

!..~ 

, 

: I -~' 

· . iA 

fl 

~1 

~ ... 

r.; 

, i 

... 

C 

4 A "¥4 ' • 

Time Linearization 

I 

! 

The nonlinear terms uw and wu are linearized by the method described in y Appendix A to yield 

where B • 1/2 corresponds to a Crank-Nicolson scheme and B = 1 to a fully 

implicit method. 

Substituting Eqs. -(B3) an'd (B4) into (Bl) we obtain 

O-'R 

(u"+I- U")/6t • 6yZ [{3U"+1 + (1-{3)u"] 

+ A~ [{3wnu

n

+

1 + {3unwn+1 + (1- 2{3)u"wn] 

+ A6[ {3W"~1 + CI-f3>w"] + A7 + {3u;wn+1 - {3wnu

y

" 

which, with some rearrangement, reduces to 

-I [I + AyZ( A3 (1-2P) ")a cr -p) j n o >..P u + P A~w ~+ -pR U 

- w n Uy "j 

where 

and Q and R are evaluated at the nth time level. 

60 

(B3) 

(B4) 

(BS) 
, 

/ 

/. 

/ 

. 

/ 

.i : 

/ ; I 

.' ' 

-_. __ ..,._ --6 

. __ -

/ 

r""'"" ." A-'. 4 ~ " ,,4 • 4 . • ,," .'c e.c ... , .. 442 '.Ie;;; ". * 'iI44P!lJ';'Q:@'l'JCl:'M;ue.A';;;PSii' 

Ii Similarly, Eq. (B2) reduces to 

! ' 

(B6) 

+ lly2Q _ B un + _7 I (1-,8) B 

/3 6 /3 - u n n Wy I 

Equations (B5) and (B6) arc solved as a coupled (2 x 2 block) system for u n+l 

n+l and w • 

Newton Iteration 

Alternatively one can linearize about a previous iteration instead of the 

previous n+B time level, to yield for the nonlinear term (wu ) for example, 

y 

(WU y)n+/3 .. [ [3'l.w* + [3C1-(3)wn] U;+I + [[32 U

y

* + 13C1-I3)u;]wn 

+ [cr-J3)2wnu; - J32 w*u

y

*] (B7) 

Substituting Eq. (B7) in Eq. (Bl) the quasilinear approximation becomes 

(Q-IR) 

CU

M1 - un)ll~t" lly2 [f3un+1 + CI-I3)un] 

+ A~[(f3W*+(I-f3)wn)f3un+1 + {f3u*+ {1-I3)un)l3wn+1 

+ IJ 2 n n f32 * * [IJ n+1 { IJ nJ (/-,..,) u w - u W +A

6 ,..,w + I-,..,)w 

+ A7 + [f3U y* +CI-(3)Uyn]f3wn+1 - [l3u; + (1-f3)u;]f3w* 

61 
, 

I 

/ 

/ 

£.~;t#!'J~~ ! --..-. . .... ,. _ ' .......... 1< .• 'T • J 

which after some manipulation reduces to 

{o[ )..~ - Lly2(AlI(,8W*+CI-,8)w"»]-R}U"+1 

Lly20 {AlI(,8U* + (1-,8)u") + As + (,8uy* + (1-,8)U;)} w"+1 

• {a [ )"/3 + Lly2 0-;) (AlIO-,8)W")] + (I~/3) R }u" 

(BB) 

{

A (1-/3) * } + Lly20 ; + -/3-A

s

W" - ,8AlIu*w - (,8u; +CI-,8}u;)w* 

Here Q and R are assumed to be evaluated at (n+a). A similar expression 

can be obtained for Eq. (B2), viz., 

{o[ :/3 - Lly2( 811(,8u~t - (I-,8)u"»] - R}W"+I 

-6y20 {SlIC,BW*+C'-,B)w") +8s +C,BWy*+C'-,8>w~)}U"+1 

.10[-'- + Lly2 (I-/3)C8 (I-,B)u")] + CI-/3) R}W" ~)./3 {3 1I /3 

{ A (/-/3) } + Lly20 f37 + p 8sU"-,8BlIU*W*- (,8Wy* + (1-,8)w;)u* 

Equations (BB) and (B9) are solved as a coupled system for u * * and w , the 

(B9) 

latest iterants. The coefficients and the Q and R operators are updated, and 

the system of equations is solved again. The process is continued until convergence is attained. The process involves: 

1. updating the coefficients A2 - A7 and B2 - B7, 

2. updating the Q and R operators, 

3. solving the resulting 2 x 2 block system, 

4. repeating steps 1-3 until convergence is attained. 

62 
" 

i 

~ ~h'i'~~J;r":'\~'~ Y.p-·f.~:'.~!j~J.¥.,h~~~~~;;.¥,~.~ :q~~;;~;i;;;~t;~~~A 

. t 

Note that Eqs. (BS), (B6), (BB) and (B9) can easily be transformed to increment 

form. However, the results in this report were obtained for the form of the 

equations given above. 

63 
\. 

I 

t 

f 

.-" 

~" 

t 'I 

t 

t 

'. 

,-" 

b-

~. 

, 

I • 

~ 

I 

, 

I 

l 

; 

r 

~ , 

REFERENCES 

1. Kendall, R. M., Bonnett, l~. S., Nardo, C. T., and Abhett, M. J.: Computational Procedure for Three Dimensional Boundary Layers on Aircraft and 

Aerospace Vehicles. Proceedings AI~~ 2nd Fluid Dynamics Conference, 

Hartford, 1975, pp. 113-121. 

2. McCroskey, W. J.: Some Current Research in Unsteady Fluid Dynamics. 

J. Fluids Engineering, Vol. 99, March 1977, pp. 8-38. 

J. Nash, J. F., and Patel, V. C.: Three Dimensional Turbulent Boundary 

Layers. Scientific and Business Consultants Inc., 1972. 

4. Krause, E.: Comment on Solution of a Three-Dimensional Boundary Layer 

with Separation, AlAA J., Vol. 7, No. 3, ~rch, 1969. 

5. Lin, T. C., and Rubin, S. G.: Viscous ~low over a Cone at Moderate 

Incidence. Part 2. Supersonic Boundary Layer., J. Fluid Mech., Vol. 59, 

Part 3, 1973, pp. 593-620. 

6. Briley, W. R.: Numerical Method for Predicting Three-Dimensional Steady 

Viscous Flow in Ducts. J. Compo Physics, Vol. 14, ~lo. I, 1974. 

7. HcDonald, It., and Briley, W. R.: Three-Dimensional Supersonic Flow of a 

Viscous or Inviscid Gas. J. of Compo Physics, Vol. 19, No.2, 1975. 

8. Blottner, F. G.: Investigation of Some Finite Difference Techniques for 

Solving the Boundary Layer Equations. Computer l1ethods in Applied 

Mechanics and Engineering. Vol. 6, 1975, pp. 1-30. 

9. Thompson, B. G. J., HacDonald, A. G. J.: The Prediction of Boundary Layer 

Profile Development for Infinite Yawed l~ings. Aero Res. Council, CP 1307. 

10. Briley, W. R., and McDonald, H.: On the Structure and Use of Linearized 

Block ADI and Related Schemes. Scientific Research Associates, Inc., 

Report No. 78-3. To be published in the Journal of Computational Physics. 

11. Peaceman, D. W., and Rachford, H. H.: 

and Elliptic Differential Equations. 

Vol. 3, 1955, pp. 28-41. 

The Numerical Solution of Parabolic 

Soc. for Indust. Appl. Math., 

12. Douglaa, J.: On the Numerical Integration of uxx + Uyy c ut by Implicit 

Methods. Soc. for Indust. Appl. Math., Vol. 3, 1955, pp. 42-65. 

13. Hitchell. A. R.: Computational Methods in Partial Differential Equations. 

Wiley, New York, 1969. 

14. Yanenko, N. N.: The Method of Fractional Steps. Springer-Verlag, 1971. 

15. Douglas, J. and Gunn, J. E.: A General Formulation of Alternating Direction Hethods. Numerische Math., Vol. 6, 1964, p. 428. 

64 
• 

& 

~ 

t 

! 

" i 

I' 

~ Q 

~ 

r 

·11 

~ il 

II 

" 11 

• /~I 1 

Ii 

-j! 

• 1 

I; 

i 

I 

: , 

i: 

i. 

t 

~ [ 

f; 

F.; 

t.: 

i:i 

~ 

~I 

f.i 

~ rl 

to; 

II 

16. Richtmyer, R. D., and Morton, K. W.: Difference Methods for Initial Value 

Problems. Second Edition. Interscience Publisher~, New York, New York, 

1967. 

17. Briley, W. R., and McDonald, H.: Solution of the Multidimensjonal 

Compressible Navicr-Stokes Equations by a Generalized Implicit Method • 

J. Compo Physics, Vol. 24, No.4, August, 1977, p. 372. 

18. Gibeling, H. J., McDonald, H., and Brilpv, W. R.: Development of a 

Three-Dimensional Combustor Flow Anal.,,,i:··. Volume I: Theoretical 

Studies. Technical Report AFAPL-TR-7'· 5:'. Vol. I, July 1975. 

19. Gibeling, H. J., McDonald, 11., and Br.i.ley, W. R.: Development of a 

Three-Dimensional Combustor Flow Analysis, . Jlume II: Theoretical 

Studies. Technical Report AFAPL-TR-75-59, Vol. II, October, 1976. 

20. Briley, H. R.: Priv.:1te communication. 

21. Rubin, S. G., Lin, T. C.: A Nu~erical Uethod for Three Dimensional 

Viscous Flow: Applic.:1tion to the Hypersonic Leading Edge. J. of Compo 

Physics, Vol. 9, 1972, p. 339. 

22. Gourlay, A. R.: 

Equation Solver. 

Hopscotch: A Fast Second-Order Parti.:1l Differential 

J. Inst. Maths. Applics. No.6, 1970, pp. 375-390. 

23. Gourlay, A. R., and McGuire, G. R.: General Hopscotch Algorithm for the 

Numerical Solution of Partial Differential Equations. J. Inst. ~mths. 

Applics. No.7, 1971, pp. 216-227. 

24. Gourlay, A. R.: Some Recent Methods for the Numerical Solution of TimeDependent Partial Differential Equations. Proc. Royal Soc. Lond. A 323, 

1971, pp. 219-233. 

25. Smith, R. E., Jr., and Kidd, A.: Comparative Study of Two Nurn~rical 

Te.chniques for the Solution of Viscous Flow in a Driven Cavity. NASA 

SP-378, 1975 • 

26. Buzbee, B. L., Golub, ·G. H., 

Solving Poisson's Equations. 

Dec., 1970, pp. 627-656. 

and Nielson, C. W.: On Direct Methods for 

SIAM J. Numerical Anal., Vol. 7, No.4, 

27. Greenberg, J. B.: Hopscotch Algorithms - A Study of Some Hybrid ..Einite':'.· 

Difference Techniques. TAE Report No. 270, Technion-Israel Institute of 

Technology, Haifa, Israel, 1976. 

28. Ciment, ~1., Leventhal, S. H., and Weinberg, B. C.: The Operator Compact 

Implicit }lethod for Parabolic Equations. J. Cornp. Physics, 28 (1978), 

pp. 135-166. 

29. Berger, A. E., Solomon, J. M., Ciment, M., Leventhal, S. It., and 

Weinberg, B. C.: Generalized Operator Compact Implicit Schemes for 

Boundary Layer Problems. Submitted for publication. 

65 
30. Allen, D. N., and Southuell, R. V.: Relaxation Methods Applied to Determine the Motion, in Two Dimensions, of a Viscous Fluid Past a Fixed 

Cylinder, Quart. J. Mech. Appl. Math., 8, 1955, pp. 129-145. 

31. I1'in, A. M.: Differencing Scheme for a Differential Equation uith a 

Small Parameter Affecting the Highest Derivative. Mat. Zametki, 6, 1969, 

pp. 237-248 - Math. Notes, 6, 1969, pp. 596-602. 

32. McDonald, H.: A Novel Finite Difference Scheme Particularly Suited for 

Fluid F10u Problems Containing Discontinuities. Scientific Research 

Associates Report P76-l, August 1976. 

33. Dennis, S. C. R.: Finite Differences Associated with Second-Order 

Differential Equations, Quart. J. Mech. Appl. Math., Vol. 23, 1960, 

pp. 487-507. 

34. El-Mistikawy, T. M., and Werle, M. J.: Numerical Methods for Inviscid, 

Viscous Fluid Flous, Report No. AFL 77-9-34, University of Cincinnati, 

March 1977. 

35. El-Histikawy, T. M., and Werle, M. J.: Numerical Method for Boundary 

Layers uith Blowing - the Exponential Box Scheme. AIM J., 16, 1978, 

pp. 749-751. 

36. Berger, A. E., Solomon, J. M., and Ciment, M.: An Analysis of a Uniformly 

Accurate Finite Difference Method for a Singular Perturbation Problem. 

Paper presented at Second International Conference on Computational Methods 

in Nonlinear Hechanics. Harch 26-30, Austin, Texas. 

37. Prenter, P. M.: Splines and Variational Methods. John Wiley & Sons, 

New York, 1975. 

38. Kendall, R. H., and Bartlett, E. P.: Nonsimilar Solution of the MultiComponent Laminar Boundary Layer by an Integral-Matrix Hethod. AIM J. 

6, 1968, p. 1089. 

39. Murphy, J. D., Presley, L. L., and Rose, H. C.: On the Calculation of 

Supersonic Separating and Attaching Flows, NASA SP-347, 1975, pp. 151-175. 

40. Ahlberg, J. H., Nilson, E. N., and \o]a1sh, J. 1.: The Theory of Splines am 

Their Applications. Academic Press, New York, 1967. 

41. Rubin, S. G., and Khosla, P. K.: Polynomial.Interpolation Methods for 

Viscous ~low Calculations, J. Compo Physics, Vol. 14, No.3, 1977, 

pp. 217-244. 

42. Rubin, S. G., and Khosla. P. K.: Turbulent Boundary Layer tnth and Without 

Mass Injection. Computers and Fluids, Vol. 5, 1977, pp. 241-259. 

43. DeBoor, C.: On Calculating uith B-Splines. J. Approximation Theory. 6, 

1972, pp. 50-62. 

66 
. -. 

~. 

!-~ . 

. .; 

:1-- - ~ , 

:,~ 

.".:::~"4 ... _J'~ 

:.aI~~ 

-. :;~r-

,. ---.. 

.:~~;; 

-', 

44. r,eDoor, C.: Package for Calculating with B-Splines. SIAM J. Numer. Anal •• 

Vol. 14, No.3, June 1977. 

45. Crandall, S. II.: Engineering Analysin. McGraw-IIUl Book Coc:pany, 1956 • 

46. Murphy, J. D.: Application of the Generalized Galerkin Method to the 

Co~putation of FluId Flows. Proc. AlAA Co~putational Fluid Dynamics 

Conference, 1973. 

47. Weinberg, B. C., L~ventha1, S. n •• Ciment, M.: The Operator Compact 

Implicit Scheme [or Vincous Flow Problems. AlAA Paper No. 77-638. 

Presented at 3rd AlAA Compo Fluid Dynamics Conference, Albuquerque NH, 

1977 • 

48. Kellogg, R. D., and Tsan, A.: Analysis of Some Difference Approximations 

for a Singular Perturbation Problem without Turning Points. To appear in 

Math. Co:p. 

49. Weinberg, B. C.,'and Glowacki, W.: lIisher Order Methods for Threenimcnsional Boundary Layers. Paper in preparation. 

50. Ames, H. F.: Numerical Hethods for Partial Differential Equations. 

Barnes & 90blc, Inc., New York, New York, 1969. 

51. von Roaenberg, D. A.: Hethods (or the NUUlerical Solution of Partial 

Di((erential F.quation5. American Elsevier Publishing Co., Inc., New York. 

New York, 1969. 

52. Douglas, J., and Jones. D. F.: On Predictor-Corrector Methods for Nonlinear Parabolic Differential Equations. Soc. for Indust. Appl. Hath., 

Vol. II, 1963, pp. 195-204. 

53. Gourlay, A. R., nnd trorris, J. L.: Finite-Difference Methods for Nonlinear 

Hyperbolic Systems. Hath. Comp., Vol. 22, 1968. pp. 28-39. 

54. Dell~an, R. E., nnd Kalaba. R. t.: Quasilinearization and Nonlinear 

Boundarv-\'al\le Problems. A!uerican Elsevier Publ. Co., Inc •• New York, 

1965. 

67 
., . 

... 

/ 

.---l 

,I 

to 

\.Ihcre 

TABLE 1. - OPERATOR COEFFICIEtITS FOR STANDARD 

ilPERATOR COMPACT IMPLICIT SCHEME 

68 

\ 

\ 

". 
~ . 

- ... -._._-

where 

TABLE II. - OPERATOR COEFFICIENTS FOR GENERALIZED 

OPERATOR COMPACT IMPLICIT SCHEME OCI-C2 

qr· 6+[P,-3]RCJ + [pZ] RCJ2 

q: · 60 +[IOP,]RCJ + [P:I]RC~ + [T"J+,P4]RC: 

qt· 6 +[p,+3]RCJ + [P,+PZ]Rc:j2 + [P4 ] RCj3 

7T, • CT"j+1 + T",-,)pz + T",+,p, + lTl 7TZ • 15 - 2pz + CCTz-llpz-3(T"j+I+ CTz) + TTZ 

7f1 .\ ~ 0", ~Ol 7f

z .\ 0 2P'-0"2~O I 

8'CT,2(10 - T"J+I- T"J-I) 0", ~O (2p,- 0"2)2/8 2p,-CTZ <O 

(T"I_' - T"lt') CTI • PI + 3 + ---''----''-

10 - T"j+l- T"j_1 

with h sufficiently small so that 

··here 

69 

" 
V·· I 

..... 

.. 

. 

I 

.. 

! 

~ 

~ ,j 

; 

• 

~ " i 

TABLE Ill. - OPERATOR COEFFICIENTS FOR GENERALIZED 

OPERATOR COMPACT UIPLICtT SCHEME OCI-Gl 

QJ-'C,+ • a -,c,+ + o-,c,+ + a -,c,+ + a -,c,+ o , 2 :5 . 

o~ .. 60 • o~ • o~ • 6 

C _ I ( OC C 0+ + I 0- -) O2 • W 2 - 202 -"6 - 2 ,p + ,p - I ,P 

0 + • 0- + _1_(OCpC - 20-p- - 20+p·) 2 2 6.' , , 

70 
T.J. .• 

--- where 

TABLE III. - CONTINUED 

w s • 0 

I I ( c c - - + +>/[ c I (. ->] }Lit • -2-40-'g=-p"="C OI.P - 0 1 P - 0 1 P P - 10 p + P 

+ 4(p· + p- + 2/)Q, - ogl 

+ t[(-20~l+otp+-1I01-P-> 

with h sufficiently small so that 

rj • rj. r; coefficients are given in TABLE I 

71 

-
TABLE IV. - OPERATOR COEFFICIENTS FOR EL-MISTlKAWY WERLE 

EXPONENTIAL BOX SCliEME 

where 

and 

rt · p+/[I-eXp(-p+>] 

rt • -(rt + rj-> 

72 

,I 

/ 

i 

~ 
TABLE V. - COMPARISON OF CONVERGENCE RATES FOR 

CONSTANT COEFFICIENT PROBLEM 

Scheme Jmax R Max error Rate 

c 

10 .100 -.576 E-07 

3.99 

OCI-G1 20 .050 -.363 E-08 

3.99 

40 .025 -.229 E-09 

10 .100 -.105 E-06 

3.93 

OCI-G2 20 .050 -.687 E-08 

3.91 

40 .025 -.438 E-09 

10 .100 .101 E-03 2.01 

CBS-COL 20 .050 .251 E-04 2.00 

40 .025 .629 E-05 

10 .100 -.163 E-07 3.97 

QBS-COL 20 .050 -.10~ E-08 3.99 

40 .025 -.656 E-10 

10 .100 -.296 E-05 3.97 

CDS-GAL 20 .050 -.189 E-06 3.99 

40 .025 -.119 E-07 

73 

;' 

II-

-II> 
I -LI' ... : .... I 

/ 

\ 

'/ 

CII' 

TABLE VI. - COMPARISON OF SPATIAL ERRORS FOR 

CONSTANT COEFFICIENT PROBLEM 

Scheme Jrnax Rc Max error Conunent 

CBS-GALl 40 2 -.158 E+Ol Oscilla tory 

QBS-COL 40 2 -.484 E-03 Oscillatory 

OCI-Gl 40 2 -.162 E-Ol Monotone 

OCI-G2 40 2 -.863 E-02 Monotone 

CBS-GALl 20 4 -.145 E+Ol Oscillatory 

QBS-COL 20 4 .612 E-Ol Oscillatory 

OCI-G1 20 4 -.586 E-Ol Monotone 

OCI-G2 20 4 -.127 E-Ol Monotone 

CBS-GAL 10 8 -.136 E+Ol Oscillatory 

QBS-COL 10 8 .288 E+02 Oscillatory 

OCI-Gl 10 8 -.561 E-Ol Monotone 

OCI-G2 10 8 -.425 E-02 Monotone 

lCollocation at x=O used. 

74 

• - 1 ~."....---"""""-- ,. 
.. 

·. / , I 

/ ,--::::~ ./ 

,r' 

! 

/ 

I 

/ 

--.-... ~-. -..-.. _ .... - .. '~.,..-.'--~--.,--/.-.~.,..~ 

TABLE VII. - SOLUTION PROFILES FOR CONSTANT 

COEFFICIENT PROBLEM - Rc • 10 

X Exact Solution Computed Solution 

(QBS-COL) (CBS-GAL) 1 

0.0 o. O. O. 

0.1 .999955 E+OO ,.374802 E+OO .118628 E+01 

0.2 .999999 E+OO .304662 E+01 .965301 E+OO 

0.3 .100000 E+01 .311375 E-01 .100646 E+01 

0.4 .100000 E+01 .251958 E+01 .998796 E+OO 

0.5 .100000 E+01 .476597 E+00 .100022 E+01 

0.6 .100030 E+01 .215299 E+01 .999958 E+OO 

0.7 .100000 E+Ol .778851 E+OO .100001 E+01 

0.8 .100000 E+01 .189431 E+01 .999999 E+OO 

0.9 .100000 E+01 .107440 E+01 .100000 E+01 

1.0 .100000 E+Ol .100000 E+01 .100000 E+Ol 

1 ' Exact derivative at x = 0 set • 

75 

(OCI-G2) 

O. 

.997532 E+OO 

.999994 E+OO 

.999999 E+OO 

.100000 E+01 

.100000 E+01 

.100000 E+01 

.100000 E+01 

.100000 E+01 

.100000 E+01 

.100000 E+01 

- 4.z . , 

• 
/f 

f 

• 

../ 

.. 

~ , . 

. I 

~ 

.f 

... 

~ 

'1 

~ .. 

~ 

.. 

e 

TABLE VIII. - COMPARISON OF RUNNING TIMES FOR 

BURGERS EQUATION - v Q 1/16 

Scheme Jmax Max error Running time1 

(st!c) 

50 .336 E-02 .016 

Allen 100 .828 E-03 .028 

Exponential 

150 .367 E-03 .043 

50 .260 E-03 .032 

CllS-GAL 100 .155 E-04 .042 

150 .297 E-03 .053 

50 .404 E-04 .078 

QBS-COL 100 -.346 E-05 .118 

150 -.679 E-06 .157 

50 .242 E-03 .025 

OCI-G1 100 .150 E-04 .044 

150 .298 E-05 .064 

50 -.696 E-04 .016 

OCI-G2 100 -.656 E-05 .031 

150 -.148 E-05 .044 

1 CPU time on C~C 7600 - excludes input/output. 

76 
, 

/ 

/ 

.. 

i 

, 

f 

r 

• 

,. 

., 

" . , 

.j 

. ;~ 

1 

,. 

" co, . ¥ , 

., 

r 

I 

¢ " ... us p 8' If.' , (.- II;c.s:;z;;;:;:;;mfiW*.:;;c' F. #1*( . :~ ..r' .. c:;; .. e ; ," Gte·' ' .. I i ] V t 

Scheme 

Allen 

El-Mistikawy 

Werle 

CBS-GAL 

QBS-COL 

OCI-G1 

OCI-G2 

Allen 

E1-MistikaW)' 

Werle 

CBS-GAL 

QBS-COL 

OCI-Gl 

OCI-G2 

TABLE IX. - COMPARISON OF SPATIAL ERRORS FOR 

BURGERS EQUATION - v • 1/24 

Jmax Rc max Max error 

25 2.4 .182 E-Ol 

25 2.4 -.213 E-01 

25 2.4 .242 E-02 

25 2.4 .243 E-02 

25 2.4 .438 E-02 

25 2.4 .872 E-03 

50 1.2 -.814 E-02 

50 1.2 -.758 E-02 

50 1.2 -.800 E-03 

50 1.2 .392 E-03 

50 1.2 .123 E-02 

50 1.2 -.279 E-03 

77 

Comment 

Monotone 

Monotone 

Oscillations 

Oscillations 

Monotone 

Monotone 

Monotone 

Monotone 

Monotone 

Monotone 

Monotone 

Mono tonI! 
>. 

~ 

(lj 

... 

10 

I:: 

.... 

't:I 

k 

0 

0 

U 

.-4 

1'3 

g 

k 

0 

z 

10 

8 

6 

4 

2 

o 

-4.0 -3.0 

ORIGINAL [.J ..... Q~ I~ 

OF POOR QUALITY 

a - B .. 1.0 

4.0 

y .. -2.0 

cS .. 0.4 

-2.0 -1.0 o 1.0 

Cross flow velocity, w 

Figure 1. - Cro~s flow velocity profiles as a function of time. 

78 

2.0 
~ 

::l 

C 

..... 

... 

0 

... 

... ... ~ 

\D X 

~ 

::<: 

,~ 

~ 

,/ ... __ .... ___ .I!--:''''~-' -~ ..... -._. -: ~ ~~=-:-.-:- _:':":~·.::"':~.::"':":"::=:V. 

0.012 

0.008 

0.004 

1.2 

-0.004 

1 cycle 

Time linearization 

(lit = .04) 

Time linearization (lit· .02) 

1.4 1.6 1.8 

Quasilinearization 

(lit = .04) 

2 iterations 

2.0 2.2 

2 cycle 

a- 1.0 

B - 4.0 

y • -2.0 

6· 0.4 

2.4 2.6 2.8 

-o.ooal, ___ -=-________ -=::::::-. ________ =:... __ ---l 

Figure 2. - Maximum error of u as a function of time. 

00 

"";;t: 

"DO 02 

0::;, 

.:0,-

.o"C 

c).:-

):0 r:;. 

r r:1 

~C 

~, 

~: 

'I 

!' 

.f 1; 

, 'l~' 

'4' f 

,Ij 

'I, 

f· 

~': 1~, 

...... ~ 

"1. ,., 

.f; ')-' '11"'- /:', .. -.,; . .'. 

1.1J' l \;; 

·t':: 

i;~ 

1,'; 

\~ 
'. \r 

... \" :!.~ ~::-iJ~~i:t2~lt.~~~~""<W~~~;;:h-V~.n""b~S('4t Io,·~r~ ... ~t.l .. '· :tJ" ',', .. ' .... ..... """'" 

,. 

/ 

/ 

,/ 

'.,. 

\ 

.~-

'\ '\ 

\\ ... '.t' 

\ 

'i -~ \ , \" 

I 

, : 

:;.} 

,'}" J.., 

~ /.' I 

.. 

i 

1 

I /, 

i , 

\' /\ 

" 

• 

:J 

c:: .... 

'" 0 

'" ... 

ClI 

e 

:l 

e 

.... 

x 

'" :<: 

:' 

0.016 

0.012 

0.008 

0.004 

0 

-0.004 

-0.008 

ORIG1N.~:' ~:;'l;~:: IJ 

OF poon Qtli-\L: N 

Reaches maximum of .0256 

at t = 1.84 

lineariza t ion 

(t.t = .04) 

-O.012~ _____________________________________ ~ ________ ~ 

Figure 3. - Maximum error of w versus time. 

80 
, " 

_.--

\ 

...... -. 

" 

, " 

--.--

. ' , 

' ....... - ..... I "-' I 

~ 

0 

~ 

~ 

CI 

El 10-2 :> 

El 

x 

r: 

:<: 

o 

.'- ._---_. "-- .-"- --. 

10 20 30 

Hesh intervals 

a· 1.0 

B. 4.0 

Y • -2.0 

6· 0.4 

Second order 

central diifcrcnces 

40 

Figure 4. - Solution of w equation, maximum error vs. mesh intervals • 

81 

50 
-. 

~ 

/ 

/ 

.. 

, 

.I 

\ \ 

10-1 

o 10 

ORlGIH;:'t. P.~, ~.'.~ ~.: 

')F POOR r.:.!;\L!"j"{ 

a a 

r. a 

Y " 

6 . 

1.0 

4.0 

-2.0 

0.4 

Sl'("ond order 

Fourth order 

OCI-Gl 

20 30 40 50 

Figure 5. - Solution of u pquation, m~ximum error VS. mush intervals • 

82 

" 
.. 

\ 

. \ 

\ 

.. ~'-\ 

'.'. 

. "', 

>,t·. 

.:.,. .... 

' . 

" 

.'. 

I 

.. ' 

': '. . '. 

,':' 

. .. .' ,-

" .. ~ . 

•... 

';', 
End of Document