fh_2000_03.pdf (3.705Mb)

Emneord - norsk:

l . Adveksjons skjema 2. Ustrukturert gitter 3. Havmodeller

Norges Forskningsråd Nordnesgaten 50 Postboks 1870 58 17 Bergen

Tlf.: 55 23 85 00 Faks: 55 23 85 31 Forskningsstasjonen Austevoll

Flodevigen havbruksstasjon havbruksstasjon 48 17 His 5392 Storebo 5984 Matredal Tlf.: 37 05 90 00 Tif.: 56 18 03 42 Tlf.: 56 36 60 40 Faks: 37 05 90 01 Faks: 56 18 03 98 Faks: 56 36 61 43 Rapport:

FISKEN OG HAVET NR. 3

-

2000

Emneord - engelsk:

1. Advection schemes 2. Unstructured grids 3. Ocean models

Tittel:

ADVECTION SCHEMES FOR COASTAL OCEAN MODELS ON UNSTRUCTURED TRIANGULAR M E S E S

Forfatter(e):

W. Paul Budgell og Morten D. Skogen

&A

Se 'onsleder

Senter:

Senter for Marint Miljø

Seksjon:

Fysisk Oseanografi

Antall sider, vedlegg inkl.:

6 1

Dato:

26.10.2000

Sammendrag:

Flere ulike adveksjonsskjema til bruk for ustrukturert gitter er testet med det formål å vurdere deres egnethet til bmk i numeriske havmodeller. De ulike skjemaene er testet på flere ulike idealiserte test problemer i både en og to dimensjoner. Til slutt er et av

skjemaene, endelig element fluks korrigert transport, implementert i en tre-dimensjonal endelig element havmodell, og testet både på et idealisert og et realistisk tilfelle hvor ferskvann fra en elv mplter saltere vann utenfor kysten.

Advection Schemes for Coastal Ocean Models on Unstructured Triangular Meshes

Norges Forskningsråd Project 13 3 109143 1

W. Paul Budgell Project Leader

and

Morten D. Skogen

Institute of Marine Research, Bergen 14 February, 2000

Summary

A variety of advection schemes suitable for use on unstructured triangular grids have been investigated for possible use in coastal ocean models based on the results of a literature survey. Idealised test cases were conducted in one and two dimensions to assess performance of the candidate schemes. One scheme, the finite element flux- corrected transport algorithm was implemented in a three-dimensional finite element ocean model and tested on an idealised and an actual estuarinelriver plume application. An article describing the results of the study is currently in preparation for submission to a refereed journal.

l. Introduction

This report describes the results of a 6-month pre-project under the BeMatA (Beregningsorientert Matematikk i Anvendelser) programrne of the Norwegian Research Council. The principal objective of the project was to provide an advection scheme suitable for use by three-dimensional coastal ocean models based on unstructured triangular grids. The sub-goals were:

to review the computational fluid dynamics literature for suitable schemes to implement and test a limited number of suitable schemes on l-D and 2 - 0 problems

to select and implement the most suitable scheme on 3-D idealised and realistic test cases

to prepare an article on the results of the study for submission to an appropri ate journal

The remainder of this report is structured for presentation of the results associated with each of the sub-goals.

2. Literature Review

It quickly became apparent in the course of conducting the literature review for this study that an enormous body of research has been conducted in the area of advection on unstructured grids over the past twenty years. Very little of this work has penetrated the field of ocean modelling, however. We consider one of the most useful contributions of this project to be the transfer of knowledge related to advection on unstructured grids gained in other branches of computational fluid mechanics (e.g., aerodynamics, groundwater hydrology, petroleum reservoirs) to ocean modelling.

Several excellent books have appeared over the past ten years on the subject of numerical treatment of advection and hyperbolic problems. Among these, Hirsch (1990), LeVeque (1992) and Toro (1999), while adressing finite difference approaches, provide valuable background on the issues involved with numerical treatment of advection problems. Finlayson (1992) and Morten (1996) include excellent overviews of finite element approaches using unstructured triangular grids.

Godlewski and Raviart (1996) and Kroner (1997) provide comprehensive descriptions of finite volume methods, including those suitable for use on triangular unstructured grids. Vreugdenhil and Koren (1993) is an extensive comparison of a wide variety of schemes for advection-diffusion problems, including finite element and finite volume approaches. These texts provided an excellent background for the extensive literature search and pursuit of specific algorithms.

Based on the literature review, the following schemes were selected for inclusion in the study:

i. Galerkin Finite Element (consistent mass matrix) ii. Lumped, Leap-frog Galerkin Finite Element

. . .

111. Streamline Upwind Petrov-Galerlun

iv. Taylor-Galerlun

v. Flux-Corrected Transport Finite Element vi. Runge-Kutta Discontinuous Galerkin vii. Cell-based Finite Volume

viii. Control Volume Finite Element Eulerian-Lagrangian Method

i) Galerkin Firzite Element

The "classic" finite element scheme is the standard Galerlun approach (see Finlayson (1992) or Morten (1996) for a derivation). If the two-dimensional advection equation is expressed as

with Dirichlet boundary conditions on the inflow part

(r,)

of the boundary

(r)

c ( x , t ) ⁼ cb ( x , t ) , X €

r1

(2) where c is concentration, u is the advecting velocity, t is time and c,, is the specified boundary forcing, then the weak formulation of the problem is

