A coastal ocean model with subgrid approximation

Accepted Manuscript

A Coastal Ocean Model with Subgrid Approximation Roy A. Walters PII: DOI: Reference:

S1463-5003(16)30010-5 10.1016/j.ocemod.2016.04.005 OCEMOD 1093

To appear in:

Ocean Modelling

Received date: Revised date: Accepted date:

3 January 2016 11 April 2016 16 April 2016

Please cite this article as: Roy A. Walters, A Coastal Ocean Model with Subgrid Approximation, Ocean Modelling (2016), doi: 10.1016/j.ocemod.2016.04.005

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT

Highlights • A semi-implicit coastal ocean model balances stability, accuracy, and efficiency.

• The methods are applied to tidal power generation.

CR IP T

• Volume averaging is used to approximate subgrid objects of various kinds.

• Thrust and power coefficients convert from free stream to volume average velocity.

AC

CE

PT

ED

M

AN US

• Three conversion methods are an iterative approach and two analytical mappings.

1

ACCEPTED MANUSCRIPT

1

A Coastal Ocean Model with Subgrid Approximation

2

Roy A. Waltersa,∗ a Ocean-River

4

Hydrodynamics, Victoria, BC, Canada

CR IP T

3

Abstract

A wide variety of coastal ocean models exist, each having attributes that reflect specific application areas. The model presented here is based on finite element methods with

AN US

unstructured grids containing triangular and quadrilateral elements. The model optimizes robustness, accuracy, and efficiency by using semi-implicit methods in time in

order to remove the most restrictive stability constraints, by using a semi-Lagrangian advection approximation to remove Courant number constraints, and by solving a wave equation at the discrete level for enhanced efficiency. An added feature is the approximation of the effects of subgrid objects. Here, the Reynolds-averaged Navier-Stokes

M

equations and the incompressibility constraint are volume averaged over one or more computational cells. This procedure gives rise to new terms which must be approx-

ED

imated as a closure problem. A study of tidal power generation is presented as an example of this method. A problem that arises is specifying appropriate thrust and power coefficients for the volume averaged velocity when they are usually referenced

PT

to free stream velocity. A new contribution here is the evaluation of three approaches to this problem: an iteration procedure and two mapping formulations. All three sets

CE

of results for thrust (form drag) and power are in reasonable agreement. Keywords: coastal ocean, numerical model, unstructured grid, subgrid, volume

6

averaging, tidal power generation, Canada, Nova Scotia, Minas Passage.

AC

5

7

1. Introduction

8

A wide range of physical phenomena can be found in the area between the continen-

9

tal slope and the landward extent of the coastal ocean which may include complicated

10

interconnecting waterways. Among other things, stratification, rotation, strong advec∗ Corresponding author Preprint submitted Elsevier Email address:[email protected] (Roy A. Walters)

April 20, 2016

ACCEPTED MANUSCRIPT

11

tion, wetting and drying, wave dispersion, scalar transport, and other types of dynamics

12

may be globally or locally important. Hence, it is not surprising that a large number of

13

coastal ocean models exist and they tend to specialize in a specific set of dynamics. The basic set of equations to solve are the well known 3-dimensional version of the

15

shallow water equations (SWE), perhaps with added features for nonhydrostatic flows,

16

dispersive waves, and scalar transport. In the early efforts to approximate these equa-

17

tions, a serious problem was encountered with spurious modes generated by coupling

18

between the pressure gradient in the momentum equations and flux term in the free

19

surface (continuity) equation. For finite difference approximations, this problem was

20

resolved with a staggered grid (C-grid) with structured grids [15] and a more recent

21

staggered grid approximation for irregular orthogonal grids [5]. For finite element ap-

22

proximations with unstructured grids, the problem led to considerably more research

23

before it was adequately resolved.

AN US

CR IP T

14

For the finite element approximations, two main strategies were used: mathemati-

25

cally manipulate the SWE to form a wave equation to replace the free surface equation

26

and use simple linear or p-order elements, or find element approximations that do not

27

support spurious modes [24]. The wave equation approach that is spectral in time pro-

28

vides an accurate and extremely efficient method to simulate tidal flows, particularly

29

in large domains. The main problem that remains is to deal with the compound tides

30

and overtides generated by the nonlinear terms and their influence on the main con-

31

stituents. The wave equation approach that uses time stepping can also be accurate and

32

efficient. The main problem encountered is approximating the advection term because

33

it is now a higher-order derivative, and treating wetting and drying with a node-based

CE

PT

ED

M

24

34

approximation. The examination of various element approximations that do not support spurious

36

modes has evolved over a number of years (see for instance [14]). Initially, mixed

AC

35

37

approximations were used such as linear approximation for sea level and a quadratic

38

approximation for velocity. This strategy required that large, global element matrices

39

be solved with the result that these methods are very inefficient and have run times

40

perhaps a thousand times longer than more recent models.

41

If we fast-forward to the present, a very useful element is the Raviart-Thomas ele3

ACCEPTED MANUSCRIPT

ment of the lowest order (RT0 ) [20]. This element uses a piecewise constant approx-

43

imation for sea level and a linear vector approximation for velocity on triangular and

44

quadrilateral elements. Many studies have shown that this element is accurate, robust,

45

and can lead to efficiencies comparable to wave equation approaches [26, 25, 30, 24].

46

Moreover, this approach has distinct advantages when approximating the advection

47

and Coriolis terms, and greatly simplifies the treatment of wetting and drying. This

48

element approximation is used in the model described in this report, RiCOM – River

49

and Coastal Ocean Model.

CR IP T

42

The next issue to consider is approximating the effects of subgrid objects. The ba-

51

sic problem is that the real world is continuous but numerical models are discretized.

52

Hence there is always some spatial and temporal scale where subgrid features are not

53

resolved. Some examples of subgrid features include submerged and emergent vegeta-

54

tion, arrays of pilings, aquaculture drop lines, and freestanding rotors for tidal power

55

generation. A common characteristic of these features is that they all can be simulated

56

using highly refined CFD models, but they are impractical to fully resolve with a larger-

57

scale coastal ocean model. Hence, some form of subgrid approximation is required. In

58

the end, the resolution is determined by a balance between computer resources and the

59

details of the problem being considered. Nonetheless, there is always some aspect of

60

the problem that cannot be resolved adequately so the question that arises is what can

61

be done about the effects.

ED

M

AN US

50

In the time domain, the problem has been largely resolved through Reynolds aver-

63

aging. The basic Navier-Stokes equations are time-averaged to arrive at the Reynolds

64

equations where the nonlinear terms give rise to averages of fluctuating components

CE

PT

62

known as Reynolds stresses. This approach then leads to a closure problem for the new

66

terms. Considerable research has been devoted to developing closures based on mixing

67

length theory, turbulence submodels, and other higher order closures.

AC

65

68

The averaging concept can be extended to the spatial domain through the use of

69

volume averages [16]. The approach leads to new terms that involve averages of spatial

70

fluctuations and requires closures similar to the time-averaging case. Although much

71

research has gone into this area, much remains to be learned.

72

The new spatially-averaged terms generally lead to some type of form drag and a 4

ACCEPTED MANUSCRIPT

73

dispersive term. The methods that can be used to approximate these terms range from

74

empirical formulations, to combined analytical-empirical formulations, and to coupled

75

large-scale and small-scale models. An example of empirical formulations is the approximation of form drag from

77

emergent and submerged vegetation using relations derived from laboratory experi-

78

ments [31, 29]. Here the drag coefficient can be formulated in terms of vegetation

79

diameter, length, spacing, and density. Reasonably accurate results were obtained in

80

comparing model results with experiments. However, the real world is characterized

81

by high spatial variability in both topography and the distribution of subgrid objects so

82

