A primitive pseudo wave equation formulation for solving the harmonic shallow water equations

A primitive pseudo wave equation formulation for solving the harmonic shallow water equations J. J. Westerink

Ocean Enyineerin9 Program, Civil Engineerin9 Department, Texas A&M University, Colleye Station, TX 77843, USA J. J. Connor and K. D. Stolzenbach

Department of Civil Engineerin9, Massachusetts Institute of Technology, Cambridye, MA 02139, USA

A finite element method formulation for solving the harmonic shallow water equations in their primitive or unmodified form is developed and analysed. The scheme, referred to as the Primitive Pseudo Wave Equation Formulation (PPWE), involves developing a weak weighted residual form of the continuity equation and furthermore forming a pseudo wave equation by substituting the discretized form of the momentum equation into the discretized form of the continuity equation. The final set of equations to be solved, the pseudo wave equation and the primitive momentum equations, deceptively resemble the discretized equations of the wave equation formulation of Lynch and Gray. Despite this resemblance, Fourier analysis indicates that the PPWE scheme is still fundamentally primitive. However, application of the PPWE scheme to a set of stringent test problems results in very good solutions with well controlled nodal oscillations. It is shown that this low degree of spurious oscillations is due to the treatment of the boundary conditions such that elevation is taken as an essential condition and normal flux is taken as a natural condition. This particular boundary condition treatment is suggested by the formation of the pseudo wave equation. Furthermore, even though the equation re-arrangement does not in itself effect the solutions, it does make the scheme much more efficient. Key Words: shallow water equations, harmonic, primitive, spurious oscillations, wave equation

1. INTRODUCTION The numerical solution of the shallow water equations has been a challenging problem of long standing interest. The use of the finite element (f.e.) method to resolve the spatial dependence in these equations has by now become a widespread solution technique. Finite difference schemes have been largely retained to resolve the time dependence although harmonic based schemes which entirely eliminate time as a variable have also come into use. However most of the proposed formulations have been severely troubled by short wavelength spatial noise which manifests itself as unphysicai node to node oscillations. The control of these spurious oscillations is vital to the integrity of both time domain and frequency domain solution procedures. In attempting to overcome this difficulty, investigators have resorted to adding artificial damping 1 or selecting numerical schemes which in of themselves introduce numerical damping 2. The problem with this of course is that not only is the short wavelength noise damped but the longer physical waves are damped as well. Accepted May 1987. Discussion closes February 1988. © 1987 Computational Mechanics Publications

188 Adv. Water Resources, 1987, Volume 10, December

More recently Gray and Lynch 3'¢ have developed a number of schemes which offer excellent control over the short wavelength spatial noise without affecting the fundamental solution. Their solution strategies share a common basis of modifying the shallow water equations prior to the application of temporal and/or spatial numerical discretizations. Their semi-implicit method [SIM) rearranges the time discretized form of the shallow water equations before spatial discretization by solving for velocities in the momentum equations and then substituting these into the flux terms in the continuity equation 3. The resulting equation, which includes a second derivative of elevation with respect to space, is now solved together with the standard momentum equations through the application of spatial discretization procedures. The SIM procedure results in a significant reduction in the amplitudes of the node to node oscillations experienced when compared to schemes based on the unmodified or so called primitive form of the shallow water equations. An even more successful class of schemes is based on the wave equation (WE) form of the shallow water equations which consists of a wave type equation and the standard primitive momentum

A primitive pseudo wave equation formulation: J. J. Westerink et al. equation 4. This wave equation is found by differentiating the primitive continuity equation with respect to time and then substituting for the resulting time differentiated flux terms by using the solution for the acceleration terms in the momentum equations. Finally, a 15ortion of the nonlinear friction term (or the entire linearized friction term) is then substituted for using the solution for velocity obtained with the continuity equation. This wave equation, which has elevation as the basic variable and includes second derivative terms of elevation with respect to both space and time, together with the primitive momentum equations are then used as the starting point of any temporal and spatial numerical discretization procedures. Lynch and Gray demonstrate that the WE formulations work very well for both time domain and harmonic approaches 4'5. It can be readily shown that a harmonic formulation which follows the equation substitution steps of the SIM approach leads to essentially the identical undiscretized equations as the harmonic WE formulation for the linearized case with constant depth and constant linear friction factor. In this paper we develop and analyse a formulation for solving the harmonic shallow water equations which in a vein similar to the WE formulation is based on the formation of a wave type equation. However, unlike the WE scheme, our formulation discretizes the shallow water equations in their primitive form and only then goes on to form a wave type equation by substituting the discretized momentum equation into the discretized continuity equation. This re-arrangement results in a final discretized set of equations consisting of a pseudo wave equation and momentum equations which closely mimic the final discretized equations obtained with Lynch's 5'° harmonic WE formulation. In fact, this pseudo wave equation includes terms which correspond to second derivatives of elevation in both time and space. However analysis shows that despite this deceptive similarity in form to the WE formulation, our method still exhibits the fundamental properties of a primitive formulation and we shall therefore refer to it as a primitive pseudo wave equation (PPWE) formulation. Nonetheless the actual performance achieved by our P P W E formulation is far superior to that of any other formulation based directly on the primitive form of the shallow water equations. Although the oscillatory behaviour of the PPWE formulation is not quite as well controlled as that of the harmonic WE formulation, the PPWE method has the advantage of, in general, involving fewer and simpler terms. This increased simplicity of the equations in conjunction with the decoupling of the solutions for elevation and velocity achieved through the formation of the pseudo wave equation can be especially beneficial to the computational efficiency of iterative harmonic formulations which solve the full nonlinear shallow water equations. Qualitatively the degree of node to node oscillation of the P P W E method is on the same order as that of the SIM formulation. It will be shown that this well controlled behaviour is due to the alternative handling of the boundary conditions in which elevation is treated as an essential boundary condition and fluxes are taken as natural boundary conditions. This particular boundary condition treatment is crucial to the success of the scheme and is suggested by the formation of the pseudo wave equation which is solved first and has elevation as its basic variable.

