- Email: [email protected]
On no-slip boundary conditions for the incompressible Navier-Stokes equations
Dynamics of Atmospheres and Oceans, 12 (1988) 47-70
47
Elsevier Science Publishers B.V., Amsterdam - Printed in The Netherlands
O N NO-SLIP B O U N D A R Y C O N D I T I O N S FOR THE INCOMPRESSIBLE N A V I E R - S T O K E S E Q U A T I O N S
M.G.G. F O R E M A N
Institute of Ocean Sciences, P.O. Box 6000, Sidney B.C., V8L 4B2 (Canada) A.F. BENNETT
College of Oceanography, Oregon State University, Corvallis, OR, 97331 (U.S.A.) (Received June 19, 1987; revised February 25, 1988; accepted March 4, 1988)
ABSTRACT Foreman, M.G.G. and Bennett, A.F., 1988. On no-slip boundary conditions for the incompressible Navier-Stokes equations, Dyn. Atmos. Oceans, 12: 47-70. Two local implementations of no-slip boundary conditions are investigated for both the vorticity-streamfunction and momentum-pressure formulations of the time-dependent planar incompressible Navier-Stokes equations, as applied to barotropic ocean circulation modelling. The objective is to determine the extent to which the local accuracy and numerical consistency of these conditions affects the global solution. The effects of a non-local implementation of no-slip conditions for the vorticity-streamfunction equations are also studied. In all cases, boundary condition effects are measured by comparing time-averaged dynamics of turbulent solutions of numerical models based on the two formulations. In the model interior, the energy and enstrophy conserving Arakawa Jacobian is used for the vorticity-streamfunction equations while an extension of the energy and potential enstrophy conserving Arakawa and Lamb finite difference scheme is used for the momentum-pressure equations. Numerical experiments performed with a non-linear model similar to Bryan's barotropic ocean reveal no significant differences between the time-averaged solutions obtained with either of the two formulations, with each using either of the two local boundary conditions. A simple one-dimensional analogue of the vorticity-streamfunction equations is solved algebraically to explain the experimental results. A similar analogue suggests that an apparent inconsistency in the no-slip boundary conditions within the Cox stratified, primitive equation, ocean circulation model should not affect the accuracy or convergence of the global solution.
1. I N T R O D U C T I O N In a numerical integration boundary conditions are used 0377-0265/88/$03.50
of the Navier-Stokes equations, no-slip to evaluate rates of diffusion. Depending
© 1988 Elsevier Science Publishers B.V.
48 upon the formulation of the dynamics, that is the use of either m o m e n t u m and pressure equations or streamfunction and vorticity equations, the diffusing quantity is m o m e n t u m or vorticity, respectively. In two dimensions, the streamfunction and vorticity equations are easier and cheaper to solve numerically and are usually the preferred formulation. However, each formulation has a different implementation of the no-slip condition. With the momentum-pressure equations, the no-slip condition is implemented locally by setting both velocity components to zero. On the other hand, with the streamfunction-vorticity formulation the no-slip condition is usually implemented locally by setting the streamfunction to zero and calculating a wall vorticity. The classical calculation of this vorticity is attributed to T h o m (1928) by Roache (1976). One objective of this study is to assess the effects of these different implementations. The local no-slip boundary condition for m o m e n t u m leads to an exact non-local formula for the wall vorticity. This formula is an area integral of the product of the vorticity field with the derivative of a Green's function. (An equivalent integral has been developed by Quartapelle and Valz-Gris (1981).) A no-slip implementation based on this formula is exactly equivalent to solving the m o m e n t u m - p r e s s u r e equations, subject to the natural, local no-slip boundary conditions for m o m e n t u m . However, alternative local methods for estimating the wall vorticity (such as the T h o m (1928) condition) are available for equations that are discrete in time. (Orszag and Israeli (1974) discuss both explicit and implicit implementations.) In this study we assess the effects of local versus non-local estimates of the wall vorticity by comparing results from the two formulations of the equations of motion, with local boundary conditions used in both cases. Although Quartapelle and Valz-Gris (1981) emphasize that local m o m e n t u m boundary conditions do not admit an exact reformulation in terms of local expressions for wall vorticity, our time-averaged numerical results suggest that the local vorticity expressions are sufficiently accurate. A second issue is the accuracy of the boundary conditions. For both formulations, the local numerical implementation of the no-slip condition often results in less accurate approximations on (or near) the boundary than in the domain interior. For example, the T h o m (1928) condition is only first-order accurate whereas the interior vorticity is usually calculated with a second-order discrete Laplacian operator. Here we study boundary condition accuracy for ocean circulation models where one might suspect that the calculations of intense flows adjacent to the western boundary may be highly dependent on the implementation of the no-slip condition. However, our experimental results indicate that the local numerical no-slip boundary conditions can be less accurate without affecting the model results. We show algebraically for a simple one-dimensional problem why this is the case.
49 A third issue is the numerical consistency (for definition, see for example, Richtmyer and Morton, 1967) of local implementations of the no-slip condition. We consider two particular local approximations to the wall vorticity for the quasi-geostrophic (henceforth QG) equations. They are the first-order Thom (1928) estimate, and a third-order estimate. Although the T h o m estimate seems to produce a numerically inconsistent approximation to Laplacian vorticity diffusion adjacent to the boundary, we show experimentally that the global solution converges and produces the same timeaveraged results as the third-order estimate. Using a simple one-dimensional analogue of the QG equations, we also show convergence algebraically. Another simple one-dimensional example suggests that an even more severe apparent numerical inconsistency in the approximation of diffusion adjacent to a no-slip boundary within the Cox (1984) stratified, primitive equation, ocean circulation model, does not affect the accuracy, nor preclude the convergence, of the global solution. Two local implementations of no-slip are also considered for the momentum and pressure formulation of the dynamics. The first implementation, though numerically inconsistent, is often employed when the primitive equations are solved on an Arakawa C lattice (e.g., the MAC method discussed in Roache (1976)). The second implementation is both consistent and locally more accurate. However, as with the QG models, both implementations seemed to converge and yield essentially the same time-averaged results. In the model interior, the Arakawa Jacobian (Arakawa, 1966) is used for the vorticity-streamfunction equations and an extension of the Arakawa and Lamb (1981) finite difference scheme is used for the m o m e n t u m - p r e s sure equations. In the absence of diffusion and time-stepping, the former conserves energy and enstrophy while the latter conserves energy and potential enstrophy. Mesinger and Arakawa (1976) show that conservation of enstrophy avoids a false cascade of energy to short wavelengths and avoids the aliasing that Phillips (1956) has demonstrated can lead to non-linear instabilities. For a domain with doubly periodic boundaries and a constant depth, it can be shown that the conservation of potential enstrophy implies the conservation of enstrophy. Even though our models do not have doubly periodic boundaries, the similar properties of these two schemes encourages a comparison of the relative merits and accuracy of the two equation formulations. Our preliminary tests with the new local boundary approximations showed little qualitative difference with the results published by Bryan (1963), and Blanford (1971), both of whom used the T h o m estimate. In fact, for the linear cases our streamfunction snap-shots at times coincident with those selected by the previous investigators were virtually identical to their results.
50 However, when we tried to replicate the more highly non-linear cases, model snap-shots did differ significantly when we employed comparable grid resolution. (Our Poisson solver did not allow complete freedom in the n u m b e r of grid intervals across the domain.) In order to gauge the precise effects of these differences, a more systematic study of the diffusion terms adjacent to a no-slip boundary was undertaken. In particular, the accuracy of the boundary conditions was measured by comparing time-averaged dynamics and monitoring the balance of terms in the equations. A modified version of a highly non-linear ocean model first solved numerically by Bryan (1963) was chosen for this study. The numerical demonstrations are carried out for two-dimensional flow but there is no obstacle, in principle, to their application to three-dimensional flow. Our interest is in numerical modelling of ocean circulation, and our dynamics include a latitude-dependent Coriolis parameter. While this affects the circulation pattern, leading to intense b o u n d a r y layers where accurate treatment of no-slip is of particular importance, the Coriolis effect itself does not influence the implementation of any of the b o u n d a r y conditions. 2. GOVERNING EQUATIONS The governing equations are the incompressible Navier-Stokes equations as applied to ocean modelling. A constant depth, a constant diffusion parameter, and a beta-plane approximation to the Coriolis parameter are assumed. The continuity equation is
v-,=0
(1)
and the non-dimensional m o m e n t u m ~v
-- +cqk×hv+cVE= ~t
equation is
E
--V2v
+ F
(2)
Re
where v(x, y, t ) = (u, v ) = velocity; E ( x , y, t ) = energy per unit mass; q( x, y, t) = potential vorticity, = [ f + ( O v / O x ) - ( O u / ~ y ) ] / h ; h = depth; f ( Y ) =f0 + BY -- non-dimensional Coriolis parameter; k = unit vector in the z direction; e = Rossby number; R e = Reynolds number; F(x, y ) = wind stress. The non-dimensional scaling for these equations is consistent with Bryan's (1963) Q G equations. In particular, the Coriolis parameter is scaled by (E/~L) -1, where L = 5000 km is the width of Bryan's ocean, and /~ = 10 -1~ cm -1 s -1. These equations are solved on a square domain with no-slip conditions u= 0 o=0
(3a) (3b)
51
on the boundaries x = 0, x = 1, and free-slip conditions v=0
(4a)
0 u/3y = 0
(4b)
on the boundaries y = 0, y = 1. Conditions 3a and 4a permit no normal flow through the boundary. They are c o m m o n l y referred to as the rigid boundary conditions. Condition 4b is a no-shear condition on the tangential velocity. In order to solve eqs. 1 and 2 numerically, eq. 1 is replaced with a Poisson equation for E. Specifically, taking X7 • eq. 2 and substituting eq. 1 yields cV 2 E = - c v - q k × h v + V ' F
(5)
N e u m a n n boundary conditions for this equation are obtained by taking n- eq. 2 and substituting v = 0. (n is the unit vector normal to the boundary.) They are cn. VE=
n.
('
Re
xy2v + F
t
(6)
The non-dimensionalized Q G equations (see Bryan (1963)) associated with eqs. 1 and 2 are -- +cJ(~,
3t
~')+--=k-v
3x
XF+--V2~
Re
"
= v 2+
(7)
(8)
where q~(x, y, t ) = transport streamfunction; ~'(x, y, t ) = vorticity; J(q~, f ) = J a c o b i a n of streamfunction and vorticity, =3+/Ox 3 f / 3 y -
3q~/ay a~/Ox. The no-slip b o u n d a r y conditions are represented as += 0
(9a)
~b/3x = 0
(9b)
whereas the free-slip b o u n d a r y conditions become + = 0 02+/0y
(10a) 2 = ~" = 0
(10b)
The test model chosen for our studies is a modification of a model first solved numerically by Bryan (1963). In order to conserve computer storage and to permit model runs with three different resolutions, our model d o m a i n was chosen to be the northern half of Bryan's rectangular domain. That is, the domain is 5000 k m x 5000 kin. Parameters Re and c are set to 100 and 0.00032, respectively. The wind stress and all other Q O parameters have the
52 same values as in Bryan (1963). For the primitive equation (PE) model, we set dimensional f0 = 0.32 × 10 -4 s -1. 3. THE NUMERICAL SOLUTIONS 3.1. The PE solution The primitive equations are solved numerically by integrating eq. 2 forward in time from a state of no motion. The Arakawa and Lamb (1981) finite difference scheme (henceforth AL) is extended to include diffusion terms, and applied to the spatial derivatives in eq. 2. The AL assumes an Arakawa C stencil for the u, v, and E variables. It was originally developed to improve the simulation of the non-linear aspects of meteorological flow over steep topography. Although no topography is included in this model, AL has other assets. Numerical experiments by Henry (personal communication, 1985) have demonstrated that AL also provides a more accurate representation of the advective terms in barotropic flow than the popular Richardson-Sielecki (Sielecki, 1968) finite difference scheme. As already mentioned, AL also conserves both energy and potential enstrophy. Equation 2 is integrated in time with a modified leapfrog method which has the diffusion terms lagged by one time step to ensure stability. In order to eliminate time splitting of the solution, a Robert (1966) filter with a parameter value of 0.01 is used. Asselin (1972) has analysed the properties of this filter. The solution procedure is to step forward eq. 2 then solve eq. 5 with a fast direct Poisson routine developed by Temperton (personal communication, 1987). This exact direct solver avoids the computational difficulties often encountered with iterative techniques (e.g., see Roache, 1976). Care is taken to maintain discrete non-divergence and to ensure that the finite difference approximation to eq. 5 is consistent with AL. This is accomplished by applying a discrete divergence operator to the AL versions of eq. 2 and substituting a centered finite difference approximation of eq. 1. The precise numerical implementation of boundary condition 6 is crucial for numerical stability (Foreman and Bennett, 1988). A successful implementation involves assuming a virtual E point half a grid interval beyond the true boundary and setting the normal velocity component to zero in the appropriate discretization of eq. 2. In effect, this means that a E / a In I is approximated on the physical boundary rather than over the outermost E points. This virtual point is then eliminated from the calculations by substituting for it into the discretization of eq. 5 adjacent to a boundary. With an Arakawa C lattice, boundaries are traditionally chosen to lie along the normal velocity components. For our models, this means that the
53 normal boundary conditions within 3 and 4 have a straightforward numerical implementation. Specifically on the eastern and western boundaries, we set u = 0, whereas on the northern and southern boundaries, we set v = 0. However, the tangential velocities are located half a grid interval in from the model boundaries thereby necessitating an indirect implementation of the tangential velocity boundary conditions. In this case, the tangential velocity conditions are implemented indirectly through the calculation of the diffusion term and the AL potential vorticities. In particular, along the northern and southern free-slip boundaries, the potential vorticity is calculated exactly as q = f / h . Denoting by Vl/2,~ and VM_I/2,j, the tangential velocities half a grid interval in from the western and eastern no-slip boundaries, the potential vorticities along these boundaries are calculated as q0o = (fj + 2 V , / z , J A x ) / h qv.j = ( ~ -
2VM_I/Z,j/Ax)/h
(lla) (llb)
respectively. Both these approximations are O(Ax) accurate. The diffusion terms in the equation for the tangential component of velocity require special treatment when they are adjacent to both a free-slip and no-slip boundary. Denoting by Ui,1/2 and Ri,u_l/2, the tangential velocities half a grid interval in from the southern and northern free-slip boundaries, the 02u/0y 2 approximations within V 2u are approximated as 2 Ui,1/2
bli,3/2 -- Ui,1/2
0y 2
Ay 2
~2bli,N-1/2
Ui,N- 3/2 -- Ui,N-1/2
0y 2
Ay 2
(12a) (12b)
respectively. Both these approximations are O(Ay) accurate and are derived using condition 4b. Adjacent to the eastern and western no-slip boundaries, care is required to avoid an inconsistent one-sided approximation to V2•; i.e., a finite difference approximation that does not approach the Laplacian as the grid interval approaches zero. A c o m m o n inconsistent approximation is derived as follows (e.g., see Roache, 1976, pp. 198-199). Consider the no-slip boundary at x = 0 and the term Ozv/Ox 2 within V =v in the v m o m e n t u m equation. Since Oo,s = 0, linear extrapolation to the virtual point v_ 1/2,j yields U 1 / 2 , j = --U1/2, j + O ( A x 2)
(13)
54 Making this substitution into the usual three point approximation to O2u1/2.j/Ox2 at v-1/2,j, namely into ~2U1/2, j -
-
~X 2
=
U3/2, j "Jr U_l/2, j -- 2U1/2, j
(14)
Ax 2
yields 2U1/2, j -
-
U3/2, j --
301/2, j
=
(15)
0X 2
AX 2
However, a Taylor expansion about va/2,j reveals that U3/2, j -- 301/2, j
Ax 2
=
U3/2, j -- 3U1/2, j + 2t)O, j
Ax 2
So eq. 15 certainly does not converge consistent approximation is -
0X 2
3 ~201/2,j 4
3X 2
+O(Ax)
(16)
8201/2,j/SX 2 as AX ~ 0. The correct
4
82Ul/2,j -
to
-
3U3/2,j -- 4U1/2, j =
+O(Ax)
AX 2
(17)
An analogous approximation is used for the 82v/Ox 2 term adjacent to the eastern boundary. Approximations 11, 12 and 17 along the no-slip and free-slip boundaries are O ( A x ) accurate whereas in the domain interior all the AL spatial differences are O(Ax 2) accurate. In a subsequent study we will measure the effects of increasing the accuracy of 11, 12 and 17 to O(Ax2). However, the focus of present PE model tests will be assess the impact of boundary conditions 17 and 15. These tests will determine the effects of an inconsistent condition that uses only 75% of the correct tangential velocity diffusion adjacent to a no-slip boundary.
3.2. The QG solution The QG equations are solved numerically in a manner similar to that employed by Bryan (1963). The major difference is the replacement of the Jacobian J+×, which only conserves energy, with the Arawaka (1966) Jacobian that conserves both energy and enstrophy. As with the PE solution, leapfrog time-stepping and a Robert (1966) time filter with a parameter value of 0.01 are employed. The Dirichlet-Poisson equation for ~k is also solved with a Temperton (personal communication, 1987) direct method. The free-slip northern and southern boundary conditions are implemented by setting both ~ = 0 and ~"= 0. Along the no-slip eastern and western boundaries, we set ~k= 0 and use 8 ~ / 8 x = 0 (condition 9b) in the
55 calculation of the boundary vorticity. This vorticity is in turn used in the calculation of the diffusion and Jacobian terms one grid interval in from the boundary. The classical local calculation of vorticity along a no-slip boundary is attributed to Thorn (1928, 1933) by Roache (1976). The method has been widely used (e.g., Bryan, 1963; Blanford, 1971) with apparent success. However, it is only O(Ax) accurate and this may lead to inaccuracies in the intense flows of the western boundary region. Beardsley (1971) developed a first-order, net vorticity conserving estimate of the boundary vorticity that is equivalent to the T h o m estimate on a Cartesian grid with Ax = Ay. He compared this estimate with the T h o m and (second-order) Pearson (1965) approximations for the linear problem of viscous decay of two-dimensional circular motion and found that the vorticity conserving scheme was most accurate. However, it remains to be seen how well the T h o m scheme performs when advection is included in the governing equations. Thom (1928) estimates the vorticity ~'w at a no-slip wall in terms of the streamfunction values qJw, at the wall, and ~w+l, one grid interval in from the wall, as 2(+w+1 -fw=
An 2
Cw)
An -
3
~)3¢,~o On 3
(18)
where An is the grid interval distance normal to the wall and ~00 is a point within 2xn of the wall. However, ~'w is only required by the conventional five point finite difference approximation to V 2fw+l and the Arakawa Jacobian J(~k~+l, ~',o+1), so it need not be calculated per se. Were it substituted directly into these approximations, both of which have An 2 as their divisor~ the truncation errors for the new estimates would become O ( A n - l ) . It therefore seems that eq. 18 leads to inconsistent numerical approximations of vorticity diffusion and advection adjacent to a no-slip boundary. It is shown later for a one-dimensional analogue of the QG equations (see Discussion) that inconsistency does occur when the Thorn estimate is combined directly with diffusion adjacent to a no-slip boundary. However, it is also shown that the global solution is convergent! Though somewhat surprising, these results do not violate the Lax Equivalence Theorem (see Richtmyer and Morton, 1967, p. 45) which says 'given a properly posed initial-value problem and a finite difference approximation to it that satisfies the consistency condition, stability is the necessary and sufficient condition for convergence'. Indeed, Bryan (1963) reports that for a linear ocean model, the Thom condition produces steady numerical solutions which apparently converge to the steady analytical solution as Ax ---, 0. However it remains to be seen if convergence still holds for approximations to time-dependent non-linear motions, and whether the first-order accuracy of fw does affect
56 the global accuracy of the solution. In order to determine this, eq. 18 will be compared with the more accurate local approximation ~'w=
6~bw+l - ~3 w + 2 + ~kw+3 2 + O(An 3) An 2
(19)
which assumes ~kw= 0 and is derived using Taylor series expansions and eq. 9. Although Roache (1975) found that second-order estimates for ~'w improved the accuracy of diffusion-dominated steady flows, such estimates were not tested here as they would also lead to an apparent (O(1)) inconsistency when combined directly with the conventional second-order finite difference diffusion operator adjacent to a no-slip boundary. The alternative to a local wall vorticity boundary condition is the exact non-local (yet difficult) condition ~4,(x)
~G(x, x I ) - [
Olnl
"D
= o
(20)
Olnl
where x lies on the boundary of the domain D, and G(x, x') is the Green's function for V 2~b= 0
(21a)
with ~b= 0
(21b)
as the boundary conditions. Quartapelle and Valz-Gris (1981) emphasize the need to satisfy an integral equation equivalent to eq. 20 in order to maintain complete equivalence with boundary conditions that are expressed in terms of the velocity. However, they do not assess the inaccuracy resulting from a violation of this strict equivalence. Our subsequent numerical experiments demonstrate that the vorticity boundary conditions, though inexact, yield essentially identical time-averaged results to those obtained with exact no-slip conditions. Boundary conditions aside, the PE and QG models should permit (or more appropriately, exclude) the same wave solutions. Specifically the non-divergence condition (eq. 1) ensures that the PE model does not permit gravity waves. And by redefining E ' = E - fq~
(22)
the Coriolis term can be removed from eq. 2. Consequently, the PE model should not allow inertial oscillations. In order that divergent waves not be permitted in the numerical solutions as well as the analytic solutions, it is important to maintain discrete non-divergence. In particular, were care not
57
taken to ensure that the finite difference approximation to eq. 5 is consistent with the m o m e n t u m equations (as described previously), some type of numerically divergent wave could be generated. It is important to note that the energy, enstrophy, and potential enstrophy conservation for both schemes is no longer valid with the inclusion of diffusion, time-stepping, and b o u n d a r y conditions. Thus, although the PE and QG models should permit exactly the same types of wave solutions, we can only say that the conservation properties of the two finite difference schemes are similar. It will therefore be instructive to see how closely the PE and QG numerical results resemble one another. 4. RESULTS OF M O D E L S I M U L A T I O N S
In order to measure the relative accuracy of the two no-slip b o u n d a r y conditions for the Q G model, three grid resolutions were chosen; namely 129, 257 and 513 points in both the x and y directions. Parameter At was set to 0.25 non-dimensional time units for the coarsest resolution and halved for each successive refinement. The first two model resolutions were started from rest and run for 6000 time units. The 513 resolution models were started at 1500 time units using linearly interpolated values calculated from the respective consistent and inconsistent 257 models and run for a further 4500 time units. In all cases, average values for each of the terms in eq. 7 were calculated over the last 4500 time units. One of these terms, the average Jacobian, was also partitioned as ~ J ( ~ , ~')) = J({~b), {~)) + {J(~b', ~ ' ) )
(23)
where { ) denotes a time average and primes denote deviations from the average values. The first component in eq. 23 is the mean advection of mean vorticity whereas the second is the eddy advection of eddy vorticity in the mean. The numerical stability of the Q G models did not seem to be affected by the two different specifications of the wall vorticity. Both 257 resolution models went unstable (i.e., caused an overflow error) at about the same time step when A t = 0.30, and both were stable for 6000 time units when A t = 0 . 2 5 . (Presumably the instability was due to a violation of the C o u r a n t - F r i e d r i c h s - L e w y condition (see Roache, 1976), but this is difficult to confirm because of the high degree of non-linearity in the problem.) Subsequent experimental runs were m a d e with At = 0.125 for the 257 resolution case and A t = 0.25, 0.0625 for the 129 and 513 resolutions, respectively, as already mentioned. The length of the averaging interval for our model diagnostics was also found to be significant. Specifically the interval must be long enough to
58
contain several eddy life cycles. Denoting b y A and B the Q G models which employ the T h o m condition (eq. 18) and the third-order condition (eq. 19), respectively, results for the 257B model were seen to differ slightly with averaging over 2500 and 4500 time units. However, there were significant differences when the averages were computed over 750 time units. F o r the 513 model, running average values for each of the terms in eq. 7 were calculated after 750, 1 5 0 0 , . . . , 4 5 0 0 time units and the largest and smallest values in the average fields were seen to settle down after - 2250 time units. The choice of At was also found to have some effect on the average values, particularly when the averaging interval is short. Runs with the 257B model again demonstrated as much difference between the average values obtained with At = 0.125 and At = 0.25 over 4500 time units as with At = 0.125 and averages over 2500 and 4500 time units, respectively. Figure 1 shows time-series plots of average kinetic energy density for the 513B model and the difference between the energies of model B and model A. The two models are practically identical for the first 1000 time units suggesting the same transient solution adjustment (via almost linear R o s s b y waves) to the initial conditions. However, dissimilar patterns after 1000 time units indicate a different life cycle of r a n d o m eddies within the model. Since the B - A curve has an approximate time average of zero and oscillations that are comparable in magnitude to those of B, we conclude that the
20
'
I
'
]
k
'
'
~
fB
f
'
10
'
5
15
W 0
__.10 Z
(.9 <
:s
I
-5
0
~
0
I
1
,
I
2
~
I
3 TIME
~
I
4
,
~
5
~ -10
6
Fig. 1. Average kinetic energy density for the 513 QG models. A and B denote the models using boundary conditions 18 and 19, respectively. B - A is the difference between model energies. Time is in thousands of non-dimensional units.
59
\
\\
\\
\
\
\
\
\
.1
/
I L
J
Fig. 2. Time-average streamfunction for QG model 513B. The dashed rectangle shows the window region for the calculation of the four diagnostic quantities and the arrow designates the row for the diagnostic profiles shown in Fig. 4.
different boundary conditions do not affect the time- and space-averaged kinetic energy density in a significant way. Figure 2 shows the average streamfunction contour field for the 513B case. The analogous field for the model A is very similar. Its maximum value is 0.5% smaller and slightly more northward than the maximum for model B. Contour plots were also computed for the four diagnostic quantities (a) -¢(J(~b', ~')), (b) - ¢ J ( ( ~ b ) , (~>), (c) (Ohb/Ox>, and (d) ¢(V2~>/Re, in the subregion of the original domain that is designated by dashed lines in Fig. 2. Average values for the finite difference approximations to each of these four terms were calculated over the last 4500 time units. The streamfunction gradient fields (diagnostic (c)) are almost identical for models A and B. The average diffusion fields (diagnostic (d)) are characterized by a ridge of high values parallel to the western boundary. The model B field has a 4% wider range of values. The fields of mean advection of mean vorticity (diagnostic (b)) have respective ridges of high and low values parallel to the western boundary. Again the patterns are quite similar for the two cases and model B has a 4% wider range. The fields of mean eddy advection of eddy vorticity (diagnostic (a)) for both models A and B
r
E
~-.
~o
~ "
~'~
~'~ ~
~ ~-.~'~ ~.
~
~
~ ~ .. ~~ ~.~.
o
~1~.
~. ~.~
F
~ ~
~.~ ~
~~.
~.~
~
@ C>
61 60
m
40
I
I
I
F
I
:"
-
~ a
w
20
-20
!
-40
0
20
40 60 GRID POINT NUMBER
80
i
100
Fig. 4. Diagnostics (a), (b), (c) and (d) for QG model 513B (solid) and 513A (dashed) across the row shown with an arrow in Fig. 2. Figure 4 shows the four diagnostic quantities for models B and A across the row in the window region that is marked by an arrow in Fig. 2. There are slight differences in diagnostics (a) and (b), but (c) and (d) are virtually identical. Within the western boundary layer of approximately 10 Ax, the balance of terms is clearly (a) plus (c) versus (b) plus (d). But within the next 30 Ax the balance is (a) plus (d) versus (b) plus (c). Contour plots of quantities (a), (b), (c) and (d) for the 129 and 257 resolution cases also show very little difference between models A and B. A comparison with the 513 results shows that the major difference is the resolution of the ridges of high and low values parallel to the western boundary. Specifically, the ridges of high values in (b) and (d) are highest for the 513 resolution case and the ridge of low values in (a) is deepest for this case. This is also seen in profiles analogous to those shown in Fig. 4. Figure 5a and b show the diagnostic (a) contour fields for models 257A and 257B. Numerical solutions of the same test problem in PE form are remarkably similar to those for the Q G solution. Owing to the larger storage requirements of the AL scheme, only two resolutions were run for the PE model; namely 129 and 257 points in both the x and £ directions. Parameter At was set to 0.25 and 0.125 non-dimensional units, respectively (the same as for the Q G model). Again all models were run for 6000 time units with average values calculated over the last 4500 time units. In order to permit comparisons with the Q G model, the vorticity, streamfunction, and diagnostic quantities (b) (the mean advection of m e a n vorticity) and (a) (the eddy
62
(a)
(b)
Fig. 5. Time-average eddy advection of eddy vorticity (diagnostic (a)) for the QG models (a) 257A and (b) 257B. Contour levels range from - 4 0 to 20 in increments of 10. H and L denote the maximum and minimum values, respectively. a d v e c t i o n o f e d d y v o r t i c i t y in the m e a n ) were c a l c u l a t e d f r o m the P E m o d e l variables u a n d v. F i g u r e 6 shows the time-series plots o f spatially a v e r a g e d kinetic e n e r g y d e n s i t y for the 257 p o i n t P E m o d e l s with d i f f u s i o n t e r m s 15 a n d 17 ( d e n o t e d
63 20
I
I
'
I
m
I
10
I
iD 15
5
UJ 0
2~ fi3 z
.~,,D-C
--4
0 ctD m
< I
-5
0 [
0
,
I
1
~
I
2
~
I
3 TIME
~
h
4
t
5
-10
,
6
Fig. 6, Average kinetic energy density for the 257 PE models, C and D denote the models using diffusion terms 15 and 17, respectively, adjacent to no-slip boundaries. D - C is the difference between the model energies. Time is in thousands of non-dimensional units.
Fig. 7. Time-average streamfunction for PE model with diffusion term 17 adjacent to the no-slip boundaries.
64 as m o d e l s C a n d D, respectively), a d j a c e n t to the n o - s l i p b o u n d a r i e s . N o t i c e t h a t b o t h P E m o d e l s s h o w m o r e distinct l o n g - p e r i o d o s c i l l a t i o n s t h a t their Q G m o d e l c o u n t e r p a r t s a n d t h a t the d i f f e r e n c e in energies b e t w e e n D a n d C
ill
( Cb
L
(a)
(b)
Fig. 8. Time-average eddy advection of eddy vorticity (diagnostic (a)) for the 257 point PE models with (a) diffusion term 15, and (b) diffusion term 17. Contour levels range from - 4 0 to 20 in increments of 10 for (a) and - 4 0 to 30 in increments of 10 for (b). H and L denote the maximum and minimum values, respectively.
65
is approximately zero for a longer time than the difference between B and A. However, the time-averaged kinetic energy density and the shorter period oscillations are approximately the same for all models. Figure 7 shows the contour field of average streamfunction for model 257D. Its maximum value is - 1% larger than the maximum for 513B in Fig. 2, and, although the contour levels are very close, there are some subtle differences. The analogous fields for model C has a maximum value that is - 3% smaller than for model D and a slightly different contour around the local low that is eastward from the western ridge of high values. The fields of eddy advection of eddy vorticity in the mean (diagnostic (a)) are shown in Fig. 8 and are again characterized by ridges of respective low and high values parallel to the western boundary and are quite similar to those for the QG model. However, the range of values for model D is 7% larger than that for model C. Also, the range of values for model D is - 3% smaller than that for the 257 point QG model B. The fields of mean advection of mean vorticity are also similar. However, model D has an 8% larger range of values than model C. 5. D I S C U S S I O N
In order to understand the similar QG model results obtained with the Thom condition and the locally more accurate and consistent condition 19, the following simple problem was studied. Consider the one-dimensional steady state QG equations 2= k
(24a)
~" = d Z ~ b / d x 2
(24b)
d2~/dx
with constant forcing k, and boundary conditions = d+/dx
(24c)
= 0
at both x -~ 0 and x = 1. The solution is +
(25a) (25b)
=
and the Thom boundary condition (eq. 18) is indeed only first-order accurate; that is (26a)
~(0) = ~k
2qJ(Ax) Ax 2
=
+ O( A x)
(26b)
66 The solution to the associated numerical problem (27a)
=k
Ax 2
(27b) ~'j =
~x 2
with boundary conditions f,0=0
(27c)
~2N = 0
(27d)
0f,] ~'0- Ax 2
(27e)
O+2N-1 ~2N -- - -
(27f)
Ax 2
for 0 a non-zero constant is
~,= '2k(j ~ x - ~)2+ ~ k ~x ~ + +j=
~4k(j A x -
½)4+
½a(j Ax-
(28a) ½)2+ fi
(28b)
with -(1/12)k{(1/2)
- [(3/2) - (3/0)1 Ax + 2 Ax 2 - Ax3[1 - ( 2 / 0 ) ] }
O~
1 -[(0-
2 ) / 0 ] Ax
(2Sc) fi=-~s
a+
(28d)
Substituting conditions 27e and 27b into eq. 27a and performing Taylor expansions about Xl, the diffusion equation adjacent to the b o u n d a r y becomes ~b3 - 4~2 + (5 + 0 ) ~ 1 AX 4
-
0
dt~l
Ax 3 dx
1 + 0 / 2 d2+1 Ax 2
+ O ( A x -1 )
(29)
dx 2
This estimate is not consistent with d4~bl/dx 4. Consequently, the T h o m estimate (0 = 2) is inconsistent when combined with diffusion adjacent to a no-slip boundary. However, comparing the numerical and analytical solutions, we see that +j and ~j are O(Ax 2) accurate when 0 = 2 and O(Ax) accurate when 0 4= 0 and 0 4= 2. So even though the T h o m condition is only first-order accurate,
67 the global numerical solution is second-order accurate. (Gustafsson (1975) proved a more general result for hyperbolic systems of equations, namely that boundary and initial approximations may be one order less accurate than the interior approximations without decreasing the overall accuracy.) Furthermore, when 0 4= 0 and 0 4= 2, conditions 27e and 27f are only O(1) accurate and hence numerically inconsistent. But the corresponding global numerical solution is O(Ax) and therefore convergent! The reason behind this somewhat surprising result is that the boundary conditions do not affect diffusion at the boundary. In particular, the boundary conditions only influence the complementary function parameters and ft. They do not affect the particular solution
~f" = ½k(j Ax-- ½)2 + ~2k Ax 2
(30)
which is second-order accurate because of the term ~ k Ax 2. Since a a n d / 3 are annihilated by the discrete diffusion operator acting on fj, only fff is diffused at the boundary. Thus the boundary conditions are irrelevant to the convergence of the diffusion approximation at the boundary. When conditions 27e and 27f are replaced by the third-order condition (eq. 19), both a and/~ become third-order accurate: otherwise eqs. 28a and 28b are unchanged. Hence ~j remains second-order accurate while ~/ becomes third-order accurate. We might then expect slightly more accurate results were eq. 19 to be used instead of eq. 18. However, results of the non-linear two-dimensional QG experiments described in the previous section suggest that there is very little difference in the mean dynamics of the models run with the two conditions. In fact it seems that any set of boundary conditions that yields an a which is at least second-order accurate will, because of the component ~ k Ax 2 within ~'7, cause the global solution ~'j to be second-order accurate. A similar simple example can be constructed to demonstrate the effects of using eq. 15 versus eq. 17 for the PE equations. In particular, the seemingly inconsistent condition 15 can be shown to produce a numerical solution that is globally second-order accurate. The widely used Cox (1984) stratified, primitive equation, ocean circulation model also has an apparent numerically inconsistent treatment of diffusion adjacent to a no-slip boundary. However, as with the Thorn condition for the QG equations, this inconsistency does not seem to affect the numerical results. In fact a simple one-dimensional version of the governing equations shows that again, the accuracy of the global solution is not affected by the apparent inconsistency at the boundary. In the Cox model, Laplacian diffusion for the depth-integrated transport leads to a biharmonic operator acting on the vertically integrated stream-
68 function, ~/,. In one dimension the finite difference approximation to this biharmonic diffusion
~4@i ~i--2 -- 4~*--1 + 61/Ji- 4~i+1 + ~i+2 V41~i - ~X 4 = mx 4
(31)
is O(Ax 2) accurate. As the Cox model employs an Arakawa B grid with the velocities lying along the boundary, the 'boundary' q, point, denoted by +1/2, actually lies ½kx in from the boundary. A no-slip boundary at point x 0 is then approximated by setting the two q, points that straddle the boundary to zero, i.e.
~1/2 = ~-1/2 = 0
(32)
Substituting 32 into eq. 31, performing Taylor expansions about x3/2, and using the true boundary conditions +0 = 0
(33a)
3q,o/aX = 0
(33b)
we see that the diffusion approximation adjacent to the boundary is
6~b3/2 -- 4~b5/2 ~ 7 /42 A+X
-- 8 Ax 2 3 0 2 ~)3 / 23 1+ O2 ( ~xl
(34)
So the numerical approximation is certainly not consistent with 04~3/2/0X 4. However, the Cox (1984) model has been used successfully in many applications and has had no reports of divergent solutions as A x - , 0. (Although to our knowledge, this limit has not been pushed; two or three points in the western boundary layer is common.) A simple example, similar to the one for the QG equations reveals that again the apparent numerical inconsistency at the boundary does not affect the convergence of the global solution. Consider the steady state problem
V4~ = k
(35)
subject to the boundary conditions 24c at both x = 0 and x = 1. The solution is given by eq. 25a and the Cox boundary conditions +(½Ax) = q~(- ½Ax) = 0
(36)
are accurate to O(Ax2). The solution to the associated numerical problem
~j+2 -- 4~j+1 + 6q,j - 4~bj_ 1 + ~j-2 =k Ax 4
(37a)
69 with boundary conditions •1/2 = ~--1/2 = + N - l / 2 = ~ N + l / 2 = 0
on a grid that staggers the boundary (i.e., j 28b but now k
-
(37b) 1 3 ...) is again given by eq. 5,5
(Ax 2 + 1)
(38a)
(2xx 2 - 1) 2
(38b)
24 13 =
k 384
So both a and/3 are accurate to O( Ax 2). Consequently the global numerical solution is also accurate to O(Ax 2) and thus convergent. 6. CONCLUSIONS The preceding analysis and numerical tests have shown little sensitivity to the numerical implementation of four no-slip boundary conditions for the incompressible Navier-Stokes equations. For the QG equations it was seen that the Thom local calculation of the wall vorticity, though seemingly inconsistent when combined with Laplacian diffusion adjacent to a no-slip boundary, not only produces a convergent global solution but also is as accurate as a locally third-order estimate. Furthermore, essentially identical time-averaged results for the QG and PE models indicate that the local implementation of the no-slip conditions for the QG equations, though not equivalent to the non-local implementation (Quartapelle and Valz-Gris, 1981), is sufficiently accurate. Our tests also confirmed that the energy and potential enstrophy conserving A r a k a w a - L a m b (1981) finite difference solution of the PE equations gives the same time-averaged results as the energy and enstrophy conserving finite difference solution of the quasi-geostrophic equations with an Arakawa (1966) Jacobian. Although an extensive analysis of no-slip boundary conditions for the widely used Cox (1984) ocean circulation model was not performed, a simple one-dimensional example suggests that here also, an apparent numerical inconsistency in the treatment of diffusion adjacent to a no-slip boundaries does not affect the accuracy nor preclude the convergence of the global solution. ACKNOWLEDGEMENTS We thank the reviewers for their constructive criticism of an earlier version of this paper.
70 REFERENCES Arakawa, A., 1966. Computational design for long-term numerical integration of equations of fluid motion: two-dimensional incompressible flow. Part 1. J. Comput. Phys., 1: 119-143. Arakawa, A. and Lamb, V.R., 1981. A potential enstrophy and energy conserving scheme for the shallow water equations. Mon. Weath. Rev., 109: 18-36. Asselin, R., 1972. Frequency filter for time integrations. Mon. Weath. Rev., 100: 487-490. Beardsley, R.C., 1971. Note on finite-difference approximations for the viscous generation of vorticity at a boundary. Mon. Weath. Rev., 99: 387-392. Blanford, R.R., 1971. Boundary conditions in homogeneous ocean models. Deep-Sea Res., 18: 739-751. Bryan, K., 1963. A numerical investigation of a nonlinear model of a wind-driven ocean. J. Atmos. Sci., 20: 594-606. Cox, M.D., 1984. A Primitive Equation Three-Dimensional Model of the Ocean. G F D L Ocean Group Technical Report No. 1, G F D L / N O A A , Princeton University, Princeton, 250 pp. Foreman, M.G.G. and Bennett, A.F., 1988. A numerical stability analysis of the two-dimensional incompressible Navier-Stokes equations. J. Comput. Phys. (in press). Gustafsson, B., 1975. The convergence rate for difference approximations to mixed initial boundary value problems. Math. Computation, 29: 396-406. Mesinger, F. and Arakawa, A., 1976. Numerical Methods Used in Atmospheric Models. Vol. 1, W M O - I C S U Joint Organizing Committee, G A R P Publication Series No. 17, Geneva, 64 pp. Orszag, S.A. and Israeli, M., 1974. Numerical simulation of viscous incompressible flows. In: M. van Dyke, W.G. Vincenti and J.V. Wehausen (Editors), Annual Review of Fluid Mechanics, Vol. 6. Annual Reviews Inc., Palo Alto, pp. 281-318. Pearson, C.E., 1965. A computational method for viscous flow problems. J. Fluid Mech., 21: 611-622. Phillips, N.A., 1956. The general circulation of the atmosphere: a numerical experiment. Q. J. R. Meteorol. Soc., 82: 123-164. Quartapelle, L. and Valz-Gris, F., 1981. Projection conditions on the vorticity in viscous incompressible flows. Int. J. Numerical Methods Fluids, 1: 129-144. Richtmyer, R.D. and Morton, K.W., 1967. Difference Methods for Initial-Value Problems. Wiley-Interscience, New York, 405 pp. Roache, P.J., 1975. The LAD, NOS and split NOS methods for the steady-state Navier-Stokes equations. Comput. Fluids, 3: 179-195. Roache, P.J., 1976. Computational Fluid Dynamics. Hermosa Publishers, Albuquerque, 446 PP. Robert, A.J. 1966. The integration of a low order spectral form of the primitive meteorological equations. J. Meteorol. Soc. Jpn., 44: 237-245. Sielecki, A., 1968. An energy-conserving difference scheme for the storm surge equations. Mon. Weath. Rev., 96: 150-156. Thom, A., 1928. An Investigation of Fluid Flow in Two Dimensions. Aerospace Research Centre, R. and M., No. 1194, United Kingdom. Thom, A., 1933. The flow past circular cylinders at low speeds. Proc. R. Soc. London, Set. A, 141: 651-666.
47
Elsevier Science Publishers B.V., Amsterdam - Printed in The Netherlands
O N NO-SLIP B O U N D A R Y C O N D I T I O N S FOR THE INCOMPRESSIBLE N A V I E R - S T O K E S E Q U A T I O N S
M.G.G. F O R E M A N
Institute of Ocean Sciences, P.O. Box 6000, Sidney B.C., V8L 4B2 (Canada) A.F. BENNETT
College of Oceanography, Oregon State University, Corvallis, OR, 97331 (U.S.A.) (Received June 19, 1987; revised February 25, 1988; accepted March 4, 1988)
ABSTRACT Foreman, M.G.G. and Bennett, A.F., 1988. On no-slip boundary conditions for the incompressible Navier-Stokes equations, Dyn. Atmos. Oceans, 12: 47-70. Two local implementations of no-slip boundary conditions are investigated for both the vorticity-streamfunction and momentum-pressure formulations of the time-dependent planar incompressible Navier-Stokes equations, as applied to barotropic ocean circulation modelling. The objective is to determine the extent to which the local accuracy and numerical consistency of these conditions affects the global solution. The effects of a non-local implementation of no-slip conditions for the vorticity-streamfunction equations are also studied. In all cases, boundary condition effects are measured by comparing time-averaged dynamics of turbulent solutions of numerical models based on the two formulations. In the model interior, the energy and enstrophy conserving Arakawa Jacobian is used for the vorticity-streamfunction equations while an extension of the energy and potential enstrophy conserving Arakawa and Lamb finite difference scheme is used for the momentum-pressure equations. Numerical experiments performed with a non-linear model similar to Bryan's barotropic ocean reveal no significant differences between the time-averaged solutions obtained with either of the two formulations, with each using either of the two local boundary conditions. A simple one-dimensional analogue of the vorticity-streamfunction equations is solved algebraically to explain the experimental results. A similar analogue suggests that an apparent inconsistency in the no-slip boundary conditions within the Cox stratified, primitive equation, ocean circulation model should not affect the accuracy or convergence of the global solution.
1. I N T R O D U C T I O N In a numerical integration boundary conditions are used 0377-0265/88/$03.50
of the Navier-Stokes equations, no-slip to evaluate rates of diffusion. Depending
© 1988 Elsevier Science Publishers B.V.
48 upon the formulation of the dynamics, that is the use of either m o m e n t u m and pressure equations or streamfunction and vorticity equations, the diffusing quantity is m o m e n t u m or vorticity, respectively. In two dimensions, the streamfunction and vorticity equations are easier and cheaper to solve numerically and are usually the preferred formulation. However, each formulation has a different implementation of the no-slip condition. With the momentum-pressure equations, the no-slip condition is implemented locally by setting both velocity components to zero. On the other hand, with the streamfunction-vorticity formulation the no-slip condition is usually implemented locally by setting the streamfunction to zero and calculating a wall vorticity. The classical calculation of this vorticity is attributed to T h o m (1928) by Roache (1976). One objective of this study is to assess the effects of these different implementations. The local no-slip boundary condition for m o m e n t u m leads to an exact non-local formula for the wall vorticity. This formula is an area integral of the product of the vorticity field with the derivative of a Green's function. (An equivalent integral has been developed by Quartapelle and Valz-Gris (1981).) A no-slip implementation based on this formula is exactly equivalent to solving the m o m e n t u m - p r e s s u r e equations, subject to the natural, local no-slip boundary conditions for m o m e n t u m . However, alternative local methods for estimating the wall vorticity (such as the T h o m (1928) condition) are available for equations that are discrete in time. (Orszag and Israeli (1974) discuss both explicit and implicit implementations.) In this study we assess the effects of local versus non-local estimates of the wall vorticity by comparing results from the two formulations of the equations of motion, with local boundary conditions used in both cases. Although Quartapelle and Valz-Gris (1981) emphasize that local m o m e n t u m boundary conditions do not admit an exact reformulation in terms of local expressions for wall vorticity, our time-averaged numerical results suggest that the local vorticity expressions are sufficiently accurate. A second issue is the accuracy of the boundary conditions. For both formulations, the local numerical implementation of the no-slip condition often results in less accurate approximations on (or near) the boundary than in the domain interior. For example, the T h o m (1928) condition is only first-order accurate whereas the interior vorticity is usually calculated with a second-order discrete Laplacian operator. Here we study boundary condition accuracy for ocean circulation models where one might suspect that the calculations of intense flows adjacent to the western boundary may be highly dependent on the implementation of the no-slip condition. However, our experimental results indicate that the local numerical no-slip boundary conditions can be less accurate without affecting the model results. We show algebraically for a simple one-dimensional problem why this is the case.
49 A third issue is the numerical consistency (for definition, see for example, Richtmyer and Morton, 1967) of local implementations of the no-slip condition. We consider two particular local approximations to the wall vorticity for the quasi-geostrophic (henceforth QG) equations. They are the first-order Thom (1928) estimate, and a third-order estimate. Although the T h o m estimate seems to produce a numerically inconsistent approximation to Laplacian vorticity diffusion adjacent to the boundary, we show experimentally that the global solution converges and produces the same timeaveraged results as the third-order estimate. Using a simple one-dimensional analogue of the QG equations, we also show convergence algebraically. Another simple one-dimensional example suggests that an even more severe apparent numerical inconsistency in the approximation of diffusion adjacent to a no-slip boundary within the Cox (1984) stratified, primitive equation, ocean circulation model, does not affect the accuracy, nor preclude the convergence, of the global solution. Two local implementations of no-slip are also considered for the momentum and pressure formulation of the dynamics. The first implementation, though numerically inconsistent, is often employed when the primitive equations are solved on an Arakawa C lattice (e.g., the MAC method discussed in Roache (1976)). The second implementation is both consistent and locally more accurate. However, as with the QG models, both implementations seemed to converge and yield essentially the same time-averaged results. In the model interior, the Arakawa Jacobian (Arakawa, 1966) is used for the vorticity-streamfunction equations and an extension of the Arakawa and Lamb (1981) finite difference scheme is used for the m o m e n t u m - p r e s sure equations. In the absence of diffusion and time-stepping, the former conserves energy and enstrophy while the latter conserves energy and potential enstrophy. Mesinger and Arakawa (1976) show that conservation of enstrophy avoids a false cascade of energy to short wavelengths and avoids the aliasing that Phillips (1956) has demonstrated can lead to non-linear instabilities. For a domain with doubly periodic boundaries and a constant depth, it can be shown that the conservation of potential enstrophy implies the conservation of enstrophy. Even though our models do not have doubly periodic boundaries, the similar properties of these two schemes encourages a comparison of the relative merits and accuracy of the two equation formulations. Our preliminary tests with the new local boundary approximations showed little qualitative difference with the results published by Bryan (1963), and Blanford (1971), both of whom used the T h o m estimate. In fact, for the linear cases our streamfunction snap-shots at times coincident with those selected by the previous investigators were virtually identical to their results.
50 However, when we tried to replicate the more highly non-linear cases, model snap-shots did differ significantly when we employed comparable grid resolution. (Our Poisson solver did not allow complete freedom in the n u m b e r of grid intervals across the domain.) In order to gauge the precise effects of these differences, a more systematic study of the diffusion terms adjacent to a no-slip boundary was undertaken. In particular, the accuracy of the boundary conditions was measured by comparing time-averaged dynamics and monitoring the balance of terms in the equations. A modified version of a highly non-linear ocean model first solved numerically by Bryan (1963) was chosen for this study. The numerical demonstrations are carried out for two-dimensional flow but there is no obstacle, in principle, to their application to three-dimensional flow. Our interest is in numerical modelling of ocean circulation, and our dynamics include a latitude-dependent Coriolis parameter. While this affects the circulation pattern, leading to intense b o u n d a r y layers where accurate treatment of no-slip is of particular importance, the Coriolis effect itself does not influence the implementation of any of the b o u n d a r y conditions. 2. GOVERNING EQUATIONS The governing equations are the incompressible Navier-Stokes equations as applied to ocean modelling. A constant depth, a constant diffusion parameter, and a beta-plane approximation to the Coriolis parameter are assumed. The continuity equation is
v-,=0
(1)
and the non-dimensional m o m e n t u m ~v
-- +cqk×hv+cVE= ~t
equation is
E
--V2v
+ F
(2)
Re
where v(x, y, t ) = (u, v ) = velocity; E ( x , y, t ) = energy per unit mass; q( x, y, t) = potential vorticity, = [ f + ( O v / O x ) - ( O u / ~ y ) ] / h ; h = depth; f ( Y ) =f0 + BY -- non-dimensional Coriolis parameter; k = unit vector in the z direction; e = Rossby number; R e = Reynolds number; F(x, y ) = wind stress. The non-dimensional scaling for these equations is consistent with Bryan's (1963) Q G equations. In particular, the Coriolis parameter is scaled by (E/~L) -1, where L = 5000 km is the width of Bryan's ocean, and /~ = 10 -1~ cm -1 s -1. These equations are solved on a square domain with no-slip conditions u= 0 o=0
(3a) (3b)
51
on the boundaries x = 0, x = 1, and free-slip conditions v=0
(4a)
0 u/3y = 0
(4b)
on the boundaries y = 0, y = 1. Conditions 3a and 4a permit no normal flow through the boundary. They are c o m m o n l y referred to as the rigid boundary conditions. Condition 4b is a no-shear condition on the tangential velocity. In order to solve eqs. 1 and 2 numerically, eq. 1 is replaced with a Poisson equation for E. Specifically, taking X7 • eq. 2 and substituting eq. 1 yields cV 2 E = - c v - q k × h v + V ' F
(5)
N e u m a n n boundary conditions for this equation are obtained by taking n- eq. 2 and substituting v = 0. (n is the unit vector normal to the boundary.) They are cn. VE=
n.
('
Re
xy2v + F
t
(6)
The non-dimensionalized Q G equations (see Bryan (1963)) associated with eqs. 1 and 2 are -- +cJ(~,
3t
~')+--=k-v
3x
XF+--V2~
Re
"
= v 2+
(7)
(8)
where q~(x, y, t ) = transport streamfunction; ~'(x, y, t ) = vorticity; J(q~, f ) = J a c o b i a n of streamfunction and vorticity, =3+/Ox 3 f / 3 y -
3q~/ay a~/Ox. The no-slip b o u n d a r y conditions are represented as += 0
(9a)
~b/3x = 0
(9b)
whereas the free-slip b o u n d a r y conditions become + = 0 02+/0y
(10a) 2 = ~" = 0
(10b)
The test model chosen for our studies is a modification of a model first solved numerically by Bryan (1963). In order to conserve computer storage and to permit model runs with three different resolutions, our model d o m a i n was chosen to be the northern half of Bryan's rectangular domain. That is, the domain is 5000 k m x 5000 kin. Parameters Re and c are set to 100 and 0.00032, respectively. The wind stress and all other Q O parameters have the
52 same values as in Bryan (1963). For the primitive equation (PE) model, we set dimensional f0 = 0.32 × 10 -4 s -1. 3. THE NUMERICAL SOLUTIONS 3.1. The PE solution The primitive equations are solved numerically by integrating eq. 2 forward in time from a state of no motion. The Arakawa and Lamb (1981) finite difference scheme (henceforth AL) is extended to include diffusion terms, and applied to the spatial derivatives in eq. 2. The AL assumes an Arakawa C stencil for the u, v, and E variables. It was originally developed to improve the simulation of the non-linear aspects of meteorological flow over steep topography. Although no topography is included in this model, AL has other assets. Numerical experiments by Henry (personal communication, 1985) have demonstrated that AL also provides a more accurate representation of the advective terms in barotropic flow than the popular Richardson-Sielecki (Sielecki, 1968) finite difference scheme. As already mentioned, AL also conserves both energy and potential enstrophy. Equation 2 is integrated in time with a modified leapfrog method which has the diffusion terms lagged by one time step to ensure stability. In order to eliminate time splitting of the solution, a Robert (1966) filter with a parameter value of 0.01 is used. Asselin (1972) has analysed the properties of this filter. The solution procedure is to step forward eq. 2 then solve eq. 5 with a fast direct Poisson routine developed by Temperton (personal communication, 1987). This exact direct solver avoids the computational difficulties often encountered with iterative techniques (e.g., see Roache, 1976). Care is taken to maintain discrete non-divergence and to ensure that the finite difference approximation to eq. 5 is consistent with AL. This is accomplished by applying a discrete divergence operator to the AL versions of eq. 2 and substituting a centered finite difference approximation of eq. 1. The precise numerical implementation of boundary condition 6 is crucial for numerical stability (Foreman and Bennett, 1988). A successful implementation involves assuming a virtual E point half a grid interval beyond the true boundary and setting the normal velocity component to zero in the appropriate discretization of eq. 2. In effect, this means that a E / a In I is approximated on the physical boundary rather than over the outermost E points. This virtual point is then eliminated from the calculations by substituting for it into the discretization of eq. 5 adjacent to a boundary. With an Arakawa C lattice, boundaries are traditionally chosen to lie along the normal velocity components. For our models, this means that the
53 normal boundary conditions within 3 and 4 have a straightforward numerical implementation. Specifically on the eastern and western boundaries, we set u = 0, whereas on the northern and southern boundaries, we set v = 0. However, the tangential velocities are located half a grid interval in from the model boundaries thereby necessitating an indirect implementation of the tangential velocity boundary conditions. In this case, the tangential velocity conditions are implemented indirectly through the calculation of the diffusion term and the AL potential vorticities. In particular, along the northern and southern free-slip boundaries, the potential vorticity is calculated exactly as q = f / h . Denoting by Vl/2,~ and VM_I/2,j, the tangential velocities half a grid interval in from the western and eastern no-slip boundaries, the potential vorticities along these boundaries are calculated as q0o = (fj + 2 V , / z , J A x ) / h qv.j = ( ~ -
2VM_I/Z,j/Ax)/h
(lla) (llb)
respectively. Both these approximations are O(Ax) accurate. The diffusion terms in the equation for the tangential component of velocity require special treatment when they are adjacent to both a free-slip and no-slip boundary. Denoting by Ui,1/2 and Ri,u_l/2, the tangential velocities half a grid interval in from the southern and northern free-slip boundaries, the 02u/0y 2 approximations within V 2u are approximated as 2 Ui,1/2
bli,3/2 -- Ui,1/2
0y 2
Ay 2
~2bli,N-1/2
Ui,N- 3/2 -- Ui,N-1/2
0y 2
Ay 2
(12a) (12b)
respectively. Both these approximations are O(Ay) accurate and are derived using condition 4b. Adjacent to the eastern and western no-slip boundaries, care is required to avoid an inconsistent one-sided approximation to V2•; i.e., a finite difference approximation that does not approach the Laplacian as the grid interval approaches zero. A c o m m o n inconsistent approximation is derived as follows (e.g., see Roache, 1976, pp. 198-199). Consider the no-slip boundary at x = 0 and the term Ozv/Ox 2 within V =v in the v m o m e n t u m equation. Since Oo,s = 0, linear extrapolation to the virtual point v_ 1/2,j yields U 1 / 2 , j = --U1/2, j + O ( A x 2)
(13)
54 Making this substitution into the usual three point approximation to O2u1/2.j/Ox2 at v-1/2,j, namely into ~2U1/2, j -
-
~X 2
=
U3/2, j "Jr U_l/2, j -- 2U1/2, j
(14)
Ax 2
yields 2U1/2, j -
-
U3/2, j --
301/2, j
=
(15)
0X 2
AX 2
However, a Taylor expansion about va/2,j reveals that U3/2, j -- 301/2, j
Ax 2
=
U3/2, j -- 3U1/2, j + 2t)O, j
Ax 2
So eq. 15 certainly does not converge consistent approximation is -
0X 2
3 ~201/2,j 4
3X 2
+O(Ax)
(16)
8201/2,j/SX 2 as AX ~ 0. The correct
4
82Ul/2,j -
to
-
3U3/2,j -- 4U1/2, j =
+O(Ax)
AX 2
(17)
An analogous approximation is used for the 82v/Ox 2 term adjacent to the eastern boundary. Approximations 11, 12 and 17 along the no-slip and free-slip boundaries are O ( A x ) accurate whereas in the domain interior all the AL spatial differences are O(Ax 2) accurate. In a subsequent study we will measure the effects of increasing the accuracy of 11, 12 and 17 to O(Ax2). However, the focus of present PE model tests will be assess the impact of boundary conditions 17 and 15. These tests will determine the effects of an inconsistent condition that uses only 75% of the correct tangential velocity diffusion adjacent to a no-slip boundary.
3.2. The QG solution The QG equations are solved numerically in a manner similar to that employed by Bryan (1963). The major difference is the replacement of the Jacobian J+×, which only conserves energy, with the Arawaka (1966) Jacobian that conserves both energy and enstrophy. As with the PE solution, leapfrog time-stepping and a Robert (1966) time filter with a parameter value of 0.01 are employed. The Dirichlet-Poisson equation for ~k is also solved with a Temperton (personal communication, 1987) direct method. The free-slip northern and southern boundary conditions are implemented by setting both ~ = 0 and ~"= 0. Along the no-slip eastern and western boundaries, we set ~k= 0 and use 8 ~ / 8 x = 0 (condition 9b) in the
55 calculation of the boundary vorticity. This vorticity is in turn used in the calculation of the diffusion and Jacobian terms one grid interval in from the boundary. The classical local calculation of vorticity along a no-slip boundary is attributed to Thorn (1928, 1933) by Roache (1976). The method has been widely used (e.g., Bryan, 1963; Blanford, 1971) with apparent success. However, it is only O(Ax) accurate and this may lead to inaccuracies in the intense flows of the western boundary region. Beardsley (1971) developed a first-order, net vorticity conserving estimate of the boundary vorticity that is equivalent to the T h o m estimate on a Cartesian grid with Ax = Ay. He compared this estimate with the T h o m and (second-order) Pearson (1965) approximations for the linear problem of viscous decay of two-dimensional circular motion and found that the vorticity conserving scheme was most accurate. However, it remains to be seen how well the T h o m scheme performs when advection is included in the governing equations. Thom (1928) estimates the vorticity ~'w at a no-slip wall in terms of the streamfunction values qJw, at the wall, and ~w+l, one grid interval in from the wall, as 2(+w+1 -fw=
An 2
Cw)
An -
3
~)3¢,~o On 3
(18)
where An is the grid interval distance normal to the wall and ~00 is a point within 2xn of the wall. However, ~'w is only required by the conventional five point finite difference approximation to V 2fw+l and the Arakawa Jacobian J(~k~+l, ~',o+1), so it need not be calculated per se. Were it substituted directly into these approximations, both of which have An 2 as their divisor~ the truncation errors for the new estimates would become O ( A n - l ) . It therefore seems that eq. 18 leads to inconsistent numerical approximations of vorticity diffusion and advection adjacent to a no-slip boundary. It is shown later for a one-dimensional analogue of the QG equations (see Discussion) that inconsistency does occur when the Thorn estimate is combined directly with diffusion adjacent to a no-slip boundary. However, it is also shown that the global solution is convergent! Though somewhat surprising, these results do not violate the Lax Equivalence Theorem (see Richtmyer and Morton, 1967, p. 45) which says 'given a properly posed initial-value problem and a finite difference approximation to it that satisfies the consistency condition, stability is the necessary and sufficient condition for convergence'. Indeed, Bryan (1963) reports that for a linear ocean model, the Thom condition produces steady numerical solutions which apparently converge to the steady analytical solution as Ax ---, 0. However it remains to be seen if convergence still holds for approximations to time-dependent non-linear motions, and whether the first-order accuracy of fw does affect
56 the global accuracy of the solution. In order to determine this, eq. 18 will be compared with the more accurate local approximation ~'w=
6~bw+l - ~3 w + 2 + ~kw+3 2 + O(An 3) An 2
(19)
which assumes ~kw= 0 and is derived using Taylor series expansions and eq. 9. Although Roache (1975) found that second-order estimates for ~'w improved the accuracy of diffusion-dominated steady flows, such estimates were not tested here as they would also lead to an apparent (O(1)) inconsistency when combined directly with the conventional second-order finite difference diffusion operator adjacent to a no-slip boundary. The alternative to a local wall vorticity boundary condition is the exact non-local (yet difficult) condition ~4,(x)
~G(x, x I ) - [
Olnl
"D
= o
(20)
Olnl
where x lies on the boundary of the domain D, and G(x, x') is the Green's function for V 2~b= 0
(21a)
with ~b= 0
(21b)
as the boundary conditions. Quartapelle and Valz-Gris (1981) emphasize the need to satisfy an integral equation equivalent to eq. 20 in order to maintain complete equivalence with boundary conditions that are expressed in terms of the velocity. However, they do not assess the inaccuracy resulting from a violation of this strict equivalence. Our subsequent numerical experiments demonstrate that the vorticity boundary conditions, though inexact, yield essentially identical time-averaged results to those obtained with exact no-slip conditions. Boundary conditions aside, the PE and QG models should permit (or more appropriately, exclude) the same wave solutions. Specifically the non-divergence condition (eq. 1) ensures that the PE model does not permit gravity waves. And by redefining E ' = E - fq~
(22)
the Coriolis term can be removed from eq. 2. Consequently, the PE model should not allow inertial oscillations. In order that divergent waves not be permitted in the numerical solutions as well as the analytic solutions, it is important to maintain discrete non-divergence. In particular, were care not
57
taken to ensure that the finite difference approximation to eq. 5 is consistent with the m o m e n t u m equations (as described previously), some type of numerically divergent wave could be generated. It is important to note that the energy, enstrophy, and potential enstrophy conservation for both schemes is no longer valid with the inclusion of diffusion, time-stepping, and b o u n d a r y conditions. Thus, although the PE and QG models should permit exactly the same types of wave solutions, we can only say that the conservation properties of the two finite difference schemes are similar. It will therefore be instructive to see how closely the PE and QG numerical results resemble one another. 4. RESULTS OF M O D E L S I M U L A T I O N S
In order to measure the relative accuracy of the two no-slip b o u n d a r y conditions for the Q G model, three grid resolutions were chosen; namely 129, 257 and 513 points in both the x and y directions. Parameter At was set to 0.25 non-dimensional time units for the coarsest resolution and halved for each successive refinement. The first two model resolutions were started from rest and run for 6000 time units. The 513 resolution models were started at 1500 time units using linearly interpolated values calculated from the respective consistent and inconsistent 257 models and run for a further 4500 time units. In all cases, average values for each of the terms in eq. 7 were calculated over the last 4500 time units. One of these terms, the average Jacobian, was also partitioned as ~ J ( ~ , ~')) = J({~b), {~)) + {J(~b', ~ ' ) )
(23)
where { ) denotes a time average and primes denote deviations from the average values. The first component in eq. 23 is the mean advection of mean vorticity whereas the second is the eddy advection of eddy vorticity in the mean. The numerical stability of the Q G models did not seem to be affected by the two different specifications of the wall vorticity. Both 257 resolution models went unstable (i.e., caused an overflow error) at about the same time step when A t = 0.30, and both were stable for 6000 time units when A t = 0 . 2 5 . (Presumably the instability was due to a violation of the C o u r a n t - F r i e d r i c h s - L e w y condition (see Roache, 1976), but this is difficult to confirm because of the high degree of non-linearity in the problem.) Subsequent experimental runs were m a d e with At = 0.125 for the 257 resolution case and A t = 0.25, 0.0625 for the 129 and 513 resolutions, respectively, as already mentioned. The length of the averaging interval for our model diagnostics was also found to be significant. Specifically the interval must be long enough to
58
contain several eddy life cycles. Denoting b y A and B the Q G models which employ the T h o m condition (eq. 18) and the third-order condition (eq. 19), respectively, results for the 257B model were seen to differ slightly with averaging over 2500 and 4500 time units. However, there were significant differences when the averages were computed over 750 time units. F o r the 513 model, running average values for each of the terms in eq. 7 were calculated after 750, 1 5 0 0 , . . . , 4 5 0 0 time units and the largest and smallest values in the average fields were seen to settle down after - 2250 time units. The choice of At was also found to have some effect on the average values, particularly when the averaging interval is short. Runs with the 257B model again demonstrated as much difference between the average values obtained with At = 0.125 and At = 0.25 over 4500 time units as with At = 0.125 and averages over 2500 and 4500 time units, respectively. Figure 1 shows time-series plots of average kinetic energy density for the 513B model and the difference between the energies of model B and model A. The two models are practically identical for the first 1000 time units suggesting the same transient solution adjustment (via almost linear R o s s b y waves) to the initial conditions. However, dissimilar patterns after 1000 time units indicate a different life cycle of r a n d o m eddies within the model. Since the B - A curve has an approximate time average of zero and oscillations that are comparable in magnitude to those of B, we conclude that the
20
'
I
'
]
k
'
'
~
fB
f
'
10
'
5
15
W 0
__.10 Z
(.9 <
:s
I
-5
0
~
0
I
1
,
I
2
~
I
3 TIME
~
I
4
,
~
5
~ -10
6
Fig. 1. Average kinetic energy density for the 513 QG models. A and B denote the models using boundary conditions 18 and 19, respectively. B - A is the difference between model energies. Time is in thousands of non-dimensional units.
59
\
\\
\\
\
\
\
\
\
.1
/
I L
J
Fig. 2. Time-average streamfunction for QG model 513B. The dashed rectangle shows the window region for the calculation of the four diagnostic quantities and the arrow designates the row for the diagnostic profiles shown in Fig. 4.
different boundary conditions do not affect the time- and space-averaged kinetic energy density in a significant way. Figure 2 shows the average streamfunction contour field for the 513B case. The analogous field for the model A is very similar. Its maximum value is 0.5% smaller and slightly more northward than the maximum for model B. Contour plots were also computed for the four diagnostic quantities (a) -¢(J(~b', ~')), (b) - ¢ J ( ( ~ b ) , (~>), (c) (Ohb/Ox>, and (d) ¢(V2~>/Re, in the subregion of the original domain that is designated by dashed lines in Fig. 2. Average values for the finite difference approximations to each of these four terms were calculated over the last 4500 time units. The streamfunction gradient fields (diagnostic (c)) are almost identical for models A and B. The average diffusion fields (diagnostic (d)) are characterized by a ridge of high values parallel to the western boundary. The model B field has a 4% wider range of values. The fields of mean advection of mean vorticity (diagnostic (b)) have respective ridges of high and low values parallel to the western boundary. Again the patterns are quite similar for the two cases and model B has a 4% wider range. The fields of mean eddy advection of eddy vorticity (diagnostic (a)) for both models A and B
r
E
~-.
~o
~ "
~'~
~'~ ~
~ ~-.~'~ ~.
~
~
~ ~ .. ~~ ~.~.
o
~1~.
~. ~.~
F
~ ~
~.~ ~
~~.
~.~
~
@ C>
61 60
m
40
I
I
I
F
I
:"
-
~ a
w
20
-20
!
-40
0
20
40 60 GRID POINT NUMBER
80
i
100
Fig. 4. Diagnostics (a), (b), (c) and (d) for QG model 513B (solid) and 513A (dashed) across the row shown with an arrow in Fig. 2. Figure 4 shows the four diagnostic quantities for models B and A across the row in the window region that is marked by an arrow in Fig. 2. There are slight differences in diagnostics (a) and (b), but (c) and (d) are virtually identical. Within the western boundary layer of approximately 10 Ax, the balance of terms is clearly (a) plus (c) versus (b) plus (d). But within the next 30 Ax the balance is (a) plus (d) versus (b) plus (c). Contour plots of quantities (a), (b), (c) and (d) for the 129 and 257 resolution cases also show very little difference between models A and B. A comparison with the 513 results shows that the major difference is the resolution of the ridges of high and low values parallel to the western boundary. Specifically, the ridges of high values in (b) and (d) are highest for the 513 resolution case and the ridge of low values in (a) is deepest for this case. This is also seen in profiles analogous to those shown in Fig. 4. Figure 5a and b show the diagnostic (a) contour fields for models 257A and 257B. Numerical solutions of the same test problem in PE form are remarkably similar to those for the Q G solution. Owing to the larger storage requirements of the AL scheme, only two resolutions were run for the PE model; namely 129 and 257 points in both the x and £ directions. Parameter At was set to 0.25 and 0.125 non-dimensional units, respectively (the same as for the Q G model). Again all models were run for 6000 time units with average values calculated over the last 4500 time units. In order to permit comparisons with the Q G model, the vorticity, streamfunction, and diagnostic quantities (b) (the mean advection of m e a n vorticity) and (a) (the eddy
62
(a)
(b)
Fig. 5. Time-average eddy advection of eddy vorticity (diagnostic (a)) for the QG models (a) 257A and (b) 257B. Contour levels range from - 4 0 to 20 in increments of 10. H and L denote the maximum and minimum values, respectively. a d v e c t i o n o f e d d y v o r t i c i t y in the m e a n ) were c a l c u l a t e d f r o m the P E m o d e l variables u a n d v. F i g u r e 6 shows the time-series plots o f spatially a v e r a g e d kinetic e n e r g y d e n s i t y for the 257 p o i n t P E m o d e l s with d i f f u s i o n t e r m s 15 a n d 17 ( d e n o t e d
63 20
I
I
'
I
m
I
10
I
iD 15
5
UJ 0
2~ fi3 z
.~,,D-C
--4
0 ctD m
< I
-5
0 [
0
,
I
1
~
I
2
~
I
3 TIME
~
h
4
t
5
-10
,
6
Fig. 6, Average kinetic energy density for the 257 PE models, C and D denote the models using diffusion terms 15 and 17, respectively, adjacent to no-slip boundaries. D - C is the difference between the model energies. Time is in thousands of non-dimensional units.
Fig. 7. Time-average streamfunction for PE model with diffusion term 17 adjacent to the no-slip boundaries.
64 as m o d e l s C a n d D, respectively), a d j a c e n t to the n o - s l i p b o u n d a r i e s . N o t i c e t h a t b o t h P E m o d e l s s h o w m o r e distinct l o n g - p e r i o d o s c i l l a t i o n s t h a t their Q G m o d e l c o u n t e r p a r t s a n d t h a t the d i f f e r e n c e in energies b e t w e e n D a n d C
ill
( Cb
L
(a)
(b)
Fig. 8. Time-average eddy advection of eddy vorticity (diagnostic (a)) for the 257 point PE models with (a) diffusion term 15, and (b) diffusion term 17. Contour levels range from - 4 0 to 20 in increments of 10 for (a) and - 4 0 to 30 in increments of 10 for (b). H and L denote the maximum and minimum values, respectively.
65
is approximately zero for a longer time than the difference between B and A. However, the time-averaged kinetic energy density and the shorter period oscillations are approximately the same for all models. Figure 7 shows the contour field of average streamfunction for model 257D. Its maximum value is - 1% larger than the maximum for 513B in Fig. 2, and, although the contour levels are very close, there are some subtle differences. The analogous fields for model C has a maximum value that is - 3% smaller than for model D and a slightly different contour around the local low that is eastward from the western ridge of high values. The fields of eddy advection of eddy vorticity in the mean (diagnostic (a)) are shown in Fig. 8 and are again characterized by ridges of respective low and high values parallel to the western boundary and are quite similar to those for the QG model. However, the range of values for model D is 7% larger than that for model C. Also, the range of values for model D is - 3% smaller than that for the 257 point QG model B. The fields of mean advection of mean vorticity are also similar. However, model D has an 8% larger range of values than model C. 5. D I S C U S S I O N
In order to understand the similar QG model results obtained with the Thom condition and the locally more accurate and consistent condition 19, the following simple problem was studied. Consider the one-dimensional steady state QG equations 2= k
(24a)
~" = d Z ~ b / d x 2
(24b)
d2~/dx
with constant forcing k, and boundary conditions = d+/dx
(24c)
= 0
at both x -~ 0 and x = 1. The solution is +
(25a) (25b)
=
and the Thom boundary condition (eq. 18) is indeed only first-order accurate; that is (26a)
~(0) = ~k
2qJ(Ax) Ax 2
=
+ O( A x)
(26b)
66 The solution to the associated numerical problem (27a)
=k
Ax 2
(27b) ~'j =
~x 2
with boundary conditions f,0=0
(27c)
~2N = 0
(27d)
0f,] ~'0- Ax 2
(27e)
O+2N-1 ~2N -- - -
(27f)
Ax 2
for 0 a non-zero constant is
~,= '2k(j ~ x - ~)2+ ~ k ~x ~ + +j=
~4k(j A x -
½)4+
½a(j Ax-
(28a) ½)2+ fi
(28b)
with -(1/12)k{(1/2)
- [(3/2) - (3/0)1 Ax + 2 Ax 2 - Ax3[1 - ( 2 / 0 ) ] }
O~
1 -[(0-
2 ) / 0 ] Ax
(2Sc) fi=-~s
a+
(28d)
Substituting conditions 27e and 27b into eq. 27a and performing Taylor expansions about Xl, the diffusion equation adjacent to the b o u n d a r y becomes ~b3 - 4~2 + (5 + 0 ) ~ 1 AX 4
-
0
dt~l
Ax 3 dx
1 + 0 / 2 d2+1 Ax 2
+ O ( A x -1 )
(29)
dx 2
This estimate is not consistent with d4~bl/dx 4. Consequently, the T h o m estimate (0 = 2) is inconsistent when combined with diffusion adjacent to a no-slip boundary. However, comparing the numerical and analytical solutions, we see that +j and ~j are O(Ax 2) accurate when 0 = 2 and O(Ax) accurate when 0 4= 0 and 0 4= 2. So even though the T h o m condition is only first-order accurate,
67 the global numerical solution is second-order accurate. (Gustafsson (1975) proved a more general result for hyperbolic systems of equations, namely that boundary and initial approximations may be one order less accurate than the interior approximations without decreasing the overall accuracy.) Furthermore, when 0 4= 0 and 0 4= 2, conditions 27e and 27f are only O(1) accurate and hence numerically inconsistent. But the corresponding global numerical solution is O(Ax) and therefore convergent! The reason behind this somewhat surprising result is that the boundary conditions do not affect diffusion at the boundary. In particular, the boundary conditions only influence the complementary function parameters and ft. They do not affect the particular solution
~f" = ½k(j Ax-- ½)2 + ~2k Ax 2
(30)
which is second-order accurate because of the term ~ k Ax 2. Since a a n d / 3 are annihilated by the discrete diffusion operator acting on fj, only fff is diffused at the boundary. Thus the boundary conditions are irrelevant to the convergence of the diffusion approximation at the boundary. When conditions 27e and 27f are replaced by the third-order condition (eq. 19), both a and/~ become third-order accurate: otherwise eqs. 28a and 28b are unchanged. Hence ~j remains second-order accurate while ~/ becomes third-order accurate. We might then expect slightly more accurate results were eq. 19 to be used instead of eq. 18. However, results of the non-linear two-dimensional QG experiments described in the previous section suggest that there is very little difference in the mean dynamics of the models run with the two conditions. In fact it seems that any set of boundary conditions that yields an a which is at least second-order accurate will, because of the component ~ k Ax 2 within ~'7, cause the global solution ~'j to be second-order accurate. A similar simple example can be constructed to demonstrate the effects of using eq. 15 versus eq. 17 for the PE equations. In particular, the seemingly inconsistent condition 15 can be shown to produce a numerical solution that is globally second-order accurate. The widely used Cox (1984) stratified, primitive equation, ocean circulation model also has an apparent numerically inconsistent treatment of diffusion adjacent to a no-slip boundary. However, as with the Thorn condition for the QG equations, this inconsistency does not seem to affect the numerical results. In fact a simple one-dimensional version of the governing equations shows that again, the accuracy of the global solution is not affected by the apparent inconsistency at the boundary. In the Cox model, Laplacian diffusion for the depth-integrated transport leads to a biharmonic operator acting on the vertically integrated stream-
68 function, ~/,. In one dimension the finite difference approximation to this biharmonic diffusion
~4@i ~i--2 -- 4~*--1 + 61/Ji- 4~i+1 + ~i+2 V41~i - ~X 4 = mx 4
(31)
is O(Ax 2) accurate. As the Cox model employs an Arakawa B grid with the velocities lying along the boundary, the 'boundary' q, point, denoted by +1/2, actually lies ½kx in from the boundary. A no-slip boundary at point x 0 is then approximated by setting the two q, points that straddle the boundary to zero, i.e.
~1/2 = ~-1/2 = 0
(32)
Substituting 32 into eq. 31, performing Taylor expansions about x3/2, and using the true boundary conditions +0 = 0
(33a)
3q,o/aX = 0
(33b)
we see that the diffusion approximation adjacent to the boundary is
6~b3/2 -- 4~b5/2 ~ 7 /42 A+X
-- 8 Ax 2 3 0 2 ~)3 / 23 1+ O2 ( ~xl
(34)
So the numerical approximation is certainly not consistent with 04~3/2/0X 4. However, the Cox (1984) model has been used successfully in many applications and has had no reports of divergent solutions as A x - , 0. (Although to our knowledge, this limit has not been pushed; two or three points in the western boundary layer is common.) A simple example, similar to the one for the QG equations reveals that again the apparent numerical inconsistency at the boundary does not affect the convergence of the global solution. Consider the steady state problem
V4~ = k
(35)
subject to the boundary conditions 24c at both x = 0 and x = 1. The solution is given by eq. 25a and the Cox boundary conditions +(½Ax) = q~(- ½Ax) = 0
(36)
are accurate to O(Ax2). The solution to the associated numerical problem
~j+2 -- 4~j+1 + 6q,j - 4~bj_ 1 + ~j-2 =k Ax 4
(37a)
69 with boundary conditions •1/2 = ~--1/2 = + N - l / 2 = ~ N + l / 2 = 0
on a grid that staggers the boundary (i.e., j 28b but now k
-
(37b) 1 3 ...) is again given by eq. 5,5
(Ax 2 + 1)
(38a)
(2xx 2 - 1) 2
(38b)
24 13 =
k 384
So both a and/3 are accurate to O( Ax 2). Consequently the global numerical solution is also accurate to O(Ax 2) and thus convergent. 6. CONCLUSIONS The preceding analysis and numerical tests have shown little sensitivity to the numerical implementation of four no-slip boundary conditions for the incompressible Navier-Stokes equations. For the QG equations it was seen that the Thom local calculation of the wall vorticity, though seemingly inconsistent when combined with Laplacian diffusion adjacent to a no-slip boundary, not only produces a convergent global solution but also is as accurate as a locally third-order estimate. Furthermore, essentially identical time-averaged results for the QG and PE models indicate that the local implementation of the no-slip conditions for the QG equations, though not equivalent to the non-local implementation (Quartapelle and Valz-Gris, 1981), is sufficiently accurate. Our tests also confirmed that the energy and potential enstrophy conserving A r a k a w a - L a m b (1981) finite difference solution of the PE equations gives the same time-averaged results as the energy and enstrophy conserving finite difference solution of the quasi-geostrophic equations with an Arakawa (1966) Jacobian. Although an extensive analysis of no-slip boundary conditions for the widely used Cox (1984) ocean circulation model was not performed, a simple one-dimensional example suggests that here also, an apparent numerical inconsistency in the treatment of diffusion adjacent to a no-slip boundaries does not affect the accuracy nor preclude the convergence of the global solution. ACKNOWLEDGEMENTS We thank the reviewers for their constructive criticism of an earlier version of this paper.
70 REFERENCES Arakawa, A., 1966. Computational design for long-term numerical integration of equations of fluid motion: two-dimensional incompressible flow. Part 1. J. Comput. Phys., 1: 119-143. Arakawa, A. and Lamb, V.R., 1981. A potential enstrophy and energy conserving scheme for the shallow water equations. Mon. Weath. Rev., 109: 18-36. Asselin, R., 1972. Frequency filter for time integrations. Mon. Weath. Rev., 100: 487-490. Beardsley, R.C., 1971. Note on finite-difference approximations for the viscous generation of vorticity at a boundary. Mon. Weath. Rev., 99: 387-392. Blanford, R.R., 1971. Boundary conditions in homogeneous ocean models. Deep-Sea Res., 18: 739-751. Bryan, K., 1963. A numerical investigation of a nonlinear model of a wind-driven ocean. J. Atmos. Sci., 20: 594-606. Cox, M.D., 1984. A Primitive Equation Three-Dimensional Model of the Ocean. G F D L Ocean Group Technical Report No. 1, G F D L / N O A A , Princeton University, Princeton, 250 pp. Foreman, M.G.G. and Bennett, A.F., 1988. A numerical stability analysis of the two-dimensional incompressible Navier-Stokes equations. J. Comput. Phys. (in press). Gustafsson, B., 1975. The convergence rate for difference approximations to mixed initial boundary value problems. Math. Computation, 29: 396-406. Mesinger, F. and Arakawa, A., 1976. Numerical Methods Used in Atmospheric Models. Vol. 1, W M O - I C S U Joint Organizing Committee, G A R P Publication Series No. 17, Geneva, 64 pp. Orszag, S.A. and Israeli, M., 1974. Numerical simulation of viscous incompressible flows. In: M. van Dyke, W.G. Vincenti and J.V. Wehausen (Editors), Annual Review of Fluid Mechanics, Vol. 6. Annual Reviews Inc., Palo Alto, pp. 281-318. Pearson, C.E., 1965. A computational method for viscous flow problems. J. Fluid Mech., 21: 611-622. Phillips, N.A., 1956. The general circulation of the atmosphere: a numerical experiment. Q. J. R. Meteorol. Soc., 82: 123-164. Quartapelle, L. and Valz-Gris, F., 1981. Projection conditions on the vorticity in viscous incompressible flows. Int. J. Numerical Methods Fluids, 1: 129-144. Richtmyer, R.D. and Morton, K.W., 1967. Difference Methods for Initial-Value Problems. Wiley-Interscience, New York, 405 pp. Roache, P.J., 1975. The LAD, NOS and split NOS methods for the steady-state Navier-Stokes equations. Comput. Fluids, 3: 179-195. Roache, P.J., 1976. Computational Fluid Dynamics. Hermosa Publishers, Albuquerque, 446 PP. Robert, A.J. 1966. The integration of a low order spectral form of the primitive meteorological equations. J. Meteorol. Soc. Jpn., 44: 237-245. Sielecki, A., 1968. An energy-conserving difference scheme for the storm surge equations. Mon. Weath. Rev., 96: 150-156. Thom, A., 1928. An Investigation of Fluid Flow in Two Dimensions. Aerospace Research Centre, R. and M., No. 1194, United Kingdom. Thom, A., 1933. The flow past circular cylinders at low speeds. Proc. R. Soc. London, Set. A, 141: 651-666.