us-ing the local strain rate of the ‡ow …eld:
e=c2SA s
@u
@x
2
+ @v
@y
2
+1 2
@u
@y + @v
@x
2
; (46)
where A = cell area; cS is the Smagorinsky coe¢ cient. Smagorinsky’s isotropic turbulence model is classically used to approximate the contribution of subgrid-scale turbulence to the resolved ‡ow subgrid-scales by large eddy simulation (LES). With the MUSCL-Hancock …nite volume method used in this work, the numerical solu-tion of the time-averaged shallow water equasolu-tions is represented by cell averages, which shares therefore some analogy with the spatially …ltered LES. Smagorinsky’s formulation takes into account the size of the resolved eddies by the introduction of the cell area which re‡ects the local level of discretisation (the original paper of Smagorinsky (1963) treats Cartesian cells, in which 2 = x y replacesA). The coe¢ cient cS is theoretically in the range 0.13. . . 0.17 (Wang, 2003), but LES of open channel ‡ow shows that the optimal value of cS varies between 0.09. . . 0.22, the lower and upper limits corresponding to rivers and lakes, respectively (Thomas and Williams, 1995).
2.6 Linearised form of the shallow water equa-tions
In Chapter 5, the numerical model is validated with analytical test problems that are based on the linearised equations. Most analytical solutions of wind-induced circulation and oscillation are indeed available only for the linearised shallow water equations, which read:
@h
@t + @p
@x + @q
@y = 0; (47)
@p
@t = gh0
@
@x + sx bx; (48)
@q
@t = gh0@
@y + sy by: (49)
This nonconservative form does not account for the advective and di¤usive term, it is therefore limited to model ‡ow on a ‡at bottom where the wave amp-litudes are much smaller than the depth, therefore the celerity of the gravity waves is constant, c0 =p
gh0, where h0 is the local still water depth. So as to keep the whole system linear, the bottom shear stress is also formulated linearly:
bx =klinp; by =klinq; (50)
where klin is a linear friction coe¢ cient [1/T]. The de…nition of the surface shear stress s is problem dependent.
In order to validate the ‡ow solver against analytical solutions, equations (47) to (49) are adapted to the ‡ux form of (1). The ‡ux vectors f1 and g1, and the vector of the source termss become
f1;lin = 2 64
p gh0h
0 3
75; g1;lin = 2 64
q 0 gh0h
3
75; slin = 2 64
0
sx bx sy by
3
75: (51)
The Riemann solver is also modi…ed to comply with the linearised shallow water equations, see Eq. (83).
Chapter III
QUADTREE MESH
This chapter is devoted to the quadtree mesh structure on which the numerical discretisation is based. An innovative numbering scheme is developed to represent cell centres, edges and corners within a common quadtree data structure. The procedures of mesh generation and manipulation are outlined.
3.1 Structured and unstructured meshing
The mesh type chosen for the discretisation of the ‡ow domain greatly determines how the numerical solution operates on that mesh. Two approaches are tradition-ally distinguished: structured grids and unstructured meshes.
A structured grid approach stores the data associated with the grid geometry and ‡ow solution in logically ordered arrays with the same number of dimensions as the ‡ow problem. The topology is regular, therefore addressing neighbours of a cell is achieved simply by varying the indices in the respective dimension which is directly translated to memory pointers. This facilitates the evaluation of …nite di¤erence stencils and the assembly of sparse coe¢ cient matrices within implicit solution strategies, resulting in computationally e¢ cient solvers that are the easi-est to vectorise. Many benchmark problems deal with rectangular domains which are straightforward to accommodate using Cartesian grids. However, arbitrary domain shapes yield stepped boundaries which are known to disturb the solution.
As an example, the regular Cartesian grid laid over Lake Balaton (Figure 4a) su¤ers from a poor discretisation of the shorelines and a crude representation of the strait. Employing a logically rectangular, body-…tted curvilinear grid (Fig-ure 4b) overcomes both problems in this case. Regardless of the method of grid generation (Thompson et al., 1985), curvilinear grids may contain highly distorted cells and highly uneven cell sizes when the geometry is complicated, unless several, piecewise structured curvilinear blocks are patched together. Unfortunately block-structuring undermines the main virtues of curvilinear grids, because neighbour indexing across blocks and especially grid generation are complicated.
In contrast, unstructured meshes o¤er more ‡exible cell types and topologies than structured grids. Typically, two-dimensional unstructured meshes are based on triangular elements that avoid hanging nodes, i.e., edges that connect to an inner point of another edge (Figure 4c). The mesh is generated so that the domain
25
(a) Regular
(b) Curvilinear
(c) Triangular
(d) Quadtree
Figure 4. Lake Balaton discretised using common mesh types. The strait of Tihany is magni…ed on the right.
3.1. STRUCTURED AND UNSTRUCTURED MESHING 27 boundaries are resolved with a string of triangle edges, and the interior is …lled with triangles according to a prescribed, possibly variable density. Meshes that allow arbitrary element types and connections are used much less frequently. Thanks to this ‡exibility, a lot of e¤ort has been invested into computational methods, mesh generation and visualisation software adapted to triangular meshes (e.g. George, 1991). Compared to structured meshes, their drawback is that indexing is replaced by connectivity lists and lookup tables, and mesh generation is typically a lengthy, user-assisted and somewhat arbitrary process. Moreover, with a …xed triangular mesh, boundary …tting is only possible if the location of the boundary is knowna priori. In the case of mildly sloped topography, water level variations may result in sawtooth boundaries unless the mesh is adapted to the solution dynamically. Local re…nement or coarsening of the initial mesh during the solution is complicated: it is either realised through subdivision followed by edge ‡ipping or by remeshing the whole domain according to a new density distribution.
Quadtree meshes are a blend of structured grids and unstructured meshes (Figure 4d). Their hierarchical structure is based on the principle of recursive decomposition of the plane and is uniquely de…ned by a linked tree structure. In spite of such a clear structure, the number of neighbours is arbitrary and a quadtree mesh can be considered an unstructured one in which the elements may have any number of sides. Originally, quadtrees and related data structures were developed for use in image processing and spatial information systems. Samet (1984) surveys several types of hierarchical data structures and their algorithms; in this work the scope is restricted to Cartesian quadtrees decomposed into equal parts on each level (which is termed regular decomposition). Mesh generation is automatic and robust, furthermore the quadtree structure is inherently designed to provide static or dynamic variability of the spatial resolution. Arbitrary boundaries are resolved by extra re…nement along those curves or more sophisticated techniques such as cut-cells and immersed boundaries, as will be outlined in §5.2. True anisotropy is di¢ cult to implement e¢ ciently on Cartesian quadtree meshes and the variation of cell sizes cannot be as smooth as on triangular or curvilinear meshes.
A generally applicable mesh should also be able to resolve an intricate reed cover, such as the one that characterises the southern part of Lake Neusiedl (see
§6.1 for the description of the lake). The conveyance of open water and reed di¤er so dramatically that the ‡ow near the reed boundary is mainly tangential, hence they should be preferably discretised like closed boundaries. Figure 5 presents the meshes generated with each of the four types. Once the conceptual model that describes zone boundary polygons, zone attributes and cell resolution is ready, generation is automatic in all cases, except for the curvilinear grid. The latter
(a) Regular (b) Curvilinear
(c) Triangular (d) Quadtree
Figure 5. The neighbourhood of the bay of Fert½orákos in the southern part of Lake Neusiedl discretised using commonly used mesh types. The shading is according to the coverage type (open water or reed).
required a lot of manual editing to accommodate zone boundaries with coordinate lines and even so, some triangles remain in the grid. The unstructured meshes are successfully concentrated to the semi-enclosed bay of Fert½orákos and to the principal strait that connects the southern and northern parts of the lake.
The merits of hierarchical meshes –‡exibility, automatic generation, easy ad-aptivity –justify their use as the basis of a lake ‡ow model. To ensure robustness and ease of implementation, a Cartesian quadtree mesh is adopted in this thesis.