idealized experimental results may just provide a first guess.

AN US

CR IP T

76

Aquaculture is an example where a combined empirical and analytical approach

84

has been successful. An interesting aspect of this application is that the canopy is

85

suspended from the surface rather than emergent from the bottom of the water column.

86

Taking into account the different physical layers in the vertical, the problem can be cast

87

in terms of a two-dimensional model [19] or in terms of a three-dimensional model

88

[17].

M

83

Coupling large and small scale models seems attractive but presents another set of

90

difficulties. Small-scale computational fluid dynamics (CFD) models and larger-scale

91

coastal ocean (CO) models operate at much different time and space scales. Because

92

of this, it is very difficult to couple them directly due to issues with interface boundary

93

conditions, sub-cycling the CFD model, computational overhead of the CFD model,

94

and other problems. In [21] we develop an alternative approach where volume aver-

95

aging is used to pass information from the CFD model to the coastal ocean model and

CE

PT

ED

89

96

bypass most of these problems. The objective of this paper is to describe the methodology used in developing this

98

coastal ocean model and the approximation for the effects of subgrid objects. The

AC

97

99

approach adopted here is to use a volume average over one or more finite elements

100

(computational cells) and formulate the resultant closures in terms of form drag and

101

spatial dispersion.

102

As an example of this approach, the estimation of power output from free-standing

103

tidal turbines is considered. Currently, turbines for this purpose are being developed 5

ACCEPTED MANUSCRIPT

and tested in several countries and can take many forms, including but not limited to

105

vertical shaft rotors, and ducted and non-ducted impellers. An important question that

106

has arisen, is how to estimate the power that can be extracted for a given site. The

107

answer depends on many factors–the geometry of the site, the tidal range, the number

108

and size of the turbines, and the efficiency of the turbines as a function of flow. The

109

field site that has been selected is Minas Passage in the Bay of Fundy which provides an

110

interesting practical perspective for this problem and represents a single tidal boundary

111

problem that has been examined in [13, 30].

CR IP T

104

In section 2, a description of the model is presented along with the treatment of

113

tidal turbines as form drag. Section 3 contains a description the field site and the model

114

setup. Following that is a presentation of the results and discussion, and conclusions.

115

2. Approach

116

2.1. Basic equations

AN US

112

The starting point in the derivation is the Navier-Stokes equations and continuity

118

equation (incompressibility condition). In the derivation, these equations are first av-

119

eraged in time using traditional Reynolds averaging, then volume-averaged in space

120

using double averaging methods (DAM) described in [16]. The continuity equation is

ED

M

117

∂ρ + ∇3 · (ρu) = 0 ∂t

PT

and the Navier-Stokes equations in a rotating frame of reference are

CE

121

(1)

Du 1 + Ω × u + ∇3 p − ν∇23 u − g = 0 Dt ρ

(2)

where D/Dt is a material derivative, u is velocity, t is time, Ω is the rotation vector

123

which is twice the Earth’s angular velocity, ρ is density which is assumed constant, p is

AC

122

124 125

pressure, ν is kinematic viscosity, g is the gravity vector, and ∇3 is the 3-dimensional gradient operator.

126

¯ + u0 , and Then decompose velocity using a Reynolds decomposition , u = u

127

¯ = h¯ a spatial decomposition, u ui + u00 where the overbar denotes a time average,

128

brackets denote a spatial average, u0 is the velocity deviation in time, and u00 is the 6

ACCEPTED MANUSCRIPT

129

deviation of the time-averaged velocity in space. In addition, let φ = Vf /V be the

130

fluid fraction, where Vf is the volume of fluid in the averaging volume V .

∂φ + ∇3 · (φh¯ ui) = 0 ∂t 132

and the momentum equation becomes

CR IP T

Averaging over time and over space (see [16]), the continuity equation becomes

131

AN US