where 52 is the spatia1 domain and v i s a sufficiently smooth weighting function for integration over the domain. In the standard Galerkin approach, (3) is approximated by finite dimensional subspaces, defined as

{gi

(x)} which form a basis of the space of functions whose derivatives are square integrable. The approximation to the solution,

c, , may be defined as

where Tsatisfies the inhomogeneous boundary condition (2). If we choose a weighting function v, to be the same as the trial function

J

{ v ,

+

^(u V

)

v

n

in the subspace spanned by the basis functions ^{@l, @?,} . . .,Gll so that each basis function

@i is substituted separately for v,, such that (5) becomes

then we have the standard Galerkin approach. The finite element method provides a means of creating the basis functions @i and the boundary function

q .

The domain L2is subdvided into elements and within each element nodal points and a polynomial approximation are chosen such that:

@ ( x , ) = 6 , for all nodal points j

@i ( X , ) is an element polynomial

The boundary function

c,

is represented by a linear combination on basis functions on the boundary

T,.

When the specific finite element polynomials are implemented in (6), a system of ordinary differential equations results:

wherec is the solution vector, M is known as the mass matrix, K is known as the stiffness matrix a n d F is the right-hand side vector. In this study, linear basis functions on triangular elements are employed. Time integration of (8) is accomplished by a centred, implicit (Crank-Nicholson) second-order scheme.

(ii) Lumped, Leap-Frog Galerkin Finite Element

This scheme is similar to (i) above, but differs in the treatment of Mand the time integration of dcldt

.

Maintaining the full, or consistent, mass matrix in (8) can be computationally expensive, since a large number of equations must be solved simultaneously through either direct or iterative means. Thus, there is a strong temptation to perform nodal quadrature, or mass lumping, to reduce M to a diagonal matrix through a sumrnation of the contributions from surrounding nodes at each node. The lumped mass treatment, when combined with leap-frog time discretisation, provides second-order accuracy in time while remaining computationally inexpensive enough for use in three-dimensional models. The leap-frog time stepping is combined with an Asselin (1972) filter in time to remove computational model oscillations.

iii) Streamline Upwind Petrov-Galerkiiz

The streamlined upwind Petrov-Galerkin (SUPG) scheme (Brooks and Hughes, 1982) applies an upwind weighting to obtain stable solutions in advection-dominated flows.

A perturbed weighting function

v7,

is employed, where

Vh

^=V,? +Ph (9)

and

replaces ( 5 ) , where v, is the Galerkin trial function and p, is a perturbation parameter introduced by Brooks and Hughes (1982) that produces upwind weighting, but only in the direction of the flow in order to minimise crosswind diffusion:

-

u*V V,,

P,, = D- liull

L$,

determines the strength of the upwinding. Several parameterisations of p, and D, are described in the literature. In this study we use ( l l ) and define

6,

as:

where Ax is the local width of the element and

In this study, since we are concerned with advection processes, the diffusion coefficient D is set to 1x10-lo, so that D becomes Ax/2 and p,, adds diffusion along streamlines by the amount incurred by first-order upwind differencing of advection.

iv) Taylor- Galerkin

The Taylor-Galerlun method (Donka, 1984) bears a close resemblance to the Lax- Wendroff schemes of the finite difference world. The advection equation (1) is expanded in a Taylor series to give:

- d'

ae

a t 2 a 3 ~

-

+--.-+---

+ o ( a t 3 )

a t at 2 at- 2 at3 ⁽¹⁴⁾

where At is the time step and n is the time level. The time derivatives are replaced by spatial derivatives through substitution of (14) into (1). Dropping higher-order terms, we obtain the second-order version of the scheme:

en+' ^-C"

+

^{u - v e =} -{v*u At ( u * v c ) } .

At 2

The added diffusion term on the right-hand side is anisotropic in the direction of the flow and, as is the case with the SUPG method in (iii), crosswind diffusion is minimised. Equation (15) is solved through the application of the Galerkin representation. Donka (1984) has als0 derived a third-order version of the Taylor- Galerlun method, which has als0 been tested in this study.

v) Flux-Corrected Transport

As we will show, second-order methods such as the standard Galerkin and Taylor- Galerkin appraoches will generate oscillatory behaviour in the absence of diffusion.

T o maintain positive, monotonic solutions in the presence of strong advection, flux- corrected transport (FCT) algorithms can be employed. Lohner et al. (1987) developed a finite element FCT method which combines a higher-order scheme (a second-order Taylor-Galerkin with a consistent mass matrix) and a low-order solution (a lumped-mass Taylor-Galerkin with added diffusion). In this study, the Lohner et al.

(1987) FCT scheme was tested with both second and third-order Taylor-Galerkin methods for the high-order scheme.

vi) Runge-Kutta Discontinuous Galerkin Method

The discontinuous Galerkin scheme for multi-dimensions (Cockburn et al., 1990) is based on the flux-conservative form of the advection equation:

where f ( c ) is the flux vector, and in this study is defined as f (c) = u c . The discontinuous Galerkin representation is obtained by obtaining a triangulation T, of the domain Q , multiplying (16) by a weighting function v, in the finite element space

V,,

integrating over the element K of the triangulation T, and replacing the exact solution c with its approximation c, E

y,

^:

c,v,dx

+ j

v*f ( r ) v,dx ⁼ O

dt K K

Integrating by parts we obtain:

