NASA Contractor Report

ICASE Report No. 89-57

F

181903

_L

ICASE

A HIGH-ORDER LAGRANGIAN-DECOUPLING METHOD

FOR THE INCOMPRESSIBLE NAVIER-STOKES EQUATIONS

L. Ho, Y. Maday,

A. Patera, E. Ronquist

Contract No. NAS1-18605

September 1989

Institute for Computer Applications in Science and Engineering

NASA Langley Research Center

Hampton, Virginia 23665-5225

Operamd by He Unive_ities Space Research Association

: (NASA-CR-IdlgO3) A HIGH-ORDER

LAGRANGIAN-DECOUPLING METHOD FOR

INCnMPRFSeI_LE NAVIPR-STOKES

Report (ICASF) 48 p

N/ SA

National Aeronautics and

Space Administration

Langley Research Cenler

Hampton, Virginia 23665-5225

TH_

FQUATIONS Final

CSCL 12A

G3/o4

NQO-I0631

Unclas

0_32363
......... L

7j

m _
A HIGH-ORDER LAGRANGIAN-DECOUPLING METHOD FOR THE

INCOMPRESSIBLE NAVIER-STOKES EQUATIONS

Lee-Wing Ho t, Yvon Maday t* 1, Anthony T. Patera t, & Einar M. Rcnquist t

tDepartment of Mechanical Engineering,

M.I.T., Cambridge, MA 02139, U.S.A.

* Laboratoire CNRS d'Analyse NumSrique,

Universit6 Pierre et Marie Curie 75252,

Paris, Cedex 05, France

ABSTRACT

In this paper we present a high-order Lagrangian-decoupling method for the unsteady

convection-diffusion and incompressible Navier-Stokes equations. The method is based upon:

Lagrangian variational forms that reduce the convection-diffusion equation to a symmetric

initial value problem; implicit high-order backward-differentation finite-difference schemes

for integration along characteristics; finite element or spectral element spatial discretizations; mesh-invariance procedures and high-order explicit time-stepping schemes for deducing

function values at convected space-time points. The method improves upon previous finite

element characteristic methods through the systematic and efficient extension to high order

accuracy, and the introduction of a simple structure-preserving characteristic-foot calculation procedure which is readily implemented on modern architectures. The new method is

significantly more efficient than explicit-convection schemes for the Navier-Stokes equations

due to the decoupling of the convection and Stokes operators and the attendant r_n

temporal stability. Numerous numerical examples are given for the convection-diffusion and

Navier-Stokes equations for the particular case of a spectral element spatial discretization.

1This research was supported in part by the National Aeronautics and Space Administration under

NASA Contract No. NAS1-18605 while the second author was in residence at the institute for Computer

Applications in Science and Engineering (ICASE), NASA Langley Research Center, Hampton, VA 23{}65.
1. Introduction

A large class of important fluid flows is described by the incompressible Navier-Stokes

equations,

um,t + up u,.,,,_ = -p,,,, +vum,ff + .f_ in

u_,q =0 in f_ ,

where z_.= (z,,,)m =1 ..... _ is tile space variable in ft of/R d, t is time, ]'(z__,t) is the prescribed

force, u(_, t) = (u,_ (z_, t))m=l,...,d is the velocity, p(_.z, t) is the pressure normalized by a

(constant) density, and v is the kinematic viscosity. We adopt Cartesian tensor indicial

u,,, = _ and summation over repea_ted subscripts, u_,q notation, with u,,,,f = 0®, , ,t ot , =

d °-Y--t Equation (la) represents conservation of momentum for a Newtonian fluid,

_=t O= w '

while equation (lb) represents conservation of mass for an essentially incompressible flow.

TILe difficulties associated with numerical solution of the incompressible Navier-Stokes

equations arise from two distinct sources: tlte pressure-divergence terms, and the convection diffusion contributions. The pressure-divergence terms require first, for optimal

convergence, that proper discrete function spaces be chosen for the velocity and pressure

[1], and second, for rapid solution, that algorithms be developed that effectively decouple

the pressure and velocity [e.g., 2-4]. The origin of the numerical difficulties associated

with the convection and convection-diffusion terms are equally well documented, but less

resolved: convection-diffusion introduces thin boundary layers for small t_; the imbalance

of convection and diffusion leads to weak spatial "stability" for small r,, in that tile ratio

of tile continuity and coercivity constants becomes large [5]; timescales for equilibration

become large for smaU _,; the convection term destroys the isotropy and symmetry of the

Stokes operator for all values of t_.

To motivate our new class of schemes we briefly summarize standard numerical approaches to the time-dependent Navier-Stokes equations. In particular, we note that_ the

unsteady Navier-Stokes equations are often treated in a semi-implicit fashion, in which

2
tile symmetric Stokesoperator is handled implicitly and tile convection term is treated

explicitly,

_+1 _ _ + jr([_,_.,,]", [u,_._,,1"-1 ...) =

At

_n +1

-e,_ + _(I.u.,,, +/_]"+_,l-u.,. + 1_]",...) in _ (2a)

u,,+l 0 in f_ (2b)

q_g

Here u_" (x) = u(__x, nat), where At is the time step, and jr and _ are appropriate explicit

and implicit time-marching schemes [6], respectively. Explicit treatment of the convective

term is motivated by the fact that solution of tile (nonlinear) implicit convection equations is difficult; no efficient (faster-than-time-like) robust iterative methods are available

for general unsymmetric anisotropic operators, and direct methods are typically memoryintensive, costly, and serial, ill particular in higher space dimensions. These arguments

for explicit convection, although by no means universally accepted (e.g., [7]), are most

easily defended for high-order spatial discretizations, in which asymmetry and bandwidth

problems are further aggravated by longer-range coupling.

The semi-implicit method of equation (2) is intended to represent a compromise of

convenience and efficiency between the symmetric Stokes operator and the unsymmetric

convection term; unfortunately, it is not an optimal compromise, as we now describe. To

begin, we assume that the semi-discretization (2) is further discretized in space, with the

spatial discretization characterized by an effective mesh spacing h. We next (plausibly)

assume that the critical time step set by the explicit convection treatment,

At.. ~O(hll__l),

is smaller than the time step required for accuracy, At,_,, and that stability titus determines

the timestep. The work required to integrate (2) to a final time T, )'r ''_" (semi-implicit),

is then

W'"" = (T/AG,))/VA, (3a)
where kVAt is the total work per timestep (e.g., ill ,clock cycles). The total work per

timestep WAr is further broken up as:

• Al_ai _A)inv WAr : avvAt,: "t- rvAt,, (3b)

Here the convective work per timestep is given by }/VAt,, = ,,vAt,,,_AJ'_'twhere _aJ,_,_,vAt,, is tile

work required to evaluate the convective terms, and a is the (order unity) number of

evaluations required for Y, and the Stokes work per timestep, }/VAt,,, is equal to the work

required to invert the Stokes system, kVAt,,.

Finally, we note that for nontrivial problems it is typically the case that w_,_,z • v/xt,c < <

Vt)_t_,,, as the solution for the pressure in incompressible simulations requires at least an

elliptic solve [8-10], and, more properly, nested elliptic solves [2-4] at each timestep. For

instance, for typical spectral-element-discretization Navier-Stokes calculations, _'_'_'_,,A,,,is

an order of magnitude less than kV_,_, for splitting Poisson schemes [11-13], and several

orders-of-magnitude less titan kV&t,, for consistent Uzawa methods [4,14,15]. (Note that

explicit treatment of the viscous terms does not significantly alter these estimates, as

the majority of the work kV_c_, derives directly from the elliptic pressure term.) These

complexity estimates are based on relatively efficient conjugate-gradient [14] and multigrid

[16,17] elliptic solvers ]4,14,15], and although kY_,_,0 will certainly be reduced by further

improvements in elliptic technology, it is clear that kyats0 will always remain significantly

larger than r rAt,c.

It is now readily seen in what sense semi-implicit methods are not entirely satisfactory.

The usual motivation for turning to an explicit rather than implicit formulation is that

the increase in number of timesteps (T/At) is more than counteracted by an associated

decrease in work per timestep, WA,. ttowever, it is clear from (3) that for the NavierStokes equations an increase in T/At from T/At,,, to T/At,, is not compensated for by

a commensurate or greater decrease in }A]A,, as the majority of the work per timestep is

due to the elliptic solves, WAr,° >> }A)At,,. A better compromise is required, in which

4
T/At << T/At,,, and _V,5,_., and VPAt., are of the same order; such a "balanced" scheme

