Oct 8, 2006 - A água do rio, correndo, cinge as múltiplas ilhas, e nesta margem direita [...] há braços lıquidos que entram pela terra e reflectem o céu.
UNIVERSIDADE DE LISBOA INSTITUTO SUPERIOR TÉCNICO
Two Dimensional Unstructured Mesh Generation applied to Shallow Water Models Andrea Mazzolari Supervisor: Doctor António Alexandre Trigo Teixeira Thesis approved in public session to obtain the PhD Degree in Civil Engineering Jury final classification: Pass with Merit Jury Chairperson: Chairman of the IST Scientific Board Members of the Committee: Doctor Piero Ruol – Full Professor, Universitá Degli Studi di Pádova, Itália; Doctor Fernando Francisco Machado Veloso Gomes, Professor Catedrático da Facultade de Engenharia, da Universidade de Porto; Doctor José Paulo Baptista Moitinho de Almeida, Professor Associado (com Agregação) do Instituto Superior Técnico, da Universidade de Lisboa; Doctor António Alexandre Trigo Teixeira, Professor Associado do Instituto Superior Técnico, da Universidade de Lisboa; Doctor Carlos Daniel Borges Coelho, Professor Auxiliar da Universidade de Aveiro; Doctor António Alberto Pires Silva, Professor Auxiliar do Instituto Superior Técnico, da Universidade de Lisboa; Doctor Maria Amélia Vieira da Costa Araújo, Investigadora Auxiliar do Instituto Superior Técnico, da Universidade de Lisboa.
Funding Institution Fundação para a Ciência e Tecnologia
2013
...uma n´evoa luminosa que o sol faz reverberar por dentro como um resplendor. A ´agua do rio, correndo, cinge as m´ ultiplas ilhas, e nesta margem direita [...] h´a bra¸cos l´ıquidos que entram pela terra e reflectem o c´eu Viagem a Portugal
Jos´e Saramago
Abstract
This research is focused on the generation of two dimensional unstructured meshes for finite element shallow water models. The work is structured over two main lines of investigation: in the first one an existing advancing front algorithm is re-activated and improved, in terms of robustness, memory efficiency, post-processing options and capacity to deal with complex geometries. An alternative implementation of the advancing front method, which introduces the concept of pseudo-island, is proposed. A graphical user friendly interface is provided for a simple and intuitive use of the algorithm itself. The second line of research consists in the definition and application of a multi-criteria strategy for calculating the node spacing function driving the domain discretization: a series of a-priori and/or a-posteriori functions, each representing a criterion of discretization, are defined and superimposed in order to calculate the target node spacing function of the prospective mesh. A-priori criteria may express spatial restrictions derived from the elevation model, geometrical features, known physical boundary conditions and internal domain forcings. A-posteriori criteria are instead related to the flow velocity or to the error estimation of an existing model solution. This strategy has been tested for the generation of meshes over idealized and real domains, for isotropic as well as anisotropic conditions, following an increasing degree of complexity. Mesh quality checks and quality enhancement procedures have been performed. Meshes resulting from the discretization of ocean, coastal and estuarine regions have been successfully used inside state-of-the art shallow water models for the resolution of hydrodynamic problems of different scales. Comparative tests between a-priori and a-posteriori criteria have been included; particular attention has been paid to the relationship between meshing criteria and mesh behaviour in presence of a wet and dry algorithm. The performance of the meshes has been evaluated both in terms of computational stability and model validation. The work done is believed to represent
v
a practical tool and a reference for the modeller in the design of unstructured meshes over shallow water real domains.
Key words: finite elements, shallow waters, unstructured mesh, node spacing function, advancing front, pseudo-island, meshing criterion, wet and dry, digital elevation model, error estimator.
vi
Acknowledgements
First of all, I would like to acknowledge Professor Antonio Trigo Teixeira for his supervision and for the attention dedicated to my work. His experience and sharing of information represent a valuable contribution for the completion of this work. I would also like to thank Dr. Amelia Ara´ ujo for interesting discussions and precious collaboration. A number of institutions, colleagues and professionals contributed to the progress of my PhD: Crystal Fulcher and Jason Fleming, for providing the ADCIRC parallel version; Jonathan Shewchuk for making available TRIANGLE; Ata Bilgili and Keston Smith for the BatTri interface; professor Lee Chi King for sending me the AAMGM source code; Cristina Monteiro for introducing me to the LIDAR dataset of the Portuguese coastline; Professor Borges Dinis for allowing the use the ABACO resources. My sincere appreciation goes to the Instituto Geogr´afico Portuguˆes (Portuguese Geographical Institute), for providing me, as a first user for research applications, the Lidar measurements of the Portuguese coastline; and to the Funda¸c˜ ao para a Ciˆencia e a Tecnologia (Science and Technology Foundation), for awarding me the doctoral grant SFRH / BD / 61161 / 2009. Many thanks to the secretary staff of the CEHIDRO and to the Apoio T´ecnico of the IST Civil Department, for the assistance in solving all kinds of practical issues a PhD student may have. To the doctoral and postdoctoral colleagues of the CEHIDRO group, in particular Maria Jo˜ ao, Nuno M., Isabel, Nuno, Nelson, Ana M., Ana, Mariana, for the help and the useful breaks had during long afternoons of work. At the end of this course, my thoughts go to my family in Italy and to Darya, for their support during difficult times; and to all the people I met in the nice country of Portugal, with whom I shared beautiful moments.
Contents List of Figures
xiii
List of Tables
xix
List of Symbols
xxi
List of Acronyms and Abbreviations
xxv
1 INTRODUCTION
1
1.1
About mesh generation and shallow water models . . . . . . . . . . . . .
1
1.2
Aims . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2
1.3
Outline of the thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
3
2 MATHEMATICAL BACKGROUND
5
2.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
2.2
Shallow water equations . . . . . . . . . . . . . . . . . . . . . . . . . . .
5
2.3
The finite element method . . . . . . . . . . . . . . . . . . . . . . . . . .
7
2.3.1
The Galerkin method . . . . . . . . . . . . . . . . . . . . . . . .
7
2.3.2
Finite elements and basis functions . . . . . . . . . . . . . . . . .
9
2.3.3
Linear basis functions for 1D finite elements . . . . . . . . . . . .
10
2.3.4
Linear basis functions for 2D finite elements . . . . . . . . . . . .
11
2.3.5
Quadratic basis functions for 2D finite elements . . . . . . . . . .
12
2.3.6
Other types of elements . . . . . . . . . . . . . . . . . . . . . . .
13
3 OVERVIEW OF THE MESH GENERATION METHODS
15
3.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
15
3.2
Classification scheme of Ho-Le (adapted) . . . . . . . . . . . . . . . . . .
15
ix
CONTENTS
3.3
3.4
3.5
3.2.1
Mesh topology first . . . . . . . . . . . . . . . . . . . . . . . . . .
16
3.2.2
Nodes first . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
16
3.2.3
Adapted mesh template . . . . . . . . . . . . . . . . . . . . . . .
17
3.2.4
Nodes and elements generated simultaneously . . . . . . . . . . .
18
Unstructured triangular mesh generation methods . . . . . . . . . . . .
19
3.3.1
The Delaunay triangulation . . . . . . . . . . . . . . . . . . . . .
19
3.3.2
The advancing front method . . . . . . . . . . . . . . . . . . . .
22
3.3.3
The quadtree method . . . . . . . . . . . . . . . . . . . . . . . .
26
3.3.4
The hybrid method . . . . . . . . . . . . . . . . . . . . . . . . . .
26
Point creation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
27
3.4.1
Independent node generation technique . . . . . . . . . . . . . .
27
3.4.2
Grid superimposition and successive subdivision . . . . . . . . .
27
3.4.3
Point creation driven by the boundary point distribution . . . .
27
3.4.4
Point creation through a node spacing function . . . . . . . . . .
28
Mesh quality evaluation and optimization . . . . . . . . . . . . . . . . .
32
4 MESHING CRITERIA FOR SHALLOW WATER MODELS
39
4.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
39
4.2
Tidal modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
40
4.3
Resolution requirements for steep bathymetric gradients . . . . . . . . .
43
4.4
Hurricane and storm surge modelling . . . . . . . . . . . . . . . . . . . .
44
4.5
Domains with complex coastline patterns . . . . . . . . . . . . . . . . .
46
4.6
In confined waters and low-lying areas . . . . . . . . . . . . . . . . . . .
50
4.7
A-posteriori meshing techniques . . . . . . . . . . . . . . . . . . . . . .
51
4.7.1
Error estimators . . . . . . . . . . . . . . . . . . . . . . . . . . .
51
4.7.2
A numerical constraint: the Courant condition . . . . . . . . . .
52
Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
53
4.8
5 UNSTRUCTURED MESH GENERATORS
55
5.1
Introducton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
55
5.2
MESHGR . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
56
5.2.1
Improvements introduced in MESHGR . . . . . . . . . . . . . . .
56
5.2.2
Structure of MESHGR . . . . . . . . . . . . . . . . . . . . . . . .
58
Other mesh generators . . . . . . . . . . . . . . . . . . . . . . . . . . . .
67
5.3
x
CONTENTS
5.3.1
TRIANGLE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
69
5.3.2
Application of TRIANGLE to a real case study . . . . . . . . . .
70
5.3.3
The SMS mesh generator . . . . . . . . . . . . . . . . . . . . . .
76
5.3.4
Application of the SMS paving algorithm to an oceanic domain .
78
5.3.5
AAMGM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
79
5.3.6
TRIGRID . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
81
5.3.7
Comparative analysis between mesh generators . . . . . . . . . .
82
6 THE MULTI-CRITERIA MESHING METHOD. APPLICATION TO IDEALIZED GEOMETRIES 85 6.1
Introducton . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
85
6.2
Test 1: rectangular domain with a steep bottom . . . . . . . . . . . . .
86
6.3
Test 2: squared domain with 2 islands . . . . . . . . . . . . . . . . . . .
96
6.4
Introduction to anisotropic meshing . . . . . . . . . . . . . . . . . . . . 101
6.5
Extension of the multi-criteria method to anisotropic meshes . . . . . . 104
6.6
Test 3: anisotropic meshing . . . . . . . . . . . . . . . . . . . . . . . . . 106
6.7
Test 4: meshing criteria for the wet and dry . . . . . . . . . . . . . . . . 108
6.8
6.7.1
Domain description and mesh generation . . . . . . . . . . . . . 112
6.7.2
Tide modelling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 113
6.7.3
Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115
7 APPLICATION OF THE MULTI-CRITERIA METHOD TO A REAL DOMAIN 119 7.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
7.2
Domain geometry and bathymetry . . . . . . . . . . . . . . . . . . . . . 120
7.3
Meshing criteria for a multi-scale domain . . . . . . . . . . . . . . . . . 121
7.4
7.3.1
Ocean, shelf break and continental shelf . . . . . . . . . . . . . . 122
7.3.2
At the coast and inside the estuary . . . . . . . . . . . . . . . . . 122
7.3.3
Numerical constraints . . . . . . . . . . . . . . . . . . . . . . . . 123
7.3.4
Final node spacing function . . . . . . . . . . . . . . . . . . . . . 124
Model calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125 7.4.1
Free surface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
7.4.2
Depth averaged velocities . . . . . . . . . . . . . . . . . . . . . . 128
xi
CONTENTS
7.4.3 7.5
7.6
Model validation: flooded area extension . . . . . . . . . . . . . . 128
Evaluation of an a-posteriori meshing criterion . . . . . . . . . . . . . . 133 7.5.1
The LTEA+CD meshing method . . . . . . . . . . . . . . . . . . 133
7.5.2
A-priori coarse estuary discretization . . . . . . . . . . . . . . . 135
7.5.3
LTEA+CD mesh convergence study . . . . . . . . . . . . . . . . 137
7.5.4
Model application . . . . . . . . . . . . . . . . . . . . . . . . . . 142
Assessment of alternative meshing criteria for a multi-scale ocean to estuary domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145
7.7
7.6.1
Mesh refinement series . . . . . . . . . . . . . . . . . . . . . . . . 146
7.6.2
Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146
Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149
8 PSEUDO-ISLANDS
153
8.1
Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 153
8.2
Test cases for idealized geometries . . . . . . . . . . . . . . . . . . . . . 154
8.3
8.4
8.5
8.2.1
Test case 1: NSF with one spike . . . . . . . . . . . . . . . . . . 154
8.2.2
Test case 2 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 159
Application of pseudo-islands to a LIDAR derived digital elevation model 160 8.3.1
Study sites . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162
8.3.2
Methodology for the mesh accuracy assessment . . . . . . . . . . 163
8.3.3
Hydrodynamic model description . . . . . . . . . . . . . . . . . . 164
Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166 8.4.1
C´ avado estuary . . . . . . . . . . . . . . . . . . . . . . . . . . . . 166
8.4.2
Ria Formosa: sites 1 and 2 . . . . . . . . . . . . . . . . . . . . . 168
8.4.3
Tidal simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . 168
Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171
9 CONCLUSIONS
173
9.1
Concluding remarks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 173
9.2
Further work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175
Bibliography
177
xii
List of Figures 2.1 2.2
Reference system adopted for the SWEs variables. . . . . . . . . . . . . e ∪ Γ. . . . . . . . . . . . . . . . . . . . . . . . . Arbitrary domain Ω = Ω
6 7
2.3
Conforming and not conforming subdivision of a 2D domain into elements 10
2.4
Elements and nodes for a 1D domain. . . . . . . . . . . . . . . . . . . .
11
2.5
Basis functions for 1D domain . . . . . . . . . . . . . . . . . . . . . . . .
11
2.6
Triangle e with nodes x¯1 , x¯2 and x¯3 . . . . . . . . . . . . . . . . . . . . .
12
2.7
Trianglular element e with 6 nodes. . . . . . . . . . . . . . . . . . . . . .
12
3.1
Classification scheme for mesh generation methods. Adapted from Ho-Le (1988). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
16
3.2
Topology decomposition of a 2D domain . . . . . . . . . . . . . . . . . .
17
3.3
A simple example of an adapted mesh template with a grid-based approach 18
3.4
Example of a mapped element approach . . . . . . . . . . . . . . . . . .
19
3.5
Mesh obtained by an algebraic system approach. . . . . . . . . . . . . .
21
3.6
Delaunay triangulation and Vorono¨ı regions for a given set of points . .
22
3.7
Steps for updating the active front in the AFM . . . . . . . . . . . . . .
24
3.8
Successive stages in the implementation of an AFM technique . . . . . .
25
3.9
Common issues during mesh generation performed by an AFM algorithm 25
3.10 Example of a quadtree recursive partition and related tree structure . .
26
3.11 Point creation driven by the boundary point distribution for a river bed domain split into sub-regions . . . . . . . . . . . . . . . . . . . . . . . . 3.12 Point distribution controlled by positioning point and line sources
28
. . .
30
3.13 Quality index qi for different triangle shapes. . . . . . . . . . . . . . . .
33
3.14 Notation used for the definition of the triangle quality indexes . . . . . .
34
3.15 Example of Laplacian smoothing . . . . . . . . . . . . . . . . . . . . . .
36
xiii
LIST OF FIGURES
3.16 Diagonal swapping for two adjacent triangles . . . . . . . . . . . . . . .
37
3.17 Diagonal swapping for obtaining triangles satisfying the Delaunay condition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1
37
Unstructured graded finite element meshes of the Western North Atlantic, Gulf of Mexico and Caribbean Sea . . . . . . . . . . . . . . . . .
42
4.2
Anisotropic mesh for an idealized bathymetric profile . . . . . . . . . . .
45
4.3
Piece-wise representation of the inverted barometer forcing function . .
46
4.4
Discretization of a hypothetical sinusoidal coastline . . . . . . . . . . . .
47
4.5
Discretization of a highly irregular and complex region in the Artic Sea
48
4.6
Example of a mesh covering the region of Lake Maracaibo . . . . . . . . ´ Obidos lagoon: generation of different meshes for different spatial scales
48
and merging algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . .
49
4.7 4.8 4.9
Southern Louisiana and estuary of the Mississippi river covered by a graded mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
51
Mesh portion covering the city of New Orleans . . . . . . . . . . . . . .
52
4.10 Gulf of Mexico discretized (a) with the LTEA approach and (b) with the wavelength to mesh size ratio criterion . . . . . . . . . . . . . . . . . . .
53
5.1
Flow chart of the GUI interfaced with MESHGR . . . . . . . . . . . . .
58
5.2
GUI provided to MESHGR . . . . . . . . . . . . . . . . . . . . . . . . .
59
5.3
Example of a non-simply connected domain formed by 3 adjacent regions with an internal opening . . . . . . . . . . . . . . . . . . . . . . . . . . .
5.4
Visual interpretation of the effect of the transformation T(x,y) between the physical and the control space. . . . . . . . . . . . . . . . . . . . . .
5.5
61 64
Identification of a prospective node N and of its node neighbours on the active front. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
66
5.6
Determination of the angles αi between AB and the nodes Qi ∈ Γ. . . .
66
5.7
Effect of Laplacian smoothing. . . . . . . . . . . . . . . . . . . . . . . .
67
5.8
Combined action of diagonal swapping and Laplacian smoothing . . . .
68
5.9
TRIANGLE: successive refinement procedure . . . . . . . . . . . . . . .
71
5.10 TRIANGLE: new series of refinements based on the wavelength to mesh size criterion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
xiv
72
LIST OF FIGURES
5.11 Computational domain, extending for 400 km along the west Iberian Coast and 615 km across the Atlantic Ocean. . . . . . . . . . . . . . . .
73
5.12 Example of “first cut” mesh for the domain of Figure 5.11. . . . . . . .
73
5.13 Unstructured triangular mesh with 41367 nodes and 80590 elements, used in Ara´ ujo et al. (2011a). . . . . . . . . . . . . . . . . . . . . . . . .
74
5.14 Comparison between the bathymetric isolines of the mesh of Figure 5.13 and of the DEM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
75
5.15 Time step requirement for a Cr = 1 for the mesh of Figure 5.13 . . . . .
76
5.16 Subdivision of a river domain into subregions . . . . . . . . . . . . . . .
77
5.17 Extension of the new ocean domain (outer polygon), compared with the mesh used in the work of Ara´ ujo et al. (2011a). . . . . . . . . . . . . . .
78
5.18 2300 Km x 1350 Km unstructured mesh generated with the SMS scalar paving method for the Atlantic domain of Figure 5.17. . . . . . . . . . .
79
5.19 Hierarchical generation principle in AAMGM . . . . . . . . . . . . . . .
80
5.20 TRIGRID clusters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
82
6.1
Test 1. Domain extension and bathymetric isolines. . . . . . . . . . . . .
87
6.2
Test 1. Background mesh. . . . . . . . . . . . . . . . . . . . . . . . . . .
87
6.3
Test 1. Mt11 generated by MESHGR before post-processing . . . . . . .
88
6.4
Test 1. Histograms of the shape quality index qi for Mt11 . . . . . . . .
88
6.5
Test 1. Mt11 before post-processing, evaluated against the criteria of maximum area change, maximum slope and node valence higher than 8
6.6
89
Test 1. Mt11 after diagonal swap and Laplacian smoothing, evaluated against the criteria of area change, maximum slope and node valence higher than 8. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
90
6.7
Test 1. Mt12 before applying post-processing operations.
91
6.8
Test 1. Mt12 quality before diagonal swapping and Laplacian smoothing,
. . . . . . . .
evaluated against the criteria of area change, maximum slope and node valence. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.9
91
Test 1. Mt12 quality after diagonal swapping and Laplacian smoothing, evaluated against the criteria of area change, maximum slope and node valence higher than 8. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
92
6.10 Test 1. Multi-criteria method . . . . . . . . . . . . . . . . . . . . . . . .
93
xv
LIST OF FIGURES
6.11 Test 1.
Mt13 before diagonal swap and Laplacian smoothing post-
processing.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
94
6.12 Test 1. Mt13 after diagonal swap and Laplacian smoothing post-processing. 94 6.13 Test 1. Domain zoning for Mt13 showing the spatial influence of each meshing criterion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
95
6.14 Test 1. Mt13. Histograms of the quality index qi before and after postprocessing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
96
6.15 Test 1. Mt13 quality before diagonal swap and Laplacian smoothing, evaluated against the criteria of area change and node valence higher than 8 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
97
6.16 Test 1. Mt13 quality after diagonal swap and Laplacian smoothing, evaluated against the criteria of area change and node valence higher than 8. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
97
6.17 Test 2. Layout of the domain and bathymetric contours. . . . . . . . . .
98
6.18 Test 2. Background mesh. . . . . . . . . . . . . . . . . . . . . . . . . . .
98
6.19 Test 2. Magnitude of the NSF ρ21 , determined as the distance from the nearest island. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
99
6.20 Test 2. Mt21 obtained by applying ρ21 , after post-processing operations
99
6.21 Test 2. NSF ρ23 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 6.22 Test 2. Mt23, obtained by applying ρ23 , after post-processing operations. 101 6.23 Test 2. Domain zoning for mesh Mt23 . . . . . . . . . . . . . . . . . . . 102 6.24 Test 2. Histograms of the quality index qi for Mt23. . . . . . . . . . . . 102 6.25 Test 2. Comparison between the Mt23 bathymetry and the scatter set bathymetry of Figure 6.17 . . . . . . . . . . . . . . . . . . . . . . . . . . 103 6.26 Anisotropic mesh . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105 6.27 Test 3. Domain extension and bathymetric contours. . . . . . . . . . . . 106 6.28 Test 3. Background mesh. . . . . . . . . . . . . . . . . . . . . . . . . . . 107 6.29 Test 3. Bathymetric gradient vectors defined over the bathymetric isolines.107 6.30 Test 3. Anisotropic mesh Mt31s . . . . . . . . . . . . . . . . . . . . . . . 109 6.31 Test 3. Anisotropic mesh Mt31c
. . . . . . . . . . . . . . . . . . . . . . 110
6.32 Test 4. Bathymetric depths for the wet and dry test case. . . . . . . . . 113 6.33 Test 4. Wet and dry test. Mesh Mt4AN with location of the reference station. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114
xvi
LIST OF FIGURES
6.34 Test 4. Water surface elevation at the reference station . . . . . . . . . . 115 6.35 Test 4. Profile of the wetted front for Mt4AN. The close-up views (a) and (b) contrast the position of the wetted front with the nearest bathymetric isoline. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116 6.36 Test 4. Comparison between the position of the wetted fronts of the anisotropic and isotropic meshes . . . . . . . . . . . . . . . . . . . . . . 117 7.1
Extension of the ocean to Lima estuary computational domain . . . . . 121
7.2
Unstructured graded mesh derived from the node spacing function ∆xf
7.3
Distribution of the element quality index qi for the post-processed mesh 125
7.4
Comparison of the observed and computed water elevations at (a) site
124
1, (b) site 3 and (c) site 4. . . . . . . . . . . . . . . . . . . . . . . . . . . 129 7.5
Comparison between the magnitude of the observed and computed depth averaged velocities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 130
7.6
Comparison between the locations of the wetted fronts taken from aerial photos and from the ADCIRC wet and dry simulation, for the Lima middle estuary. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 131
7.7
Water surface extension in the Lima estuary calculated along the day 8th October 2006 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132
7.8
Mesh M70 for the ocean and estuarine domain . . . . . . . . . . . . . . 136
7.9
Local truncation error calculated for meshes M70, M70 1 and M70 2 . . 138
7.10 Areas of refinement for M70 1 with respect to M70 . . . . . . . . . . . . 140 7.11 Areas of refinement for M70 2 with respect to M70 1 . . . . . . . . . . . 141 7.12 Comparison between measured and computed depth averaged velocities for meshes M70 and M70 D . . . . . . . . . . . . . . . . . . . . . . . . . 145 7.13 Evaluation of the area gap between the wetted boundary extracted from aerial photos and the computed wetted boundary of mesh M1 LT, for September 10th 2010, at 10.30am. . . . . . . . . . . . . . . . . . . . . . . 150 8.1
Test 1. NSF defined over a squared domain . . . . . . . . . . . . . . . . 155
8.2
Test 1. Mesh T1 A, produced with the circulating option and mesh T1 B, produced with the front searching option of the standard AFM implemented in MESHGR, for the NSF of Figure 8.1. . . . . . . . . . . 155
8.3
Successive stages of the pseudo-island front searching algorithm . . . . . 157
xvii
LIST OF FIGURES
8.4
Test 1. Pseudo-island front searching algorithm applied to the NSF of Figure 8.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158
8.5
Test 1. qi distribution into classes for T1 C and T1 D . . . . . . . . . . 159
8.6
Test 2. NSF defined over a rectangular domain with one linear and one curvilinear high gradient areas. . . . . . . . . . . . . . . . . . . . . . . . 160
8.7
Test 2. Meshes obtained for the NSF of Figure 8.6 . . . . . . . . . . . . 161
8.8
Test 2. qi distribution into classes for T2 B (left side) and T2 C (right side) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161
8.9
C´avado estuary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162
8.10 View of the Ria Formosa, with the delimitations of the selected study sites. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163 8.11 Extension of the Ria Formosa domain used for a wet and dry tide model application. The location of site 1 is evidenced as well. . . . . . . . . . . 165 8.12 C´avado estuary bathymetry and test meshes . . . . . . . . . . . . . . . . 167 8.13 Comparison between meshes and the respective wetted areas for the study site 1 of the Ria Formosa during flood tide . . . . . . . . . . . . . 170
xviii
List of Tables 3.1
Advantages and disadvantages of the mesh generation methods. . . . . .
20
5.1
Some general properties of the tested mesh generators . . . . . . . . . .
83
6.1
Quality measures for the meshes of test 1, before and after post-processing. 95
6.2
Meshes parameters for the wet and dry modelling test . . . . . . . . . . 115
7.1
Manning values associated to the estuary surface cover types . . . . . . 126
7.2
Mesh properties for the LTEA+CD convergence study. . . . . . . . . . . 137
7.3
Performance indices for meshes M70, M70 1, M70 2 and M70 D, at sites 1, 3 and 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143
7.4
Properties of the meshes corresponding to different meshing criteria . . 146
7.5
Performance indices for meshes M1, M1 LT, M70 and M10, at sites 1, 3 and 4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 147
7.6
Assessment of the wet and dry area extension for alternative meshing criteria. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149
8.1
Mesh statistics for the discretization tests of the C´avado estuary. . . . . 167
8.2
Main properties of the meshes for the selected sites of the Ria Formosa. 169
8.3
Measured and computed tidal constituent values for the Faro-Olh˜ao gauge171
xix
List of Symbols ai , bi , ci
side lengths of a triangle
b
bottom depth (positive if above MSL)
c
proportionality coefficient
d
bottom depth (positive if below MSL)
ei
mesh element
f
1) known term of a PDE; 2) Coriolis parameter
g
module of the gravitational field
h
water column depth
i, j, k
indices
m
meters
n
Manning coefficient
~n
normal vector
p
fluid pressure
pa
barometric pressure at the free surface
qi
element equilaterality index
r
Earth radius
r1
in-circle radius of a triangle
r2
circum-circle radius of a triangle
s
seconds
si
mesh anisotropic coefficient
xxi
LIST OF SYMBOLS
t
time
u
1) component of the mesh element orientation vector along xi ; 2) exact solution of a PDE
u ˜
approximate solution of a PDE
v
component of the mesh element orientation vector along xj
v¯i
matrix eigenvector in the xi direction
w
polynomial function
xi
1) axis in the plane orthogonal to xj ; 2) mesh node in 1D
x¯i
mesh node in 2D or more
xj
axis in the plane orthogonal to xi
z
vertical axis in a Cartesian coordinate system
A
1) mesh element area; 2) stiffness matrix
Ai
area of the ith mesh element
Aij
element ij of the stiffness matrix
AB
segment of the advancing front
Cf
bottom drag coefficient
Cr
Courant number
Fi
Coriolis force in the xi direction
I
identity matrix
Lh
topographic length scale
M
metric tensor
Pi
point of a scatter set
Qi
1) point of a scatter set; 2) point belonging to the advancing front
R
wavelength to grid size ratio
Ri
ith subregion of a domain
T
wave period
T
transformation between the physical and the dual space
xxii
LIST OF SYMBOLS
Ui
depth averaged velocity along the xi direction
Uj
depth averaged velocity along the xj direction
V
flow velocity module
Vi
Vorono¨ı region
VN
sub-space of a functional space
V0
functional space
α
1) scalar coefficient; 2) Earth elasticity factor
αi
1) scalar coefficient; 2) unknown variable in a system of equations
α ¯
N-dimension column vector
α ˆ i , βˆi , γˆi
triangle interior angles
γi
in-circle to circum-circle ratio
γw
fluid specific weight
δij
Kronecker’s delta
η
free surface elevation (relative to MSL)
θ
degrees latitude
λ
1) wavelength; 2) degrees longitude
λi
1) basis function; 2) ith eigenvalue of a matrix
µh
horizontal eddy viscosity coefficient, under the Boussinesq hypothesis
µi
target side deviation parameter
νT
depth average horizontal eddy viscosity
ρ
target element side
ρi
target element side for the ith meshing criterion
ρw
fluid density
ςi
edge length quality index
τi
edge ratio
τsi
surface shear stress along the xi direction
τb
bottom shear stress module
xxiii
LIST OF SYMBOLS
τbi
bottom shear stress along the xi direction
+ τˆM E
local truncation error
χ
Newtonian equilibrium tidal potential
Γ
1) domain boundary(ies); 2) advancing front for t >0
Γ0
initial advancing front
∆t
time step
∆x
mesh element size
Θ
blending function ranging from 0 to 1
Φi
basis function in the Galerkin method
Ω
computational domain
e Ω
computational domain without boundary(ies)
B
background mesh
R
set of the real numbers
∇h
module of the bathymetric gradient
~ ∇h
bathymetric gradient
xxiv
List of Acronyms and Abbreviations 1D
One Dimensional
2D
Two Dimensional
3D
Three Dimensional
AAMGM
Automatic Adaptive Mesh Generator Method
ADCIRC
ADvanced CIRCulation model
ADCP
Acoustic Doppler Current meter Profiler
AFM
Advancing Front Method
BVP
Boundary Value Problem
CDT
Conforming Delaunay Triangulation
DEM
Digital Elevation Model
DT
Delaunay Triangulation
FEM
Finite Element Method
GUI
Graphical User Interface
LIDAR
LIght Detection And Ranging
LTEA
Local Truncation Error Analysis
LTEA-CD
Local Truncation Error Analysis with Complex Derivatives
LTEA+CD
Local Truncation Error Analysis with Complex Derivatives, including bottom stresses and the Coriolis force
MAE
Mean Average Error
xxv
LIST OF ACRONYMS AND ABBREVIATIONS
MSL
Mean Sea Level
NBCS
Number of Boundary Curves/Segments of the MESHGR domain
NELEM
Number of ELEMents of the MESHGR background mesh
NBV
Number of Boundary Vertices of the MESHGR domain
NPOIG
Number of POInts of the MESHGR backGround mesh
NREG
Number of REGions of the MESHGR domain
NSE
Nash-Sutcliffe Efficiency coefficient
NSF
Node Spacing Function
PDE
Partial Differential Equation
PSLG
Planar Straight Line Graph
RMSE
Root Mean Square Error
SMS
Surface water Modelling System
SWE
Shallow Water Equation
TIN
Triangulated Irregular Network
xxvi
1
INTRODUCTION 1.1
About mesh generation and shallow water models
Mesh generation is the area of the mathematical sciences which provides solutions for the discretization of a computational domain. In the engineering practice, complex real case problems are normally described in terms of partial differential equations, for which no analytical resolving expression is known. The finite difference, the finite element and the finite volume methods are nowadays well established techniques for replacing the continuous form of the differential equations into a system of discrete algebraic equations. The way this discretization is accomplished produces different mesh types: the finite difference method requires a uniform spacing between points, giving origin to a structured mesh (or grid). The finite element and finite volume methods allow a varying spacing, and are associated to what is known as an unstructured mesh. For clarity the term “grid” is used within the text when referring to a structured mesh, and “mesh” just to unstructured meshes, although in some literature, especially of American origin, these two terms are used interchangeably. In the last decades major advances in modelling research and computational capabilities provided robust and efficient mesh generation codes, with an increasing degree of automation and powerful computer-aided design tools, both freeware and commercial. On the other side, mesh generation may still be considered an “art”, in the sense that there does not exist something like a best meshing technique or ideal mesh. On the contrary, every modelling application presents peculiar necessities and constraints. The present research is focused on the problem of mesh generation related to shallow
1
1. INTRODUCTION
water models. Unstructured meshes have received attention from developers and practitioners for their flexibility in handling complex domains, as estuaries, deltas, coastal zones and the open ocean. They are often used in storm surges and inundation studies. A mesh should, in general, meet these conditions: a) discretize the domain at a level of refinement which does not influence the model output. b) avoid computational overload. The two conditions are in conflict, and the task of the modeller is to find a good balance between them. Property a) is fulfilled when the mesh sizing adjusts to the physical gradients of the variables involved and reproduces the geometrical domain singularities. Also the mesh quality, evaluated in terms of geometrical parameters of the elements, should be considered. Property b), on the contrary, is achieved when there is not over-refinement. The overall accuracy of a model depends, among other factors, on a), while its computational performance on b). For what said, it is clear that the mesh generation process is a non-trivial step in the build-up of a computational model. However, often in the modelling practice, in the published literature and inside the model user manuals, the problem of spatial discretization is or poorly treated or completely overlooked, leaving uncertainty about the level of accuracy of the results.
1.2
Aims
Despite the extensive amount of work published, the process of mesh generation for the vast majority of shallow water modelling studies is not described: these studies normally lack a discussion about the adopted meshing approach and its influence on the model results. Mesh quality analysis are rare too, while common experience suggests that the generated meshes present localized low quality features, in terms of element shape and connectivity, poorly resolved areas, small interior angles, which need to be corrected by a painstaking and error prone manual editing. The present research tries to give an insight into these issues by setting the following aims:
2
1.3 Outline of the thesis
• write a catalogue of the common meshing criteria, which define the node placement inside the mesh, applied to shallow water modelling studies, and their respective areas of validity; • define a comprehensive meshing strategy, based on the multi-criteria concept, according to the information acquired in the previous point;
• compare the effectiveness and efficiency of a-priori meshing criteria with established a-posteriori criteria;
• improve and develop an existing mesh generation algorithm (MESHGR) and its
relative graphical user interface, for the creation of unstructured meshes. The algorithm, based on the advancing front method, should be able to perform the spatial discretization of real multiply-connected domains according to a pre-defined node spacing function, to the desired level of refinement; to apply isotropic and anisotropic conditions to the element shape; to control the mesh topology by quality enhancement options; and to export the final mesh to available finite element models in a ready-to-use format.
1.3
Outline of the thesis
The thesis is structured into three parts: the first one, after this Introduction (chapter 1), defines the background over which the research develops: the shallow water equations and the generation of finite elements (chapter 2); a comprehensive classification of the mesh generation methods (chapter 3) and a review of the meshing criteria encountered in the literature (chapter 4) follow. The second part (chapter 5) describes the development of the MESHGR algorithm performed during this research and proposes a comparative analysis between MESHGR and other tested mesh generators. The third part presents a series of applications to test cases with an increasing order of complexity: idealized geometries are resolved in chapter 6, while the discretization of real domains for tide and storm surge model applications are included in chapter 7. Chapter 8 introduces an alternative version of the advancing front meshing method. The final chapter gives an evaluation of the work done and proposes further lines of research.
3
2
MATHEMATICAL BACKGROUND 2.1
Introduction
This chapter gives an outline of the shallow water equations and of the finite element method, from which the necessity of designing unstructured meshes in coastal and ocean models originates. The essential information is given here, as these topics are related but do not strictly represent the subject of the present work. The reader looking for a more exhaustive presentation is asked to follow the given references.
2.2
Shallow water equations
The shallow water equations (SWEs) are a set of partial differential equations, obtained, under the assumption that the vertical component of the fluid motion is negligible if compared with the horizontal components, by integrating over the vertical the mass and momentum conservation laws. The SWEs are applicable to a variety of hydrodynamic problems, wherever the horizontal length scale is much greater than the vertical length scale. In this paragraph the conservative form of the SWEs is reported, as expressed in Trigo-Teixeira (1994a), to which the reader is addressed for their full mathematical derivation. Under the mentioned condition, and considering the reference system of Figure 2.1, the vertical momentum equation is reduced to the hydrostatic pressure distribution:
5
2. MATHEMATICAL BACKGROUND
Figure 2.1: Reference system adopted for the SWEs variables: pa is the atmospheric pressure, MSL is the mean sea level, b(xi , xj , t) the bottom elevation, defined negative if below MSL, η(xi , xj , t) the free surface elevation and h the water column depth.
∂p = −gρw ∂z
(2.1)
where p is the pressure, z the vertical direction, g the gravity acceleration and ρw the fluid density. The SWEs can be written as: ∂ 1 ∂(Ui h) + ρw [Uj Ui h + δij g(h2 − d2 )] = ρw Fi + ∂t ∂xj 2 ∂d ∂pa ∂ ∂Ui ∂Uj +γw (h − d) −h + [hµh ( + )] + τsi − τbi ∂xi ∂xi ∂xj ∂xj ∂xi ρw
(2.2)
where the indices i and j correspond to the planar directions xi and xj , t is the time, Ui and Uj are the depth averaged velocities oriented along directions xi and xj , h is the water depth, δij is the Kronecker delta, d = d(xi , xj , t) is the bottom elevation, Fi is the term accounting for the Coriolis effect, γw = gρw is the fluid specific weight, pa is the barometric pressure at the free surface, µh is the horizontal eddy viscosity coefficient coming from the Boussinesq turbulent viscosity hypothesis for the Reynolds stresses, τsi is the surface shear stress and τbi the bottom stress in the xi direction. The SWEs, of which (2.2) is one of the many existing versions appearing in the literature, represent the fundamental equations adopted in numerical schemes for the simulation of tides, storm surges, ocean and estuarine currents.
6
2.3 The finite element method
2.3
e ∪ Γ. Figure 2.2: Arbitrary domain Ω = Ω
The finite element method
In most cases a partial differential equation (PDE), associated to a boundary value problem (BVP), does not have an analytical solution and needs to be solved in an approximate way. The Galerkin method is a general procedure for converting a continuous PDE into a discrete formulation of easier solution, and is reported here to illustrate the reason why a computational domain, associated to a PDE problem, has to be discretized into simple geometrical elements. The literature sources of this paragraph are the books of Cuvelier et al. (1986), Lapidus and Pinder (1982) and Rappaz et al. (1998).
2.3.1
The Galerkin method
A general PDE associated to its BVP can be written as: e in Ω
Lu = f u = uo
(2.3)
along Γ
where the operator L expresses the terms of a second order linear PDE in the unknown u = u(¯ x), being x ¯ ∈ Ω (Figure 2.2). f is the known term of the equation. u0 is the
7
2. MATHEMATICAL BACKGROUND
value that u takes along the domain boundary Γ. u is a hypothetical function satisfying (2.3), and is supposed to belong to the functional space V0 composed by all possible linear combinations w of basis functions Φi : w=
∞ X
α i Φi
(2.4)
i=1
where w ∈ V0 , and αi are unknown scalars. u can as well be written as a polynomial function, which expression, substituted into (2.3), allows to reformulate the initial problem in these terms: find u ∈ V0 such as: Z
a(u, Φi ) =
∀ Φi ∈ V 0
f Φi dΩ Ω
(2.5)
being a(u, Φi ) the term obtained by substitution of (2.4) into the term Lu of (2.3). (2.5) is called the variational (or weak) formulation of the initial problem. The Galerkin method supposes the existence of a sub-space VN ⊂ V0 composed by N basis functions (also called shape or interpolating functions) Φj in which an approximate formulation of (2.5) can be defined. In other terms, this is equivalent to say that there is a u ˜ ∈ VN such that: a(˜ u , Φj ) =
Z
f Φj dΩ
j = 1, . . . , N
(2.6)
Ω
where u ˜ takes the form: u ˜=
N X
α j Φj
(2.7)
j=1
(2.6) is a system of N equations in the unknowns αj , which is also possible to express in a matrix form: ¯ A¯ α=F
(2.8)
where the element Ai,j = Ai,j (Φi , Φj ). A is called stiffness matrix. α ¯ is the N dimension ¯ the N dimension column vector with Fj = column vector of the unknowns αj and F R Ω f Φj dΩ. The finite element method (FEM), which is a particular application of the
Galerkin method, provides a numerical technique for generating the basis functions Φj .
8
2.3 The finite element method
2.3.2
Finite elements and basis functions
The FEM applied to a Galerkin formulation of a BVP consists in constructing the basis functions Φj provided a spatial discretization of Ω into sub-regions, called elements. It is required that the discretization satisfies the following properties: i) there is a finite number of elements ek , k = 1, . . . , K; ii) given two elements ek and ej , with k 6= j, then either ek ∩ ej = ∅ or: ii.a) ek and ej have a common point for 1D domains; ii.b) ek and ej have a common edge or a common vertex for 2D domains; ii.c) ek and ej have a common face, a common edge or a common vertex for 3D domains; iii) the union
SK
k=1 ek
= Ω.
Figure 2.3(a) shows, for a 2D domain, a discretization which is meeting the requirements i), ii) and iii) and Figure 2.3(b) one which is not. The vertices x ¯j ∈ Ω of each element are called nodal points, or in a simpler way nodes. The set of elements in
which a region is subdivided, its nodes and the connections between nodes, constitute a mesh. The basis functions Φj (¯ xj ) : Ω → R, j = 1, . . . , N , should satisfy the following proper-
ties:
a) Φj is equal to 1 in x ¯j and vanishes at all other nodal points: Φj (¯ xi ) = δij
(2.9)
b) Φj has a prescribed behaviour (linear, quadratic, cubic,..) on Ω. The number of nodes per element depends on this behaviour, as demonstrated later; c) Φj is continuous on Ω. In order to illustrate how the FEM determines the basic functions, examples for simple discretizations for 1D and 2D domains are presented hereafter.
9
2. MATHEMATICAL BACKGROUND
(a)
(b)
Figure 2.3: (a) Conforming and (b) not conforming subdivision of a 2D domain into elements.
2.3.3
Linear basis functions for 1D finite elements
Considering a partition of a linear segment into N elements ej , with N + 1 nodal points xj (Figure 2.4), the wanted basis functions for this domain should be linear (condition b) and also meet conditions a) and c). Considering the following functions defined on the element ej = [xj−1 , xj ]: Φj−1 (x) = Φj (x) =
xj − x xj − xj−1
x − xj−1 xj − xj−1
(2.10) (2.11)
It is easy to demonstrate that, along ej , Φj−1 and Φj satisfy conditions a) and c) and correspond to acceptable basis function for the given discretization (Figure 2.5).
10
2.3 The finite element method
Figure 2.4: Elements and nodes for a 1D domain.
(a)
(b)
Figure 2.5: Basis functions (a) Φj−1 and (b) Φj defined for the interval ej = [xj−1 , xj ].
2.3.4
Linear basis functions for 2D finite elements
Considering a partition of a 2D domain Ω into triangular elements, the vertices x ¯j of each triangle are considered as nodes where to define a linear function (Figure 2.6). Therefore a candidate basis function for Φj (¯ xj ) assumes the form: λj (¯ xj ) = λ(x1 , x2 )j = α1 x1 + α2 x2 + α3
(2.12)
which is linear and continuous on the generic triangle e. (x1 , x2 )j are the coordinates of node j. In order to find the corresponding linear shape functions Φj (¯ xj ), j = 1, 2, 3, the coefficients α1 , α2 and α3 are determined according to condition b): λj (¯ xj ) = δij
i,j = 1, 2, 3
11
(2.13)
2. MATHEMATICAL BACKGROUND
Figure 2.6: Triangle e with nodes x¯1 , x¯2 and x¯3 .
Figure 2.7: Trianglular element e with 6 nodes.
which represents, for each function, a system of 3 equations in 3 unknowns. After calculating the solutions α1∗ , α2∗ and α3∗ , and substituting them into (2.12), it follows that Φ1 = λ1 , Φ2 = λ2 and Φ3 = λ3 . Again, u ˜ is obtained by substituting the found basis functions in (2.7) and resolving system (2.8).
2.3.5
Quadratic basis functions for 2D finite elements
In case quadratic basis functions are wanted, their corresponding general form is expressed as: Φj (x1 , x2 ) = α1 x21 + α2 x1 x2 + α3 x22 + α4 x1 + α5 x2 + α6
(2.14)
where the 6 unknown parameters αj need a system of 6 equations for each basis function, which implies the identification of 6 nodes per element (Figure 2.7). The shape functions Φj (¯ xj ) should be continuous, and, in order to meet (2.9), should be set equal to 1 in x ¯j and to 0 in x ¯i , with i 6= j. It is possible to demonstrate that: Φ1 = λ1 (2λ1 − 1)
12
2.3 The finite element method
Φ2 = λ2 (2λ2 − 1) Φ3 = λ3 (2λ3 − 1) Φ12 = 4λ1 λ2
Φ13 = 4λ1 λ3 Φ23 = 4λ2 λ3 where λ1 , λ2 and λ3 are the linear basis functions for the 2D finite elements.
2.3.6
Other types of elements
The FEM uses the same approach for finding the basis functions of other types of elements, both in 2D and 3D. The derivation of basis functions for these more complex cases can be found in Cuvelier et al. (1986).
13
3
OVERVIEW OF THE MESH GENERATION METHODS 3.1
Introduction
In the first part of this chapter a classification scheme of 2D mesh generation methods is presented. Several classifications or surveys are available in the literature (among many: Thacker, 1980; Ho-Le, 1988; Mavriplis, 1995a; Bossen and Heckbert, 1996; Owen, 1998; Frey and George, 2000). Due to the large amount of work done on the subject, any attempt of subdividing mesh generators into classes may result arbitrary and incomplete. The article of Ho-Le (1988), with integration of the more recent reviews of Lo (2002) and of Thompson et al. (1999), will be taken as a reference in this study, as it offers a systematic approach based on the temporal order of creation of nodes and elements. In the second part of the chapter the most successful mesh generation algorithms, namely the Delaunay triangulation (DT), the advancing front method (AFM) and the quadtree approach, are briefly presented; the third part is dedicated to the problem of node placement, while the last part deals with mesh quality measures and quality enhancement techniques.
3.2
Classification scheme of Ho-Le (adapted)
While most surveys are usually a simple listing of the most popular meshing algorithms, the taxonomic criterion here applied and adapted from Ho-Le (1988) defines classes and sub-classes according to the order of creation of nodes and elements. This scheme
15
3. OVERVIEW OF THE MESH GENERATION METHODS
Figure 3.1: Classification scheme for mesh generation methods. Adapted from Ho-Le (1988).
seems to be the most systematic, as it surveys structured as well as unstructured mesh generators, with triangular or other shaped elements. New algorithms can be included in the sub-classes. Particular regard is given to the methods for the triangular meshing of 2D domains, which is the focus of the present research, while the discretization in three dimensions is not considered. Four main classes, each of them with one or more sub-classes, are devised (Figure 3.1).
3.2.1
Mesh topology first
The set of elements, which represents the mesh topology, is created, followed by smoothing techniques for adjusting the node position. The identification of this class is merely nominal, as an algorithm for creating a mesh topology has necessarily to rely on other external mesh generator.
3.2.2
Nodes first
The nodes of the prospective mesh are created first, and then connected to form the elements. Two cases are given: a) Topology decomposition approach: a domain is first decomposed into coarse elements, i.e. by connecting its boundary vertices. Further subdivisions can be performed to satisfy the mesh density requirements (Figure 3.2). b) Node connection approach: it is a widely applied technique because it is conceptually simple. It comprehends two steps: the node generation and the following
16
3.2 Classification scheme of Ho-Le (adapted)
nodes connection into elements. The DT and the AFM in the formulation of Lo (1985) belong to this class. A description of different node generation methods is given in paragraph 3.4.
(a)
(b)
Figure 3.2: Topology decomposition of a 2D domain: (a) the first subdivision uses only the boundary polygon vertices; (b) the second subdivision inserts internal nodes.
3.2.3
Adapted mesh template
A regular, pre-generated grid template is created first and then fit into the domain polygon. Three sub-cases follow: a) Grid-based approach: a pre-generated rectangular/quadrilateral grid is superimposed to the domain (Figure 3.3). An example of this application is the quadtree approach (Baehmann et al., 1987). b) Mapped element approach (or transfinite mapping): a domain with three or four sides is first divided in sub-regions; a parametric grid template is then transferred, through a blending function, to each sub-region. The final mesh results from the union of the separately generated meshes (Figure 3.4). The method is analytically expressed through the projection theory, and may have multiple applications. In the algebraic generation, a sub-class of this case (see chapter VIII of Thompson et al., 1985), Lagrange polynomials generate curvilinear coordinate systems which can be fit to the domain boundary and to internal boundary conforming curves. The intersection of these lines generates the internal mesh points (Figure 3.5). In the
17
3. OVERVIEW OF THE MESH GENERATION METHODS
Figure 3.3: A simple example of an adapted mesh template with a grid-based approach (from Ho-Le, 1988). In this case triangular and quadrilateral elements coexist in the final mesh.
discrete transfinite mapping, the boundary curve is approximated as a sequence of points on a list, expressed by a unique coordinate. Most curves can be concatenated. An example is the use of cubic or higher order splines. c) Conformal mapping approach: it is a generalization of the mapped element approach, dealing with regions of more than four sides. Given the 2D simply connected straight sided polygon P, a second polygon Q is defined in the control space (see paragraph 3.4.4) with the same number of vertices of P and meshed with ideally shaped elements. The mesh of Q is then mapped back into P.
3.2.4
Nodes and elements generated simultaneously
Two different solutions have been developed: a) The geometry decomposition approach: a complex region can be recursively split, with insertion of nodes along the new inserted boundary, according to specified node density requirements (Bykat, 1983). b) The advancing front method, which is also called paving method (Peraire et al., 1987), developed the original formulation of Lo (1985), in which points are regularly
18
3.3 Unstructured triangular mesh generation methods
Figure 3.4: Example of a mapped element approach: a) a region is subdivided into convex elements, b) a unit mesh template is mapped into each element in order to form c) the final mesh (from Ho-Le, 1988).
placed inside the domain beforehand, by adopting a simultaneous generation of points and elements. See paragraph 3.3.2 for a complete presentation of this method. Table 3.1 compares the mentioned methods in terms of advantages and disadvantages.
3.3
Unstructured triangular mesh generation methods
Among the many methods of triangular mesh generation reported in the literature, two of them will be described here in more detail: the DT and the AFM. This choice comes from the fact that most mesh generators adopted in commercial and open source codes are derived from these two concepts (see Owen, 1998), and particularly in the field of coastal and ocean modelling practice. Brief details of other methods, namely the quadtree and the hybrid techniques, are given in sections 3.3.3 and 3.3.4.
3.3.1
The Delaunay triangulation
The mesh generation based on the Delaunay triangulation (DT) (Delaunay,1934) uses a simple geometry criterion for connecting a set of already existing points, in order to form a triangulation of non-overlapping elements. Strictly speaking, the DT only deals with node connection, assuming that the points to connect are somehow already
19
3. OVERVIEW OF THE MESH GENERATION METHODS
Table 3.1: Advantages and disadvantages of the mesh generation methods.
20
3.3 Unstructured triangular mesh generation methods
Figure 3.5: Mesh obtained by an algebraic system approach. The curves crossing the domain are expressed in terms of linear Lagrangian polynomials (from Trigo-Teixeira, 1994b).
present. For this reason, and in order to keep the explanation simpler, the creation of points is treated separately. For a given set of N points Pi , belonging to a convex planar region Ω, a set of sub-regions Vi can be drawn, where: Ω=
N [
Vi
(3.1)
i=1
in such a way that each Vi is the space closer to Pi than any other point. These convex, non-overlapping regions are called Vorono¨ı regions (Vorono¨ı, 1908). If every couple of points which belong to nearby Vorono¨ı regions, thus having a boundary segment in common, is connected by straight lines, a tessellation of Ω in triangles is achieved, called DT (Figure 3.6). The DT represents the dual graph of the Vorono¨ı tessellation for a given set of points. Among the many properties of the DT, two are particularly important for planar meshes: - a circle circumscribing any Delaunay triangle does not contain any other input points in its interior (in-circle criterion). This property is used to build algorithms for the triangulation; - as a consequence of the in-circle criterion, the DT maximizes the minimum internal
21
3. OVERVIEW OF THE MESH GENERATION METHODS
Figure 3.6: Delaunay triangulation and Vorono¨ı regions for a given set of points (from Thompson et al., 1999, page 1-13).
angle of the resulting triangles, tending to generate well shaped elements. When a planar region is to be meshed, different approaches can be followed: given a planar straight line graph (PSLG), that is a collection of boundary vertices connected by straight segments, a constrained DT of a PSLG is defined when each PSLG segment is forced to be present as a triangle edge in the DT. The triangles may result in not meeting all the Delaunay properties. A conforming DT (CDT) of a PSLG is a triangulation which maintains the Delaunay properties, where additional vertices, called Steiner points, are inserted in the original PSLG segments. Steiner points allow all the segments to belong to the final mesh and all triangles to maintain the Delaunay properties (Shewchuk, 1996). As already said, the DT itself is a triangulation method, which does not impose any requirement on the node placement. This is why DT algorithms usually include their specific node generation, insertion and deletion procedures. A popular approach is to first mesh the set of nodes belonging to the PSLG. As this mesh is usually coarse and irregular, new nodes are inserted incrementally into the existing mesh, and a DT is performed again. This process can continue until the mesh density meets the required refinement. The BatTri triangulator (Bilgili et al., 2005) is a typical example of triangulation for successive refinements.
3.3.2
The advancing front method
The advancing front method (AFM) uses the concept of active front, composed by a series of line segments on which elements are generated. The method was first applied
22
3.3 Unstructured triangular mesh generation methods
to the triangulation of an already existing set of nodes (Lo, 1985). In later versions, nodes and elements started to be generated simultaneously: in particular elements may be created one by one (Lee, 1999; Lee and Hobbs, 1999; L¨ohner and Parikh, 1988), or several at the same time (L¨ ohner and Cebral, 1999). Although several AFM versions have been developed, the basic structure of the algorithm consists of these common steps: a) after the domain boundary is discretized into line segments and boundary nodes are identified, the starting generation front is initialized. Nodes of internal boundaries are included as well. It is recommended to orient the front in such a way to leave the area to be meshed on the left hand side of the front; b) a new node is placed inside the domain and the element is formed; c) the active front is updated; d) steps b) and c) are repeated until the domain is completely meshed. Steps a) to c) are graphically shown in Figure 3.7, while step d) is shown in Figure 3.8. The methods for determining the node placement are various: paragraph 3.4.4 gives an overall description of frequent approaches for node generation and distribution, with particular regard to the concepts of node spacing function, meshing criterion, background mesh and metric tensor. The optimal offset of a new node from the segment belonging to the active front is a function of the algorithm input, and may be achieved by successive trials. As the AFM approach generates nodes and elements simultaneously, potential issues as the poor element area grading and the intersection of frontal active fronts (Figure 3.9) have to be taken into consideration. A complete description of an AFM algorithm is offered in chapter 5. A variant of the standard approach is the layered AFM (Pirzadeh, 1996; Garimella and Shephard, 2000), where a layer of elements is generated at once from a portion of the active front. The DT and AFM can also be combined together in the same algorithm. Mavriplis (1995b) suggested a mesh generator in which internal elements are generated through an AFM, and then connected by a DT. In Horritt (2000) the mesh generation is carried out by an AFM algorithm for the first layers bounding the domain, and then switches to a Delaunay triangulator for covering the unmeshed interior domain. The combined approach seems to be more time efficient than the AFM by itself, as more elements
23
3. OVERVIEW OF THE MESH GENERATION METHODS
Figure 3.7: Steps for updating the active front in the AFM: a) Initial active front; b) Identification of a new node; c) Front updated (from Thompson et al., 1999, page 17-5).
24
3.3 Unstructured triangular mesh generation methods
Figure 3.8: Successive stages in the implementation of an AFM technique (from Thompson et al., 1999, page 17-6).
Figure 3.9: Common issues during mesh generation performed by an AFM algorithm: connection/collision of frontal fronts (left); connection of elements of different sizes (right). From El-Hamalawi (2004).
25
3. OVERVIEW OF THE MESH GENERATION METHODS
Figure 3.10: Example of a quadtree recursive partition and related tree structure (from Thompson et al., 1999, page 15-3).
may be generated at once, to create well shaped elements, and to avoid the collision of two frontal fronts with different element sizes (El-Hamalawi, 2004).
3.3.3
The quadtree method
This method consists in superimposing a grid over the domain, recursively subdivided into an elementary shape with varying size in order to be consistent with a predefined spacing requirement. This process is based on a tree interrogation procedure, where the tree root corresponds to the starting undivided element, i.e a square. The tree sons are, at the first level, the four equal squares obtained by bisecting each side of the original square. Each of the four sons is interrogated and eventually subdivided into four other sons if the specific subdivision criterion is met. Graded meshes are obtained by placing nodes in cells of different levels of the tree structure. The partition ends when the subdivision criterion is met all over the domain. A simple quadtree algorithm example is given in Figure 3.10.
3.3.4
The hybrid method
Hybrid meshes are the result of the combination of different types of meshes. The aim is to achieve a better adaptation to complex geometries to be discretized, especially
26
3.4 Point creation
for surfaces and volumes. In two dimensions, unstructured with structured meshes, unstructured triangular and quadratic meshes, or unstructured triangular meshes produced with different methods (Trigo-Teixeira, 1994b) are common solutions. GFGEN (2001) is a renown mesh generator supporting the design of hybrid meshes.
3.4
Point creation
The set of points to which a triangulation method is applied can be generated in several ways. Here a summary of the most common methods is presented. It has to be mentioned that the task of node creation and element generation have been presented here separately, for the sake of simplicity, although in many algorithms the two operations are concurrent (as in the AFM), or a loop for node and element formation is run until a condition, like a spacing function, is satisfied (as in the meshing for successive refinements).
3.4.1
Independent node generation technique
A technique external to the triangulation method is applied. For example, if a circular region has to be meshed, a polar grid can be generated and adapted to the circle, and the points then connected to form the mesh.
3.4.2
Grid superimposition and successive subdivision
A regular grid is superimposed to the domain. This technique is called quadtree (see paragraph 3.3.3 for a more detailed explanation).
3.4.3
Point creation driven by the boundary point distribution
After the domain boundary is discretized into segments, whose vertices reflect geometrical and/or physical properties, points are created inside the domain extending the boundary point distribution with a smooth gradation (Frey, 1987; George and Seveno, 1994). In Figure 3.11 a river bed domain is split into sub-regions with different boundary point densities. The point generation algorithm distributes the nodes inside the domain in a uniform grading placement between opposite segments.
27
3. OVERVIEW OF THE MESH GENERATION METHODS
(a)
(b)
Figure 3.11: Point creation driven by the boundary point distribution for a river bed domain split into sub-regions. In (b) the point creation is driven by a bi-directional interpolation between two frontal boundary segments of each sub-region.
3.4.4
Point creation through a node spacing function
The node spacing function (NSF), also called size function, defines the target distance between nodes in the prospective mesh. This concept is widely used in the coastal modelling practice, as it enables to calculate the distance between nodes as a function, for example, of the bathymetry or flow conditions. Among the first publications explicitly focusing on the concept of NSF, Frey (1987) develops an algorithm for mesh refinement, where, after a DT of the boundary nodes is done, additional nodes are added one-at-a-time inside the domain: the insertion is accepted if maintains a valid DT and meets a specified NSF. The NSF can be either explicitly or implicitly defined. The explicit form is expressed as a function of the position, in the form: ρ = ρ(x, y)
(3.2)
where ρ is the target spacing and x and y the domain coordinates. ρ must be consistent with the distribution of the boundary nodes. In the implicit form, after the boundary nodes are placed in sufficient number to approximate well the domain boundary position and a node spacing factor is attributed to each of them, the NSF is constructed over the domain from the information stored on the boundary nodes. The definition of a NSF relies on the combined use of the following tools: a) one or several meshing criteria; b) a background mesh;
28
3.4 Point creation
c) a Riemann metric tensor, commonly used for anisotropic meshing. a) meshing criteria Node placement for domains over which the problem unknowns are highly variable requires different levels of mesh resolution. In order to meet a specific sizing requirement, a meshing criterion is used. For isotropic conditions, a meshing criterion corresponds to a scalar function to define the desirable element size throughout the domain in terms of distance between nodes. However a criterion may also express other geometrical element characteristics, as the shape elongation and the direction of elongations in anisotropic conditions (chapter 6). Two main types of criteria need to be distinguished: a-priori criteria, dependent on fields known before a problem solution is calculated, as the bathymetry/topography, the geometry, physical forcings or other physical domain properties, and a-posteriori criteria, derived from an error estimator, calculated from a previous model solution (Franca et al., 1995), or from the discretized governing equations (Hagen et al., 2000). In the last case a-posteriori techniques are always related to adaptive meshing modelling, where the nodes of a previously designed mesh are redistributed and/or new ones inserted in areas with the highest estimated errors, in order to minimize the error of the flow problem solution. Examples of recurrent meshing criteria for shallow water finite element model applications are listed in chapter 4. A particular definition of scalar NSF is when point creation is controlled by point and lines sources, which represent a way to locally enhance the mesh refinement. A point source is a punctual location inside the domain from where the node density decreases in the outward radial direction. A line source can be considered as a sequence of point sources. For a set of N point sources j with position P¯j , the local point spacing ρi can be defined as: ¯
Bj |x¯i −Pj | } ρi = minN j=1 {Aj e
(3.3)
where Aj and Bj are amplification and decay spacing parameters, |x¯i − P¯j | is the scalar
distance between the generic domain point i the point source Pj . Figure 3.12 exempli-
fies the result of point and line sources.
29
3. OVERVIEW OF THE MESH GENERATION METHODS
Figure 3.12: Point distribution controlled by point and line sources (Thompson et al., 1999, p.1-21).
b) Background mesh The node placement can be controlled by a background mesh underlying the domain to be meshed (Lee and Hobbs, 1999): this mesh is composed by elements at whose vertices a NSF, a size and stretching coefficients can be attributed. Small spacing is set for example in correspondence of zones where high gradients of the problem variables are expected. The size and shape of a mesh element are calculated by interpolation from the spacing coefficients of the background element it belongs to. The background mesh may also be represented by a previous modelling solution in an adaptive meshing process (Frey, 1987). c) the Riemann metric tensor In many computational applications, problems are highly dependent on preferential directions and may require elongated elements (see chapter 6). In the anisotropic mesh generation, the element characteristics cannot be fully addressed only by a scalar function. A 2 x 2 symmetric matrix, called in differential geometry the Riemann metric tensor, is therefore used. In this section the basis of the metric tensor approach are presented for the 2D case. An exhaustive description of a metric approach for a DT is found in Borouchaki et al. (1997a; 1997b).
30
3.4 Point creation
Given a domain Ω ⊂ R2 and a point x ¯ ∈ Ω, the metric tensor M (¯ x) is given by the
symmetric positive matrix
M (¯ x) =
a(¯ x) b(¯ x) c(¯ x) d(¯ x)
(3.4)
where a(¯ x), c(¯ x) > 0 and a(¯ x)c(¯ x) − b(¯ x) 2 > 0
(3.5)
The metric field M (¯ x) defines a Riemann structure on Ω. For any two points x¯1 , x ¯2 ∈ Ω, a real function L(M, x¯1 , x ¯2 ) is defined as: L(M, x¯1 , x ¯2 ) =
Z
1q
Q(x¯1 , x ¯2 , ε)T M (Q(x¯1 , x ¯2 , ε))Q(x¯1 , x ¯2 , ε)dε
0
(3.6)
where Q(x¯1 , x ¯2 , ε) = x¯1 + ε(x¯2 − x¯1 )
06ε61
(3.7)
In general L(M, x¯1 , x ¯2 ) is the distance between x¯1 and x¯2 with respect to the metric tensor M (¯ x). In case M (¯ x) = I, then L(M, x¯1 , x ¯2 ) reduces to the Euclidean distance x¯2 − x¯1 . For (3.4) and (3.5) M (¯ x) has positive real eigenvalues λ1 and λ2 which can be expressed as:
M (¯ x) = v¯1
T λ1 0 v¯1 v¯2 v¯2 0 λ2
(3.8)
where v¯1 and v¯2 are the eigenvectors of M (¯ x), with v¯1 ⊥ v¯2 . Since λ1 , λ2 > 0, it is possible to write:
λ1 =
1 l12
(3.9)
λ2 =
1 l22
(3.10)
and
Supposing that the length of the vector u ¯ = x¯2 − x¯1 is l1 in the v¯1 direction and l2
in the v¯2 direction, the length of u ¯ with respect to the metric will be, substituting in (3.6): L(M, x¯1 , x ¯2 ) = L(M, u ¯) = 1
(3.11)
Therefore the metric can be defined in Ω in such a way that controls the element characteristics of a mesh. For defining an anisotropic triangle, three parameters are needed: one versor (v¯1 or v¯2 ), and two scalars: l1 and l2 . The versor represents the privileged direction of stretching, and the two scalar parameters the element size along
31
3. OVERVIEW OF THE MESH GENERATION METHODS
the privileged direction and the direction normal to it. If l1 = l2 , the mesh will be isotropic. For (3.11), meshing a domain Ω over which a metric is defined corresponds to generating an ideal (= equilateral) mesh in a dual space, with the edge lengths tending to 1. The metric, performing the transformation of the point coordinates from the physical space to a control space, constructs an ideal size element and maps its coordinates back into the physical space. Application examples are given for DTs (Legrand et al., 2007) or AFMs (Lee, 1999).
3.5
Mesh quality evaluation and optimization
The main purpose of a meshing algorithm is to generate a good quality mesh, where the term “good” can be interpreted as able to limit the discretization errors of the governing model equations to a determined degree of approximation, while at the same time avoid computational overload. The concept of mesh quality is thus related to a specific modelling application. Hagen et al. (2000; 2001) demonstrated that, for triangular meshes applied to a 2D shallow water model, the localized truncation error decreases for equilateral triangles and limited area ratios between adjacent elements, as the odd-order truncation terms become less significant. On the other side, flow simulations with pronounced directional characteristics may be computationally more efficient if anisotropy conditions are imposed (Jansen et al., 2001). In this section the most recurrent triangle quality measures are reported. These may derive from purely geometry quantities, like the element angles, shape, size, and area grading, or from the level of reproduction of domain features and algorithm target size functions: - the interior angles are used as indicators of the element shape. Usually the extremal angles (maximum and minimum) are considered. For isotropy conditions, the angle divergence from the angle of the equilateral triangle is commonly used: δα ˆi = |
π −α ˆi| 3
(3.12)
where α ˆ i is the interior angle. In the literature shape factors derived from δ α ˆ i have been defined (El-Hamalawi, 2004). Some triangulation algorithm, as the DT, assures the maximization of the minimum interior angle;
32
3.5 Mesh quality evaluation and optimization
Figure 3.13: Quality index qi for different triangle shapes.
- the element equilaterality: it indicates the deviation in shape of an element from the equilateral triangle. It is defined as (Lee and Hobbs, 1999): √ 4 3Ai qi = 2 ai + b2i + c2i
(3.13)
where Ai is the area and ai , bi and ci the length sides of the element. qi measures how much a triangle is skewed from the isotropy condition: for an equilateral triangle qi = 1 (Figure 3.13). In case of a graded isotropic mesh, the mesh generator should be targeted to maximize qi ; - the in-circle to circum-circle radius ratio (Geuzaine and Remacle, 2009; P´ebay and Baker, 2003): γi = (
r1 sinˆ αi sinβˆi sinˆ γi )i = 2 r2 sinˆ αi + sinβˆi + sinˆ γi
(3.14)
being r1 the radius of the inscribed circle, r2 the radius of the circumcircle, α ˆ i , βˆi and γˆi the three inner angles of the triangle i (Figure 3.14). This parameter evaluates the isotropy/anisotropy grade of each element : an equilateral triangle has γi = 1, while for a degenerated triangle (collinear vertices) γi = 0. - the edge ratio (P´ebay and Baker, 2003): τi =
max{ai , bi , ci } min{ai , bi , ci }
(3.15)
P´ebay and Baker (2003) also report a series of other quality measures derived from
33
3. OVERVIEW OF THE MESH GENERATION METHODS
Figure 3.14: Notation used for the definition of the triangle quality indexes
possible combination of the elementary triangle geometric parameters. - the element area change: the area transition between one element and its neighbours is measured. A maximum threshold is necessary for limiting the estimation of the error truncation (Hagen et al., 2000); - the valence or degree number of a node: it counts the numbers of elements to which a mesh node belongs. An optimal valence of 6 is reported to increase the quality of the mesh topology (Field and Frey, 1995; Fortunato et al., 2011). Numerous quality indexes accounting for a deviation of the element from predetermined target functions, i.e. the size and shape of an ideal element, have been defined. Between them: - the maximum slope: the bathymetry gradient of an element is contrasted against a threshold value; - the quality parameter µi defines the triangle quality in relation to its deviation from the target size ρi (Lee and Hobbs, 1999): µi = where ρ˜i is the actual element size; the target element size is attained;
ρ˜i ρ˜i (2 − )qi ρi ρi
(3.16)
ρ˜i ρ˜i (2 − ) is a factor reaching a value of 1 when ρi ρi
34
3.5 Mesh quality evaluation and optimization
- similar to the previous one, the edge length quality ςi (Legrand et al., 2006) indicates the deviation of the length of the element side ρ˜i from the target length ρi : ςi = 1 − |
ρ˜i − 1| ρi
(3.17)
For an ideal case ςi = 1. Several other expressions of quality metrics calculating an element shape and size discrepancy from an ideal element can be found (Lee, 1999; Gorman et al., 2006; Power et al., 2006). If the generated meshes contain anisotropy elements, these metrics are usually derived with respect to a metric in a parametric space, where the target element corresponds to an equilateral triangle of edge length equal to 1: in these cases the quality measures defined for isotropic meshes are extended to Riemannian spaces (Dompierre et al., 2005). After the quality parameter is defined, an optimization process can be performed in order to enhance the mesh quality. This may include node connection changes, node movement, node insertion/depletion, or a combination of them, and can be usually related to purely geometric or adaptive processes. As the mesh quality depends on the particular application, the subsequent mesh enhancement options will be defined and applied accordingly to the given measure of element quality. Many optimization algorithms, especially the ones belonging to successive refinement approaches, are actual node generation/insertion methods. Adaptive re-meshing obtained from an a-posteriori error function also fall within this category. Common topology modification procedures are: - the Laplacian smoothing (called in this way because the formula for repositioning nodes derives from a finite difference approximation of the Laplace equation): each node not located along the boundary is moved toward the centroid of the polygon formed by the adjacent triangles (Figure 3.15): N 1 X x ¯P i x ¯P = N
(3.18)
i=1
where x ¯P is the new position of point P and x ¯P i the positions of its neighbouring points. (3.18) calculates the coordinates of the new node position simply as the average of the coordinates of the neighbouring nodes. A variation of (3.18) introduces
35
3. OVERVIEW OF THE MESH GENERATION METHODS
a weighted coefficient of relaxation (e.g.: George and Seveno, 1994). Other expressions reposition the internal nodes in order to increase the minimum internal angle or maximize the worst element quality (Ollivier-gooch and Freitag, 1997);
(a)
(b)
Figure 3.15: Example of Laplacian smoothing: a) node P is visualized in its original position with all neighbour elements. b) node P is repositioned to the centroid of the polygon;
- the diagonal swapping: the common edge of two neighbour triangles is swapped (Figure 3.16); The criterion defining whether or not a diagonal should be swapped may be different: the optimization of the node valence to 6 (see chapter 5), the conformity with the Delaunay condition (Figure 3.17), or the improvement of the equilateral quality measure of the involved elements (Geuzaine and Remacle, 2009); - the deletion of elements with bad aspect ratio: this is specially applied to bad shaped elements formed along curved boundaries. A node reconnection procedure follows; - the local mesh refinement (also known as mesh enrichment): in areas where pronounced gradients of the flow variables are expected, elements may be split in two or four, or divided by the insertion of breaklines. As an example, after the bathymetry gradient is evaluated, an element-by element refinement can be applied in order to meet the requirement of the maximum allowed element slope.
36
3.5 Mesh quality evaluation and optimization
Figure 3.16: Diagonal swapping for two adjacent triangles
Figure 3.17: Diagonal swapping for obtaining triangles satisfying the Delaunay condition. (from Legrand et al., 2000).
37
4
MESHING CRITERIA FOR SHALLOW WATER MODELS 4.1
Introduction
In this chapter a review of the meshing criteria encountered in the literature of ocean and coastal modelling is presented. The review is limited to 2D unstructured finite element model meshes, for both isotropic and anisotropic conditions. The mesh generation and mesh editing are a critical step in the set-up of a numerical model. The definition of a meshing criterion and of the element shape/sizes does not have an absolute validity, but has to be related to the particular case under investigation and to other model settings. For example, the problem of coastline resolution deals also with the type of coastal boundary conditions, and the proper definition of both elements is necessary (Jones and Davies, 2005). Keeping in mind these considerations, it seems a good starting point to understand how modellers resolved frequent meshing issues related to ocean and coastal domains. The domain dimensions of the given examples vary over several orders of magnitude, from river estuaries and urban areas to the ocean scale. Criteria which include the effect of the Earth curvature on the domain resolution requirements have not been reported, as in those cases the condition of planar surface is not respected. This effect should however be considered for polar regions and for the modelling of global ocean dynamics.
39
4. MESHING CRITERIA FOR SHALLOW WATER MODELS
4.2
Tidal modelling
1D early attempts to discretize spatially the SWEs for tidal waves propagating over constant depth (Gray and Lynch, 1977; Le Provost and Vincent, 1986) identified a criterion based on the constant ratio between wavelength (λ) and mesh element size (∆x): R=
λ ∆x
(4.1)
with R ranging form 15 to 25. The tidal wavelength can be expressed, using the inviscid long wave theory formulation, by the shallow water relation: λ=T
p gh
(4.2)
where h is the water depth, g the gravity acceleration and T the tidal period. As the wavelength decreases for a wave moving into shallower waters, the element size has √ to decrease as h to maintain a constant λ\∆x ratio. In this way a graded mesh is generated, with an increasing resolution in correspondence of the shallower domain areas, where major wave modifications are expected. However later studies (Westerink et al., 1992a; Westerink et al., 1992b; Westerink et al., 1994; Westerink et al., 1995), carried out on a 2D finite element tidal model, pointed out that the wavelength to mesh size ratio criterion is insufficient by itself to provide a convergent and optimally graded mesh. The hydrodynamic computations were performed using the generalized wave-continuity equation (Lynch and Gray, 1979) for the finite element solution of the SWEs, included in the ADCIRC code (Luettich et al., 1992). The determination of an optimal mesh was done through convergence tests, in which the computations obtained with each mesh are compared with a “true” solution given by a very finely resolved uniform grid. The goal was to determine a graded mesh for which the flow computational solution is comparable with the finest grid, but with a significantly smaller number of nodes. This research pointed out major shortcomings for the wavelength to mesh size ratio approach: - the nature of the criterion is mainly one-dimensional, that is for monochromatic waves over constant depth, hence not suited to capture the bi-linear characteristics of tidal flows, i.e. amphidroms or tidal currents;
40
4.2 Tidal modelling
- a resolution proportional to the wavelength gives a lower node spacing than the one required to properly resolve steep slopes, as across the shelf break, canyons or seamounts; - shallow regions, as the continental shelf of the Gulf of Mexico and the Caribbean Sea, need a higher refinement level than deeper waters in order to maintain the same level of numerical precision; on the other side, if the λ\∆x ratio is set to a number which ensures a sufficient resolution for the shallower portions of the domain, deep waters will subsequently be over-resolved, leading to an unnecessary computational cost; - a refined representation of the coastal boundary may positively influence the computed tide levels. The convergence studies showed that the location of the highest errors were localized near the Atlantic coast (for the absolute amplitude) and near amphidromes (for the component phase), that is where high gradients of the mentioned variables are found. An overall good refinement is proposed by setting a λ\∆x ratio of not less than 25 and up to 150 in deep water regions. Additional resolution, in the range 50 ≤ λ\∆x ≤ 100
(Westerink et al., 1994) or 250 ≤ λ\∆x ≤ 500 (Westerink et al., 1995; Luettich and Westerink, 1995) should be achieved along the shelf break and shelf slope. Examples
of meshes for the Western North Atlantic, Gulf of Mexico and Caribbean Sea are presented in Figure 4.1. The studies mentioned so far deal with mesh convergence trends for a given bathymetry, supposing that the most detailed available bathymetric dataset is used. Luettich and Westerink (1995) separated the effect of an increased mesh resolution from the effect of an increased bathymetry resolution. Two sequences of convergence tests were performed: the first series assessed the convergence of the model results for uniform grids of increasing resolution interpolated with the same bathymetry scatter set (75 Km resolution). Then the most refined grid of the first sequence was tested with bathymetry scatter sets of increasing resolution. It was concluded that the bathymetry resolution is critical for achieving convergent numerical results, as an improved representation of steep slopes and local bottom features influence the model accuracy. This conclusion is however valid only if the mesh spacing was small enough to catch these features.
41
4. MESHING CRITERIA FOR SHALLOW WATER MODELS
(a)
(b)
(c)
Figure 4.1: Unstructured graded finite element meshes of the Western North Atlantic, Gulf of Mexico and Caribbean Sea based on a λ\∆x ≥ 25: a) coarse representation of the coastal boundary; b) fine representation of the coastal boundary; c) coarse representation of the coastal boundary and additional resolution in the vicinity of the shelf break (from Westerink et al., 1994).
42
4.3 Resolution requirements for steep bathymetric gradients
4.3
Resolution requirements for steep bathymetric gradients
Considering the additional mesh refinement required along the shelf break, it is clear that other meshing criteria should be introduced in order to resolve properly areas of steep topography: for ocean domains, these include, apart from the shelf break and slope, seamounts, islands, shoals and bank sides. Local circulation patterns and internal tides are expected around these features (Xing and Davies, 1998), which may compromise the performance of a modelling encompassing bigger spatial scales. Hannah and Wright (1995) analysed the source of errors between an analytical and a numerical model solution for wind-driven flows. They pointed out that both unresolved physical (i.e. the pressure gradients or the eddy viscosity layer for the vertical column) and topographic variables (steep slopes) introduce errors and that all length scales should be properly discretized. A discretization driven by the topographic length Lh is suggested: ∆x ≤ αLh
(4.3)
being Lh ≃
h ∆h/∆x
(4.4)
where ∆x is the mesh spacing, h the water depth, ∆h the vertical depth change per element and α the refinement coefficient. Lynch et al. (1995) applied the topographic length scale refinement for modelling the sub-tidal circulation of a shallow sand bank, adopting for α the range of values [0.1 - 1]. Gorman et al. (2006) suggested a re-meshing technique, valid for isotropic and anisotropic meshes, where an optimization algorithm is based on the estimation of a bathymetry derived error function. Legrand et al. (2007) developed an anisotropy meshing approach tested for an ocean domain, encompassing continental shelf, shelf break and ocean shelf. Two refinement criteria, one based on the hydrodynamics and one on the bathymetry, are combined together: the hydrostatic consistency condition for steep bottom slopes, and the tide celerity for gentle slope gradients. A metric tensor defines the three parameters needed for the anisotropic meshing:
43
4. MESHING CRITERIA FOR SHALLOW WATER MODELS
- a size function L// proportional to the shallow water wavelength: L// (x, y) = L∗
r
h(x, y) h∗
(4.5)
where h(x, y) is the water depth at the (x, y) location and L∗ a user specified target element size for a given column depth h∗ . - a size function L⊥ , normal to L// , and proportional to the slope gradient: L⊥ (x, y) = σ|
~ ∇h(x, y) −1 | h(x, y)
(4.6)
where σ is a coefficient related to the vertical grid discretization of the momentum ~ equation, and ∇h(x, y) the bathymetry gradient. - a privileged direction v¯(x, y) set along the direction normal to the bathymetry gradient. The orientation of this vector can be set also along other directions, i.e. the major axes of the tide ellipse. Minimum and maximum bounds are set to avoid over-refining in the shallowest areas or under-refining where the slope gradient is low. In this way in areas where the slope is steep the node density is increased along the bathymetry gradient direction, while on the perpendicular direction the discretization is driven by the wavelength function; in areas of limited slope, the mesh tends to be isotropic with size based on the local wavelength (Figure 4.2). Suggestions of other meshing criteria are proposed: for example, the element target size could also be set proportional to the semi-axis of the tide ellipses, in order to resolve the circulation patterns along narrow channels.
4.4
Hurricane and storm surge modelling
Blain et al. (1998) investigated how hurricane modelling, in terms of storm scale, path, forward velocity and coastline refinement, drives the mesh requirements. The major finding was that the inverted barometer function, reproducing the pressure field of the moving storm, must be properly reproduced (Figure 4.3) in deep waters. In shallower waters the interaction with the bottom increases and big improvements in the simulation results rely on the enhancement of refinement of the coastline and of the shallowest
44
4.4 Hurricane and storm surge modelling
Figure 4.2: Anisotropic mesh for an idealized bathymetric profile including three regions that may be regarded, from the left to the right of the top figure, as the continental shelf, the continental slope and the ocean shelf. The mesh element sizing considers the requirements of the tide wave propagation and of the relative slope gradient (from Legrand et al., 2007). The variable D for the water depth appearing in this figure, taken from the original article, corresponds to the water depth h used in the present work.
45
4. MESHING CRITERIA FOR SHALLOW WATER MODELS
Figure 4.3: Piece-wise representation of the inverted barometer forcing function over a discretized 1D domain (from Blain et al., 1998).
part of the domain. A critical area is around the storm landfall point, where major errors have been found.
4.5
Domains with complex coastline patterns
Irregular coastline geometries or domains with several internal discontinuities like archipelagos or reefs require specific attention when defining the level of discretization of the boundaries and of the adjacent areas. The hydrodynamic interaction between flow and boundary geometry and shape may determine processes which propagate over a wide range of spatial scales (Wolanski et al., 2003). Although the proper level of refinement depends on the purposes of the application, a better coastline discretization corresponds to a more refined representation of the shallowest bathymetry and to improved simulation results (Blain et al., 1998) (Figure 4.4), in conjunction with proper boundary conditions (Jones and Davies, 2005). Kliem and Greenberg (2003) found that the accurate representation of channels and straights connecting the Arctic Sea with the North Atlantic, in both width and cross-section profiles, is critical for the correct estimation of flow transport between the two oceans (Figure 4.5). In Figure 4.6 a BatTri application shows a series of bathymetry based successive refinements of the inlet connecting Lake Maracaibo to the Gulf of Venezuela. An island as well appears in the middle of the inlet. In this case a suitable mesh refinement criterion would set a NSF proportional to the narrowness of the strait/inlet. Until now however the author
46
4.5 Domains with complex coastline patterns
Figure 4.4: Discretization of a hypothetical sinusoidal coastline using (a) five, (b) nine and (c) 17 points per wavelength (from Blain et al., 1998).
does not know of any code which implements this option automatically. When two water bodies are connected by a narrow geometry, as an inlet, a hybrid mesh generation, which splits the computational domain into geometrically simpler regions which are meshed separately, is sometimes applied. Trigo-Teixeira (1994b) applies this ´ technique to a multi-scale domain of the Obidos lagoon (Portugal): an AFM generator is adopted to mesh the lagoon and the nearshore ocean region, where the NSF is defined by hand on a background coarse mesh. The interconnecting channel is meshed trough an algebraic interpolation function, in order to obtain a regular node location pattern following the tide streamlines (Figure 4.7). A specific algorithm is applied to connect the separately generated meshes into one. The combined use of different algorithms represents in this case a successful solution to address specific meshing requirements to different sub-regions, although it adds the extra operation of linking the disjoint meshes. Around islands the flow motion may assume jet or eddy-like patterns. Legrand et al. (2006) defined a Delaunay triangulation algorithm in which the target element function ∆x is a linear combination of two criteria: the shallow water wavelength and a blending
47
4. MESHING CRITERIA FOR SHALLOW WATER MODELS
Figure 4.5: Discretization of a highly irregular and complex region in the Artic Sea. The resolution ranges from 1.1 km in the narrow straits to 53 km in Baffin Bay (from Greenberg et al., 2007).
Figure 4.6: Example of a mesh covering the region of Lake Maracaibo (on the left). The square on the upper right zooms the inlet of the mesh on the left, and the bottom square represents an example of a possible local refinement.
48
4.5 Domains with complex coastline patterns
(a)
(b)
´ Figure 4.7: Obidos lagoon: (a) generation of different meshes for different spatial scales and (b) merging algorithm. From Trigo-Teixeira (1994b).
49
4. MESHING CRITERIA FOR SHALLOW WATER MODELS
function proportional to the distance d(x, y) from the nearest island or reef: ∆x(x, y) = ∆xM IN + Θ(d(x, y))
s
h(x, y) (∆xM AX − ∆xM IN ) hM AX
(4.7)
where ∆xM AX and ∆xM IN set respectively the maximum and minimum element sizes. The blending function Θ(d(x, y)) is expressed in terms of a piece-wise cubic polynomial, and is equal to 0 in proximity of the shoreline, and to 1 for distances larger than the recirculation eddies scale lengths.
4.6
In confined waters and low-lying areas
When the flow propagates over very shallow areas and interacts with the coastline morphological features, like estuaries, inlets, lagoons, wetlands and urbanized low-lying areas, the actual topography must be reproduced at a level which allows capturing the physics of the process and the local details of the flow. At this scale, meshing criteria for open regions, dependent on the wavelength, on the topographic length scale, or on the combination of both, are not suitable any more. A meshing approach based on multiple requirements is needed, based on the peculiar needs of each case: Jones and Richards (1992) propose the manual insertion of new nodes and the use of diagonal swapping for the inclusion of hydraulic structures in the mesh. In coastal and river inundation studies, given a high resolution digital elevation model (DEM) providing water and land coverage, the proper reproduction of small scale features likely to influence the water propagation, such as channels, levees, roads, railways, buildings, etc., is recommended (Blanton and Luettich, 2008; Westerink et al., 2008; Shen et al., 2006; Bunya et al., 2010). Xu et al. (2010) use a hybrid triangular and quadrangular mesh where elements covering the low lying coastal areas coincide with the DEM grid cells. In order to improve the representation of the modelled flood event, specific algorithms have been developed to extract significant points from a LIDAR dataset for the identification of vegetation and man-made features (Cobby et al., 2003; Bates et al., 2003). Figure 4.8 shows a discretization of the Mississippi estuary, while Figure 4.9 shows a highly detailed discretization, including the channel and embankment network, for the city of New Orleans.
50
4.7 A-posteriori meshing techniques
Figure 4.8: Southern Louisiana and estuary of the Mississippi river covered by a graded mesh. The element size is shown on the right (from Bunya et al., 2010). Resolution on the continental shelf is determined by the wavelength, the topographic length and storm size criteria. The mesh applies localized refinement to major hydraulic features, down to 50 m element size for critical channels and conveyances.
4.7
A-posteriori meshing techniques
A-posteriori meshing schemes derive the NSF from the model solution obtained with a previous attempt to discretize the domain. Error estimators and the Courant condition fall into this category.
4.7.1
Error estimators
Classical error estimators of the Navier-Stokes equations identify the residual of the flow solution on an element-per element basis (Oden et al., 1990; Franca et al., 1995; Carey, 1995). In the localized truncation error analysis, called for brevity LTEA, (Hagen et al., 2000; Hagen et al., 2001; Hagen and Parrish, 2004) the error is associated to the discretization of the model governing equations, and the node placement is optimized by setting an uniform distribution of the truncation error all over the domain. In this way the physics of the shallow water flow are incorporated directly into the algorithm of mesh generation. One example of LTEA application is shown in Figure 4.10, where the same domain, corresponding to the Gulf of Mexico, is discretized by two different criteria: the mesh on the left is generated by the LTEA method, while the mesh on the right by keeping a constant wavelength to mesh size ratio. Although the two meshes
51
4. MESHING CRITERIA FOR SHALLOW WATER MODELS
Figure 4.9: Mesh portion covering the city of New Orleans for an ADCIRC application (modified from Greenberg et al., 2007). The meandering feature crossing the map horizontally is the Mississippi river.
have the same number of nodes, the LTEA derived mesh produces a better agreement with tide data.
4.7.2
A numerical constraint: the Courant condition
Spatial discretization and numerical stability of a model which solves PDEs in an explicit time scheme are interrelated: it must be ensured that the model is able to “catch” the flow transient between two adjacent points, all over the mesh, otherwise computational instability issues may rise. For example, in the case of a wave propagating in a 2D domain, the model time step should be a fraction of the time required to the wave to cross whatever cell, and, as a consequence, the computational time step is controlled by the smallest element of the mesh. This requirement is defined by the Courant-Lewy-Friedrichs condition (Courant et al., 1967), expressed in one dimension as: ∆t ≤
∆x V
(4.8)
and the corresponding Courant number Cr : Cr =
∆t ∆x V
52
≤1
(4.9)
4.8 Conclusions
(a)
(b)
Figure 4.10: Gulf of Mexico discretized (a) with the LTEA approach and (b) with the wavelength to mesh size ratio criterion. The mesh on the left showed better tide model predictions (from Hagen et al., 2001).
where ∆t is the time step (in s), ∆x the local mesh length scale (in m), and V the estimated flow velocity (in m/s). When ∆x is lowered, ∆t must decrease as well in order to satisfy (4.8). For non-linear applications in deep ocean flows, a practical upper limit of Cr is 0.5; however, in case of more demanding simulations including high flow gradients or wetting and drying schemes, an even tighter constraint (Cr < 0.1) must be imposed (Blain et al., 2010).
4.8
Conclusions
In this chapter a review of a-priori and a-posteriori meshing criteria for finite element shallow water models has been presented. The a-posteriori approach relies on a preexisting discretization, which is coarsened/refined according to the prescription of an error estimator or of the flow velocity module. For a-priori criteria systematic mesh convergence studies are carried out in order to quantify the improvement of the flow solution for meshes with different local or global refinements. No indications are usually given about the “absolute” best element size to use across the domain: although the criterion is the same, the associated level of refinement is devised case by case, according to the computational model used, the applied forcing functions, the boundary conditions and the bathymetric features. In general, these conclusions emerge:
53
4. MESHING CRITERIA FOR SHALLOW WATER MODELS
a) for computational domains covering multiple spatial and hydrodynamic scales, unstructured meshes are preferable to structured ones because more flexible in meeting varying resolution requirements; b) a resolution proportional to the wavelength gives an insufficient node spacing to properly resolve steep slopes, as across the shelf break, canyons or seamounts; c) critical refinement zones with respect to the accuracy of the solution are identified in correspondence of high bathymetry changes (shelf break and shelf slope), features such as inlets, rivers, navigation channels, levee systems and emerged structures likely to drive the flow in wet and dry studies, and complex coastline geometry; d) higher refinement is needed where the 2D flow propagation pattern is pronounced (amphidromes, shallow waters); e) the pressure forcing field needs to be properly reproduced along the storm trajectory and at the land fall point.
54
5
UNSTRUCTURED MESH GENERATORS 5.1
Introducton
The fundamental approaches for the discretization of 2D planar surfaces in triangular elements, as accounted in chapter 3, are conceptually simple, based on elementary geometrical properties. However a surprising flourishing of meshing algorithms, both commercial and open source, has been developed (see the reviews of Owen, 1998, and Lo, 2002). Looking at the literature of the past decades, algorithms have moved towards higher levels of automation and flexibility, going hand by hand with the available increased computational power and graphical capabilities: operations such as handling of non-trivial domain shapes, constrained triangulations, complex boundaries, node placement/replacement control options, post-processing operations for the mesh quality enhancement, etc., have been gradually evolved, made robust and included in user friendly graphical user interfaces (GUI). The range of available mesh generators is therefore huge and at first sight disorientating, especially for those not acquainted with the subject. The intention of this chapter is to retrace the investigation carried out in the first phase of this research, when several softwares were analysed. Priority is given to the description of the new version of MESHGR, an existing advancing front mesh generator (Trigo-Teixeira, 1994a). The presentation of other unstructured mesh generators, which have been tested as potential alternative meshing tools, is included. A final comparative analysis is made in order to identify a finite element discretization
55
5. UNSTRUCTURED MESH GENERATORS
tool suited for the specific requirements of the present investigation.
5.2
MESHGR
MESHGR is an AFM algorithm which generates 2D triangular unstructured meshes for finite element applications. In the version presented here the code, which has been written in Fortran and has been developed since the 1980s, has been improved in order to cope with multi-scale domains, complex boundary geometries and NSFs derived from multiple criteria, as often encountered in the ocean and coastal modelling practice. Documented past applications of the algorithm (Trigo-Teixeira, 1994a; Trigo-Teixeira, 1994b) were dealing with simple geometries and rough boundary discretizations, a maximum number of elements limited to 5000, reduced variation in the element sizes and the use of a NSF defined manually and not derived from geometric or bathymetric criteria. When the domain extended over different spatial scales, a series of meshes had to be created separately and then merged together by an auxiliary external program. A series of improvements have been introduced to MESHGR, in terms of memory allocation, user-friendly GUI, resolution of internal bugs, pre-processing of the input data files, post-processing routines, and export options of the output. The general idea which addresses the current research effort is to develop and make available to the scientific community a user-friendly tool which produces high quality, ready-to-use unstructured triangular meshes, with a special target for oceanic, coastal and estuarine domains, but which is also applicable to other fields of the FEM. The improved algorithm has shown to overcome the limitations of the previous version: wide multi-scale real domains have been successfully discretized, and the produced meshes have been applied into finite element model applications (Mazzolari et al., 2013; Flor et al., 2013). In the current chapter an overall algorithm description is given, while meshing tests for idealized and real domains are reported in chapters 6 and 7 respectively. An alternative implementation of the AFM, which introduces the concept of pseudo-islands, is presented in chapter 8.
5.2.1
Improvements introduced in MESHGR
The MESHGR input is a boundary polyline for the definition of the computational domain and an underlying background mesh for the definition of the NSF. A number
56
5.2 MESHGR
of alterations to the original MESHGR code became necessary in order to obtain a more robust, flexible and efficient tool able to cope with complex multi-scale domains and boundary geometries, as often encountered in the current ocean and coastal modelling practice. These changes have been summarized in the following list: - change from a sequential code structure to a modular structure: the variable initialization, the main body of the code, the graphical subroutines, dialogues and messages were separated into different modules; - passage from a static memory allocation to a dynamic memory allocation. This change allows to increase the input and output mesh dimension from the previous 5 000 elements to over 106 elements. However, as the number of nodes and elements of the prospective mesh is not known beforehand, and the Fortran programming language requires the array size to be declared somewhere inside the code, a default maximum dimension of 2 ∗ 106 elements and 106 nodes is set for the output mesh; this parameter may be increased according to the operating system and machine memory storage capabilities; - separation of the background mesh input file from the domain input file, in order to prepare or modify the files independently; - introduction of a pre-processing program which automates the preparation of the MESHGR input files in the Surface-water modelling system SMS (Aquaveo, 2000) environment: in particular an ADCIRC (Luettich et al., 1992) unstructured mesh is re-written into a background mesh file and a SMS boundary polyline into a domain file; - iterative loop for applying post-processing quality and checking options, with possibility to export the mesh at each step in ready-to-use formats of some widely used models: ADCIRC, RMA2 (Donnel et al., 2006) and CGWAVE (Demirbilek and Panchang, 1996); - elimination of unused variables and parts of the code; - modification of the algorithm default parameters, in order to increase the range of their admissible values. In this way the algorithm demonstrated greater flexibility in handling complex boundaries and highly variable NSFs;
57
5. UNSTRUCTURED MESH GENERATORS
Figure 5.1: Flow chart of the GUI interfaced with MESHGR
- introduction of the pseudo-island option (see chapter 8); - substitution of the subroutines of the GPACK GUI (Rohlf et al., 1987) with the Windows-based Winteracter (2009) library. A flow chart visualizing the use of MESHGR interfaced the new GUI is represented in Figure 5.1, while the layout of the new GUI is shown in Figure 5.2.
5.2.2
Structure of MESHGR
The mesh generation process follows a hierarchical framework: this means that, after the definition of the domain and of its sub-regions, the starting and ending points of the boundary curves (called vertices) are saved for first in the prospective mesh, followed by the nodes interpolated along the boundary curves and at last by the nodes inside the domain. These successive steps are identified: a) definition of the background mesh;
58
5.2 MESHGR
Figure 5.2: GUI provided to MESHGR. The left sub-window shows the background mesh and its numbered nodes, the right sub-window the generated mesh. The box in the middle asks to the user which option to perform.
59
5. UNSTRUCTURED MESH GENERATORS
b) definition of the domain geometry; c) discretization of the domain boundaries according to the background mesh specifications; d) discretization of the domain; e) mesh quality enhancement options. In this section each of the mentioned points is explicitly treated, with the most significant parts of the algorithm written in pseudo-code. a) Background mesh The first part of the input dataset is represented by the definition of a triangular element background mesh B which covers the whole domain extension. B addresses the element characteristics of the prospective mesh, in terms of element side length and isotropy/anisotropy conditions. Every background mesh node i is accompanied by a series of attributes: the Cartesian location on the plane (xi , yi ), a scalar corresponding to the target NSF (ρi ), an anisotropy coefficient (si ), and two values (ui and vi ) defining the element orientation on the plane in case an anisotropic condition is set. ρi may derive from a NSF representing one or a series of meshing criteria defined by the user (Mazzolari et al., 2012). The NSF is interpolated externally to MESHGR to the nodes of B. si , ui and vi determine whether a mesh is isotropic or not: si defines the element shape (si =1 for an equilateral triangle); ui and vi represent the components of the elongation vector in case the element is not equilateral (si < 1). B is read in this sequence: a.1) declaration of the number of nodes (NPOIG) and elements (NELEM); a.2) for i =1, NELEM, definition of the mesh topology; a.3) for i =1, NPOIG, declaration of xi , yi , ρi , si , ui and vi .
b) Definition of the domain The second part of the input dataset contains the geometrical description of Ω (Figure
60
5.2 MESHGR
Figure 5.3: Example of a non-simply connected domain formed by 3 adjacent regions(R1 , R2 and R3 ) and an internal opening into R2 . The arrows indicate the segment orientation, set in a way that keeps the domain on the left side of the segment. Segment and arc extremities are evidenced as well.
2.2), in terms of domain node locations and boundaries. Ω may be simply or non-simply connected, formed by n ≥ 1 non-overlapping regions Ri : Ω=
n [
Ri
(5.1)
i=1
A set of boundary curves defines Ω and its regions. Each curve is a closed polyline, formed by segments and/or arcs of circle, not self-intersecting and not intersecting any other boundary curves. Adjacent regions may share the same boundary curve. A curve is defined by a starting and an ending point, called vertices, oriented in a way that respects the convention of the left side: this means that, if the curve is run from the starting point to the ending point, the domain remains to the left. It follows that external boundaries are oriented anti-clockwise and internal boundaries (i.e. internal openings) clockwise. An example of an idealized non-simply connected domain formed by 3 regions is presented in Figure 5.3. Ω is read inside MESHGR in this order: b.1) declaration of the number of regions (NREG), vertices (NBV) and boundaries (NBCS) delimiting Ω;
61
5. UNSTRUCTURED MESH GENERATORS
b.2) for i =1, NBV, declaration of the coordinates of the boundary segments/arcs vertices; b.3) for i =1, NBCS, definition of the boundary type (0=linear, 1= curvilinear) and of the vertex number at its extremities. If the segment is curvilinear, the coordinates of the arc center and the radius are declared as well; b.4) for i = 1, NREG, declaration of the number of domain segments delimiting the ith region and of their sequential number.
c) Discretization of the domain boundaries Before the AFM is applied for the discretization of the interior domain, the discretization of the domain boundaries is performed according to the background mesh specifications. The interpolated points, together with the vertices of the domain boundary segments, form the initial advancing front Γ0 . The following steps are performed: c.1) the segment vertices are saved as the first points lying on Γ0 ; c.2) for i = 1, NBCS: c.2.1) read the initial vertex. Call it j ; c.2.2) detect the background mesh element where j lies; c.2.3) interpolate the meshing parameters ρj ,sj , uj and vj ; c.2.4) according to the meshing parameters of c.2.3), find the next interpolated node j+1 along the segment and save it in Γ0 ; c.2.5) j=j+1 ; c.2.6) if the end of the segment is reached go to the next segment; if not go to c.2.3). At the end of this loop Γ0 is initialized. Even if curvilinear arcs are declared in the domain initialization, the process of boundary discretization will generate just linear segments in closed loop(s). The series of nodes belonging to Γ0 represents the node string located along the boundaries of the final mesh.
62
5.2 MESHGR
d) domain discretization After the initialization of Γ0 , a method to determine the departure zone, that is the portion of the active front from which an element is created, is set (see point d.1). For supporting the creation of anisotropic meshes, a bi-univocal correspondence ′
′
(x , y ) = T(x, y) between the physical space and a control space (Borouchaki and Frey, 1998; Legrand et al., 2007) is defined as a function of the meshing parameters. This transformation maps the physical space (x and y coordinates) into a normalized space ′
′
(x and y coordinates) where the target elements are equilateral triangles: x′ = u(ux + vy) + v(vx − uy) s T(x, y) = y ′ = v(ux + vy) − u(vx − uy) s
(5.2)
In this way the problem of creating anisotropic elements/meshes in the physical space is reduced to the problem of creating equilateral elements in the control space. In the literature the triangles of the normalized space have unit side length, while MESHGR just requires the triangles to be equilateral, with the size defined by the NSF in the physical space. The inverse transformation T−1 maps back the node coordinates for the final mesh in the physical space. Figure 5.4 gives a visual interpretation of T(x, y). ′
′
For the isotropic case it can be easily shown that x = x and y = y. The domain discretization in MESHGR is performed according to the following procedure (the case of a domain formed by a single region is presented (NREG=1) for the sake of clarity, although the description can be extended to n regions): d.1) choose the method for implementing the evolution of the advancing front: - front circulating: the order in which elements are formed follows the spatial sequence of the segments along the active front. The front is run in the direction that keeps the domain on the left side; - front searching: at each step the element is formed from the shortest segment (in the normalized space) of the active front; d.2) find the shortest segment AB ∈ Γ0 ; d.3) calculate the meshing parameters of AB: sAB = min(sA ; sB )
63
(5.3)
5. UNSTRUCTURED MESH GENERATORS
Figure 5.4: Visual interpretation of the effect of the transformation T(x,y) between the physical and the control space.
ρA + ρB 2 uA + uB = 2 vA + vB = 2
ρAB =
(5.4)
uAB
(5.5)
vAB
(5.6)
2 uAB and vAB are normalized in order to have u2AB + vAB = 1;
d.4) map AB into the control space using relation (5.2); d.5) calculate the target NSF corresponding to the side length for the prospective element as: ρN
ρN = 2AB 0.55AB
if 0.55AB ≤ ρAB ≤ 2AB if ρAB > 2AB if ρAB < 0.55AB. ′
(5.7)
d.6) calculate the coordinates of the prospective element node N (in the control space); d.7) map back into the physical space. Find N (Figure 5.5); d.8) form the searching circle C˜ with radius equal to 0.8ρN (Figure 5.5);
64
5.2 MESHGR
d.9) detect all Qi ∈ Γ internal to C˜ and call this list C(Qi ) (Γ is the advancing front
and coincides with Γ0 only when the first mesh element is created). Order C(Qi ) from
the nearest to N to the furthest; d.10) loop over C(Qi ) to check if the connection of Qi with AB is possible. A connection is possible if these two conditions are satisfied: \ i ) no other nodes Qj ∈ Γ, j 6= i, are inside the triangle AQ i B; ii ) no segments Qj−1 Qj and Qj Qj+1 ∈ Γ intersect either AQi and BQi . Hence priority is given to the nodes belonging to Γ instead of point N. If the connection is possible go to d.13); d.11) if no connections are possible ∀Qi ∈ C(Qi ), or C(Qi ) = ∅, check if the connection
of N with AB is possible, according to the conditions set in d.10). If so go to d.13); d.12) if no nodes are satisfying d.10) nor d.11):
d.12.1) loop over Qi ∈ Γ: provided the connection of Qi with AB is possible, save the 30 nodes which give the highest α ˆ i , being α ˆ i > 0 (Figure 5.6). Order
them from the one that gives the maximum α ˆ i to the minimum α ˆ i . Call this set of nodes Φ(Qα ); \ d.12.2) loop over Φ(Qα ): select the node with the minimum area of AQ α B. Then ¯ ¯ \ go to d.13). If the minimum area of AQ α B > A, being A a default value, go to d.14); d.13) create a new element. Update Γ by deleting AB and inserting the other two sides of the new element; d.14) select the next segment according to the chosen AFM implementation. Label again this segment AB; d.15) go back to d.3).
e) Mesh quality enhancement options After Ω has been discretized, MESHGR provides two post-processing routines for the enhancement of the mesh quality: the diagonal swapping and the Laplacian smoothing,
65
5. UNSTRUCTURED MESH GENERATORS
Figure 5.5: Identification of a prospective node N and of its node neighbours on the active front.
Figure 5.6: Determination of the angles αi between AB and the nodes Qi ∈ Γ.
66
5.3 Other mesh generators
Figure 5.7: Effect of Laplacian smoothing.
both acting on the whole set of elements forming the mesh. The Laplacian smoothing improves the element equilateral shape by repositioning a node at the barycenter of the polygon formed by its neighbour elements (Figure 5.7). The diagonal swapping procedure modifies the mesh in order to obtain a number of connections per node as close as possible to an optimal topology: 6 neighbouring elements for interior nodes, 3 for nodes along the domain boundary segment and 1 if the node is at the extremity of a boundary segment. Atkinson et al. (2004) showed the numerical improvements of solving the SWEs on an equilateral, 6 valence, mesh configuration. An example of a mesh topology modified according to the mentioned post-processing routines, represented in Figure 5.8, demonstrates how diagonal swapping can provide an optimal and uniform distribution of node connections, targeted to 6 elements per node. The functionality of the two options is however limited to a topological improvement of the mesh (diagonal swapping) or to repositioning the internal nodes (Laplacian smoothing), while node insertion or deletion operations are not included.
5.3
Other mesh generators
In the following paragraphs a series of alternative mesh generators, some of which widely used in the coastal engineering community, are presented. Where possible, some practical meshing studies are documented in order to retrace the trials carried out when looking for a reference algorithm for the present research.
67
5. UNSTRUCTURED MESH GENERATORS
Figure 5.8: Combined action of diagonal swapping and Laplacian smoothing: (a) an odd number of node connection is transformed into (b) an even one by diagonal swapping and (c) a higher degree of equilaterality is reached with Laplacian smoothing.
68
5.3 Other mesh generators
5.3.1
TRIANGLE
TRIANGLE (Shewchuk, 1996) is a mesh generator available with the BatTri matlab interface (Bilgili et al., 2005) which performs conforming and constrained DTs by means of an iterative node insertion algorithm. The initial input are the bathymetry scatter set and a PSLG (see paragraph 3.3.1), expressed as a sequence of boundary nodes and corresponding segments. Islands can be specified as internal holes in the domain. After the first rough mesh is generated, successive refinements can be performed where, at each step, the user is asked to choose from a series of bathymetry dependent criteria. Diagnostic plots enable to accept the new mesh or to go back to the previous one. The mesh refinement scheme is based on an element-by-element check: if the criterion is not met, the algorithm inserts a new node inside the triangle and re-triangulates the mesh. If the criterion is respected, the element is left unaltered. The available criteria are listed here, where h indicates the average element depth, ∇h the module of
the bathymetric gradient evaluated at the vertices of the element, α is a user defined coefficient and A the area extension which the element can have:
- h/∇h refinement: the maximum area which the element can have is given by the relation: A≤
h ∇h α
(5.8)
- h refinement: the maximum element area is linearly related to a scalar factor: A≤
h α
(5.9)
- 1/∇h refinement: the maximum element area is related to the inverse of the bathymetric gradient. This means that a steeper slope imposes a smaller area than a flat bathymetry: A≤
1 ∇h α
(5.10)
- maximum slope refinement: all elements which exceed a threshold bathymetric gradient, set by the user, are refined. This option is useful to better represent sand waves and other small scale bottom features; - tide wavelength to mesh size ratio refinement: this relationship expresses the CourantFriedrichs-Levy condition (Formula 4.8) in terms of the tidal wavelength to mesh size
69
5. UNSTRUCTURED MESH GENERATORS
ratio R wanted by the modeller. The area that an element can have is given by: A≤
ghT 2 R2
(5.11)
- constant maximum area refinement: a constant bathymetry independent maximum area extension is imposed on all the elements; - delta(h)/h refinement: the following condition is checked: hM AX − hM IN delta(h) ≤α = h hM IN
(5.12)
where hM AX and hM IN are the depths of the deepest and the shallowest nodes of the element. If the condition is not met, the element is refined until its area is less than another user defined area threshold, called tarea. An example of a refinement procedure for successive steps is presented in Figure 5.9, where a simple squared domain with an idealized bathymetry is refined according to the local bathymetry gradient: more than producing a mesh for a tide propagation modelling, this example aims at exploring the possibilities offered by TRIANGLE in constructing well graded triangular meshes. The series of successive refinements of Figure 5.9 demonstrates how the steeper portion of the domain, represented by the darker zone of Figure 5.9(b), can be successfully discretized, until a desired grade of detail. An unstructured mesh suitable for modelling applications should meet a series of geometric properties, among which an area change transition between neighbour elements not exceeding a reference value (i.e. 0.5). Starting from the mesh of Figure 5.9(f), further discretizations are imposed in order to achieve this requirement. This time the tide wavelength to mesh size criterion is applied, for a hypothetical tide wave propagating over the domain. The meshes resulting from the second series of refinements (Figure 5.10) better resolve the shallower part of the domain, with acceptable element area transitions.
5.3.2
Application of TRIANGLE to a real case study
The unstructured mesh used for the ADCIRC storm simulations reported in Ara´ ujo et al. (2011a) was created with TRIANGLE and its interface BatTri. The chosen
70
5.3 Other mesh generators
(a)
(b)
(c)
(d)
(e)
(f )
Figure 5.9: TRIANGLE: successive refinement procedure: (a) domain bathymetry; (b) spatial extent of the inverse of the bathymetry gradient; (c) “first cut” mesh and (d) to (f) successive refinements.
71
5. UNSTRUCTURED MESH GENERATORS
(a)
(b)
Figure 5.10: TRIANGLE: new series of refinements based on the wavelength to mesh size criterion.
domain extends from 39.4◦ to 43.1◦ N and 8.67◦ to 16.19◦ W (Figure 5.11), in order to place the open ocean boundary in deep water and far enough from the coast to avoid potential interference issues (Blain et al., 1994). The coastline was located at a depth of 4 m below MSL. The DEM used for the mesh generation had different bathymetry sources: at the Portuguese/Spanish shelf, hydrographic charts (resolutions of 1:150 000 and 1:1000 000) published by the Portuguese Navy (Instituto Hidrografico) were used. The contours of these charts were digitized from the coastline up to depths of 4000 m. The bathymetry for deeper waters was obtained from the World database of the Institute of Geophysics and Planetary Physics (Smith and Sandwell, 1997), with a grid resolution of 1 min. The meshing methodology consisted in a series of successive refinements of an initial coarse mesh, applying each time the h/A ratio criterion (formula 5.9) for a specific range of depths. Following Aretxabaleta (2003), the domain was split into three zones: the first one, with depths ranging from −12 m to
−180 m (BatTri considers negative depths for the seabed), the second, encompassing the shelf break and shelf slope, until a depth of −2100 m, and the third zone from −2100 m to the exterior ocean boundary. A “first cut” mesh was then created by
imposing a maximum element size of 2000 m, 3000 m and 15000 m for the three zones respectively (Figure 5.12). At this point a series of “try and error” refinements was performed: at each step, a user defined α coefficient is entered and refinement is applied to all elements not satisfying relation (5.9). The constraint may be applied to a specific
72
5.3 Other mesh generators
Figure 5.11: Computational domain, extending for 400 km along the west Iberian Coast and 615 km across the Atlantic Ocean.
Figure 5.12: Example of “first cut” mesh for the domain of Figure 5.11.
73
5. UNSTRUCTURED MESH GENERATORS
Figure 5.13: Unstructured triangular mesh with 41367 nodes and 80590 elements, used in Ara´ ujo et al. (2011a).
bathymetric range or to the whole domain. The repeated application of (5.9) produced the highest refinement near the coast and in the presence of offshore sand banks (e.g. Galicia seamount). This refinement criterion was chosen because it is related to the wavelength to mesh size ratio. In fact, considering that λ=T and that ∆x ≈
p gh
s√
3 A 4
(5.13)
(5.14)
being T the wave period and ∆x the side length of an equilateral triangle with the same area A of the mesh triangle, it results that: λ 2 h ∝( ) A ∆x
(5.15)
The final mesh (Figure 5.13) has 41367 nodes and 80590 elements, with a spatial resolution ranging from approximately 17 km in the open ocean to 250 m near the coast. With this depth-dependent refinement procedure, the final mesh was able to reproduce
74
5.3 Other mesh generators
Figure 5.14: Comparison between the bathymetric isolines of the mesh of Figure 5.13 (darker lines) and of the DEM (brighter lines) for the Galicia seamount and for the shelf break, where the major bathymetry gradients are found.
with a good approximation the isolines of the nautical charts. Attention was taken in assuring a smooth element grading from deep to shallow waters and in enhancing the resolution of the shelf break and continental slope, where higher bathymetry gradients are located. Laplacian smoothing was applied to the final mesh to achieve a higher element quality index. A visual check of the adequacy of the mesh resolution can be done by plotting the depth contours of the mesh against the isolines of the DEM (Figure 5.14). For the shown region, which encompasses a complete cross section of the shelf break, the concordance is high; some visible discrepancies appear just when the bathymetric isolines present sudden changes. BatTri also provides an option to visualize the computational time step required by the Courant-Lewy-Friedrichs condition for the present mesh, and the corresponding Courant number Cr . The required time step is shorter in the shallower domain, where
75
5. UNSTRUCTURED MESH GENERATORS
Figure 5.15: Time step requirement for a Cr = 1 for the mesh of Figure 5.13, visualized near the coast, where smaller elements are generated. The minimum time step to satisfy the constraint is around 500 s. The bar scale is expressed in seconds.
the element size decreases substantially. The minimum time step which may be used is around 500 s for a Cr = 1; however, in order to prevent any model instability, a stricter Cr should be considered (Blain et al., 2010). In the related ADCIRC simulations, a time step of 30 s was chosen.
5.3.3
The SMS mesh generator
The Surface-water Modeling System SMS is a commercial analysis tool interfacing a wide range of numerical models for supporting hydrodynamic computations. It is provided with a mesh generator which makes available the following series of discretization methods: - patching: after a four sided domain is defined, together with a node distribution along the domain boundaries, the location of the interior nodes is interpolated from the node location along opposite sides. It can be classified as an algebraic interpolation method, recommended for domain with a preferential flow direction (channels, rivers, etc.); - triangulation: the nodes of a scatter set are saved as nodes of the prospective mesh, while the elements are generated according to a DT;
76
5.3 Other mesh generators
(a)
(b)
Figure 5.16: Subdivision of a river domain into subregions: (a) one subregion is selected and attributed to a specific meshing method; (b) resulting composite mesh, including both triangles and rectangles. The used meshing algorithm is the boundary spaced paving for main river course, the boundary spaced adaptive tessellation for the right and left banks.
- adaptive tessellation: this method overlays a quadtree structure on the domain, and splits recursively the quad until the local NSF is met. Groups of interior segments and points can be forced into the final mesh. The NSF can be determined in two ways: in the boundary spaced adaptive tessellation the element size is a function of the varying node spacing along the domain perimeter; in the scalar adaptive tessellation, the NSF is calculated from the bathymetry scatter set. As sharp discontinuities of the element size may arise when passing form a quadtree level to the upper/lower, a transition parameter called bias can be adjusted by the user; - paving: this option is based on the AFM: similarly to the adaptive tessellation, the boundary spacing and the scalar paving density modes are available. Triangular and/or quadrilateral elements can be formed according to the numerical model in use (see GFGEN, 2001). The listed methods can also be combined together if a domain is subdivided into subregions, and a different method is attributed to each of them. An example of a combination of different mesh types is applied to a meandering river course in Figure 5.16. SMS also provides the LTEA method (Hagen et al., 2000; Hagen et al., 2002; Hagen et al., 2006) and the LTEA+CD method (Parrish and Hagen, 2009), which are currently available for the ADCIRC module. These are two versions of an a-posteriori adaptive criterion, which relocates the nodes of an initial mesh according
77
5. UNSTRUCTURED MESH GENERATORS
Figure 5.17: Extension of the new ocean domain (outer polygon), compared with the mesh used in the work of Ara´ ujo et al. (2011a).
to the estimation of an error associated with a the harmonic analysis of a M2 tide solution. The methodology will be tested in chapter 7 in order to assess and integrate a series of a-priori meshing criteria. SMS, which is a commercial software, provides several mesh visualization and post-processing tools, as spring relaxation (another term for Laplacian smoothing used inside SMS), local and global mesh refinement, mesh editing and quality checking. SMS has been extensively used as a support for the input pre-processing and the visualization of the MESHGR output (see the flowchart 5.1).
5.3.4
Application of the SMS paving algorithm to an oceanic domain
SMS provides a robust and flexible meshing tool applicable to large ocean domains, as for the case presented hereafter. The discretized domain covers a significant part of the North Atlantic Ocean, including the Mid Atlantic Ridge (Figure 5.17). The shallow water wavelength to mesh size ratio criterion was considered for the meshing (formula 4.1). The period T was set equal to 12h 25min, corresponding to the M2 tidal constituent. A NSF equivalent to a wavelength to mesh size ratio R = 1200 was
78
5.3 Other mesh generators
Figure 5.18: 2300 Km x 1350 Km unstructured mesh generated with the SMS scalar paving method for the Atlantic domain of Figure 5.17.
generated for the bathymetry range 4 m - 2500 m, assuring a fair discretization of the continental shelf and shelf slope, while a R = 100 condition was imposed for deeper bathymetries. The smooth junction of the two functions was assured by a maximum variation of the adjacent element area of a factor of 0.6. The scalar paving algorithm was chosen to perform the meshing, The final mesh resulted in 60766 nodes and 118656 elements, with an element resolution ranging from 80 Km along the open ocean to 220 m along the Portuguese coastline. The generated mesh (Figure 5.18) has been later tested in Marujo da Silva (2011).
5.3.5
AAMGM
The Automatic Adaptive Mesh Generation (AAMGM) program (Lee, 1999; Lee and Hobbs, 1999) implements an AFM algorithm, where a NSF, similarly to MESHGR, is interpolated from a background mesh. The order in which the mesh nodes are created follows a hierarchical generation principle (Figure 5.19): first the endpoints at the extremities of the boundary curves are saved as mesh nodes, then the nodes interpolated along the boundary curves, and at last the points generated inside the domain are added to the prospective mesh. AAMGM supports the creation of isotropic and anisotropic, triangular and quadrilateral elements, and is provided of the quality enhancement procedures of Laplacian smoothing and diagonal swapping. AAMGM was taken into consideration as a possible alternative to MESHGR, being the similarities
79
5. UNSTRUCTURED MESH GENERATORS
Figure 5.19: Hierarchical generation principle in AAMGM: (a) the extremities of the boundary curves delimiting the domain (and its subregions if present) are saved as mesh nodes; the mesh is then completed with the (b) nodes generated along the boundaries and (c) the interior nodes. A NSF defines the mesh element size. From Lee and Hobbs (1999).
80
5.3 Other mesh generators
between the two codes numerous. Moreover, AAMGM is provided of a cycle of refinements, where the final mesh at one level is used as the background mesh of the successive level. In practice the convergence of the final mesh to the NSF values is obtained within 3 cycles. This internal loop may overcome one of the major limitations of MESHGR, where the representation of the final NSF heavily depends on the design of the initial background mesh. On the other side the input NSF is expressed in analytical terms and not by piecewise functions depending on the bathymetry or other model variables, representing a serious drawback for applications to real domain. Due to the impossibility of performing real case tests, AAMGM was not considered appropriate for discretizing ocean and coastal domains.
5.3.6
TRIGRID
TRIGRID (Henry and Walters, 1993) is an open source 2D unstructured mesh generation package (http://trigrid.sourceforge.net/), targeted to finite element applications in the coastal environment. It performs constrained DT, with equilateral triangles as target elements. The triangle area is related to a scalar density function, which is proportional to specific flow, geometry or bathymetry characteristics. This algorithm has been considered because of the novel cluster-based method applied for the subdivision of the computational domain, depending on a predefined NSF. The algorithm is implemented as follows: - after a scalar NSF is set, the domain is covered by a Cartesian grid, where the element size is half the smallest triangle side of the final model mesh (although in the referred article it is not explained how to set this length); - cluster elements are formed by incorporating the squares of the Cartesian grid (Figure 5.20), starting from the square with the deepest depth. The cluster shape is controlled in order to be not distorted. The addition of single elements to one cluster stops when the cluster area/node function meets a predefined limit value. Some inter-spaces between neighbour clusters are admitted; - the centroids of the clusters form a set of nodes suitable to be the mesh nodes. A DT is then performed to form the mesh;
81
5. UNSTRUCTURED MESH GENERATORS
Figure 5.20: TRIGRID clusters: (a) initial formation of a cluster; (b) formation of a cluster in presence of already existing clusters; (c) badly shaped cluster. The underlying grid is the Cartesian grid whose elements have half size of the smallest triangle in the final model mesh (from Henry and Walters, 1993).
- the quality of the mesh can be enhanced with several post-processing operations (mesh splitting or mesh joining; bathymetry re-interpolation; element shape adjustment with smoothing algorithms, like the Laplacian relaxation, etc..). The code accepts a scalar function as a criterion for driving the node generation: the authors propose that node distance is proportional to the inverse of the water depth, in order to maintain a constant resolution over a wavelength in shallow waters: ρ∝
1 h
(5.16)
where ρ is the target NSF and h the bathymetric depth. TRIGRID has been applied to support the meshing of wide coastal domains (Walters et al., 2001). The package is downloadable at http://sourceforge.net/projects/trigrid/, and provided with a graphicbased editor, but the example tests could not be entirely run, maybe due to some missing code subroutines. The available contacts at the page seem not to be active any more. For these reasons the software was dropped off from further research.
5.3.7
Comparative analysis between mesh generators
A conclusive paragraph is dedicated to the comparison of the characteristics of the aforementioned mesh generators (Table 5.1). An absolute criterion for the determination of a “best” mesh generator does not exist, as each code has been written for applications of a specific modelling field, in which these codes are documented to be successful. In the present research an open source discretization tool able to implement the NSF approach for finite element ocean and coastal models was needed. Delaunay triangulators were excluded because the meshing of the domain often passed through a
82
5.3 Other mesh generators
Table 5.1: Some general properties of the tested mesh generators
series of successive refinement steps and not as a only step discretization, according to a unique final NSF reflecting multiple criteria. At the end the choice fell on MESHGR, which, among the tested and reviewed algorithms, is an AFM generator which meet the mentioned requirements. The code, which had been applied only to simple domain geometries, manually defined NSF and generated meshes with a limited number of elements, was considered a good opportunity for development and applications to more demanding conditions.
83
6
THE MULTI-CRITERIA MESHING METHOD. APPLICATION TO IDEALIZED GEOMETRIES 6.1
Introducton
A critical issue in the current modelling practice is still represented by the definition of the efficient mesh node placement, in order to assure the convergence of the model results without overburdening the computation of an unnecessary degree of refinement. A useful tool for addressing the control of the element characteristics in a mesh generator algorithm is the NSF, introduced in chapter 3. Isotropic meshes derive from scalar NSFs, while anisotropic meshes need to be defined by means of vector NSFs, where, apart from the element size, other fields, i.e. the element skewness and the direction of elongation, must be specified. The function may be determined from a-priori or a-posteriori criteria. A-priori criteria derive from known physical and hydrodynamic factors of the relevant processes being studied, such as the bathymetry, the topographic length scale, and expected characteristics of the flow motion. A-posteriori criteria are used in adaptive meshing techniques (Borouchaki et al., 1997a, 1997b; Hagen et al., 2000; Bacopoulos et al., 2011) where a “first try” mesh is modified in order to minimize or drive to an uniform value an error estimation of the flow problem solution.
85
6. THE MULTI-CRITERIA MESHING METHOD. APPLICATION TO IDEALIZED GEOMETRIES
The recent improvements of the model capabilities in the prediction and hindcast of complex natural phenomena have been achieved, among other factors, by including bigger spatial and temporal scales, by increasing the spatial resolution of the input datasets and by considering all the relevant physical processes involved. While in the past the model input parameters and boundary conditions were used as tuning tools, which made the model case dependent, the today mainstream approach relies on the underlying physics of the studied process: all relevant forcing mechanisms contributing to the event must be correctly reproduced. The evolution of the storm surge modelling, where different computational models are coupled together, clearly follows this trend (Cialone et al., 2002; Ebersole et al., 2010; Dietrich et al., 2011; Honghai et al., 2013). The multi-criteria meshing method (Mazzolari et al., 2012) formally translates this approach to the task of mesh generation: it requires that a series of a-priori and/or a-posteriori, scalar and/or vector meshing criteria is identified (from here the term multi-criteria comes). Each criterion is expressed in terms of a NSF, covering the whole domain or part of it. The final NSF driving the meshing process is obtained by merging the requirements of the previously defined constraints in a unique piecewise linear function. In the process of exploring the capacities (and limits) presented by MESHGR, a sequence of applications of the multi-criteria method combined with the AFM are developed here for idealized bathymetries, following an increasing level of complexity. In each example the bathymetry and the coastline are digitized in a way that reproduces features often encountered in real domains. Simply connected and multiply connected 2D regions are first considered for isotropic mesh generation; at a second stage anisotropic conditions are imposed as well. The main purpose of these tests is to evaluate the capacity and robustness of MESHGR in generating good quality meshes according to target NSFs derived from multiple criteria.
6.2
Test 1: rectangular domain with a steep bottom
In order to show how the multi-criteria method operates, two a-priori criteria have been applied to a rectangular domain. The domain bathymetry has isobaths fairly parallel to the coastline, with the presence of a steep feature representing an idealized
86
6.2 Test 1: rectangular domain with a steep bottom
Figure 6.1: Test 1. Domain extension and bathymetric isolines.
Figure 6.2: Test 1. Background mesh. The superimposed rectangle represents the domain boundaries.
continental shelf slope (Figure 6.1). The background mesh of Figure 6.2, formed by 46 nodes and 67 elements, will be used for the definition of the NSFs. The first meshing criterion is the well-known wavelength to mesh size ratio: T λ = ρ11 (x, y) = R
p gh(x, y) 200
(6.1)
where ρ11 is the target element size for test 1 and criterion 1, λ the shallow water tide wavelength, T is the period of the M2 semidiurnal tide constituent, considered equal to 12h 25min, and h the bottom depth. A R equal 200 has been chosen as a reference value according to Westerink et al. (1995). For (6.1), MESHGR produces, for a
87
6. THE MULTI-CRITERIA MESHING METHOD. APPLICATION TO IDEALIZED GEOMETRIES
Figure 6.3: Test 1. Mt11 generated by MESHGR before applying post-processing operations.
Figure 6.4: Test 1. Histograms of the shape quality index qi for Mt11. Only the range of values [0.5-1] is represented.
front searching implementation of the advancing front, the mesh shown in Figure 6.3. This mesh is named Mt11, because it belongs to test 1 and is derived from the first chosen criterion. From here on, all test meshes will be named after this convention. Mt11 has 628 nodes and 1122 elements, and has not been submitted yet to the quality improvement subroutines provided by MESHGR, which include Laplacian smoothing and diagonal swapping (see paragraph 5.2.2). A MATLAB subroutine (Hill, 2007) is used to visualize, in form of histograms, the shape quality index qi , defined by formula (3.13). Figure 6.4 shows that almost 500 elements over 1122 are equilateral and many others are almost attaining an equilateral shape, while “bad” shaped elements, where “bad” means with a qi below 0.6, are 6.
88
6.2 Test 1: rectangular domain with a steep bottom
Figure 6.5: Test 1. Mt11 before post-processing, evaluated against the criteria of maximum area change, maximum slope and node valence higher than 8. None of the elements has a slope higher than 0.1.
Mt11 has been subsequently post-processed with the diagonal swapping and the Laplacian smoothing routines. This makes it possible to compare Mt11 before and after post processing against the quality criteria of element area change, maximum slope and node valence. The element area change evidences neighbour elements with a difference in the area extension of more than 50% between them. The maximum slope indicator calculates, after the mesh is interpolated to the bathymetry, the maximum inclination among the sides of an element, defined as the ratio between the depths of the two extremities and their horizontal distance. It evidences all elements with at least one side exceeding a threshold value, set in this case to 0.1. It is known to modellers that rapid changes in slope may origin computational errors (Hannah and Wright, 1995). Nodes with a valence higher than 8, where the term valence indicates the number of edges attached to a node, are evidenced if present. The comparison of Figure 6.5 and Figure 6.6 reveals the usefulness of the quality enhancement process: nodes with a valence higher than 8 disappear, while the number of adjacent elements with an excessive area change is consistently reduced. This result is achieved despite the high number of nodes and sides, with respect to the total, being constrained to the domain boundaries and therefore fixed. For the present bathymetry no issues are related to the mesh slope. A peculiarity of the unprocessed mesh is that all nodes with a valence higher than 8 are connected with elements which exceeds the area change limit of 50%.
89
6. THE MULTI-CRITERIA MESHING METHOD. APPLICATION TO IDEALIZED GEOMETRIES
Figure 6.6: Test 1. Mt11 after diagonal swap and Laplacian smoothing, evaluated against the criteria of area change, maximum slope and node valence higher than 8. None of the elements has a slope higher than 0.1.
The second meshing criterion is obtained as the inverse of the bathymetric gradient, according to: ρ12 (x, y) = c
1 ~ |∇h|
(6.2)
~ is the bathymetric gradient calculated in SMS where c is a coefficient set to 2000 and ∇h at each mesh node by averaging the normals of all the planar faces of the triangulated bathymetry attached to that node. In the space the normal vector to a surface, often simply called the “normal”, is a vector perpendicular to the surface itself. If the surface is given as a scalar function F = F (x, y, z), defined over a domain Ω ⊆ R3 , then the normal ~nP to a point P ∈ Ω is expressed as: ˆ P ~ P = ( ∂F ˆi + ∂F ˆj + ∂F k) ~nP = ∇F ∂x ∂y ∂z
(6.3)
~ P denotes the gradient of F and ˆi, ˆj and kˆ the unit vectors along the x, y where ∇F and z directions. In SMS, which is used as support interface, the bathymetry surface is usually represented as a linearly triangulated regular or irregular network, which can be written as: h(x, y) = ax + by + d
90
(6.4)
6.2 Test 1: rectangular domain with a steep bottom
Figure 6.7: Test 1. Mt12 before applying post-processing operations.
Figure 6.8: Test 1. Mt12 quality before diagonal swapping and Laplacian smoothing, evaluated against the criteria of area change, maximum slope and node valence.
where the depth h is interpolated from the nodes of the bathymetry scatter set. Assuming that P coincides with a node, the normal to P is given by: ~ P = (aˆi + bˆj)P ~nP = ∇h
(6.5)
~ P is perpendicular to the bathymetric isolines. It has to be reminded that ∇h
The mesh obtained from criterion (6.2) with the front searching option is called Mt12, and has 942 nodes and 1813 elements (Figure 6.7). The following figures show the distribution of some attributes of Mt12 (area change, maximum slope and node valence) before (Figure 6.8) and after (Figure 6.9) post-processing.
Their comparison shows
again the good performance of the MESHGR post-processing options, where in the final mesh all nodes have a valence not higher than 8, and the number of abrupt area
91
6. THE MULTI-CRITERIA MESHING METHOD. APPLICATION TO IDEALIZED GEOMETRIES
Figure 6.9: Test 1. Mt12 quality after diagonal swapping and Laplacian smoothing, evaluated against the criteria of area change, maximum slope and node valence higher than 8.
transitions decrease from 88 to 5. The new node locations in the post-processed mesh require re-interpolating the bathymetry in order to check the maximum element slope, but, as for the unprocessed mesh, no slope warnings are detected. The application of the multi-criteria meshing method calculates the final NSF as: ρ13 = min{ρ11 (x, y); ρ12 (x, y)}
(6.6)
which generates a more stringent requirement than the application of the single criteria separately. Figure 6.10 illustrates how the multi-criteria approach operates in the definition of the final element target size from the single NSFs: the method, applied here just to two criteria, can be extended to as many as needed. It has to be outlined that the interpolation of the minimum size function may result in a final value which is different than every single NSF, as a consequence of interpolating piecewise functions over a coarse background mesh. In Figure 6.10 this is evident for the interval 20 km to 45 km along the horizontal axis. A smaller interpolation error would be obtained if the spacing of the background mesh is reduced. The discretization derived from ρ13 produces, adopting the front searching option as done before, Mt13 (Figure 6.11), with 1624 nodes and 3087 elements, while for ρ11 the mesh had 628 nodes and 1122 elements, and for ρ12 942 nodes and 1813 elements. Also for Mt13 the diagonal swapping and Laplacian smoothing procedures are applied (Figure 6.12) in order to produce a better graded mesh.
92
6.2 Test 1: rectangular domain with a steep bottom
Figure 6.10: Test 1. Multi-criteria method: (a) depth profile and (b) magnitude of the single NSFs defined over the back ground mesh, for a domain cross section profile.
93
6. THE MULTI-CRITERIA MESHING METHOD. APPLICATION TO IDEALIZED GEOMETRIES
Figure 6.11: processing.
Test 1.
Mt13 before diagonal swap and Laplacian smoothing post-
Figure 6.12: Test 1. Mt13 after diagonal swap and Laplacian smoothing post-processing.
Figure 6.13 represents a zoning of the domain based on the leading criterion: the shaded part corresponds to the elements whose size is derived from ρ12 , the rest being derived from ρ11 : the steep slope, located horizontally in the center of the domain, falls in the zone meshed by the NSF dependent on the bathymetric gradient. The problem of under-refinement due to a constant wavelength to mesh size ratio function is resolved here, as the level of discretization is made dependent on the bathymetric gradient and proportional to a user defined coefficient. The qi distribution is represented before (Figure 6.14a) and after (Figure 6.14b) the post-processing, which evidences the effects of Laplacian smoothing and diagonal swapping on the element shape: bad shaped elements (qi below 0.6) are reduced from 15 to 1, while, although the absolute number of elements in the highest quality class decreases,
94
6.2 Test 1: rectangular domain with a steep bottom
Figure 6.13: Test 1. Domain zoning for Mt13 showing the spatial influence of each meshing criterion. The evidenced part corresponds to the area where the final density is determined by ρ12 . The meshing of the complementary area is driven by ρ11 .
the dimension of the histograms located in the interval 0.95 < qi ≤ 1 tends to increase in the post-processed mesh. After exporting Mt13 into SMS, it is possible to evaluate the maximum area change, maximum node connection and element slope (Figure 6.15 and Figure 6.16). As already happened for Mt11 and Mt12, the post-processing operations increase the mesh quality, by reducing the element couples with an area change over 50% from 88 to 8 and lowering the maximum node valence from 9 to 8. No slopes steeper than 0.1 are registered in both the unprocessed and processed Mt13. Table 6.1 reports an overview of the discussed quality parameters of the meshes for test 1, underlying again the effect of post-processing operations.
Table 6.1: Quality measures for the meshes of test 1, before and after post-processing.
95
6. THE MULTI-CRITERIA MESHING METHOD. APPLICATION TO IDEALIZED GEOMETRIES
(a)
(b)
Figure 6.14: Test 1. Mt13. Histograms of the quality index qi : (a) before post-processing, (b) after post-processing. Only the range of values 0.55-1 is represented.
6.3
Test 2: squared domain with 2 islands
In this example MESHGR is tested for a squared domain of 100 km side length, with the presence of two islands which are not included in the discretization (Figure 6.17). The domain is multiply connected. The depth contours and the extension of the domain to be meshed have been digitized manually. From one side the bathymetry scatter set reproduces the two islands on a gently inclined sea floor, and on the other side it has been kept as simple as possible. The background mesh drawn in Figure 6.18 has nodes which coincide, for simplicity, with the points of the bathymetry scatter set: it follows that the point distribution is denser in correspondence of the curvilinear internal boundaries. This observation may be extended to more complex cases: in general the refinement of the background mesh should be able to reflect, as far as possible, the spatial characteristics of the domain geometry (coastline, islands, bathymetry...) or of a physical variable. If this prerequisite is not met, local discrepancies between original
96
6.3 Test 2: squared domain with 2 islands
Figure 6.15: Test 1. Mt13 quality before diagonal swap and Laplacian smoothing, evaluated against the criteria of area change and valence higher than 8. No slope gradients higher that 0.1 are detected.
Figure 6.16: Test 1. Mt13 quality after diagonal swap and Laplacian smoothing, evaluated against the criteria of area change and node valence higher than 8. No slope gradients higher that 0.1 are detected.
and interpolated features may arise. As for test 1, the target mesh element size is obtained by applying a series of intermediate criteria, derived from geometric and bathymetric functions. A first criterion, related to the domain geometry, calculates the distance from the nearest island: ρ21 (x, y) = min{disland1 (x, y); disland2 (x, y)}
(6.7)
where dislandi (x, y) is the distance of a point of coordinate (x,y) from the centre of the ith island. The magnitude of ρ21 is shown in Figure 6.19.
The application of ρ21 as
input to MESHGR gives, for the circulating front option, Mt21, represented in Figure
97
6. THE MULTI-CRITERIA MESHING METHOD. APPLICATION TO IDEALIZED GEOMETRIES
Figure 6.17: Test 2. Layout of the domain and bathymetric contours. The coastline is represented by the last contour on the right.
Figure 6.18: Test 2. Background mesh, containing 37 nodes and 60 elements, superimposed to the domain boundaries.
98
6.3 Test 2: squared domain with 2 islands
Figure 6.19: Test 2. Magnitude of the NSF ρ21 , determined as the distance from the nearest island.
Figure 6.20: Test 2. Mt21 obtained by applying ρ21 , after post-processing operations: elements which exceed the maximum element area change condition of 50% and nodes with a valence higher than 7 are evidenced.
99
6. THE MULTI-CRITERIA MESHING METHOD. APPLICATION TO IDEALIZED GEOMETRIES
6.20. The front circulating implementation of the AFM is chosen here for testing the quality of the output meshes in a similar way to what done for the searching front option. Mt21, obtained after a diagonal swap and a Laplacian smoothing processing, has 3582 nodes and 6890 elements, 3 of which do not respect the area change transition limitation of 50%, and are all located along the domain boundaries. The element side lengths go from 1000 m along the island boundaries to 8979 m on the top left and bottom right domain corners. It is clear that a higher grade of refinement is reached in the vicinity of the islands, while the shallower waters near the coastline still lack a good discretization. For this reason the wavelength to mesh size criterion (formula 6.1) is selected as an extra refinement solution, called here ρ22 . The final NSF, ρ23 , is eventually derived as the minimum between ρ21 and ρ22 (Figure 6.21).
Figure 6.21: Test 2. NSF ρ23 .
The resulting mesh, Mt23 (Figure 6.22), obtained by selecting the circulating front option, has 4724 nodes, 9105 elements, and a node spacing ranging from 810 m to 5811 m. It can be appreciated the additional refinement brought to the shallowest zone (right part of the domain). A few elements show an area change transition of over 50%,
100
6.4 Introduction to anisotropic meshing
Figure 6.22: Test 2. Mt23, obtained by applying ρ23 , after post-processing operations. Neighbour elements with an area transition higher that 50% and nodes with a valence higher that 7 are evidenced.
while 13 nodes have a valence equal or higher than 7. Figure 6.23 evidences the zone where the wavelength to mesh size criterion drives the meshing process. The rest of the domain is meshed according to the geometry derived criterion. Figure 6.24 plots the distribution of qi for Mt23, confirming the capacity of MESHGR to produce meshes with an overall high grade of well-shaped elements. As a last quality check, the mesh representation of the underlying bathymetry is visualized in Figure 6.25. Although this type of test is visual and therefore only qualitative, a good indication of the exact reproduction of the bathymetry is given. Here the agreement between the mesh and the scatter set contours is optimal along straight segments, while the only detectable differences are localized in the vicinity of singular points.
6.4
Introduction to anisotropic meshing
Unstructured meshes may include anisotropic characteristics, required by physical and/or numerical conditions. Anisotropic applications are usually related to flow problems with strong directional patterns: boundary layers, shear layers, shock-fronts, airfoils
101
6. THE MULTI-CRITERIA MESHING METHOD. APPLICATION TO IDEALIZED GEOMETRIES
Figure 6.23: Test 2. Domain zoning for Mt23 showing the area where the node placement is driven by the wavelength to mesh size ratio function.
Figure 6.24: Test 2. Histogram of the quality index qi for Mt23. Only the range of values [0.55-1] is represented.
102
6.4 Introduction to anisotropic meshing
Figure 6.25: Test 2. Comparison between the Mt23 bathymetry and the scatter set bathymetry of Figure 6.17. The mesh contours are represented with a darker line.
and other geometries subject to high Mach and Reynolds numbers. The creation of anisotropic elements in 2D domains requires the definition of extra parameters in addition to the element size, as the element shape and the direction of elongation. The triangulation methods already reviewed for isotropic meshes (Chapter 3) are applied as well for anisotropic meshes: the DT (Vallet et al., 1991; Bossen and Heckbert, 1996; Borouchaki et al., 1997a, 1997b; Castro-Diaz and Hecht, 1995); a modified version of the AFM, where a layer of elements is generated from the active front (Pirzadeh, 1994; Hassan et al., 1995; Garimella and Shephard, 2000), or hybrid techniques (Cougny and Shepard, 1999). Also for the anisotropic meshing these algorithms provide a set of options for local topology modification: the most common are node insertion/deletion, node movement, edge swapping and area smooth grading. Central to the insertion of anisotropic conditions is the use of a metric tensor (see section 3.4.4c), derived from a-priori or a-posteriori criteria: the metric introduces a mapping of the physical space into a control space, where the target elements have ideal isotropic shapes, as equilateral triangles or squares of edge length equal to unity. Often anisotropic and isotropic requirements coexist in domains with different spatial scale regions, which require mixed meshes (Jansen et al., 2001). A-priori anisotropic discretizations have been applied, in oceanic domains, to steep bottom features, as the shelf break and seamounts, or tide ellipses (Legrand et al.,
103
6. THE MULTI-CRITERIA MESHING METHOD. APPLICATION TO IDEALIZED GEOMETRIES
2007). Applications of anisotropic meshes to hydrodynamic models are mainly associated to a-posteriori adaptive techniques, where the re-meshing is targeted to capture the local anisotropic flow characteristics and at the same time to optimize the computational costs (Remacle et al., 2006; Piggott et al., 2009; Porta et al., 2012). Cobby et al. (2003) modelled river floods with stretched elements along the river bed, in order to align the element sides with the channel banks. Hasan et al. (2012) demonstrated that quadratic elements aligned with the bathymetric contours improve the predicted tides in an estuarine domain. As a drawback, the use of elongated triangles may be associated with rapidly varying element size, which results in higher order truncation error terms (Hagen et al., 2000; 2001). The anisotropic tests presented in this section are meant to show how anisotropic conditions can be included in the multi-criteria method and be inserted into MESHGR. MESHGR documented applications have dealt so far only with isotropic meshes, while its performance for anisotropic discretizations is unexplored. Moreover, elongated elements orientated along the bathymetric isolines may prove to provide a more accurate and realistic representation of the flood propagation and retreat. In order to validate this hypothesis a wet and dry tide study test is performed for an idealized geometry with ADCIRC. The numerical accuracy and computational load will be evaluated as well. This investigation may contribute to give insight into a not yet well explored area of applied modelling.
6.5
Extension of the multi-criteria method to anisotropic meshes
The introduction of an anisotropic requirement to the mesh generation process complies with the conceptual framework of the multi-criteria meshing strategy: the element anisotropy is introduced as an extra discretization criterion, which may be superimposed to the existing isotropic criteria. In MESHGR the anisotropic information is defined by 3 attributes for each background mesh node i : the anisotropy coefficient si , which controls the element shape, and the 2 components of the anisotropy vector, ui and vi , setting the direction of elongation. The smaller the anisotropy coefficient is,
104
6.5 Extension of the multi-criteria method to anisotropic meshes
the more the elements are stretched. In order to show how MESHGR performs, a rectangular domain is meshed with the front searching method, for a constant anisotropic condition si = 0.2 and a constant ρ = 50 m. (Figure 6.26). The direction of elongation is set parallel to the horizontal axis (ui = 1 and vi = 0). Figure 6.26 demonstrates that
(a)
(b)
Figure 6.26: (a) Anisotropic mesh defined by a constant anisotropy coefficient s = 0.2 and a constant target element side length equal to 50 m; (b) particular of the mesh showing the element dimensions. The actual anisotropic ratio of the element is 0.23.
the actual stretching coefficient may slightly differ from the target one (in the close up view the actual si is 0.23 against a target value of 0.2). The discrepancy is even more evident for the elements generated along the boundary, which, according to the front searching method, are most probably the last ones to be generated and hence have to fill the gap between the already meshed area and the boundary.
105
6. THE MULTI-CRITERIA MESHING METHOD. APPLICATION TO IDEALIZED GEOMETRIES
Figure 6.27: Test 3. Domain extension and bathymetric contours.
6.6
Test 3: anisotropic meshing
This test case illustrates how the anisotropic coefficient and the components of the element orientation vector are calculated. The main aim is the generation of elongated elements which adjust to the orientation of the bathymetric isolines. Elements oriented in such way may enhance the realistic rendering of the wetted front in inundation studies. The bottom is represented by an idealized bathymetry (Figure 6.27) whose depths range from -4 m to 15 m (sea depths are considered positive downwards), with the presence of two emerged features, one on the left and one on the right side of the domain. This example reproduces in a simplified way a coastal geometry where the wet and dry process involves a significant domain extension. A background mesh (Figure 6.28) has been digitized manually with 59 nodes and 110 elements. The background nodes, which in this case do not correspond to the points of the bathymetry scatter set, are denser where the curvature of the bathymetric isolines is higher. The NSF for the element size has been set constant and equal to 10 m. The anisotropic attributes si , ui and vi need to be calculated for every ith background node. si is set proportional ~ −1 . ∇h, ~ calculated from equation (6.5), is represented in Figure 6.29, and to |∇h|
shows, as expected, a direction perpendicular to the bathymetric isolines. In order to ~ −1 is projected linearly between avoid too skewed shapes, the range of values of |∇h|
106
6.6 Test 3: anisotropic meshing
Figure 6.28: Test 3. Background mesh. The mesh boundaries coincide with the domain sides. Higher refinement corresponds to areas with a higher curvature of the bathymetric isolines.
Figure 6.29: Test 3. Bathymetric gradient vectors defined over the bathymetric isolines.
107
6. THE MULTI-CRITERIA MESHING METHOD. APPLICATION TO IDEALIZED GEOMETRIES
0.5 (maximum allowed elongation) and 1 (equilateral condition). The direction of ~ ~ elongation is set perpendicular to |∇h|: consequently the two components of |∇h|,
∇hi and ∇hj , are interpolated over the background mesh in order to define ui and vi
respectively.
Given these input specifications, two meshes are generated: Mt31s with the searching front option (Figure 6.30), and Mt31c with the circulating front option (Figure 6.31). Both are showed before post-processing operations in order to evaluate the original node location determined by the anisotropic conditions. Each figure is accompanied by close-up views in order to compare the element orientation with the bathymetric contours. Mt31s has 1350 nodes and 2578 elements, Mt31c 1344 nodes and 2566 elements. In both meshes, equilateral elements are formed along the domain edges, where the bottom is flat, while anisotropic shapes appear in the central part, in correspondence of higher bathymetric gradients. The circulating front option produces a more regular grading of the element areas and orientations than the front searching option, which generates abrupt area transition and highly skewed elements. The element elongation is generally oriented according to the bathymetric isolines when the curvature of the contour is wide (external isolines of Figures 6.30b, 6.30c, 6.31b and 6.31c), while for the inner curves the direction of elongation becomes less controlled. Limitations in generating good orientated elements in correspondence of pronounced curvatures may be overcome by setting a sizing criterion depending on the local curvature of the bathymetric isoline: smaller element sizes are needed where the radius of the curvature is smaller. This option has however not been investigated further. As already mentioned, in order to increase the level of control on the final mesh and improve the accordance between element orientation and bathymetry contours, the background node positioning should be done with care: higher refinement is needed where more control is wanted on the element size, shape and orientation, especially where the bathymetric contours are irregular.
6.7
Test 4: meshing criteria for the wet and dry
A wet and dry algorithm allows the classification of mesh elements as either wet or dry, according to predetermined criteria in order to simulate a flow inundation process.
108
6.7 Test 4: meshing criteria for the wet and dry
(a)
(b)
(c)
Figure 6.30: Test 3. (a) Anisotropic mesh Mt31s, generated by the searching front option with a target element side length of 10 m. The elongation of the elements is set perpendicular to the bathymetric gradient; (b) and (c): comparison between the element orientation and the bathymetric isolines.
109
6. THE MULTI-CRITERIA MESHING METHOD. APPLICATION TO IDEALIZED GEOMETRIES
(a)
(b)
(c)
Figure 6.31: Test 3. (a) Mesh Mt31c, generated by the circulating front option with a target element side length of 10 m. The element elongation is set perpendicular to the bathymetric gradient; (b) and (c): comparison between the element orientation and the bathymetric isolines.
110
6.7 Test 4: meshing criteria for the wet and dry
The insertion of a wet and dry scheme in a computational model has a twofold facet: from a numerical point of view, the algorithm should be efficient, numerically stable, and should not introduce spurious water oscillations. On the other side the physics of the process need to be considered, as the the governing equations have to be solved in constrained conditions, as a reduced water depth. A comprehensive and updated review of existing wet and dry algorithms for circulation models based on the SWEs is given by Medeiros and Hagen (2013), where these 4 classes are identified: - thin layer algorithm: a thin layer is superimposed to the whole domain in order to include all mesh nodes in the computation. A subsequent numerical condition is applied for distinguishing wet and dry elements; - element removal : an element is classified as either dry, wet, or partially wet/dry, according to model specific conditions and checks. If an element is dry it is excluded from the computation. Also the treatment of partially wet/partially dry elements is much dependent on the model; - depth extrapolation: the water depth is extrapolated from the wetted elements to nearby dry areas, according to specific conditions (i.e. bottom friction); - negative depth: the free water surface extends also underground, and the governing equations are calculated all over the mesh. Elements are classified as dry if the water elevation is negative. The ADCIRC wet and dry scheme, to which the following tests are targeted, belongs to the element removal class: the element state is hence determined by a series of conditions based on the local flow properties and on the definition of a water depth threshold: the capacity of the mesh in reproducing the DEM isolines is fundamental, as the temporal evolution of the wetted front and the estimation of the flooded area extension are critical factors in the assessment of the flooding risk. In these cases, when designing a mesh for wet and dry studies, a modeller should consider these, sometime contrasting, guidelines (Donnel et al., 2006): - keep the smoothness and continuity of the DEM isolines; - maintain a reduced depth difference between nodes of the same element;
111
6. THE MULTI-CRITERIA MESHING METHOD. APPLICATION TO IDEALIZED GEOMETRIES
- minimize the number of elements changing state at each time step; - avoid the “saw tooth” effect (a rugged shape of the wetted front); - avoid too distorted elements. Moreover, wet and dry schemes usually burden the simulation by imposing strict Courant conditions (Blain et al., 2010). In the following numerical experiment an anisotropic factor is introduced in the mesh design in order to align the element shape along the bathymetric isolines. The mesh is tested in a tidal simulation over an idealized bathymetry. Potential improvements are assessed, in comparison with an identical simulation for an isotropic mesh, in terms of shape and evolution of the wetted front, free surface elevation, model stability and computation costs.
6.7.1
Domain description and mesh generation
The idealized bathymetry used for the definition of the NSFs is represented in Figure 6.32. It corresponds to a semi-enclosed basin with a central island, connected to the open sea by an inlet. Two scalar criteria are used for sizing the mesh elements: ρ41 is the inverse of the bathymetric gradient, according to the relation: ρ41 (x, y) =
1 ~ |∇h(x, y)|
(6.8)
The second criterion is targeted to increase the mesh refinement for depths close to 0: ρ42 (x, y) = max{(c|h(x, y)|; 10}
(6.9)
where |h(x, y)| is the absolute value of the water depth and c = 10 a proportionality
coefficient. ρ42 imposes a lower bound of 10 m to the element size for avoiding overrefinement near MSL . The final NSF is set to the minimum between the previously
defined functions. The MESHGR coefficients of anisotropy are calculated according to the procedure of the previous anisotropic test. Here si varies between 0.4 and 1. The obtained mesh, Mt4AN (Figure 6.33), has 10012 nodes and 19649 elements. Close-up views of the mesh demonstrate how the element configuration tends to align with the underlying bathymetric isolines. A second mesh, called Mt4IS and not shown here, is generated with the same final NSF of Mt4AN for isotropic conditions (si = 1), resulting in 5706 nodes and 11061 elements.
112
6.7 Test 4: meshing criteria for the wet and dry
Figure 6.32: Test 4. Bathymetric depths for the wet and dry test case.
6.7.2
Tide modelling
A periodic free surface elevation, with an amplitude of 1 m and a semi-diurnal period, is applied to the ocean boundary of both meshes. The simulation is performed with ADCIRC, where a wet and dry scheme sets whether an element is wet or dry according to empirical rules (Luettich and Westerink, 1995; Dietrich et al., 2004; Blain et al., 2010). A ∆t = 0.5 s and a linear parametrization of the bottom friction, with a drag coefficient equal to 0.0025, are adopted. The simulation with isotropic conditions did not present instability problems. For the anisotropic mesh, the set of input parameters used in the isotropic case originated instability issues in correspondence of a few elongated elements (the diagonal term of the corresponding nodes became negative). Therefore a Laplacian smoothing was applied to the entire mesh for relaxing the original anisotropic discretization and the time step was lowered to 0.25 s for achieving a stable run. The characteristics of the isotropic and anisotropic meshes are reported in Table 6.2.
113
6. THE MULTI-CRITERIA MESHING METHOD. APPLICATION TO IDEALIZED GEOMETRIES
Figure 6.33: Test 4. Wet and dry test. Mesh Mt4AN with location of the reference station.
6.7.3
Results
The water elevation time series computed for the meshes Mt4AN and Mt4IS at the reference station depicted in Figure 6.33, are shown if Figure 6.34, while the profiles of the wetted front are presented in Figures 6.35 and 6.36. Comparing the results between isotropic and anisotropic conditions, it can be stated that, for the present case: - the presence of distorted elements inside the mesh does not influence the computed water elevation; - the wetted front for Mt4AN does have a smoother profile than for Mt4IS, making it more similar to how it should look like in a real case; - the evolution of the wetted front is however coherent in both cases; when the wetted
114
6.8 Conclusions
Table 6.2: Parameters of the isotropic and anisotropic meshes for the wet and dry modelling test
Figure 6.34: Test 4. Water surface elevation at the reference station
fronts are superimposed (Figure 6.36) their relative distance is within the limits of the local mesh element size; - the anisotropic conditions impose a higher computational cost as a consequence of a higher node number and a smaller time step in order to get a stable run. It can be concluded that the alignment of the mesh elements to the bathymetric contours offers a more realistic visualization of the wetted front profile. However, if an isotropic mesh already provides a discretization where the hydrodynamic and bathymetric gradients are resolved, the introduction of an anisotropic orientation of the elements does not provide an essential improvement in the estimation of the flooded area extent, nor in the computed water elevations.
6.8
Conclusions
In relation to the exposed idealized examples of isotropic and anisotropic meshing for different domains and size functions, these conclusions can be drawn:
115
6. THE MULTI-CRITERIA MESHING METHOD. APPLICATION TO IDEALIZED GEOMETRIES
Figure 6.35: Test 4. Profile of the wetted front for Mt4AN. The close-up views (a) and (b) contrast the position of the wetted front with the nearest bathymetric isoline.
116
6.8 Conclusions
Figure 6.36: Test 4. Comparison between the position of the wetted fronts of the anisotropic and isotropic meshes, at the same time step during the simulation.
- MESHGR handles multi-connected domains properly: test 2 suggests that as many “holes” as required can be introduced, provided that the condition of a not-disjoint domain is preserved; - MESHGR demonstrates to be fully adapted to receive a background information derived from the multi-criteria method; - post processing subroutines enhance substantially the quality of the isotropic mesh. Specifically, diagonal swapping is suited to lower the maximum node valence; Laplacian smoothing, combined with the previous subroutine, decreases the number of abrupt area changes between adjacent elements; - a critical factor on which the accordance between the NSF and the actual size/shape
117
6. THE MULTI-CRITERIA MESHING METHOD. APPLICATION TO IDEALIZED GEOMETRIES
of the final mesh elements depends is the refinement of the background mesh; - anisotropic meshing is successfully performed in MESHGR, within the limits expressed in the previous point; - the multi-criteria methodology for meshing domains with highly varying bathymetries may overcome the poor level of discretization often encountered when just one single mesh criterion or a localized refinement are applied; - the introduction of anisotropic conditions in a wet a dry study gives a more realistic rendering of the wetted front profile, although the improvements in the estimation of the flooded area extension are limited.
118
7
APPLICATION OF THE MULTI-CRITERIA METHOD TO A REAL DOMAIN 7.1
Introduction
Shallow water models are often required to simulate processes extending over several spatial and hydrodynamic scales, from the wide ocean to coastal bays, estuaries and low lying areas exposed to inundation. Unstructured meshes are appropriate tools for meeting the resolution requirements of these models, where a coarse discretization in deep waters coexist with localized refinements wherever the involved phenomena present high gradients, as in the case of physical forcings, steep bottom features, detailed coastlines, etc.. This chapter describes the application of the multi-criteria meshing method to a domain going from the Western Atlantic Ocean Shelf in front of the Iberian Peninsula to the Lima River estuary. The discretization is performed by MESHGR, which demonstrates its capability to handle complex real case geometries and to generate ready-to-use meshes with element sizes varying over three orders of magnitude. The obtained mesh is then successfully tested in a wet and dry tide modelling study. In the second part of the chapter a mesh convergence study is performed, for the mentioned domain, on the basis of the a-posteriori local truncation error analysis with complex derivatives LTEA+CD (Parrish, 2007; Parrish and Hagen, 2009).
119
7. APPLICATION OF THE MULTI-CRITERIA METHOD TO A REAL DOMAIN
In the last part, the performance of the meshes derived from the multi-criteria method, from the LTEA+CD, from a uniform coarse (50 m) and a uniform fine (10 m) discretization of the Lima estuary, are compared in a wet and dry ADCIRC tide simulation.
7.2
Domain geometry and bathymetry
The domain under investigation extends from the North Atlantic Ocean to the coast of the Western Iberian Peninsula, including the estuary of the Lima River (Figure 7.1). This estuary has been object of several studies, with a focus on salinity structure and salt water intrusion (Pinho and Vieira, 2007; Vale and Dias, 2011a), or on the pollutant tracking in the lower estuary (Vale and Dias, 2011b). Rebord˜ao and TrigoTeixeira (2009) performed a wet and dry simulation for a mesh of constant element size inside the estuary: however, no quantitative results of the flooding process are reported. None of the mentioned studies includes ocean scales, nor has a specific focus on the mesh generation. The estuary represents a non-trivial task in terms of spatial discretization: the lower part has a channel-like shape with the presence of a harbour. The middle estuary presents a complex morphology with a network of channels and sedimentary banks covered by brackish vegetation, exposed to the covering and uncovering of the tide propagation. The tides are semidiurnal, with a range varying from 1.1 m on neap tides to 3.7 m on spring tides. For the determination of the domain boundary location, the river floodplains adjacent to the river course are included as well up to the 5 m bathymetric contour above MSL, in order to assure a buffer zone between the land boundary and the high water spring level. The domain is extended upstream until Ponte de Lima, at 24 km distance from the river mouth, beyond the influence of the tide propagation. Ponte de Lima is considered the furthest inland location of the tidal wave (Trigo-Teixeira, 2008). The position of the ocean boundary is similar to a previous modelling study (Ara´ ujo et al., 2011b) which gave satisfactory results in reproducing the tide levels at the tide gauge in the Viana do Castelo harbour. The DEM is obtained from different data sources: for the ocean waters, depths up to 4000 m are digitized from the 1: 150 000 and 1: 1 000 000 nautical Portuguese charts, while deeper bathymetry is provided by the 1 grid resolution global seafloor database
120
7.3 Meshing criteria for a multi-scale domain
from the Institute of Geography and Planetary Physics (Smith and Sandwell, 1997). The estuary bathymetry is made up by a series of field surveys, specifically a 5 m resolution regular grid for the lower estuary, a 25 m x 10 m resolution scatter set for the middle estuary and a series of cross-sections for the upper estuary and the Lima course. The flood plain topography is digitized from the 1: 25 000 topography map of the Portuguese Army Geographical Institute. The final scatter set is obtained in the form of a linearly interpolated triangulated irregular network (TIN).
Figure 7.1: Extension of the computational domain: (a) Ocean boundary location and bathymetry; (b) Lima estuary extension and bathymetry (negative depths correspond to points above mean sea level) and (c) middle estuary, with the urban area of Viana do Castelo and the position of the Eiffel bridge. The main and secondary river channels are visible. Numbers correspond to the deployment sites of 4 acoustic doppler current meter profilers.
7.3
Meshing criteria for a multi-scale domain
Published ocean and coastal numerical applications often omit the description of the mesh generation process and of the used sizing functions, leaving the reader uncertain about the efficiency of the proposed discretization. However, criteria for shallow water
121
7. APPLICATION OF THE MULTI-CRITERIA METHOD TO A REAL DOMAIN
models are well documented and show recurrent formulations (see chapter 4). The purpose of this paragraph is to explicitly describe the steps involved in the creation of a mesh, for the multi-scale domain previously described, from the definition of the meshing criteria to the post-processing operations.
7.3.1
Ocean, shelf break and continental shelf
One of the most recurrent NSF for the ocean and continental shelf is the wavelength (λ) to mesh size (∆x) ratio (Le Provost and Vincent, 1986): T λ ∆x1 (x, y) = = R
p gh(x, y) 750
(7.1)
where T is in this case the M2 tide period and R has been chosen equal to 750 after a series of preliminary tests, in order to obtain, for the oceanic part of the domain, the same maximum and minimum element sizes of the mesh of Figure ??, successfully tested in Ara´ ujo et al. (2011a). Mesh convergence studies (Westerink et al., 1994; Luettich and Westerink, 1995) however pointed out several shortcomings of (7.1), as the inadequate resolution of steep slope features, or the impossibility to adapt the mesh to 2D flow patterns as amphidromes (see chapter 3). Some of these limitations may be overcome by a NSF derived from the topographic length scale Lh (Huthnance, 1995, Legrand et al., 2007): ∆x2 (x, y) = c1 Lh = c1
h(x, y) ~ |∇h(x, y)|
(7.2)
~ is the bathymetric gradient and c1 a coefficient of proportionality equal in where ∇h this case to 0.25. (7.2) is superimposed upon (7.1) for increasing the mesh refinement in areas of steep bathymetry.
7.3.2
At the coast and inside the estuary
In near shore and estuarine regions the element size should not only reflect the decrease in the wavelength but also be related to the local bottom features, as the flow interaction with the bathymetry is enhanced. The criteria proposed in open regions are not valid any more and must be integrated or substituted with local ones, taking into account specific constraints. The first criterion applied to the estuary relates the mesh resolution
122
7.3 Meshing criteria for a multi-scale domain
to the module of the inverse of the bathymetry gradient: ∆x3 (x, y) =
1 ~ |∇h(x, y)|
(7.3)
A lower bound of 15 m is applied to avoid the definition of too reduced sizes when the local gradient is steep. According to (7.3) the most refined estuarine areas are expected along the river banks or the slopes of the emerged bedforms, in order to better reproduce close bathymetry isolines and represent the temporal evolution of the wetted front. The following criterion, ∆x4 , addresses the resolution of the horizontal length scales of the estuary, represented by the width of the harbour navigation channel, the river main and secondary channels, the river banks and adjacent floodplains. A minimum number of nodes per cross-section must be assured, as under-resolution may significantly decrease the system conveyance (Westerink et al., 2008). Therefore the highest refinement is attributed to the navigation and main river channels (∆x4 = 20 m), in order to guarantee at least 3 elements in the narrowest cross-section; this level of refinement is relaxed for the river secondary channels (∆x4 ranging from 30 to 35 m); the refinement for banks and marshes is set to vary from 30 to 60 m accordingly to the distance from the river longitudinal axis. A radial function, ∆x5 , is applied from the river mouth towards the ocean, up to a distance of 2 km, for assuring a smooth size grading between the NSFs defined in the estuary and in the ocean: ∆x5 (x, y) = c2
p (x − xC ) + (y − yC )
(7.4)
with c2 being a coefficient of proportionality equal to 2000, and (xC , yC ) the coordinates of the center of the river mouth.
7.3.3
Numerical constraints
Spatial discretization and numerical stability are interrelated. For ADCIRC, applications with the wet and dry option should respect a severe Courant condition (Blain et al., 2010). The related NSF can be expressed as: ∆x6 (x, y) =
123
VM AX ∆t 0.1
(7.5)
7. APPLICATION OF THE MULTI-CRITERIA METHOD TO A REAL DOMAIN
where VM AX is the module of the maximum velocity, calculated for the same domain extension and input conditions of the present case, at the nodes of a “first try” mesh, and ∆t the computational time step. In practice ∆x6 acts as a lower bound over the whole domain, assuring a Cr < 0.1.
7.3.4
Final node spacing function
The final NSF, to be read in MESHGR through a background unstructured mesh, is calculated as: ∆xf (x, y) = max{min (∆xi )5i=1 ; ∆x6 }
(7.6)
where ∆x1 to ∆x5 operate as minimization conditions on the final target size function, while ∆x6 sets a spatially variable lower bound on the mesh element size. MESHGR discretizes the domain by the advancing front searching method. Laplacian smoothing is applied for assuring better area size transitions between adjacent elements. The obtained mesh, called M1 (Figure 7.2), has 63 516 nodes and 124 035 elements. The element size varies throughout the domain from 17.6 km in deep waters to 8 m inside the estuary. Areas of refinement in the open ocean are observed in correspondence of seamounts and along the shelf break and continental slope. Inside the estuary the highest level of refinement is attained along the main river channel, with spot-like node concentrations due to elevated local bathymetry gradients. Figure 7.3 shows the distribution of the element shape quality index qi (formula 3.13).
Figure 7.2: Unstructured graded mesh M1 derived from the node spacing function ∆xf . Localized refinement is seen (a) in correspondence of seamounts and of the shelf break, and, for a closer estuary view (b), along the main river channel and where the bathymetry gradient is pronounced. Numbers refer to the ADCPs locations.
124
7.4 Model calibration
Figure 7.3: Distribution of the element quality index qi for the post-processed mesh
7.4
Model calibration
A wet and dry tide modeling study is performed by means of ADCIRC, the 2D version of the finite element ADvanced CIRCulation coastal ocean model (Luettich et al., 1992), which solves the shallow water depth integrated barotropic mass and momentum equations under the hypothesis of incompressibility, Boussinesq and hydrostatic pressure. The governing equations, where the free surface stresses and atmospheric pressure gradients are neglected as not considered in this study, are (Westerink et al., 2008): ∂Ui h ∂(Uj hcosθ) + =0 ∂λ ∂θ Ui ∂Ui Uj ∂Ui tanθ ∂Ui + + − Ui + f Uj = ∂t rcosθ ∂λ r ∂θ r g ∂(η − αχ) νT ∂ ∂Ui h ∂Ui h − + + − τb U i rcosθ ∂λ h ∂λ ∂λ ∂θ ∂η 1 + ∂t rcosθ
∂Uj Uj ∂Uj Ui ∂Uj tanθ + + + Ui + f Ui = ∂t rcosθ ∂λ r ∂θ r g ∂(η − αχ) νT ∂ ∂Uj h ∂Uj h − − τb U j + + r ∂θ h ∂θ ∂λ ∂θ
(7.7)
(7.8)
(7.9)
where η is the free surface elevation relative to the MSL, h = η + d the total water depth, with d being the bathymetric depth relative to the MSL, λ and θ the latitude and longitude coordinates, Ui and Uj the depth averaged horizontal velocities, r the
125
7. APPLICATION OF THE MULTI-CRITERIA METHOD TO A REAL DOMAIN
Earth radius, g the gravity acceleration, f the Coriolis parameter, χ the Newtonian equilibrium tidal potential, α the effective earth elasticity factor, τb the bottom stress, dependent on the averaged velocity by a quadratic formulation, and νT the horizontal eddy viscosity. ADCIRC provides a wetting and drying option tested for a number of idealized and real geometries (Blain et al., 2010). The bottom friction is expressed in its two components by a quadratic horizontal velocity dependent law: τbi = ρw Cf τbj = ρw Cf
p
p
U 2 + V 2 Ui
(7.10)
U 2 + V 2 Uj
(7.11)
with ρw being the fluid density and Cf the drag coefficient. For the present case Cf has been expressed in terms of the Manning roughness coefficient n (Luettich et al., 1992): Cf =
gn2 (d + η)1/3
(7.12)
n is assumed to be space variable according to the land cover type detected from available aerial photos of the estuary. For the present simulations, the Manning coefficients of Table 7.1 associated to each land type are taken from Chow (1981) and Bunya et al. (2010). As input parameters, a νT equal to 5 m2 s−1 and a ∆t of 0.5 s are used, which Table 7.1: Manning values associated to the estuary surface cover types, according Chow (1981) and Bunya et al. (2010)
were able to stabilize the present model run. A simulation of the astronomical tide is carried out in order to assess the model capacity to reproduce, for the given domain discretization, the hydrodynamics of the estuary. The modelled tides are compared with measurements of a field campaign which took
126
7.4 Model calibration
place between the 5th and the 14th October 2006, providing time series of water elevation and velocities of 4 bottom anchored Nortek 2 MHz acoustic Doppler current profilers (ADCPs), deployed along the main river channel (Figure 7.1): sites 1 and 2 were located respectively downstream and upstream of the Eiffel Bridge; site 3 in the middle estuary and site 4 along the river course, in an area still strongly affected by the tide propagation. The free surface levels were corrected from the influence of atmospheric pressure variations. The ADCP at site 2 was re-deployed after a local fisherman accidentally collected it: the measures of this instrument, being potentially subject to errors, are not used for the model validation. The model is forced by freshwater inflow at the upstream river section, with daily averaged flow rate values registered at Ponte da Barca, which is the nearest hydrographic station with available data for the analysed period, located at 17 km from the domain upstream boundary. The river record is made available by the Portuguese Hydraulic Resources National Information System. Consequently to the presence of the upstream dams of Alto Lindoso and Touvedo, the river regime in normal conditions corresponds to the minimum environmental flow rate. The flow rate for the present simulation time varies between 2.73 and 15.15 m3 /s. The ocean boundary nodes have been forced with amplitude and phase values of 13 astronomy harmonic constituents (K1, K2, L2, M2, MU2, N2, NU2, 2N2, O1, P1, Q1, S2 and T2), extracted from the Le Provost global tide dataset (Le Provost et al., 1998). In addition, the tide potential functions of the mentioned constituents have been applied to the mesh interior nodes. Free lateral slip is allowed at the remaining land boundaries and at the wet and dry interface.
7.4.1
Free surface
The model shows a good reproduction of the tide ranges and phases, for the sites 1, 3 and 4 (Figure 7.4). Results for site 2 are not reported as the instrument, which suffered severe manipulation, may show a bias in the measures. The model reproduces the asymmetry of the tide curve registered at site 4, with a steeper rising curve and a delayed falling water profile. Also the damping of the tide range along the estuary is well reproduced, with an average decrease between sites 1 and 4 of 0.8 m at spring tide, while the measured difference is 0.65 m. The maximum calculated phase difference between the same sites occurs at low water and is around 2 hours, while the averaged registered lag is 1.5 hours. An underestimation of the water elevation at high tide, not
127
7. APPLICATION OF THE MULTI-CRITERIA METHOD TO A REAL DOMAIN
present for site 1, appears evident at site 4. This difference remains unaltered even if a twofold refined mesh is used, suggesting that the committed error has not to be attributed to the adopted mesh. A reason of the model inaccuracy may rely in the mass balance errors of the wet and dry scheme, which become relevant when the intertidal areas occupy an extensive part of the estuary (Blain et al., 2010).
7.4.2
Depth averaged velocities
The magnitude of the predicted and observed depth averaged velocities during the whole tide cycle is compared in Figure 7.5. At site 1 the measured series has a strong flood flow registered on the 9th October, which might be of meteorological origin. The predicted velocities reproduce generally well the observed time series, with an optimal phase accordance. From processing the data of ADCP at site 1, the observed velocity component along the Latitude direction shows a cut-off during the flood phase of the spring tide, which is not noticed in the Longitudinal component, and disappears at neap tide. This behaviour may be related to some errors in the data acquisition and explains the overestimation of the flood current, while the ebb tide intensity is in agreement with the observations. At site 3 the peak flood velocity is underestimated of around 0.2 m/s, while the difference in the ebb phase is again minimal. At site 4 the slight asymmetry between flood and ebb tide is well reflected in the modelled currents: the computed current profile is able to reproduce the deformation of the tidal curve during ebb tides, with however a damping of the peak current of around 0.2 − 0.3 m/s.
7.4.3
Model validation: flooded area extension
The capacity of the model in reproducing the wet and dry in the estuary is assessed with the support of time referenced aerial orthophotos. The calibrated model, forced by the same 13 astronomy constituents of the previous simulation and at the upstream section by the river flow rate, is run between 10th and 14th September 2010, in order to include the moment when the photos were taken. Figure 7.6 represents the position of the modelled wetted boundary with respect to the real limits of the inundated area, at 10:30am of September 13th 2010. The extension of the flooded area in the orthophotos has been digitized manually. Wherever the visual detection of the water extent was uncertain in the photos, due to the reduced depth of the water column, the bathymetry contour 0.84 m below MSL, corresponding to the tide level at the harbour
128
7.4 Model calibration
(a)
(b)
(c)
Figure 7.4: Comparison of the observed and computed water elevations at (a) site 1, (b) site 3 and (c) site 4.
129
7. APPLICATION OF THE MULTI-CRITERIA METHOD TO A REAL DOMAIN
(a)
(b)
(c)
Figure 7.5: Comparison between the magnitude of the observed and computed depth averaged velocities at (a) site 1, (b) site 3 and (c) site 4. Positive values correspond to the flood phase.
130
7.4 Model calibration
gauge, is followed. The parts of the estuary which are flooded in the photos, the main and secondary channels, result flooded also in the wet and dry simulation. The distance between the wetted boundaries is in general lower than 50 m, a value which is comparable with the local mesh element size. It can be concluded that the used mesh represents accurately the underlying DEM and that the flooding process during a tidal cycle is realistically simulated. The following series of figures visualize the extension of
Figure 7.6: Comparison between the locations of the wetted fronts taken from aerial photos and from the ADCIRC wet and dry simulation, for the Lima middle estuary. Discrepancies are comparable with the local mesh element size. Small polygons belonging to the modeled wetted boundary represent cells which have remained flooded from previous tide cycles.
the flooded area along a spring tide cycle. Figures 7.7 (a) to (d) allow to appreciate the extent of the intertidal area along different phases of a spring cycle: at spring tide high water the river sediment banks and secondary channels are almost completely flooded, while at spring tide low water only the main river channel remains wet. The elevated extension of the intertidal area in the middle estuary, which for the simulated tide cycle is 63% of the inundated area at high tide, may explain the presence of mass balance errors and be related to the vertical off-set between water level observations and predictions.
131
7. APPLICATION OF THE MULTI-CRITERIA METHOD TO A REAL DOMAIN
(a)
(b)
(c)
(d)
Figure 7.7: Water surface extension in the Lima estuary calculated along the day 8th October 2006: a) low tide at 9am; b) increasing tide at 12pm; c) high tide at 15pm; d) low tide at 22pm.
132
7.5 Evaluation of an a-posteriori meshing criterion
7.5
Evaluation of an a-posteriori meshing criterion
The multi-criteria meshing method described in the previous paragraph proved to generate a good quality mesh that was successfully applied to a tidal simulation. As a conceptually different approach, the a-posteriori LTEA+CD (Parrish, 2007; Parrish and Hagen, 2009; Bacopoulos et al., 2011) is applied here to the multi-scale domain previously defined. The LTEA+CD is a relatively recent discretization methodology formulated and applied so far only to full tidal zones and not yet to DEM exposed to flooding. A mesh convergence study is developed here in order to assess if the algorithm is capable of producing meaningful discretizations and ready to use meshes for wet and dry applications as well. In this sense, a series of model performance indices will be discussed.
7.5.1
The LTEA+CD meshing method
The multi-criteria methodology presented in the previous section relies on a mix of a-priori and a-posteriori criteria combined together to define the final NSF for the meshing. The LTEA+CD, specifically targeted to oceanic, coastal and estuarine domains, is based on the principle of redistributing the nodes of an existing mesh in order to get the local truncation error as uniform as possible all over the computational domain. The truncation error expression is obtained as the difference between the continuous and the discrete non-conservative shallow water momentum equation, in its linearized harmonic form, using complex derivatives and without the inclusion of the advective terms. The LTEA+CD represents an evolution of the previous LTEA (Hagen et al., 2000; 2001; 2002; 2006; Hagen and Parrish, 2004) and LTEA-CD (Parrish and Hagen, 2007) methods, as, in the error expression, the estimation of the velocity from a non-linear simulation, the spatially variable bottom stress and the Coriolis force, not included in the LTEA and LTEA-CD, are now present. The truncation error estimator + τˆM E is expressed as (Parrish and Hagen, 2009):
∆x6 ˆ [ω(iv0 − u ˆ0 ) + (τb vˆ)0 + ˆi(τb u ˆ0 ) − ˆi(f vˆ)0 + (f u ˆ)0 ](6) (7.13) 1440 √ where ω is the frequency of the single harmonic constituent, ˆi = −1, ∆x the local + τˆM E =
element size (for hypothesis the triangles are supposed to be equilateral), u and v the
133
7. APPLICATION OF THE MULTI-CRITERIA METHOD TO A REAL DOMAIN
velocities in the x and y direction, u ˆ and vˆ the harmonic velocities in the x and y direction. The suffix 0 refers to the central node of the “valence cell”, that is a group of 6 equilateral triangles superimposed to the mesh and having in common the mesh node 0. τb is the spatially variable bottom stress, calculated here on a per node basis by a quadratic formulation (formulas 7.10 and 7.11) using the velocities derived from the harmonic solution and averaging in the time, and f is the Coriolis term. The superscript (6) refers to the order of complex differentiation. From (7.13) it is therefore + possible to set a target τˆM E (i.e. a uniform maximum allowed value) and solve for ∆x,
which becomes the target element size function of the adapted mesh. The estimation of the bottom stresses make the method particularly attractive for the meshing of very limited water depths, where interactions between the flow and the bathymetry are enhanced. The method has demonstrated to improve the model efficiency by creating meshes with less nodes/elements but with the same or even higher performance of the existing mesh. Here the LTEA+CD is applied, for the North Atlantic Ocean to Lima estuary domain extension of Figure 7.1, to an existing a-priori mesh (paragraph 7.5.2). As the LTEA+CD algorithm has been implemented only for deep and shallow waters constantly covered by the tide, the aim of this application is to explore the algorithm behaviour when the re-meshing includes areas above low tide: this is achieved by running the model in a wet mode: intertidal and supratidal depths are fictitiously decreased to the minimum low tide level, in order to assure that all elements remain wet during the simulation. The procedure, fully automatized within the SMS mesh generation toolbox, consists of the following steps: 1) a simplified linear ADCIRC simulation is run with the existing mesh, by forcing the domain with the M2 tide constituent. Although more constituents may be selected, only the M2 astronomical component has been chosen here as dominant at the estuary entrance (Mazzolari et al., 2012). The model is run in a wet mode, meaning that the intertidal and supratidal zones were attributed a minimum depth in order to remain always wet; 2) a harmonic analysis of the elevation and velocity output is performed on a per node basis;
134
7.5 Evaluation of an a-posteriori meshing criterion
3) the LTEA+CD uses the results of step 2) to calculate the localized error and, as a function of it, a target NSF, in a way that the prospective mesh is designed so as to redistribute the truncation error uniformly; 4) the target NSF is used by the SMS mesh generator to create the new mesh. A scale factor is applied to the NSF for generating a mesh with the specific number of nodes selected by the user. The nodes redistribution is allowed along the domain boundaries too. A smoothing function is applied at the end to limit the maximum area change transition between adjacent elements. A mesh convergence study is presented hereafter: a series of adapted meshes are generated, starting from an initial, a-priori mesh with a coarse estuary discretization, and tested in a tide modelling study.
7.5.2
A-priori coarse estuary discretization
An initial mesh is designed according to the a-priori criteria of the wavelength to mesh side ratio (formula 7.1) and of the topological length scale (formula 7.2) for the ocean domain, and with a constant 70 m element side along the estuary. A recurrent empirical discretization approach to domains dominated by channel-like features imposes a minimum number of elements per channel cross section (Parrish, 2007; Westerink et al., 2008; Salisbury et al., 2011), as a rough representation of the cross section bathymetry may reduce the wetted area and consequently the discharge rates (Feyen et al., 2002). The minimum number of required elements varies according to the modeller’s needs and application options (with or without wet and dry; with or without slip boundary condition). For the Lima estuary, considering that the navigation channel and main river channel minimum width is around 230 m, a 70 m constant element size is chosen to guarantee at least 3 elements per cross section. This grade of resolution is extended in the intertidal areas and overland until the land boundaries. It has to be noted that some sections of the upper river course, of the estuary tributary and secondary channels remain under-refined. The mesh (Figure 7.8), generated by MESHGR, is called M70, and has 31784 nodes (-50% with respect to M1) and 61502 elements.
135
7. APPLICATION OF THE MULTI-CRITERIA METHOD TO A REAL DOMAIN
(a)
(b)
Figure 7.8: Mesh M70 for the (a) ocean and (b) estuarine domain, with dark lines representing the isoline 0 m with respect to MSL. It is possible to appreciate that at least 3 elements are guaranteed across the narrowest main river channel section.
136
7.5 Evaluation of an a-posteriori meshing criterion
Table 7.2: Mesh properties for the LTEA+CD convergence study.
7.5.3
LTEA+CD mesh convergence study
In this mesh convergence study the LTEA+CD is iteratively applied to M70. The newly created meshes are tested, under the same model conditions, for establishing the simulation performance and computational gains of the new model configuration against the M70 benchmark. At each refinement, the velocity amplitude and phase, the water surface elevation amplitude and phase, the bottom shear stress and the Coriolis force, calculated on a per node basis from the linear run, are used as input terms for the error estimation, according to the 4 steps of paragraph 7.5.1. Two meshes are created by setting a target node number equal to the nodes of M70: M70 1, derived from M70, results in 31763 nodes and 60596 elements; M70 2, derived by M70 1, has 31594 nodes and 59931 elements. A third mesh, M70 D, is generated from M70 by imposing the double of the original node number, and results in 60075 nodes and 116023 elements. All discretizations are performed by the SMS paving method. Their main properties are reported in Table 7.2. Each refinement, for either a constant or an increased node number, reduces the overall mesh quality, even if a smoothing of the LTEA+CD output is applied. This property may be explained by the enhanced variations that the LTEA+CD method introduces to the original element size distribution. The time step used for the wet and dry simulations is controlled by the smallest element size. Figure 7.9 depicts the localized truncation error for meshes M70, M70 1 and M70 2 along a domain transect going, from left to right, from the Western Ocean boundary to the upstream section of the Lima River: a clear decreasing trend of the local error is visible, along the ocean shelf, for successive re-meshing attempts. This behaviour is in line with the LTEA+CD documented literature results. The spike on the right side of Figure 7.9a corresponds to the steep slope of the shelf break: in this area the
137
7. APPLICATION OF THE MULTI-CRITERIA METHOD TO A REAL DOMAIN
(a)
(b)
Figure 7.9: Local truncation error calculated for meshes M70, M70 1 and M70 2 along a domain transect going, from left to right, (a) from the Western ocean boundary to the Lima river mouth, where the isolated peak corresponds to the shelf break, and (b) from the toe of the continental slope to the upstream boundary of the estuary domain.
138
7.5 Evaluation of an a-posteriori meshing criterion
elements undergo rapid size changes due to the transition from the deep ocean shelf to the continental shelf. It has been demonstrated that excessive area changes, often associated to not completely equilateral elements, cause the truncation error to increase (Hagen et al., 2001). The sudden rise of the error in all meshes when approaching the estuary mouth corresponds to the increase in the velocity gradients, as shallower depths are entered and the flow has to adapt to the harbour and estuary constrained geometry. A close-up view of the error profile in this area and inside the estuary (Figure 7.9b) demonstrates that the error measure varies with sharp changes in very limited depths, as a consequence of the high gradients of the velocity fields encountered here. The error estimation inside the estuary, on the contrary of what observed for the deep water region, does not only decrease or become uniform, but seems to increase for successive LTEA+CD applications. This behaviour is probably directly related to the physics of the problem: as the mesh is more refined, the bathymetry interpolated to the model becomes more detailed and more irregular (smaller-scale features are represented), while the model reflects the enhanced irregularity with more pronounced gradients in the flow fields, which on their turn increase the error. The error estimation for M70 D, not shown in the figures, is in general higher than for all the other meshes. Again, this may be explained by considering that doubling the number of nodes corresponds to a more detailed bathymetry representation. A series of figures are added to describe some refinement properties of the NSFs introduced by the LTEA+CD method. Figures 7.10 and 7.11 show that, for the oceanic domain, areas of mesh refinement are mainly located along the seamount slopes, the continental shelf break and shelf slope, and along the coastline for very shallow waters. These locations are all characterized by an increased localized error than nearby areas, as illustrated in Figure 7.9. Most of the estuary is subject to refinement, with bigger element size variations occurring when M70 1 is obtained from M70 than when M70 2 is obtained from M70 1. However, in both cases, a higher grade of refinement is located along the navigation channel and the main river channel; white areas, indicating where the mesh become coarser, correspond to the intertidal flat, sandbanks or minor estuary channels. The harmonic analysis of the linear run output shows that the water surface elevation and phase have a smooth gradient throughout the estuary, with a damping of the tidal signal in the upstream direction. On the contrary, the velocity field has more
139
7. APPLICATION OF THE MULTI-CRITERIA METHOD TO A REAL DOMAIN
(a)
(b)
Figure 7.10: (a) Areas of refinement in the ocean for M70 1 with respect to M70. Areas in white correspond to element relaxation; (b) differences between the element side lengths of mesh 70 1 with respect to mesh M70 for the central part of the Lima estuary. Areas in white correspond to element relaxation.
140
7.5 Evaluation of an a-posteriori meshing criterion
(a)
(b)
Figure 7.11: (a)Areas of refinement for M70 2 with respect to M70 1. Areas in white correspond to element relaxation; (b) differences between the element side lengths of mesh 70 2 with respect to mesh M70 1 for the central part of the Lima estuary. Areas in white correspond to element relaxation.
141
7. APPLICATION OF THE MULTI-CRITERIA METHOD TO A REAL DOMAIN
pronounced changes, and it is easy to observe that an higher amount of mesh refinement occurs where the flow velocity gradient is higher. This occurs especially along the main river course and around the singularity points of the estuary land boundary, as these features introduce local disturbances in the flow field which affect the LTEA+CD output. A more uniform flow field is reproduced over the intertidal and supratidal zones, most probably because of the artificial, constant minimum depth imposed to all points potentially becoming dry during a tide cycle, with the result of flattening the bathymetry.
7.5.4
Model application
A series of astronomical tide simulations has been performed, with the wet and dry option turned on, by maintaining the set of parameter values and input boundary conditions which gave a stable run for M1 (see paragraph 7.4). In each simulation only the mesh changes, with depths interpolated linearly from the underlying DEM. The validation period goes from 15:30 of 5th October 2006 to 00:00 of 14th October 2006, for which the ADCP water elevations and currents at sites 1, 3 and 4 are available for comparison. An appropriate set of statistical indices is chosen for the quantification of the model errors with respect to the measurements. The chosen indices are the mean average error (MAE), the root mean square error (RMSE) and the Nash-Sutcliffe efficiency coefficient (NSE, Nash and Sutcliffe, 1970), defined as: M AE =
N 1 X obs |yi − yicomp | N
(7.14)
i=1
RM SE =
N 1 X obs (yi − yicomp )2 N i=1
N SE =
PN
obs i=1 (yi
!0.5
P obs − y comp )2 − y¯)2 − N i=1 (yi i PN obs − y 2 (y ¯ ) i=1 i
(7.15)
(7.16)
where N is the sample dimension, equal here to 1207, yiobs the ith measurement and yicomp the ith computed value of y. y¯ is the average of the measured values. MAE and RMSE provide a statistical error estimation in the variable units, while NSE is a normalized performance index, with an efficiency of 1 for a perfect reproduction of the event. Table 7.3 reports the water elevation and depth average velocity performance
142
7.5 Evaluation of an a-posteriori meshing criterion
Table 7.3: Performance indices for meshes M70, M70 1, M70 2 and M70 D, at sites 1, 3 and 4.
indices, at sites 1, 3 and 4, for the meshes M70, M70 1, M70 2 and M70 D. As a general overview of the model efficiency, the mesh performance is higher in the reproduction of the free water surface than the currents, while the lower range of variability of the indices for the water surface than for the velocities demonstrates that the model is more sensitive to the current fields than to the free surface. These results are commonly observed for depth averaged hydrodynamic models (French, 2010). The model efficiency also diminishes along the estuary upstream direction: the water elevation maximum NSE at site 1 is 0.989 against 0.825 at site 4; for the currents the NSE decreases from 0.801 at site 1 to 0.763 at site 4. At site 3 the velocity performance is better than at site 1: this result can be explained by an anomaly in the measurements at site 1, where the Latitude component shows a cut-off during the flood spring tide phase, which disappears at neap tide and which is not observed at the other sites. It has to be noticed that, going upstream, also the MAE and RMSE are increasing, even if the ranges of variation of the water surface and of the currents are lowering. This evidences that the model accumulates errors, as the transport process associated to a non-conservative
143
7. APPLICATION OF THE MULTI-CRITERIA METHOD TO A REAL DOMAIN
mass algorithm for the wet and dry propagates inside the estuary. More refined bathymetric sources, for both the ocean and the estuarine areas, may clarify if the model results can be improved. If the performance of the different discretizations is evaluated, at site 1 the water elevation indices remain fairly constant, while an improvement in the NSE is observed for the velocities for the LTEA+CD derived meshes. At site 3 the results are somehow contrasting, where no mesh clearly performs better for all indices: quite surprisingly all the re-meshing attempts perform the same or slightly worse than M70 in the water elevation; in the depth average velocity however M70 2 performs better. At site 4 the trend is clearer, with M70 D outperforming the other meshes in all fields, apart the water elevation MAE. It has to be noted here the improvement made by the LTEA+CD discretizations in the current prediction at site 4, where the absolute errors were substantially decreased and the NSE increased with respect to M70. This is the only site where the a-posteriori meshing approach introduces a substantial improvement in the model performance, suggesting that the effects of the adapted discretization are propagated throughout the whole domain. Considering that the performance values usually differ only for the second or third decimal case, it is likely that all the proposed discretizations, apart from M70, achieve an acceptable convergence for practical modelling applications. In order to visually represent the magnitude of the differences between runs with significant performance discrepancy, the time series of the computed velocities for M70 and M70 D are compared with the measured time series at sites 3 and 4 (Figure 7.12). The correspondent plots of the water surface elevations show very limited changes and therefore the reference time series remain those of the validated model (see figure 7.4). Both simulations show a damped tidal signal in comparison with the measurements; however, significant differences are noticeable: at site 4 M70 overestimates the ebb peak velocity, while M70 D underestimates it. At site 4 the peak velocity values for M70 D are closer to the measured ones than for M70 for both ebb and flood phases, along the whole tidal cycle.
144
7.6 Assessment of alternative meshing criteria for a multi-scale ocean to estuary domain
(a)
(b)
Figure 7.12: Comparison between measured and computed depth averaged velocities for meshes M70 and M70 D at (a) site 3 and (b) site 4.
7.6
Assessment of alternative meshing criteria for a multiscale ocean to estuary domain
Here a mesh convergence study is made on the basis of a series of meshes, corresponding to different meshing criteria and approaches, applied to the previous tide simulation problem. The purpose of this test is dual: from one side the level of accuracy and efficiency between different discretizations are assessed, in order to understand which meshing criteria guarantee an acceptable performance with the fewest possible nodes and elements. From the other side this analysis explores the influence of the chosen criteria in rendering the temporal wet and dry front evolution.
145
7. APPLICATION OF THE MULTI-CRITERIA METHOD TO A REAL DOMAIN
Table 7.4: Properties of the meshes corresponding to different meshing criteria
7.6.1
Mesh refinement series
Four different meshes are tested: M1 (Figure 7.2); M1 LT, obtained by applying the LTEA+CD procedure to M1 with an aproximately constant node number and a minimum element size of 10 m; M70 (Figure 7.8); and M10, which imposes a constant element size of 10 m for the Lima estuary up to the domain land boundaries, while the resolution for the open sea is unaltered with respect to M1. The fine 10 m resolution is adopted in order to obtain a sufficient representation of all secondary estuary channels (at least 3 nodes per cross-section). M10 is generated by MESHGR and has 289523 nodes (4.2 times more than M1) and 572011 elements. The mesh resolution is just the double of the 5 m resolution of the surveyed bathymetry. Table 7.4 summarizes the main properties of these meshes. It is to be noticed the general high quality of the element shapes. The time step used for the wet and dry tide modelling simulations is related to the smallest mesh element size.
7.6.2
Results
An ADCIRC astronomical tide simulation is performed for each mesh, adopting the model set-up and input conditions described in paragraph 7.4. For the assessment of the model validation on a quantitative basis, the performance indices MAE, RMSE and NSE are evaluated at sites 1, 3 and 4 (Table 7.5). The considerations expressed for Table 7.3 are valid also for this series of refinements: the errors and the variability of the performance indices tend to increase if following the estuary upstream direction, while the velocity field is more sensitive to the spatial discretization than to the water elevation. The indices do not identify univocally an optimal discretization: if analysing the water surface results, at site 1 the performance is almost identical for all meshes, while at sites 3 and 4 some limited discrepancies appear, with M10 having the best performance as a consequence of the very refined discretization along the estuary. M1 LT
146
7.6 Assessment of alternative meshing criteria for a multi-scale ocean to estuary domain
Table 7.5: Performance indices for meshes M1, M1 LT, M70 and M10, at sites 1, 3 and 4.
147
7. APPLICATION OF THE MULTI-CRITERIA METHOD TO A REAL DOMAIN
performs as M1 at sites 1 and 3, while it has the worst output of all meshes at site 4. Concerning the depth averaged velocity indices, M1 and M1 LT have the best values at sites 1 and 4, while M1 LT outperforms the other meshes at site 3. It has to be noted, for the two most upstream sites, the poor performance of M70. M10, although reproducing the estuary bathymetry with the highest resolution, does not show a clear improvement in reproducing the local current intensities. This result is in contrast with the analysis of the water surface elevation index and with the logical perception of the modeller, who introduces more refinement where more accuracy is desired. It is believed that 3D phenomena as stratification along the water column, registered during the measurement campaign (Compass Hydrographic Services, 2006) and reported also in other studies (Vale and Dias, 2011a) are not accounted in a 2D depth averaged model, and may introduce a flaw in the performance analysis. On the contrary, the values obtained for the water surface elevation, that is less dependent than the velocity on local flow patterns, attribute the highest performance to M10. M1 LT does not show any clear performance improvements with respect to the mesh it originates from. Apart from the local 3D disturbances, the apparently scarce effectiveness of the LTEA+CD algorithm may be explained by the fact that the re-meshing is targeted to an optimal node redistribution based on the M2 derived error only, while the simulation uses 13 input astronomical constituents, whose interaction increases in the very shallow waters of the estuary. Moreover, a sub-optimal node placement is most probably occurring inside the estuary for the presence of boundary singularity points, which introduce fictitious distortions in the velocity field. The capacity of the proposed discretizations in reproducing the position of the wetted front and the flooded area extension inside the estuary is assessed as well: similarly to what did for M1, the modelled wetted front and the limit of the estuary submerged areas digitized from aerial photographs are compared for 13th September 2010, at 10.30am, for all meshes. The areas bounded by the real and the computed fronts are summed up for the portion of the estuary occupied by two sandbanks exposed to wet and dry (Figure 7.13). The lower the shaded area, the more faithful the flooding simulation is. The results for each mesh are reported in Table 7.6. The series of meshes of the first convergence study are included as well. Table 7.6 indicates that the LTEA+CD method increases the accuracy in reproducing the spatial extension of the wetted area, even if the mesh node number is maintained constant for the re-meshing: this occur
148
7.7 Conclusions
Table 7.6: Assessment of the wet and dry area extension for alternative meshing criteria.
for M1 LT with respect to M1; for M70 2 with respect to M70 1, and for M70 1 with respect to M70. As previously observed, the re-meshing tends to localize the refinement where the flow has higher velocity gradients, as occurring along the channel sides and sandbanks. If a severe flooding event is reproduced, with the water surface level exceeding the normal tidal range, the wet and dry front would extend much more overland, and a new evaluation of the LTEA+CD capacity to refine more extended areas should be made. M70 has the lowest accuracy, being the coarsest mesh inside the estuary and being designed independently of the bedforms, while M10 has the highest accuracy. The average distance between boundaries for M10 is 9.4 m, which corresponds fairly well to the element size and which represents a validation of the accuracy of the followed methodology.
7.7
Conclusions
In this chapter several meshing methods have been tested for the discretization of a multi-scale real domain, going from the Western North Atlantic Ocean to the Lima River estuary, up to the limits of tidal propagation. All generated meshes allow to simulate a wet and dry process within the estuary, with a varying model accuracy. In particular: - the MESHGR algorithm has proved to be robust and flexible enough to handle complex boundary shapes and to accept multiple requirements coming from the multi-
149
7. APPLICATION OF THE MULTI-CRITERIA METHOD TO A REAL DOMAIN
Figure 7.13: Evaluation of the area gap between the wetted boundary extracted from aerial photos and the computed wetted boundary of mesh M1 LT, for September 10th 2010, at 10.30am.
criteria method applied to a real, multi-scale domain; - the ADCIRC model configuration found in the model validation has resulted in being stable also for all the successive discretizations; - the LTEA+CD method has been applied with success to a domain exposed to wet and dry; the re-meshing has demonstrated to be more sensitive to the gradient of the velocity field and to the boundary singular points, which should be avoided for maximizing the efficiency of the method, than to the water elevation; - the a-priori criteria of placing 2 to 3 elements in the lower and middle estuary channel cross-sections may lead to inaccuracies in reproducing the velocity field of the upper part of the estuary;
150
7.7 Conclusions
- the accuracy in reproducing the wet and dry front position varies between half of the local element size for the coarse uniform mesh to the full element size for the finest mesh; the LTEA+CD derived meshes outperform the original mesh in reproducing the position of the wetted fronts. More tests should be carried out in the case of storm surges.
151
8
PSEUDO-ISLANDS 8.1
Introduction
In river and coastal flood modelling applications, given an appropriate NSF describing the problem geometry, the capacity of the mesh generator to resolve sudden elevation changes and hydraulically relevant objects is fundamental for the resolution of the associated partial differential equations describing the flow and for the rendering of the flooded area extension: The unstructured mesh interpolating the DEM should ideally retain all the natural and man-made small scale features, as channels, roads, railways, levees and vegetation obstructions, which are likely to influence the water levels and flood propagation pattern, as this may improve the accuracy of the model results (Yu and Lane, 2006; Shen et al., 2006; Schubert et al., 2008). The standard advancing front method, as exemplified by MESHGR, follows this implementation order: the domain boundaries are discretized; the advancing front is initialized; a departure zone is set according to a predefined geometrical condition, and internal nodes and elements are created as long as the front is not empty. However, tests with idealized geometries presented hereafter demonstrate that, in presence of a NSF with locally high gradients, reproducing cases such as shock fronts, hydraulic jumps, bores or steep bathymetries, the standard procedure may lead to an inaccurate discretization or to a poor mesh quality. In order to overcome this drawback, a new advancing front initialization procedure is proposed, where equilateral triangles, called pseudo-islands, are introduced inside the domain in correspondence of a NSF derived condition: this may be represented by threshold value associated to the NSF gradient
153
8. PSEUDO-ISLANDS
or to NSF local minima. The triangles have the barycentre in the point where the condition is met, the side lengths equal to the local NSF value and the side segments orientated according to the left hand side rule. The sides of the pseudo-islands are saved together with the segments derived from the usual discretization of the domain boundaries in the advancing front. If the departure zone is set to the shortest segment of the active front (front searching option), the formation of elements will radiate from the sides of the pseudo-islands, in an inside to outside direction of propagation of the front, until the domain external boundaries are reached. The element size is supposed to increase for increasing distances from the nearest pseudo-island. The priority given to the formation of the smallest elements improves the ability of the algorithm to resolve the remaining unmeshed zone (George and Seveno, 1994). The following section describes how the concept of pseudo-island has been implemented in the MESHGR standard AFM algorithm; two test cases for idealized geometries and meshing criteria are successfully resolved with the introduction of pseudo-islands; in the following paragraph different locations of a LIDAR derived DEM of the Portuguese coast are selected to perform a series of discretization tests and a methodology is defined to assess the accuracy of a mesh with respect to the underlying bathymetry scatter set. The results of the mesh accuracy tests with and without pseudo-islands are presented together with an application of a pseudo-island derived mesh to a wet and dry simulation. Finally, some conclusions are drawn in the last section.
8.2 8.2.1
Test cases for idealized geometries Test case 1: NSF with one spike
An idealized NSF, defined over a squared domain, is characterized by a constant ρ equal to 22 000 m and with a spike corresponding to a target value of 2500 m, as visualized in Figure 8.1. The standard MESHGR algorithm produces the meshes represented in Figure 8.2. The localized high gradient of the NSF represents a critical feature, which the standard AFM version fails to represent appropriately: the area where the NSF decreases rapidly is discretized without any refinement. Pseudo-islands are therefore introduced and added to the advancing front. In case of sudden NSF variations, it is beneficial for a correct discretization to place pseudo-islands in correspondence of the local minima. When the front searching option is chosen (point d.1 of the MESHGR
154
8.2 Test cases for idealized geometries
Figure 8.1: Test 1. (a) NSF defined over a squared domain with a background value of 22000 m and a sudden decrease in the center to a value of 2500 m; (b) a longitudinal domain transect of the NSF.
Figure 8.2: Test 1. Mesh T1 A, produced with the circulating option and mesh T1 B, produced with the front searching option of the standard AFM implemented in MESHGR, for the NSF of Figure 8.1.
155
8. PSEUDO-ISLANDS
pseudo-code), the elements will start to be created around the pseudo-islands, thus giving priority to the most refined parts of the mesh. The following pseudo-code is a possible, straightforward solution for the introduction of pseudo-islands. The sequence of instructions of the main algorithm is unaltered until the end of the boundary discretization (point c.2.5), where Γ0 is initialized: p.1) define a threshold condition based on the NSF. Call it ρM IN ; p.2) for i = 1, NPOIG (number of nodes of B) check if: p.2.1) i falls inside the domain, at a distance from the boundary higher than 2ρi (the coefficient 2 creates a buffer distance between potential pseudo-islands and the domain boundaries); p.2.2) ρi is lower than αρM IN , with α ≤ 1; If conditions p.2.1) and p.2.2) are met, save node i in the dedicated array [POINTSIN]; p.3) order [POINTSIN] in increasing order of ρi ; p.4) for i = 1, dimension [POINTSIN] check if, for all points j for which a pseudoisland already exists, the distance ij > k(ρi + ρj ). The coefficient k ≥ 1 creates a buffer distance between neighbour pseudo-islands. If true:
p.4.1) a pseudo-island is created with side length equal to ρi and with barycentre the background mesh node i; p.4.2) the sides and vertices of the pseudo-island are added to Γ0 ; p.5) choose a departure zone option (front searching or circulating): for more than one pseudo-islands, the searching option is recommended; p.6) follows the sequence of instruction d.2) to d.15). p.7) at the end of the meshing process, the triangles corresponding to the pseudoislands are included as elements and nodes of the final mesh. Figure 8.3 visualizes successive stages in the domain discretization according to the pseudo-island front searching option. All pseudo-islands are identified before the actual spatial discretization starts. The circulating option is suggested when a unique
156
8.2 Test cases for idealized geometries
(a)
(b)
(c)
(d)
(e)
(f )
Figure 8.3: Successive stages of the pseudo-island front searching algorithm: (a) discretization of the boundary domain and identification of the NSF minima; (b) construction of the pseudo-islands and generation of the first element; (c) to (e) built-up of the mesh around the pseudo-islands and merging of the internal fronts; (f) completion of the domain discretization and insertion of the pseudo-islands in the final mesh.
pseudo-island is created, while the front searching option is suitable to make the sides of all pseudo-islands active. With this pseudo-code implemented inside MESHGR, the discretization of the NSF of Figure 8.1 produces the mesh T1 C (Figure 8.4). T1 C has 106 nodes and 178 elements. The algorithms identifies 4 pseudo-islands, 1 corresponding to the points where ρi < ρM IN . The 4 pseudo-islands are added at last 2 to the element list. A major limitation of T1 C is however the presence of distorted elements and abrupt area changes, especially between elements located at the transition between the high and low NSF gradients. In fact the relation controlling the element side is equation (5.7), which produces 29 couples of elements with an area ratio higher than 0.5, while 3 nodes have a valence higher than 8. It has been found that decreasing the upper limit regulating the side transition between adjacent elements, as given by the following formula, would produce a smoother element area grading (mesh T1 D in
157
8. PSEUDO-ISLANDS
Figure 8.4: Test 1. Pseudo-island front searching algorithm applied to the NSF of Figure 8.1. Mesh T1 C, on the left, obtained for minimum/maximum element side transition coefficients of 0.55 and 2; and Mesh T1 D, on the right, for minimum/maximum element side transition coefficients of 0.55 and 1.15. 4 pseudo-islands are evidenced.
Figure 8.4): ρN
ρN = 1.15AB 0.55AB
if 0.55AB ≤ ρAB ≤ 1.15AB if ρAB > 1.15AB if ρAB < 0.55AB.
(8.1)
Including the 4 pseudo-islands, mesh T1 D has 300 nodes and 570 elements, around three times more than T1 C. While the number of pseudo-islands is the same, here the number of nodes with a valence higher than 8 is reduced to 0 and the couples of elements with an abrupt area transition (area ratio higher than 0.5) from 29 to 2, 1 of which having a node belonging to the domain boundary and then set prior to the calculation of the pseudo-islands. A useful measure of the geometrical quality of the mesh elements is the index qi (formula 3.13), which calculates how equilateral is a triangle. The distribution of qi for T1 C and T1 D (Figure 8.5) demonstrates that a more controlled area grading in correspondence of high NSF variations prevents the formation of highly distorted elements (7 for T1 C and none for T1 D), while assuring nearly equilateral shapes for the majority of the elements.
158
8.2 Test cases for idealized geometries
(a)
(b)
Figure 8.5: Test 1. qi distribution into classes for T1 C (left side) and T1 D (right side). Only the range of values [0.6 - 1] is represented.
8.2.2
Test case 2
A second test case presents a target NSF with linear and curvilinear shock-like features (Figure 8.6). Also in this situation the standard MESHGR version is not able to reproduce the rapidly changing NSF values (mesh T2 A in Figure 8.7a). On the contrary, the pseudo-island algorithm, after detecting the positions of the NSF minima, can associate pseudo-islands to each of them, from where elements are formed (mesh T2 B and mesh T2 C in Figures 8.7b and 8.7c). Mesh T2 B is generated according to the side transition relation (5.7), while mesh T2 C with the stricter side limitation of formula (8.1).
Mesh T2 B has 472 nodes and 910 elements, 31 of which are added to the
final mesh as pseudo-islands. Mesh T2 C has 1129 nodes and 2224 elements, with the same number of pseudo-islands of T2 B. As observed for test case 1, the pseudo-island algorithm with the original side change condition between adjacent elements does not assure an acceptable area grading and well-shaped elements: T2 B contains 182 couples of adjacent elements with an area ratio higher than 0.5, and 6 nodes with a valence number higher than 7; for T2 C instead the abrupt area transitions are reduced to 38, 30 of which have at least 1 node belonging to the boundary (and thus with a position independent of the discretization of the inner domain), while the number of nodes with a valence higher than 7 reduces to 0. The distribution of the element shape index qi (Figure 8.8) confirms an enhanced mesh quality for T2 C : the 11 elements with a qi < 0.6 are located along the domain boundary, where the node location is fixed at the initialization of the front and remains unchanged during the meshing.
159
8. PSEUDO-ISLANDS
Figure 8.6: Test 2. NSF defined over a rectangular domain with one linear and one curvilinear high gradient areas.
8.3
Application of pseudo-islands to a LIDAR derived digital elevation model
In this section the potential improvements brought by the pseudo-island concept are quantitatively evaluated for real case applications. A 2 m resolution LIDAR survey of the Portuguese coastline provided by the Portuguese Geographic Institute (IGP), retaining information of the bare Earth surface without inclusion of buildings and vegetation, is used as the reference DEM. LIDAR data have been increasingly used in the current coastal modelling practice, as they incorporate features to an unprecedented level of detail. This property becomes decisive when the flow interaction with the bottom and with significant hydraulic structures is modelled for areas exposed to tide and storm surge flooding. However, elevation error arise whenever the DEM is discretized into a mesh: a proper discretization should hence be targeted to minimize the elevation errors, while at the same time maintain the element quality and sizing properties which make the mesh usable in finite element applications.
160
8.3 Application of pseudo-islands to a LIDAR derived digital elevation model
Figure 8.7: Test 2. Meshes obtained for the NSF of Figure 8.6: (a) Mesh T2 A, created without pseudo-islands; (b) Mesh T2 B, obtained with the pseudo-island option and minimum/maximum side transition coefficients of 0.55 and 2; (c) Mesh T2 C, obtained with the pseudo-island option and minimum/maximum side transition coefficients of 0.55 and 1.15.
(a)
(b)
Figure 8.8: Test 2. qi distribution into classes for T2 B (left side) and T2 C (right side). Only the range of values [0.4 - 1] is represented.
161
8. PSEUDO-ISLANDS
Figure 8.9: C´ avado estuary, characterized by a breakwater separating the river course from the low-lying marsh.
8.3.1
Study sites
Two areas of the LIDAR data presenting some typical features of coastal environments are identified as study cases. The first one is located on the estuary of the C´avado river, near the town of Esposende, in the North of Portugal (Figure 8.9). The site is characterized by an extended breakwater delimiting the right side of the main river course from the backside low lying intertidal marsh, with a series of openings regulating the mass exchange between the two water bodies. The structure was designed to dissipate high energetic waves entering the river mouth. The bathymetric gradients of the breakwater, the river and the marsh bottom irregularities present a recurrent discretizazion challenge for the coastal modeller. The second area is the region of the Ria Formosa, in the South of Portugal, an extended barrier lagoon connected to the ocean by a series of morphologically active inlets (Figure 8.10a). The interior of the lagoon is a marsh with a complex network of tidal channels and sand-mud banks,
162
8.3 Application of pseudo-islands to a LIDAR derived digital elevation model
Figure 8.10: View of the Ria Formosa, with the delimitations of the selected study sites.
resulting in local high bathymetric gradients. In particular, two sites have been selected as representative of the lagoon morphology (Figure 8.10b): site 1 has a triangular shape and encompasses a major meandering tidal channel connected to a series of minor channels. Abrupt local bathymetry changes are observed between the fairly flat bank surfaces and the channel bottom. A good description of theses irregularities is essential for assuring a stable and reliable simulation of the wet and dry process. Site 2 includes the extremities of two barrier islands, the inlet between them and part of the salt marsh behind.
8.3.2
Methodology for the mesh accuracy assessment
In order to assess the accuracy of a mesh in reproducing the underlying bathymetry, expressed as a scatter set of points, the LIDAR dataset is subdivided randomly into a training and a test scatter set, containing respectively the 90% and the 10% of the points. It is believed that, with relation to the 2 m DEM resolution, the training scatter set is dense enough not to deteriorate the bathymetric/topographic information, while
163
8. PSEUDO-ISLANDS
the remaining 10% of the points represent a sufficient population for a sound statistical analysis. After one or a series of meshing criteria are defined analytically, a final NSF is calculated and read in the AFM generator for performing the discretization. The training scatter set is interpolated to the mesh, which in its turn is interpolated to the points of the test scatter set. The elevation error is defined as the absolute difference between zitest (x, y), elevation of the test point, and zimesh (x, y), the interpolated elevation of the mesh in the same point. All interpolations are linear. The test points are assumed to be the “true” values of the terrain elevation. For the statistical description of the error, two estimators are considered: the mean average error (MAE), calculated as: N 1 X test M AE = |zi (x, y) − zimesh (x, y)| N
(8.2)
i=1
and the root mean square error (RMSE): "
N 2 1 X test mesh RM SE = zi (x, y) − zi (x, y) N i=1
#0.5
(8.3)
where N is the total number of test points.
8.3.3
Hydrodynamic model description
The meshes RF1 noPI and RF1 PI, generated respectively without and with pseudoislands for site 1 of the Ria Formosa (see section 8.4.2), are incorporated into a wider mesh encompassing a larger part of the lagoon, of the barrier islands and of the nearby coastal zone (Figure 8.11). A wet and dry tide modelling of the Ria is carried out with the aim of assessing to which extent the meshes influence the reproduction of the wetted areas. The outer mesh, generated through the SMS interface (Surface-water Modeling System, 2013), is generated according to the wavelength to mesh size ratio (Le Provost and Vincent, 1986; Westerink et al., 1992a) for the ocean part; inside the Ria two NSFs are applied: the first is the module of the inverse of the bathymetry/topography gradient (formula 7.3); the second is based on the channel geometry with an element side length ranging from 50 m along the axis of the main channels to 300 m at a distance of 700 m or more from the channel axis.
164
8.3 Application of pseudo-islands to a LIDAR derived digital elevation model
Figure 8.11: Extension of the Ria Formosa domain used for a wet and dry tide model application. The location of site 1 is evidenced as well.
Out of the study site 1, the bathymetry has been digitized from the 1: 15 000 and 1: 150 000 nautical charts. As nautical charts may not provide a sufficient characterization of the lagoon morphology and bathymetry, the main objective of this model application cannot be a validation study of the hydrodynamics of the lagoon, but rather a test of the behaviour of a wet and dry algorithm, under the same model conditions, for different AFM discretization approaches. The 2D finite element circulation model ADCIRC (Luettich and Westerink, 2004) is used, which solves the depth integrated shallow water equations, under the hydrostatic pressure approximation, according to the generalized wave continuity equation formulation (Luettich et al., 1992). The governing equations are written in formulas (7.7-7.9). The quadratic parametrization of the bottom friction is adopted (formula 7.10), with a constant drag coefficient of 2.5 ∗ 10−3 , a value chosen
in accordance with previous tide modelling applications (Blain and Rogers, 1998). 13 astronomical components (K1, K2, L2, M2, MU2, N2, 2N2, NU2, O1, P1, Q1, S2 and T2), extracted from the tidal database of Le Provost et al. (1998), have been imposed to the domain ocean boundary, while the potential functions of the same components have been applied to the interior nodes. A wet and dry algorithm (Luettich and Westerink, 1995; Luettich and Westerink, 1999) is included, which calculates the extent of the flooded areas on an element-by-element basis, according to the local water elevation
165
8. PSEUDO-ISLANDS
and water elevation gradient. An 8 day long simulation, from the 2nd January to the 10th January 2013, is successfully performed for assessing the computational stability of the proposed discretization and for representing the extension of the wet and dry areas during different tidal phases. The same simulation is then prolonged to a total of 30 days for including a complete spring-neap tide cycle, from which a harmonic analysis of the water elevation results is extracted for comparison with the measured astronomical constituents at the Faro-Olh˜ ao tide gauge (Figure 8.10).
8.4 8.4.1
Results and discussion C´ avado estuary
Considering the continuous slope variations between river bed, breakwater and adjacent areas, a target NSF is set equal to the inverse of the module of the bathymetry gradient ∇h:
∆x = (∇h)−1
(8.4)
Relashionship (8.4) assures the placement of smaller elements where the bathymetric variations are higher, while elements in flat zones remain coarse. Keeping the same NSF condition and the element sizing transition formula (8.1), three different discretizations are performed by means of MESHGR: mesh CA PIa, where pseudo-islands are inserted if the local bathymetric slope exceeds 15◦ and with the coefficient of minimum distance between neighbouring pseudo-islands k set to 1.5 (point p.4 of the pseudo-code); mesh CA PIb, again with a threshold slope of 15◦ and k equal to 3; and mesh CA noPI, without pseudo-islands. Meshes CA PIb and CA noPI are visualized in Figure 8.12. Both discretizations evidence the presence of the breakwater and the location of its internal opening. Zones of CA PIb with a denser discretization than CA noPI, along the breakwater and on the estuary bottom, are attributable to the insertion of pseudoislands, which improve the local terrain representation. Distorted elements along the domain boundary are attributed to the fact that MESHGR creates the nodes along the boundaries before the placement of the pseudo-islands and the internal nodes: this sequence may lead to undesired shapes when the area transition coefficient is small (1.15 in this case) and the domain has a limited number of elements in comparison
166
8.4 Results and discussion
Figure 8.12: (a) C´ avado estuary bathymetry and discretization domain; (b) mesh CA PIb, created with pseudo-islands; (c) mesh CA noPI, created without pseudo-islands.
with its total dimension. Table 8.1 summarizes the accuracy properties of the three meshes. Table 8.1: Mesh statistics for the discretization tests of the C´avado estuary.
The statistics of the absolute vertical elevation error are sensitive to the number of pseudo-islands introduced into the mesh: for a higher number of pseudo-islands, the errors decrease. It can be inferred that the opening of internal fronts, in correspondence of the local NSF minima, is beneficial to the discretization process as it increases the accuracy of the mesh in reproducing the underlying DEM. This result can have a twofold interpretation: locally, the presence of a pseudo-island forces the algorithm to place an element precisely where the dicretization condition (in this case the bathymetry gradient) is more demanding, and not in its neighbourhood as in a standard implementation
167
8. PSEUDO-ISLANDS
of the advancing front (an exaggeration of this effect has been demonstrated in the idealized test cases 1 and 2). From a global point of view, the priority in meshing denser areas inside the domain, instead of starting the discretization near the boundaries, increases the total number of nodes and elements, and, as a consequence, the underlying surface is in general better represented.
8.4.2
Ria Formosa: sites 1 and 2
For site 1 two meshes are produced: RF1 noPI and RF1 PI, without and with pseudoislands respectively. Both discretizations are performed according to the following NSF: ∆x = 3(∇h)−1
(8.5)
The empirical coefficient of relationship (8.5) has been chosen since, after a series of preliminary tests, it demonstrated to produce a mesh with a suitable range of element sizes to be included in the wet and dry tide modelling of the Ria Formosa. For site 2, four different algorithm configurations have been tested, with emphasis in reproducing sudden bathymetric changes according to relashionship (8.4). As a benchmark case, the first created mesh, RF2 noPI, does not include pseudo-islands. The second attempt introduces pseudo-islands in correspondence of bathymetric slopes more inclined than 20◦ (mesh RF2 PIa) and setting the buffer coefficient k equal to 3. In the third case the threshold gradient was decreased to 10◦ and k was kept equal to 3, producing mesh RF2 PIb. In the fourth test the threshold gradient is 10◦ and k is decreased to 1.5 (mesh RF2 PIc), in order to increase the pseudo-island distribution. Table 8.2 summarizes the accuracy properties of the meshes for sites 1 and 2. The observations drawn for the series of meshes of the C´avado estuary are confirmed also for sites 1 ans 2 of Ria Formosa. The most evident result is the improvement of the overall mesh accuracy when pseudo-islands are inserted into the discretization: the higher is the number of pseudo-islands, the lower is the elevation MAE and RMSE.
8.4.3
Tidal simulation
Table 8.3 compares the amplitude and phase of the measured tidal constituents for the Faro-Olh˜ ao gauge with the computed values of the 30 day long tide simulation at the same site. The harmonic analysis of the modelled results show the same output
168
8.4 Results and discussion
Table 8.2: Main properties of the meshes for the selected sites of the Ria Formosa.
values independently of using RF1 noPI or RF1 PI into the wider mesh. The harmonic analysis of the measured values is provided by the tidal charts of the Portuguese Hydrographical Institute (Instituto Hidrogr´afico, 2013). The model performance in reproducing the tidal amplitudes at the principal inlet of the lagoon is optimal, and acceptable values are achieved for the phases. Although no other nearby tidal gauges are available, it can be inferred that, considering the relatively limited domain extension and the similar hydrodynamic conditions, the present model configuration is able to reproduce the tide elevation conditions regulating the exchanges of water between the inlets and the open ocean. The reproduction of the flooding and drying cycle inside the Ria is highly dependent on the quality and resolution of the DEM interpolated to the mesh. As part of the lagoon bathymetry has been digitized and may lack the resolution needed, the analysis of the wet and dry process will be limited to site 1, for which the LIDAR data has been used. Figure 8.13 compares the wetted areas when using mesh RF1 noPI and RF1 PI, for the same computational time step corresponding to a flood phase. It is noticeable that, in the discretization with pseudo-islands, the wetted area assumes a more continuous and precise pattern: in particular the medium and small tidal channels, whose sides are evidenced by the depth 0 relative to MSL, are connected to the bigger channel on the right side of the picture, on the contrary of what occurs for RF1 noPI, and, as a result, the flooding phase is reproduced more faithfully. The barrier island, on the bottom left of the picture, is instead well delineated in both meshes due to its bigger scale size.
169
8. PSEUDO-ISLANDS
Figure 8.13: Comparison between meshes and the wetted areas for the study site 1 of the Ria Formosa during flood tide, at 4.30am of January 2013: (a) mesh RF1 noPI, created without pseudo-islands; (b) wetted area for RF1 noPI; (c) mesh RF1 PI, created with pseudo-islands; (d) wetted area for RF1 PI . The dark lines correspond to the bathymetric depth 0 with respect to MSL.
170
8.5 Conclusions
Table 8.3: Measured and computed tidal constituent values for the Faro-Olh˜ao gauge
8.5
Conclusions
This chapter has presented an alternative implementation of the advancing front method for the meshing of 2D domains. The new version introduces internal openings, called pseudo-islands, which are added to the boundary active front. The placement of the pseudo-islands in correspondence of the specific conditions derived from the node spacing function, i.e. the local minima or a threshold gradient, gives priority to the discretization of the zones with the highest resolution requirements, while the front propagates in an inside to outside direction. When creating a new node, a stricter side length condition between the new and the existing element must be introduced for obtaining smooth area transitions and fairly equilateral triangles. Laplacian smoothing and diagonal swapping post-processing routines are applied as well for achieving an overall good mesh quality. Discretization tests, performed on idealized geometries and on real domains described by a high resolution digital elevation model, demonstrated the advantages of using pseudo-islands, both in terms of discretization performance than of modelling of inundation processes. The higher is the number of pseudo-islands introduced in the meshing process, the higher is the number of active fronts and the larger is the final mesh. The increased discretization accuracy favours the reproduction of irregular bathymetries, small scale hydraulic features and the realistic modelling of flooding processes. The method may include, in presence of pseudo-islands, the redistribution the boundary nodes in order to improve the element shape along the domain boundary and the application of the pseudo-islands to an a-posteriori, error-based node spacing function.
171
9
CONCLUSIONS 9.1
Concluding remarks
In this research an advancing front method algorithm (MESHGR) has been improved and applied to the meshing of idealized and real shallow water domains. The passage from a static to a dynamic memory allocation, the automation of the creation of the input files, the resolution of a series of bugs, the addition of a basic graphical user interface, have made the code more robust and able to cope with multi-connected, multi-scale geometries and with node spacing functions derived from multiple criteria. A series of export options have been added for writing the final mesh into the mesh format of finite element models widely used in the ocean and coastal research. It has been demonstrated the effectiveness of the post-processing routines in enhancing the quality of the isotropic mesh: diagonal swapping is specifically suited to lower the overall maximum node valence; Laplacian smoothing, combined with the previous subroutine, decreases the number of abrupt area changes between adjacent elements. A critical factor on which the accordance between the NSF and the actual size/shape of the final mesh elements depends is the level of refinement of the background mesh. The background node positioning should be done with care: higher refinement is needed where more control is wanted on the element size, shape and orientation, especially where the node spacing function contours are irregular. Various a-priori and a-posteriori criteria, for isotropic and anisotropic conditions, have been tested. The introduction of a distortion factor for a more realistic reproduction of the wetted front has turned out to be feasible, although limited improvements in
173
9. CONCLUSIONS
the rendering of the front have been obtained with respect to an isotropic mesh which already assure a sufficient refinement of the intertidal areas. Moreover, anisotropic elements exceeding a certain elongation factor expose the model to computational instabilities. The multi-criteria meshing method has demonstrated to be a meaningful and comprehensive approach for meshing domains varying over several spatial scales and with several input forcing functions, as often present in shallow water models and storm surge inundation studies. The simultaneous consideration of carefully chosen criteria, based on the physics and the geometry of the flow problem and on the computational constraints of the model, overcome the limitations presented by applying just a single criterion, or the necessity of refinements by hand. The modeller’s sensitivity and previous experience of similar study cases are however necessary for an effective choice of the criteria, while the algorithm capacity to design well graded meshes with high quality properties is a prerequisite of a good model application. Within the framework of the multi-criteria method, alternative approaches have been assessed for tide modelling with inclusion of the wet and dry process. In particular, the a-posteriori LTEA+CD method has been applied to such domains in order to study the grade of convergence of a-priori meshes. Finally, the pseudo-island option introduced in the standard structure of the advancing front method has shown good results in the discretization of both idealized and real test cases, with node spacing functions derived from bathymetries with pronounced gradients. The introduction of internal fronts, combined with a more controlled area transition between elements, has allowed a more accurate representation of irregular bathymetries. Comparative tests demonstrate that, for the same input node spacing function, the presence of pseudo-islands placed over the highest gradients decrease the mesh average vertical elevation error, in comparison with a mesh obtained by the standard advancing front method. A pseudo-island derived mesh applied to a wet and dry study of the Ria Formosa has shown a more realistic representation of the wetted area extension than a mesh created, for the same node spacing function, without pseudo-islands.
174
9.2 Further work
9.2
Further work
Several aspects of this research offer possibilities of further development, both on the programming and on the modelling side. In relation to the MESHGR algorithm, the memory storage allocated to the output mesh, currently set to 106 nodes and 2 ∗ 106
elements, should be increased by at least a factor of 5, in order to make the code potentially able to reproduce the largest meshes reported until now, in the author’s knowledge, in 2D shallow water modelling studies. More effort should be spent in the development of the graphical user interface, to be re-converted to an open-source library, and provided with other essential options, as the possibility to design the background mesh and the domain input boundary, and to calculate the node spacing functions (in the present version this is done with the help of the SMS interface). It is desirable to couple the multi-criteria approach with a-posteriori algorithms, as
the LTEA+CD, for addressing the issue of mesh optimization and de-refinement. This aspect is essential for the implementation of quasi-real time storm surge forecasting systems. As the presented modelling applications always relied on ADCIRC, it would be useful to include alternative 2D shallow water models and wet and dry algorithms. A comparative analysis between different models would give insight into the nature of the systematic errors observed, for example, along the Lima estuary, and to separate the errors caused by the discretization and by the model itself in reproducing field data. Experiences of this type are already being undertaken in collaboration with other universities. The pseudo-island algorithm has proven to bring advantages in the discretization of complex, high resolution real bathymetries. However other tests are needed: application to a broader number of domains, bathymetries, meshing criteria and threshold conditions will corroborate the validity of the method. Also, new versions of the pseudo-island algorithm may receive attention, where, for example, the pseudo-islands are created during the discretization process and not all together before the formation of the first mesh element.
175
Bibliography [1] Aquaveo, 2000. SMS: Surface-Water Modeling System, Version 10.1, Reference manual. Accessed 1.4.2011. URL http://www.xmswiki.com/xms/SMS:SMS [2] Ara´ ujo, M. A. V. C., Trigo-Teixeira, A., Mazzolari, A., 2011a. Simulation of storm surge events at the Portuguese coast (Viana do Castelo). Littoral 2010 - Adapting to Global Change at the Coast: Leadership, Innovation, and Investment, Article number 120040. URL http://dx.doi.org/10.1051/litt/201112004 [3] Ara´ ujo, M. A. V. C., Trigo-Teixeira, A., Mazzolari, A., 2011b. Wave set-up in the modeling of storm surge at Viana do Castelo (Portugal). Proceedings of the 11th International Coastal Symposium. Journal of Coastal Research SI 64(1), 971–975. [4] Aretxabaleta, A., 2003. Pelosab: Finite element mesh generation tutorial: using BatTri. Numerical Methods Laboratory Report NML-03-14, Dartmouth College, NH. [5] Atkinson, J. H., Westerink, J. J., Luettich, R. A., 2004. Two-dimensional dispersion analyses of finite element approximations to the shallow water equations. International Journal for Numerical Methods in Fluids 45(20), 715–749. [6] Bacopoulos, P., Parrish, D. M., Hagen, S. C., 2011. Unstructured mesh assessment for tidal model of the South Atlantic Bight and its estuaries. Journal of Hydraulic Research 49, 487–502. [7] Baehmann, P. L., Wittchen, S. L., Shephard, M. S., Grice, K. R., Yerry, M. A., 1987. Robust, geometrically based, automatic two-dimensional mesh generation. International Journal for Numerical Methods in Engineering 24(8), 1043–1078. [8] Bates, P. D., Marks, K. J., Horritt, M. S., 2003. Optimal use of high-resolution topographic data in flood inundation models. Hydrological Processes 17, 537–557. [9] Bilgili, A., Smith, K. W., Lynch, D. R., 2005. BatTri: A 2D unstructured grid generator for finite element circulation modeling. Computers & Geosciences 32(5), 632–642.
177
BIBLIOGRAPHY
[10] Blain, C. A., Linzell, R. S., Chu, P., Massey, C., 2010. Validation test report for the Advanced CIRCulation Model (ADCIRC) v45.11. NRL Memorandum Report, NRL/MR/7320–10-9205, Naval Research Laboratory, Washington, DC. [11] Blain, C. A., Rogers, W. E., 1998. Coastal tide prediction using the ADCIRC-2DDI hydrodynamic finite element model: model validation and sensitivity analyses in the Southern North Sea/English Channel. Naval Research Laboratory, Oceanography Division, Stennis Space Center, MS. Report number NRL/FR/7322–98-9682. [12] Blain, C. A., Westerink, J. J., Luettich, R. A. J., 1994. The influence of domain size on the response characteristics of a hurricane storm surge model. Journal of Geophysical Research 99, 18467–18479. [13] Blain, C. A., Westerink, J. J., Luettich, R. A. J., 1998. Grid convergence studies for the prediction of hurricane storm surge. International Journal for Numerical Methods in Fluids 26, 369–401. [14] Blanton, B. O., Luettich, R. A., 2008. North Carolina coastal flood analysis system model grid generation. RENCI Technical Report Series TR-08-05, 18p. [15] Borouchaki, H., Frey, P. J., 1998. Adaptive triangular-quadrilateral mesh generation. International Journal for Numerical Methods in Engineering 41, 915–934. [16] Borouchaki, H., George, P. L., Hecht, F., Laug, P., Saltel, E., 1997a. Delaunay mesh generation governed by metric specifications. Part I: Algorithms. Finite Elements in Analysis and Design 25, 61–83. [17] Borouchaki, H., George, P. L., Mohammadi, B., 1997b. Delaunay mesh generation governed by metric specifications. Part II: Application examples. Finite Elements in Analysis and Design 25, 85109. [18] Bossen, F. J., Heckbert, P. S., 1996. A pliant method for anisotropic mesh generation. 5th International Meshing Roundtable, Sandia National Laboratories, 63–76. [19] Bunya, S., Dietrich, J. C., Westerink, J. J., Ebersole, B. A., Smith, J. M., Atkinson, J. H., Jensen, R., Resio, D. T., et al., 2010. A high-resolution coupled riverine flow, tide, wind, wind wave, and storm surge model for Southern Louisiana and Mississippi. Part I: model development and validation. Monthly Weather Review 138, 345–377. [20] Bykat, A., 1983. Design of a recursive, shape controlling mesh generator. International Journal for Numerical Methods in Engineering 19, 1375–1390. [21] Carey, G., 1995. A posteriori error estimation and mesh refinement. In: Lynch, D. R., Davies, A. M. (Eds.), Quantitative Skill Assessment for Coastal Ocean Models. Coastal and Estuarine Series. American Geophysical Union, Washington DC, 47, 15–29.
178
BIBLIOGRAPHY
[22] Castro-Diaz, M. J., Hecht, F., 1995. Anisotropic surface mesh generation. INRIA Research Report. RR-2672. [23] Chow, V. T., 1981. Development of uniform flow and its formulas. Mc Graw Hill (Ed.) Open Channel Hydraulics. International Student edition, 17th Printing, 88-123. [24] Cialone, M. A., Militello, A., Brown, M. E., Kraus, N. C., 2002. Coupling of wave and circulation numerical models at Grays Harbor entrance, Washington, USA. Proceedings of the 28th International Conference on Coastal Engineering, 1279–1291. [25] Cobby, D. M., Mason, D. C., Horritt, M. S., Bates, P. D., 2003. Two-dimensional hydraulic flood modelling using a finite-element mesh decomposed according to vegetation and topographic features derived from airborne scanning laser altimetry. Hydrological Processes 17, 1979–2000. [26] Compass Hydrographic Services, 2006. Measurements of current and water levels in the Lima River, Portugal. For Hidrodata Lda. Rep Number 06/11/10048/1, 17pp. [27] Cougny, H. L. D., Shepard, M. S., 1999. Parallel unstructured grid generation. In: Handbook of Grid Generation. CRC Press, Boca Raton, FL, 24.10-24.16. [28] Courant, R., Friedrichs, K., Lewy, H., 1967. On the partial difference equations of mathematical physics. IBM Journal of Research and Development 11, 215–234. [29] Cuvelier, C., Segal, A., van Steenhoven, A. A., 1986. Finite element method and NavierStokes equations. Reidel Publishing Company, Dordrecht, 483pp. [30] Delaunay, B. N., 1934. Sur la sphere vide. Bulletin de l’Acad´emie des Sciences de l’URSS. Classe des sciences math´ematiques et naturelle 7, 793–800. [31] Demirbilek, Z., Panchang, V. G., 1996. CGWAVE: A coastal surface water wave model of the mild slope equation, Technical Report CHL-98-26. US Army Corps of Engineers, Waterways Experiment Station, Vicksburg, MS 39180. [32] Dietrich, J. C., Kolar, R. L., Luettich, R. A., 2004. Assessment of ADCIRCs wetting and drying algorithm. In: Proceedings 15th International Conference on Computional Methods in Water Resources Vol. 2, Miller C. T. et al. (Eds.) Elsevier, Amsterdam. [33] Dietrich, J. C., Westerink, J. J., Kennedy, A. B., Smith, J. M., Jensen, R. E., Zijlema, M., Holthuijsen, L. H., Dawson, C., Jr., R. A. L., Powell, M. D., Cardone, V. J., Cox, A. T., Stone, G. W., Pourtaheri, H., Hope, M. E., Tanaka, S., Westerink, L. G., Westerink, H. J., Cobell, Z., 2011. Hurricane Gustav (2008) waves and storm surge: hindcast, synoptic analysis, and validation in Southern Louisiana. Monthly Weather Review 139, 2488–2522. [34] Dompierre, J., Vallet, M. G., Labbe, P., Guibault, F., 2005. An analysis of simplex shape measures for anisotropic meshes. Computer Methods in Applied Mechanics and Engineering 194, 4895–4914.
179
BIBLIOGRAPHY
[35] Donnel, B. P., Letter, J. V., McAnally, W. H., et al., 2006. Users guide to RMA2 WES Version 4.5. U.S. Army Engineering Research and Development Center, Waterways Experiment Station, Coastal and Hydraulics Laboratory. [36] Ebersole, B. A., Westerink, J. J., Bunya, S., Dietrich, J. C., Cialone, M. A., 2010. Development of storm surge which led to flooding in St. Bernard Polder during Hurricane Katrina. Ocean Engineering 37, 91–103. [37] El-Hamalawi, A., 2004. A 2D combined advancing front- Delaunay mesh generation scheme. Finite Elements in Analysis and Design 40, 967–989. [38] Feyen, J. C., Atkinson, J. H., Westerink, J. J., 2002. GWCE-based shallow water equation simulations of the Lake PontchartrainLake Borgne tidal system. Computational Methods in Water Resources, II. Vol. 47. S. M. Hassanizadeh et al., (Eds.), Developments in Water Science, Elsevier. [39] Field, D. A., Frey, W. H., 1995. Structural improvement of planar triangulations: some constraints and practical issues. Communications in Numerical Methods in Engineering 11, 191–198. [40] Flor, A. P., Mazzolari, A., Gon¸calves, A. B., Ara´ ujo, M., Trigo-Teixeira, A., 2013. Influence of elevation modelling on hydrodynamic simulations of a tidally-dominated estuary. Journal of Hydrology 497, 152–164. [41] Fortunato, A. B., Bruneau, N., Azevedo, A., Ara´ ujo, M. A. V. C., Oliveira, A., 2011. Automatic improvement of unstructured grids for coastal simulations. Proceedings of the 11th International Coastal Symposyum. Journal of Coastal Research SI 64, 1028–1032. [42] Franca, A. S., Haghighi, K., Oliveira, L. S., 1995. Error estimation and adaptivity in finite element analysis of transport problems. In R.W. Johnson and E.D. Hughes (Eds.). Quantification of Uncertainty in Computational Fluid Dynamics 213, 19-24. [43] French, J. R., 2010. Critical perspectives on the evaluation and optimization of complex numerical models of estuary hydrodynamics and sediment dynamics. Earth Surface Processes and Landforms 35, 174–189. [44] Frey, P. J., George, P. L., 2000. Mesh generation: application to finite elements. Hermes Science. [45] Frey, W. H., 1987. Selective refinement: a new strategy for automatic node placement in graded triangular meshes. International Journal for Numerical Methods in Engineering 24, 2183–2200. [46] Garimella, R. V., Shephard, M. S., 2000. Boundary layer mesh generation for viscous flow simulations. International Journal for Numerical Methods in Engineering 49, 193–218.
180
BIBLIOGRAPHY
[47] George, P. L., Seveno, E., 1994. The advancing-front mesh generation revisited. International Journal for Numerical Methods in Engineering 37, 3605–3619. [48] Geuzaine, C., Remacle, J. F., 2009. Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities. International Journal for Numerical Methods in Engineering 79, 1309–1331. [49] GFGEN, 2001. User guide version 4.35. US Army Engineer Research and Development Center. Waterways Experiment Station Coastal and Hydraulics Laboratory. [50] Gorman, G. J., Piggott, M. D., Pain, C. C., de Oliveira, C. R. E., Umpleby, A. P., Goddard, A. J. H., 2006. Optimisation based bathymetry approximation through constrained unstructured mesh adaptivity. Ocean Modelling 12, 436–452. [51] Gray, W., Lynch, D., 1977. Time-stepping schemes for finite element tidal model computations. Advances in Water Resources 1, 83–95. [52] Greenberg, D. A., Dupont, F., Lyard, F. H., Lynch, D. R., Werner, F. E., 2007. Resolution issues in numerical models of oceanic and coastal circulation. Continental Shelf Research 27, 1317–1343. [53] Hagen, S. C., Horstmann, O., Bennett, R. J., 2002. An unstructured mesh generation algorithm for shallow water modeling. International Journal of Computational Fluid Dynamics 16, 83–91. [54] Hagen, S. C., Parrish, D. M., 2004. Unstructured mesh generation for the western North Atlantic tidal model domain. Engineering with Computers 20, 136–146. [55] Hagen, S. C., Westerink, J. J., Kolar, R. L., 2000. One-dimensional finite element grids based on a localized truncation error analysis. International Journal for Numerical Methods in Fluids 32, 241–262. [56] Hagen, S. C., Westerink, J. J., Kolar, R. L., Horstmann, O., 2001. Two-dimensional, unstructured mesh generation for tidal models. International Journal for Numerical Methods in Fluids 35, 669–686. [57] Hagen, S. C., Zundel, A. K., Kojima, S., 2006. Automatic, unstructured mesh generation for tidal calculations in a large domain. International Journal of Computational Fluid Dynamics 20, 593–608. [58] Hannah, C. G., Wright, D. G., 1995. Depth dependent analytical and numerical solutions for wind-driven flow in the coastal ocean. In: Lynch, D. R., Davies, A. M. (Eds.), Quantitative Skill Assessment for Coastal Ocean Models. Vol. 47 of Coastal Estuarine Studies. American Geophysical Union, Washington D.C., 125-152.
181
BIBLIOGRAPHY
[59] Hasan, G. M. J., van Maren, D. S., Cheong, H. F., 2012. Improving hydrodynamic modeling of an estuary in a mixed tidal regime by grid refining and aligning. Ocean Dynamics 62, 395–409. [60] Hassan, O., Probert, E. J., Morgan, K., Peraire, J., 1995. Mesh generation and adaptivity for the solution of compressible viscous high speed flows. International Journal for Numerical Methods in Engineering 38, 1123–1148. [61] Henry, R. F., Walters, R. A., 1993. Geometrically based, automatic generator for irregular triangular networks. Communications in Numerical Methods in Engineering 9, 555–566. [62] Hill, D. F., 2007. Tidal modeling of Glacier Bay, Alaska. Accessed 5.5.2010. URL http://water.engr.psu.edu/hill/research/glba/software.html [63] Ho-Le, K., 1988. Finite element mesh generation methods: a review and classification. Computer Aided Design 20(1), 27–38. [64] Honghai, L., Lihwa, L., Burks-Copes, K. A., 2013. Modeling of coastal inundation, storm surge, and relative sea-level rise at naval station Norfolk, Norfolk, Virginia, U.S.A. Journal of Coastal Research 29(1), 18–30. [65] Horritt, M. S., 2000. Development of physically based meshes for two dimensional models of meandering channel flow. International Journal of Numerical Methods in Engineering 47, 2019–2037. [66] Huthnance, J., 1995. Circulation, exchange and water masses at the ocean margin: the role of physical processes at the shelf edge. Progress in Oceanography 35, 353–431. [67] Instituto Hidrogr´ afico, 2013. Tabelas de mar´e. Accessed 7.3.2013. URL http://www.hidrografico.pt [68] Jansen, K. E., Shephard, M. S., Beall, M. W., 2001. On anisotropic mesh generation and quality control in complex flow problems. 10th International Meshing Roundtable, Newport Beach, California. Sandia National Laboratories, 341–349. [69] Jones, J. E., Davies, A. M., 2005. An intercomparison between finite difference and finite element (TELEMAC) approaches to modelling west coast of Britain tides. Ocean Dynamics 55, 178–198. [70] Jones, N. L., Richards, D. R., 1992. Mesh generation for estuarine flow modeling. Journal of Waterway Port Coastal and Ocean Engineering 118, 599–614. [71] Kliem, N., Greenberg, D. A., 2003. Diagnostic simulations of the summer circulation in the Canadian arctic archipelago. Atmosphere-Ocean 41, 273–289. [72] Lapidus, L., Pinder, G. F., 1982. Numerical solution of partial differential equations in science and engineering. John Wiley and Sons, New York, 667pp.
182
BIBLIOGRAPHY
[73] Le Provost, C., Lyard, F., Molines, J. M., Genco, M. L., Rabilloud, F., 1998. A hydrodynamic ocean tide model improved by assimilating a satellite altimeter-derived data set. Journal of Geophysical Research 103, 5513–5529. [74] Le Provost, C., Vincent, P., 1986. Some tests of precision for a finite element model of ocean tides. Journal of Computational Physics 65, 273–291. [75] Lee, C. K., 1999. Automatic adaptive mesh generation using metric advancing front approach. Engineering Computations 16(2), 230263. [76] Lee, C. K., Hobbs, R. E., 1999. Automatic adaptive finite element mesh generation over arbitrary two-dimensional domain using the advancing front technique. Computers & Structures 71, 9–34. [77] Legrand, S., Deleersnijder, E., Delhez, E., Legat, V., 2007. Unstructured, anisotropic mesh generation for the nw european continental shelf, the continental slope and neighboring ocean. Continental Shelf Research 27, 1344–1356. [78] Legrand, S., Deleersnijder, E., Hanert, E., Legat, V., Wolanski, E., 2006. High-resolution unstructured meshes for hydrodynamic models of the great barrier reef, australia. Estuarine, Coastal and Shelf Science 68, 36–46. [79] Legrand, S., Legat, V., Deleersnijder, E., 2000. Delaunay mesh generation for an unstructured-grid ocean general circulation model. Ocean Modelling 2, 17–28. [80] Lo, S. H., 1985. A new mesh generation scheme for arbitrary planar domains. International Journal for Numerical Methods in Engineering 21, 1403–1426. [81] Lo, S. H., 2002. Finite element mesh generation and adaptive meshing. Progress in Structural Engineering and Materials 4, 381–399. [82] L¨ ohner, R., Cebral, J. R., 1999. Parallel advancing front grid generation. In: International Meshing Roundtable. pp. 67–74. [83] L¨ ohner, R., Parikh, P., 1988. Generation of three-dimensional unstructured grids by the advancing-front method. International Journal for Numerical Methods in Fluids 8, 1135– 1149. [84] Luettich, R. A., Westerink, J. J., 1995. Continental shelf scale convergence studies with a barotropic tidal model. Vol. 47. In: Quantitative Skill Assessment for Coastal Ocean Models. Coastal and Estuarine Studies, Lynch and Davies (Eds.). Washington DC. [85] Luettich, R. A., Westerink, J. J., 1995. Implementation and testing of elemental flooding and drying in the adcirc hydrodynamic model. Final Report for the Department of the Army, Contract Number DACW39-94-M-5869.
183
BIBLIOGRAPHY
[86] Luettich, R. A., Westerink, J. J., 1999. Elemental wetting and drying in the ADCIRC hydrodynamic model: upgrades and documentation for ADCIRC version 34.XX. Contractors Report, Department of the Army, US Army Corps of Engineers, Waterways Experiment Station, Vicksburg, MS., March 1999, 8pp. [87] Luettich, R. A., Westerink, J. J., 2004. Formulation and numerical implementation of the 2D/3D ADCIRC finite element model version 44.XX. [88] Luettich, R. A. J., Westerink, J. J., Scheffner, N. W., 1992. ADCIRC: an advanced three-dimensional circulation model for shelves, coasts and estuaries, report 1: theory and methodology of ADCIRC-2DDI and ADCIRC-3DL. Technical report DRP-92-6. U.S. Army Engineers Waterways Experiment Station, Vicksburg, MS, 137p. [89] Lynch, D., Justin, T. C. I., Naimie, C. E., Werner, F. E., 1995. Convergence studies of tidally rectified circulation on georges bank. In: Quantitative Skill Assessment for Coastal Ocean Models. Coastal and Estuarine Studies 47, 153–174. [90] Lynch, D. R., Gray, W. G., 1979. A wave equation model for finite element tidal computations. Computers & Fluids 7, 207–228. [91] Marujo da Silva, N. R. C., 2011. Storm surge hydrodynamic modelling. MSc thesis. Instituto Superior Tecnico. URL https://dspace.ist.utl.pt/bitstream/2295/1066510/1/dissertacao.pdf [92] Mavriplis, D. J., 1995a. Unstructured mesh generation and adaptivity. Technical Report ICASE 95-26, NASA Langley, Hampton VA. [93] Mavriplis, D. J., 1995b. An anvancing front Delaunay triangulation algorithm designed for robustness. Journal of Computational Physics 117(1), 90–101. [94] Mazzolari, A., Trigo-Teixeira, A., Ara´ ujo, M., 2012. Two dimensional unstructured mesh generation for shallow water models based on the multi-criteria strategy. In: 2as Jornadas de Engenharia Hidrografica. Lisboa, 20-22 junho de 2012. [95] Mazzolari, A., Trigo-Teixeira, A., Ara´ ujo, M., 2013. A multi-criteria meshing method applied to a shallow water model. Proceedings of the 12th International Coastal Symposium (Plymouth, England). Journal of Coastal Research, SI 65, 1170–1175. [96] Medeiros, S. C., Hagen, S. C., 2013. Review of wetting and drying algorithms for numerical tidal flow models. International journal for numerical methods in fluids, 71, 473–487. [97] Nash, J. E., Sutcliffe, J. V., 1970. River flow forecasting through conceptual models part I. A discussion of principles. Journal of Hydrology 10, 282–290. [98] Oden, J., Demkowicz, L., Rachowicz, W., Westermann, T., 1990. A posteriori error analysis in finite elements: the element residual method for symmetrizable problems with
184
BIBLIOGRAPHY
applications to compressible euler and Navier-stokes equations. Computer Methods in Applied Mechanics and Engineering 82, 183–203. [99] Ollivier-gooch, C., Freitag, L. A., 1997. Tetrahedral mesh improvement using swapping and smoothing. International Journal for Numerical Methods in Engineering 40, 3979– 4002. [100] Owen, S. J., 1998. A survey of unstructured mesh generation technology. Proceedings of the 7th International Meshing Roundtable, Sandia National Laboratories, 239–267. [101] Parrish, D. M., 2007. Target element sizes for finite element tidal models from a domain wide, localized truncation error analysis incorporating bottom stress and coriolis force. PhD Thesis. Department of Civil and Environmental Engineering, University of Central Florida. [102] Parrish, D. M., Hagen, S. C., 2007. 2D unstructured mesh generation for oceanic and coastal tidal models from a localized truncation error analysis with complex derivatives. International Journal of Computational Fluid Dynamics 21, 277–296. [103] Parrish, D. M., Hagen, S. C., 2009. Incorporating spatially variable bottom stress and Coriolis force into 2D, a posteriori , unstructured mesh generation for shallow water models. International Journal for Numerical Methods in Fluids 60, 237–261. [104] P´ebay, P. P., Baker, T. J., 2003. Analysis of triangle quality measures. Mathematics of Computation 72 (244), 1817–1839. URL http://dx.doi.org/10.1090/S0025-5718-03-01485-6 [105] Peraire, J., Vahdati, M., Morgan, K., Zienkiewicz, O. C., 1987. Adaptive remeshing for compressible flow computations. Journal of Computational Physics 72, 449–466. [106] Piggott, M. D., Farrell, P. E., Wilson, C. R., Gorman, G. J., Pain, C. C., 2009. Anisotropic mesh adaptivity for multi-scale ocean modelling. Philosophical Transactions of The Royal Society A: Mathematical, Physical and Engineering Sciences 367, 4591–4611. [107] Pinho, J. L. S., Vieira, J. M. P., 2007. Mathematical modelling of salt water intrusion in a northern portuguese estuary. Water in Celtic Countries: Quantity, Quality and Climate Variability. IAHS Publications 310, 277–287. [108] Pirzadeh, S., 1994. Viscous unstructured three-dimensional grids by the advancing-layers method. In: In proceedings of the 32nd Aerospace Sciences Meeting and Exhibit, AIAA94-0417, Reno, NV, January 1994. [109] Pirzadeh, S., 1996. Three-dimensional unstructured viscous grids by the advancing-layers method. AIAA Journal 34(1), 43–49.
185
BIBLIOGRAPHY
[110] Porta, G. M., Perotto, S., Ballio, F., 2012. Anisotropic mesh adaptation driven by a recovery-based error estimator for shallow water flow modeling. International Journal for Numerical Methods in Fluids 70 (3), 269–299. URL http://dx.doi.org/10.1002/fld.2688 [111] Power, P. W., Piggott, M. D., Fang, F., Gorman, G. J., Pain, C. C., Marshall, D., Goddard, A. J. H., Navon, I. M., 2006. Adjoint goal-based error norms for adaptive mesh ocean modelling. Ocean Modelling 15, 3–38. [112] Rappaz, M., Bellet, M., Deville, M. O., 1998. Numerical modeling in materials sciences and engineering. Springer Series in Computational Mathematics, 2. [113] Rebord˜ ao, I., Trigo-Teixeira, A. A., 2009. Tidal propagation in the Lima Estuary. Proceedings of the 10th International Coastal Symposium. Journal of Coastal Research SI 56, 1400–1404. [114] Remacle, J. F., Frazo, S. S., Li, X., Shephard, M. S., 2006. An adaptive discretization of shallow-water equations based on discontinuous Galerkin methods. International Journal for Numerical Methods in Fluids 52, 903–923. [115] Rohlf, F. J., Lopes, J. M. B., Delgado, L., 1987. GPACK: Graphics Package Users manual. Release 3.1. Dept. of Mechanical Engineering and Centro de Informatica. Instituto Superior Tecnico. Technical University of Lisbon, 232pp. [116] Salisbury, M. B., Hagen, S. C., Coggin, D., Bacopoulos, P., Atkinson, J., Roberts, H., 2011. Unstructured mesh development for the Big Bend region (Florida). Florida Watershed Journal 4, 11–14. [117] Schubert, J. E., Sanders, B. F., Smith, M. J., Wright, N. G., 2008. Unstructured mesh generation and landcover-based resistance for hydrodynamic modeling of urban flooding. Advances in Water Resources 31, 1603–1621. [118] Shen, J., Zhang, K., Xiao, C., Gong, W., 2006. Improved prediction of storm surge inundation with a high-resolution unstructured grid model. Journal of Coastal Research 22, 1309–1319. [119] Shewchuk, J., 1996. Triangle: engineering a 2D quality mesh generator and Delaunay triangulator. In: Proceedings of the First Workshop on Applied Computational Geometry. Association for Computing Machinery, Philadelphia, PA, 124133. [120] Smith, W. H. F., Sandwell, D. T., 1997. Global sea floor topography from satellite altimetry and ship depth soundings. Science 277, 1956–1962. [121] Surface-water Modeling System, 2013. SMS User Manual 11.1. Accessed 10.2.2013.
186
BIBLIOGRAPHY
[122] Thacker, W. C., 1980. A brief review of techniques for generating irregular computational grids. International Journal for Numerical Methods in Engineering 15, 1335–1341. [123] Thompson, J. F., Soni, B. K., Weatherill, N. G., 1999. Handbook of Grid Generation. CRC Press. [124] Thompson, J. F., Warsi, Z. U. A., Mastin, C. W., 1985. Numerical grid generation: foundations and applications. URL www.Hpc.msstate.edu/publications/gridbook/ [125] Trigo-Teixeira, A., 1994a. Finite element modeling of hydrodynamics in coastal zones. PhD Thesis. Department of Civil Engineering, University College of Swansea. [126] Trigo-Teixeira, A., 1994b. Hybrid mesh generation using algebraic interpolation and the advancing front concept. seminrio luso francs de modelao em hydrulica martima. coimbra, 11-13-1994. Journal de Estudos de Engenharia Civil SI 6, 216–222. [127] Trigo-Teixeira, A. A., 2008. Dredging works in the Port of Viana do Castelo and the impact on the Eiffel Bridge. A case study using finite element hydrodynamic modeling. Proceedings of the PIANC Mediterranean Days of Coastal and Port Engineering, Palermo. [128] Vale, L. M., Dias, J. M., 2011a. The effect of tidal regime and river flow on the hydrodynamics and salinity structure of the lima estuary: use of a numerical model to assist on estuary classification SI 64(2), 1604–1608. [129] Vale, L. M., Dias, J. M., 2011b. Coupling of a lagrangian particle tracking module to a numerical hydrodynamic model: simulation of pollution events inside an estuarine port area. Proceedings of the 11th International Coastal Symposium. Journal of Coastal Research SI 64(2), 1609–1613. [130] Vallet, M. G., Hecht, F., Mantel, B., 1991. Anisotropic control of mesh generation based upon Voronoi type methods. In: Numerical Grid Generation in Computational Fluid Dynamics and Related Fields. A. S. Arcilla, J. H¨ auser, P. R. Eiseman and J. F. Thompson (Eds.), Elsevier Science Publisher B.V. Amsterdam, 93-103. [131] Vorono¨ı, G., 1908. Nouvelles applications des param`etres continus `a la th´eorie des formes quadratiques. Journal f¨ ur die reine und angewandte Mathematik 134, 198–287. [132] Walters, R. A., Goring, D. G., Bell, R. G., 2001. Ocean tides around New Zealand. New Zealand Journal of Marine and Freshwater Research 35, 567–579. [133] Westerink, J. J., Luettich, R., Blain, C., Hagen, S., 1995. Surface elevation and circulation in continental margin waters. In: Carey, G. (Eds.), Finite Element Modeling of Environmental Problems. Wiley, Chichester, 39-60.
187
BIBLIOGRAPHY
[134] Westerink, J. J., Luettich, R. A., Feyen, J. C., Atkinson, J. H., Dawson, C., Roberts, H. J., Powell, M. D., Dunion, J. P., Kubatko, E. J., Pourtaheri, H., 2008. A basin to channel-scale unstructured grid hurricane storm surge model applied to Southern Louisiana. Monthly Weather Review 136, 833–864. [135] Westerink, J. J., Luettich, R. A., Muccino, J. C., 1994. Modelling tides in the western North Atlantic using unstructured graded grids. Tellus Series 46A, 178–199. [136] Westerink, J. J., Muccino, J. C., Luettich, R. A. J., 1992a. Resolution requirements for a tidal model of the Western North Atlantic and Gulf of Mexico. Proceedings of the 9th International Conference on Computational methods in Water Resources, Denver, CO. Russell et al. (Eds.), Computational Mechanics Publications, Southampton, UK. [137] Westerink, J. J., Muccino, J. C., Luettich, R. A. J., 1992b. Tidal and hurricane storm surge computations for the Western North Atlantic and Gulf of Mexico. Estuarine and Coastal Modeling II, M. Spaulding (Ed.), ASCE, 538-550. [138] Winteracter, Jun. 2009. The Fortran 9x GUI toolset. URL http://www.winteracter.com/index.htm [139] Wolanski, E., Brinkman, R., Spagnol, S., McAllister, F., Steinberg, C., Skirving, W., Deleersnijder, E., 2003. Merging scales in models of water circulation: perspectives from the great barrier reef. Advances in Coastal Modeling. Elsevier, 67, 411-429. [140] Xing, J., Davies, A. M., 1998. A three-dimensional model of internal tides on the MalinHebrides shelf and shelf edge. Journal of Geophysical Research 1032, 27821–27848. [141] Xu, H., Zhang, K., Shen, J., Li, Y., 2010. Storm surge simulation along the U.S. East and Gulf Coasts using a multi-scale numerical model approach. Ocean Dynamics 60, 1597–1619. [142] Yu, D., Lane, S. N., 2006. Urban fluvial flood modelling using a two-dimensional diffusionwave treatment, part 1: mesh resolution effects. Hydrological Processes 20, 1541–1565.
188