where e is an edge of the element and n,,, is an outward-directed unit normal vector at e . Replacing f (c)*neXK with the outward-directed discontinuous flux at e :

Replacing the integrals with quadrature rules:

jK

f (c)*vv,dx =

x

M ^{o,f (c)*v v,,}

IKI

]=l

where o, are quadrature coefficients, we obtain the weak formulation

The equation (21) is solved over all the elements, using a lumped-mass time derivative and employing a second order Runge-Kutta scheme to produce the Runge- Kutta Discontinuous Galerkin (RKDG) method (Cockburn and Shu, 1998). In this study triangular elements with linear basis functions are used, resulting in a second- order scheme in time and space. The rnethod outlined above will not be free of oscillations. Slope limiting is applied to deal with this problem. In this study, the minmod and total variance bounded in the means (TVBM) slope limiters described by Cockburn et al. (1990) are employed.

vii) Cell-Based Finite Volurne

MUSCL (Monotone Upstream-Centred Schemes for Conservation Laws) (van Leer, 1979) type schemes for cell-centred, triangular, finite volume approaches (Hubbard, 1999) were investigated in this study. The governing equation is the flux form of the

advection equation as given in (16) and is integrated over a control volume K , applying the divergence theorem:

where n is outward-directed unit vector, normal to the boundary aK of the control volume. Approximation of the boundary integral in (22) leads to the finite volume discretisation

where 'Fis the average value of c in the control volume K, L is the number of edges of the control volume and n, is the outward-directed unit vector normal to the 1" edge, flm is the approximation to the flux at the l" edge. Assuming 'F is piecewise constant within each cell, a first-order upwind scheme is obtained by defining

1 1

-

f'(co, c, )*n, = -(fo ₂

+

f, )-n, --li-n, ₂

I

^{(c, -co)}

where co is the value of c in the cell under consideration, c, is the concentration in an adjacent cell and is an appropriate average of the advection velocity vector

i = af/ac evaluated from the solution values of co and c,.

To obtain a higher-order scheme, higher-order reconstruction can be applied. To obtain second-order spatial accuracy, a linear approximation to the solution across the cell is required. A linear reconstruction of c within a cell can be expressed as

c = F + r * G (25)

where r is a position vector relative to the centroid of the cell and G is a gradient operator, to be defined. G i s selected so as to maximise the magnitude of the concentration gradient in the cell without generating any new extrema at the midpoints of the cell edges. Thus, the numerical flux function of (24) at a cell edge can now be written as

f'(c,,c,) =f'(co +r,,-Go,c, +r,,-G,) (26) where r, is the vector from the centroid of cell i to the midpoint of the edge between cells i and j , G, is the gradient of the reconstructed solution in cell i . c, is the interior reconstructed solution value relative to the cell under consideration and c, is the corresponding exterior value, taken from the adjacent cell (which is generally different).

The gradient limiters examined in this study are based on gradients computed from neighbouring triangles. In Fig 1, the gradients can be computed from

F,

, k = 0,1,2,3.

The triangles defined by connecting the centroids are designated A,, . The Limited Central Difference is based on the gradient from A,,, , and is designated V(A,,3).

The gradient G is determined from a limiting operation where

max (c, - c, , O ) _.

if ro,* G > max (c, - c,, O)

otherwise and the LCD limited gradient is obtained from

G = a v ( A l , 3 ) = ( ^1=1,2,3 mlnu.

')

V ( A 1 2 3 )

Fig 1. Configuration of cell centred, triangular finite volume scheme

The Maximum Limited Gradient (MLG) scheme of Batten et al. (1996) determines a reconstruction based on the four gradient operators

v ( A 0 2 3 ) 3 v ( A 1 0 3 ) , v ( A 1 2 0 )

which are limited in the same manner as in (27) and (28) to produce the limited gradients

G,, = a o V ( A 1 2 3 ) , G, = a , V ( A 0 2 3 ) ,

(29)

'

2 =%v(A101)9 G 3 % V ( A ~ 2 ~ )

The MLG gradient operator is then taken to be the G, of (29) with the largest slope

I G ~ / .

In one dimension, the MLG reduces to the Superbee limiter (Sweby, 1984). The

Durlofsky et al. (1992) limiter consists of selecting the gradient with the maximum slope from

{ v ( A 0 2 ) v ( A 1 0 3 ) v ( A 1 2 0 ) O)

which does not create new extrema at the midpoints of the cell edges.

Time integration is accomplished through a second-order Runge-Kutta scheme.

viii) Control Volt~me Finite Element Eulerian-Lagrangian Method

Eulerian-Lagrangian Methods (ELM) have become increasingly popular in a variety of fields (Sorek, 1988; Binning and Celia, 1996; Rasch and Williamson, 1990;

Staniforth and Cote, 1991), combining the accuracy of the Lagrangian approach with the convenience of a fixed grid. ELMs, however, can have serious mass conservation errors (Baptista, 1987; Russell, 1989). To overcome mass conservation problems, the Eulerian-Lagrangian Localized Adjoint Methods (ELLAMs) have been developed (Celia et al., 1990; Russell, 1989; Binning and Celia, 1996). We wished to include an ELLAM scheme in the comparison study, but the duration of this pre-project was too limited to permit the development and implementation of such complex schemes.

