- Email: [email protected]
Quiet: A reduced noise finite element model for tidal circulation
QUIET: A reduced noise finite element model for tidal circulation WILLIAM G. GRAY and INGEMAR P. E. KINNMARK Water Resources Program, Department o f CtvdEngmeermg, Prmceton, New Jersey 08544, USA
The model described herein uses the vertically averaged wave equation and m o m e n t u m equations to simulate two-dimensional tidal flow. The model uses quadratic isoparametric elements in space and a new explicit in time scheme to compute surface elevations and velocities which are essentially noise-free. The model is shown to be efficient in that it requires the solution o f no matrices. A sample application to the southern North Sea demonstrates the utility of the model.
y-direction 3V aV OV --+U--+V--+g-3t ~x ay
a(H-- h) + fU + rV= Ay]H ay
(2b)
or
conservative x-dtrectlon 3(HU)
a(HUU) t - - - - t
ax
~t
D(HUV) -
-
a(H-- h) + g H - -
ay
~x
(3a)
-- f H V + rHU = A x
INTRODUCTION A QUadratic Isoparametnc Exphclt in Tmae (QUIET) finite element model has been developed and documented The model is based on the wave equation formulation of the continuity equaUon descnbed by Lynch and Gray ~ coupled with the freedom to select various time weighting schemes for the momentum equations The QUIET model is noteworthy in that computational noise is suppressed wtule the physical characteristics are not damped. The model also does not require the solution of any matnces. In this report, the computational features of QUIET wtuch account for its speed and versatility are developed Data requirements are also dtscussed, and an example simulation of the southern part of the North Sea is presented.
EQUATION FORMULATION The continuity and momentum equations used when modelling vertacaUy well-mLxed essentially honzontal tidal circulation, are obtained by integrating the point equations through the depth of the flow. A clear derivation of the equatmns is gwen by Prltchard 2 Following ttus work we obtmn Continuity 3H --
+
3t
a(HU)+ a(UV) = ax ay
0
(1)
y-direction a(HUV) + a(HVV)
3(HV) ~t
ax
a ( H - h)
~ gH - -
~y
ay
(3b)
+ f H U + r H V = Ay where
is positive eastward is posttwe northward Y is tnne t u , v are the vertically averaged velocities in the x- and y-directions respectively is the total depth of flow H h xs bathymetry is the Conohs parameter f ts the non-hnear bottom friction 7 g is gravity, and Ax, Ay are the atmospheric wind stresses. X
Most finite element models of estuarlne flow solve equation (1) In conjunction with either the non-conservatwe momentum equations (2) or the conservative forms (3). Lynch and Gray 1 have proposed using a wave equation in place of equation (1) because of some advantageous numerical properties. The wave equation is derived by first differentiating equation (1) with respect to tame to obtain"
Momentum equations
+
non-conservatxve x-direction bU 8U --+U--+ 3t ~x
~U V--+g--
ay
a(H-- h) ax
Next equations (3) are used to replace. f V + rU = A x / H
(2a)
Received August 1982 Written dlscussmn closes September 1983.
130
(4)
Adv Eng Software 1983, Vol 5, No 3
a~u) bt
and
a(HV) Ot
0141-1195/83/030130-07 $2 00 ©1983CMLPubhcatlons
where
and obtain:
b2n at 2
el' (r,b)= IIr.b
b [ b(HUU)+ b(HUI0 t- gH b(H -- h) bx [ bx by bx
-my
+
by
J
(5)
Finally equation (1) is used to replace:
and ~2 is the areal regmn of interest. To facilitate the computations, mtegralaon by parts (Green's theorem) is applied to the third bracketed term m equation (9) and the integral around the boundary of the domain ,s re-expressed using the conservatwe momentum equation. The result of these operations is
<
02H
--r
+
by
J
by
\
aH + f x H v --A], V~b,) --
so that the final form of the wave equation is.
a~H
--+~ at 2
OZ
O rO(HU03 a(HUV)
at
ax
I--
[
Ox
O(H -- h) I - - - + g H - by Ox
y[°7)+,.q
+
aT
-- HV--
ay
h)
z= 1,M
b2H
OH
0t 2
0t
< ( ~ + rv + f x v), ep,)+ ( g V ( , - h ) , ~t)
(6)
The computer model developed solves equations (2) and (6) for H, U and V To make their mampulatton easier in the subsequent portmn of the paper, we can re-write these equations using vector notation:
---I- r----V
. IV. (Hw) + gItV(H--h)+ f x h v - - A ]
-- H r . V r = 0
(7)
av A --+ v.Vv+gV(H--h) +fxv+ rv----= 0 (8) 0t H Dlscretisation of these equalaons is accomphshed by applying the fimte element method in the spatial domam along with a finite chfference marching procedure in the time domain. These aspects of the model will be discussed in some detail so that the reasons for its speed and low computer core requirements are revealed.
+<(v.vv--A),
Applicataon of the finite element method to the wave and momentum equations has been presented previously by Lynch and Gray. 1 Here we will bnefly demonstrate their results and concentrate on the selection of the bass function type, quadrature scheme, and time stepping algonthm. Use of the fimte element method in space requires that the differential equation be orthogonal to each of a set of M basis functions ~bz(x, y):
Ot 2 '
--
i:l,M
#,>= 0
(12)
To complete the formulation of the problem, a time stepping scheme must be selected and the dependent variables must be expanded in terms of the i~asas functions.
Time stepping The model makes use of three time levels in discretising equations (11) and (12). For simphcity we introduce the notation: 5P =
Pt+At--Pt--at
where the subscnpts refer to the tune at which the function p is evaluated. In the model, the first equation solved is wave equation (11) wluch is discretlsed as 1
DISCRETISATION OF THE G O V E R N I N G EQUATIONS
i = 1, M
(11)
where P is the boundary of the region under study. In lake manner, apphcation of the finite element method to momentum equation (8) yaelds:
] Or t- fHU-- Ay -- HU ~x
= 0
• n~bt ds = 0
F
- m V - Ax ] - 5 [ °(Mvv) + b(HVV) l by [ Ox by + gH a(ItO~ ?-
(10)
~2
zv- Ax] - b [ b ( n v v ) + o(nvv ~y I bx by
+ gH O(H-- h) + fHU + rHV -- Ay] = 0
dA
1
2
At----2 (6H, ¢,)+ 2 At
+
Mv]. n#,dS t F
i = 1,M
(13)
All terms on the fight side of equation (13) are known because they are evaluated at ttrne t or t--At. The surface integral creates no additional problems because it is nonzero only when the node is a boundary node. At the boundary nodes, either the stage or the normal discharge must be specified. If the stage is known, equation (13) is not formulated at that node. If the discharge is specified, then all terms m the surface mtegral will be known. It should be noted that equation (13) is an explicit, centered
Adv Eng Software. 1983, Vol. 5, No. 3
131
approxlmatlon and is rdentrcal to that previously reported by Lynch and Gray 1 For the momentum equations, time differencmg is performed as follows.
The transformation from the x-y domam to the ,$+I domam is accomphshed element-wise by expandmg the x and y co-ordinates m terms of the basis functions such that
$ @iv, &> + cY((T& + fx6v), &)
(19a)
(28 - 1)
Y
=~((V,-Vt--D.t)r~~)-P((rtVt+fXVt),~~)
=
f YJ@jl(& 77) j=1
VW
Then the Jacobian matrix for the transformation is readrly obtamed as
J=
(20)
(14) where OL, 0, y, E, f3 are parameters to be selected by the user, and known quantities appear on the right side of equatron (14) Note that because the wave equation is solved prior to the continuity equation, Ht+At is known for purposes of solvmg the momentum equation. Note that a centred leapfrog scheme is obtamed if 0 = f, cu = 4, /3 = 0, y = 0, E= 1. A two-level explicit scheme which 1s centred m time except for the convective terms results if t9 = 1, a! = /3 = Y=E=& Spatzal dzscretzsatzorz In the model, a Gale&m approach is adopted wherein Lagrangran quadratic rsoparametrrc basis functions are used to drscretise the spatial domain, and as welghtmg functions. A general 9-node quadrilateral element 1s depicted m Fig. la This element is transformed mto a twoumt square m the g-r, coordinate system (Fig lb) for purposes of mtegratron and expansion of the dependent varrable m terms of the basis functrons. The basis functions for node z of the element m Fig. lb may be expressed rn the general form @,(E, r?) = Wk (0 a, (rl)
(15)
This matrix contams all the mformation needed to evaluate derivatrves of dependent vanables with respect to x and y in terms of derivatives with respect l and n Furthermore, the determmant of J 1s the Jacobian which 1s used to transform the mtegrations m equations (13) and (14), which are 111 the x-y plane, to mtegratrons over the 2 x 2 squares m the t-11 plane Dependent vanables are expanded m terms of the basrs functions such that p= F
(21)
PJ(t)@J(x,y)
J=l
where PJ AS the approximation to P at node I. In equation (21) P represents variables such as H, U, V and h but also mcludes products of variables such as HU, HV, rH, HUU, etc. Insertion of expansions such as equation (2 1) mto wave equation (13) yields
5 (< (1 +: At
6HJ@J, @I> = 2((HfJ
--Ht-~t,) $11, @t)
J=1
where ak(t) = f [2 -
3@ [(I +
q (v> = a [2 - 37?;
ttk) - (ik + #I
1 [Cl + m?J) - (q + VI21
z=(tkt2)+3h,+1)
(16)
-Ant2
HJVJVJ
v#J +gHJ c”
(17) (18)
<@‘#k- hk) V+k)
k=l
+fx+& -AJ@J] ,v@I> +
At’(H,Vj@J
and & and q, are the values of .$ and 77 at the node for whtch the basis function is desrred.
z=l,M (22)
Momentum equation (14) becomes Mi g
(
( e + or$At + crAtfxb@J, $r) = (20 - 1)
x ((vr, -
-hat,) $11, 4,) -@t%y,vtJ + fxbJ) $1, @‘r)
(1 - 0) At ((rrJvt -A 5 + fx vr-
At> $1, 41)
-~A~(([~~HJ~~H~+(~-E)H~-A~,IV~J~~~
lb la Fzgure 1. Nzne node Lagrangzun quadratic zsoparametrzc element zn x-y and 5~ coordznates
132
Adv Eng. Software 1983, Vol. 5, No 3
- (hJv$J, $,)>
- A t c” ([email protected]~J’ @I) k=l
+
2lHt+zxt!
+
\2
/ Ht]J
in equations (22) and (23) may be evaluated directly and appear m Tables 1 and 2 respectwely. Inspectton of the terms m the second columns of these tables,provades an indicahon of the amount of work that must be done by the program. It should be noted that the summatmn over e, the element numbers, xs done only over elements which contain node t. The coefficients, a(~ e, ~/e) are equal to 1/9, 4[9 or 16/9 depending on whether the node is a corner, mid-side, or msd-element node respectwely. Evaluatmn of the Jacoblan for element e at node t, je, xs accomplished taking advantage of the fact that a¢~/a~ and aCz/a~/are non-zero at no more than three nodes m an element containing node t. The boundary Jacobian for element e at node t, J~, is determined by transformmg the curved element side m x-y space to a straight side of twoumt length m ~-space.
+,, +,)
= 1,M
(23)
Spattal integratlon The next step m the solutmn process is the evaluatmn of the integrals which appear in equations (22) and (23) Integration is carried out using Sirnpson's rule because this quadrature formula mmtmises the required computational effort over Lagrangian quadratic elements. 3 All mtegratmns are performed m the local ~ coordinate system. For an arbitrary functmn p defined on the m t e r v a l - - I ~ < ~ 1 , -- 1 ~
1
;f --1 --1
p(~,~)d~dr~--- y,
''
~ a(i,DpO,/) (24)
Eq uatto n solution
t=-- 1 1=--1
Combination of the Lagrangaan, quadratic, Isoparametrlc elements with Stmpson's role and the described time stepping scheme has ymlded a set of matrix equatmns for which the matrix is diagonal. Thus the solution of the equations at each node can be made with minimal effort. The solution of the wave equation is accomplished by solving.
where a(t, I) = (-~-- i~) (~ --12). However, because ¢~ Is non-zero at only one of the integration points, 1
1
I- f p(~' ~1)¢' d~ dr/= a(k,l) p(k,1)
(25)
[~ea(~,li~)Je]{O+r'tA--~-f)aH,}=r,
--I --I
where 1 ~
where r~ xs composed of the terms in equatmn (22) which are known at time level t or t - A t . Then Ht+at t IS obtained as
is used,
t= (k +'2) + 3(l + l)
(26)
-- l <.k,l<. l
Note that the nine-point Smapson's formula is thus reduced to a one-point formula. Integratmns of p multxplled by a¢~/Ox or O¢,/by reqmre that the mtegrand be evaluated at five of the rune quadrature points Thus the integrals
Ht+~ti = Ht_gati + ri/{[~e a(~e, 71e) J~ ] [ 1 + rtiAt/2]} (27)
Table 1 Quantltaes appeanng m equation (22)
• (1
Slmpson's Rule mtegratmn of terms in equation (22)
rt~t \
rtlAt\
1=1" + - ' ~ ) 6 H I (¢1, Ck>
e=l
M
Z (H~1-Ht-atl) <¢j, ¢,>
N
(Hti-H,_a, t) y. a(g, he) s e
]=1
e=l
M
N
Htlvtjvt1: 07¢1, V¢i) 1=1 M
M
e=l]=l M
]=1 k = l
M
E ~ ~ Htlvtlvtl: (V¢l)~(VCt)~a(~, r/~) J~ k=l
N M M
Z E g,j(%-hk) <¢/V¢~,V¢,>
~ ~ Htk(Ht/- hi) (V¢l)~ (V¢t)~a(~ek,rl~) J~ e=l]=l N
]
N
1 + ~ - ) a g i ~ a(~ e, r/e) j e
k=l
M
(f x/-/t]v9 -- At/). (¢/, V¢i)
ea e e
A e=lk=l
N M
ZZ Z~v~j~. <~jv~, ~,>
y~ y~ H,~v,~. (V~)ga(~ e, ,7e) je
1 k
e=lk=l
N F
e =1\
\
e2
e
dt
Adv Eng Software. 1983, Vol 5, No 3
133
Table 2. Quantities appearing in momentum equation (23)
Sampson's Rule integration of terms m momentum equahon (23)
M
N
(0 + are At) 8v]<@l,0,>
(0 + ~'rtlAt ) 8v, E a ( ~ , r/e) j e
i=1
e=l
N ~AtfxSv, Z a( ~je, r/e) j~
M
~AtfxSvl<~P], ~> I=I
e=l
M
N
(v,,-v~_a,? y a(~ e, r/e) je
Y~ ( % - vt- ,, t/) <¢,/, 0,>
e=l
I=1
N
M
(rtzVtt+ fxvt) ~ a(~je, r/e) j e
Z (zqvtl + fx vt~)
e=l
M
N
(%v,_a,, + f×v,_a,,) y~ a(~ e, ne) je
(1"tlVt_ A tl + f x vt_ &tl)(q~l, ~b,)
e=l
1=1
M
N
y (r~H, + cH,~ + (1 - e ) ~ , _ a,-h) I=I M
N
M
M
e=lk=l
l=lk=l
~= At+Atl <~I'~P*)
~At+At~ - - Z ~ - a ( ~_, ~e) j f e=l
I
Note that if H is specified at a boundary, this fact may be applied after equation (27) has been applied at all nodes. The momentum equations are solved for U and V from equation (23) and the terms in Table 2 as" 1
[ e~ a(~, 7/e) Je[{(0 +
At) 8U, --fodXt8 Vl} = s~,
J
(28a)
[~ a(~f, 77e) Jr] {(0 + ~t, At) 8 Vz+ faAtSU,}= Sy, L
("/SHk + eHtk + (1 -- e ) H t_ zxtt¢ -- h) (Vq~k)ea (~e, r/e) j e
Z Z v,kv,," (v*k)~a(~f, r/e) je
Y'. 5". vqvo,. <¢kVOj, ¢,>
1 Hi+At
M
2E e=l k=l
-I
(2Sb) Solution of equations (28a) and (28b) for U and V at the new tune level yields
Ut+At, = lit -zXt, + [(0 + ~%At) Sx, + (f~At) Sy,]
(29a) v,+at, = v,-zxtt + [(0 + ~%At)
Sy t - (y~At)
t+At l
BOUNDARY CONDITIONS Because of the explicit nature of the model wherein a solution of equations (27) or (29) at one node is uncoupled from the solution at other nodes, the apphcahon of the boundary condmons is strmghtforward. The model is capable of treating the prescribed conditions as detailed below
Spectftcatton of stage When the stage is spemfied at nodal point t on the boundary, equation (27) for Ht+zxtz is replaced by the known value.
Speclficatton of normal veloctty or normal volumetrw flow The normal velocity or volumetnc flow may be spemfied to be zero, constant or tame varying at node t located on the boundary This speclficatlon affects both the wave equation and the momentum equation solution. If the normal velocity, Vm, as specified at node l, the surface integral term m the wave equation (the last item m Table 1) is chscretlsed in tune to obtain
sx, l \ ~t +rHv .n~idS F (29b)
With these equations, as with the depth, Dmchlet boundary condmons on velocity may be specified after apphcatIon of equataon (29) at all nodes.
134 'Adv Eng Software 1983, Vol 5, No. 3
8 V,,, Htt--Ht-att + H,, -~tt + Tt H, Vn,,) = ~e ( vntz At × ( ~ - ~~,e 2 ~~ J~e
(30)
Note that m the interest of simphcity, the tmae derivative of depth terms has been treated exphclfly and is not centred m tane. It is possible to modify equation (26) such that the surface integral (30) is centred in time. However, the additional computational effort required to incorporate tins feature is significant but does not sxgmficantly unprove the accuracy of the model. When the volumentric flow, H i v t. n = Qnt is specified at node z, the surface integral in the wave equation may be &scretlsed to obtain.
f\-~t r
TQ) " n $ t d S = e ( 2 A t
.'
__,~)
e (31)
Here all approxtmations are centred in tune. Of course specification of the normal velocity or flow must also be incorporated into the momentum equations. However, this boundary condltton may be applied after equation (29) is solved to obtam imtlal estunates of Ut+A t and Vt+~tr The normal direction at node i is obtained directly when the slope of the element boundary is contlnuous at node t or m an average sense when the slope is discontinuous foUowmg the procedure of Gray. 3 Thus:
/I,
sm¢
cos ~JtVx,
where ~ is the angle between the x-axis and the normal chrection and V m and Vii are the normal and tangential velocities respectively. Thus specification of the normal velocity is accomplished in three steps" (1) obtain Vx . from the esttrnates of x. and y.velocity computed froffi equation (29) (denoted U~+at~ and V*+atl) by rotating these velocities; (2) specify Vn t + A t ~ •; (3) rotate Vxt+at i and Vnt+Att oack into x - y co-ordlnatbs. Thus if me normaJ velocity IS specllied at node t: •
rr
,
t + A t ;
-
Ut +A h = Vnt +A tt COSO+ U~+A tz sin2~b -- V*+ A fi sln~ cos~b
(33a) Vt +~ h = Vnt +A h sm~b-- U*+Afi sm~b c o s $ + gT+ A tt cosZ~b
(33b) Specification o f both veloclty or volumetric f l o w components
Within the model, velocity or volumetnc flow components may be specified either in x - y co-ordinates or in normal-tangential co-ordinates. The normal component, either specified darecfly or obtmned from rotation of the x - y components, is used in an evaluation of the surface integral in the wave equations as described by equations (30) or (31) Specified x - y components, or rotated normaltangential components, are used m lieu of equation (29). Of course if the boundary conditions are in terms of flow, they must be divided by H t +ts t to obtain the velocities. I/O FEATURES OF THE MODEL The model has been designed to try to minlmlse the consequences of user error in mput data. Computation may be performed m either the FPS or SI system of units. Checks are made wltinn the model to verify that the expected number of data cards for the various quantities have been
read m, that no dependent variable at any node is specified twice, and that the boundary condition data is read m at the expected number of nodes as well as for all nodes where it is expected. Failure to satisfy these criterion results in spectfic error messages followed by the termmataon of the code prior to the start of computations. Data to be provided IS discussed below. B o t t o m friction spectficatton
Either the Manning or Chezy formula may be used to compute r or a tune-mvariant form for r may be specified. Appropriate values of Manmng's n, Chezy's C, or tnnelnvariant, must be specified at all comer nodes but at side or rind-element nodes only as desired• These parameters will be mterpolated to side or mid-element nodes where no values have been specified. A message is issued indicating winch nodes required interpolation. Latttude
Latitude of the model region must be provided so that the Corlohs acceleration may be properly accounted for. Bathymetry
Depth of water below a reference datum must be provided at all corner nodes, but at side and mid-element nodes only as desired. Values will be mterpolated for side and mid-element nodes where no value has been specified should the user desire. Geometry
Node numbers must be provided for all eight nodes on the boundary of an element in counter-clockwise order. Locations must be specified for all corner nodes• Locations will be interpolated at side nodes where no location has been specified Centre node numbers and locations wall be determined by the computer if not specified by the user. Checks are made to ensure that node numbering is counterclockwise and that side nodes appearmg m two elements lie between the same corner nodes in both elements. lnlttal condtttons
Values of stage and velocity must be specified. At comer nodes where values are not read in, ( H - - h ) , U o r Vwtll be defaulted to zero• Side and mid-element values may be computed by interpolation as desired. Atmospheric forcing
Velocity components of the wind or wind stress must be supplied at each node. Because it may be desirable to change ttus specification as computations proceed, a sub-routine exists which the user may have to alter. Boundary condittons
Because of the variety of forms of tidal forcing that may appear m a model, specification of H, U, V or Vn at a boundary is done through a sub-routine that the user must alter for each application. Both river and sea type flow boundaries may be accommodated• Mtscellan eou s
The model checks to be sure that dimensioned variables have been allocated enough space. It provides the form a new DIMENSION statement should take if desired. Tune weighting parameters, tune step, number of tune steps
Adv. Eng. Software 1983, Vol 5, No. 3
135
7
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
,._._f. ......... .
:/
,
//
t
[. -. _. + . . . . . . . . . ,,~
t
--
~*-
_
_. . . . . . . . . ,,
.,:---~........... ,.
"?.."
,.-.. /-z~'a-.
"~ ..... ....
T , " ~ ; ". •. I
.,,
_/ --"- - - . .-.-.~. - - ~ . ~
. . . .
,"
~
/
,"
- ...... ,S".....
,/
,
.7l==" " ? "
""
::>"
,,t ~,. ";'"
t:,,,, -~----~.L
..... ~"! . . . . . ,/ .... . . . . . i ........
..... l t.. ,
-os
_____~
"~'.
"
," . . . . . Z~:--==7~L-~:~~
.~'--4-~,' ....... ,<......
, .............
t /"
,,
..t... ,'.-3,....!.,~ /
/
d
. . . . . . . . .
".:--.
'-:-.:: " . ~........ " ""-r"
....• " ,'
- - - - : " -~-" ;xl
.... ':.-
,,
,,.
1../;-"....
F~gure 2. Fimte element grid and computed low eo-ttdal and co-range lines for the southern North Sea from QUIET
stratlon of the utility o f the model and not a fine-tuned simulation The North Sea was discretlsed into 45 elements and 211 nodes as in Fig. 2. Bathymetry data was obtmned from sea charts and tidal forcing was interpolated and extrapolated from the data of Leendertse 4 for 12 and 13 September 1958. Atmospheric forcing was neglected and because data for the river inflows from the Thames, Maas and Rhine was not readily avatlable, these inputs were neglected Manmng's n was specified to be constant and equal to 0.035 at all nodes Time weighting parameters selected for equation (23) were 0 = 1.0 and a =/3 = 7 = e = 0.5 A time step of 155 seconds was selected and the model was run for 2880 time steps. Low co-tidal and corange hnes computed during the fifth 24,8 h cycle are presented in Fig. 2. Comparison with the tuned results of Leendertse in Fig 3, indicate that the model does a reasonable j o b o f representing the stage m the North Sea.
CONCLUSION
."~" ~-. i
"--.-.A \
ll-t I
ll\ i%
OOIlr I
I \ i
\-2.00 \
\
\\
,,<,<>,. # ',~.~o
\
~ --
tt - - \ ~I \ ~ 'l 8 0 0 h r ~ , . / ' -
x ..~\
~,,
/ ',,
-
~ %
mmm
\,
/
~ll
,
The features of the QUIET model which make it a simple user-oriented tidal simulation tool have been enumerated. The facts that no matrices need be solved, that the wave equation is used m place o f the prumhve continuity equation, and that ttrne weighting may be applied to the m o m e n t u m equations are all features which contribute to the speed, accuracy and stability of the model. Copies of the documentation are avatlable from the authors. ACKNOWLEDGEMENT
Figure 3. Computed low co-ttdal and co-range lines for the southern North Sea from Leendertse (1967)
This work was supported in part by the National Science Foundation through grant No. CME-7921076.01.
REF ERENCES or simulation period are specified by the user. The model does not include a graphics package for output as requirements vary from system to system. The program contains 1850 statements including comment cards.
EXAMPLE SIMULATION The QUIET model has been applied to the southern portion of the North Sea No adjustment of parameters was attempted so the simulation is merely a first-cut demon-
136
Adv Eng. Software. 1983, Vol 5, No 3
1 Lynch, D. R and Gray, W. G. A wave equation model for fimte element tidal computations, Computers and Fluzds, 1979, 7, 3,207 2 Pritchaxd, D. W. Two-dimensional models, m G. H. Ward and W H. Espey (eds), Estuarme Modehng An Assessment, NTIS No PB206-807, Water Quality Office, U.S.E.P.A., 1971 3 Gray, W G An efficient finite element scheme for twodimensional surface water computation, Proc. Ftrst IntL Conf. on Fmtte Elements m WaterResources, Pentech Press, 1977 4 Leendertse, J. J Aspects of a computational model for longperiod water-wave propagation, Rand Memorandum RM-5294PR, Santa Momca, Cahforma, 1967
The model described herein uses the vertically averaged wave equation and m o m e n t u m equations to simulate two-dimensional tidal flow. The model uses quadratic isoparametric elements in space and a new explicit in time scheme to compute surface elevations and velocities which are essentially noise-free. The model is shown to be efficient in that it requires the solution o f no matrices. A sample application to the southern North Sea demonstrates the utility of the model.
y-direction 3V aV OV --+U--+V--+g-3t ~x ay
a(H-- h) + fU + rV= Ay]H ay
(2b)
or
conservative x-dtrectlon 3(HU)
a(HUU) t - - - - t
ax
~t
D(HUV) -
-
a(H-- h) + g H - -
ay
~x
(3a)
-- f H V + rHU = A x
INTRODUCTION A QUadratic Isoparametnc Exphclt in Tmae (QUIET) finite element model has been developed and documented The model is based on the wave equation formulation of the continuity equaUon descnbed by Lynch and Gray ~ coupled with the freedom to select various time weighting schemes for the momentum equations The QUIET model is noteworthy in that computational noise is suppressed wtule the physical characteristics are not damped. The model also does not require the solution of any matnces. In this report, the computational features of QUIET wtuch account for its speed and versatility are developed Data requirements are also dtscussed, and an example simulation of the southern part of the North Sea is presented.
EQUATION FORMULATION The continuity and momentum equations used when modelling vertacaUy well-mLxed essentially honzontal tidal circulation, are obtained by integrating the point equations through the depth of the flow. A clear derivation of the equatmns is gwen by Prltchard 2 Following ttus work we obtmn Continuity 3H --
+
3t
a(HU)+ a(UV) = ax ay
0
(1)
y-direction a(HUV) + a(HVV)
3(HV) ~t
ax
a ( H - h)
~ gH - -
~y
ay
(3b)
+ f H U + r H V = Ay where
is positive eastward is posttwe northward Y is tnne t u , v are the vertically averaged velocities in the x- and y-directions respectively is the total depth of flow H h xs bathymetry is the Conohs parameter f ts the non-hnear bottom friction 7 g is gravity, and Ax, Ay are the atmospheric wind stresses. X
Most finite element models of estuarlne flow solve equation (1) In conjunction with either the non-conservatwe momentum equations (2) or the conservative forms (3). Lynch and Gray 1 have proposed using a wave equation in place of equation (1) because of some advantageous numerical properties. The wave equation is derived by first differentiating equation (1) with respect to tame to obtain"
Momentum equations
+
non-conservatxve x-direction bU 8U --+U--+ 3t ~x
~U V--+g--
ay
a(H-- h) ax
Next equations (3) are used to replace. f V + rU = A x / H
(2a)
Received August 1982 Written dlscussmn closes September 1983.
130
(4)
Adv Eng Software 1983, Vol 5, No 3
a~u) bt
and
a(HV) Ot
0141-1195/83/030130-07 $2 00 ©1983CMLPubhcatlons
where
and obtain:
b2n at 2
el' (r,b)= IIr.b
b [ b(HUU)+ b(HUI0 t- gH b(H -- h) bx [ bx by bx
-my
+
by
J
(5)
Finally equation (1) is used to replace:
and ~2 is the areal regmn of interest. To facilitate the computations, mtegralaon by parts (Green's theorem) is applied to the third bracketed term m equation (9) and the integral around the boundary of the domain ,s re-expressed using the conservatwe momentum equation. The result of these operations is
<
02H
--r
+
by
J
by
\
aH + f x H v --A], V~b,) --
so that the final form of the wave equation is.
a~H
--+~ at 2
OZ
O rO(HU03 a(HUV)
at
ax
I--
[
Ox
O(H -- h) I - - - + g H - by Ox
y[°7)+,.q
+
aT
-- HV--
ay
h)
z= 1,M
b2H
OH
0t 2
0t
< ( ~ + rv + f x v), ep,)+ ( g V ( , - h ) , ~t)
(6)
The computer model developed solves equations (2) and (6) for H, U and V To make their mampulatton easier in the subsequent portmn of the paper, we can re-write these equations using vector notation:
---I- r----V
. IV. (Hw) + gItV(H--h)+ f x h v - - A ]
-- H r . V r = 0
(7)
av A --+ v.Vv+gV(H--h) +fxv+ rv----= 0 (8) 0t H Dlscretisation of these equalaons is accomphshed by applying the fimte element method in the spatial domam along with a finite chfference marching procedure in the time domain. These aspects of the model will be discussed in some detail so that the reasons for its speed and low computer core requirements are revealed.
+<(v.vv--A),
Applicataon of the finite element method to the wave and momentum equations has been presented previously by Lynch and Gray. 1 Here we will bnefly demonstrate their results and concentrate on the selection of the bass function type, quadrature scheme, and time stepping algonthm. Use of the fimte element method in space requires that the differential equation be orthogonal to each of a set of M basis functions ~bz(x, y):
Ot 2 '
--
i:l,M
#,>= 0
(12)
To complete the formulation of the problem, a time stepping scheme must be selected and the dependent variables must be expanded in terms of the i~asas functions.
Time stepping The model makes use of three time levels in discretising equations (11) and (12). For simphcity we introduce the notation: 5P =
Pt+At--Pt--at
where the subscnpts refer to the tune at which the function p is evaluated. In the model, the first equation solved is wave equation (11) wluch is discretlsed as 1
DISCRETISATION OF THE G O V E R N I N G EQUATIONS
i = 1, M
(11)
where P is the boundary of the region under study. In lake manner, apphcation of the finite element method to momentum equation (8) yaelds:
] Or t- fHU-- Ay -- HU ~x
= 0
• n~bt ds = 0
F
- m V - Ax ] - 5 [ °(Mvv) + b(HVV) l by [ Ox by + gH a(ItO~ ?-
(10)
~2
zv- Ax] - b [ b ( n v v ) + o(nvv ~y I bx by
+ gH O(H-- h) + fHU + rHV -- Ay] = 0
dA
1
2
At----2 (6H, ¢,)+ 2 At
+
Mv]. n#,dS t F
i = 1,M
(13)
All terms on the fight side of equation (13) are known because they are evaluated at ttrne t or t--At. The surface integral creates no additional problems because it is nonzero only when the node is a boundary node. At the boundary nodes, either the stage or the normal discharge must be specified. If the stage is known, equation (13) is not formulated at that node. If the discharge is specified, then all terms m the surface mtegral will be known. It should be noted that equation (13) is an explicit, centered
Adv Eng Software. 1983, Vol. 5, No. 3
131
approxlmatlon and is rdentrcal to that previously reported by Lynch and Gray 1 For the momentum equations, time differencmg is performed as follows.
The transformation from the x-y domam to the ,$+I domam is accomphshed element-wise by expandmg the x and y co-ordinates m terms of the basis functions such that
$ @iv, &> + cY((T& + fx6v), &)
(19a)
(28 - 1)
Y
=~((V,-Vt--D.t)r~~)-P((rtVt+fXVt),~~)
=
f YJ@jl(& 77) j=1
VW
Then the Jacobian matrix for the transformation is readrly obtamed as
J=
(20)
(14) where OL, 0, y, E, f3 are parameters to be selected by the user, and known quantities appear on the right side of equatron (14) Note that because the wave equation is solved prior to the continuity equation, Ht+At is known for purposes of solvmg the momentum equation. Note that a centred leapfrog scheme is obtamed if 0 = f, cu = 4, /3 = 0, y = 0, E= 1. A two-level explicit scheme which 1s centred m time except for the convective terms results if t9 = 1, a! = /3 = Y=E=& Spatzal dzscretzsatzorz In the model, a Gale&m approach is adopted wherein Lagrangran quadratic rsoparametrrc basis functions are used to drscretise the spatial domain, and as welghtmg functions. A general 9-node quadrilateral element 1s depicted m Fig. la This element is transformed mto a twoumt square m the g-r, coordinate system (Fig lb) for purposes of mtegratron and expansion of the dependent varrable m terms of the basis functrons. The basis functions for node z of the element m Fig. lb may be expressed rn the general form @,(E, r?) = Wk (0 a, (rl)
(15)
This matrix contams all the mformation needed to evaluate derivatrves of dependent vanables with respect to x and y in terms of derivatives with respect l and n Furthermore, the determmant of J 1s the Jacobian which 1s used to transform the mtegrations m equations (13) and (14), which are 111 the x-y plane, to mtegratrons over the 2 x 2 squares m the t-11 plane Dependent vanables are expanded m terms of the basrs functions such that p= F
(21)
PJ(t)@J(x,y)
J=l
where PJ AS the approximation to P at node I. In equation (21) P represents variables such as H, U, V and h but also mcludes products of variables such as HU, HV, rH, HUU, etc. Insertion of expansions such as equation (2 1) mto wave equation (13) yields
5 (< (1 +: At
6HJ@J, @I> = 2((HfJ
--Ht-~t,) $11, @t)
J=1
where ak(t) = f [2 -
3@ [(I +
q (v> = a [2 - 37?;
ttk) - (ik + #I
1 [Cl + m?J) - (q + VI21
z=(tkt2)+3h,+1)
(16)
-Ant2
HJVJVJ
v#J +gHJ c”
(17) (18)
<@‘#k- hk) V+k)
k=l
+fx+& -AJ@J] ,v@I> +
At’(H,Vj@J
and & and q, are the values of .$ and 77 at the node for whtch the basis function is desrred.
z=l,M (22)
Momentum equation (14) becomes Mi g
(
( e + or$At + crAtfxb@J, $r) = (20 - 1)
x ((vr, -
-hat,) $11, 4,) -@t%y,vtJ + fxbJ) $1, @‘r)
(1 - 0) At ((rrJvt -A 5 + fx vr-
At> $1, 41)
-~A~(([~~HJ~~H~+(~-E)H~-A~,IV~J~~~
lb la Fzgure 1. Nzne node Lagrangzun quadratic zsoparametrzc element zn x-y and 5~ coordznates
132
Adv Eng. Software 1983, Vol. 5, No 3
- (hJv$J, $,)>
- A t c” ([email protected]~J’ @I) k=l
+
2lHt+zxt!
+
\2
/ Ht]J
in equations (22) and (23) may be evaluated directly and appear m Tables 1 and 2 respectwely. Inspectton of the terms m the second columns of these tables,provades an indicahon of the amount of work that must be done by the program. It should be noted that the summatmn over e, the element numbers, xs done only over elements which contain node t. The coefficients, a(~ e, ~/e) are equal to 1/9, 4[9 or 16/9 depending on whether the node is a corner, mid-side, or msd-element node respectwely. Evaluatmn of the Jacoblan for element e at node t, je, xs accomplished taking advantage of the fact that a¢~/a~ and aCz/a~/are non-zero at no more than three nodes m an element containing node t. The boundary Jacobian for element e at node t, J~, is determined by transformmg the curved element side m x-y space to a straight side of twoumt length m ~-space.
+,, +,)
= 1,M
(23)
Spattal integratlon The next step m the solutmn process is the evaluatmn of the integrals which appear in equations (22) and (23) Integration is carried out using Sirnpson's rule because this quadrature formula mmtmises the required computational effort over Lagrangian quadratic elements. 3 All mtegratmns are performed m the local ~ coordinate system. For an arbitrary functmn p defined on the m t e r v a l - - I ~ < ~ 1 , -- 1 ~
1
;f --1 --1
p(~,~)d~dr~--- y,
''
~ a(i,DpO,/) (24)
Eq uatto n solution
t=-- 1 1=--1
Combination of the Lagrangaan, quadratic, Isoparametrlc elements with Stmpson's role and the described time stepping scheme has ymlded a set of matrix equatmns for which the matrix is diagonal. Thus the solution of the equations at each node can be made with minimal effort. The solution of the wave equation is accomplished by solving.
where a(t, I) = (-~-- i~) (~ --12). However, because ¢~ Is non-zero at only one of the integration points, 1
1
I- f p(~' ~1)¢' d~ dr/= a(k,l) p(k,1)
(25)
[~ea(~,li~)Je]{O+r'tA--~-f)aH,}=r,
--I --I
where 1 ~
where r~ xs composed of the terms in equatmn (22) which are known at time level t or t - A t . Then Ht+at t IS obtained as
is used,
t= (k +'2) + 3(l + l)
(26)
-- l <.k,l<. l
Note that the nine-point Smapson's formula is thus reduced to a one-point formula. Integratmns of p multxplled by a¢~/Ox or O¢,/by reqmre that the mtegrand be evaluated at five of the rune quadrature points Thus the integrals
Ht+~ti = Ht_gati + ri/{[~e a(~e, 71e) J~ ] [ 1 + rtiAt/2]} (27)
Table 1 Quantltaes appeanng m equation (22)
• (1
Slmpson's Rule mtegratmn of terms in equation (22)
rt~t \
rtlAt\
1=1" + - ' ~ ) 6 H I (¢1, Ck>
e=l
M
Z (H~1-Ht-atl) <¢j, ¢,>
N
(Hti-H,_a, t) y. a(g, he) s e
]=1
e=l
M
N
Htlvtjvt1: 07¢1, V¢i) 1=1 M
M
e=l]=l M
]=1 k = l
M
E ~ ~ Htlvtlvtl: (V¢l)~(VCt)~a(~, r/~) J~ k=l
N M M
Z E g,j(%-hk) <¢/V¢~,V¢,>
~ ~ Htk(Ht/- hi) (V¢l)~ (V¢t)~a(~ek,rl~) J~ e=l]=l N
]
N
1 + ~ - ) a g i ~ a(~ e, r/e) j e
k=l
M
(f x/-/t]v9 -- At/). (¢/, V¢i)
ea e e
A e=lk=l
N M
ZZ Z~v~j~. <~jv~, ~,>
y~ y~ H,~v,~. (V~)ga(~ e, ,7e) je
1 k
e=lk=l
N F
e =1\
\
e2
e
dt
Adv Eng Software. 1983, Vol 5, No 3
133
Table 2. Quantities appearing in momentum equation (23)
Sampson's Rule integration of terms m momentum equahon (23)
M
N
(0 + are At) 8v]<@l,0,>
(0 + ~'rtlAt ) 8v, E a ( ~ , r/e) j e
i=1
e=l
N ~AtfxSv, Z a( ~je, r/e) j~
M
~AtfxSvl<~P], ~> I=I
e=l
M
N
(v,,-v~_a,? y a(~ e, r/e) je
Y~ ( % - vt- ,, t/) <¢,/, 0,>
e=l
I=1
N
M
(rtzVtt+ fxvt) ~ a(~je, r/e) j e
Z (zqvtl + fx vt~)
e=l
M
N
(%v,_a,, + f×v,_a,,) y~ a(~ e, ne) je
(1"tlVt_ A tl + f x vt_ &tl)(q~l, ~b,)
e=l
1=1
M
N
y (r~H, + cH,~ + (1 - e ) ~ , _ a,-h)
N
M
M
e=lk=l
l=lk=l
~= At+Atl <~I'~P*)
~At+At~ - - Z ~ - a ( ~_, ~e) j f e=l
I
Note that if H is specified at a boundary, this fact may be applied after equation (27) has been applied at all nodes. The momentum equations are solved for U and V from equation (23) and the terms in Table 2 as" 1
[ e~ a(~, 7/e) Je[{(0 +
At) 8U, --fodXt8 Vl} = s~,
J
(28a)
[~ a(~f, 77e) Jr] {(0 + ~t, At) 8 Vz+ faAtSU,}= Sy, L
("/SHk + eHtk + (1 -- e ) H t_ zxtt¢ -- h) (Vq~k)ea (~e, r/e) j e
Z Z v,kv,," (v*k)~a(~f, r/e) je
Y'. 5". vqvo,. <¢kVOj, ¢,>
1 Hi+At
M
2E e=l k=l
-I
(2Sb) Solution of equations (28a) and (28b) for U and V at the new tune level yields
Ut+At, = lit -zXt, + [(0 + ~%At) Sx, + (f~At) Sy,]
(29a) v,+at, = v,-zxtt + [(0 + ~%At)
Sy t - (y~At)
t+At l
BOUNDARY CONDITIONS Because of the explicit nature of the model wherein a solution of equations (27) or (29) at one node is uncoupled from the solution at other nodes, the apphcahon of the boundary condmons is strmghtforward. The model is capable of treating the prescribed conditions as detailed below
Spectftcatton of stage When the stage is spemfied at nodal point t on the boundary, equation (27) for Ht+zxtz is replaced by the known value.
Speclficatton of normal veloctty or normal volumetrw flow The normal velocity or volumetnc flow may be spemfied to be zero, constant or tame varying at node t located on the boundary This speclficatlon affects both the wave equation and the momentum equation solution. If the normal velocity, Vm, as specified at node l, the surface integral term m the wave equation (the last item m Table 1) is chscretlsed in tune to obtain
sx, l \ ~t +rHv .n~idS F (29b)
With these equations, as with the depth, Dmchlet boundary condmons on velocity may be specified after apphcatIon of equataon (29) at all nodes.
134 'Adv Eng Software 1983, Vol 5, No. 3
8 V,,, Htt--Ht-att + H,, -~tt + Tt H, Vn,,) = ~e ( vntz At × ( ~ - ~~,e 2 ~~ J~e
(30)
Note that m the interest of simphcity, the tmae derivative of depth terms has been treated exphclfly and is not centred m tane. It is possible to modify equation (26) such that the surface integral (30) is centred in time. However, the additional computational effort required to incorporate tins feature is significant but does not sxgmficantly unprove the accuracy of the model. When the volumentric flow, H i v t. n = Qnt is specified at node z, the surface integral in the wave equation may be &scretlsed to obtain.
f\-~t r
TQ) " n $ t d S = e ( 2 A t
.'
__,~)
e (31)
Here all approxtmations are centred in tune. Of course specification of the normal velocity or flow must also be incorporated into the momentum equations. However, this boundary condltton may be applied after equation (29) is solved to obtam imtlal estunates of Ut+A t and Vt+~tr The normal direction at node i is obtained directly when the slope of the element boundary is contlnuous at node t or m an average sense when the slope is discontinuous foUowmg the procedure of Gray. 3 Thus:
/I,
sm¢
cos ~JtVx,
where ~ is the angle between the x-axis and the normal chrection and V m and Vii are the normal and tangential velocities respectively. Thus specification of the normal velocity is accomplished in three steps" (1) obtain Vx . from the esttrnates of x. and y.velocity computed froffi equation (29) (denoted U~+at~ and V*+atl) by rotating these velocities; (2) specify Vn t + A t ~ •; (3) rotate Vxt+at i and Vnt+Att oack into x - y co-ordlnatbs. Thus if me normaJ velocity IS specllied at node t: •
rr
,
t + A t ;
-
Ut +A h = Vnt +A tt COSO+ U~+A tz sin2~b -- V*+ A fi sln~ cos~b
(33a) Vt +~ h = Vnt +A h sm~b-- U*+Afi sm~b c o s $ + gT+ A tt cosZ~b
(33b) Specification o f both veloclty or volumetric f l o w components
Within the model, velocity or volumetnc flow components may be specified either in x - y co-ordinates or in normal-tangential co-ordinates. The normal component, either specified darecfly or obtmned from rotation of the x - y components, is used in an evaluation of the surface integral in the wave equations as described by equations (30) or (31) Specified x - y components, or rotated normaltangential components, are used m lieu of equation (29). Of course if the boundary conditions are in terms of flow, they must be divided by H t +ts t to obtain the velocities. I/O FEATURES OF THE MODEL The model has been designed to try to minlmlse the consequences of user error in mput data. Computation may be performed m either the FPS or SI system of units. Checks are made wltinn the model to verify that the expected number of data cards for the various quantities have been
read m, that no dependent variable at any node is specified twice, and that the boundary condition data is read m at the expected number of nodes as well as for all nodes where it is expected. Failure to satisfy these criterion results in spectfic error messages followed by the termmataon of the code prior to the start of computations. Data to be provided IS discussed below. B o t t o m friction spectficatton
Either the Manning or Chezy formula may be used to compute r or a tune-mvariant form for r may be specified. Appropriate values of Manmng's n, Chezy's C, or tnnelnvariant, must be specified at all comer nodes but at side or rind-element nodes only as desired• These parameters will be mterpolated to side or mid-element nodes where no values have been specified. A message is issued indicating winch nodes required interpolation. Latttude
Latitude of the model region must be provided so that the Corlohs acceleration may be properly accounted for. Bathymetry
Depth of water below a reference datum must be provided at all corner nodes, but at side and mid-element nodes only as desired. Values will be mterpolated for side and mid-element nodes where no value has been specified should the user desire. Geometry
Node numbers must be provided for all eight nodes on the boundary of an element in counter-clockwise order. Locations must be specified for all corner nodes• Locations will be interpolated at side nodes where no location has been specified Centre node numbers and locations wall be determined by the computer if not specified by the user. Checks are made to ensure that node numbering is counterclockwise and that side nodes appearmg m two elements lie between the same corner nodes in both elements. lnlttal condtttons
Values of stage and velocity must be specified. At comer nodes where values are not read in, ( H - - h ) , U o r Vwtll be defaulted to zero• Side and mid-element values may be computed by interpolation as desired. Atmospheric forcing
Velocity components of the wind or wind stress must be supplied at each node. Because it may be desirable to change ttus specification as computations proceed, a sub-routine exists which the user may have to alter. Boundary condittons
Because of the variety of forms of tidal forcing that may appear m a model, specification of H, U, V or Vn at a boundary is done through a sub-routine that the user must alter for each application. Both river and sea type flow boundaries may be accommodated• Mtscellan eou s
The model checks to be sure that dimensioned variables have been allocated enough space. It provides the form a new DIMENSION statement should take if desired. Tune weighting parameters, tune step, number of tune steps
Adv. Eng. Software 1983, Vol 5, No. 3
135
7
. . . . . . . . . . . . . . . . . . . . . . . . . . . .
,._._f. ......... .
:/
,
//
t
[. -. _. + . . . . . . . . . ,,~
t
--
~*-
_
_. . . . . . . . . ,,
.,:---~........... ,.
"?.."
,.-.. /-z~'a-.
"~ ..... ....
T , " ~ ; ". •. I
.,,
_/ --"- - - . .-.-.~. - - ~ . ~
. . . .
,"
~
/
,"
- ...... ,S".....
,/
,
.7l==" " ? "
""
::>"
,,t ~,. ";'"
t:,,,, -~----~.L
..... ~"! . . . . . ,/ .... . . . . . i ........
..... l t.. ,
-os
_____~
"~'.
"
," . . . . . Z~:--==7~L-~:~~
.~'--4-~,' ....... ,<......
, .............
t /"
,,
..t... ,'.-3,....!.,~ /
/
d
. . . . . . . . .
".:--.
'-:-.:: " . ~........ " ""-r"
....• " ,'
- - - - : " -~-" ;xl
.... ':.-
,,
,,.
1../;-"....
F~gure 2. Fimte element grid and computed low eo-ttdal and co-range lines for the southern North Sea from QUIET
stratlon of the utility o f the model and not a fine-tuned simulation The North Sea was discretlsed into 45 elements and 211 nodes as in Fig. 2. Bathymetry data was obtmned from sea charts and tidal forcing was interpolated and extrapolated from the data of Leendertse 4 for 12 and 13 September 1958. Atmospheric forcing was neglected and because data for the river inflows from the Thames, Maas and Rhine was not readily avatlable, these inputs were neglected Manmng's n was specified to be constant and equal to 0.035 at all nodes Time weighting parameters selected for equation (23) were 0 = 1.0 and a =/3 = 7 = e = 0.5 A time step of 155 seconds was selected and the model was run for 2880 time steps. Low co-tidal and corange hnes computed during the fifth 24,8 h cycle are presented in Fig. 2. Comparison with the tuned results of Leendertse in Fig 3, indicate that the model does a reasonable j o b o f representing the stage m the North Sea.
CONCLUSION
."~" ~-. i
"--.-.A \
ll-t I
ll\ i%
OOIlr I
I \ i
\-2.00 \
\
\\
,,<,<>,. # ',~.~o
\
~ --
tt - - \ ~I \ ~ 'l 8 0 0 h r ~ , . / ' -
x ..~\
~,,
/ ',,
-
~ %
mmm
\,
/
~ll
,
The features of the QUIET model which make it a simple user-oriented tidal simulation tool have been enumerated. The facts that no matrices need be solved, that the wave equation is used m place o f the prumhve continuity equation, and that ttrne weighting may be applied to the m o m e n t u m equations are all features which contribute to the speed, accuracy and stability of the model. Copies of the documentation are avatlable from the authors. ACKNOWLEDGEMENT
Figure 3. Computed low co-ttdal and co-range lines for the southern North Sea from Leendertse (1967)
This work was supported in part by the National Science Foundation through grant No. CME-7921076.01.
REF ERENCES or simulation period are specified by the user. The model does not include a graphics package for output as requirements vary from system to system. The program contains 1850 statements including comment cards.
EXAMPLE SIMULATION The QUIET model has been applied to the southern portion of the North Sea No adjustment of parameters was attempted so the simulation is merely a first-cut demon-
136
Adv Eng. Software. 1983, Vol 5, No 3
1 Lynch, D. R and Gray, W. G. A wave equation model for fimte element tidal computations, Computers and Fluzds, 1979, 7, 3,207 2 Pritchaxd, D. W. Two-dimensional models, m G. H. Ward and W H. Espey (eds), Estuarme Modehng An Assessment, NTIS No PB206-807, Water Quality Office, U.S.E.P.A., 1971 3 Gray, W G An efficient finite element scheme for twodimensional surface water computation, Proc. Ftrst IntL Conf. on Fmtte Elements m WaterResources, Pentech Press, 1977 4 Leendertse, J. J Aspects of a computational model for longperiod water-wave propagation, Rand Memorandum RM-5294PR, Santa Momca, Cahforma, 1967