Dh¯ ui 1 1 ¯ i) − g + Ω × h¯ ui + ∇3 (φh¯ pi) − ∇3 · (φνh∇3 u Dt ρφ φ 1 1 + ∇3 · (φhu0 u0 i) + ∇3 · (φhu00 u00 i) φ φ I I 1 1 ¯) · n + p¯n ˆ dS − (ν∇3 u ˆ dS = 0 ρVf Vf

(3)

(4)

Some new terms have been created by averaging the nonlinear terms and are similar

134

in concept to the derivation of Reynolds stresses through time averaging. In eq. 4,

135

the sixth term is the spatially averaged Reynolds stress, the seventh term is a spatial

136

dispersion term, and the last two terms are form drag and skin friction. Because form

137

drag generally dominates skin friction, the last two terms can be written as a drag force

ED

M

133

fF D = (1/2)ρAf Cd |ur |ur

(5)

where Af is frontal area of an object measured in the plane that is normal to the di-

139

rection of flow, Cd is the non-dimensional drag coefficient, and ur is the reference

140

velocity. The properties of the spatial dispersion term are not known in any detail but

141

appear to control the wake mixing behind the subgrid object. These issues are discussed

142

in the results section.

CE

PT

138

Next, the equations are cast in the form of the 3-dimensional shallow water equa-

AC

143 144

145

tions by using the hydrostatic assumption. Rewriting the continuity equation as ∂Vf + ∂t

I

h¯ uidS = 0

and using the hydrostatic approximation in the momentum equation

7

(6)

ACCEPTED MANUSCRIPT

+ − +

146

1 g∇(φh¯ η i) φ ¯ 1 ∂ ∂u 1 (φhAv i) − ∇ · (φhAh ∇¯ ui) φ ∂z ∂z φ 1 ∇3 · (φhu00 u00 i) − FD = 0 φ

ˆ × h¯ fz ui +

CR IP T

Dh¯ ui Dt

(7)

where ∇ is the horizontal gradient operator, Reynolds stresses are written using the

147

eddy viscosity concept with coefficients Av and Ah , and FD contains the form drag

148

and friction terms in the last line of eq. 4.

In the computational domain where there are no subgrid objects, φ = 1 and the

150

equations reduce to the standard SWE with additional terms from the volume average.

151

In addition, for the types of coastal ocean problems considered here, φ ≈ 1 and is

AN US

149

152

constant so the equations can again be cast into a more traditional form where the

153

continuity equation is

∂w =0 ∂z

(8)

and the momentum equation expressed in nonconservative form is

ED

154

M

∇·u+

∂ Du ∂u ˆ × u + g∇η − + fz (Av ) − ∇ · (Ah ∇u) + F = 0 Dt ∂z ∂z

(9)

where all variables are time and space averaged so the overbar and brackets can be

156

dropped, the coordinate directions (x, y, z) are aligned in the east, north, and ver-

157

tical directions; u(x, y, z, t) is now the horizontal velocity with components (u, v);

158

ˆ is the upward w(x, y, z, t) is the vertical velocity; f (x, y) is the Coriolis parameter; z

159

unit vector; η(x, y, t) is the distance from the reference elevation to the free surface;

160

g is the gravitational acceleration; Av (x, y, z, t) and Ah (x, y) are the kinematic verti-

AC

CE

PT

155

161

cal and horizontal viscosity, respectively; ∇ is the horizontal gradient operator (∂/∂x,

162

∂/∂y); and F contains body forces that include form drag, skin friction, and spatial

163

dispersion terms derived above.

164

The Reynolds stress terms have been written in terms of an eddy viscosity coeffi-

165

cient and are incorporated into Av and Ah . These terms must be approximated through 8

ACCEPTED MANUSCRIPT

166

a turbulence closure submodel or by empirical methods. Note that because of volume

167

averaging, the stress terms should also be volume averaged. This is an active research

168

subject and more information may be found in [7]. For the most part, then, the time and volume averaged equations reduce to the

170

standard SWE with additional terms that account for the effects of the subgrid object

171

on the flow. Locally where the fluid volume fraction, φ, varies significantly, the more

172

exact averaged equations can be used.

The free surface equation is derived by vertically-integrating the continuity equa-

173 174

CR IP T

169

tion and using the kinematic free surface and bottom boundary conditions [18].

AN US

Z η ∂η +∇·( udz) = 0 ∂t h

(10)

Boundary conditions for Eqs. 8 to 10 generally fall into two categories: conditions

176

at open (sea) boundaries and conditions at solid (land) boundaries. At open bound-

177

aries sea level η, radiation conditions, or a combination of these are generally set. In

178

addition, discharge may be specified for river or other inflow. At land boundaries, the

179 180

M

175

b ) = 0 where n b is the unit normal. normal component of velocity vanishes so that (u · n If the flow is viscous, then stress conditions also need to be specified. Either a stress

condition is specified by Eq. 11 or u|h = 0 and the bottom boundary layer needs to be

182

resolved by the vertical grid placement.

ED

181

The 2-dimensional vertically-averaged form of Eq. 9 is similar with u replaced by

184

the depth-averaged value and the vertical component of stress (fourth term) replaced

185

by (τ s − τ b )/ρH where τ b is the bottom stress, τ s is the surface stress which is

neglected here, ρ is a reference density, and H(x, y, t) is the total water depth given by

CE

186

PT

183

187

H = η − h where h(x, y) is the bottom elevation measured from a reference elevation.

For the derivation, see the development of the depth-averaged equations in [18] or

189

many other publications.

AC

188

190

191

The bottom stress τ b is given by τb = Cb |u|u = γu ρ

(z = h),

where Cb is a bottom friction coefficient and γ is defined by this equation. 9

(11)

ACCEPTED MANUSCRIPT

192

2.2. Form drag Form drag that is modeled as a subscale process (eq. 5) can be included into the

194

momentum equation in a general sense with the body force term F in Eq. 9. Then F

195

includes Fd which is the form drag per unit volume which is given by

Fd = (1/2)λf Cd |u|u

CR IP T

193

(12)

where λf is the frontal area Af of the object divided by the averaging volume. Note that

197

both Af and the averaging volume are important for scaling from the discrete turbine

198

scale to the larger-scale model. Moreover, some care must be taken as the volume-

199

averaged model velocity in eq. 12 will not be the same as the reference velocity in eq.

200

5.

AN US

196

From a model perspective, the obvious averaging volume is the volume of a compu-

202

tational cell- an element volume when using finite element methods. In the following,

203

A and V will refer to the horizontal area and volume of the computational cell.

204

2.2.1. 3D case

M

201

Expand Af = W ∆z and V = A∆z in the expression for λf so that

205

206

where W (z) is width of the object and A is element area at depth z.

207

2.2.2. 2D case

The 2-dimensional equations are derived by depth integrating the 3-dimensional

PT

208

¯ d is the depth-averaged form drag approximated as equations with the result that F

CE

209

211

212 213

¯ d = (1/2)λf Cd |U|U = (cf d /H)|U|U F

is non-dimensional. Hence, the form drag coefficient and bottom friction coefficient in a 2-dimensional

formulation can be merged to give τ /(ρH) = [(Cb + cf d )/H]|U|U

214

(14)

where U is the depth-averaged velocity and cf d = (1/2)(Af /A)Cd by definition and

AC

210

(13)

ED

Fd = (1/2)(W/A)Cd |u|u

where Cb is the non-dimensional bottom friction coefficient. 10

(15)

ACCEPTED MANUSCRIPT

215

2.3. Discretization and Solution Methods The objective of the time and space discretization methods is to optimize the model

217

for accuracy, stability, and efficiency. The momentum and free surface equations are

218

discretized in time using semi-implicit methods in order to remove the most restrictive

219

stability constraints on gravity waves and vertical stress [2, 4]. The material derivative

220

in the momentum equation is discretized using semi-Lagrangian methods to remove

221

stability constraints on advection [2, 27]. The Coriolis term is discretized with a 3rd-

222

order Adams-Bashforth scheme to provide stability [28]. In addition, the spatial dis-

223

cretization uses an unstructured grid with either triangular or quadrilateral elements for

224

flexibility in representing irregular multiply-connected domains [25]. More details are

225

provided in the subsections that follow.

226

2.3.1. Time Discretization

Using the θ-method, the discretized equations for the free surface and momentum

227

and

 Z η  Z η η n+1 − η n n+1 n + ∇ · θ( u dz) + (1 − θ)( u dz) = 0 ∆t h h

(16)

ED

229

can be written as [26]

M

228

AN US

CR IP T

216

PT

un+1 − u∗ ∆t

ˆ × un + g∇(θη n+1 + (1 − θ)η n ) + fz −

∂ ∂un+1 (Av ) − ∇ · (Ah ∇un ) + Fn+1 = 0 ∂z ∂z

(17)

where n denotes the time level tn = n∆t, u∗ is the velocity at the foot of the La-

231

grangian trajectory, the vertical component of stress and form drag are treated implic-

232

itly, and quadratic stress terms in bottom friction and form drag are approximated as

CE

230

(|u|u)n+1 = |un |un+1 .

AC 233

234

For θ > 1/2, the gravity wave terms are stable [4]. The remaining terms with

235

stability constraints are the horizontal stress and the Coriolis terms. In practice, the

236

horizontal stress is only used when resolved horizontal boundary layers are present.

237

In that case, the term is sub-stepped in time if necessary in order to satisfy a stability

11

ACCEPTED MANUSCRIPT

constraint [12]. As pointed out in [11], the stability of the Coriolis term has 2 aspects –

239

one relating to the time discretization and one to the spatial discretization (next section).

240

The standard explicit treatment of this term is not stable [28] and the instability is

241

masked when significant dissipation is present, but becomes evident once beyond the

242

shelf break. For our purposes, a stable 3rd-order Adams-Bashforth scheme is used

243

[28]. Other models [6] have had success with a 2nd-order Adams-Bashforth scheme

244

with modified coefficients.

245

2.3.2. Spatial Discretization

CR IP T

238

In each finite element, the dependent variables for sea level, η, and normal velocity

247

on an edge, un , are approximated as a constant and a linear function, respectively.

248

This particular element is known as the Raviart-Thomas element of lowest order [20].

249

Hence sea level is a piecewise constant globally and velocity is a piecewise linear

250

function globally. The approximation for velocity is somewhat unusual in that un is a

251

scalar whereas the approximation function is a vector.

AN US

246

This particular approximation method has a number of important advantages.

253

• With piecewise constant η, the free surface and continuity equations become

255

finite volume expressions that conserve mass exactly.

ED

254

M

252

• Normal velocity is continuous on each edge which simplifies the calculation of fluxes and Lagrangian trajectories. There is no need to evaluate a discontinuous

257

normal velocity on an edge such as with Discontinuous Galerkin Methods.

258

• All terms in the equations can be calculated directly from the velocity approximation. In particular, the Coriolis term is calculated directly from a vector prod-

CE

259

PT

256

260

uct between the unit vector in the vertical and the vector basis functions for

261

velocity (Fig. 2 in [26]). Thus, there is no need to reconstruct the full velocity field such as with irregular grid finite difference methods that use circumcenters

263

in an orthogonal grid.

AC 262

264

• There is no constraint on element shape except to satisfy accuracy requirements

265

(equilateral elements are most accurate). There is no limit on element corner

266

angles. 12

ACCEPTED MANUSCRIPT

267

On the other hand, there is a serious shortcoming to this approach. A direct ap-

268

plication of the Galerkin finite element methodology leads to element mass matrices

269

that are fully populated nv × nv matrices, where nv is the number of vertices in an element. After assembly, a global sparse matrix for velocity must be solved which re-

271

sults in a very inefficient model. Fortunately, there are at least 2 ways to diagonalize

272

the matrix without reducing the accuracy significantly. One of these ways, suggested

273

in [1], results in a diagonal matrix with the distance between adjacent circumcenters

274

on the diagonal; the other way results in a diagonal matrix with the normal distance

275

to an edge between adjacent centroids on the diagonal. The latter is similar in concept

276

to the Integrated Compartment Method proposed long ago [32] and reduces to the cir-

277

cumcenter formulation for equilateral and square elements. This centroid approach is

278

used in the model described here. A more detailed description of these approaches is

279

found in [26].

AN US

CR IP T

270

In the horizontal plane, the model is discretized as described above. In the ver-

281

tical, the element edges are defined by vertical lines extending from the nodes at the

282

free surface to the bottom and the elements are brick- or pie-shaped depending on

283

whether quadrilateral or triangular elements are used in the horizontal. The vertical

284

discretization uses terrain-following σ coordinates. Vertical velocity is approximated

285

as a piecewise constant on the upper and lower faces of each element. Horizontal ve-

286

locity can be discretized at the midpoint on the horizontal edges of each element or

287

the face-centers of the layers. The former uses linear basis functions so the momentum

288

equation reduces to a standard finite element tridiagonal system along each vertical line

289

connecting the midpoint of the edges. In the case of layers, velocity is constant on each

CE

PT

ED

M

280

290

face and a finite difference approach is used to form a tridiagonal system. In this study,

291

layers are chosen since they are more compatible with volume averaging. Using these methods, the discretized equation for the free surface can be written as

AC

292 293

294

[26]

Ae

Z Z Nes Nes X X ηen+1 − ηen = −θ li ( un+1 dz) − (1 − θ) l ( unj(i) dz) i j(i) n n ∆t H H i=1 i=1 j(i) j(i)

(18)

where subscript e = 1, . . . , Ne is the element index, Ne is the number of elements, 13

ACCEPTED MANUSCRIPT

Ae is the area of element e in the x-y plane, li is the edge length of side i, Nes is the

296

number of element sides, uj(i) is the normal velocity on edge j defined by side i of

297

element e, H is the water depth, and a suitable convention is used to determine normal

298

direction uniquely [5]. This equation states that the volume change at the surface is the

299

net volumetric flux through the sides.

CR IP T

295

For this element approximation, there is one degree of freedom for sea level η on

301

each element and one degree of freedom for the normal component of velocity uj on

302

each edge j. Hence it is more convenient to express the momentum equation in terms

303

of the normal velocity component. Since (17) is invariant under solid rotation in the

304

horizontal plane, the momentum equation for the normal velocity component can be

305

written

AN US

300

un+1 − u∗j ∂ ∂η n+θ ∂uj n+1 j − (Av ) + Fjn+1 = −g + Cjn + ∇ · (Ah ∇u)nj (19) ∆t ∂z ∂z ∂n where uj is the normal component of velocity positive outwards on side j, the body

307

force F n+1 is evaluated as γ|un |un+1 where γ are the coefficients used, Cjn is the

M

306

Coriolis term evaluated on side j, the superscript denotes the time level, η n+θ =

309

θη n+1 + (1 − θ)η n , and ∂/∂n is the gradient operator normal to the cell edge. Note

ED

308

that horizontal friction is generally neglected because it is usually small compared to

311

vertical friction. If it is included, the discretization follows standard methods used in

312

FEM. For the example presented later, horizontal friction is neglected.

PT

310

The Coriolis term is evaluated with the 3rd-order Adams-Bashforth scheme as

CE

313

314

316

317

318

f (23vjn − 16vjn−1 + 5vjn−2 ) 12

(20)

where vjn ≡ (ˆ z × u)nj is defined from the corresponding term in eq. 17 and can

be thought of as the tangential component of velocity positive in a counter-clockwise

AC

315

Cjn =

sense around each element. After assembly of all the element contributions, the discretized momentum equa-

tion for the normal velocity on a side can be written as [26] Aun n+1 = G − gθ∆t∆η n+1 Z 14

(21)

ACCEPTED MANUSCRIPT

where matrix A contains all the implicit contributions on the left hand side of (19), un

320

is a vector of the normal velocity components uj on the edges, G is a vector containing

321

the explicit terms in (19), ∆η n+1 is the difference between element values adjacent to

322

an edge, and Z is a vector of depth increments. These terms are defined in more detail

323

in [26]. In general, A is a large sparse matrix that is block tridiagonal at each velocity

324

node in the horizontal grid. This matrix is inverted by diagonalizing in the horizontal

325

as described earlier and solving the tridiagonal system in the vertical at each node [26].

326

2.3.3. Solution Method

CR IP T

319

In general, the system of equations given by (18) and (19) are solved by first solving

328

for un by solving the tridiagonal system in (19) and substituting this result into (18) to

329

form a discrete wave equation that involves only η at the n + 1 time level. This equa-

330

tion is positive definite and diagonally dominant, and is solved by standard conjugate

331

gradient methods.

AN US

327

However, the treatment of wetting and drying can be improved by writing the free

333

surface equation (Eq. 18) in finite volume form and adopting the method proposed by

334

Casulli [3]. Then the momentum equation is solved for un n+1 and substituted into the

335

finite volume free surface equation, so that the resulting discrete wave equation can be

336

written in compact vector form as [3]

V(ζ) + Tζ = b

(22)

PT

ED

M

332

337

where ζ is a vector of sea level with components ζi = ηin+1 , V is a vector of volumes

338

with components Vi (ηin+1 ), and T is the flux matrix arising from the second term in (18). As shown by Casulli [3], an efficient Newton-type algorithm for solving Eq. 22

340

is given by

AC

CE 339

341

342

343

[P(ζ m ) + T]∆ζ m+1 = [V(ζ m ) + Tζ m − b]

where m is the iteration index, P is the wet area, and ∆ζ m+1 = ζ m+1 − ζ m .

(23)

Once ζ is known, the values for η n+1 are backsubstituted into eq. 21 to solve for

un n+1 .

15

ACCEPTED MANUSCRIPT

344

3. Tidal Power Example From the development thus far, the tidal power problem has been reduced to finding

346

the appropriate values for the drag coefficient Cd defined by a known reference veloc-

347

ity, and a suitable closure for the spatial dispersion term in (4). From the perspective of

348

tidal turbines used for power generation, two coefficients are necessary: a thrust coef-

349

ficient CT (which is equivalent to a drag coefficient) and a power coefficient Cp , along

350

with a useable reference velocity. The double-averaged velocity produced by the model

351

is the obvious choice for reference velocity. What remains then is how to determine

352

the values for the coefficients based on this reference velocity. The spatial dispersion

353

term appears to be mainly significant in the wake region and is not considered in this

354

example.

In the context of turbines, the drag coefficient, Cd , is written as the thrust coeffi-

355 356

AN US

CR IP T

345

cient, CT . Using this notation, the drag force and power can be written (24)

357

where At is the frontal area of the turbine normal to the direction of flow, and ur is the

358

free stream reference velocity, and

359

ED

M

F = (1/2)ρAt CT |ur |ur

where ur is the corresponding free stream reference speed. Drag from the structure

360

supporting the turbine as well as low speed cut-in and high speed cut-off are neglected.

361

For an actuator disk in a free stream flow, values for the drag and power coeffi-

(25)

CE

PT

P = (1/2)ρAt Cp ur 3

cients can be derived using conservation of mass and momentum. The coefficients

363

at maximum energy extraction were derived independently by Lanchester (1915), by

364

Betz (1920), and by Zhukowsky (1920) and are generally known as the Lanchester-

AC

362

365

Betz limit. For this case, CT = 8/9 and Cp = 16/27 where the reference velocity

366

is the free stream velocity. Garrett and Cummins [9] also derive these results and in

367

addition derive a correction term for flows constrained by lateral boundaries. For that

368

case, the power can be larger due to the increased pressure drop across the turbines.

369

In practical terms, CT and Cp are defined for various turbines by their manufacturers. 16

ACCEPTED MANUSCRIPT

Typical values for CT range from 0.4 to 0.8 and for Cp range from 0.3 to 0.5 based

371

on free stream velocity as the reference velocity. For a simple straight channel with

372

constant depth and width and discharge, the free stream velocity is defined upstream.

373

However, for more complicated geometry and transient flows, defining a free stream

374

velocity is not possible in general. Additional complications arise with shear flows

375

such as the vertical shear in tidal flows and the determination of averaging volume for

376

multi-turbine deployments.

CR IP T

370

In most cases, the drag coefficient for various geometric objects and thrust coeffi-

378

cients for turbines use free stream velocity as the reference. Then the relation between

379

free-stream and volume-averaged velocity needs to be determined and the coefficients

380

adjusted accordingly. We can define a new set of coefficients CT∗ and Cp∗ as

AN US

377

F = (1/2)ρAt CT∗ |U|U 381

where U is the volume-averaged velocity, and

M

P = (1/2)ρAt Cp∗ U 3

382

(26)

(27)

where U is the volume-average speed.

Note that the model velocity is discretized on the edges so care must be exercised

384

to ensure that the volume average of the edge-based equations is equal to the volume-

385

averaged equations. For the linear terms this is a straightforward average. The ad-

386

vective fluxes through the edges are defined the same in both cases. The only term

387

needing special attention is the form drag term in eq. 26. Approximating the term as

388

(1/2)ρAt CT∗ |U |u on each edge results in an equivalent approximation for F.

CE

PT

ED

383

389

In a previous study [21], some of the issues with volume-averaged velocities were

390

addressed by using a detailed CFD model to calculate the flow around a turbine and then volume-averaging the resultant velocity to determine the parameters to be used in

392

a larger-scale coastal ocean model (COM). A few of the more important results that

393

relate to the present study are

AC 391

394

1. Both CT∗ and Cp∗ were robust parameters based on using volume-averaged veloc-

395

ity as the reference velocity. That is, they were essentially constant as a function 17

ACCEPTED MANUSCRIPT

396

of velocity when operating at maximum power generation. Hence, these param-

397

eters could be determined for the maximum velocity and used over the entire

398

range of tidal velocity. 2. The averaging volume must be at least twice the turbine diameter on each edge

400

to maintain accuracy between the small-scale CFD model and large-scale COM.

401

In a physical sense, the volume needs to be larger than the detailed flow varia-

402

tions generated around the turbine because the latter model cannot resolve them

403

nor does it contain the right physics at the turbine scale. One implication is

404

that the vertical variation in velocity for the COM must be averaged over the

405

total averaging-volume because the grid spacing in the vertical is usually much

406

smaller that the turbine diameter. Toward this end, local coefficients in each layer

407

are defined from the local velocity in each layer such that the thrust and power

408

sum to the correct total for the turbine.

AN US

CR IP T

399

3. Ultimately, the consistency between the CFD model and COM was dependent on

410

how accurately the COM could produce the same volume-averaged velocity as

411

the CFD model. For single turbines in a simple channel geometry, the agreement

412

was good. For more complicated field-scale geometry, differences in physics and

413

turbulence closures between the models led to greater errors.

ED

M

409

In the following section which describes a field-scale problem, other options are

415

considered for determining coefficients when detailed CFD calculations may or may

416

not be available. This example uses the same grid and sites as were used in [21]. The

417

location is Minas Passage in the upper Bay of Fundy.

PT

414

3.1. Minas Passage

CE 418

Minas Passage connects Minas Basin with the upper Bay of Fundy (Fig. 1) and is

420

well known for its large tidal currents and range. This area has considerable potential

AC

419

421

for tidal power generation and has been the subject of many studies from a theoretical

422

and practical perspective.

423

The model domain includes the Gulf of Maine, Bay of Fundy, Minas Passage,

424

and Minas Basin (Fig.1). This area is located on the east coast of North America

425

between approximately 44◦ and 46◦ North latitude. The resonant period of this system 18

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

ED

Figure 1: Geometry and bathymetry of Bay of Fundy. AO is Atlantic Ocean, GB is Georges Bank, GoM is Gulf of Maine, BoF is Bay of Fundy, MP is Minas Passage, MB is Minas Basin, B is Boston, and CC is Cape Cod.

is slightly longer than the period of the dominant tidal constituent, M2 , [8, 10] and this

427

area is well known for the large tide range in the upper bay. The node point for the

428

oscillation is near the shelf break east of Georges Bank; however, this is complicated

CE

PT

426

429

by other resonances parallel to the shelf that extend down to Boston and Cape Cod. The dominant tidal constituent is the M2 and the amplitude varies from 0.5 m at

431

the open boundary to over 6 m in Minas Basin. In a previous study [30], 8 constituents

AC

430

432

are used in the open boundary conditions– M2 , S2 , N2 , K2 , K1 , O1 , P1 , and Q1 .

433

Most of these constituents are small compared to the M2 (particularly the diurnal con-

434

stituents) but were retained in order to assess the accuracy of the model against tide and

435

current observations. In general, the model reproduced the observed tidal harmonics

436

and ADCP current meter observations so that the model is sufficiently accurate for the 19

ACCEPTED MANUSCRIPT

437

evaluation of tidal power potential [30]. During flood tide, a flow separation develops at Cape Split at the entrance to Minas

439

Passage with a high velocity core in the northern part of the passage and a counter eddy

440

in the southern part (Fig. 2). Peak velocity approaches 4m/s in the deep part of the

441

passage outside the eddy. The reattachment length of the eddy and its strength increase

442

until peak flood.

PT

ED

M

AN US

CR IP T

438

CE

Figure 2: Minas Passage at flood tide. CS identifies Cape Split.  denote the turbine deployment sites. Velocity vectors have constant length and are decimated by 5.

During ebb tide, the velocity is more uniform across the passage and an eddy devel-

AC

443

444

ops southwest of Cape Split. At the start of the ebb tide, the eddies that were created

445

during the flood are advected out Minas Passage. The passage of these eddies is re-

446

peatable over the tidal cycles and can be seen in the current meter data as well as the

447

model results.Overall, the flood/ebb behavior is consistent with other published results

448

[10, 22]. 20

ACCEPTED MANUSCRIPT

The original model grid was constructed from unstructured triangles that vary in

450

edge length from 12 km on the open boundary to 70 m in Minas Passage. This grid

451

was modified by inserting four 200 m by 200 m square grids at the simulated turbine

452

deployment sites. The square grids are composed of either a 3x3 or 5x5 array of square

453

elements (66.7 m or 40 m edge length). The turbine was placed at the center of the

454

square grids (Fig. 3). This particular grid configuration was created so that CFD cal-

455

culations could be made on the square grids. This same grid is used in the results

456

presented here.

CR IP T

449

A single turbine was placed at each test site. All the turbines were 16 m in di-

458

ameter and placed 15 m from the bottom of the water column. Note that each set of

459

square grids was rotated to align with the maximum flood velocity to simplify the CFD

460

calculations in a previous study.

AC

CE

PT

ED

M

AN US

457

Figure 3: Grid around the turbine sites in Minas Passage.

21

ACCEPTED MANUSCRIPT

For the present example, the small constituents are neglected and only the M2

462

is retained. In total, there are 69583 vertices and 135101 triangular and quadrilateral

463

elements. On a desktop computer using openmp with 5 cores, the model runs 270 times

464

faster (2D) and 60 times faster (3D) than simulated time. The calculations presented

465

for this study are 2D and a set of 3D calculations can be found in [21].

466

3.2. Estimation of coefficients

CR IP T

461

The estimation procedures adopted here use the velocity in the absence of a turbine

468

at the turbine site as an approximation for the free stream velocity. Then the coefficients

469

CT and Cp for the Lanchester-Betz limit are used to calculate thrust and power where

470

the reference velocity is the approximation for free stream velocity. From the model

471

perspective, the volume-averaged velocity is smaller than the free stream velocity when

472

form drag is applied so the coefficients must be adjusted to recover the same thrust and

473

power values. This adjustment is made in two ways: (1) the coefficients are adjusted

474

iteratively until the thrust and power from the initial estimate is reached, and (2) the

475

coefficients are mapped from a free stream reference to a cell averaged velocity using

476

two analytical approaches.

M

AN US

467

The general procedure is to spin up the model without form drag until it reaches a

478

periodic steady state. At a time near maximum flood currents at the turbine sites, a ref-

479

erence thrust and power are calculated based on the approximate free stream velocity

480

(case 0). In case 1, form drag from the turbines is introduced and the model is spun

481

up to the same time as before. Thrust and power are calculated based on the volume-

482

averaged velocity which gives values that are smaller than the free stream calculations.

483

The coefficients CT and Cp are then modified iteratively to estimate CT∗ = CT u2f s /u2m

484

and Cp∗ = CT u3f s /u3m where uf s is the free stream speed and um is the volume aver-

485

aged speed at iteration m. Iterations were continued until the errors in thrust and power

CE

PT

ED

477

were small (< 1%). Typically one or two iterations are sufficient. In case 2 and 3, the

487

coefficients CT and Cp are mapped into CT∗ and Cp∗ and thrust and power calculations

488

are compared to the reference calculations.

AC 486

489

In this study, the turbines are assumed to be omnidirectional; that is, the incident

490

flow is always normal to the face of the turbine. For a strongly bi-directional flow such

22

ACCEPTED MANUSCRIPT

491

as encountered here, the differences with fixed turbines are small and occur mainly

492

when the tide is turning and the currents are weak. In the first case using the Lanchester-Betz limit, CT = 8/9 and Cp = 16/27 are

494

referenced to the free stream velocity. The resulting values for thrust and power at the

495

4 sites is listed in Table 1. The values in the table are taken at 88.5 hrs from the start of

496

the simulation, a point near maximum flood currents.

CR IP T

493

Table 1: Coefficients, thrust, and power for case 1. The values for CT , Cp , T , and P are based on free stream velocity. CT∗ and Cp∗ are referenced to the volume-averaged velocity.

CT 0.8889 0.8889 0.8889 0.8889

Cp 0.5926 0.5926 0.5926 0.5926

CT∗ 1.0064 0.9222 0.9056 0.9904

Cp∗ 0.7139 0.6262 0.6094 0.6969

T (kN) 829 717 693 759

P (kW) 1678 1354 1287 1472

AN US

site 1 2 3 4

Generally a monotonic convergence would be expected because of negative feed-

498

back between changes in the coefficients and thrust and power. However, that was not

499

the case at site 1, indicating that there is a weak interaction between the turbines even

500

with the large spacing between them. The two upstream turbines (2,3) change the flow

501

slightly and affect the two downstream turbines (1,4). Here, the upstream direction is

502

taken at maximum flood velocity which is the time when maximum power is gener-

503

ated. In retrospect, a better estimate of free stream velocity at a particular site would

504

be when the other three turbines were in place. Recalculations result in larger thrust

505

and power at site 3 (≈ 3%) and 4 (≈ 4%), and smaller changes at sites 1 and 2 (Table

506

2). For the present case, this new estimate results in marginal changes compared to

507

the first estimate but is used as the reference thrust and power in this analysis. Note

508

that as the number of turbines increases, this procedure becomes more time consuming

CE

PT

ED

M

497

and tedious since a separate simulation is needed to define the approximate free stream

510

speed at each site.

AC 509

511

The values for CT∗ and Cp∗ follow a trend that has been defined previously in sim-

512

ulations with a single turbine in a simple channel flow. For a specific device with a

513

known thrust and power, the volume-averaged velocity becomes smaller as the averag-

23

ACCEPTED MANUSCRIPT

Table 2: Coefficients, thrust, and power for case 1. The values for CT , Cp , T , and P are based on free stream velocity with the other 3 turbines in place. CT∗ and Cp∗ are referenced to the volume-averaged velocity.

CT 0.8889 0.8889 0.8889 0.8889

CT∗ 1.0021 0.9221 0.9235 1.0266

Cp 0.5926 0.5926 0.5926 0.5926

Cp∗ 0.7094 0.6262 0.6309 0.7355

T (kN) 822 717 700 771

P (kW) 1659 1353 1307 1507

CR IP T

site 1 2 3 4

ing volume becomes smaller. Hence for the same thrust and power, the coefficients CT∗

515

and Cp∗ are larger with smaller averaging volume. As a rough measure, the coefficients

516

increase by approximately 4 % when a square element with edge length of 4 times the

517

turbine diameter is reduced to 2 times the diameter. There is a similar trend when depth

518

is reduced. In the present case, site 1 and 4 are contained in elements that are 40 m

519

on an edge and in water depths of 36 and 37 m. On the other hand, sites 2 and 3 are

520

contained in elements that are 66.7 m on an edge and in water depths of 55 and 48 m,

521

respectively. Using the rough scaling above, an increase of 8 % in the coefficients from

522

sites 2 and 3 to sites 1 and 4 is reasonable.

M

AN US

514

In the second case, the coefficients CT and Cp are mapped to CT∗ and Cp∗ using a

523

525

ED

CT = 4a(1 − a)

(28)

where a = (u0 − ut )/u0 is the axial induction factor, u0 is free stream speed and ut is the speed at the turbine disk. To relate u0 to ut , this expression can be written as

CE

526

method from [23]. Using actuator disk theory, CT is written as

PT

524

u0 =

(1 +

√

2 ut 1 − CT )

(29)

Waldman et al. [23] assume that this relation also holds for relating the free stream

AC

527

528

velocity to the area averaged velocity u in the grid cell in the plane of the turbine so

529

that

u0 =

(1 +

√

2 u 1 − CT )

24

(30)

ACCEPTED MANUSCRIPT

where  = At /Ac corrects for the fact that the disk only occupies a part of the cross

531

section. Ac is the cross section area of the computational cell in the plane of the turbine.

532

Then eq. 30 can be substituted into eqs. 24 and 25 where u0 = ur . Comparing with

533

eqs. 26 and 27, the coefficients can be corrected to use the average velocity, resulting

534

in

Cp∗ = Cp



and

535

2 √ 1 + 1 − CT

2

2 1 − CT

3

√

AN US

= CT



CT∗

CR IP T

530

1+

(31)

(32)

536

These equations apply to the case where  is small so that Ac is sufficiently large

537

that it contains all the flow variations induced by the presence of the turbine. In a

538

previous study [21], we found that accuracy deteriorated when a edge of the averaging

539

cube was less than twice the turbine diameter.

The calculated coefficients, thrust, and power for case 2 are listed in Table 3. Com-

541

paring this table with the results after numerical convergence in Table 2, the values for

542

the coefficients as well as thrust and power are all smaller. Sites 2 and 3 have the largest

543

elements that are 66.7 m on an edge or 4.1 times the turbine diameter and hence have

544

the largest Ac and smallest . At these sites the results are in reasonable agreement with

545

the other results with differences less than 1 or 2 %. Sites 1 and 4 have the smallest

546

elements that are 40 m on an edge or 2.5 times the turbine diameter and have smaller

547

values for Ac and larger . The differences at these sites are larger, approximately 5 %

548

for thrust and 8 % for power. Past experience would indicate that the differences would

549

increase significantly at even smaller averaging volumes.

CE

PT

ED

M

540

In the third case, another mapping is derived from continuity requirements between

AC

550 551

the upstream flow and the sum of the flow through and bypassing the turbine. From

552

continuity

Ac u = (Ac − At )ub + At ut

25

(33)

ACCEPTED MANUSCRIPT

Table 3: Coefficients, thrust, and power for case 2. The values for T and P are based on free stream velocity. CT∗ , Cp∗ , T ∗ , and P ∗ are referenced to the volume-averaged velocity.

CT∗ 0.9438 0.9108 0.9141 0.9458

Cp∗ 0.6483 0.6147 0.6180 0.6504

T ∗ (kN) 784 709 699 730

P ∗ (kW) 1543 1332 1303 1388

T (kN) 822 717 700 771

P (kW) 1659 1353 1307 1507

CR IP T

site 1 2 3 4

In the limit as  becomes small, the bypass flow ub can be approximated as the free

554

stream flow. Then a relation between the cell averaged and free stream flow can be

555

derived as

AN US

553

−1  p  u u0 = 1 − (1 − 1 − CT ) 2

(34)

−2  p  CT∗ = CT 1 − (1 − 1 − CT ) 2

(35)

−3  p  Cp∗ = Cp 1 − (1 − 1 − CT ) 2

(36)

and

ED

557

M

and the coefficients are corrected as

556

These corrections result in larger coefficients and hence larger thrust and power

559

(Table 4). In general, these results are in better agreement with the reference thrust and

560

power from case 1. The differences in thrust range from 3 % at sites 1 and 4 with 40 m

561

edge lengths to less than 1 % at sites 2 and 3 with 66.7 m edge length. Differences in

CE

PT

558

562

power range from 5 % to less than 1% at the same sites. In comparing the time series of power at the 4 sites, 2 of the sites (1 and 4) have

564

relatively large differences between the cases and 2 of the sites (2 and 3) have relatively

AC

563

565

small differences between the cases. For comparison, the reference calculation (case 0),

566

the results after numerical convergence (case 1) and the 2 results after mapping (case 2

567

and 3) are compared at sites 1 and 3 (Figure 4). Power at site 1 (the upper set of lines)

568

displays the largest differences in the 4 results, whereas power at site 3 (the lower set of

569

lines) displays the smallest differences. The results at site 4 are similar to those at site 26

ACCEPTED MANUSCRIPT

Table 4: Coefficients, thrust, and power for case 3. The values for T and P are based on free stream velocity. CT∗ , Cp∗ , T ∗ , and P ∗ are referenced to the volume-averaged velocity.

CT∗ 0.9703 0.9217 0.9264 0.9735

Cp∗ 0.6759 0.6257 0.6305 0.6792

T ∗ (kN) 801 716 704 739

P ∗ (kW) 1594 1349 1318 1415

T (kN) 822 717 700 771

P (kW) 1659 1353 1307 1507

CR IP T

site 1 2 3 4

1 but with slightly smaller amplitude. The results at site 2 are similar to those at site 3

571

but with larger values for power. The difference in results at the various sites probably

572

gives a good sense of the variability between different approximation methods. What

573

is really needed is a rigorous comparison between data from an installed turbine and

574

numerical model simulations using the different approximation methods. I am unaware

575

of any such data that is available.

576

4. Concluding Remarks

M

AN US

570

In the formulation of this coastal ocean model, an attempt is made to optimize the

578

model for accuracy, efficiency, and robustness. Toward this end, a semi-implicit time

579

integration removes all the major stability constraints, a semi-Lagrangian (SLM/ELM)

580

advection approximation removes Courant number restrictions, a finite element dis-

581

cretization gives a compact stencil, and solution of a wave equation at the discrete level

582

results in good efficiency. The advection approximation uses quadratic interpolation at

583

the foot of the trajectory which has negligible dispersion compared to linear interpola-

584

tion methods.

CE

PT

ED

577

585

The model also employs a volume averaging method to approximate subgrid ob-

586

jects as body forces. For a coastal scale model, the shallow water equations are modified to include a form drag and spatial dispersion term. The methods are applied to

588

tidal power generation from turbines that are subgrid in scale at a study site in Minas

589

Passage at the upper end of Bay of Fundy. Some options are presented for estimating

590

the thrust and power coefficients, including a numerical convergence procedure and 2

591

methods of mapping from a free stream velocity reference to a volume averaged veloc-

AC 587

27

M

AN US

CR IP T

ACCEPTED MANUSCRIPT

ED

Figure 4: Variation in power between different approximation methods. The thick line denotes site 1 and the thin line denotes site 3. A different color is used for each of the 4 cases. In the legend, C denotes the case (C0 is case 0, etc.) and S denotes the site (S1 is site 1, etc.).

ity reference. Although the two mappings are ostensibly for the same conditions where

593

 is small and the flow bypassing the turbine is approximately the free stream veloc-

594

ity, the two sets of calculated coefficients differ. All sets of results are in reasonable

595

agreement, but the smallest difference is between the the coefficients from numerical

CE

PT

592

596

iteration (case 1) and the mapping based on continuity (case 3). In order to simplify the analysis, only the dominant tidal constituent (M2 ) was used

598

in the boundary forcing. Using the full set of constituents would introduce fortnightly

599

and monthly modulations in the tides and lead to greater variability in power output.

600

The results presented here would still apply under the same assumptions. For a more

601

complete analysis, the assumption that the turbines operate with coefficients given by

602

the Lanchester-Betz limit would have to be dropped and field-scale data used instead.

AC

597

28

ACCEPTED MANUSCRIPT

For this analysis, the locations of the turbines are at sites where tidal turbines will be

604

tested in the future. Choosing other sites depends on a number of factors including but

605

not limited to the maximum velocity, the water depth, navigation hazards, turbulence

606

intensity, accessibility for maintenance, access to the power grid, etc. This evaluation

607

is beyond the scope of this study.

CR IP T

603

In these examples, the spatial dispersion term was neglected, mainly because it is

609

not clear how to calculate this term except by small-scale CFD simulations. Intuitively,

610

this term should have a role in modifying the wake and the flow that bypasses the tur-

611

bine disk. When the averaging volume is of the order of the size of the turbine (a blade

612

diameter on an edge), then the flow in the laterally adjacent elements is accelerated

613

by the spatial dispersion in addition to the surface pressure gradient. Clearly, much

614

more work needs to be done to try and understand the magnitude and importance of

615

this term.

616

Acknowledgments

AN US

608

M

The author wishes to acknowledge helpful review comments by Michael Shives,

617

Phillip Gillibrand, and Vincenzo Casulli.

619

References

ED

618

[1] Baranger, J., Maitre, J.F., Oudin, F., 1996. Connection between finite volume and

621

mixed finite element methods. Mathematical Modelling and Numerical Analysis

622

30, 445-465.

PT

620

[2] Casulli, V., 1987. Eulerian-Lagrangian methods for hyperbolic and convection

624

dominated parabolic problems. In: Taylor, C., D.R Owen, Hinton, E., (Eds),

625

Computational Methods for Non-linear Problems, Pineridge Press, Swansea, pp.

AC

CE 623

626

239–268.

627

[3] Casulli, V., 2009. A high-resolution wetting and drying algorithm for free-surface

628

hydrodynamics. Int. J. Numer. Meth. Fluids 60, 391–408. DOI: 10.1002/fld.1896

29

ACCEPTED MANUSCRIPT

629

[4] Casulli V., Cattani E., 1994. Stability, accuracy, and efficiency of a semi-implicit

630

method for three-dimensional shallow water flow. Computers Math. Applications

631

27(4), 99-112. [5] Casulli, V., Walters, R. A., 2000. An unstructured grid, three-dimensional model

633

based on the shallow water equations. International Journal for Numerical Meth-

634

ods in Fluids 32 (3), 331-348.

CR IP T

632

[6] Fringer, O.B, Gerrirsen, M., Street, R.L., 2006. An unstructured grid, finite vol-

636

ume, nonhydrostatic, parallel coastal ocean simulator. Ocean Modelling 14, 139-

637

173.

AN US

635

638

[7] Finnigan, J.J., Shaw, R.H., 2008. Double-averaging methodology and its applica-

639

tion to turbulent flow in and above vegetation canopies. Acta Geophysica 56(3),

640

534–561.

[8] Garrett, C., 1972. Tidal resonance in the Bay of Fundy and Gulf of Maine. Nature 238, 441–443.

642

643

[9] Garrett, C., Cummins, P., 2007. The efficiency of a turbine in a tidal channel. J. Fluid Mechanics 588, 243-251. doi:10.1017/S0022112007007781

ED

644

645

M

641

[10] Greenberg, D.A., 1979. A numerical model investigation of tidal phenomena in the Bay of Fundy and Gulf of Maine. Marine Geodesy 2(2), 161–187.

PT

646

[11] Ham, D.A, Kramer, S. C., Stelling, G. S., Pietrzak, J. D., 2007. The symmetry and

648

stability of unstructured mesh C-grid shallow water models under the influence

CE

647

of Coriolis. Ocean Modelling 16, 47-60.

649

650

[12] Jankowski, J.A., (2009). Parallel implementation of a non-hydrostatic model for free surface flows with semi-Lagrangian advection treatment. Int. J. Numer. Meth.

652

Fluids 59, 1157?1179.

AC 651

653

[13] Karsten, R.H., McMillan J.M., Lickley, M.J., Haynes, R.D., 2008. Assessment of

654

tidal current energy in the Minas Passage, Bay of Fundy. Proc. of the Institution

655

of Mechanical Engineers, Part A: Journal of Power and Energy. 30

ACCEPTED MANUSCRIPT

656

[14] Le Roux D.Y., Rostand V., Pouliot B., 2007. Analysis of numerically induced os-

657

cillations in 2D finite-element shallow-water models Part I: Inertia-gravity waves.

658

SIAM Journal of Scientific Computing 29, 331-360. [15] Leendertse, J.J, 1967. Aspects of a computational model for long period water

660

wave propagation. Tech Report RM-5294-PR, Santa Monica, Rand Memoran-

661

dum.

CR IP T

659

[16] Nikora, V., McLean, S., Coleman, S.,Pokrajac, D., Walters, R., 2007. Double-

663

averaging concept for rough-bed open-channel and overland flows: Theoretical

664

background. J. Hydraul. Eng. ASCE 133, 873–883.

AN US

662

665

[17] O’Donncha,F., Hartnett, M., Plew, D.R., 2015. Parameterizing suspended canopy

666

effects in a three-dimensional hydrodynamic model. Journal of Hydraulic Re-

667

search, DOI: 10.1080/00221686.2015.1093036

668

[18] Pinder, G.F., Gray, W.G., 1977. Finite Elements in Surface and Subsurface Hydrology, Academic Press.

M

669

[19] Plew, D.R., 2011. Shellfish farm-induced changes to tidal circulation in an em-

671

bayment, and implications for seston depletion Aquaculture Environment Inter-

672

actions 1, 201–214. doi: 10.3354/aei00020

ED

670

[20] Raviart,P. A., Thomas, J. M., 1977. A mixed Finite Element Method for 2nd Order

674

Elliptic Problems. Mathematical Aspects of the Finite Element Method, Lecture

675

Notes in Math. Springer Berlin.

CE

PT

673

[21] Shives, M., Crawford, C., Hiles, C., Walters, R., 2013. Combining Numerical

677

Methods for Basin and Turbine Scales for Improved Modelling of in-situ Turbine

678

Arrays. Proc. 10th European Wave and Tidal Energy Conf., Aalborg, Denmark.

AC

676

679

[22] Sucsy, P.V., Pearce, B.R., Panchang, V.G., 2006. Comparison of two-and three-

680

dimensional model simulation of the effect of a tidal barrier on the Gulf of Maine

681

tides. Journal of Physical Oceanography 23(6), 1231–1248.

31

ACCEPTED MANUSCRIPT

[23] Waldman, S., Genet, G., Bast, S., Side, J., 2015. Correcting for mesh size de-

683

pendency in a regional models representation of tidal turbines. Proceedings of

684

the 11th European Wave and Tidal Energy Conference 6-11th Sept 2015, Nantes,

685

France.

686

CR IP T

682

[24] Walters, R.A., 2005, Coastal ocean models: two useful finite element methods, Continental Shelf Revsearch 25, 775-793.

687

[25] Walters, R. A., Gillibrand, P., Bell, R., Lane, E.M., 2010. A study of tides and

689

currents in Cook Strait, New Zealand. Ocean Dynamics 60, 1559–1580. DOI

690

10.1007/s10236-010-0353-8

AN US

688

691

[26] Walters, R.A., Hanert, E., Pietrzak, J.D., Le Roux, D.Y., 2009. Comparison of un-

692

structured, staggered grid methods for the shallow water equations. Ocean Mod-

693

elling 28, 106–117.

694

[27] Walters, R.A., Lane, E.M., Henry, R.F., 2008. Semi-Lagrangian methods for a finite element coastal ocean model. Ocean Modelling 19, 112–124.

M

695

[28] Walters, R.A., Lane, E.M., Hanert, E., 2009. Useful time-stepping methods

697

for the Coriolis term in a shallow water model. Ocean Modelling 28, 66–74.

698

doi:10.1016/j.ocemod.2008.10.004

699

ED

696

[29] Walters, R.A., Plew, D.R., 2008. Numerical modelling of environmental flows using DAM: Some preliminary results. Acta Geophysica 56(3), 918–934.

701

PT

700

[30] Walters, R.A., Tarbotton, M.R, Hiles, C.E., 2013. Estimation of tidal power potential. Renewable Energy 51, 255–262.

CE

702

[31] Wu, W., Shields Jr., F. F., Bennett, S. J,. and Wang, S. S. Y., 2005. A depth-

704

averaged two-dimensional model for flow, sediment transport, and bed topogra-

AC

703

705

phy in curved channels with riparian vegetation. Water Resources Research 41,

706

W03015.

707

[32] Yeh, G-T, 1981. Numerical solutions of Navier-Stokes equations with an inte-

708

grated compartment method (ICM). International Journal for Numerical Methods

709

in Fluids 1(3), 207-223. 32

ACCEPTED MANUSCRIPT

Figure Captions

710 711

Figure 1: Geometry and bathymetry of Bay of Fundy. AO is Atlantic Ocean, GoM

713

is Gulf of Maine, BoF is Bay of Fundy, MP is Minas Passage, MB is Minas Basin, B

714

is Boston, and CC is Cape Cod.

715

Figure 2: Minas Passage at flood tide. CS identifies Cape Split.  denote the tur-

716 717

CR IP T

712

bine deployment sites. Velocity vectors have constant length and are decimated by 5.

718

AN US

Figure 3: Grid around the turbine sites in Minas Passage.

719 720

Figure 4: Variation in power between different approximation methods. The thick

722

line denotes site 1 and the thin line denotes site 3. A different color is used for each of

723

the 4 cases. In the legend, C denotes the case (C0 is case 0, etc.) and S denotes the site

724

(S1 is site 1, etc.).

AC

CE

PT

ED

M

721

33