Thus, we were very pleased to co-operate with Dr. Anabela Oliveira of the National Civil Engineering Laboratory in Lisbon, Portugal, who has recently developed a new, advanced scheme, a control volume finite element Eulerian-Lagrangian method, designated as VELA (Oliveira, 1997). VELA combines the mass-conservation properties of the control volume finite element schemes with a new high-accuracy quadrature integration technique for Eulerian-Lagrangian calculations.

References

Asselin, R., Frequency filters for time integrations, Mon. Weather Rev., 100: 487-490, 1972.

Baptista, A.M., Solution of Advectiorz-Dominated Transport by Eulerian-Lagrangian Methods using the Backwards Method of Characteristics, Ph.D. Dissertation, Massachusetts Institute of Technology, Cambridge, 1987.

Batten, P., C. Lambert and D.M. Causon, Positively conservative high-resolution convection schemes for unstructured elements, Int. J. Numer. Methods Fluids, 39: 1821-1838, 1996.

Binning, P. and M.A. Celia, A finite-volume Eulerian-Lagrangian localized adjoint method for solution of the contaminant transport equations in two-dimensional multiphase flow systems, Water Resources Research, 32: 103-1 14, 1996.

Brooks, A.N. and T.J.R. Hughes, Stream-line upwindfPetrov-Galerkin formulation for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations, Conzput. Methods Appl. Mech. Engrg., 32: 199-259,

1982.

Celia, M.A., T.F. Russell, I. Herrera and R.E. Ewing, An Eulerian-Lagrangian localized adjoint method for the advection-diffusion equation, Adv. Wat.

Resourc., 13: 187-206, 1990.

Cockburn, B., S. Hou and C.W. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws IV: The multidimensional case, Math. Comp., 54: 545-581, 1990.

Cockburn, B., and C.W. Shu, The Runge-Kutta discontinuous Galerkin finite element method for conservation lasw V: Multidimensional systems, J. Comput. Phys., 141: 199-224, 1998.

Donea, J., A Taylor-Galerkin method for convective transport problems, Int. J.

Numer. Math. Engrg., 20: 101-1 19, 1984.

Durlofsky, L.J., B. Enquist and S. Osher, Triangle based adaptive stencils for the solution of hyperbolic conservation laws, J. Comput. Phys., 98: 64-73, 1992.

Finlayson, B.A., Numerical Methods for Problems with Moving Fronts, Ravenna Park Publishing, Seattle, 605 p, 1992.

Godlewski, E. and P.A. Raviart, Numerical Approximation of Hyperbolic Systems of Consewation Laws, Springer-Verlag, New York, 509 p, 1996.

Hirsch, C., Numerical Computation of Intemal and External Flows, Vo1.2, Computational Methods for Inviscid and Viscous Flows, John Wiley and Sons, New York, 691 p, 1990.

Hubbard, M.E., Multidimensional slope limiters for MUSCL-type finite volume schemes on unstructured grids, J. Comput. Plzys., 155: 54-74, 1999.

Kroner, D., Numerical Schemes for Consewation Laws, Wiley-Teubner, Chichester, 508 p, 1997.

Lohner, R., K. Morgan, J. Peraire and M. Vahdati, Finite element flux-corrected transport (FEM-FCT) for the Euler and Navier-Stokes equations, Int. J. Num.

Methods Fluids, 7: 1093-1 109, 1987.

LeVeque, R.J., Numerical Methods for Coasewation Laws, Birkhauser Verlag, Basel, 214 p, 1992.

Lynch, D.R. and F.E. Werner, Three-dimensional hydrodynamics on finite elements.

part

II:

Non-linear time-stepping model, Int. J. Numer. Methods Fluids, 12:

507-533,1991.

Morton, K.W., Numerical Solution of Convection-Difision Problems, Chapman and Hall, London, 372 p, 1996.

Oliveira, A.P., Eulerian-Lagrangian analysis of Transport and Residence Times in Estuaries and Coasts, Ph.D. Dissertation, Oregon Graduate Institute of Science and Technology, Portland, 1997.

Rasch, P.J., and D.L. Williamson, On shape-preserving interpolation and semi- Lagrangian transport, SIAM J. Sci. Stat. Cornput., 11: 656-687, 1990.

Russell, T.F., Eulerian-Lagrangian localized adjoint methods for advection-dominated problems, in Numerical Analysis 1989, Pitma~z Research Notes, D.F. Griffiths and G.A. Watson (editors), Longman Scientific and Technical, 228: 206-228, 1989.

Segal, A, Finite element methods for advection-diffusion equations, in Numerical Methods for Advection-Difision Problems, Vreugdenhil, C.B. and B. Koren (editors), Vieweg, Braunschweig, Wiesbaden, pp 195-214, 1993.

Sorek, S., Eulerian-Lagrangian method for solving transport in aquifers, Adv. Wat.

Resourc.,ll: 67-73, 1988.

Staniforth, A., and J. Cote, Semi-Lagrangian integration schemes for atmospheric models - a review, Mon. Wea. Rev., 119: 2206-2223, 199 1.

Sweby, P.K., High resolution schemes using flux limiters for hyperbolic conservation laws, SIAM J. Numer. Anal., 21: 995-1011, 1984.

Toro, E.F., Riemann Solvers and Numerical Methods for Fluid Mechanics, Springer, Heidelberg, 624 p, 1999.

Vreugdenhil, C.B. and B. Koren (eds.), Numerical Methods for Advection-Difision Problems, Vieweg, Braunschweig, Wiesbaden, 370 p, 1993.