would clearly result in improvements in computational efficiency, in particular for highorder spatial discretizations. These comments are especially relevant when (2) is viewed a

method for the solution of the steady-state equations.

We present in this paper a new high-order Lagrangian-decoupling method for the

convection-diffusion and Navier-Stokes equations based on high-order finite-differences in

"time" and variational (finite or spectral element) discretizations in "space". The scheme

preserves symmetry as regards the diffusion and Stokes operators, thereby allowing for fast

and robust iterative solution, yet provides for much better temporal stability for the full

Navier-Stokes equations than explicit-convection alternatives such as (2); in essence, tile

new method exchanges a significant increase in At above At,, for an easily accomodated

increase in the convective work that brings }'VAt,, into balance with W,',t.0. Tile direct ancestors of the method are the Lagrangian and arbitrary-Lagrangian-Eulerian finite element

techniques [18,19], however the technique can also be thought of as a rigorous subcycling

approach [20], or as a high-order non-dissipative upwinding scheme [21,22]. The essential

differences between our method and its closest predecessor, the characteristic scheme of

Pironneau [18], are the systematic and efficient extension to high-order accuracy, and a

simple structure-preserving Lagrangian remeshing scheme (characteristic-foot calculation)

which is readily implemented on modern architectures.

The structure of the paper is as follows. In Section 2 we introduce the strong and

Lagrangian variational forms of the convection-diffusion equation. In Section 3 we describe

the nodal function spaces for finite [23] or spectral [14,24] element spatial discretization; we

present the new high-order finite-difference-in-"time" Lagrangian-decoupling method; and

we summarize the stability, accuracy, and computational complexity of the new technique.

In Section 4 we briefly describe an extension of the scheme to the full unsteady nonlinear

Navier-Stokes equations. Lastly, in Section 5, numerous numerical examples are given of

tlte new scheme for the particular case of a spectral element spatial discretization.
2. Convection-Diffusion Equation

2.1 Strong Form

In this paper we shall adopt the common practice of using the convection-diffusion

equation as a simple model for the Navier-Stokes equations (1). As our scheme is particularly tied to the structure of the incompressible Navier-Stokes equations (the inequality

]/Vat,, << )'VAt,, is only occasionally true for convection-diffusion alone), it is imperative

that we return to the parent equation to demonstrate that our method is, indeed, applicable to the full Navier-Stokes system. This extension is described and illustrated in Sections

4 and 5, respectively.

Our simple convection-diffusion equation is given by

¢,, +u,¢,, = re,,, + a Vt e (0,T), V_xe fz , (4)

where ¢(x_, t) is a passive scalar, u is the prescribed convection velocity, x is space, t is

time, v is diffusivity, g(x, t) is a prescribed source, and fl is the domain in/R d. In addition

to equation (4) we require initial conditions

¢(z,t=0)=0 v__ef/ , (5)

and boundary conditions

¢(_,t) I,,r_= o vt _ (o,:r) , (6)

where Of/ denotes tile boundary of f_, and _h will denote the unit outward normal on Of/.

The simple initial conditions (5) and homogeneous essential boundary conditions (6) are

chosen for clarity of presentation; the method, in fact, applies to general non-homogeneous

and non-essential bouILdary conditions [25].

For the passive scalar equation (4) the convection velocity u_ is prescribed. The field

u.U_is required to be solenoidal,

7L_,_ = o (7)
to vanish on tile domain boundary,

u leo = 9_ (8)

and to be independent of time,

u = _(__) (9)

The requirement (8) is closely coupled to (6), and can be relaxed to include inflow and

outflow conditions [25]. Time-dependent Lt will be considered in Section 4 in tile context

of the Navier-Stokes equations.

2.2 Weak Forms

Our numerical formulation rests on two Lagrangian variational forms, ttle first for

the convection-diffusion equation (4), tile second for a seemingly trivial equation which

we shall denote the "invariance" equation. To better motivate these variational forms,

we remark that our numerical scheme, like previous finite-element characteristic methods

[18,19], corresponds to finite-difference methods in space-time along u-characteristics, and

variational methods in space. The purpose of the Lagrangian variational forms is first,

to (Lagrangian) transform "space" into "time" as regards the convection operator, and

second, to allow for (variational) finite or spectral element spatial discretization of the

diffusion operator.

a. Convection-Diffusion Equation

To pose tile variational problems we first define the classical Sobolev spaces L2(I2)