2. P R I M I T I V E P S E U D O WAVE E Q U A T I O N FORMULATION

a. Governing equations The shallow water equations consist of a set of depth averaged momentum equations and a depth averaged continuity equation. They describe tidal wave propagation and many other depth averaged type free surface circulation problems. In their linearized form these equations are expressed as: q,, + (uh),x + (vh) ~,= 0

(1)

u., + gtl, x -- fi, + ~ u = 0

(2a)

2 v., +gq,r + f u + ~ v = O

(2b)

tl

where t x,y u,v

h g

f

2

time cartesian coordinates depth averaged components of velocity in the x,y coordinate directions surface elevation relative to mean sea level (MSL) depth to MSL acceleration due to gravity Coriolis factor linearized bottom friction coefficient

The associated boundary conditions are the prescription of elevation, q*, and the prescription of normal flux, Q*, which are respectively expressed as:

~l(x,y,t)=~l*(x,y, t) on F,

(3a)

Q,(x,y, t)=Q*(x, y, t) on F~

(3b)

For the linearized equations normal flux is computed as:

Q, = Ot,xuh + ~,fl,h

(4)

where ~,x and ~,~, are the direction cosines on the boundary. For periodic boundary and/or other periodic forcings, these linearized equations greatly simplify by assuming that the responses in both elevation and velocities are also periodic and of the form:

q(x, y, t) = Re{~(x, y) e i'''',

(5a)

u(x, y, t)= Re{fi(x, y)e i''''',

(5b)

v(x, y, t)= Re{ b(x, y) e/'''']

(5c)

where ~) fi,b i ~o

complex amplitude of elevation complex amplitudes of velocities ( - 1) 1/2 frequency

Substitution into (1) through (3) leads to the harmonic form of the linearized shallow water equations which are no longer time dependent and are expressed as: io9~ + (fih),~,+ (bh), r = 0

(6) (7a)

Adv. Water Resources, 1987, Volume 10, December

189

A primitive pseudo wave equation formulation: J. J. Westerink et al. i(~+gt}.r+ffi+2?~=O~

(7b)

~/(x,y)=~/*(x,y) on F,

(8a)

Q.(x,y)=Q*(x,y) on F o

(8b)

This set of harmonic equations can be used to compute responses for dominant astronomical tidal components in regions where the nonlinearities are not important or they can be used as the core of fully nonlinear harmonic schemes which are based on either perturbation analyses v or iterative type schemes 8.

b. Weighted residual formulation The primitive form of the harmonic shallow water equations given by equations (6) through (8) represent the starting point of our numerical formulation. In establishing our weighted residual formulation, we must keep in mind that ultimately we will substitute the discretized momentum equations into the discretized continuity equation, resulting in a final system of algebraic equations with elevation representing the unknown variable. This then suggests that elevation be taken as the essential boundary condition. This is turn implies that the continuity equation be used in order to establish the symmetrical weak weighted residual form such that elevation indeed represents the essential boundary condition and furthermore normal flux represents the natural boundary condition. Using Galerkin's method, the error in the continuity equation is weighted by the variation in elevation, fit},and is integrated over the interior domain, ~. The error on the natural boundary is then taken into consideration by weighting the error in the normal boundary flux with fit} and integrating over the flux prescribed boundary, Vq. These combined integrated and weighted interior and natural boundary errors are then required to vanish resulting in the following expression:

{ io.)~l+ (~h).. + (~h),,,1~t} dn c,

+

[{-0,,+0*}f;TdF=0]

(9)

ro

Integrating the terms involving spatial derivatives in equation (9) by parts, considering relationship (4) and accounting for the variation in elevation vanishing on the essential boundaries yields the symmetrical weighted residual weak form for the continuity equation:

(2*@

i

The momentum equations (7) are now multiplied through by depth h such that the application of the f.e. method will yield the same derivative matrices in both the continuity and momentum equations. The weighted residual form for the momentum equations is then derived by weighting each of these equations with a residual velocity and integrating over the interior domain

Adv. Water Resources, 1987, Volume 10, December

{icohh +ght} x - /hb + ).fil fftdf~=O

(lla)

{io~h{,+ght}.,.+lhz]+).b~b~dn=O

(llb)

c. Finite element method formulation The final weighted residual form of the shallow water equations requires that the variables be represented by functions with at least C o continuity. Using identical C o bases for both elevation and velocities in addition to the depth and the linearized friction coefficient, ,;~,allows each of these variables to be represented within the elements by an expression of the form:

a(x, y)=dp(x, y)a !"~

(12)

where

a(x, y) representative variable which stands for h, {', t}, h, ), and the variations 3h, fit and 3t} 4,(x, y) elemental vector of C o interpolating polynomials a ~"~ elemental vector of nodal values for the representative variable Substitution of the elemental expansions of the form given by equation (12) for the variables and parameters into the final weighted residual form of the governing equations (10) and (11), summing over the elements and taking into consideration inter-element functional continuity and the arbitrary variation of 6fi, 6~ and fit} readily leads to the global f.e. form of the continuity and momentum equations: itnn~ - D/.7= -/~.

(13)

itoMvlfr+MFl[r+McC+gDrtl = 0

(14)

where elevation amplitude vector velocity amplitude vector incorporating both x and y components M, continuity equation coefficient matrix Mu momentum equation mass matrix My linearized friction distribution matrix Mc Coriolis matrix P, loading vector for flux prescribed boundaries U

These global matrices and vectors are assembled using the following elemental matrices and vectors:

(15)

dr (10)

190

as follows:

M~v ")=

I

dpr~h~"¥p d~ '

I~;~=

0

~-

0

jr

tI

'fC I"~dpdfl tt ± ',

'ff

~gTdph(n)¢~d

1

(16)

0 (17) ~rq~2¢"'~d

A primitive pseudo wave equation formulation: J. J. Westerink et al.

,j I-

0

M(c,)=

4.-

f

~-~h{"}~d ~ t

~T~bh~"~bdfl '

(18)

e

I

~

e

J

(19)

(20) We have now established the discretized continuity and momentum equations which will be used in order to form our pseudo wave equation. We have implemented our PPWE method using linear interpolating polynomials over triangular elements. This choice of bases leads to elemental systems of three equations for continuity and six equations for momentum. The corresponding gobal system of N equations for continuity and 2N equations for momentum, where N equals the number of nodes, are banded and sparse. Furthermore the matrices M,, Me, and Mr are all symmetric, while the Coriolis matrix, Mc, is skew symmetric.

d. Development of the pseudo wave equation and final solution procedure We now go on to form the pseudo wave equation as was assumed in the development of the weighted residual formulation for our PPWE scheme. This procedure is straightforward and starts by solving the discretized harmonic momentum equation (14) for velocity: t?=

-

gM~JTDTq

(21)

structure for MTor consisting of 2 x 2 submatrices lying along the main diagonal with all other terms equal to zero. The inversion of 191ToT then only involves independently inverting the 2 x 2 submatrices and yields an inverted matrix M~o~rwith exactly the same structure as MroT with all nonzero entries contained in three diagonal bands. With the use of these lumping procedures, the system matrix of (23) is banded with a bandwidth which is double of that of the continuity equation coefficient matrix, M,, due to the product DMTVo~TDT. There are therefore no apparent advantages from a computational efficiency viewpoint to lump M,. Finally this system matrix is always complex and is nonsymmetric when Coriolis effects are taken into consideration. Our PPWE scheme therefore has successfully decoupled the solutions for elevation and velocity with tj being solved for with the pseudo wave equation which has a small and bangled system matrix and with [r being solved for with the primitive momentum equation in an extremely efficient manner when 1VIToris lumped. We note at this point that harmonic formulations which form a pseudo wave equation with velocity as the basic variable are not plausible due to problems relating to the ill-conditioning of the final system matrix produced with this type of scheme. A velocity pseudo wave equation is obtained by substituting for elevation into the discretized momentum equation. Solving for elevation with the discretized continuity equation (13) yields:

0 =.1 M,~-~(DC-/~,,)

(24)

I(t)

Substituting equation (24) into equation (14) and rearranging leads to: ( -- to2 Mr: + iOMF + itoMc

+ gDTM,~- t D)(r = gDTM,; tl~,

where

(25)

~/[TOT-~- (koMt, + Mr + Mc)

(22)

Substituting for /~" into the discretized continuity equation (13) leads to: (/toM, + 9Dr~/I~olTDT)~= -- e ,

(23)

Equation (23) represents our pseudo wave equation and elevation prescribed boundary conditions are directly substituted into this system of equations whereas the natural flux prescribed boundary conditions are accounted for in the load vector/~,. The pseudo wave equation is now solved for elevation independently of the momentum equation. After solving for tj with (23), we substitute for tj into equation (21) and directly solve for/J. The system matrix of the pseudo wave equation (23) is of sufficient rank N to solve for the N unknown elevations. Furthermore for a consistent formulation this system matrix is fully populated due to the inversion of MTOT' Therefore in order for the PPWE scheme to be practical it is necessary that IMror be made into a naturally partitional matrix through lumping procedures applied to the mass matrix, Me,, the frictional distribution matrix, MF, and the Coriolis matrix, Mc. The lumping for the symmetrical mass and frictional distribution matrices is achieved by combining rows onto the diagonal and for the skew-symmetric Coriolis matrix by combining rows onto adjacent off-diagonal positions. This leads to a

Hence equation (25) is solved for velocity/J which may then be backsubstituted into equation (24) in order to obtain elevation. This strategy, however, fails due to a rank problem with the total left hand side system matrix generated in (25). Due to the toz factor multiplying the momentum mass matrix, Mu, and the to factor multiplying the linearized friction and Coriolis matrices, Mr and Me, the contribution of the 9DTM~- 1D matrix dominates the sum of these matrices. The contribution of Mr, MF and M c matrices becomes especially small for larger wavelengths (i.e., lower frequencies) and may disappear entirely due to the roundoff accuracy of the computer. Hence even though the rank of M U, MF and Mc is 2N, the rank of the dominating component of the left hand side system matrix, gDTM,;- 1D, is only N which makes the system matrix in general ill-conditioned and in the case of very long waves equal to rank N. Clearly we can not solve for 2N velocities with a system of 2N equations of only rank N. The implementation of harmonic schemes which substitute the continuity equation into the momentum equation unmistakably demonstrate the severe numerical problems associated with this approach. This therefore indicates that it is only appropriate to form harmonic pseudo wave equations with ~j as the basic variable as with our PPWE formulation.

Adv. Water Resources, 1987, Volume 10, December

191

A primitive pseudo wave equation formulation: J. J. Westerink et 3. ANALYSIS OF T H E P P W E F O R M U L A T I O N

al.

(26) and (27) yields:

Let us now compare the PPWE formulation to primitive and WE formulations and furthermore investigate this new scheme's fundamental behaviour. The final equations that result in the PPWE scheme are the peseudo wave equation and the primitive momentum equation. Implementation of the PPWE scheme to a one dimensional constant depth problem with constant nodal spacing Ax and using a semi-lumped scheme with 191voT but not M, being lumped leads to the following set of discretized equations at each interior node j:

I

2gh ~(1 1 -cos(2aAx))} (ko) 2 + io9 ~). + (2Ax)2

0

1

ti,o+~}J

i{igSin(aAx)~

E,o~o_j __EOol

(29)

where for the semi-lumped PPWE scheme: 1 -+-4rlj 1 +) qJ+ 6

A = 1 (2 cos(aAx) + 4)

(30a)

and for the fully lumped PPWE scheme: 0

(i°°+~)OJ+g(~lJ+ x-~lJ-~-)

(26)

(27)

A= 1

For nontrivial solutions the determinant of the matrix in equation (29) must vanish and thus: {(io9)2 +i{o~+

For the fully lumped PPWE scheme which lumps both lfelror and M, the leading (~j_ 1 + 4~j + ~i+ 1) term in the pseudo wave equation (26) reduces to 6r/j. These discretized equations are very nearly identical to those produced by the harmonic WE scheme with linear elements. The major difference between the equations produced by the PPWE and the WE schemes is that the second spatial derivative in the pseudo wave equation (26) appears as a less accurate centered approximation which uses points two nodes away from n o d e j instead of points directly adjacent to node j as in the WE formulation. The second difference is that the acceleration and friction terms in the discretized momentum equation (27) are in their lumped form whereas these terms are consistently handled in the standard WE formulation. When the fully lumped PPWE scheme is compared to a fully lumped WE formulation, this second difference disappears. We now apply a Fourier analysis scheme to study the PPWE formulation in more detail, to properly classify it and to look into the question of whether a fully or partially lumped scheme will yield better accuracy. We shall follow the propagation factor analysis developed for f.e. schemes by Pinder and Gray 9 and applied by Lynch 6 to study a variety of primitive and WE schemes which use both time domain discretizations and harmonic formulations. The general solution to the one-dimensional linear harmonic shallow water equations consists of series solutions with components of the following form: ~(x) = ~/o e~,X

(28a)

U(x) = U ° e i'x

(28b)

where a r/°,U °

spatial frequency spatially and temporally independent series coefficients

Due to the linearity of the equations only one set of components need be considered. Substituting (28) into

192

Adv. Water Resources, 1987, Volume 10, December

(30b)

2oh 1 (1-cos(2aAx))

i~o+

(2Ax) 2 A

=0 (31)

Like the WE formulation three roots for o) result from the requirement that the determinant vanish. In fact one of these roots is identical for both the harmonic WE and PPWE formulations and corresponds to o) = i2/h which is associated with a simple decay with no propagation. The other two roots are the physical roots which correspond to a progressive and a regressive wave and after some rearranging may be expressed as:

~o=a,Ii2

/gh(sinaAx) 2 (2 ~2 }2-~-a-+X/-A-\~) -\2/~-mm)}

(32)

This expression must be compared with that obtained by performing the same substitution and subsequent analysis for the undiscretized shallow water equations (6) and (7) as was just applied to the discretized PPWE equations (21) through (23). This results in the exact analytical expression relating o) to a:

i2

°t +10"

)

2

(33)

Since as the spatial grid is refined aAx--~ 0 and since in this limit equation (32) reduces to equation (33), it is clear that the PPWE formulation yields consistent solutions. Assuming that the part of equation (32) under the square root sign remains positive, comparison of equation (32) to (33) reveals that the ratio of the amplitude of the numerically computed waves to that of the analytical waves equals unity. Therefore the PPWE scheme exhibits perfect analytical damping for the waves which propagate. A comparison of the computed to analytic wave celerity is possible by examining the ratio of these two wave speeds which are respectively given by the real parts of equations (32) and (33) and dividing each by

A primitive pseudo wave equation formulation." J. J. Westerink et al. >t-¢..)

a. This celerity ratio equals:

o

x/l(sinaAx]2_(

2

.._I U.l > l,.iJ

)2

1.2

72

(34) I _J

t2_ OE

-72

O3 >el

This celerity ratio is plotted in Fig. 1 for all propagating waves for both the semi-lumped and fully lumped formulations. Both show a phase lag, although the accuracy of the semi-lumped scheme is superior. Consideration of the shortest possible waves, those of the grid wavelength with L = 2Ax or aAx = ~, reduces one of the roots of equation (32) to co = 0. Hence the shortest waves are neither propagated nor damped. This behaviour is the same as that experienced by the primitive equation formulations. In fact equation (32) and thus (34) are almost identical to those obtained for primitive harmonic formulations and are the same when comparing the lumped primitive scheme to the lumped PPWE formulation. Therefore despite the similarity in appearances between the resulting discretized equations for the PPWE and WE formulations, it is clear that the PPWE is basically a primitive formulation. The underlying reason for this is the less accurate centered spatial approximation of the second derivative which resulted in the pseudo wave equation. A factor which significantly influences the structure of the second spatial derivative term is the lumping of the I~¢lroT matrix which is used for both the semi- and fully lumped ^schemes. A consistent scheme which lumps neither Mror nor M, leads to an approximation for this second derivative which is distributed over an entire range of points adjacent to nodej. In fact for a very long grid with no significant boundary influence at node j, the pseudo wave equation and momentum equation are expressed as:

6

(/co) 2 + ico

/

9h(ao;li+~:~tak(~;-k+i'lj+k))=O

i,>

,

,,,

0.'~

0.2

0 (->

l'rl

I i I"/," ff

I~J" //

1

C

S_L

~

IVy// 2

X

__0.0

- -

6 I0 20 60 PARTS PER WAVELENGTH (L/Ax)

g -216

-288

I00

Fig. 1. Ratio of numerical wave speed to analytical wave speed for the PPWE formulation

where •"

.

2/ka6x5

B = 2 aksm [ ~ )

(38)

k=l

Comparison to equation (33) shows that the consistent PPWE scheme also exhibits perfect analytical damping for propagating waves. Furthermore, when a sufficient number of terms are considered in the expansion for B, one of the two roots of equation (37) will again be zero for grid scale wavelengths. Thus the consistent PPWE scheme is also classified as a primitive scheme like the semi- and fully lumped PPWE schemes despite the weak involvement of the j + l and other nodes in the approximation of the second spatial derivative in the pseudo wave equation. The celerity ratio for the consistent PPWE scheme is expressed as:

(35) (o'Ax) 2

(io)+~)(Oj-,+40j+Uj+,]_

CONSISTENT SEIMLUMp_E_D__ OLEY-L UMPED

/

(39)

/~j+ 1 -- 0j- 1_)=

6

)+g~

2Ax

0

(36) where ao = -3.21539, a1=0.43078, a2=1.49227, a3= - 0.39985, a4 = 0.107142, a 5 = - 0.028710, a 6 = 0.007692, a 7 = - 0 . 0 0 2 0 6 1 , as=0.000552, a9=0.000149, etc. The approximation for the second spatial derivative still mainly involves nodes j and j + 2, however nodes j + 1, j___3, j + 4 , etc. also come into play. Repeating the propagation factor analysis for this case gives the following physical roots:

{i). /9t7

B

co=o" 2~Ta---~A (o'Ax)2

2

}

(2ha) 2

(37)

This celerity ratio is plotted in Fig. I and exhibits superior behaviour over both the semi-lumped and fully lumped schemes. A consistent scheme which does not lump Mro7 is however not feasible due to the fact that a fully populated matrix would have to be dealt with when solving (23). Finally we note that using a consistent form of Mror to solve for the primitive momentum equation while using the lumped form of Mror to form the pseudo wave equation should have no effect on the results. It will however make solving for (' a more expensive proposition. • We therefore conclude that within the limits of the Fourier analysis performed, the PPWE formulation appears to be a primitive formulation. However as will

Adv. Water Resources, 1987, Volume 10, December

193

A primitive pseudo wave equation formulation: J. J. Westerink et al.

",~--- 0 PEN OCEAN \

\

BOUNDARY \ \ \

\

\

\ \

\ \ \

½ __!

K

Fig. 2c. r2

Fig. 2a.

Geometry of linear test problem

Fig. 2b.

Coarse grid for linear" test problem

become apparent in the next section, the performance characteristics of the P P W E scheme are more like those of a WE scheme than of a primitive scheme.

4. E V A L U A T I O N O F T H E P E R F O R M A N C E O F THE PPWE FORMULATION In order to verify the accuracy of our P P W E scheme and investigate its numerical behaviour we shall apply it to a number of test cases for which analytical solutions are available. We have selected test cases from a stringent series of test problems developed by Lynch and Gray 3"4"~°. The geometry used for these test problems consists of a quarter of an annulus enclosed with land

194

Refined grid .lot linear test problem

~

Adv. Water Resources, 1987, Volume 10, December

boundaries on three sides and open to the ocean on the outer edge as shown in Fig. 2a. Bathymetry can be constant or vary linearly or quadratically for these radially symmetrical cases and linearized bottom friction values may be specified. These two dimensional depth varying test problems give considerable insight into a numerical scheme's grid wavelength oscillatory behaviour and capability to propagate the longer physical waves. The specifications for the test problems were selected to be identical to those used by Lynch and Gray such that a comparison to their results with both primitive schemes and their SIM and WE formulations is possible. The dimensions used for the geometry are an inner radius of r 1= 2 x 106 ft and an outer radius ofr z = 5 × 106 ft. A tidal forcing period of 12.4 hours is used and the linear bottom friction coefficient )~ is specified such that ), = ~h where r = 0.0001 sec-~. The Coriolis factor f is set to zero. Three cases are considered, the first with constant bathymetry at h = 30 ft and the latter two with bathymetry varying with radial distance in a linear fashion between h = 20 ft at r~ to h = 50 ft at r 2 and quadratically between h = 10 ft at r I and h = 62.5 ft at r 2. Two f.e. grids are used for these test runs; the coarse grid shown in Fig. 2b which is identical to that used by Lynch and Gray 3 5 and the refined grid shown in Fig. 2c. The coarse grid is discretized in a radial manner using A r = 5 x 104ft and AO=n/16 whereas the refined grid uses A r = 2 . 5 x 104ft and A0=~/32. Finally zero normal flux boundary conditions are included in the P P W E formulation as natural boundary conditions in the load vector/~,. The ocean boundary conditions will be specified as essential boundary conditions equal to ~)*= 0.10ft at all nodes on this boundary. Application of both the fully and semi-lumped P P W E schemes to the three test cases indeed confirms the Fourier Analysis predictions in that the semi-lumped scheme shows both better accuracy and more controlled oscillations than the fully lumped scheme. Results obtained using the semi-lumped P P W E formulation with the coarse grid for all three cases are shown in Figs 3 through 5 and are plotted with exact analytical solutions

A primitive pseudo wave equation formulation: J. J. Westerink et al. 0.16 --

formulations 5. The degree of oscillatory behaviour experienced with the P P W E formulation is more similar to that of the SIM scheme of Gray and Lynch 3. The spurious oscillations of the P P W E scheme are in general a bit larger than those of the SIM scheme for the constant bathymetry case while they are slightly smaller for the quadratic bathymetry case. The accuracy of the P P W E solutions tends to be somewhat better than that of the SIM scheme for all cases since the SIM method somewhat underpredicts the exact solution while the P P W E method does not. Both the P P W E and SIM schemes exhibit oscillatory behaviour at a significant proportion of the nodes. The harmonic WE scheme on the other hand typically exhibits pronounced oscillations at only one or two nodal points. Comparison of P P W E to Lynch's harmonic WE results 5 indicates that the WE formulation

0.14 Z 0

X

I--

<[ > w d ILl 0.12

• • O.IC

I

NUMERICAL NUMERICAL I

2

MAX MIN

i

\ I

3

X

I

%L

4

O.16 -

R(xlOSft)

Fig. 3a. Elevation solution for constant bathymetry for coarse grid 0.09

O.IZ Z

o_ bILl ..3 W

0.06

•

0.12

U>.(,.)

o,

Ld

>

0.03

I

0.10

I 3

2

I

I 4

1

5

R(xlOSft) • 0

~ 2

NUMERICAL

i 3

i

R(xl05

Fig. 3b.

i 4

Fig. 4a. Elet'ation solution for linear bathymetry Jbr coarse grid

MIN i

] 5

0.06

ft)

-

Velocity solution for constant bathymetry .lbr

coarse grid generated from formulas developed by Lynch and Gray 1°. Velocity values indicated in the figures are absolute values of velocity. The predicted numerical solution at any value of 0 is contained between the m a x i m u m and minimum points plotted. The separation between the m a x i m u m and minimum values at any value of r is representative of the m a x i m u m degree of node to node oscillation. The results shown indicate that our P P W E formulation does not exhibit the severe oscillatory behaviour normally associated with schemes which are based on the primitive shallow water equations. In fact the amplitudes of the node to node oscillations are typically at least an order of magnitude less than those experienced with standard primitive formulations. This can be verified by comparison to results presented for primitive formulations which use both typical time stepping schemes T M ~ and linear harmonic

"~

0.04

Q)

S UJ >

002

/ 0

: i

2

J

l

3

J 4

i

I 5

R(xlOSft)

Fig. 4b.

Velocity solution

for

linear bathymetry ./'or

coarse grid

Adv. Water Resources, 1987, Volume 10, December

195

A primitive pseudo wave equation formulation: J. J. Westerink et al. 0'16 I

~ O.1~ Z

o_ I.--> W ...I W

O.12

,,

• O102

NUMERICAL MAX

4

3

5

R(xlOSft) Fig. 5a. Elevation solution for quadratic bathymetry for coarse grid 0.06 -

"s 0.04 u--

v >.I'-£..)

q

I

W

>

0.02 • • 02

,,

,

3

NUMERICAL MAX

NUMERICAL MIN ,

,

4

,

I 5

R(xlO-~ ftl

Fig. 5b. Velocity solution for quadratic bathymetry fiyr coarse grid is indeed superior. When comparing for corresponding cases the worst oscillation resulting in the P P W E computations to the worst oscillation experienced with the harmonic WE formulation, the amplitude of the elevation oscillation is between 1 and 4 times larger for the P P W E scheme and the amplitude of the velocity oscillations are between 2 and 3 times larger. It is noted from Figs 3 through 5 that velocities on the inner boundary at r=r~ do not identically equal zero. This is due to our formulation's treatment of normal flux as a natural boundary condition. As was described in Section 2, our formulation combines the integrated weighted interior error of the continuity equation with the integrated weighted error on the normal flux boundary and forces this sum to vanish. Hence there will be some error on elevation and velocity both in the interior domain (Q
196

Adv. Water Resources, 1987, Volume 10, December

boundary (at r = r 1 and 0 - 0, 0 = n/2). The degree of error in velocity for the constant depth case is the same on the inner boundary (at r = rl ) as in the interior domain while this error seems somewhat larger than interior errors for the linearly and quadratically varying bottom cases, However it is noted that all three cases exhibit approximately the same level of inner boundary error in terms of flux, which is the pertinent variable involved in the natural boundary error term. Furthermore the degree of flux error on the inner boundary is consistent with the degree of flux error exhibited in the interior domain and on the elevation prescribed boundary. The inner boundary error is in all cases the worst at the corner nodes. The degree of residual error and of the spurious oscillations on this boundary significantly improves at nodes other than these corner nodes. Finally the residual error of the velocities/fluxes at the corner nodes and at other nodes on the inner boundary are not entirely radially oriented. entirely radially oriented. Thus the quality of our solution is identical on both the natural boundary and in the interior domain. The exact solution will be progressively better satisfied both in the interior and on the natural boundary as grid size is refined. This is indeed shown in Figs 6 through 8 which represent the three cases simulated with the more refined grid. Both flux and velocity residual errors on the natural boundary are now identical to those in the interior domain. The solutions are demonstrated to be consistent and the degree of node to node oscillation is even further reduced from the coarse grid computations. The worst oscillations experienced with the refined grid for both elevation and velocity are reduced by a factor of between approximately 2 to 4 when compared to coarse grid results. This reduction factor is progressively higher for the more rapidly varying depth cases. This is in fact the explanation for the improvement in computed velocity errors on the natural boundaries such that they are now comparable to interior domain velocity errors. Furthermore when comparing all nodes instead of simply comparing the nodes with the worst oscillations, the average reduction in elevation oscillations is approximately equal to 3 while the average reduction in velocity oscillations is well over a factor of 5 for all depth varying cases. A further consequence of treating normal flux as a natural boundary condition is that the solutions on alternating radial rows are not identical as is the case with schemes which use the more traditional approach of satisfying normal flux exactly (i.e., solutions at 0 = 0, 2A0, 4A0 etc. are all identical to each other as are solutions at 0 = A 0 , 3AO, 5A0 etc.). This is due to the fact that the introduction of elevation prescribed essential boundary conditions into our final system of equations (23) does not introduce a sufficient number of constraints to force this behaviour. Figs 3 through 8 therefore indicate the m a x i m u m and minimum solution values at the nodes over the range of 0 values at a given value of r. Solutions are however symmetrical about 0=~z/4. A result of the solutions not having identical values at alternating radial rows is that velocities are not perfectly radial, although the direction of our computed velocities typically matches the direction of the radial line on which the given node lies to an average of less than 0.5 degrees for the coarse grid and less than 5 × 10 -5 degrees for the refined grid. These

A primitive pseudo wave equation formulation." J. J. Westerink et al. 0.16

Oscillations for velocity remain at about the same level as nonclamped P P W E formulations. This then suggests than an important part of the success of the PPWE scheme lies in its treatment of flux prescribed boundary conditions as natural boundary conditions. In fact, implementation of the entirely standard solution procedure of solving the continuity equation (13) and the momentum equation (14) simultaneously (i.e., without the formation of the pseudo wave equation) with elevation and velocity being strictly enforced as essential boundary conditions in each respective equation leads to the identical poor quality solution as the previously considered clamped flux arrangement. Thus this makes it clear that the PPWE scheme's control over oscillations stems entirely from the boundary condition treatment and not from the formation of the pseudo wave equation.

- -

0.14 Z 0 I.< hi J LIJ

0.12

• •

NUMERICAL MAX NUMERICAL MIN

"=, \

0.10

I

I

I

I

3

2

1

"~

4

5

0.16--

R(xlOSft)

Fig. 6a. Elevation solution for constant bathymetry for refined grid O.IZ

0.09 Z 0 I-W J hi

"O

0.12

0.06 >I--

•

q UJ >

,.

O.IC

0.03

2

NUMERICAL MAX

3

%.

\

4

5

R(xlOSft) ,l~

02

• A I

I

3

Fig. 7a. Elevation solution Jor lineal" bathymetry for refined grid

NUMERICAL MAX NUMERICAL MIN I

R(xlOSft)

I

4

I

I

5

0.06 -

Fig. 6b. Velocity solution for constant bathymetry jor refined grid

averages include the inner boundary nodes where deviations from the radial direction of the residual error velocities/fluxes were the most severe. We can override the specification of the flux boundary conditions through the load vector P , by enforcing them strictly in the traditional manner by simply replacing equations in system (14) by the exact specified land boundary velocities after having performed the necessary rotations. Now the substitution and solution steps of the P P W E scheme are implemented as previously described. Even though clamping the velocities on the land boundaries in this fashion yields the exact solution there, it does not in general lead to very satisfactory results. Nodal oscillations for elevation in both the radial and angular (0) directions become much worse and in fact reminiscent of standard primitive formulation results.

0.04 >.I-£.)

q W

>

0.02

/

• NUMERICMAX AL

¢

•

1

I

3

NUMERICAL MtN i

i

R(xlOSft)

i

4

I

5

Fig. 7b. Velocity solution for linear bathymetry for refined grid

Adv. Water Resources, 1987, Volume 10, December

197

A primitive pseudo wave equation formulation: J. J. Westerink et al. 016 I

~

0.I

0.12

~

-

I

~

O, IC

EXACT

I

~ MAX

I 3

I

I 4

I

R(xlOSft) Fig. 8a. Elevation solution for quadratic bathymetry for refined grid 0.06 --

"0 0.04

>I£3

S >

0.02

I

EXACT

• • O

J

2

I

NUMERICAL NUMERICAL I

3

[ 4

MAX MtN I

I

5

R(xlOSft)

Fig. 8b. Velocity solution for quadratic bathymetry for refined grid However the formation of the pseudo wave equation is intimately related to this particular boundary condition treatment and furthermore is responsible for making the overall solution much more economical. 5. C O N C L U S I O N S The development of the PPWE scheme indicates that it is possible to obtain good solutions exhibiting low spurious oscillations with numerical formulations based on the discretized form of the primitive harmonic shallow water equations. Rearranging the primitive discretized equations leads to a final set of discretized equations which at first sight are deceptively similar to those produced by the harmonic WE formulation. Despite this similarity, Fourier analysis indicates that the P P W E formulation is still fundamentally primitive in nature due

198

Adv. Water Resources, 1987, Volume 10, December

to the less accurate approximation of the second spatial derivative of elevation which results in the pseudo wave equation. This approximation involves nodes which are two nodes removed from the node at which the derivatives are evaluated instead of the adjacent nodes. However the performance of the PPWE scheme is distinctly superior to that of other primitive formulations. The P P W E scheme exhibits excellent control over short wavelength noise while not damping the fundamental solution. The extent of the spurious oscillations is similar to that of the SIM method but is not quite as well controlled as that of the WE formulation. Theequations for the PPWE scheme are however in general much simpler than those for the harmonic WE scheme. The PPWE scheme's similarity in oscillatory behaviour to the standard temporal SIM scheme does not stem from the P P W E scheme being any more closely related to the SIM scheme than to the WE scheme since the latter two schemes would be essentially the same for a harmonic formulation. In fact for the standard time domain SIM formulation, control over spurious oscillations is based solely on the second spatial derivative and it is exactly in this regard that the PPWE scheme fails to match the harmonic WE and hypothetical harmonic SIM schemes. In fact the PPWE scheme's equation re-arrangement in of itself is not a requirement for its better controlled oscillatory behaviour although it does significantly increase the efficiency of the solution procedure. Numerical tests indicate that relaxing the velocity specified boundary conditions and treating them as natural boundary conditions is responsible for the well controlled behaviour of the PPWE formulation. This boundary condition treatment goes hand in hand with the formation of the pseudo wave equation which has elevation as its basic variable. Even though these boundary conditions do lead to residual flux/velocity errors on the natural boundaries, the extent of error in terms of flux is on the same order as in the interior domain and on the essential boundaries. Furthermore this error on the natural boundary diminishes with interior errors as the grid is refined. Nonzero normal fluxes on land boundaries, of course, can be potentially troublesome when a circulation model is used in conjunction with a transport model. This potential difficulty can be readily overcome by simply setting these residual nonzero normal land boundary velocities/fluxes to zero only after the computations have been completed. (Recall that trying to enforce zero normal flux as part of the solution resulted in oscillations in elevation.t Doing so will not introduce additional errors locally since these residual normal fluxes are consistent with the flux errors throughout the domain. Neither will this zeroing out of computed residual velocities on land boundaries cause global mass balance problems since over a tidal cycle the nonzero harmonic boundary velocities/fluxes always integrate out to zero. Of course when considering the nonlinear shallow water equations with nonlinear finite amplitude terms included as part of the flux terms, a slight change would be seen in global mass conservation due to this land boundary flux zeroing procedure. However, in general we do not expect the boundary condition treatment of the PPWE scheme to lead to problems when the resulting circulation model is coupled to a transport model, as long as the land boundary fluxes are zeroed after completing the computations.

A primitive pseudo wave equation formulation: J. J. Westerink et al. T h e w o r k i n g v e r s i o n of t h e P P W E s c h e m e is b a s e d o n t h e s e m i - l u m p e d f o r m u l a t i o n w h i c h t a k e s a d v a n t a g e of b e t t e r p h a s e b e h a v i o u r t h a n t h e fully l u m p e d v e r s i o n b u t still m a i n t a i n s b a n d e d a n d / o r d i a g o n a l m a t r i c e s for t h e e q u a t i o n s w h i c h m u s t be s o l v e d . T h i s p a r t i a l l u m p i n g a n d the fact t h a t t h e s o l u t i o n s for e l e v a t i o n a n d v e l o c i t y h a v e been decoupled through the equation re-arrangements makes the PPWE scheme very economical.

4 5 6 7 8

REFERENCES Wang, J. D. and Connor, J. J. Mathematical modeling of near coastal circulation, R.M. Parsons Laboratory Report No. 200, Massachusetts Institute of Technology, Cambridge, MA, 1975 Kawahara. M., Takeuchi, N. and Yoshida, T. Two step explicit finite element method for tsunami wave propagation analysis, Int. J. Num. Meth. Eng., 1978, 12, 331-351 Gray, W. G. and Lynch, D. R. On the control of noise in finite element tidal computations: a semi-implicit approach, Computers and Fluids, 1979, 7, 47-67

9 10 11

Lynch, D. R. and Gray, W. G. A wave equation model for finite element tidal computations. Computers and Fluids. 1979.7,207 228 Lynch, D. Comparison of spectral and time-stepping approaches for finite element modeling of tidal circulation, Oceans 81, Boston, Massachusetts, 1981, 810-814 Lynch, D. R. Finite element solution of the shallot' water equations, PhD Thesis, Department of Civil Engineering, Princeton University, 1978 Le Provost. C. and Rougier, G. Numerical modeling of the harmonic constituents of the tides, with application to the English Channel, J. Phys. Oceanoyr., 1981, I1, 1123-1138 Westerink. J. J.. Connor, J. J. and Stolzenbach, K. D. A frequency-time domain finite element model for tidal circulation based on the least squares harmonic analysis method, Int. J. Num. Meth. Fluids, 1987 tin press) Pinder, G. F. and Gray, W. G. Finite element simulation in surface and subsurface hydrology, Academic Press. N e t York. 1977 Lynch, D. R. and Gray, W. G. Analytical solutions for computer flow model testing, J. Hydraulics Dit'., ASCE, 1978, HYI0, 1409 1428 Gray, W. G. Some inadequacies of finite element models as simulators of two-dimensional circulation, Adv. 14"~ter Resources. 1982.5, 171-177

Adv. Water Resources, 1987, Volume 10, December

199