3. One-Dimensional Test Cases

The behaviour of the schemes in one-dimension was examined with idealised test cases in which a rectangular box and a triangular distribution of concentration is advected from left to right with a velocity of 1. The grid size,Ax, was 0.02 on a domain from x ⁼ O to x = 2 . Both the box and triangle had an initial width of 0.4, and an amplitude of 1.0, centred at x ⁼ 0.5 at t ⁼ O. The features were then advected for the time interval t = O to t = l . The time step, A t , was in all cases 0.005, for a Courant number C, = 0.25 . On the upstream ( x = O) boundary the Dirichlet condition c = O was imposed, while on the downstream boundary ( x = 2) the Neumann ac/ax = O condition was specified.

Behaviour of the Galerkin finite element, leap-frog lumped-mass Galerkin, streamline upwind Petrov-Galerkin, Taylor-Galerkin, flux-corrected transport, Runge-Kutta discontinuous Galerkin and weighted average flux schemes is shown in Figs 2-12.

The weighted average flux (WAF) scheme with Superbee limiter (Toro, 1999) is included, since it is the one-dimensional equivalent of the MUSCL finite volume scheme with a MLG limiter. In all cases, the exact solution is shown as a solid line, while the numerical solution is a dashed line.

Shown in Fig 2 are the results for the standard Galerlun finite element method. The box case shows strong oscillatory behaviour and negative values, but the steepness of the box sides is well-represented. Slight oscillatory behaviour is evident in the triangle case, but the accuracy overall is extrernely good.

It has been reported in the literature, e.g., Segal (1993), that if one maintains the consistent mass matrix in the Galerkin finite element scheme, well-behaved results can be acheived with mesh Peclet numbers = uAx/D, where D is a diffusion coefficient, as large as 15, whereas a lumped-mass treatment would require = 2 for stability. To test this possibility, diffusion was added the standard Galerlun finite element model to produce a mesh Peclet number of 15. The results are shown in Fig 3. The oscillations in the box case have nearly disappeared, while the amplitude of the solution in the triangular test case has only been reduced by 15 percent.

The approach commonly used in three-dimensional simulations (e.g., Lynch and Werner, 1991), is to use a leap-frog time-stepping, lumped mass matrix approximation to the Galerkin finite element method. The results of this approach on the one-dimensional tests are shown in Fig 4. It can be seen that the box simulation is plagued by severe dispersive ripples, while the scheme als0 appears to be dissipative.

Phase errors are evident in the triangular test, where the numerical solution is leading the exact concentration distribution.

The streamline upwind Petrov-Galerlun results are shown in Fig 5. For these test cases the scheme reduces to the upwind advection scheme. Thus, while the results are positive and monotonic, the scheme is very diffusive.

Shown in Fig 6 are the results for the second-order Taylor-Galerkin scheme.

Oscillatory behaviour occurs near the sides of the box, but is much less severe than in the standard Galerkin approach. The treatment of the triangular test case is very accurate, with negligible oscillatory behaviour.

The third-order Taylor-Galerkin results are shown in Fig 7. The box results show the essentially the same degree of oscillation as the second-order case, but the oscillation pattern is reversed. The triangle test case results are even more accurate than in the second-order case.

The finite element flux-corrected transport (FCT) scheme based on the second-order Taylor-Galerlun method as the high-order scheme are shown in Fig 8. The results are positive and monotonic. The treament of the box case is extremely accurate. The amplitude of the triangle is well-reproduced, but the FCT scheme develops a discontinuous step structure on the sides of the triangle. Simulations were also conducted with a third-order Taylor-Galerkin method as the high-order scheme, but the results are indistinguishable from those with the second-order scheme shown in Fig 8.

The results from the Runge-Kutta discontinuous Galerlun (RKDG) scheme with no slope limiting are shown in Fig 9. These results are nearly identical to the second- order Taylor-Galerlun results of Fig 6. This is not too surprising since both schemes are second order.

The RKDG results using a TVBM slope limiter are shown in Fig 10. Oscillations are nearly eliminated, with only slightly negative values apparent on the leading edges of the concentration features.

The minmod limiter eliminates all oscillations in the RKDG results and all values are positive in Fig l l. An asymmetry develops in the box results between the leading and trailing sides of the box. There is also greater attenuation of the amplitude in the triangle case.

The results from the weighted average flux (WAF) method with a Superbee limiter are shown in Fig 12. This method was included to emulate the one-dimensional behaviour of the MUSCL finite volume scheme with MLG (maximum limited gradient). The box case is well-represented, with symrnetrical behaviour of the leading and trailing edges. The triangle peak is somewhat attenuated and spread as in the RKDG minmod results.

Galerkin

Galerkin

Fig 2. Galerkin finite element on box and triangle test cases.

Galerkin Pe=l 5

Galerkin Pe=15

Fig 3. Galerlun finite element with

P,

= 15 on box and triangle test cases.

Leap-Frog, Lumped Mass Matrix Galerkin

Leap-Frog, Lumped Mass Matrix Galerkin

Fig 4. Leap-frog lumped mass Galerkin on box and triangle test cases

Fig

Petrov-Galerkin

Petrov-Galerkin

Taylor-Galerkin 2nd Order

Taylor-Galerkin 2nd Order

Fig 6. Second-order Taylor-Galerkin on box and triangle test cases.

Taylor-Galerkin 3rd Order