and HI(a)[z61,

L2(a) = {v(x) measurable in fl, f v2dgl < c_}

f_

H01(fl) = {v(x) e L2(fl); v,, e L2(fi), q = 1,... ,d, v 10t_ = 0}

(lo.)

(lob)
In what follows we shall denote the space Hl(fl) as Y. We shall also require the temporalspatial space Z defined by

Z = {v measurable in

T g

q--1

v(__,t)• n_(n) .,_. t • (O,T)} (11)

where T is the final time of integration.

denoted by II" IIL,, II. Ilu,, and II. Ilz, respectively.

The Lagrangian variational form of the equations (4)-(6) is then:

that ¢(x, t = 0) = 0, and Vt • (0, T)

_{ fnC¢_n) =-. fn¢,,¢,,dn +fn_gdn

for all ¢ • Z such that 6(x,s = 0) = ¢0(x), and Vs = T - t • (0, T)

Tile norms associated With these spaces will be

Find (;b • Z such

(12)

6. -_,6,, = 0 Vx• n (13)

for some Co(x_) • Hol(f_). Here ¢(x,s) = C(x,T - s). Equation (12) is readily derived

from the strong form (4)-(6) by integration by parts in space and time, use of "Reynolds

transport theorem", and substitution of equation (13) [27l; note that the particularly

simple form of the first term in equation (12) requires a solenoidal convecting velocity field

u. We now suppose that the domain _ is the image of some reference domain _ through a

smooth, invertible mapping A: _ D. The equation (13) for the test function ¢(x_, s) can

then be written in Lagrangian form,

{ ax---ia s) = -u(X_(a,s)) Vs • (0, T), Va • _ (I4a)

_, X--,

_(x(__,_),_) = 0

X(a,, ^ = 0) =

¢(z__,s= 0) ¢0(x_)

where X__ is a standard Lagrangian spatial variable.

V_a • _ (14b)

Vz_.• D ,

The variable X__.will remain smooth

for smooth u.
Note that equations (12) and (14) are the variational restatement of the simple physical

fact that the rate of change of ¢ in space-time along a characteristic of the convection

equation is equal to the divergence of the diffusion fluxes. In our numerical scheme the

convection initial value problem will be treated by high-order finite-differences, and the

diffusion boundary value problem will be treated by finite or spectral elements.

b. Invariance Equation

As with all Lagrangian schemes, some form of remeshing or characteristic-foot calculation is required for tlle scheme to proceed for long times. A key element of our high-order

scheme is a remeshing strategy which is both efficient and accurate; the method is based on

an invariance procedure, in which we associate to any continuous function q(x__) E H01(f_)

a new function q(x, s) and related (trivial) evolution equation

q(__,_= 0) = _(_), _q(z,_) = 0 v_ e _ (15)

Equation (15) can be written in a u-Lagrangian variational form as: Find q(x, S) E Z such

that q(x,s = 0) = q(x), and Vs E (0, S)

._d {/pq} + f 1,u,_q,,_ =0 (16)

t2

for all p E Z such that p(x, s = 0) = p0(x_), and Vs E (0, S)

p,, - u,,,p,,,, = 0 Vx e s2 (17)

for some p0(x) E HI(o). Equation (17) for p(x_, s) can be written in Lagrangian form,

{ _-(_,s)=-__,(x(._, s))

_p(X_(.,_), _)= o

{ x(_,s = o) = __ p(_, _ = 0)= v0(_z)

with _ defined as for equation (14).

w e (0,s), Va_e fi (18a)

(18b)
The equivalence of ((16), (18)) and (15) can be understood by noting that the u--

motion of the "mesh" through (18) is cancelled in (16) by the u-convection term. Although

it may appear overly complicated to re-introduce Lagrangian convection by the u--field

into an equation (15) which is essentially Eulerian, this will allow us to accurately remesh

variables while preserving the underlying structure of the spatial discretization. Indeed, we

shall see that the Lagrangian simplification of (4) in (12) in exchange for the Lagrangian

complication of (15) ill (16) is equivalent to striking a better balance between At and

W ,.,,,./ WAt,, .

Remark on Inflow/Outflow Boundary Conditions: We note that although we have

implicitly written tile hyperbolic equation (16) as having Dirichlet boundary conditions

on 0f_ (i.e., q(x, s) E Z), this is, in fact, equivalent to the proper choice of no boundary

conditions (recall u._ fir,, = 0), given that q(x) C H01. In treating inflow/outflow problems

(in which we impose essential boundary conditions on ¢ at inflow, u,,,h._ < 0, and natural

boundary conditions on ¢ at outflow, u._ fi,,, > 0), we require boundary conditions on q in

(16) only at inflow, with this inflow value obtained through extensiofi of the inflow profile

upstream.

10
3. Discrete Equations

3.1 Spatial Discretization

a. Discrete Function Spaces

We first need to generate ttle discrete function spaces ("nodal finite element spaces")

associated with tile space Y = H01(f_), after which our method will follow directly from tile

Lagrangian variational forms of Section 2.2 and high-order finite-difference discretization

of the "time" terms in equation (12) and (16). We begin by assuming a conforming,

non-overlapping domain decomposition [28] :D,

_= U _ , (19)

in which the 12_ are defined with respect to reference volumes _/k (e.g., a segment in/R 1, a

triangle or square in/R 2, a tetrahedron or cube in/R 3) by an invertible elemental mapping

(in fact, .4 of equation (14)),

c _ _:LL _-)_,_ e fi_ v,. e {1,... ,d}. (2oh)

Tile decomposition is geometrically conforming in that we require that the intersection

F = (gf__ 17 (9121 be either null, or an entire face, edge, or vertex of both f_h and _'lz.

Within each reference domain _, k = {1,..., K}, we introduce a polynomial space

/P_(fi_), a unisolvent set of N _ basis points _'A_, m E {x,...,d},i e {1,...,N'}, and a set

of associated Lagrangian interpolants Qi,_,

Q_,_ e/P_ (fit)

Qi,_(_i,_)=_,i Vi,jE{1,...,Nk} 2

, (21)

11
where 6_i is tile Kronecker delta symbol. For functional conformity we require that tile

/pk(fi_) and F '_ be chosen such that for all k E {1,...,K} and for any v({_) E /P_(fi_)

and any x E F = 0{2 _ O Of_ t such that P is non-null there exists a w E h°l(_ t) such

that w(F-!'a(x)) = v(F -1 '_(x)). With our nodal Lagrangian-interpolant representation

we now make precise tlle particular isoparametric map FA _ of interest in our Lagrangian

characteristic technique,

N j,

FAJ'(__)-__, X_,_O"'(_) Vm e{1,...,d} Vke {1,...,K}, (22)

i=l

in which the X _'_ are the images of the _'_.

i

In addition to the conforming assumption on the /P_(_) and F 'A we require: the

_,k be chosen such that those _i._ on O_ _ are unisolvent with respect to the trace space

of/P_(_); the usual "local-global" coincidence condition on the X_._'_ (e.g., for each _i'_

on an edge shared by elements k and l, there must exist a unique _i.I which maps to the

same global point),

x_ ,_ = x__' w e {_,... ,,_'}, v(d,k) e s, (23)

Here A/" is the number of global nodes in the system, and Si is the set of (local node,

element) couples associated with the unique global node X._.. _ Note that although we

introduce here the standard finite element global identifier S¢ to describe inter-element

connectivity and continuity, this construct serves purely for purposes of succint presentation; in practice, our schemes are based entirely on element-local constructs (see Section 5)

[13,14,29]. Summarizing our superscript notation: X refers to a globally defined quantity;

X '_ to an elementally defined quantity; X i'h to a nodally/elementally defined quantity;

and X i to a nodally/globally defined quantity.

We can now define our discrete polynomial subspace Y'h(X') C Hl(f_) = Y as

(

0,_)_s,

12
where Vk _ {1,...,h'} and Vj E {1,,..,N _}

Qj,_(£-I ,_(x_)) v_ e (24b) 0_"(_-)= 0 V__ ,

where the set {1,...,.Af}/s excludes those i E {1,...,.N'} for which X._._ is on Off. The

space Yh(X _) is thus characterized by the domain decomposition _D, by h, symbolic for the

/pi(_i), k = 1,..., K, and by the the global nodes X___. The basis we choose to represent

va E Yh is the nodal basis used to describe the space in (24)

_f

v,(_) = _ v_ _ 0J,_(_) (2_)

i=l 0,*)Es,

V 'h - o ie {i ,... ,W}_s , (25b)

' = _'h(X'). where the set {1, ..., A/'}_ s is the comlflement of {1, ...,]g'}/s in {1, ...,A/}, and v h

For " and explicitly write q(X i) to denote the nodal values simplicity we shall write v h = v',

of a function q not a priori in Yh (X_).

b. Discrete Inner Products

Having defined our spaces in Section 3.1a, we now complete the spatial discretization

by defining the discrete inner products associated with the continuous inner products

required in the variational forms defined in 2.2. To begin, we remark that the convectiondiffusion equation (12) requires two inner products, the L 2 inner product and an H 1

bilinear form. Their discrete forms with respect to our space Yp, are

W 2¢

V_b',¢ e (ff_Ar) 2 E E _biSii(.Xl)dP 1 =--

i=l j:l

E E E _ _ _ , (26)

i:1 ./:1 P,,.J' "*'

{e.b)es i

and

v_', ¢J e (_'¢)2

._" j¢

i:l j =1

13
1¢ Ar

i=l j=l ,.,.b,.,. *

(r,k)E$,;

, (27)

respectively. Here J ,k (_) and G_h_,,(__) are tile Jacobian and tile derivative transformation

matrix associated with tile change of variables (22),

J'*(_) = det F '_ (_) (28a)

and

j ,*)-1F,,,({) ,' = (2sb)

respectively. The invariance equation (16) requires the L 2 inner product (26), as well as a

convection operator,

V¢',¢ i E (/R_¢) 2 _ __, ¢'DO(X')¢ j =_

i=.l j-=l

/v" /d

E E ¢'¢J E fi'_ u.Q'"G,_O"'d_

i'-1 i=1 ,.,.b,.,.

[.J,)E$i

(29)

In (26)-(28) _ refers either to exact quadrature or numerical quadrature, depending on

the particular spatial discretization chosen (see Section 5).

Lastly, we require the linear form

i--1

()'p,*'es, " - -

(30)

for the inhomogeneous term in (12).

3.2 Full Discretization

Our fully discrete convection-diffusion equation now corresponds to insertion of, the

discrete linear _md bilinear forms Of (26)-(30) into the variational form (12), followed by

14
discretization of tile "space-time" derivative associatedwith the left-hand side of (12) by

the Qth order backward-differentiation implicit finite-difference formula [6]: Find _;,_+1 =

¢i;.+1 = ¢/(t.+l = (n + 1)At) such that

Z &B'J (x';"+1-')¢ j;"+l-' =

i=I g-:-0

fl'{_--_ A'J(X_;'*÷I)¢ J';"+l + _(X____/;"+l)g;"+l} , Vi E {1,...,J_f}/a (31a)

i--1

¢,;.+1 : 0 vi _ {1,...,N}_, (31b)

Here X z;'_+l = X_.t;* refers to the the base ("undeformed") discretization of _, and thus

the Aii(XI;'_+I),gi(X _;"+1)_ __ are known; it can also be shown from the incompressibility

of u and (14a) that B_/(X I;'+1-_) is, in fact, independent of X_J, and can thus be taken

as B_i(XJ;'*+I). Note that any nodal variable without tilde refers to nodal values at the

base-discretization points X_ i;'. The backward differentiation coefficients are given in [6]

with fl, .-_ O(1/At), and 8' " O(1).

It remains now to find the _i;.+1-, for q > 0; _;.+l-q represents the value of ¢ at

time t "+l-q at tile "foot" of the characteristic whose "head" at time t "+1 is at X__J;" [18].

To this end, we use our invariance procedure (16)-(18) to determine _;.+1-_ from the

values ¢/;.+1-,, the latter being known from previous timesteps. To wit, at each timestep

t "+1 we perform the following subproblem: Vq 6 {1,..., Q}

_i;n+l-e;m:O : ¢i;.+1- f )

___,;;._=0 = _(x__';')

Xi;;m=O __ X_; •

Vi e {1, ..., A/'} (32)

15
vm = o,...,qatl_ -1 ,vi a {1,...,X}

1 B'J (X._l;;m"t'l)_ j;"-'t-l-';m'l'l - Z Bq(X---l;;'a)'t_ j;''l-l-tt;m "-

As

j :0 i =0

2

p :==0i..- 1

I (_ _ij(__XJ;;m+l)_j;;m+l _ _ij(Xz;;m)_j;;m)

j "-0 / =0

2 .,V"

p ===:0,_=--1

(x_;;m+l - x _;;')/A_ = - _ "r,_;;"-' ,

p-.=O

(33a)

(33b)

(33c)

(34)

We have introduced here an explicit third order Adams-Bashforth scheme in s (with

coefficients % [6]) to solve the invariance problem. The subproblem timestep is As < At,

with _- assumed integer; the notation ¢_;'_ refers to timestep n, while _;";'_ refers to

sub-timestep m for a given t". Equations (33a) and (33c) correspond to discretization

of (16) and (18) respectively; in conjunction with (34), (33a) and (33c) represent the

Lagrangian requirement (14). Determination of the _ requires application of invariance to

the velocity as well, (33b), as the velocity enters into the convective operator D _i, (29),

as well as the mesh update equation for X _. Note that, in practice, we do not assume

Bi_(X ';;'') = B_(X t;;') in (33a) and (33b), as actually calculating Bq (X_.';;_) from (26)

appears to improve stability.

The complete scheme comprises: Lagrangian variational forms that reduce the convection-diffusion equation to a symmetric initial value problem along characteristics; highorder, implicit finite-difference sclmmes for integration along characteristics; standard

16
finite-element-like discretizations for the symmetric boundary-value-problem spatial operators; mesh invariance and high-order (explicit) time-stepping for deducing the function

values at previous space-time points. We now turn to a summary of the characteristics of

the method; the next section is a heuristic anticipation of the numerical experiments of

Section 5, and of theoretical analyses treated in [25].

3.3 Stability, Accuracy and Computational Complexity

a. Stability

For tile backward-differentiation scheme (31) unconditional absolute stability can be

intuited in several ways. First, for tile Q = 1 scheme, energy arguments are readily

constructed [18]. Second, for any Q _< 6 tile region of absolute stability in tile complex

hat plane for the model problem ut = /_u for the backward differentiation formula includes

the entire negative real axis, implying that for diffusion partial differential equations these

schemes are absolutely stable. As (31) is a diffusion equation along characteristics, we

expect unconditional stability. Indeed, it is readily shown from von-Neumann analysis

that for the simple case of Fourier spatial discretization the modulus of the backwards

differentiation "temporal" growth factors for space-time integration of convection diffusion

are identical to those for time integration of standard diffusion; as expected, convection

effects only the phase of the propagator.

As regards the application of the third order Adams Bashforth technique to the invariance equations, (33), it is clear that the usual Courant condition will apply, giving a

condition of the form As < O(h/I u_l). We note that our Lagrangian-decoupling method

bears a resemblance to earlier subcycle proposals [20]; an essential difference is that (33) rigorously and stably decouples the convection and elliptic (or Stokes) contributions, whereas

subcycle approaches involve repeated application of the Convection operator in splitting

fashion.

17
b. Accuracy

Convergence of ¢h to ¢ in H 1 at a fixed time T (or in the Z norm) is achieved as

At --, 0 in (31), As _ 0 in (33), and as Yh ---* Y. The error due to space-time integration

by the backward differentiation formula in (31) is expected to be of order (At) Q. The

ability to obtain virtually arbitrary-order accuracy constitutes a major difference between

the current technique and past methods [18}; high-order Adams-Moulton schemes would

not be stable, and would also require multiple Laplacian evaluations on deformed (not

X _;') domains. As concerns the temporal accuracy of the Adams-Bashforth treatment of

tile invariance equation, (33), we expect O((As) 3) errors. TILe fact that the error can

be controlled through As and a high-order time-stepping scheme represents a significant

improvement over earlier remeshing schemes which are typically O(At) [19]; low-order

remeshing methods can result in strongly accuracy-limited timesteps.

The spatial errors are proportional to the error in the best-fit approximation of the

solution ¢ by the underlying discrete space Yh [14,18,25]. In particular, for h-type finite

element refinement [23] (K _ oo, N k fixed) we expect algebraic convergence, and for

p-type [30,31], or spectral element refinement [14,15,32] (K fixed, N k _ oo) we expect

spectral convergence (e.g., exponenti',d convergence for infinitely smooth functions). Note

that the spatial discretization does not only effect the treatment of the diffusion operator in

(31): spatial discretization controls tile numerical diffusion and dispersion in the invariance

procedure (33); the allowable At in (31) is (accuracy, not stability) limited by the degree

of geometric distortion that can be represented by the spatial approximation [25,33].

Since the characteristic procedure transforms space into time in the convection operator, a "time"-stepping scheme for (31) limited to low-order approximation would potentially preclude the use of a high-order method in the remaining spatial operators. The

fact that the current scheme is high-order in time allows methods that are high-order in

t

space to once again be considered in a characteristic framework. Although it might appear

18
that a variational (e.g., spectral) method in spacewould suggesta matching variational

(spectral) approximation in time [24], the initial-value-problem nature of the convection

operator indicates that finite-difference approximation will be more efficient, with none

of the attendant "complex geometry" objections typically raised as regards high-order

finite-difference approximations in space.

Lastly, it should be remarked that even tile low-order (Q = 1) scheme can be quite

accurate. In particular, we note that if the invariance procedure is exact (e.g., As small),

then the temporal error in (31) is due solely to inexact integration of the diffusion terms

along a characteristic. This contrasts sharply with conventional upwinding schemes [21],

in which the majority of the numerical diffusion derives from mesh-induced inaccurate

representation of the characteristic foot. As a result, low-order forms of (31)-(34) can

achieve high accuracy if: gradients in ¢ are small, and the diffusivity t_ is small; gradients

in ¢ are primarily in the cross-stream direction, that is, orthogonal to u [5]; gradients in ¢

are large in the flow direction, but the flow velocity (and hence effective spatial increment

Ax = l uU.[ At) is "small" (such as in stagnation-flow boundary layers).

c. Computational Complexity and Generality

The beauty of tile characteristic scheme [18] is that it reduces each timestep to a

symmetric elliptic solve (in the case of Navier-Stokes, a symmetric Stokes solve) that can

be very rapidly solved iteratively [4,14,16,17] ; that is, the asymmetrizing and destabilizing

effect of the convection operator is eliminated. For our particular implementation, the price

to be paid is the solution of a "bare" (non-pressure-corrected) convection operator (33).

The key to the computational efficiency of the scheme, as described in detail in equation

(3) of the Introduction, is that the increase in work associated with (33), WAr,,, is more

than compensated for by the large At allowed in (31). In particular, we notethat the work

19
for the characteristic scheme,YY 'h, is given by

T T

wch I'L A vgv a| i_,¢ = _ _WA,,. , (35) At,, a vvAt,o + At'.:,

where tile first term represents (33), and tile second term (31). Here At,, is as in (3),

a' is inflated from (3) by virtue of tile additional work associated with Q > 1 in (33),

and At',,, is tile new accuracy limit that reflects tile effect of geometric distortion and the

accuracy of tile Qth order backward differentiation scheme. Although many factors must

be considered to determine whether, at fixed accuracy, )4/"h is much less than )IV°'_', the

potential clearly exists for great savings. In particular, the two terms in (35) are better

balanced than in (3), and the dominant work per timestep, _'Y_u, is effected much less

frequently, T/At',_, << T/Ate,. In summary, the stability-limiting and work intensive

parts of the calculation have been Lagrangian-decoupled.

A central aspect of the efficiency of the method is the invariance procedure. Previously proposed schemes are low-order [18,19], or use a nonlinear search and solution

algorithm to find the ¢ directly. In the latter approach, [(33a), (33b)] are replaced by a

"nonlinear solver" to obtain the _ directly from the X _ given by (33c) [18]. Although the

"nonlinear solver" may appear faster that the time-stepping approach (33), in practice the

opposite will typically be true: first, even for low-order approximations, the "nonlinear

solver" approach requires non-elemental, unstructured calculations, leading to inefficient

vectorization and parallelization - the time-evolution solver preserves topology and structure; second, as tile "nonlinear solver" approach destroys discretization structure, tensor

product representations will not be preserved, leading to orders-of-magnitude increases

in work for high-order approximations ([14,34], see Section 5) - the time-evolution solver

maintains an initially tensor-product representation for all "time" s.

Lastly, we briefly comment on generality. Iterative solvers, unlike direct solvers, are

typically non-robust with respect to modifications to the equations and physics, as their

convergence (rate) is dependent not only on discrete-equation structure, but also on tile

20
structure of tile spectrum. It is usually tile case, however, that increased physical complexity in incompressible continua enters into the "divergence of the flux tensor" portion

of the equation, not the acceleration term, and thus the Lagrangian-decoupling method

should be generally applicable to a large class of problems.

21
4. Navier-Stokes Equations

In this section we describe tile extension of the convection-diffusion scheme of Section

3 to the full nonlinear unsteady Navier-Stokes equations (1). There are many ways in

which to deal with the nonlinearity; we propose here a linearization technique which is

both simple to implement and illustrative of tile general procedure.

To begin, we consider a semi-discretization in time, in which u,, (x) -um (x, n'At)

are known for 0 < n' < n, with u_,,,, = 0 Vn' E {O,...,n}. To advance the solution,

that is, find u_+l(x), we write u,_(x,t) = u,_ (x)+ 6u,_(x,t), and choose the particular

linearization of (1) that yields

u,,,,t + uq um,q - -p,_ + r'u,,,,_ + f,_ in [2 (36a)

uq,q = 0 in ft (36b)

for t near t". The line arization (36) is preferrable over other possible distributions of 6u,_ in

that stability is readily shown, _ ff_ u,_ u,_ < 0. We now apply our Lagrangian-decoupling

method and variational spatial operators exactly as for convection-diffusion, except that

now the velocity appearing in the convection operator D ii is frozen at u_.."during the substep

(32)-(34), and the Laplacian operator in (31) is replaced with an appropriate symmetric

Stokes discretization [14,15]. Note that, unlike the compressible flow case, (33a) has only

a single "characteristic" direction, _A, shared by all components of the velocity.

As regards the accuracy of the resulting scheme, it is clear that the method (36) is

O(At) in time due to the linearization in 6u,,, ; this does not, however, preclude use of

a O((At) q) scheme in space-time in the Lagrangian-decoupling procedure. In essence,

the former relates to the temporal accuracy of the integration along the characteristic,

whereas the latter relates to the spatial accuracy of the integration; note that when a

steady-state is achieved, the O(At) errors vanish, whereas the spatial errors remain. For

time dependent problems, variants of (36) which are more accurate in time are readily

22
generated by appropriate high-order linearizations [25]. The spatial errors due to the

Stokesdiscretization are standard [1,14,35,36],and will not be discussedfurther.

23
5. Spectral Element Implementation

In this section we consider the spectral element method as a particular example of

the general variational spatial discretization described in Section 3.1. We give numerical

results which demonstrate the convergence rate of the spectral-element-based Lagrangiandecoupling method for the convection-diffusion of a passive scalar. We also demonstrate

the viability of the iacthod for the full Navier-Stokes equations, ill particular as regards

increased stability and associated computational savings.

5.1 Description of Spectral Element Spatial Discretization

The spectral element method is a high-order variational method for the spatial discretization of partial differential equations [14,15,32]. As in (19), tile computational domain

is first broken up into K subdomains (spectral elements) O k , k = {1, ..., K}. For reasons

of efficiency (i.e., tile availability of tensor product sum factorization [34]) the elements

12_ are taken to be line segments in 1121, quadrilaterals in /R 2, and hexahedra (bricks) in

_3, defined with respect to the reference volumes I_* =] - 1, 1[J in J_ttJ by an invertible

mapping F ,k (20). Although geometrically nonconforming spectral element methods have

been developed [37,38 l, we restrict ourselves here to the conforming case.

Within each reference volume _* the discrete polynomial subspace/P* is taken to be

the set of all polynomials of degree < N in each spatial direction [14]; the spectral element

spatial discretization call be characterized by the discretization pair h = (K, N). Although

the particular choice of basis points _'* does not effect the error estimates, it greatly effects

the conditioning and sparsity of the resulting set of algebraic equations, and is critical for

the efficiency of parallel iterative solution procedures [13,29]. In order to take advantage

of efficient sum-factorization techniques [34], the N h = (N + 1) a basis points _+'_ for _*

are taken to be the d-tensor product of the N + 1 Gauss-Lobatto Legendre points [39];

the corresponding Lagrangian interpolants Qi,_ of (21) are therefore the product of d N th

24
order one-dimensionalLagrangian interpolants I14].

In the spectral element method the discrete inner products defined in Section 3.1b

are based on tensor-product Gauss-Lobatto Legendre integration, that is, _ is equivalent

to Gaussian quadrature. TILe quadrature points are the same as the basis points; with

N + 1 quadrature points we can integrate exactly polynomials of degree < 2N - 1. This

particular choice of numerical integration (denoted "consistent" integration) can be shown

to be sufficient and optimal for most problems [14,33]. In particular, consistent integration is sufficient for complex geometries, an important consideration when considering the

invariance equation (32), in which tile base geometry may deform significantly.

TILe spectral element method combines the geometric flexibilty of a low-order method

with the rapid convergence rate associated with spectral techniques. Considering only

spatial errors, it can be shown that tile spectral element solution to the convection-diffusion

equation (4) or to tile full Navier-Stokes equations (1) converges spectrally fast to the exact

solution for K fixed, N _ oo, with exponential convergence obtaining for locally analytic

geometry, data and solution. This rapid convergence rate derives from the good stability

and approximation properties of the polynomial spaces /P_, and the accuracy associated

with Gauss-Lobatto Legendre quadrature and interpolation [14,39]. Moreover, the discrete

solution suffers from minimal numerical dispersion and diffusion, a fact which is important

in the solution of tile invariance equation (I6).

Having defined our spectral element spatial discretization, we need to solve the set of

algebraic equations (31)-(34) corresponding to tile convection-diffusion of a passive scalar.

For the implicitly treated symmetric positive-definite elliptic kernel (diffusion) in (31) we

employ preconditioned conjugate gradient iteration or rapidly-convergent intra-element

multigrid techniques [16,17]; for the Stokes operator of (36) we employ splitting schemes

[11,12] or preconditioned Uzawa methods [4,14] to reduce the procedure to elliptic subproblems. The key to minimal work per iteration is the use of tensor product elements,

spaces, bases, and (consistent) quadratures, allowing for sum-factorization and efficient

25
matrix-vector product evaluations.

We note that the solution of (31) is stable for all time steps At (Q < 6). However, in

order to evaluate ¢ from previous timesteps at the foot of the relevant characteristics, we

need to solve a pure convection problezn fore and l_L.The details of this evaluation are given

in (33), in which we have chosen to use an explicit third-order Adams-Bashforth integration

scheme; the associated spectral element Courant condition is At,, < O(1/KlldN2) [15].

It is important to note that our approach to solving the invariance equation (16) preserves

tile underlying tensor-product structure of the spatial discretization, and hence no costly

re-interpolation is necessary. Indeed, if referencing were done only to the base geometry

X_I;', interpolation would require O(KN 2"_) operations, compared to O(KN d+l) for our

tensor-product scheme [34].

5.2 One-dimensional results for a passive scalar

In this section we consider tile solution of tile convection-diffusion equation (4) in one

space dimension. In the first test problem we choose 12 =]0,1[, v = 0.01, u = 1, with

qb(x,t = O) - sin(2_rx), and periodic boundary conditions. The exact solution is given by

¢(z,t) = e-_(2")_tsin(27r(x -ut)). First, we demonstrate the temporal accuracy of the

new Lagrangian-decoupling method. Figure 1 shows the L °°-error at a fixed time t = 2 as

a function of the time step At for a fixed spatial discretization K = 2 and N = 15; this

spatial discretization assures that the temporal error is always dominating. The results

clearly demonstrate that the temporal error is O((At) Q), where Q is the order of the

backward differentiation formula in (31). Since As is typically much smaller than At,

the temporal error resulting from solving (33) will, in general, be much smaller than the

temporal error resulting from solving (31), at least for Q < 3. We note that for pure

convection (1., = 0), the only temporal error is due to the integration of the invariance

equation (32)-(34); for this smooth problem, the accuracy improves as v decreases (in

contrast to "mesh-upwinding"). Next, we demonstrate the spatial error of the scheme. In

26
Figure 2 we plot tile L _ -error at a fixed time t = 2 as a function of tile polynomial degree

N, for fixed K = 2 elements, with Q = 3, At = 0.01, and AtAs = 100. As expected, we

obtain exponential convergence as long as N < 10; for N > 10 the temporal error becomes

dominating given our particular choice of Q and At.

As our second one-dimensional test problem we consider the steady state solution

to (4) in fl =]0, 1[, _, = 0.1, u = 1, and essential boundary conditions ¢(x = 0, t) = 0,

"Iv-1 which has a

¢(x = 1,t) = 1. The exact solution as t --* oo is given by 4) - ;,'7";':-T,

boundary layer of thickness _, at z = 1. The numerical steady solution is obtained by

integrating the unsteady convection-diffusion equation to a time t = 10, starting from an

initial condition _b(x, t = 0) = x. Figure 3 shows the L =-error as a function of the time

step At for Q = 1,2,3,4, for a fixed spatial discretization K = 2 and N = 11. Note that

in this case tile time step At has to be smaUer than O(_) before the expected convergence

rate is observed. Also note that, unlike the standard Eulerian convection scheme, ill which

tile steady state error is purely spatial, the error in our new scheme is contaminated by the

temporal discretization due to the transformation of space into time. In Figure 4 we plot

tile L _ -error as a function of tile polynomial degree N for fixed K = 2, Q = 3, At = 10- 3,

and AtAs = 10. The results demonstrate exponential convergence, at least for N < 10;

for N > 10 the temporal error again becomes dominating.

We note that it is a drawback to tile current method that At must be the same at all

points on the domain, as this does not allow for local refinement near important spatial

structures. Fortunately, the flow "through the wall" (u,_ fir,, _ 0) associated with the

normal boundary layer of Fig. 3 is not typical; if u,,, h_ = 0, the characteristics scheme is,

in some sense, self-adaptive, as Ax (x f u dr'.

5.3 Two-dimensional results for a passive scalar

The first two-dimensional test problem we consider is pure convection (t_ = 0) of

a passive scalar _ in f_ =]0,1[ 2. The boundary conditions are periodic conditions at

27
Xl = 0 and Xl = 1, and homogeneous Neumann conditions at x2 = 0 and x2 = 1. Tile

prescribed solenoidal velocity field is Ul = 1, u2 = 0, and the exact solution is given by

¢(Xl, x2, t) = sin(27r(Xl - ult))cos(27rx2). In Figure 5 we plot the error in the H I seminorm at a fixed time t = 5 as a function of the polynomial degree N for fixed K = 4

spectral elements, Q = 1 and At = 0.1. The results clearly demonstrates that exponential

convergence is achieved due to the smooth nature of the solution. Note that for this pure

convection problem (for which our new scheme is not of practical interest, as )4;_xt_0 is

effectively zero), the only temporal error is the error due to integrating the invariance

equation (32), and thus Q = 1, At = .1 is sufficient. The subtimestep As is sufficiently

small here that the spatial error is dominating.

As the second test problem we consider the steady state solution to the convectiondiffusion equation (4) in the two-dimensional domain 12 : Xl e] - 1/2, 1/2[, x2 El0, 1[. The

given velocity field is tile stagnation potential-flow Ul = Xl, u2 = -x2. With boundary

conditions ¢(Xl, x2 = 0, t) = 1 and ¢(Xh x2 = 1, t) = 0 the steady state solution depends

only on x2, and is given as ¢ = 1 - f[2 ,-_,/,v_f

f_ ,_t,/2vd _ In Figure 6 we plot the error in the

H 1 semi-norm at a fixed (effectively infinite) time t = 5 as a function of the time step At

for Q = 1,2,3 for a fixed spatial discretization consisting of K = 4 spectral elements of

polynomial degree N = 10. Note that we impose the exact solution as essential boundary

conditions on the whole domain boundary. The expected convergence rate O((At) Q) is

clearly demonstrated. For Q = 1 we repeat the experiment, but now imposing essential

boundary conditions only along x2 = 0 and x2 = 1, with homogeneous natural (Neumann)

boundary conditions imposed along Xl = -1/2 and Xl = 1/2; the results demonstrate that

the scheme works equally well for outflow boundary conditions. In Figure 7 we plot the

error in the semi-norm at a fixed time t = 5 as a function of the polynomial degree N for

fixed K = 4, Q = 2 or 3, At = 10-3, and 5 _< AtAs _< 10. Again, the results verify the

spectral accuracy of the scheme.

Our stagnation flow results are different from the one-dimensional normal boundary

28
layer in two aspects. First, tile boundary layer scale only as t/112 here, as opposed to

the atypical _, for normal layers. Second, and more importantly, we note that there is no

evidence in Fig. 6 of a boundary layer threshold effect in At, as u,_ fi,,, = 0 and the effective

Ax near the wall is therefore small. In most real viscous flows, in which no-slip is applied,

u2 would in fact vary as -x 2, yielding even better accuracy than the slip u2 = -x2 case

studied here.

5.4 Two-dimensional Navier-Stokes results

We now return to the unsteady Navier-Stokes system, which was the original motivation for the new scheme due to to the potential savings in computational cost. As a

test problem we consider flow in a xt-periodic grooved channel, in which the flow is driven

by a constant pressure gradient in the xl-direction. These first results are intended only

to demonstrate the stability of the method, and to hint at the potential computational

savings possible.

The tests described below correspond to the geometry shown in Fig. 8a for a Reynolds

number R = _ = 25, where f is the constant force (pressure gradient normalized by

is 3

density), and H is the half-channel width. Only steady state results are presented. In

all cases the spectral element mesh (X I;') is based on the K = 6 elements shown in Fig.

8a, with N = 6 in each element. In order to determine the accuracy and efficiency of

the characteristic-decoupling scheme we first solve the problem by using a standard semiimplicit "Eulerian" scheme of the type (2), in which tile convective term is treated with

a third-order Adams-Bashforth scheme. A steady-state solution (t = 4) is reached after

834 timesteps, with the solution shown in Figures 8b and 8c. The maximum velocity,

normalized by H2f/t/, is found to be U,_, = 1.0.

Next, we repeat this test using the new Lagrangian-decoupling method. We use the

first order Q = 1 scheme (only steady-state results are examined), with a timestep At

which is 50 times larger than that for the Eulerian scheme (note that for tlle Eulerian

29
schemethe Navier-Stokes timestep At is chosen from stability considerations, At = At,, ;

for the characteristic scheme only As in (33) is limited by At,,). The final time t = 4 is

now reached in 16 timesteps, giving the results shown in Figures 9a and 9b; in terms of

computational cost, tile new scheme obtained the t = 4 solution 16 times faster titan the

standard approach. As regards accuracy, tlle two solutions look qualitatively the same,

however, the maximum non-dimensional velocity in tile characteristic case is U_,,, = 0.88,

indicating a 10% error with respect to tile gulerian scheme. Although this discrepancy

may appear to be due to tile severe distortion of the mesh, shown in Fig. 9c as a plot

of X _;;At/A', tests on plane Poiseuille flow (in which only distortion effects are present)

indicate that this degree of distortion is readily handled by the h = (K = 6, N = 6)

approximation. Reduction of At or increasing Q reveals that the error is, in fact temporal;

systematic numerical convergence results for Navier-Stokes are given in [25].

The results of the these first Navier-Stokes calculations are by no means conclusive;

tests to determine computational savings at fixed accuracy are currently underway I251.

Of particular importance is understanding the dependence of At',,, and _'V"_ on temporal

order Q, the solution-induced deformation (in particular in boundary layers and near

singularities), the spectral (or finite) element discretization, incomplete iteration, and the

Stokes solver efficiency. It is clear, nevertheless, that the technique holds great promise in

striking an efficient compromise between temporal stability and equation structure.

ACKNOWLEDGEMENTS:

This work was supported by the ONR and DARPA under Grant N00014-89-J-1610,

by the ONR under Contract N00014-88-K-0188, by the NSF under Grant 8806925-ASC,

and by NASA under NAS1-18605.

30
References

[1] F. Brezzi, On the existence, uniqueness and approximation of saddle point problems

arising from Lagrangian multipliers, R.A.I.R.O. Numer. Anal. 8 (197_t) 129-151.

[2] M.O. Bristeau, R. Glowinski and J. Periaux, Numerical methods for the Navier-Stokes

equations. Applications to the simulation of compressible and incompressible viscous flows,

Computer Physics Report, to appear.

[3] J. Cahouet and J.P. Chabard, Some fast three-dimensional finite element solvers for

the generalized Stokes problem, Int. J. Numer. Meth. in Fluids 8 (1988) 869-895.

[4] Y. Maday, D. Meiron, Patera, A.T. and R0nquist, E.M., Analysis of iterative methods

for the steady and unsteady Stokes problem: Application to spectral element discretizations, J. Comput. Phys., submitted.

[5] A. Brooks and T.J.R. Hughes, Streamline upwind/Petrov Galerkin formulation for

convection-dominated flows witlt particular empltasis on the incompressible Navier-Stokes

equations, Comp. Meth. Appl. Mech. Eng. 32 (1982)199.

[6] C.W. Gear, Numerical Initial Value Problems in Ordinary Differential Equations (Prentice -Hall, Englewood Cliffs, 1971).

[7] Y. Yamaguchi, C.J. Chang and R.A. Brown, Multiple bouyancy-driven flows in a vertical

cylinder heated from below, Phil. Trans, Roy. Soc. London A (1984) 312-519.

[8] A.J. Chorin, Numerical solution of incompressible flow problems, in: J.M. Ortega and

W.C. Reinboldt, eds., Studies in Numerical Analysis 2, (SIAM, 1970) 64-71.

[9] R. Temam, Navier-Stokes Equations. Theory and Numerical Analysis (North-Holland,

Amsterdam, 1977).

[10] S.A. Orszag and L.C. Kells, Transition to turbulence in plane Poiseuille flow and plane

Couette flow, J. Fluid Mech. 96 (1980) 159.

31
[11] K.Z. Korczak and A.T. Patera, An isoparametric spectral element method for solution

of the Navier-Stokes equations in complex geometry, J. Comput. Phys. 62 (1986) 361-382.

[12] P.F. Fischer, E.M. RCnquist and A.T. Patera, Parallel spectral element methods for

viscous flows, ill: Proceedings of the Workshop on Methods and Algorithms for Partial

Differential Equations on Advanced Processors, Austin, Texas, 1988, to appear.

[13] P.F, Fischer, Parallel spectral element methods for the incompressible Navier-Stokes

equations (Ph.D. Thesis, Massachusetts Institute of Technology, 1989).

[14] Y. Maday and A.T. Patera, Spectral element methods for the Navier-Stokes equations,

in: A.K.Noor and J.T. Oden, eds., State of tile Art Surveys in Computational Mechanics

(ASME, New York, 1989) 71-144.

I15]E.M. ROnquist, Optimal spectral element methods for the unsteady three-dimensional

Navier-Stokes equations, (Ph.D. Thesis, Massachusetts Institute of Technology, 1988).

[16] E.M. Rcnquist and A.T. Patera, Spectral element multigrid. I. Formulation and

numerical results, J. Sci. Comput. 2(4) (1987) 389-402.

[17] Y. Maday and M. Munoz, Spectral element multigrid. II. Theoretical justification, J.

Sci. Comput., to appear.

[18] O. Pironneau, On tile transport-diffusion algorithm and its applications to tile NavierStokes equations, Numer. Math. 38 (1982) 309-:]32.

[19] B. Ramaswamy, Numerical simulation of unsteady viscous free surface flow, J. Comput.

Phys., submitted.

[20] P. Gresho, S. Chan, R. Lee and C. Upson, A modified F.E.M. for solving the timedependent incompressible Navier-Stokes equations, Int. J. Num. Meth. F1. 4 (1984)

557-598 (Part 1) and 619-640 (Part 2).

[21] P.3. Roache, Computational Fluid Dynamics (Hermosa Publishers, Albuquerque,

1982).

32
[22]K. Kuwahara, Developmentof high-Reynolds-number-flowcomputations, in: K. Kuwahara, R. Mendezand S.A. Orszag,eds., Lecture Notes in Engineering 24 (Springer-Verlag,

1985)28-40.

I231 G. Strang and d. Fix, An Analysis of the Finite Element Method (Prentice-Hal/,

Englewood Cliffs, 1973).

[24] C. Canuto, M. Hussaini, A. Quarteroni and T. Zang, Spectral Methods in Fluid

Dynamics (Springer-Verlag, 1987).

[25] Y. Maday, A.T. Patera and E.M. R0nquist, Analysis of a high-order Lagrangiandecoupling method for the incompressibl e Navier-Stokes equations, ICASE Report 89-53

(in preparation).

[26] R.A. Adams, Sobolev Spaces (Academic Press, 1975).

[27] L.W. Ho, A spectral element stress formulation of the Navier-Stokes equations (Ph.D.

Thesis, Massachusetts Institute of Technology, I989).

[28] R. Glowinski, G. Golub, G. Meurant and J. Periaux, J., eds., Proceedings of tile

First International Conference on Domain Decomposition Methods for Partial Differential

Equations (SIAM, Philadelphia, 1988).

[29] P.F. Fischer and A.T. Patera, Parallel spectral element solution of tile Stokes problem,

J. Comput. Phys, to appear.

[30] I. Babuska and M.R. Dorr, Error estimates for tile combined h- and p-version of the

finite element method, Numer. Math. 37 (1981) 257.

[31] B.A. Szabo, Mesh design for the p-version of the finite element method, Comp. Meth.

in Appl. Math. Engr. 55 (1986) 181-197.

[32] A.T. Patera, A spectral element method for Fluid Dynamics: Laminar flow in a channel

expansion, J. Comput. Phys. 54 (1984) 468-488.

33
[33]Y. Maday andE.M. Rcnquist, Optimal error analysis of spectral methods with emphasis on noneonstant coefficients and deformed geometries, in Proe. Int. Conf. on Spectral

and High-Order Meth. for Partial Differential Equations, Italy, June 1989, to appear.

[34] S.A. Orszag, Spectral methods for problems in complex geometries, J. Comput. Phys.,

37 (1980) 37-70.

[35] V. Girault and P.A. Raviart, P.A., Finite Element Approximation of the Navier-Stokes

Equations (Springer-Verlag, Berlin, 1986).

[36] Y. Maday, A.T. Patera and E.M. R0nquist, A well-posed optimal spectral element

approximation for the Stokes problem, SIAM J. Numer. Anal., to appear.

[37] C. Bernardi, Y. Maday and A.T. Patera, A new non-conforming approach to domain

decomposition: the mortar element method, in: H. Brezis and J.L. Lions, eds., Nonlinear

Partial Differential Equations and their Applications, to appear.

[38] Y. Maday, C. Mavriplis and A.T. Patera, Non-conforming mortar element methods:

application to spectral discretizations, in: T. Chan, ed., Proceedings of the 7th International Conference on Domain Decomposition Techniques for Partial Differential Equations

(SIAM, Philadelphia, 1988).

[39] P.J. Davis and P. Rabonowitz, Methods of Numerical Integration (Academic Press,

1985).

34
Figure Captions

Figure 1. A plot of tile discretization error II ¢ - Ch [1_- at a fixed time t = 2 as a

function of the timestep At for a fixed (accurate) spatial discretization, when solving the

convection-diffusion equation (4) using the new Lagrangian-decoupling method. The exact

solution is given as ¢(z,t) = e-V(2")_'sin(2rr(z -ut) in the (periodic) domain f/ = [0,1],

with v = 0.01 and u = 1. The results show the (dominating) temporal error of the Q = 1

(C)), Q = 2 (D), and Q = 3 (A) backward differentiation schemes.

Figure 2. A plot of the discretization error II ¢-¢h I1_- at a fixed time t = 2 as a function

of the polynomial degree N in the It" = 2 fixed elements, when solving tile convectiondiffusion equation (4) using tile new Lagrangian-decoupling method. The timestep is set

small to minimize temporal pollution. TILe exact one-dimensional periodic solution is given

as ¢(x,t) = e-v(e")_'sin(27r(x -ut) for x E _ = [0, 1], with v = 0.01 and u = 1. Spectral

convergence is achieved for N < 10; for N > 10 the (fixed) temporal error becomes

dominating.

Figure 3. A plot of the steady-state discretization error II ¢- Ch IlL- as a function of tile

timestep At for a fixed (accurate) spatial discretization, when solving the convectiondiffusion equation (4) using the new Lagrangian-decoupling method. The exact onedimensional steady-state boundary layer solution is given by ¢ = _ for x E _ = [0, 1]

(v = 0.1). The results show the (dominating) temporal error of the Q = 1 (O), Q = 2 (e),

Q -- 3 (A), and Q = 4 (A) backward differentiation schemes.

Figure 4. A plot of the steady-state discretization error ]] ¢ - Ch ILL*"as a function of

the polynomial degree N in K -- 2 fixed elements, when solving the convection-diffusion

equation (4) using the new Lagrangian-decoupling method. The timestep is set small

to minimize temporal pollution. The exact one-dimensional steady-state boundary layer

,-'_-1 for x E fi -- [0,1] (v = 0.1). Spectral convergence is solution is given by ¢ = ;,-7";'=-_1

achieved for N < 10; for N > 10 tile (fixed) temporal error becomes dominating.

35
Figure 5. A plot of tile discretization error in tile H 1 semi-norm, i ¢ - Ch Jl, at a fixed

time t = 5 as a function of the polynomial degree N in K = 4 fixed elements, when solving

the convection equation (4) (v = 0) using the new Lagrangian-decoupling method. The

exact two-dimensional solution is given by ¢(Xl, X2, t) = sin(2_r(xl -ult)cos(27rx2) for

(xl, x2) C fi = [0, 1] 2, with ul = 1, u2 = 0. Periodic boundary conditions are imposed at

xt = 0 and Xl = 1, with lmmogeneous Neumann conditions imposed at x2 = 0 and x2 = 1.

Spectral accuracy is achieved.

Figure 6. A plot of the steady-state discretization error in the H I semi-norm, I ¢ -

Ch I1, as a function of the timestep At for a fixed (accurate) spatial discretization, when

solving the two-dimensional convection-diffusion equation (4) using the new Lagrangiandecoupling method. The exact steady-state solution is given as ¢ = 1 - 'f_" '-c_/_'ae for

:rt E [-1/2,1/2], x2 E [0, 1], for an imposed stagnation flow ut = xl, u2 = -x2. With

imposed Dirichlet boundary conditions for ¢, the results show the (dominating) temporal

error for backward differentiation schemes of order Q = 1 (O), Q = 2 (0), and Q = 3 (A).

For comparison we also plot the the error for the case Q = 1 with homogeneous Neumann

conditions imposed along x t : +1/2 (outflow) (.).

Figure 7. A plot of the steady-state discretization error in the H I semi-norm, I ¢-¢J, Jl, as

a function of the polynomial degree N in K = 4 fixed spectral elements, when solving the

two-dimensional convection-diffusion equation (4) using the new Lagrangian-decoupling

method. The timestep is set small to minimize temporal pollution. The exact solution is

given by ¢ = 1 - f_'_ ,-¢_'_'d_ for Xl e [-1/2,1/2], z2 E [0, 1], for an imposed stagnation

fo t,,-ts/3vd_

flow Ul = Xl, u2 = -x2. Spectral accuracy is achieved.

gradient in the Xl direction.

is reached at a time t = 4.

Figure 8. Two-dimensional Navier-Stokes solution in a xl-periodic grooved channel at

a Reynolds number R = _ 25, in which the flow is driven by a constant pressure

It 3 ---

Starting from rest, an approximately steady state solution

Here the numerical solution is obtained using a standard

36
(Eulerian) semi-implicit schemedescribedby equation (2), basedon tile spectral element

mesh shown in Fig. 8a. In Fig. 8b and 8c we plot the solution in terms of streamlines and

velocity vectors, respectively.

Figure 9. Two-dimensional Navier-Stokes solution in a zl-periodic grooved channel at

a Reynolds number R = _ = 25, in which the flow is driven by a constant pressure

gradient in the zl direction. Starting from rest, an approximately steady state solution is

reached at a time t = 4. Here the numeric',d solution is obtained using tile new Lagrangiandecoupling method. In Fig. 9b and 9c we plot the solution in terms of streamlines and

velocity vectors, respectively, and in Fig. 9c we indicate the mesh distortion that occurs in

the invariance-equation integration. Note that although the characteristic method allows

for large At, and hence rapid convergence to the steady-state, the method does not reduce

the parametric stiffness associated with long timescales as R _ co.

37
I0

°

IIIIII1

IIIIII 'T:

|O "I

10-2

10-3

8 I0 -_ ,4I O. s

_.11

0 "e

I0 "']

10 -8

10 -9

lO-Z

0OA

0

El

00

z_

0

E] &

:: I_lI

lO-I

_t

00&

00A

I _II I I 1 I

10o

-

FIGURE1

38
0

0

0

S I0 -'i

I lO-S

0

0

I0 -6

0

0

0 0 0

I 0 -9

2 3 4 5 B 7

N

8 9 L0 IL 12 13

FIGURE 2

39
-100

I0-!

I0-s

0

I | I I I I I i

0

o

A

A

A

0

A

0

i i i IIlr]

©

0

&

&

I _ 1 l 11|

lO-i

0

&

i ; l t i i_

©

l I I l I I J

10 o

FIGURE 3

40
10 °

I0 °|

10 .2

10 -3

8

,._ 10 -'_

I

-_ i0 -s

10-a

10 -7

10 -8

10 -9

0

_

0

3

0

4

0

5

0

6

0

7

N

0

8

0

9

0

I0

0

11 12 13

FIGURE 4

41
10 °

i0-1

I0-2

I I I

I

6

l I I

?

N

FIGURE 5

&

I

8

&

l

9

!

10

i

1

42
I0

° I I ! | I I I I I I I I i I I II I I I I I I I II I I I I I 1 I I

s,-i

10 -|

10 -2 __--

10-s

•_ 0_4 I

lO-S_

10-tl

10 -_

F

10-e i

10-4

®O I 1 i I iill

10-3

@OA i , ! ® ® @ O 0 D i i i i i,l i

10 -z

At

FIGURE6

®

@

DA

[]A i I I I i Ill

10 -I

®

El Ai I I I I

e

I

I0o

43
1"00 = ,

lO-I

lO-Z

10-3

+=,.

i0-'i__-

#

_--- lO-S

i0-8

i0-?

lO-a

10-9 - ,

I

I l I i I

O

! • 1 l !

2 3 q

0

I ' I ' I ' I ' I

0

0

I , l i I

5 6 '7

N

O

O

I i I

8 9

O

O

1 J 1 l_

IO II

FIGURE 7

44
12H "I

6H

FIGURE 8 (a)

FIGURE 8 (b)

FIGURE 8 (c)

45
FIGURE 9 (a)

FIGURE 9 (b)

FIGURE 9 (c)

46
_,_,._,c.._ Report Documentation Page

1. Report No 2. Government Accession No. 3. Recipient's Catalog No.

NASA CR- 181903

ICASE Report No. 89-57

4. Title"and Subtitle ....

A HIGH-ORDER LAGRANGIAN-DECOUPLING METHOD FOR THE

INCOMPRESSIBLE NAVIER-STOKES EQUATIONS

7. Author(s)

L. Ho, Y. Maday, A. Patera, E. R_nquist

9. Performing Organization Name and Address

Institute for Computer Applications in Science and

Engineering

Mail Stop 132C, NASA Langley Research Center

Hampton_ VA 23665-5225

12. Sponsoring Agency Name and Addre_

National Aeronautics and Space Administration

Langley Research Center

Hampton, VA 23665-5225

5. Report Date

September 1989

6. Performing Organization Code

8. Performing Organization Report No.

89-57

10_ Work Unit No.

505-90-21-01

11, Contract or Grant No.

NASI-18605

13. Type of Report and I_eriod Covered

Contractor Report

14. Sponsoring Agency Code

15. Supplementaw Notes

Langley Technical Monitor:

Richard W. Barnwell

Proceedings of ICOSAMOM Meeting 1989

Final Report

16, Abstract

In this :paper we present a high-order Lagranglan-decoupling method for the

unsteady convection-diffusion and incompressible Navier-Stokes equations. The

method is based upon: Lagrangian variational forms that reduce the convectiondiffusion equation to a symmetric initial value problem; implicit hlgh-order

backward-differentation finite-difference schemes for integration along characteristics; finite element or spectral element spatial discretizations; mesh-invariance procedures and high-order explicit time-stepplng schemes for deducing

function values at convected space-time points. The method improves upon previous

finite element characteristic methods through the systematic and efficient

extension to high order accuracy, and the introduction of a simple structurepreserving characteristic-foot calculation procedure which is readily implemented

on modern architectures. The new method is significantly more efficient than

explicit-convection schemes for the Navier-Stokes equations due to the decoupling

of the convection and Stokes operators and the attendant increase in temporal

stability. Numerous numerical examples are given for the convection-diffusion

and Navier-Stokes equations for the particular case of a spectral element spatial

discretizaiton.

17. Key Wor_ (Sugg_tedbyAu_or(s))

characteristic, imcompressible NavierStokes, spectral methods

19. Security Classif. (of this report)

18. Distribution S'tatement

64 - Numerical Analysis

Unclassified - Unlimited

_. Secur_ C_ssif. (of this _ga) 21.'No. of _ges 22. Price

A0 3

NASA-Langley, 198g

Unclassified

NASA FORM 1626 OCT

Unclassifie4 47