Taylor-Galerkin 3rd Order

Fig 7. Third-order Taylor-Galerkin on box and triangle test cases.

FCT

FCT

Fig 8. Flux-Corrected Transport with 2nd-order Taylor-Galerkin

RKDG with no limiter

RKDG with no limiter

Fig 9. Runge-Kutta Discontinuous Galerkin with no limiter on box and triangle test cases.

RKDG with TVBM limiter

RKDG with TVBM limiter

Fig 10. Runge-Kutta Discontinuous Galerkin with TVBM limiter on box and triangle test cases.

RKDG with MINMOD limiter

RKDG with MINMOD limiter

Fig l l. Runge-Kutta Discontinuous Galerkin with minmod limiter on box and triangle test cases.

WAF with SUPERBEE limiter

WAF with SUPERBEE limiter

Fig 12. Weighted Average Flux with Superbee limiter (MUSCL finite volume with MLG limiter) on box and triangle test cases.

4. Two-Dimensional Test Cases

Two-dimensional test cases were constructed for the schemes in this study using the triangular mesh of Fig 13. The mesh contains 4225 nodes and 8192 triangles.

Fig 13. Triangular mesh for 2-D test cases.

The square domain extends from -1 to l in both x and y . The velocity field is solid body rotation, such that V = (-2ny, 2 n ~ ) ~ . For the cylinder test case, the initial conditions are (see Fig 14):

l for r 5 0.25 O for r > 0.25

and for the cone test case, the initial conditions are (see Fig 15):

cos2 ( 2 n r ) for r 5 0.25 O for r > 0.25

where r' = ( x - 0 . 5 ) ~

+

^y'. The initial fields are rotated for one full revolution (from t = O until t = l ) , where the initial distribution of c should be unchanged from the initial conditions. The solution is set to zero on the inflow boundaries. For all schemes except VELA, a time step of 2 . 9 1 8 ~ l o 4 is used, which corresponds to C,. 5 113, a stability requirement for schemes employing Runge-Kutta time integration. For the VELA method, which has no time-step stability restriction, a time step of 0.02 was used.

The results from the tests on the various schemes are summarised in Tables 1 and 2 and are presented in Figs 16-40. The contour interval for concentration on all the contour plots is 0.05.

1

^Scheme

_I

_,,c

_l

_cm,

_I

_L-

Standard Galerkin finite element Galerlun FEM with

P,

⁼ 15 Leap-frop. lumped mass Galerkin LF lumped Galerkin with

P,

= 2 Streamline u ~ w i n d Petrov-Galerlun

-0.306 - 3 . 3 4 ~ -0.493

2"" Order Taylor-Galerkin 3rd Order Tavlor-Galerkin

O.

-0.19 1

FCT with 2nd order TG FCT with 3rd order TG MUSCL with no limiter

1.333 0.900 1.579

-0.218 -1.081

MUSCL with LCD MUSCL with MLG

0.747 0.632 0.902 0.286

1.182

-3.68 X lo-' -4.13 X 10.' -0.0754

MUSCL with Durlofsky et al.

RKDG with no limiter

I

^{W L A}

1

^-0.0705

1

^1.092

1

^0.538

1

0.846 0.574 1.252

1.880

O.

- 2 . 9 4 ~ 10'"

Table l. Summary of results from cylinder test case.

0.607 1.081 0.994

0.999 1.109

-0.0980 -0.0597

0.553 0.570 RKDG with TVBM

RKDG with MINMOD

/

^Scheme

_l

_cmin

₁

_cm,

_{l r,}

0.570 0.576 0.590 0.951

1 .O00

0.644 0.65 1 0.801

1.084 -0.0401

-4.97 X lo-''

0.678 0.549 1.058

0.999

Standard Galerkin finite element Galerlun FEM with

P,

⁼ 15 Leap-frog lumped mass Galerlun LF lumped Galerkin with

P,

⁼ 2

-0.0250 - 4 . 1 5 ~ 10‘'

Streamline upwind Petrov-Galerlun 2nd Order Tavlor-Galerkin

-0.232 O.

3ra Order Taylor-Galerkin FCT with 2"d order TG

MUSCL with Durlofsky et al.

I

^O.

1

^0.338

1

^0.663

0.984 0.424

-0.01 18 -0.0213

0.030 0.582 0.932

0.091

-0.0632 - 3 . 2 5 ~ 10-l2 FCT with 3rd order TG

MUSCL with no limiter MUSCL with LCD MUSCL with MLG

0.618 0.911 0.922

0.984

0.559 0.959 0.539 0.835 - 2 . 4 6 ~ lo-''

-0.0207 O.

- 7 . 7 5 ~ 10-l' RKDG with no limiter

RKDG with TVBM

Table 2. Summary of results from cone test case.

0.076 0.024 1.052

0.539

0.435 0.056 0.475 0.191

RKDG with MLNMOD VELA

0.091 0.455

-0.0102 -0.0099 - 8 . 2 6 ~ -0.0070

0 . 9 7 5 0.870

0.0%

0.130 0.752

0.987

0.251 0.015

It should be noted that the RKDG and MUSCL schemes in this study are element or cell-based. Since there are roughly twice as many triangles as nodes, the effective grid size, or cell length, for these schemes is approximately l/& the cell length of the rest of the methods.

It can be seen (e.g., Figs 27,32,40) that if the concentration distribution is smooth and well-resolved spatially, then the second-order schemes Taylor-Galerlun, RKDG and MUSCL, without limiters, VELA and SUPG provide excellent results. If sharp discontinuities are encountered, or if it is essential that positivity be maintained, then the RKDG scheme with a minmod limiter and the MUSCL scheme with a MLG limiter produce the best results.

If the scheme must be vextex (or node) based and be suitable for inclusion in a three- dimensional model, then the second-order Taylor Galerlun scheme and the flux- corrected transport method based upon the second-order Taylor-Galerkin can be considered. The FCT scheme is robust (Fig 29), although its clipping action will tend to reduce peak concentrations (Fig 30). It is because of its robustness and compatibility with an existing three-dimensional finite element model that the FCT scheme with second-order Taylor-Galerkin is implemented for the three-dimensional simulations.

Fig 14. Exact solution for cylinder test case.

Fig 15. Exact solution for cone test case.

Galerkin, Consistent Mass Matrix

Fig 16. Results from standard Galerkin finite element on cylinder test case.

Galerkin, Consistent Mass Matrix

Galerkin, Consistent Mass Matrix

Fig 17. Results from standard Galerkin finite element on cone test case.

Galerkin, Consistent Mass Matrix 1

b

^P,=15

Fig 18. Results from Galerkin finite element with

c

⁼ 15 on cylinder test case.

Galerkin, Consistent Mass Matrix

k

Fig 19. Results from Galerlun finite element with = 15 on cone test case.

Galerkin, Lumped Mass Matrix Leap-Frog

Fig 20. Results from leap-frog lumped mass Galerkin on cylinder test case.

Galerkin, Lumped Mass Matrix Leap-Frog

X

Fig 21. Results from leap-frog lumped mass Galerkin on cone test case.

Galerkin, Lumped Mass Matrix Leap-Frog, P,=2

Fig 22. Results from leap-frog lumped mass Galerkin on cylinder test case, with P, = 2 .

ss Matrix

Galerkin, Lumpeel Mass Matrix Leap-Frog. P,=2

Fig 23. Results from leap-frog lumped mass Galerkin on cone test case, with

P,

⁼ 2 .

SUPG

SUPG

Fig 24. Results from streamline upwind Petrov-Galerlun on cylinder test case.

SUPG

SUPG

Fig 25. Results from streamline upwind Petrov-Galerlun on cone test case.

Taylor-Galerkin 2nd Order

0.5

. . - -

Fig 26. Results from second-order Taylor-Galerkin on cylinder test case.

Fig 27. Results from second-order Taylor-Galerkin on cone test case.

1

0.5

*

^{0 -}

-0.5

-1 -1

-

Taylor-Galerkin 2nd Order

- -

- -

-

-

- -

- -

-

- -

-

-

- - -

' " " I I i ' i l l l l l l l i I I

-0.5 O 0.5 1

X

Taylor-Galerkin 3rd Order

Fig 28. Results from third-order Taylor-Galerkin on cylinder test case.

Tayior-Galerkin 2nd Order with FCT

0.5

Fig 29. Results from flux-corrected transport with 2nd order Taylor-Galerkin on cylinder test case.

Taylor-Galerkin 2nd Order with FCT

Taylor-Galerkin 2nd Order with FCT

Fig 30. Results from flux-corrected transport with 2nd order Taylor-Galerlun on cone test case.

Local Discontinuous Galerkin

' I

No Limiter

Fig 3 1. Results from Runge-Kutta discontinuous Galerkin (RKDG) with no slope lirniting on cylinder test case.

Local Discontinuous Galerkin

Fig 32. Results from Runge-Kutta discontinuous Galerkin (RKDG) with no slope limiting on cone test case.

Local Discontinuous Galerkin TVBM

Fig 33. Results from Runge-Kutta discontinuous Galerlun (RKDG) with TVBM slope limiter on cylinder test case.

1

0.5

0 -

-0.5

-1 -

Local Discontinuouc Galerkin

-

^TVBM

- - - -

-

- -

- -

- - - -

-

-

-

- - - -

l l l l l l l l l l l ' l l i l lI ll

1 -0.5 O 0.5 1

X

Local Discontinuous Galerkin l

Fig 34. Results from Runge-Kutta discontinuous Galerkin (RKDG) with TVBM slope limiter on cone test case.

Locai Discontinuous Galerkin MINMOD

0.5

Fig 35. Results from Runge-Kutta discontinuous Galerkin (RKDG) with minmod slope limiter on cylinder test case.

Local Discontinuous Galerkin

k

^M'NMoD

Fig 36. Results from Runge-Kutta discontinuous Galerkin (RKDG) with minmod slope limiter on cone test case.

MUSCL M LG

Fig 37. MUSCL-type finite volume with MLG slope limiter on cylinder test case.

1

0.5

0 -

-0.5

-1 -1

MUSCL

- MLG

- -

-

- -

- -

-

-

- - -

-

-

- - - -

-

" " ' " " " " " ' " ~ ~

-0.5 O 0.5 1

X

Fig 38. MUSCL-type finite volume with MLG slope limiter on cone test case.

1

0.5

*

^{0 -}

-0.5

-1 -

-

^MUSCL

- MLG

-

-

-

-

- - - - - - - -

-

-

- - -

' " ' l l l " l l l l ' l ' l i i l

1 -0.5 O 0.5 1

X

Fig 39. Results from VELA on cylinder test case.

Fig 40. Results from VELA on cone test case.

5. Three-Dimensional Test Cases

As discussed in the previous section, the flux-corrected transport method based upon the second-order Taylor-Galerkin scheme was implemented for use in a three- dimensional finite element model. The FCT scheme was implemented in Quoddy (Lynch and Werner, 1991), to compute horizontal advective fluxes of momentum, temperature and salinity. The vertical diffusive and advective fluxes were computed as before, with the vertical diffusion treated irnplicitly using a Crank-Nicholson scheme. Because of the positive, monotonic nature of the FCT scheme, it was possible to run the model with zero horizontal viscosity and diffusivity.

The consistent mass matrix was used in the FCT Taylor-Galerkin scheme. This necessitated a matrix inversion. However, since the horizontal geometry and mesh characteristics were the same at each vertical (sigma-coordinate) level, it was only necessary to perform one LU decomposition at the beginning of the simulation, then apply the back-substitution step with the constant LU matrices at each of the 21 sigma-levels and at each time step. Thus, a directional splitting was accomplished between the horizontal, with an unstructured triangular mesh, and the vertical with a structured, vertically-stretched, sigma coordinate system.

Fig 41. Mesh and model domain for the idealised estuary/river plume case. Distances are in km.

Two three-dimensional test cases were used in the study. The first is an idealised riverlestuary supplying a (nearly) freshwater flux of 80.000 m3s-l, with a salinity of 1 psu and temperature of 8OC, entering a shelf sea with salinity 34.8 psu and temperature of -lS°C. The depth of the domain is constant at 10 m. These conditions are roughly representative of river discharge on to the Siberian shelves in the Arctic.

The mesh and the model domain are shown in Fig 41. The typical mesh size is 5 km, with 4450 nodes and 8490 triangles.

The results after 120 days are shown in Figs 42 and 43.

Fig 42. Computed surface salinity at day 100 from the idealised estuarylriver plume test

The low-salinity water enters the shelf region primarily as a rightward-directed coastal jet (Fig 42). An arrested, leftward plume is also evident, consistent with the

theory of river plume dynarnics on shallow shelf seas. In Fig 43 is shown a close-up of the river mouth area with the surface salinity and velocity. Even though no horizontal viscosity of diffusivity are used, the results are smooth and well-behaved while preserving a realistic frontal structure and strong velocity shear across the front.

Previously, it was not possible to obtain a numerically stable to this problem using Quoddy, with its lumped-mass, leap-frog time-stepping advection.

0.5 mi'

Fig 43. Results of the idealised estuary/river plume simulation in the vicinity of the river mouth.

A second test case represents a more realistic application, with actual topography from the Kara Sea region and realistic buoyancy forcing. Shown in Figs 44 and 45 are the finite element mesh and bottom topography for the Kara Sea portion of the model domain. A total of 6663 nodes and 12415 elements were used in the simulation, which included the Barents Sea, as well. River inflow with the same salinity and temperature as in the previous case ( l psu, 8OC) were supplied at the Ob River (60.000 m3s-l) and Yenisei River (80.000 m3s-l), with the river water entering an (initially) homogeneous shelf sea with a salinity of 34.8 psu and a temperature of -1.5 OC

.

Fig 44. Finite element mesh for the Kara Sea portion of the model domain.

Both the second-order Taylor-Galerkin and FCT based on second-order Taylor Galerkin schemes were used to advect heat, salt and momentum in this application.

Using the Taylor-Galerkin scheme alone resulted in numerically-unstable results, with negative salinities being encountered after 20 days. Adding horizontal diffusion improved the results but, with

4

= 15, negative salinities were still encountered after 30 days. Using the FCT approach eliminated the need for horizontal diffusion and viscosity and no numerical difficulties were encountered during the course of the 180 day simulation. Shown in Fig 46 are the results from day 120 of the simulation. Fig 47 shows the detailed surface salinity and velocity fields in the vicinity of the river mouths. As was the case for the idealised estuarylriver plume sirnulation, the FCT algorithm produces well-behaved tracer and velocity fields with a weil-defined frontal structure and strong velocity shear in the frontal zone.

Fig 45. Bathyrnetry for the Kara Sea portion of the model domain.

fig 46. Surface salinity from day 120 of the sirnulation.

Fig 47. Surface salinity and velocity near the Ob and Yenisei River mouth at day 120 of the simulation.

6. Conclusions

Several schemes examined in this study perfonned well and would prove useful for coastal ocean transport problem. Based on the results of the idealised one and two- dimensional tests, we recornmend that the MUSCL finite volume, RKDG and VELA methods be combined with slope limiting schemes and be used to examine systems of equations. These schemes are particularly suitable for multi-dimensional systems of equations, since they do not require the maintenance of a consistent mass matrix for high accuracy results. The schemes should first be applied to the two-dimensional shallow water equations on a rotating reference frame, then to the full three- dimensional primitive equations.

Acknowledgements

We wish to thank Dr. Anabela Oliveira for conducting the two-dimensional test case simulations with VELA and making the results available to us. A more comprehensive description of VELA will be included in a joint paper with Dr.

Oliveira. We als0 wish to thank Professor Magne Espedal and Dr. Helge Dahle and Kenneth Hvistendahl Karlsen at the Department of Mathematics, University of Bergen for pointing us to recent literature and their recornrnendations concerning ELLAM, discontinuous Galerkin and finite volume schemes.